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
- 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.)
- 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)
