vtk表面提取参数研究

最新推荐文章于 2024-05-26 00:15:00 发布

pkrobbie

最新推荐文章于 2024-05-26 00:15:00 发布

阅读量2.4k

点赞数 4

分类专栏： C/C++ Image processing 文章标签：医学影像

本文链接：https://blog.csdn.net/pkrobbie/article/details/122029730

版权

C/C++ 同时被 2 个专栏收录

23 篇文章 1 订阅

订阅专栏

Image processing

9 篇文章 0 订阅

订阅专栏

marching cubes是三维图形处理中常见的算法。实际使用中，对一个影像数据做表面提取，然后平滑得到一个组织表面是一个常用功能。下面这段代码是参考3d slicer的分割流程的代码。

    boneExtractor = vtkMarchingCubes::New();
    boneExtractor->ComputeGradientsOff();
    boneExtractor->ComputeNormalsOff(); 
    boneExtractor->SetValue(0, 254);
    boneExtractor->ComputeScalarsOff();

    smooth = vtkWindowedSincPolyDataFilter::New();
    smooth->SetInputData(boneExtractor->GetOutput());
    smooth->SetNumberOfIterations(20); 
    double passBand = 0.01;
    smooth->SetPassBand(passBand);
    smooth->BoundarySmoothingOff();
    smooth->FeatureEdgeSmoothingOff();
    smooth->NonManifoldSmoothingOn();
    smooth->NormalizeCoordinatesOn();

    normalGenerator = vtkPolyDataNormals::New();
    normalGenerator->ConsistencyOn();
    normalGenerator->SplittingOff();

    normalGenerator->SetInputData(smooth->GetOutput());

多数情况下，这个代码都可以很好工作，但是，个别CT数据用这个流程处理后会留下一些奇怪的线。从切面来开，似乎是存在一些无法平滑的面。使用vtk8.1.2时会发生这个问题。

在这里插入图片描述

于是，对同一个CT数据做了一些对比测试，发现一些有趣的现象：

vtk 9.0使用同样的pipeline, 没问题。
vtk 8.1中使用vtkFlyingEdges3d filter替换vtkMarchingCubes也没有问题
vtk8.1.2，如果不设置vtkImageData的Spaing/Origin，也没有这种现象，但是，生成的polydata的坐标不正确。

第一个推论是这种现象可能和原始数据的origin和spacing有关。查了一下数据属性，vtkImageData的origin是[-202.7, -201.5, -1495.75], spacing是[0.787109 0.787109 0.625]。而vtkImageData默认的origin是[0, 0, 0], spacing是[1, 1, 1]。看起来有点像浮点数计算误差为导致MarchingCubes输出的polydata发生的变化影响了平滑。

分析vtk代码发现，vtkMarchingCubes可以通过参数控制判断两个点重合的阈值。在上述代码中增加了这几行进行测试，然后改变SetTolerance的参数进行测试，看vtkMarchingCubes生成的polydata中顶点和面片的数量变化。

    vtkNew<vtkPointLocator> locator;
    boneExtractor->SetLocator(locator);
    locator->SetTolerance(0.01);

同时也记录了9.0的vtkMarchingCubes和8.1的vtkFlyingEdges3d的输出结果特征。结果见下表

算法和参数	顶点数量	多边形数量
8.1 vtkFlyingEdges3D	155112	310224
9.0 vtkMarchingCubes + 默认locator(vtkMergePoints)	155112	310224
8.1 vtkMarchingCubes + 默认locator(vtkMergePoints)	161006	310224
8.1 vtkMarchingCubes + PointLocator tolerance 0.00001	181783	310224
8.1 vtkMarchingCubes + PointLocator tolerance 0.00002	155112	310224
8.1 vtkMarchingCubes + PointLocator tolerance 0.0001	155112	310224
8.1 vtkMarchingCubes + PointLocator tolerance 0.003	155112	310224
8.1 vtkMarchingCubes + PointLocator tolerance 0.004	126826	253325
8.1 vtkMarchingCubes + PointLocator tolerance 0.005	93231	186451
8.1 vtkMarchingCubes + PointLocator tolerance 0.01	93113	186226
8.1 vtkMarchingCubes + PointLocator tolerance 0.02	93113	186226
8.1 vtkMarchingCubes + PointLocator tolerance 0.05	93113	186226

可以看到一下几个规律：

8.1 vtkFlyingEdges3D和9.0 vtkMarchingCubes生成的顶点数量相同，少于8.1 vtkMarchingCubes(默认参数)的顶点数量。三个算法生成的面片数量相同。估计这些多出来的点造成的拓扑结构变化影响了某些位置的平滑效果。
继续增大8.1 vtkMarchingCubes使用的tolerance会发现，实际上真正需要的顶点和面片的数量可以进一步降低。对于这个数据tolerance达到0.01之后，面片和顶点数量就不再增加了。也就是说8.1 vtkFlyingEdges3D和9.0 vtkMarchingCubes仍然会生成一些多余的点和面片。

从Marching Cubes的原理看，所有顶点都位于原来的voxel顶点或者每个边的中点。理论上，tolerance只要小于0.5*min(spacing)就是合理的。不同算法再不同条件上多出来的顶点是由于用浮点数直接判断相等带来的。

回到最初的问题，对于这个测试数据把tolerance生成0.0001或更大就可以保证平滑效果。所以，考虑CT数据特征，最终代码选择·SetTolerance(0.01)。

在这里插入图片描述

如果使用vtkFlyingEdges3D/或者不希望直接给vtkMarchingCubes设置locator，可以考虑在vtkFlyingEdges3D/vtkMarchingCubes之后先经过一个vtkCleanPolyData再进行平滑，猜测可以达到类似的效果。时间有限，没有做实际测试。请大家自行测试。

回头再来分析不同算法之间的区别：

vtk 8.1的vtkMarchingCubes先给每个点都乘以spacing，加上origin，所以，在spacing有效位数很多的时候，更容易出现计算误差带来多余的点。

vtk 9.0的vtkMarchingCubes修改了流程：先不考虑spacing和orgin来计算等值面，最后再对生成的polydata每个点做空间坐标变换。这样就减少了有效位数带来的计算误差，减少了不必要的特征点。

vtkFlyingEdges3D的计算过程类似，算法分成4个步骤，前3步都不使用spacing/origin属性，直到最后一步才使用Origin、Spacing生成计算最终点的坐标。