VTK Python 中 PolyData 与 ImageData 的体素体积分数为何偏移？

结论

代码中有明确的半体素坐标差异：stencil 在 i * voxel_size 对应的位置分类，Monte Carlo 则以 (i + 0.5) * voxel_size 为中心采样。网格切片的 z 坐标、图像显示坐标和采样中心也没有统一。这些差异足以造成所描述的偏移，应先统一坐标，再检查距离符号和网格质量。

这些坐标不一致可以从代码中直接确认，但仅凭提供的代码，无法判断实际数据中所有误差的来源。当前材料也没有可直接检查的附图或 TIFF 数据。

半体素偏移从哪里产生

当前代码将标量写入 GetPointData()。对于这里从零开始的 extent 和默认方向矩阵，vtkImageData 的采样点位置为：

[ \mathbf{p}_{ijk}=\mathbf{o}+(i h,j h,k h) ]

其中，(\mathbf{o}) 是 origin，(h) 是 spacing。SetOrigin(0, 0, 0) 表示索引零对应的采样点位于原点，不能同时把原点理解为围绕该采样点的体素小立方体的最小角点。

原代码采用了两套位置：

因此，inside[i, j, k] 与 centers2[i, j, k] 对应不同的物理位置，三个轴都相差半个 spacing。

还需要区分 VTK 的点和单元。SetDimensions(nx, ny, nz) 创建的是这些数量的采样点；对于三维网格，相邻点之间的单元数量为 (nx - 1, ny - 1, nz - 1)。当前算法可以把每个点标量解释为以该点为中心的体素值，但采样、stencil 和绘图都必须沿用这个约定。

统一输出图像、stencil 和采样中心

沿用原代码中体素中心位于 0.5 * voxel_size 的设计，可以将输出图像的 origin 设置为半个 spacing，再用同一组参数生成 centers2。

下面的代码替换输出网格定义和 pol2stenc 设置部分，后续的 vtkImageStencil 处理可以继续使用：

h = float(voxel_size)
origin = np.array([0.5 * h, 0.5 * h, 0.5 * h])

image = vtk.vtkImageData()
image.SetDimensions(nx, ny, nz)
image.SetOrigin(*origin)
image.SetSpacing(h, h, h)
image.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 1)
image.GetPointData().GetScalars().Fill(1)

x = origin[0] + np.arange(nx) * h
y = origin[1] + np.arange(ny) * h
z = origin[2] + np.arange(nz) * h
xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')
centers2 = np.stack([xx, yy, zz], axis=-1)

pol2stenc = vtk.vtkPolyDataToImageStencil()
pol2stenc.SetInputData(poly)
pol2stenc.SetOutputOrigin(*image.GetOrigin())
pol2stenc.SetOutputSpacing(*image.GetSpacing())
pol2stenc.SetOutputWholeExtent(image.GetExtent())
pol2stenc.Update()

也可以保留输出图像的零 origin，并把 centers2 改成从零开始。两种约定都可以使用，但同一索引必须始终指向同一个物理位置。

输入图像的 origin 需要单独判断。getMeshes() 使用 img.SetOrigin(0, 0, 0)，所以生成的网格已经处于这套输入坐标中。代码中另行计算的 centers1 没有传入该函数，不会改变网格位置。如果原始数据约定首个像素中心位于 0.5 * resol，应在生成网格前设置对应的输入 origin。这个选择取决于原始数据的空间定义，不能为了让图看起来对齐而随意平移网格。

修正距离符号和体积分数定义

对于封闭、朝外定向且内外定义明确的表面，vtkImplicitPolyDataDistance 的负值对应内部，正值对应外部。因此，原函数中的：

sum(di > 0 for di in d) / N

计算的是外部体积分数，注释 fraction INSIDE the mesh 写反了。

外部体积分数的估计为：

[ \hat f_{\mathrm{out}}=\frac{1}{N}\sum_{n=1}^{N}\mathbf{1}[d(\mathbf{p}_n)>0] ]

忽略恰好落在表面上的零测度情形，内部体积分数为 1 - f_out。这与外部体积除以内部体积的比值不同；后者在内部体积为零时没有有限值。

可以将函数写成：

rng = np.random.default_rng(42)

def voxel_fraction_mc(center, voxel_size, implicit, N=2000):
    offsets = rng.uniform(-0.5, 0.5, size=(N, 3)) * voxel_size
    sample_points = np.asarray(center) + offsets
    distances = np.fromiter(
        (implicit.EvaluateFunction(p) for p in sample_points),
        dtype=float,
        count=N,
    )
    return np.mean(distances > 0.0)

frac = np.array([
    voxel_fraction_mc(center, h, implicit)
    for center in centers2[intersection]
])

fraction_out = np.zeros((nx, ny, nz), dtype=float)
fraction_out[outside] = 1.0
fraction_out[intersection] = frac

这里的 N=2000 只是示例采样数，应根据精度和计算成本调整。原代码的 N=50 在真实体积分数接近 0.5 时，估计值的标准误差约为 0.071，足以产生明显的随机起伏。增加采样数可以减小噪声，坐标偏移仍需通过统一空间定义解决。

原代码还使用了 2 * frac，并将外部标记为 2，所以绘制的 centers2_label 范围是 [0, 2]。如果希望颜色直接表示 [0, 1] 的外部体积分数，应绘制上面的 fraction_out。

让切片和绘图使用同一套物理坐标

修正输出 origin 后，第 k 层中心的 z 坐标为：

z_slice = image.GetOrigin()[2] + k * h

原代码的切平面使用 k * h，与原 centers2 的第 k 层中心相差半个体素。

绘图也存在类似问题：pcolormesh(array, shading='nearest') 默认把数组元素中心放在整数索引处，而 points / h 没有扣除图像 origin。直接使用物理坐标绘制，可以减少这种转换错误：

k = nz // 2
z_slice = origin[2] + k * h

plane = vtk.vtkPlane()
plane.SetOrigin(0.0, 0.0, z_slice)
plane.SetNormal(0.0, 0.0, 1.0)

cutter = vtk.vtkCutter()
cutter.SetInputData(poly)
cutter.SetCutFunction(plane)
cutter.Update()

x_edges = origin[0] + (np.arange(nx + 1) - 0.5) * h
y_edges = origin[1] + (np.arange(ny + 1) - 0.5) * h

fig, ax = plt.subplots()
# 数组按 [x, y, z] 存储：绘图横轴为 y，纵轴为 x。
plot = ax.pcolormesh(
    y_edges,
    x_edges,
    fraction_out[:, :, k],
    shading='flat',
    vmin=0.0,
    vmax=1.0,
)

cut_poly = cutter.GetOutput()
if cut_poly.GetNumberOfPoints() > 0:
    cut_points = vtk_to_numpy(cut_poly.GetPoints().GetData())
    ax.scatter(cut_points[:, 1], cut_points[:, 0], c='r', s=1)

ax.set_xlabel('y (mm)')
ax.set_ylabel('x (mm)')
ax.set_aspect('equal')
fig.colorbar(plot, ax=ax, label='Outside volume fraction')

如果继续使用整数索引绘图，网格点应转换为 (points - origin) / h，再按当前数组布局交换用于横轴、纵轴的分量。

体积分数统计的是整个三维体素，而红色轮廓只表示体素中心平面上的截线。即使坐标完全一致，红线也不必穿过所有混合体素，或恰好沿着体积分数为 0.5 的位置。

边界候选体素还需要单独检查

原来的 intersection 来自内外掩码膨胀后的交集，只能作为边界附近的候选区域，既不能证明其中每个体素都与网格相交，也不能保证所有相交体素都会被选中。例如，一个小封闭结构完全落在体素内部，却没有覆盖任何采样中心时，仅凭中心分类可能发现不了它。

对于边长为 h 的立方体，中心到任意位置的最大距离为：

[ r=\frac{\sqrt{3}}{2}h ]

若体素与网格表面相交，其中心到表面的最近距离必然不大于 r。在距离计算有效的前提下，可以用这个条件构造保守候选集合，替换原来的膨胀筛选：

flat_centers = centers2.reshape(-1, 3)
d_center = np.fromiter(
    (implicit.EvaluateFunction(p) for p in flat_centers),
    dtype=float,
    count=len(flat_centers),
).reshape(nx, ny, nz)

half_diagonal = 0.5 * np.sqrt(3.0) * h
intersection = np.abs(d_center) <= half_diagonal

fraction_out = (d_center > 0.0).astype(float)
for i, j, k in np.argwhere(intersection):
    fraction_out[i, j, k] = voxel_fraction_mc(
        centers2[i, j, k], h, implicit
    )

这个条件会包含部分实际不相交的体素，但适合用来排查候选区域遗漏。它需要遍历所有中心，大数据量下应评估成本。Monte Carlo 本身也可能漏采很小的内外体积部分，不能把一次采样得到的 0 或 1 当作严格的无相交证明。

最后，用一个明确位于内部的点和一个明确位于外部的点检查距离符号。若结果异常，再检查网格是否封闭、是否有非流形边或自相交，以及多个材料界面是否具有一致的内外定义。AutoOrientNormalsOn() 和孔洞填充步骤不能单独证明这些条件已经满足。

备注：内容仅供参考。