diff --git a/examples/01-filter/image-fft-perlin-noise.py b/examples/01-filter/image-fft-perlin-noise.py new file mode 100644 index 0000000000..b83a0799db --- /dev/null +++ b/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() ` to +sample :func:`pyvista.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 +# `_ 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 `. + +# 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. diff --git a/examples/01-filter/image-fft.py b/examples/01-filter/image-fft.py new file mode 100644 index 0000000000..165ba9b6cc --- /dev/null +++ b/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 +`_. + +This example was inspired by `Image denoising by FFT +`_. + +""" + +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') diff --git a/examples/01-filter/sampling_functions_2d.py b/examples/01-filter/sampling_functions_2d.py index 27df2b8549..4e4f94e09e 100644 --- a/examples/01-filter/sampling_functions_2d.py +++ b/examples/01-filter/sampling_functions_2d.py @@ -1,4 +1,6 @@ """ +.. _perlin_noise_2d_example: + Sample Function: Perlin Noise in 2D ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Here we use :func:`pyvista.core.imaging.sample_function` to sample @@ -6,8 +8,7 @@ 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 `_ The development of Perlin Noise has allowed computer graphics artists to better represent the complexity of natural phenomena in visual diff --git a/pyvista/_vtk.py b/pyvista/_vtk.py index d44b65dcf3..e358ef9c18 100644 --- a/pyvista/_vtk.py +++ b/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, diff --git a/pyvista/core/datasetattributes.py b/pyvista/core/datasetattributes.py index 67fcdeef7e..68a55a8db3 100644 --- a/pyvista/core/datasetattributes.py +++ b/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) diff --git a/pyvista/core/filters/uniform_grid.py b/pyvista/core/filters/uniform_grid.py index 7163757108..31c6722b25 100644 --- a/pyvista/core/filters/uniform_grid.py +++ b/pyvista/core/filters/uniform_grid.py @@ -1,4 +1,4 @@ -"""Filters module with a class to manage filters/algorithms for uniform grid datasets.""" +"""Filters with a class to manage filters/algorithms for uniform grid datasets.""" import collections.abc import numpy as np @@ -7,6 +7,7 @@ from pyvista import _vtk, abstract_class from pyvista.core.filters import _get_output, _update_alg from pyvista.core.filters.data_set import DataSetFilters +from pyvista.errors import AmbiguousDataError, MissingDataError @abstract_class @@ -414,3 +415,346 @@ def image_threshold( # run the algorithm _update_alg(alg, progress_bar, 'Performing Image Thresholding') return _get_output(alg) + + def fft(self, output_scalars_name=None, progress_bar=False): + """Apply a fast Fourier transform (FFT) to the active scalars. + + The input can be real or complex data, but the output is always + :attr:`numpy.complex128`. The filter is fastest for images that have + power of two sizes. + + The filter uses a butterfly diagram for each prime factor of the + dimension. This makes images with prime number dimensions (i.e. 17x17) + much slower to compute. FFTs of multidimensional meshes (i.e volumes) + are decomposed so that each axis executes serially. + + The frequencies of the output assume standard order: along each axis + first positive frequencies are assumed from 0 to the maximum, then + negative frequencies are listed from the largest absolute value to + smallest. This implies that the corners of the grid correspond to low + frequencies, while the center of the grid corresponds to high + frequencies. + + Parameters + ---------- + output_scalars_name : str, optional + The name of the output scalars. By default, this is the same as the + active scalars of the dataset. + + progress_bar : bool, optional + Display a progress bar to indicate progress. + + Returns + ------- + pyvista.UniformGrid + :class:`pyvista.UniformGrid` with applied FFT. + + See Also + -------- + rfft: The reverse transform. + low_pass: Low-pass filtering of FFT output. + high_pass: High-pass filtering of FFT output. + + Examples + -------- + Apply FFT to an example image. + + >>> from pyvista import examples + >>> image = examples.download_moonlanding_image() + >>> fft_image = image.fft() + >>> fft_image.point_data # doctest:+SKIP + pyvista DataSetAttributes + Association : POINT + Active Scalars : PNGImage + Active Vectors : None + Active Texture : None + Active Normals : None + Contains arrays : + PNGImage complex128 (298620,) SCALARS + + See :ref:`image_fft_example` for a full example using this filter. + + """ + # check for active scalars, otherwise risk of segfault + if self.point_data.active_scalars_name is None: + try: + pyvista.set_default_active_scalars(self) + except MissingDataError: + raise MissingDataError('FFT filter requires point scalars.') from None + + # possible only cell scalars were made active + if self.point_data.active_scalars_name is None: + raise MissingDataError('FFT filter requires point scalars.') + + alg = _vtk.vtkImageFFT() + alg.SetInputDataObject(self) + _update_alg(alg, progress_bar, 'Performing Fast Fourier Transform') + output = _get_output(alg) + self._change_fft_output_scalars( + output, self.point_data.active_scalars_name, output_scalars_name + ) + return output + + def rfft(self, output_scalars_name=None, progress_bar=False): + """Apply a reverse fast Fourier transform (RFFT) to the active scalars. + + The input can be real or complex data, but the output is always + :attr:`numpy.complex128`. The filter is fastest for images that have power + of two sizes. + + The filter uses a butterfly diagram for each prime factor of the + dimension. This makes images with prime number dimensions (i.e. 17x17) + much slower to compute. FFTs of multidimensional meshes (i.e volumes) + are decomposed so that each axis executes serially. + + The frequencies of the input assume standard order: along each axis + first positive frequencies are assumed from 0 to the maximum, then + negative frequencies are listed from the largest absolute value to + smallest. This implies that the corners of the grid correspond to low + frequencies, while the center of the grid corresponds to high + frequencies. + + Parameters + ---------- + output_scalars_name : str, optional + The name of the output scalars. By default, this is the same as the + active scalars of the dataset. + + progress_bar : bool, optional + Display a progress bar to indicate progress. + + Returns + ------- + pyvista.UniformGrid + :class:`pyvista.UniformGrid` with the applied reverse FFT. + + See Also + -------- + fft: The direct transform. + low_pass: Low-pass filtering of FFT output. + high_pass: High-pass filtering of FFT output. + + Examples + -------- + Apply reverse FFT to an example image. + + >>> from pyvista import examples + >>> image = examples.download_moonlanding_image() + >>> fft_image = image.fft() + >>> image_again = fft_image.rfft() + >>> image_again.point_data # doctest:+SKIP + pyvista DataSetAttributes + Association : POINT + Active Scalars : PNGImage + Active Vectors : None + Active Texture : None + Active Normals : None + Contains arrays : + PNGImage complex128 (298620,) SCALARS + + See :ref:`image_fft_example` for a full example using this filter. + + """ + self._check_fft_scalars() + alg = _vtk.vtkImageRFFT() + alg.SetInputDataObject(self) + _update_alg(alg, progress_bar, 'Performing Reverse Fast Fourier Transform.') + output = _get_output(alg) + self._change_fft_output_scalars( + output, self.point_data.active_scalars_name, output_scalars_name + ) + return output + + def low_pass( + self, + x_cutoff, + y_cutoff, + z_cutoff, + order=1, + output_scalars_name=None, + progress_bar=False, + ): + """Perform a Butterworth low pass filter in the frequency domain. + + This filter requires that the :class:`UniformGrid` have a complex point + scalars, usually generated after the :class:`UniformGrid` has been + converted to the frequency domain by a :func:`UniformGridFilters.fft` + filter. + + A :func:`UniformGridFilters.rfft` filter can be used to convert the + output back into the spatial domain. This filter attenuates high + frequency components. Input and output are complex arrays with + datatype :attr:`numpy.complex128`. + + The frequencies of the input assume standard order: along each axis + first positive frequencies are assumed from 0 to the maximum, then + negative frequencies are listed from the largest absolute value to + smallest. This implies that the corners of the grid correspond to low + frequencies, while the center of the grid corresponds to high + frequencies. + + Parameters + ---------- + x_cutoff : double + The cutoff frequency for the x axis. + + y_cutoff : double + The cutoff frequency for the y axis. + + z_cutoff : double + The cutoff frequency for the z axis. + + order : int, optional + The order of the cutoff curve. Given from the equation + ``1 + (cutoff/freq(i, j))**(2*order)``. + + output_scalars_name : str, optional + The name of the output scalars. By default, this is the same as the + active scalars of the dataset. + + progress_bar : bool, optional + Display a progress bar to indicate progress. + + Returns + ------- + pyvista.UniformGrid + :class:`pyvista.UniformGrid` with the applied low pass filter. + + See Also + -------- + fft: Direct fast Fourier transform. + rfft: Reverse fast Fourier transform. + high_pass: High-pass filtering of FFT output. + + Examples + -------- + See :ref:`image_fft_perlin_example` for a full example using this filter. + + """ + self._check_fft_scalars() + alg = _vtk.vtkImageButterworthLowPass() + alg.SetInputDataObject(self) + alg.SetCutOff(x_cutoff, y_cutoff, z_cutoff) + alg.SetOrder(order) + _update_alg(alg, progress_bar, 'Performing Low Pass Filter') + output = _get_output(alg) + self._change_fft_output_scalars( + output, self.point_data.active_scalars_name, output_scalars_name + ) + return output + + def high_pass( + self, + x_cutoff, + y_cutoff, + z_cutoff, + order=1, + output_scalars_name=None, + progress_bar=False, + ): + """Perform a Butterworth high pass filter in the frequency domain. + + This filter requires that the :class:`UniformGrid` have a complex point + scalars, usually generated after the :class:`UniformGrid` has been + converted to the frequency domain by a :func:`UniformGridFilters.fft` + filter. + + A :func:`UniformGridFilters.rfft` filter can be used to convert the + output back into the spatial domain. This filter attenuates low + frequency components. Input and output are complex arrays with + datatype :attr:`numpy.complex128`. + + The frequencies of the input assume standard order: along each axis + first positive frequencies are assumed from 0 to the maximum, then + negative frequencies are listed from the largest absolute value to + smallest. This implies that the corners of the grid correspond to low + frequencies, while the center of the grid corresponds to high + frequencies. + + Parameters + ---------- + x_cutoff : double + The cutoff frequency for the x axis. + + y_cutoff : double + The cutoff frequency for the y axis. + + z_cutoff : double + The cutoff frequency for the z axis. + + order : int, optional + The order of the cutoff curve. Given from the equation + ``1/(1 + (cutoff/freq(i, j))**(2*order))``. + + output_scalars_name : str, optional + The name of the output scalars. By default, this is the same as the + active scalars of the dataset. + + progress_bar : bool, optional + Display a progress bar to indicate progress. + + Returns + ------- + pyvista.UniformGrid + :class:`pyvista.UniformGrid` with the applied high pass filter. + + See Also + -------- + fft: Direct fast Fourier transform. + rfft: Reverse fast Fourier transform. + low_pass: Low-pass filtering of FFT output. + + Examples + -------- + See :ref:`image_fft_perlin_example` for a full example using this filter. + + """ + self._check_fft_scalars() + alg = _vtk.vtkImageButterworthHighPass() + alg.SetInputDataObject(self) + alg.SetCutOff(x_cutoff, y_cutoff, z_cutoff) + alg.SetOrder(order) + _update_alg(alg, progress_bar, 'Performing High Pass Filter') + output = _get_output(alg) + self._change_fft_output_scalars( + output, self.point_data.active_scalars_name, output_scalars_name + ) + return output + + def _change_fft_output_scalars(self, dataset, orig_name, out_name): + """Modify the name and dtype of the output scalars for an FFT filter.""" + name = orig_name if out_name is None else out_name + pdata = dataset.point_data + if pdata.active_scalars_name != name: + pdata[name] = pdata.pop(pdata.active_scalars_name) + + # always view the datatype of the point_data as complex128 + dataset._association_complex_names['POINT'].add(name) + + def _check_fft_scalars(self): + """Check for complex active scalars. + + This is necessary for rfft, low_pass, and high_pass filters. + + """ + # check for complex active point scalars, otherwise the risk of segfault + if self.point_data.active_scalars_name is None: + possible_scalars = self.point_data.keys() + if len(possible_scalars) == 1: + self.set_active_scalars(possible_scalars[0], preference='point') + elif len(possible_scalars) > 1: + raise AmbiguousDataError( + 'There are multiple point scalars available. Set one to be ' + 'active with `point_data.active_scalars_name = `' + ) + else: + raise MissingDataError('FFT filters require point scalars.') + + scalars = self.point_data.active_scalars + + meets_req = scalars.ndim == 2 and scalars.shape[1] == 2 + if not meets_req: + raise ValueError( + 'Active scalars must be complex data for this filter, represented ' + 'as an array with `dtype=numpy.complex128`.' + ) diff --git a/pyvista/core/imaging.py b/pyvista/core/imaging.py index 6bd04bfd55..2daa90468f 100644 --- a/pyvista/core/imaging.py +++ b/pyvista/core/imaging.py @@ -111,6 +111,8 @@ def sample_function( >>> surf = pyvista.sample_function(noise, dim=(200, 200, 1)) >>> surf.plot() + See :ref:`perlin_noise_2d_example` for a full example using this function. + """ samp = _vtk.vtkSampleFunction() samp.SetImplicitFunction(function) diff --git a/pyvista/errors.py b/pyvista/errors.py index 6474eb87b5..d1ebe32ee1 100644 --- a/pyvista/errors.py +++ b/pyvista/errors.py @@ -4,6 +4,14 @@ class MissingDataError(ValueError): """Exception when data is missing, e.g. no active scalars can be set.""" + def __init__(self, message='No data available.'): + """Call the base class constructor with the custom message.""" + super().__init__(message) + class AmbiguousDataError(ValueError): """Exception when data is ambiguous, e.g. multiple active scalars can be set.""" + + def __init__(self, message="Multiple data available."): + """Call the base class constructor with the custom message.""" + super().__init__(message) diff --git a/pyvista/examples/downloads.py b/pyvista/examples/downloads.py index dafd28e4b0..a6b2e677fc 100644 --- a/pyvista/examples/downloads.py +++ b/pyvista/examples/downloads.py @@ -4088,3 +4088,37 @@ def download_cells_nd(load=True): # pragma: no cover """ return _download_and_read("cellsnd.ascii.inp", load=load) + + +def download_moonlanding_image(load=True): # pragma: no cover + """Download the Moon landing image. + + This is a noisy image originally obtained from `Scipy Lecture Notes + `_ and can be used to demonstrate a + low pass filter. + + See the `scipy-lectures license + `_ for more details + regarding this image's use and distribution. + + Parameters + ---------- + load : bool, optional + Load the dataset after downloading it when ``True``. Set this + to ``False`` and only the filename will be returned. + + Returns + ------- + pyvista.UniformGrid or str + ``DataSet`` or filename depending on ``load``. + + Examples + -------- + >>> from pyvista import examples + >>> dataset = examples.download_moonlanding_image() + >>> dataset.plot(cpos='xy', cmap='gray', background='w', show_scalar_bar=False) + + See :ref:`image_fft_example` for a full example using this dataset. + + """ + return _download_and_read('moonlanding.png', load=load) diff --git a/pyvista/plotting/plotting.py b/pyvista/plotting/plotting.py index 7402c39af7..48139a206d 100644 --- a/pyvista/plotting/plotting.py +++ b/pyvista/plotting/plotting.py @@ -7,7 +7,6 @@ import os import pathlib import platform -import re import textwrap from threading import Thread import time @@ -272,6 +271,7 @@ def __init__( self.last_image = None self._has_background_layer = False self._added_scalars = [] + self._prev_active_scalars = {} # set hidden line removal based on theme if self.theme.hidden_line_removal: @@ -2382,6 +2382,14 @@ def add_mesh( # Scalars formatting ================================================== added_scalar_info = None if scalars is not None: + # track the previous active scalars + if mesh.memory_address not in self._prev_active_scalars: + self._prev_active_scalars[mesh.memory_address] = ( + mesh, + mesh.point_data.active_scalars_name, + mesh.cell_data.active_scalars_name, + ) + show_scalar_bar, n_colors, clim, added_scalar_info = self.mapper.set_scalars( mesh, scalars, @@ -4632,30 +4640,20 @@ def _remove_added_scalars(self): as active scalars. """ - - def remove_and_reactivate_prior_scalars(dsattr, name): - """Remove ``name`` and reactivate prior scalars if applicable.""" - # reactivate prior active scalars - if name.endswith('-normed'): - orig_name = name[:-7] - if orig_name in dsattr: - dsattr.active_scalars_name = orig_name - elif name.endswith('-real'): - orig_name = name[:-5] - if orig_name in dsattr: - dsattr.active_scalars_name = orig_name - elif re.findall('-[0-9]+$', name): - # component - orig_scalars = re.sub('-[0-9]+$', '', name) - if orig_scalars in dsattr: - dsattr.active_scalars_name = orig_scalars - - dsattr.pop(name, None) - + # remove the added scalars for mesh, (name, assoc) in self._added_scalars: dsattr = mesh.point_data if assoc == 'point' else mesh.cell_data - remove_and_reactivate_prior_scalars(dsattr, name) + dsattr.pop(name, None) + + # reactivate prior active scalars + for mesh, point_name, cell_name in self._prev_active_scalars.values(): + if point_name is not None and mesh.point_data.active_scalars_name != point_name: + mesh.point_data.active_scalars_name = point_name + if cell_name is not None and mesh.cell_data.active_scalars_name != cell_name: + mesh.cell_data.active_scalars_name = cell_name + self._added_scalars = [] + self._prev_active_scalars = {} def __del__(self): """Delete the plotter.""" diff --git a/pyvista/plotting/renderer.py b/pyvista/plotting/renderer.py index 700b3cb48a..0c3a7d3e1d 100644 --- a/pyvista/plotting/renderer.py +++ b/pyvista/plotting/renderer.py @@ -1116,15 +1116,21 @@ def show_bounds( Bolds axis labels and numbers. Default ``True``. font_size : float, optional - Sets the size of the label font. Defaults to 16. + Sets the size of the label font. Defaults to + :attr:`pyvista.global_theme.font.size + `. font_family : str, optional Font family. Must be either ``'courier'``, ``'times'``, - or ``'arial'``. + or ``'arial'``. Defaults to :attr:`pyvista.global_theme.font.family + `. color : color_like, optional - Color of all labels and axis titles. Default white. - Either a string, rgb list, or hex color string. For + Color of all labels and axis titles. Defaults to + :attr:`pyvista.global_theme.font.color + `. + + Either a string, RGB list, or hex color string. For example: * ``color='white'`` @@ -1374,16 +1380,25 @@ def show_bounds( # set font font_family = parse_font_family(font_family) - for i in range(3): - cube_axes_actor.GetTitleTextProperty(i).SetFontSize(font_size) - cube_axes_actor.GetTitleTextProperty(i).SetColor(color.float_rgb) - cube_axes_actor.GetTitleTextProperty(i).SetFontFamily(font_family) - cube_axes_actor.GetTitleTextProperty(i).SetBold(bold) - - cube_axes_actor.GetLabelTextProperty(i).SetFontSize(font_size) - cube_axes_actor.GetLabelTextProperty(i).SetColor(color.float_rgb) - cube_axes_actor.GetLabelTextProperty(i).SetFontFamily(font_family) - cube_axes_actor.GetLabelTextProperty(i).SetBold(bold) + props = [ + cube_axes_actor.GetTitleTextProperty(0), + cube_axes_actor.GetTitleTextProperty(1), + cube_axes_actor.GetTitleTextProperty(2), + cube_axes_actor.GetLabelTextProperty(0), + cube_axes_actor.GetLabelTextProperty(1), + cube_axes_actor.GetLabelTextProperty(2), + ] + + for prop in props: + prop.SetColor(color.float_rgb) + prop.SetFontFamily(font_family) + prop.SetBold(bold) + + # Note: font_size does nothing as a property, use SetScreenSize instead + # Here, we normalize relative to 12 to give the user an illusion of + # just changing the font size relative to a font size of 12. 10 is used + # here since it's the default "screen size". + cube_axes_actor.SetScreenSize(font_size / 12 * 10.0) self.add_actor(cube_axes_actor, reset_camera=False, pickable=False, render=render) self.cube_axes_actor = cube_axes_actor diff --git a/tests/conftest.py b/tests/conftest.py index 1455a3f7d4..182ca2bf64 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -103,6 +103,13 @@ def pointset(): return pyvista.PointSet(points) +@fixture() +def noise_2d(): + freq = [10, 5, 0] + noise = pyvista.perlin_noise(1, freq, (0, 0, 0)) + return pyvista.sample_function(noise, bounds=(0, 10, 0, 10, 0, 10), dim=(2**4, 2**4, 1)) + + def pytest_addoption(parser): parser.addoption("--reset_image_cache", action='store_true', default=False) parser.addoption("--ignore_image_cache", action='store_true', default=False) diff --git a/tests/test_datasetattributes.py b/tests/test_datasetattributes.py index 67e1006302..e7b0835220 100644 --- a/tests/test_datasetattributes.py +++ b/tests/test_datasetattributes.py @@ -128,6 +128,9 @@ def test_active_scalars_name(sphere): sphere.point_data[key] = range(sphere.n_points) assert sphere.point_data.active_scalars_name == key + sphere.point_data.active_scalars_name = None + assert sphere.point_data.active_scalars_name is None + def test_set_scalars(sphere): scalars = np.array(sphere.n_points) @@ -189,6 +192,10 @@ def test_set_vectors(hexbeam): hexbeam.point_data.set_vectors(vectors, 'my-vectors') assert np.allclose(hexbeam.point_data.active_vectors, vectors) + # check clearing + hexbeam.point_data.active_vectors_name = None + assert hexbeam.point_data.active_vectors_name is None + def test_set_invalid_vectors(hexbeam): # verify non-vector data does not become active vectors @@ -197,6 +204,17 @@ def test_set_invalid_vectors(hexbeam): hexbeam.point_data.set_vectors(not_vectors, 'my-vectors') +def test_set_tcoords_name(): + mesh = pyvista.Cube() + old_name = mesh.point_data.active_t_coords_name + assert mesh.point_data.active_t_coords_name is not None + mesh.point_data.active_t_coords_name = None + assert mesh.point_data.active_t_coords_name is None + + mesh.point_data.active_t_coords_name = old_name + assert mesh.point_data.active_t_coords_name == old_name + + def test_set_bitarray(hexbeam): """Test bitarrays are properly loaded and represented in datasetattributes.""" hexbeam.clear_data() @@ -459,6 +477,9 @@ def test_normals_get(plane): plane_w_normals = plane.compute_normals() assert np.array_equal(plane_w_normals.point_data.active_normals, plane_w_normals.point_normals) + plane.point_data.active_normals_name = None + assert plane.point_data.active_normals_name is None + def test_normals_set(): plane = pyvista.Plane(i_resolution=1, j_resolution=1) diff --git a/tests/test_grid.py b/tests/test_grid.py index b08f01a71b..e3d72cd22f 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -9,6 +9,7 @@ import pyvista from pyvista import examples from pyvista._vtk import VTK9 +from pyvista.errors import AmbiguousDataError, MissingDataError from pyvista.plotting import system_supports_plotting from pyvista.utilities.misc import PyvistaDeprecationWarning @@ -920,6 +921,71 @@ def test_cast_uniform_to_rectilinear(): assert rectilinear.bounds == grid.bounds +def test_fft_and_rfft(noise_2d): + grid = pyvista.UniformGrid(dims=(10, 10, 1)) + with pytest.raises(MissingDataError, match='FFT filter requires point scalars'): + grid.fft() + + grid['cell_data'] = np.arange(grid.n_cells) + with pytest.raises(MissingDataError, match='FFT filter requires point scalars'): + grid.fft() + + name = noise_2d.active_scalars_name + noise_fft = noise_2d.fft() + assert noise_fft[name].dtype == np.complex128 + + full_pass = noise_2d.fft().rfft() + assert full_pass[name].dtype == np.complex128 + + # expect FFT and and RFFT to transform from time --> freq --> time domain + assert np.allclose(noise_2d['scalars'], full_pass[name].real) + assert np.allclose(full_pass[name].imag, 0) + + output_scalars_name = 'out_scalars' + # also, disable active scalars to check if it will be automatically set + noise_2d.active_scalars_name = None + noise_fft = noise_2d.fft(output_scalars_name=output_scalars_name) + assert output_scalars_name in noise_fft.point_data + + noise_fft = noise_2d.fft() + noise_fft_inactive_scalars = noise_fft.copy() + noise_fft_inactive_scalars.active_scalars_name = None + full_pass = noise_fft_inactive_scalars.rfft() + assert np.allclose(full_pass.active_scalars, noise_fft.rfft().active_scalars) + + +def test_fft_low_pass(noise_2d): + name = noise_2d.active_scalars_name + noise_no_scalars = noise_2d.copy() + noise_no_scalars.clear_data() + with pytest.raises(MissingDataError, match='FFT filters require point scalars'): + noise_no_scalars.low_pass(1, 1, 1) + + noise_too_many_scalars = noise_no_scalars.copy() + noise_too_many_scalars.point_data.set_array(np.arange(noise_2d.n_points), 'a') + noise_too_many_scalars.point_data.set_array(np.arange(noise_2d.n_points), 'b') + with pytest.raises(AmbiguousDataError, match='There are multiple point scalars available'): + noise_too_many_scalars.low_pass(1, 1, 1) + + with pytest.raises(ValueError, match='must be complex data'): + noise_2d.low_pass(1, 1, 1) + + out_zeros = noise_2d.fft().low_pass(0, 0, 0) + assert np.allclose(out_zeros[name][1:], 0) + + out = noise_2d.fft().low_pass(1, 1, 1) + assert not np.allclose(out[name][1:], 0) + + +def test_fft_high_pass(noise_2d): + name = noise_2d.active_scalars_name + out_zeros = noise_2d.fft().high_pass(100000, 100000, 100000) + assert np.allclose(out_zeros[name], 0) + + out = noise_2d.fft().high_pass(10, 10, 10) + assert not np.allclose(out[name], 0) + + @pytest.mark.parametrize('binary', [True, False]) @pytest.mark.parametrize('extension', ['.vtk', '.vtr']) def test_save_rectilinear(extension, binary, tmpdir): diff --git a/tests/test_plotter.py b/tests/test_plotter.py index b3876d268b..24e3fc0863 100644 --- a/tests/test_plotter.py +++ b/tests/test_plotter.py @@ -13,7 +13,7 @@ def test_plotter_image(): plotter = pyvista.Plotter() - with pytest.raises(AttributeError, match='not yet been set up'): + with pytest.raises(AttributeError, match="not yet been set up"): plotter.image @@ -72,17 +72,17 @@ def test_prepare_smooth_shading_texture(globe): """Test edge cases for smooth shading""" mesh, scalars = _plotting.prepare_smooth_shading(globe, None, True, True, False, None) assert scalars is None - assert 'Normals' in mesh.point_data - assert 'Texture Coordinates' in mesh.point_data + assert "Normals" in mesh.point_data + assert "Texture Coordinates" in mesh.point_data def test_prepare_smooth_shading_not_poly(hexbeam): """Test edge cases for smooth shading""" - scalars_name = 'sample_point_scalars' + scalars_name = "sample_point_scalars" scalars = hexbeam.point_data[scalars_name] mesh, scalars = _plotting.prepare_smooth_shading(hexbeam, scalars, False, True, True, None) - assert 'Normals' in mesh.point_data + assert "Normals" in mesh.point_data expected_mesh = hexbeam.extract_surface().compute_normals( cell_normals=False, @@ -115,11 +115,11 @@ def test_remove_scalars_single(sphere, hexbeam): assert len(pl._added_scalars) == 2 for mesh, (name, assoc) in pl._added_scalars: - assert name == 'Data' + assert name == "Data" if mesh is sphere: - assert assoc == 'point' + assert assoc == "point" else: - assert assoc == 'cell' + assert assoc == "cell" pl.close() assert pl._added_scalars == [] @@ -132,7 +132,7 @@ def test_remove_scalars_complex(sphere): """Test plotting complex data.""" data = np.arange(sphere.n_points, dtype=np.complex128) data += np.linspace(0, 1, sphere.n_points) * -1j - point_data_name = 'data' + point_data_name = "data" sphere[point_data_name] = data pl = pyvista.Plotter() @@ -155,7 +155,7 @@ def test_remove_scalars_normalized(sphere): assert sphere.n_arrays == 0 # test scalars are removed for normalized multi-component - point_data_name = 'data' + point_data_name = "data" sphere[point_data_name] = np.random.random((sphere.n_points, 3)) pl = pyvista.Plotter() pl.add_mesh(sphere, scalars=point_data_name) @@ -165,7 +165,7 @@ def test_remove_scalars_normalized(sphere): def test_remove_scalars_component(sphere): - point_data_name = 'data' + point_data_name = "data" sphere[point_data_name] = np.random.random((sphere.n_points, 3)) pl = pyvista.Plotter() pl.add_mesh(sphere, scalars=point_data_name, component=0) @@ -181,7 +181,7 @@ def test_remove_scalars_component(sphere): def test_remove_scalars_rgba(sphere): - point_data_name = 'data' + point_data_name = "data" sphere[point_data_name] = np.random.random(sphere.n_points) pl = pyvista.Plotter() pl.add_mesh(sphere, scalars=point_data_name, opacity=point_data_name) @@ -196,9 +196,39 @@ def test_remove_scalars_rgba(sphere): # assert sphere.point_data.active_scalars_name == point_data_name +def test_active_scalars_remain(sphere, hexbeam): + """Ensure active scalars remain active despite plotting different scalars.""" + point_data_name = "point_data" + cell_data_name = "cell_data" + sphere[point_data_name] = np.random.random(sphere.n_points) + hexbeam[cell_data_name] = np.random.random(hexbeam.n_cells) + assert sphere.point_data.active_scalars_name == point_data_name + assert hexbeam.cell_data.active_scalars_name == cell_data_name + + pl = pyvista.Plotter() + pl.add_mesh(sphere, scalars=range(sphere.n_points)) + pl.add_mesh(hexbeam, scalars=range(hexbeam.n_cells)) + pl.close() + + assert sphere.point_data.active_scalars_name == point_data_name + assert hexbeam.cell_data.active_scalars_name == cell_data_name + + def test_no_added_with_scalar_bar(sphere): - point_data_name = 'data' + point_data_name = "data" sphere[point_data_name] = np.random.random(sphere.n_points) pl = pyvista.Plotter() pl.add_mesh(sphere, scalar_bar_args={"title": "some_title"}) assert sphere.n_arrays == 1 + + +def test_add_multiple(sphere): + point_data_name = 'data' + sphere[point_data_name] = np.random.random(sphere.n_points) + pl = pyvista.Plotter() + pl.add_mesh(sphere) + pl.add_mesh(sphere, scalars=np.arange(sphere.n_points)) + pl.add_mesh(sphere, scalars=np.arange(sphere.n_cells)) + pl.add_mesh(sphere, scalars='data') + pl.show() + assert sphere.n_arrays == 1