search

pyvista / pyvista

3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK)
https://docs.pyvista.org
MIT License
2.64k stars 482 forks source link

UniformGridFilters.gaussian_smooth() does not correctly set the scalar according to argument #1870

Open jweine opened 2 years ago

jweine commented 2 years ago

Describe the bug, what's wrong, and what you expect: Found behaviour:
When calling the dataset filter gaussian_smooth on a UniformGrid instance while specifying a non-active scalar as keyword argument, the smoothing is applied to the active scalar and saved to a data array with the specified name.

Expected behaviour:
When specifying the keyword argument 'scalars' the correct scalar field should be set as active and the gaussian smoothing should be applied accordingly. After applying the filter, the active scalar of the UniformGrid instance should be the same as before the application.

To Reproduce

The example below, produces a 3D gaussian random distribution and creates a 3D histogram, which is saved as cell data to a UniformGrid whose dimensions match the histogram bins (Image 1). Then the density is estimated by applying the gaussian_smooth filter, which produces the expected result (Image 2).
When creating an arbitrary additional data array and thereby changing the active scalar of the grid before filtering, the filter uses the dummy field (image 3).

import pyvista
import numpy as np

random_normal = np.random.normal(0., 1.5, size=(100000, 3))
uniform_grid = pyvista.UniformGrid((61, 61, 61), (0.1, 0.1, 0.1), (-3., -3., -3.))

extend = np.reshape(uniform_grid.extent, [3, 2])
bin_edges = [np.arange(0., uniform_grid.extent[2*i+1] + 1) * uniform_grid.spacing[i] + uniform_grid.origin[i] for i in range(3)]

histogram = np.histogramdd(random_normal, bins=bin_edges)

uniform_grid["histogram"] = histogram[0].reshape(-1)
uniform_grid["dummy"] = np.zeros_like(histogram[0].reshape(-1))
uniform_grid = uniform_grid.cell_data_to_point_data()

uniform_grid.set_active_scalars("histogram")  
correct_grid = uniform_grid.gaussian_smooth(std_dev=4, scalars="histogram")

uniform_grid.set_active_scalars("dummy")  
incorrect_grid = uniform_grid.gaussian_smooth(std_dev=4, scalars="histogram")

pyvista.close_all()
plotter = pyvista.Plotter(shape=(1, 3))
plotter.add_mesh(uniform_grid.slice_orthogonal(), scalars="histogram")
plotter.subplot(0, 1)
plotter.add_mesh(correct_grid.slice_orthogonal(), scalars="histogram")
plotter.subplot(0, 2)
plotter.add_mesh(incorrect_grid.slice_orthogonal(), scalars="histogram")
plotter.show()

Screenshots grafik

System Information:


--------------------------------------------------------------------------------
  Date: Sat Nov 27 17:49:30 2021 Mitteleuropäische Zeit

                OS : Windows
            CPU(s) : 12
           Machine : AMD64
      Architecture : 64bit
       Environment : Jupyter
        GPU Vendor : NVIDIA Corporation
      GPU Renderer : Quadro T2000/PCIe/SSE2
       GPU Version : 4.5.0 NVIDIA 471.41

  Python 3.8.11 (default, Aug  6 2021, 09:57:55) [MSC v.1916 64 bit (AMD64)]

           pyvista : 0.32.1
               vtk : 9.0.3
             numpy : 1.19.5
           imageio : 2.9.0
           appdirs : 1.4.4
            scooby : 0.5.7
            meshio : 4.4.6
        matplotlib : 3.4.3
         pyvistaqt : 0.5.0
             PyQt5 : 5.15.4
           IPython : 7.27.0
          colorcet : 1.0.0
        ipyvtklink : 0.2.1
             scipy : 1.7.1
        itkwidgets : 0.32.1
              tqdm : 4.62.3
--------------------------------------------------------------------------------
banesullivan commented 2 years ago

As far as I am aware, we are "correctly" setting up the underlying VTK algorithm to process on the specified scalars array with this line:

https://github.com/pyvista/pyvista/blob/d957ba67c2da56295e109ee72c65277c142dbf2f/pyvista/core/filters/uniform_grid.py#L49

which makes me wonder if the underlying vtkImageGaussianSmooth filter has a bug where it processes the active scalars rather than the set array to process. After some surface level digging, this does appear to be the case at this line in the vtkImageGaussianSmooth filter: https://github.com/Kitware/VTK/blob/6a9c565da01bcd6295d0bcbb66c0a9d0d2eaa69e/Imaging/General/vtkImageGaussianSmooth.cxx#L324

The filter is fetching a pointer to the active array, not what is specified to process by the algorithm -- I'm pretty sure this is a bug and would reccomend opening an issue with VTK over it.

To work around this in PyVista, we could have some sort of context manager to set the active scalars to what is specified in PyVista's gaussian_smooth method then reset when complete (because processing on an array should not set it as active). We may have to do some testing to see if this has any side effects -- I do know that if the source uniform grid is plotter somewhere, like in the pyvistaqt.BackgroundPlotter, setting and resetting the active scalars will have visual side effects.

jweine commented 2 years ago

I would be willing to contribute, and the Idea of locally setting and resetting the active scalar seems relatively straight forward. If welcome, I could look into seeting up some test cases.

banesullivan commented 2 years ago

Any contributions you could provide for this would be greatly appreciated!