{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fast Fourier Transform with Perlin Noise {#image_fft_perlin_example}\n\nThis example shows how to apply a Fast Fourier Transform (FFT) to a\n`pyvista.ImageData`{.interpreted-text role=\"class\"} using\n`pyvista.ImageDataFilters.fft`{.interpreted-text role=\"func\"} filter.\n\nHere, we demonstrate FFT usage by first generating Perlin noise using\n`pyvista.sample_function() <pyvista.core.utilities.features.sample_function>`{.interpreted-text\nrole=\"func\"} to sample\n`pyvista.perlin_noise <pyvista.core.utilities.features.perlin_noise>`{.interpreted-text\nrole=\"func\"}, and then performing FFT of the sampled noise to show the\nfrequency content of that noise.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n\nimport pyvista as pv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Generate Perlin Noise\n\nStart by generating some [Perlin\nNoise](https://en.wikipedia.org/wiki/Perlin_noise) as in\n`perlin_noise_2d_example`{.interpreted-text role=\"ref\"} example.\n\nNote that we are generating it in a flat plane and using a frequency of\n10 in the x direction and 5 in the y direction. The unit of frequency is\n`1/pixel`.\n\nAlso note that the dimensions of the image are powers of 2. This is\nbecause the FFT is much more efficient for arrays sized as a power of 2.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"freq = [10, 5, 0]\nnoise = pv.perlin_noise(1, freq, (0, 0, 0))\nxdim, ydim = (2**9, 2**9)\nsampled = pv.sample_function(noise, bounds=(0, 10, 0, 10, 0, 10), dim=(xdim, ydim, 1))\n\n# warp and plot the sampled noise\nwarped_noise = sampled.warp_by_scalar()\nwarped_noise.plot(show_scalar_bar=False, text='Perlin Noise', lighting=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Perform FFT of Perlin Noise\n\nNext, perform an FFT of the noise and plot the frequency content. For\nthe sake of simplicity, we will only plot the content in the first\nquadrant.\n\nNote the usage of `numpy.fft.fftfreq`{.interpreted-text role=\"func\"} to\nget the frequencies.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"sampled_fft = sampled.fft()\nfreq = np.fft.fftfreq(sampled.dimensions[0], sampled.spacing[0])\nmax_freq = freq.max()\n\n# only show the first quadrant\nsubset = sampled_fft.extract_subset((0, xdim // 2, 0, ydim // 2, 0, 0))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Plot the Frequency Domain\n\nNow, plot the noise in the frequency domain. Note how there is more high\nfrequency content in the x direction and this matches the frequencies\ngiven to\n`pyvista.perlin_noise <pyvista.core.utilities.features.perlin_noise>`{.interpreted-text\nrole=\"func\"}.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# scale to make the plot viewable\nsubset['scalars'] = np.abs(subset.active_scalars)\nwarped_subset = subset.warp_by_scalar(factor=0.0001)\n\npl = pv.Plotter(lighting='three lights')\npl.add_mesh(warped_subset, cmap='blues', show_scalar_bar=False)\npl.show_bounds(\n axes_ranges=(0, max_freq, 0, max_freq, 0, warped_subset.bounds[-1]),\n xtitle='X Frequency',\n ytitle='Y Frequency',\n ztitle='Amplitude',\n show_zlabels=False,\n color='k',\n font_size=26,\n)\npl.add_text('Frequency Domain of the Perlin Noise')\npl.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Low Pass Filter\n\nLet\\'s perform a low pass filter on the frequency content and then\nconvert it back into the space (pixel) domain by immediately applying a\nreverse FFT.\n\nWhen converting back, keep only the real content. The imaginary content\nhas no physical meaning in the physical domain. PyVista will drop the\nimaginary content, but will warn you of it.\n\nAs expected, we only see low frequency noise.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"low_pass = sampled_fft.low_pass(1.0, 1.0, 1.0).rfft()\nlow_pass['scalars'] = np.real(low_pass.active_scalars)\nwarped_low_pass = low_pass.warp_by_scalar()\nwarped_low_pass.plot(show_scalar_bar=False, text='Low Pass of the Perlin Noise', lighting=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# High Pass Filter\n\nThis time, let\\'s perform a high pass filter on the frequency content\nand then convert it back into the space (pixel) domain by immediately\napplying a reverse FFT.\n\nWhen converting back, keep only the real content. The imaginary content\nhas no physical meaning in the pixel domain.\n\nAs expected, we only see the high frequency noise content as the low\nfrequency noise has been attenuated.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"high_pass = sampled_fft.high_pass(1.0, 1.0, 1.0).rfft()\nhigh_pass['scalars'] = np.real(high_pass.active_scalars)\nwarped_high_pass = high_pass.warp_by_scalar()\nwarped_high_pass.plot(show_scalar_bar=False, text='High Pass of the Perlin Noise', lighting=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sum Low and High Pass\n\nShow that the sum of the low and high passes equals the original noise.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"grid = pv.ImageData(dimensions=sampled.dimensions, spacing=sampled.spacing)\ngrid['scalars'] = high_pass['scalars'] + low_pass['scalars']\n\nprint(\n 'Low and High Pass identical to the original:',\n np.allclose(grid['scalars'], sampled['scalars']),\n)\n\npl = pv.Plotter(shape=(1, 2))\npl.add_mesh(sampled.warp_by_scalar(), show_scalar_bar=False, lighting=False)\npl.add_text('Original Dataset')\npl.subplot(0, 1)\npl.add_mesh(grid.warp_by_scalar(), show_scalar_bar=False, lighting=False)\npl.add_text('Sum of the Low and High Passes')\npl.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Animate\n\nAnimate the variation of the cutoff frequency.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def warp_low_pass_noise(cfreq, scalar_ptp=np.ptp(sampled['scalars'])):\n \"\"\"Process the sampled FFT and warp by scalars.\"\"\"\n output = sampled_fft.low_pass(cfreq, cfreq, cfreq).rfft()\n\n # on the left: raw FFT magnitude\n output['scalars'] = output.active_scalars.real\n warped_raw = output.warp_by_scalar()\n\n # on the right: scale to fixed warped height\n output_scaled = output.translate((-11, 11, 0), inplace=False)\n output_scaled['scalars_warp'] = output['scalars'] / np.ptp(output['scalars']) * scalar_ptp\n warped_scaled = output_scaled.warp_by_scalar('scalars_warp')\n warped_scaled.active_scalars_name = 'scalars'\n # push center back to xy plane due to peaks near 0 frequency\n warped_scaled.translate((0, 0, -warped_scaled.center[-1]), inplace=True)\n\n return warped_raw + warped_scaled\n\n\n# Initialize the plotter and plot off-screen to save the animation as a GIF.\nplotter = pv.Plotter(notebook=False, off_screen=True)\nplotter.open_gif(\"low_pass.gif\", fps=8)\n\n# add the initial mesh\ninit_mesh = warp_low_pass_noise(1e-2)\nplotter.add_mesh(init_mesh, show_scalar_bar=False, lighting=False, n_colors=128)\nplotter.camera.zoom(1.3)\n\nfor freq in np.geomspace(1e-2, 10, 25):\n plotter.clear()\n mesh = warp_low_pass_noise(freq)\n plotter.add_mesh(mesh, show_scalar_bar=False, lighting=False, n_colors=128)\n plotter.add_text(f\"Cutoff Frequency: {freq:.2f}\", color=\"black\")\n plotter.write_frame()\n\n# write the last frame a few times to \"pause\" the gif\nfor _ in range(10):\n plotter.write_frame()\n\nplotter.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The left mesh in the above animation warps based on the raw values of\nthe FFT amplitude. This emphasizes how taking into account more and more\nfrequencies as the animation progresses, we recover a gradually larger\nproportion of the full noise sample. This is why the mesh starts\n\\\"flat\\\" and grows larger as the frequency cutoff is increased.\n\nIn contrast, the right mesh is always warped to the same visible height,\nirrespective of the cutoff frequency. This highlights how the typical\nwavelength (size of the features) of the Perlin noise decreases as the\nfrequency cutoff is increased since wavelength and frequency are\ninversely proportional.\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}