How do you voxelize an isosurface of a 3D scalar field in Python?

Question

In a nutshell, I am interested in voxelizing a 2D object embedded in 3D with Python. At the end, I'd like to manipulate the voxelized data such as their coordinates. Along the way, I ended up using VTK but their little documentation threw me off. I was able to voxelize a surface but cannot retrieve their data.

Questions

1. Where is the voxel data (coordinates etc.) stored in the pyvista.core.pointset.UnstructuredGrid Class or vtkUnstructuredGrid Class? (pyvista is just a package built on VTK.)
2. Is there an easy way to voxelize a 2D object embedded in a 3D data? (Set 1 if the 2D object intersects with a voxel, 0 otherwise.)

---

In case you are interested, I provide a sample code that turns an isosurface(a mesh) of a scalar function to voxels. (The problem is that I don't know where the info about the voxels are stored in the class (pyvista.core.pointset.UnstructuredGrid Class which is supposed to be a fancy version of vtkUnstructuredGrid Class)

import numpy as np
import matplotlib.pyplot as plt
import pyvista
from skimage import measure

# Define positional grids (x, y, z)
x, y, z = np.linspace(-10, 10, 50), np.linspace(-10, 10, 50), np.linspace(-10, 10, 50)
xxx, yyy, zzz = np.meshgrid(x, y, z)
dx, dy, dz = x[1] - x[0], y[1] - y[0], z[1] - z[0] 

# Define sample data
m=6. # m=2: sphere, higher m: more cubic
rrr = np.abs((xxx**m + yyy**m + zzz**m) ** (1/m))
data = 5000 * np.exp(-0.5 * rrr)

# Find a isosurface with value 
isovalue = 1680
verts, faces_skimg, normals, values = measure.marching_cubes_lewiner(data, isovalue, spacing=(dx, dy, dz))

# Plot the isosurface
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:, 1], faces_skimg, verts[:, 2], cmap='Spectral', lw=1)

# Define some function for later
def convertFaces_skimg2vtk(faces_skimg):
    """This function fixes the data format of scikit-image to the one of vtk/pyvista."""
    nf, _ = faces_skimg.shape
    faces = []
    for i in range(nf):
        faces += [3]
        faces += list(faces_skimg[i, :])
    return np.asarray(faces)
# sckit-image and vtk use different formats to describe faces. 
faces_pyvista = convertFaces_skimg2vtk(faces_skimg)

# Create a mesh, then voxelize
mesh = pyvista.PolyData(verts, faces_pyvista)
voxels = pyvista.voxelize(mesh, density=mesh.length/200)

# Show the created mesh and voxels
p = pyvista.Plotter()
# p.add_mesh(mesh, color=True, show_edges=False, opacity=1) # show a mesh with VTK
p.add_mesh(voxels, color=True, show_edges=True, opacity=0.5) # show voxels with VTK
p.show()

This code gives the isosurface visualized using matplotlib and vtk (voxelized).

matplotlib- isosurface(a mesh)

vtk- isosurface(voxels)

Answer 1

I discovered that VTK actually has a class called vtkExtractSurface which works like a charm See this as well. Based on the documents, this method could be slower than what @Jens Munk suggests as it requires the voxelized data first. (It seems that I could directly use this method without voxelizing the isosurface but I am not sure how to do it.)

I leave the updated code just in case somebody needs it.

import numpy as np
import pyvista
from skimage import measure

# Define positional grids (x, y, z)
x, y, z = np.linspace(-10, 10, 50, endpoint=True), np.linspace(-10, 10, 50, endpoint=True), np.linspace(-10, 10, 50, endpoint=True)
xxx, yyy, zzz = np.meshgrid(x, y, z)
dx, dy, dz = x[1] - x[0], y[1] - y[0], z[1] - z[0] 

# Define sample data
m=2. # m=2: sphere, higher m: more cubic
rrr = np.abs((xxx**m + yyy**m + zzz**m) ** (1/m))
data = 5000 * np.exp(-0.5 * rrr)

# Find a isosurface with value 
value = 1680
verts, faces_skimg, normals, values = measure.marching_cubes_lewiner(data, value, spacing=(dx, dy, dz))
verts[:, 0] += np.min(yyy)
verts[:, 1] += np.min(xxx)
verts[:, 2] += np.min(zzz)

# Define a useful function for later.
def convertFaces_skimg2vtk(faces_skimg):
    """This function fixes the data format of scikit-image to the one of vtk/pyvista."""
    nf, _ = faces_skimg.shape
    faces = []
    for i in range(nf):
        faces += [3]
        faces += list(faces_skimg[i, :])
    return np.asarray(faces)

faces_pyvista = convertFaces_skimg2vtk(faces_skimg)

# Create a mesh, then voxelize
mesh = pyvista.PolyData(verts, faces_pyvista)
voxels = pyvista.voxelize(mesh, density=mesh.length/200)
slices = voxels.slice_orthogonal(x=0, y=0, z=0)

# solution
voxels['cell_ids'] = np.arange(voxels.n_cells)
cell_ids = voxels.extract_surface()['cell_ids']
surf_vox = voxels.extract_cells(cell_ids)
slices = surf_vox.slice_orthogonal(x=0, y=0, z=0)

# Show the created mesh and voxels
p = pyvista.Plotter()
# p.add_mesh(mesh, color=True, show_edges=False, opacity=1) # show a mesh with VTK
# p.add_mesh(voxels, color=True, show_edges=True, opacity=0.1) # show voxels with VTK
p.add_mesh(slices, color='skyblue')

p.add_bounding_box()
p.show()

Answer 2

I have created an example of what you are doing without the use of pyvista. It is orders of magnitude faster than your solution. It should not be hard to modify this to use two stencils to only address the volume in between two surfaces. Anyway, it runs really fast, so there is no need to not voxelize the volume.

Here goes

import numpy as np
import matplotlib.pyplot as plt
from skimage import measure

# Define positional grids (x, y, z)
x, y, z = np.linspace(-10, 10, 50), np.linspace(-10, 10, 50), np.linspace(-10, 10, 50)
xxx, yyy, zzz = np.meshgrid(x, y, z)
dx, dy, dz = x[1] - x[0], y[1] - y[0], z[1] - z[0] 

# Define sample data
m=6. # m=2: sphere, higher m: more cubic
rrr = np.abs((xxx**m + yyy**m + zzz**m) ** (1/m))
data = 5000 * np.exp(-0.5 * rrr)

# Find a isosurface with value 
isovalue = 1680
verts, faces_skimg, normals, values = measure.marching_cubes_lewiner(data, isovalue, spacing=(dx, dy, dz))

# Plot the isosurface
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(verts[:, 0], verts[:, 1], faces_skimg, verts[:, 2], cmap='Spectral', lw=1)

# Here goes VTK

import vtk

dxyz = 0.1

# Convenience function for displaying results side-by-side    
def vtkSubfigs(nrows=1, ncols=1, sharecamera=False):
  renderWindow = vtk.vtkRenderWindow()
  renderers = []
  camera = None
  if sharecamera:
    camera = vtk.vtkCamera()
  for irow in range(nrows):
    for icol in range(ncols):
      renderer = vtk.vtkRenderer()
      renderer.SetViewport(0.0 + icol*1.0/ncols,0.0 + irow*1.0/nrows,
                           (icol+1)*1.0/ncols,(irow+1)*1.0/nrows)
      if camera is not None:
        renderer.SetActiveCamera(camera)
      renderWindow.AddRenderer(renderer)
      renderers.append(renderer)
  return renderWindow, renderers

points = vtk.vtkPoints()

# Create the topology of the point (a vertex)
vertices = vtk.vtkCellArray()

# Add points and vertives
for i in range(0, len(verts)):
    p = verts[i].tolist()
    point_id = points.InsertNextPoint(p)
    vertices.InsertNextCell(1)
    vertices.InsertCellPoint(point_id)

# Create id list for a single face
def mkVtkIdList(it):
    vil = vtk.vtkIdList()
    for i in it:
        vil.InsertNextId(int(i))
    return vil
    
# Faces
polys = vtk.vtkCellArray()
for i in range(0, len(faces_skimg)):
  polys.InsertNextCell(mkVtkIdList(faces_skimg[i, :]))
    
# Create a poly data object
polydata = vtk.vtkPolyData()

# Set the points and vertices we created as the geometry and topology of the polydata
polydata.SetPoints(points)
polydata.SetVerts(vertices)
polydata.SetPolys(polys)

polydata.Modified()

# Visualize surface
mapper1 = vtk.vtkPolyDataMapper()
mapper1.SetInputData(polydata)

actor1 = vtk.vtkActor()
actor1.SetMapper(mapper1)

# White image to apply stencils to
whiteImage = vtk.vtkImageData()
bounds = polydata.GetBounds()
spacing = np.ones(3) * dxyz
whiteImage.SetSpacing(spacing)

dim = np.zeros(3,dtype=np.int32)
for i in range(3):
  dim[i] = np.ceil((bounds[i * 2 + 1] - bounds[i * 2]) / spacing[i]).astype(np.int32)
  dim[i] = dim[i] + 3
whiteImage.SetDimensions(dim)
whiteImage.SetExtent(0, dim[0] - 1, 0, dim[1] - 1, 0, dim[2] - 1)

origin = np.zeros(3, dtype=np.float64)
origin[0] = bounds[0] - 1.5*spacing[0]
origin[1] = bounds[2] - 1.5*spacing[1]
origin[2] = bounds[4] - 1.5*spacing[2]
whiteImage.SetOrigin(origin)
whiteImage.AllocateScalars(vtk.VTK_UNSIGNED_CHAR,1)

# fill the image with foreground voxels:
inval = 1
outval = 0
count = whiteImage.GetNumberOfPoints()

# Fast way of setting scalar data
data = np.ones(count,dtype='B')
data[:] = inval
a = vtk.vtkUnsignedCharArray()
a.SetArray(data, data.size, True)
whiteImage.GetPointData().SetScalars(a) # Don't delete data parameter

# polygonal data --> image stencil:
pol2stenc = vtk.vtkPolyDataToImageStencil()
pol2stenc.SetInputData(polydata)

pol2stenc.SetOutputOrigin(origin)
pol2stenc.SetOutputSpacing(spacing)
pol2stenc.SetOutputWholeExtent(whiteImage.GetExtent())
pol2stenc.Update()

# Cut the corresponding white image and set the background:
imgstenc = vtk.vtkImageStencil()
imgstenc.SetInputData(whiteImage)
imgstenc.SetStencilConnection(pol2stenc.GetOutputPort())

imgstenc.ReverseStencilOff()
imgstenc.SetBackgroundValue(outval)
imgstenc.Update()

liverImage = imgstenc.GetOutput()

startLabel = 1
endLabel = 1

# Pad the volume so that we can change the point data into cell data
extent = liverImage.GetExtent()

pad = vtk.vtkImageWrapPad()
pad.SetInputData(liverImage)
pad.SetOutputWholeExtent(extent[0], extent[1] + 1,
                         extent[2], extent[3] + 1,
                         extent[4], extent[5] + 1)
pad.Update()

# Copy the scalar point data of the volume into the scalar cell data
pad.GetOutput().GetCellData().SetScalars(
    imgstenc.GetOutput().GetPointData().GetScalars())

selector = vtk.vtkThreshold()
selector.SetInputArrayToProcess(0, 0, 0,
                                vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS,
                                vtk.vtkDataSetAttributes.SCALARS)
selector.SetInputConnection(pad.GetOutputPort())
selector.ThresholdBetween(startLabel, endLabel)
selector.Update()

# Shift the geometry by 1/2
transform = vtk.vtkTransform()
transform.Translate(-.5, -.5, -.5)

transformModel = vtk.vtkTransformFilter()
transformModel.SetTransform(transform)
transformModel.SetInputConnection(selector.GetOutputPort())

geometry = vtk.vtkGeometryFilter()
geometry.SetInputConnection(transformModel.GetOutputPort())

mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(geometry.GetOutputPort())
mapper.SetScalarRange(startLabel, endLabel)
mapper.SetScalarModeToUseCellData()
mapper.SetColorModeToMapScalars()

actor = vtk.vtkActor()
actor.SetMapper(mapper)

# Create two renderes inside window sharing camera
renderWindow, renderers = vtkSubfigs(ncols=2, sharecamera=True)
renderWindow.SetSize(800,800)

renderWindowInteractor = vtk.vtkRenderWindowInteractor()
renderWindowInteractor.SetRenderWindow(renderWindow)

renderers[0].AddActor(actor1)
prop = actor1.GetProperty()
prop.SetColor(vtk.vtkColor3d(1,0,0))
prop.SetOpacity(0.65)

renderers[0].SetBackground(1,1,1)

renderers[1].AddActor(actor)
renderers[1].SetBackground(1,1,1)
renderers[1].ResetCamera()

renderWindowInteractor.Initialize()

renderWindow.Render()
renderWindowInteractor.Start()