如何将VTK生成的颅骨3D表面模型转换为实心模型？

问题描述 #

我使用VTK官方示例代码创建颅骨3D模型，模型运行正常，但仅生成空心的表面结构，需要将内部填充，得到完整的实心3D模型。原代码如下：

#!/usr/bin/env python

# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.vtkCommonCore import (
    VTK_VERSION_NUMBER,
    vtkVersion
)
from vtkmodules.vtkCommonDataModel import vtkMergePoints
from vtkmodules.vtkFiltersCore import (
    vtkFlyingEdges3D,
    vtkMarchingCubes
)
from vtkmodules.vtkFiltersModeling import vtkOutlineFilter
from vtkmodules.vtkIOImage import vtkMetaImageReader
from vtkmodules.vtkRenderingCore import (
    vtkActor,
    vtkPolyDataMapper,
    vtkRenderWindow,
    vtkRenderWindowInteractor,
    vtkRenderer
)


def main():
    # vtkFlyingEdges3D was introduced in VTK >= 8.2
    use_flying_edges = vtk_version_ok(8, 2, 0)

    file_name = get_program_parameters()

    colors = vtkNamedColors()

    # Create the RenderWindow, Renderer and Interactor.
    ren = vtkRenderer()
    ren_win = vtkRenderWindow()
    ren_win.AddRenderer(ren)
    iren = vtkRenderWindowInteractor()
    iren.SetRenderWindow(ren_win)

    # Create the pipeline.
    reader = vtkMetaImageReader()
    reader.SetFileName(file_name)
    reader.Update()

    locator = vtkMergePoints()
    locator.SetDivisions(64, 64, 92)
    locator.SetNumberOfPointsPerBucket(2)
    locator.AutomaticOff()

    if use_flying_edges:
        try:
            using_marching_cubes = False
            iso = vtkFlyingEdges3D()
        except AttributeError:
            using_marching_cubes = True
            iso = vtkMarchingCubes()
    else:
        using_marching_cubes = True
        iso = vtkMarchingCubes()
    iso.SetInputConnection(reader.GetOutputPort())
    iso.ComputeGradientsOn()
    iso.ComputeScalarsOff()
    iso.SetValue(0, 1150)
    if using_marching_cubes:
        iso.SetLocator(locator)

    iso_mapper = vtkPolyDataMapper()
    iso_mapper.SetInputConnection(iso.GetOutputPort())
    iso_mapper.ScalarVisibilityOff()

    iso_actor = vtkActor()
    iso_actor.SetMapper(iso_mapper)
    iso_actor.GetProperty().SetColor(colors.GetColor3d('Ivory'))

    outline = vtkOutlineFilter()
    outline.SetInputConnection(reader.GetOutputPort())

    outline_mapper = vtkPolyDataMapper()
    outline_mapper.SetInputConnection(outline.GetOutputPort())

    outline_actor = vtkActor()
    outline_actor.SetMapper(outline_mapper)

    # Add the actors to the renderer, set the background and size.
    #
    ren.AddActor(outline_actor)
    ren.AddActor(iso_actor)
    ren.SetBackground(colors.GetColor3d('SlateGray'))
    ren.GetActiveCamera().SetFocalPoint(0, 0, 0)
    ren.GetActiveCamera().SetPosition(0, -1, 0)
    ren.GetActiveCamera().SetViewUp(0, 0, -1)
    ren.ResetCamera()
    ren.GetActiveCamera().Dolly(1.5)
    ren.ResetCameraClippingRange()

    ren_win.SetSize(640, 480)
    ren_win.SetWindowName('HeadBone')

    ren_win.Render()
    iren.Start()


def get_program_parameters():
    import argparse
    description = 'Marching cubes surface of human bone.'
    epilogue = '''
    '''
    parser = argparse.ArgumentParser(description=description, epilog=epilogue,
                                     formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument('filename', help='FullHead.mhd.')
    args = parser.parse_args()
    return args.filename


def vtk_version_ok(major, minor, build):
    """
    Check the VTK version.

    :param major: Major version.
    :param minor: Minor version.
    :param build: Build version.
    :return: True if the requested VTK version is greater or equal to the actual VTK version.
    """
    needed_version = 10000000000 * int(major) + 100000000 * int(minor) + int(build)
    try:
        vtk_version_number = VTK_VERSION_NUMBER
    except AttributeError:  # as error:
        ver = vtkVersion()
        vtk_version_number = 10000000000 * ver.GetVTKMajorVersion() + 100000000 * ver.GetVTKMinorVersion() \
                             + ver.GetVTKBuildVersion()
    if vtk_version_number >= needed_version:
        return True
    else:
        return False


if __name__ == '__main__':
    main()

解决方案 #

Marching Cubes/Flying Edges这类等值面算法本质是提取表面网格（PolyData），不会生成实心体。要得到实心模型，可采用以下两种方案：

方案1：基于体数据直接生成实心体（适合3D打印/有限元分析）

直接从原始体数据中过滤出骨组织区域，转为实心体网格：

1. 导入额外的VTK模块：vtkThreshold（过滤体素）、vtkDataSetMapper（渲染体网格）
2. 替换原代码中等值面提取的部分，改为阈值过滤+体网格转换

修改后的核心代码：

# 替换原等值面提取的代码段
# --------------------------
# 使用vtkThreshold过滤出骨组织体素
threshold = vtkThreshold()
threshold.SetInputConnection(reader.GetOutputPort())
threshold.ThresholdByUpper(1150)  # 保留值>=1150的体素（骨组织）
threshold.Update()

# 将过滤后的体数据转为UnstructuredGrid（实心体）
grid_mapper = vtkDataSetMapper()
grid_mapper.SetInputConnection(threshold.GetOutputPort())
grid_mapper.ScalarVisibilityOff()

solid_actor = vtkActor()
solid_actor.SetMapper(grid_mapper)
solid_actor.GetProperty().SetColor(colors.GetColor3d('Ivory'))
# --------------------------

# 替换原添加iso_actor的代码，改为添加solid_actor
ren.AddActor(outline_actor)
ren.AddActor(solid_actor)

方案2：闭合表面后生成实心四面体网格（保留表面精度）

如果需要保留等值面的高精度，先确保表面闭合，再生成实心体：

1. 导入额外模块：vtkFillHolesFilter（补洞）、vtkDelaunay3D（生成四面体网格）
2. 对原等值面进行补洞处理，再生成实心体

核心代码示例：

# 在原iso提取代码后添加以下步骤
# --------------------------
# 补全表面孔洞，确保网格闭合
fill_holes = vtkFillHolesFilter()
fill_holes.SetInputConnection(iso.GetOutputPort())
fill_holes.SetHoleSize(1000.0)  # 根据实际情况调整孔洞大小阈值

# 基于闭合表面生成实心四面体网格
delaunay = vtkDelaunay3D()
delaunay.SetInputConnection(fill_holes.GetOutputPort())
delaunay.SetTolerance(0.01)

# 渲染实心体
solid_mapper = vtkDataSetMapper()
solid_mapper.SetInputConnection(delaunay.GetOutputPort())
solid_mapper.ScalarVisibilityOff()

solid_actor = vtkActor()
solid_actor.SetMapper(solid_mapper)
solid_actor.GetProperty().SetColor(colors.GetColor3d('Ivory'))
# --------------------------

# 替换原添加iso_actor的代码，改为添加solid_actor
ren.AddActor(outline_actor)
ren.AddActor(solid_actor)

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