Voxelize a Surface Mesh

Create a voxel model (like legos) of a closed surface or volumetric mesh.

This example also demonstrates how to compute an implicit distance from a bounding pyvista.PolyData surface.

# sphinx_gallery_thumbnail_number = 2
from pyvista import examples
import pyvista as pv
import numpy as np

# Load a surface to voxelize
surface = examples.download_foot_bones()
surface

PolyData	Information
N Cells	4204
N Points	2154
X Bounds	-5.633e+00, 5.633e+00
Y Bounds	-1.860e+00, 1.860e+00
Z Bounds	-2.125e+00, 2.126e+00
N Arrays	0

cpos = [(7.656346967151718, -9.802071079151158, -11.021236183314311),
 (0.2224512272564101, -0.4594554282112895, 0.5549738359311297),
 (-0.6279216753504941, -0.7513057097368635, 0.20311105371647392)]

surface.plot(cpos=cpos, opacity=0.75)

Out:

[(7.656346967151718, -9.802071079151158, -11.021236183314311),
 (0.2224512272564101, -0.4594554282112895, 0.5549738359311297),
 (-0.6279216753504943, -0.7513057097368636, 0.20311105371647395)]

Create a voxel model of the bounding surface

voxels = pv.voxelize(surface, density=surface.length/200)

p = pv.Plotter()
p.add_mesh(voxels, color=True, show_edges=True, opacity=0.5)
p.add_mesh(surface, color="lightblue", opacity=0.5)
p.show(cpos=cpos)

Out:

[(7.656346967151718, -9.802071079151158, -11.021236183314311),
 (0.2224512272564101, -0.4594554282112895, 0.5549738359311297),
 (-0.6279216753504943, -0.7513057097368636, 0.20311105371647395)]

We could even add a scalar field to that new voxel model in case we wanted to create grids for modelling. In this case, let’s add a scalar field for bone density noting:

voxels["density"] = np.full(voxels.n_cells, 3.65) # g/cc
voxels

Header	Data Arrays
UnstructuredGrid Information N Cells 93007 N Points 113156 X Bounds -5.633e+00, 5.584e+00 Y Bounds -1.860e+00, 1.858e+00 Z Bounds -2.125e+00, 2.097e+00 N Arrays 3	Name Field Type N Comp Min Max vtkOriginalPointIds Points int64 1 3.685e+03 7.283e+05 vtkOriginalCellIds Cells int64 1 3.624e+03 7.016e+05 density Cells float64 1 3.650e+00 3.650e+00

voxels.plot(scalars="density", cpos=cpos)

Out:

[(7.656346967151718, -9.802071079151158, -11.021236183314311),
 (0.2224512272564101, -0.4594554282112895, 0.5549738359311297),
 (-0.6279216753504943, -0.7513057097368636, 0.20311105371647395)]

A constant scalar field is kind of boring, so let’s get a little fancier by added a scalar field that varies by the distance from the bounding surface.

voxels.compute_implicit_distance(surface, inplace=True)
voxels

Header	Data Arrays
UnstructuredGrid Information N Cells 93007 N Points 113156 X Bounds -5.633e+00, 5.584e+00 Y Bounds -1.860e+00, 1.858e+00 Z Bounds -2.125e+00, 2.097e+00 N Arrays 4	Name Field Type N Comp Min Max vtkOriginalPointIds Points int64 1 3.685e+03 7.283e+05 implicit_distance Points float64 1 -6.951e-01 4.148e-01 vtkOriginalCellIds Cells int64 1 3.624e+03 7.016e+05 density Cells float64 1 3.650e+00 3.650e+00

contours = voxels.contour(6, scalars="implicit_distance")

p = pv.Plotter()
p.add_mesh(voxels, opacity=0.25, scalars="implicit_distance")
p.add_mesh(contours, opacity=0.5, scalars="implicit_distance")
p.show(cpos=cpos)

Out:

[(7.656346967151718, -9.802071079151158, -11.021236183314311),
 (0.2224512272564101, -0.4594554282112895, 0.5549738359311297),
 (-0.6279216753504943, -0.7513057097368636, 0.20311105371647395)]