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.

from __future__ import annotations

import numpy as np
import pyvista as pv

from pyvista import examples

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

PolyData (0x705009453280)
  N Cells:    4204
  N Points:   2154
  N Strips:   0
  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

PolyData 2,154 points 4,204 cells 158 KiB

Bounds⧉X [-5.633e+00, 5.633e+00]Y [-1.860e+00, 1.860e+00]Z [-2.125e+00, 2.126e+00]

Cellsfaces 4,204

cpos = pv.CameraPosition(
    position=(7.656346967151718, -9.802071079151158, -11.021236183314311),
    focal_point=(0.2224512272564101, -0.4594554282112895, 0.5549738359311297),
    viewup=(-0.6279216753504941, -0.7513057097368635, 0.20311105371647392),
)

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

Static Scene

Interactive Scene

Create a voxel model of the bounding surface

voxels = surface.voxelize()

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

Static Scene

Interactive Scene

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

UnstructuredGrid (0x705009452440)
  N Cells:    10478
  N Points:   15818
  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:   1

UnstructuredGrid 15,818 points 10,478 cells 1016 KiB

Bounds⧉X [-5.633e+00, 5.633e+00]Y [-1.860e+00, 1.860e+00]Z [-2.125e+00, 2.126e+00]

• Cell Data: (1)
  density active
  • density⧉
    scalar
    float64
    [3.650e+00, 3.650e+00]
    active

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

Static Scene

Interactive Scene

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

UnstructuredGrid (0x705009452440)
  N Cells:    10478
  N Points:   15818
  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:   2

UnstructuredGrid 15,818 points 10,478 cells 1.1 MiB

Bounds⧉X [-5.633e+00, 5.633e+00]Y [-1.860e+00, 1.860e+00]Z [-2.125e+00, 2.126e+00]

• Point Data: (1)
  • implicit_distance⧉
    scalar
    float64
    [-6.912e-01, 3.280e-01]
• Cell Data: (1)
  density active
  • density⧉
    scalar
    float64
    [3.650e+00, 3.650e+00]
    active

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

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

Static Scene

Interactive Scene

Tags: filter