Add image_fft and image_rfft

examples/01-filter/image-fft-perlin-noise.py

@@ -0,0 +1,206 @@
+"""
+.. _image_fft_perlin_example:
+
+Fast Fourier Transform with Perlin Noise
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This example shows how to apply a Fast Fourier Transform (FFT) to a
+:class:`pyvista.UniformGrid` using :func:`pyvista.UniformGridFilters.fft`
+filter.
+
+Here, we demonstrate FFT usage by first generating Perlin noise using
+:func:`pyvista.sample_function() <pyvista.core.imaging.sample_function>` to
+sample :func:`pyvista.perlin_noise <pyvista.core.common_data.perlin_noise>`,
+and then performing FFT of the sampled noise to show the frequency content of
+that noise.
+"""
+
+import numpy as np
+
+import pyvista as pv
+
+###############################################################################
+# Generate Perlin Noise
+# ~~~~~~~~~~~~~~~~~~~~~
+# Start by generating some `Perlin Noise
+# <https://en.wikipedia.org/wiki/Perlin_noise>`_ as in
+# :ref:`perlin_noise_2d_example` example.
+#
+# Note that we are generating it in a flat plane and using a frequency of 10 in
+# the x direction and 5 in the y direction. The unit of frequency is
+# ``1/pixel``.
+#
+# Also note that the dimensions of the image are powers of 2. This is because
+# the FFT is much more efficient for arrays sized as a power of 2.
+
+freq = [10, 5, 0]
+noise = pv.perlin_noise(1, freq, (0, 0, 0))
+xdim, ydim = (2**9, 2**9)
+sampled = pv.sample_function(noise, bounds=(0, 10, 0, 10, 0, 10), dim=(xdim, ydim, 1))
+
+# warp and plot the sampled noise
+warped_noise = sampled.warp_by_scalar()
+warped_noise.plot(show_scalar_bar=False, text='Perlin Noise', lighting=False)
+
+
+###############################################################################
+# Perform FFT of Perlin Noise
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Next, perform an FFT of the noise and plot the frequency content.
+# For the sake of simplicity, we will only plot the content in the first
+# quadrant.
+#
+# Note the usage of :func:`numpy.fft.fftfreq` to get the frequencies.
+
+sampled_fft = sampled.fft()
+freq = np.fft.fftfreq(sampled.dimensions[0], sampled.spacing[0])
+max_freq = freq.max()
+
+# only show the first quadrant
+subset = sampled_fft.extract_subset((0, xdim // 2, 0, ydim // 2, 0, 0))
+
+
+###############################################################################
+# Plot the Frequency Domain
+# ~~~~~~~~~~~~~~~~~~~~~~~~~
+# Now, plot the noise in the frequency domain. Note how there is more high
+# frequency content in the x direction and this matches the frequencies given
+# to :func:`pyvista.perlin_noise <pyvista.core.common_data.perlin_noise>`.
+
+# scale to make the plot viewable
+subset['scalars'] = np.abs(subset.active_scalars)
+warped_subset = subset.warp_by_scalar(factor=0.0001)
+
+pl = pv.Plotter(lighting='three lights')
+pl.add_mesh(warped_subset, cmap='blues', show_scalar_bar=False)
+pl.show_bounds(
+    axes_ranges=(0, max_freq, 0, max_freq, 0, warped_subset.bounds[-1]),
+    xlabel='X Frequency',
+    ylabel='Y Frequency',
+    zlabel='Amplitude',
+    show_zlabels=False,
+    color='k',
+    font_size=26,
+)
+pl.add_text('Frequency Domain of the Perlin Noise')
+pl.show()
+
+
+###############################################################################
+# Low Pass Filter
+# ~~~~~~~~~~~~~~~
+# Let's perform a low pass filter on the frequency content and then convert it
+# back into the space (pixel) domain by immediately applying a reverse FFT.
+#
+# When converting back, keep only the real content. The imaginary content has
+# no physical meaning in the physical domain. PyVista will drop the imaginary
+# content, but will warn you of it.
+#
+# As expected, we only see low frequency noise.
+
+low_pass = sampled_fft.low_pass(1.0, 1.0, 1.0).rfft()
+low_pass['scalars'] = np.real(low_pass.active_scalars)
+warped_low_pass = low_pass.warp_by_scalar()
+warped_low_pass.plot(show_scalar_bar=False, text='Low Pass of the Perlin Noise', lighting=False)
+
+
+###############################################################################
+# High Pass Filter
+# ~~~~~~~~~~~~~~~~
+# This time, let's perform a high pass filter on the frequency content and then
+# convert it back into the space (pixel) domain by immediately applying a
+# reverse FFT.
+#
+# When converting back, keep only the real content. The imaginary content has no
+# physical meaning in the pixel domain.
+#
+# As expected, we only see the high frequency noise content as the low
+# frequency noise has been attenuated.
+
+high_pass = sampled_fft.high_pass(1.0, 1.0, 1.0).rfft()
+high_pass['scalars'] = np.real(high_pass.active_scalars)
+warped_high_pass = high_pass.warp_by_scalar()
+warped_high_pass.plot(show_scalar_bar=False, text='High Pass of the Perlin Noise', lighting=False)
+
+
+###############################################################################
+# Sum Low and High Pass
+# ~~~~~~~~~~~~~~~~~~~~~
+# Show that the sum of the low and high passes equals the original noise.
+
+grid = pv.UniformGrid(dims=sampled.dimensions, spacing=sampled.spacing)
+grid['scalars'] = high_pass['scalars'] + low_pass['scalars']
+
+print(
+    'Low and High Pass identical to the original:', np.allclose(grid['scalars'], sampled['scalars'])
+)
+
+pl = pv.Plotter(shape=(1, 2))
+pl.add_mesh(sampled.warp_by_scalar(), show_scalar_bar=False, lighting=False)
+pl.add_text('Original Dataset')
+pl.subplot(0, 1)
+pl.add_mesh(grid.warp_by_scalar(), show_scalar_bar=False, lighting=False)
+pl.add_text('Sum of the Low and High Passes')
+pl.show()
+
+
+###############################################################################
+# Animate
+# ~~~~~~~
+# Animate the variation of the cutoff frequency.
+
+
+def warp_low_pass_noise(cfreq, scalar_ptp=sampled['scalars'].ptp()):
+    """Process the sampled FFT and warp by scalars."""
+    output = sampled_fft.low_pass(cfreq, cfreq, cfreq).rfft()
+
+    # on the left: raw FFT magnitude
+    output['scalars'] = output.active_scalars.real
+    warped_raw = output.warp_by_scalar()
+
+    # on the right: scale to fixed warped height
+    output_scaled = output.translate((-11, 11, 0), inplace=False)
+    output_scaled['scalars_warp'] = output['scalars'] / output['scalars'].ptp() * scalar_ptp
+    warped_scaled = output_scaled.warp_by_scalar('scalars_warp')
+    warped_scaled.active_scalars_name = 'scalars'
+    # push center back to xy plane due to peaks near 0 frequency
+    warped_scaled.translate((0, 0, -warped_scaled.center[-1]), inplace=True)
+
+    return warped_raw + warped_scaled
+
+
+# Initialize the plotter and plot off-screen to save the animation as a GIF.
+plotter = pv.Plotter(notebook=False, off_screen=True)
+plotter.open_gif("low_pass.gif", fps=8)
+
+# add the initial mesh
+init_mesh = warp_low_pass_noise(1e-2)
+plotter.add_mesh(init_mesh, show_scalar_bar=False, lighting=False, n_colors=128)
+plotter.camera.zoom(1.3)
+
+for freq in np.geomspace(1e-2, 10, 25):
+    plotter.clear()
+    mesh = warp_low_pass_noise(freq)
+    plotter.add_mesh(mesh, show_scalar_bar=False, lighting=False, n_colors=128)
+    plotter.add_text(f"Cutoff Frequency: {freq:.2f}", color="black")
+    plotter.write_frame()
+
+# write the last frame a few times to "pause" the gif
+for _ in range(10):
+    plotter.write_frame()
+
+plotter.close()
+
+
+###############################################################################
+# The left mesh in the above animation warps based on the raw values of the FFT
+# amplitude. This emphasizes how taking into account more and more frequencies
+# as the animation progresses, we recover a gradually larger proportion of the
+# full noise sample. This is why the mesh starts "flat" and grows larger as the
+# frequency cutoff is increased.
+#
+# In contrast, the right mesh is always warped to the same visible height,
+# irrespective of the cutoff frequency. This highlights how the typical
+# wavelength (size of the features) of the Perlin noise decreases as the
+# frequency cutoff is increased since wavelength and frequency are inversely
+# proportional.

examples/01-filter/image-fft.py

@@ -0,0 +1,109 @@
+"""
+.. _image_fft_example:
+
+Fast Fourier Transform
+~~~~~~~~~~~~~~~~~~~~~~
+
+This example shows how to apply a Fast Fourier Transform (FFT) to a
+:class:`pyvista.UniformGrid` using :func:`pyvista.UniformGridFilters.fft`
+filter.
+
+Here, we demonstrate FFT usage by denoising an image, effectively removing any
+"high frequency" content by performing a `low pass filter
+<https://en.wikipedia.org/wiki/Low-pass_filter>`_.
+
+This example was inspired by `Image denoising by FFT
+<https://scipy-lectures.org/intro/scipy/auto_examples/solutions/plot_fft_image_denoise.html>`_.
+
+"""
+
+import numpy as np
+
+import pyvista as pv
+from pyvista import examples
+
+###############################################################################
+# Load the example Moon landing image and plot it.
+
+image = examples.download_moonlanding_image()
+print(image.point_data)
+
+# Create a theme that we can reuse when plotting the image
+grey_theme = pv.themes.DocumentTheme()
+grey_theme.cmap = 'gray'
+grey_theme.show_scalar_bar = False
+grey_theme.axes.show = False
+image.plot(theme=grey_theme, cpos='xy', text='Unprocessed Moon Landing Image')
+
+
+###############################################################################
+# Apply FFT to the image
+# ~~~~~~~~~~~~~~~~~~~~~~
+# FFT will be applied to the active scalars, ``'PNGImage'``, the default
+# scalars name when loading a PNG image.
+#
+# The output from the filter is a complex array stored by the same name unless
+# specified using ``output_scalars_name``.
+
+fft_image = image.fft()
+fft_image.point_data
+
+
+###############################################################################
+# Plot the FFT of the image
+# ~~~~~~~~~~~~~~~~~~~~~~~~~
+# Plot the absolute value of the FFT of the image.
+#
+# Note that we are effectively viewing the "frequency" of the data in this
+# image, where the four corners contain the low frequency content of the image,
+# and the middle is the high frequency content of the image.
+
+fft_image.plot(
+    scalars=np.abs(fft_image.point_data['PNGImage']),
+    cpos="xy",
+    theme=grey_theme,
+    log_scale=True,
+    text='Moon Landing Image FFT',
+)
+
+
+###############################################################################
+# Remove the noise from the ``fft_image``
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Effectively, we want to remove high frequency (noisy) data from our image.
+# First, let's reshape by the size of the image.
+#
+# Next, perform a low pass filter by removing the middle 80% of the content of
+# the image. Note that the high frequency content is in the middle of the array.
+#
+# .. note::
+#    It is easier and more efficient to use the existing
+#    :func:`pyvista.UniformGridFilters.low_pass` filter. This section is here
+#    for demonstration purposes.
+
+ratio_to_keep = 0.10
+
+# modify the fft_image data
+width, height, _ = fft_image.dimensions
+data = fft_image['PNGImage'].reshape(height, width)  # note: axes flipped
+data[int(height * ratio_to_keep) : -int(height * ratio_to_keep)] = 0
+data[:, int(width * ratio_to_keep) : -int(width * ratio_to_keep)] = 0
+
+fft_image.plot(
+    scalars=np.abs(data),
+    cpos="xy",
+    theme=grey_theme,
+    log_scale=True,
+    text='Moon Landing Image FFT with Noise Removed',
+)
+
+
+###############################################################################
+# Convert to the spatial domain using reverse FFT
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Finally, convert the image data back to the "spatial" domain and plot it.
+
+
+rfft = fft_image.rfft()
+rfft['PNGImage'] = np.real(rfft['PNGImage'])
+rfft.plot(cpos="xy", theme=grey_theme, text='Processed Moon Landing Image')

examples/01-filter/sampling_functions_2d.py

@@ -1,13 +1,14 @@
 """
+.. _perlin_noise_2d_example:
+
 Sample Function: Perlin Noise in 2D
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Here we use :func:`pyvista.core.imaging.sample_function` to sample
 Perlin noise over a region to generate random terrain.
 
 Perlin noise is atype of gradient noise often used by visual effects
 artists to increase the appearance of realism in computer graphics.
-Source:
-https://en.wikipedia.org/wiki/Perlin_noise
+Source: `Perlin Noise Wikipedia <https://en.wikipedia.org/wiki/Perlin_noise>`_
 
 The development of Perlin Noise has allowed computer graphics artists
 to better represent the complexity of natural phenomena in visual

pyvista/_vtk.py

@@ -399,6 +399,12 @@
     except ImportError:  # pragma: no cover
         pass
 
+    from vtkmodules.vtkImagingFourier import (
+        vtkImageButterworthHighPass,
+        vtkImageButterworthLowPass,
+        vtkImageFFT,
+        vtkImageRFFT,
+    )
     from vtkmodules.vtkRenderingCore import (
         vtkActor,
         vtkActor2D,

pyvista/core/datasetattributes.py

@@ -476,11 +476,15 @@ def active_t_coords_name(self) -> Optional[str]:
 
     @active_t_coords_name.setter
     def active_t_coords_name(self, name: str) -> None:
+        if name is None:
+            self.SetActiveTCoords(None)
+            return
+
         self._raise_no_t_coords()
         dtype = self[name].dtype
         # only vtkDataArray subclasses can be set as active attributes
         if np.issubdtype(dtype, np.number) or dtype == bool:
-            self.SetActiveScalars(name)
+            self.SetActiveTCoords(name)
 
     def get_array(self, key: Union[str, int]) -> pyvista_ndarray:
         """Get an array in this object.
@@ -1150,6 +1154,10 @@ def active_scalars_name(self) -> Optional[str]:
 
     @active_scalars_name.setter
     def active_scalars_name(self, name: str) -> None:
+        # permit setting no active scalars
+        if name is None:
+            self.SetActiveScalars(None)
+            return
         self._raise_field_data_no_scalars_vectors()
         dtype = self[name].dtype
         # only vtkDataArray subclasses can be set as active attributes
@@ -1177,6 +1185,10 @@ def active_vectors_name(self) -> Optional[str]:
 
     @active_vectors_name.setter
     def active_vectors_name(self, name: str) -> None:
+        # permit setting no active
+        if name is None:
+            self.SetActiveVectors(None)
+            return
         self._raise_field_data_no_scalars_vectors()
         if name not in self:
             raise KeyError(f'DataSetAttribute does not contain "{name}"')
@@ -1313,6 +1325,10 @@ def active_normals_name(self) -> Optional[str]:
 
     @active_normals_name.setter
     def active_normals_name(self, name: str) -> None:
+        # permit setting no active
+        if name is None:
+            self.SetActiveNormals(None)
+            return
         self._raise_no_normals()
         self.SetActiveNormals(name)