如何用SimpleITK/VTK/PyVista提取MRI NRRD的皮肤表面并处理PolyData？

核心需求 #

从脑部MRI的3D重建PolyData中仅保留最大的皮肤表面，需解决孔洞填充及最大表面提取的问题。

已尝试方案 #

• 方案1：使用SimpleITK生成掩码，计划通过形态学闭运算填充孔洞后提取边界，代码片段：

filter = sitk.BinaryThresholdImageFilter()
filter.SetLowerThreshold()
filter.SetUpperThreshold()
mask = filter.Execute(img) 

• 方案2：通过PyVista处理DICOM转换的NRRD文件，阈值处理后调用extract_surface()和fill_holes()，但fill_holes()未改善结果。

可行解决步骤 #

方法一：基于PyVista的完整流程 #

1. 表面填充+最大组件提取

先提取表面，再针对性填充孔洞，最后筛选最大连通表面：

import pyvista as pv

# 读取NRRD文件
volume = pv.read("brain_mri.nrrd")
# 阈值处理：根据皮肤实际HU值调整范围（示例为-1000到0，需匹配你的数据）
skin_volume = volume.threshold([-1000, 0])
# 提取初始表面
skin_surface = skin_volume.extract_surface()
# 填充孔洞：设置足够大的hole_size（比如1000）覆盖皮肤区域的孔洞
skin_surface_filled = skin_surface.fill_holes(hole_size=1000)
# 提取最大连通组件（仅保留皮肤表面）
largest_skin = skin_surface_filled.extract_largest_component()
# 可视化验证
largest_skin.plot()

注意：hole_size参数需根据数据调整，太小无法填充大孔洞，太大可能误填充有效区域。

2. 体数据预处理优化

如果表面填充效果不佳，先对阈值后的体数据做形态学闭运算，再提取表面：

# 对体数据做形态学闭运算，填充内部小空洞
structuring_element = pv.Cube()  # 也可选用Sphere，按需调整
closed_volume = skin_volume.morphological_operation("close", kernel=structuring_element, iterations=2)
# 提取表面并筛选最大组件
closed_surface = closed_volume.extract_surface()
final_skin = closed_surface.extract_largest_component()

方法二：基于SimpleITK的流程 #

补全掩码处理与边界提取的完整步骤：

import SimpleITK as sitk
import numpy as np
import pyvista as pv

# 读取图像（假设img已加载）
filter = sitk.BinaryThresholdImageFilter()
filter.SetLowerThreshold(-1000)  # 替换为皮肤实际HU下限
filter.SetUpperThreshold(0)       # 替换为皮肤实际HU上限
mask = filter.Execute(img)

# 形态学闭运算填充孔洞：调整KernelRadius适配孔洞大小
closing_filter = sitk.BinaryMorphologicalClosingImageFilter()
closing_filter.SetKernelRadius(2)  # 对应5x5x5结构元（半径2，边长=2*2+1）
closed_mask = closing_filter.Execute(mask)

# 提取边界图像
contour_filter = sitk.BinaryContourImageFilter()
contour_filter.SetFullyConnected(True)
contour_image = contour_filter.Execute(closed_mask)

# 转换为PyVista PolyData并提取最大表面
contour_array = sitk.GetArrayFromImage(contour_image)
# 转换数组顺序：SimpleITK为(z,y,x)，PyVista需要(x,y,z)
contour_array = np.transpose(contour_array, (2,1,0))
# 创建UniformGrid
grid = pv.UniformGrid()
grid.dimensions = contour_array.shape
grid.spacing = img.GetSpacing()
grid.origin = img.GetOrigin()
grid.point_data["values"] = contour_array.flatten()
# 提取表面并筛选最大组件
skin_surface = grid.threshold(1).extract_surface()
final_skin = skin_surface.extract_largest_component()

关键注意事项 #

• 阈值校准：皮肤HU值通常在-1000到0之间，但需根据你的MRI数据直方图调整，确保仅选中皮肤区域。
• 参数调试：孔洞填充、形态学运算的参数需反复调试，匹配数据中的孔洞大小。
• 连通组件提取：这一步是确保仅保留皮肤表面的核心，因为脑部其他组织的表面远小于皮肤。

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