如何在 Python 中对 3D 标量场的等值面进行体素化？

如何解决如何在 Python 中对 3D 标量场的等值面进行体素化？

简而言之，我对使用 Python 对嵌入 3D 中的 2D 对象进行体素化感兴趣。最后，我想操纵体素化数据，例如它们的坐标。一路上，我最终使用了 VTK，但他们的小文档让我失望。我能够对表面进行体素化，但无法检索其数据。

问题

1. pyvista.core.pointset.UnstructuredGrid 类或 vtkUnstructuredGrid 类中存储的体素数据（坐标等）在哪里？ （pyvista 只是一个构建在 VTK 上的包。）
2. 是否有一种简单的方法可以对嵌入在 3D 数据中的 2D 对象进行体素化？ （如果 2D 对象与体素相交，则设置为 1，否则设置为 0。）

---

如果您有兴趣，我提供了一个示例代码，可将标量函数的等值面（网格）转换为体素。 （问题是我不知道关于体素的信息在类中存储在哪里（pyvista.core.pointset.UnstructuredGrid Class，它应该是 vtkUnstructuredGrid Class 的一个奇特版本）

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

# Define positional grids (x,y,z)
x,z = np.linspace(-10,10,50),np.linspace(-10,50)
xxx,yyy,zzz = np.meshgrid(x,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,dz))

# Plot the isosurface
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot_trisurf(verts[:,0],verts[:,1],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,show_edges=True,opacity=0.5) # show voxels with VTK
p.show()

这段代码给出了使用 matplotlib 和 vtk（体素化）可视化的等值面。

ma​​tplotlib-等值面（网格）

vtk-等值面（体素）

解决方法

我已经创建了一个示例，说明您在不使用 pyvista 的情况下所做的事情。它比您的解决方案快几个数量级。修改它以使用两个模板来仅处理两个表面之间的体积应该不难。无论如何，它运行得非常快，所以没有必要不对音量进行体素化。

来了

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

# Define positional grids (x,y,z)
x,z = np.linspace(-10,10,50),np.linspace(-10,50)
xxx,yyy,zzz = np.meshgrid(x,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,dz))

# Plot the isosurface
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot_trisurf(verts[:,0],verts[:,1],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,dim[1] - 1,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,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))
prop.SetOpacity(0.65)

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

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

renderWindowInteractor.Initialize()

renderWindow.Render()
renderWindowInteractor.Start()

,

我发现 VTK 实际上有一个名为 vtkExtractSurface 的类，它的作用就像一个魅力 参见 this。根据文档，这种方法可能比@Jens Munk 建议的要慢，因为它首先需要体素化数据。 （似乎我可以直接使用这种方法而无需对等值面进行体素化，但我不知道该怎么做。）

我留下更新的代码，以防有人需要它。

import numpy as np
import pyvista
from skimage import measure

# Define positional grids (x,50,endpoint=True),endpoint=True)
xxx,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,value,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,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,show_edges=True,opacity=0.1) # show voxels with VTK
p.add_mesh(slices,color='skyblue')

p.add_bounding_box()
p.show()