VTK Python中PolyData与ImageData相交问题：复合体素比例计算偏移异常排查

我的核心目标是对分割后的体数据完成两项关键操作：

• 对分割体数据生成的网格进行平滑处理
• 将平滑后的网格重新体素化为更高分辨率的体数据

这么做的最终目的是：

1. 计算等值面的法线
2. 通过识别新体数据中体素相对于平滑网格的外部/内部比例，找出“复合体素”

具体执行步骤如下：

1. 读取存储在vtkImageData中的分割图像栈
2. 使用vtkSurfaceNets3D生成等值面网格，同时完成孔洞填充、平滑等后处理
3. 将处理后的网格转换为模板，生成新的体素化数据，并记录每个体素相对于SurfaceNets网格的内外状态

我当前卡在了体素与网格相交判断及内外比例计算这一步：
我采用蒙特卡洛方法，在每个体素内生成N个采样点，通过vtkImplicitPolyDataDistance计算这些点到网格的距离，以此得到体素的内外比例。但计算结果出现了明显的偏移异常：

在切片可视化图中，比例从黄色（1=完全外部）到紫色（0=完全内部）渐变，红色点标记该切片处的网格轮廓，两者位置明显不匹配。

我无法确定偏移的根源：是网格与原始图像的原点/空间间距不匹配导致的？还是vtkImplicitPolyDataDistance本身的特性引发的问题？

以下是相关实现代码，其中体素内外比例由voxel_fraction_mc函数中的frac变量计算得到：

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

内容的提问来源于stack exchange，提问作者Daniel Bichou