VTK Python中 PolyData和 ImageData的交集

编程语言 2026-07-10

我的目标是对分割后的体积进行网格光滑，并在更高分辨率下对生成的网格进行体素化。这样做的原因是我希望能够在等值面处计算法线，并通过确定相对于所计算的光滑网格的“外/内”比率，在新生成的体积中确定“组合体素”。以下是实现这些步骤的方法：

我有一个分割的图像堆栈，存储在一个 vtkImageData 中
我用vtkSurfaceNets3D生成等值面的网格（经平滑处理、孔洞被填充等）
我把它做成一个模板，以创建一个新体素化的体积，同时追踪哪些体素在SurfaceNets网格的内部和外部
现在关键的部分来了：我尝试确定哪些体素实际被网格相交，并试图计算它们相对的外/内比例。为此，我为每个体素考虑一个蒙特卡洛方法，在体素内取N 个点并用 vtkImplicitPolyDataDistance 计算它们到网格的距离。

但是，分数结果很奇怪，因为似乎存在某种不匹配。图像中所示的分数从黄色（1=在外部）变化到紫色（0=在内部），红点显示绘制切片处的网格轮廓。似乎存在某种偏移。

我无法确定这个偏移是由网格和图像之间的原点/间距不匹配引起，还是由于 vtkImplicitPolyDataDistance 方法的特定行为所致。

下面是我的参考代码。计算得到的分数值是 frac，由函数 voxel_fraction_mc 定义。

import tifffile
import vtk
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk
plt.close('all')

#%% Parametres
filename = 'sphere_stack.tiff'
resol = 0.02 # taille physique d'un pixel en mm

add_mat = False # add material around the image stack
matId_matiere = 2
# material thickness to add in all 3 directions (symetric)
ep_z = 2
ep_y = 2
ep_x = 2

convert2VTK = False

#%% Fetch data
img_stack = (tifffile.imread(filename)-1)
if add_mat:
    img_stack = np.pad(img_stack,((int(ep_z/resol),int(ep_z/resol)),(int(ep_y/resol),int(ep_y/resol)),(int(ep_x/resol),int(ep_x/resol))),'constant',constant_values=matId_matiere)
img_stack = img_stack.transpose((2,1,0))
# visualitation
id_img = img_stack.shape[2]//2
single_img = img_stack[:,:,id_img].astype(np.int8)

x1 = np.linspace(0.5*resol, img_stack.shape[0]*resol-0.5*resol, img_stack.shape[0])
y1 = np.linspace(0.5*resol, img_stack.shape[1]*resol-0.5*resol, img_stack.shape[1])
z1 = np.linspace(0.5*resol, img_stack.shape[2]*resol-0.5*resol, img_stack.shape[2])
xx1, yy1, zz1 = np.meshgrid(x1,y1,z1, indexing='ij')
centers1 = np.stack([xx1, yy1, zz1], axis=-1)

#%% SurfaceNets3D
from vtk.numpy_interface import dataset_adapter as dsa
def getMeshes(stack,resol,hole_fill_size=10*resol):
    """
        Uses VTK SurfaceNets3D to generate pixel coordinate meshes.

        stack: 3 dimensional labeled array. Expecting z, y, x dimensions, but
               shouldn't matter.
    """
    #VTK seems to use x as the first index.
    nx,ny,nz = stack.shape
    img = vtk.vtkImageData();
    img.SetDimensions(nx,ny,nz)
    img.SetSpacing(resol,resol,resol)
    img.SetOrigin(0,0,0)
    flat = stack.ravel(order='F')
    vtk_array = numpy_to_vtk(flat, deep=True)
    img.GetPointData().SetScalars(vtk_array)
    snets = vtk.vtkSurfaceNets3D()
    snets.SetInputData(img)
    snets.SetOutputMeshTypeToTriangles()
    # snets.SmoothingOff()
    snets.Update()
    # Hole fill
    fill = vtk.vtkFillHolesFilter()
    fill.SetInputData(snets.GetOutput())
    fill.SetHoleSize(hole_fill_size)
    fill.Update()
    # Smoothing
    smooth = vtk.vtkSmoothPolyDataFilter()
    smooth.SetInputData(fill.GetOutput())
    smooth.SetNumberOfIterations(20)
    smooth.Update()
    #The normals are not all the same direction.
    nrm = vtk.vtkPolyDataNormals()
    nrm.ConsistencyOn()
    nrm.AutoOrientNormalsOn()
    nrm.SetSplitting(False)
    nrm.SetInputDataObject(smooth.GetOutputDataObject(0) )
    nrm.Update()
    polys = nrm.GetOutput()
    #polys = snets.GetOutput()
    pda = dsa.WrapDataObject(polys)
    points = np.array(pda.GetPoints())
    return polys, pda, points, nrm

def vtk_to_triangles(pda):
    cells = pda.GetPolys()
    cell_array = cells.GetData()
    polys = vtk_to_numpy(cell_array)
    faces = []
    i = 0
    while i < len(polys):
        n = polys[i]
        face = polys[i+1:i+1+n]
        faces.append(face)
        i += n + 1
    return np.array(faces)
poly, pda, points, normals = getMeshes(img_stack,resol)  
#%% Voxelize
# Definition of the new voxelization
voxel_size = 0.8*resol
nx = int(img_stack.shape[0]*resol / voxel_size)
ny = int(img_stack.shape[1]*resol / voxel_size)
nz = int(img_stack.shape[2]*resol / voxel_size)
x = np.linspace(0.5*voxel_size, nx*voxel_size-0.5*voxel_size, nx)
y = np.linspace(0.5*voxel_size, ny*voxel_size-0.5*voxel_size, ny)
z = np.linspace(0.5*voxel_size, nz*voxel_size-0.5*voxel_size, nz)
xx, yy, zz = np.meshgrid(x,y,z, indexing='ij')
centers2 = np.stack([xx, yy, zz], axis=-1)
# Identification of voxels "inside", "outside" and at the frotier
image = vtk.vtkImageData()
image.SetDimensions(nx, ny, nz)
image.SetOrigin(0,0,0)
image.SetSpacing(voxel_size, voxel_size, voxel_size)
image.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 1)
image.GetPointData().GetScalars().Fill(1)

pol2stenc = vtk.vtkPolyDataToImageStencil()
pol2stenc.SetInputData(poly)
pol2stenc.SetOutputWholeExtent(image.GetExtent())
pol2stenc.SetOutputSpacing(voxel_size, voxel_size, voxel_size)
pol2stenc.Update()
inside = vtk.vtkImageStencil()
inside.SetInputData(image)
inside.SetStencilConnection(pol2stenc.GetOutputPort())
inside.ReverseStencilOff()
inside.SetBackgroundValue(0)  # extérieur = 0
inside.Update()
insideoutput = inside.GetOutput()
inside = insideoutput.GetPointData().GetScalars()
inside = vtk_to_numpy(inside)
inside = inside.reshape((nz, ny, nx)).transpose(2,1,0).astype(bool)
outside = vtk.vtkImageStencil()
outside.SetInputData(image)
outside.SetStencilConnection(pol2stenc.GetOutputPort())
outside.ReverseStencilOn()
outside.SetBackgroundValue(0)  # extérieur = 0
outside.Update()
outsideoutput = outside.GetOutput()
outside = outsideoutput.GetPointData().GetScalars()
outside = vtk_to_numpy(outside)
outside = outside.reshape((nz, ny, nx)).transpose(2,1,0).astype(bool)

implicit = vtk.vtkImplicitPolyDataDistance()
implicit.SetInput(poly)

dilate = vtk.vtkImageDilateErode3D()
dilate.SetInputData(insideoutput)
dilate.SetDilateValue(1)
dilate.SetErodeValue(0)
dilate.SetKernelSize(3,3,3)
dilate.Update()

dilate2 = vtk.vtkImageDilateErode3D()
dilate2.SetInputData(outsideoutput)
dilate2.SetDilateValue(1)
dilate2.SetErodeValue(0)
dilate2.SetKernelSize(3,3,3)
dilate2.Update()

erode = vtk.vtkImageDilateErode3D()
erode.SetInputData(insideoutput)
erode.SetDilateValue(0)
erode.SetErodeValue(1)
erode.SetKernelSize(3,3,3)
erode.Update()

boundary = vtk.vtkImageMathematics()
boundary.SetOperationToMultiply()
boundary.SetInput1Data(dilate2.GetOutput())
boundary.SetInput2Data(dilate.GetOutput())
boundary.Update()
intersection = vtk_to_numpy(boundary.GetOutput().GetPointData().GetScalars()).reshape((nz, ny, nx)).transpose(2,1,0).astype(bool)

#%% Voxel labels
# Volumic fraction of composite voxels
def voxel_fraction_mc(center, voxel_size, implicit, N=50): # fraction volumique monte carlo
    pts = np.random.uniform(-0.5,0.5,(N,3))*voxel_size + center
    d = [implicit.EvaluateFunction(p) for p in pts]
    return sum(di > 0 for di in d)/N # fraction INSIDE the mesh

frac = np.array([voxel_fraction_mc(pt, voxel_size, implicit) for pt in centers2[intersection].reshape(-1,3)])
centers2_label = np.ones(centers2.shape[:-1])
centers2_label[inside] = 0
centers2_label[outside] = 2
centers2_label[intersection] = 2*frac

#%% Visualization
# 2D
plane = vtk.vtkPlane()
plane.SetOrigin(0,0,centers2_label.shape[-1]//2*voxel_size)
plane.SetNormal(0, 0, 1)

# --- mesh slice ---
cutter = vtk.vtkCutter()
cutter.SetInputData(poly)  # ton mesh
cutter.SetCutFunction(plane)
cutter.Update()

# --- fetch point coordinates ---
cut_poly = cutter.GetOutput()
points_vtk = cut_poly.GetPoints()

points = vtk_to_numpy(points_vtk.GetData())

fig,axs = plt.subplots(1,2,figsize=(12,8))
axs[0].pcolormesh(single_img, edgecolors='black', linewidth=0.5, shading='nearest')
axs[0].set_aspect('equal')
axs[1].pcolormesh(centers2_label[:,:,centers2_label.shape[2]//2], edgecolors='black', linewidth=0.5, shading='nearest')
axs[1].set_aspect('equal')
axs[1].scatter(points[:,1]/voxel_size,points[:,0]/voxel_size,color='r',s=1)
axs[0].set_xticks(np.arange(-0.5, img_stack.shape[1], 1), minor=True)
axs[0].set_yticks(np.arange(-0.5, img_stack.shape[0], 1), minor=True)
axs[0].grid(which='minor', color='black', linestyle='-', linewidth=1,alpha=0.5)
axs[0].tick_params(which='both', bottom=False, left=False, labelbottom=False, labelleft=False)
axs[1].scatter(points[:,1]/voxel_size,points[:,0]/voxel_size,color='r',s=1)
axs[1].set_xticks(np.arange(-0.5, centers2_label.shape[1], 1), minor=True)
axs[1].set_yticks(np.arange(-0.5, centers2_label.shape[0], 1), minor=True)
axs[1].grid(which='minor', color='black', linestyle='-', linewidth=1,alpha=0.5)
axs[1].tick_params(which='both', bottom=False, left=False, labelbottom=False, labelleft=False)