pyvista / pyvista

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

Include end bound in `pv.voxelize()` and `pv.voxelize_volume()` to avoid cropping the mesh #6609

Closed GuillaumeBroggi closed 1 month ago

GuillaumeBroggi commented 2 months ago

Describe the feature you would like to be added.

Voxelizing a mesh with pv.voxelize() or pv.voxelize_volume() will crop the mesh.

For instance, for a simple cube with the default parameters, the mesh bounds are [-0.5, 0.5]:

import pyvista as pv
cube = pv.Cube()
cube
HeaderData Arrays
PolyDataInformation
N Cells6
N Points8
N Strips0
X Bounds-5.000e-01, 5.000e-01
Y Bounds-5.000e-01, 5.000e-01
Z Bounds-5.000e-01, 5.000e-01
N Arrays3
NameFieldTypeN CompMinMax
NormalsPointsfloat323-1.000e+001.000e+00
TCoordsPointsfloat322-1.000e+001.000e+00
FaceIndexCellsint6410.000e+005.000e+00

After voxelization, the mesh bounds are [-0.5, 0.3], the missing 0.2 corresponding to a missing row of cells:

pv.voxelize(cube, density=0.2)
HeaderData Arrays
UnstructuredGridInformation
N Cells64
N Points125
X Bounds-5.000e-01, 3.000e-01
Y Bounds-5.000e-01, 3.000e-01
Z Bounds-5.000e-01, 3.000e-01
N Arrays2
NameFieldTypeN CompMinMax
vtkOriginalPointIdsPointsint6410.000e+001.240e+02
vtkOriginalCellIdsCellsint6410.000e+006.300e+01

One could expect the cube mesh to remain whole. Even if cropping is the intended behavior, there is no straightforward method to obtain a complete voxelized mesh of the cube.

As mentioned in https://github.com/pyvista/pyvista/issues/4312#issuecomment-1516994440, using np.linspace() instead of np.arange() would implement the proposed feature. Would a PR implementing this be welcome, or would it bring unwanted side effects?

Note

https://github.com/pyvista/pyvista/pull/6514 does not implement this behavior, instead the mesh is padded to enclose the cube:

pv.voxelize(cube, density=0.2, enclosed=True)
HeaderData Arrays
UnstructuredGridInformation
N Cells216
N Points343
X Bounds-6.000e-01, 6.000e-01
Y Bounds-6.000e-01, 6.000e-01
Z Bounds-6.000e-01, 6.000e-01
N Arrays2
NameFieldTypeN CompMinMax
vtkOriginalPointIdsPointsint6410.000e+003.420e+02
vtkOriginalCellIdsCellsint6410.000e+002.150e+02
tkoyama010 commented 2 months ago

Thank you for your suggestion. How about add it as an new option first? If the new method is better, we can make anothor breaking-change and make it the default. Pull requests are always welcome.

GuillaumeBroggi commented 2 months ago

In that case, I can start with the following signature that will not break anything:

voxelize(mesh, density=None, check_surface=True, enclosed=False, include_end_bound=False)

Maybe it could be made into a more generic function in the future, something like (just writing down some ideas):

voxelize(mesh, density=None, check_surface=True, range=["include", "exclude"])

where the behavior are defined for each bound:

The current behavior would be achieved with range=["include", "exclude"]. Using range=[float, float] or range=["include", "include"] would update the density if required to accommodate both bounds.

tkoyama010 commented 2 months ago

Awesome. Let's try it!

GuillaumeBroggi commented 1 month ago

Implemented in both #6610 and #6611