李丹 彭海欣

摘 要：MarchingCubes算法是医学图像三维重建获取等值面常用的方法，该方法原理简单，但在遍历立方体的时候会遍历许多空立方体，浪费时间。本文提出一种基于顶点状态获取相邻非空体元的方法，根据顶点的状态判断某一面是否含有等值边，含有等值边的那面相邻体元有很大可能含有等值边，减少空立方体的遍历。此外，本文还将算法在基于GPU的基础上进行并行化处理。

关键词：MarchingCubes；医学图像；相邻体元；GPU

中图分类号：TP391.41 文献标识码：A 文章编号：1003-5168（2017）03-0057-04

Abstract： MarchingCubes algorithm is a commonly used method for 3D reconstruction of medical images. The method is simple， but a lot of empty cubes will be traversed when traversing the cube， which wastes time. This paper presented a method based on vertex state to obtain adjacent non - empty body elements and determine whether a side contains the equivalent edge according to the state of the vertex to. It is very likely that the adjacent surface element with the same edge can contain the equivalent edge， which reduces the ergodic of the empty cube. In addition， the algorithm was parallelized based on GPU.

Keywords： MarchingCubes；medical image；adjacent voxel；GPU

近年来，三维重建技术迅速发展，其应用领域十分广泛，其中在医学方面有着重要的学术意义及实践意义。移动立方体MarchingCubes算法[1]（简称MC算法）原理简单，实现容易，逼近精度较高[2]。但是其仍存在不足之处，遍历立方体时会存在空立方体，浪费大量精力，避免对空立方体的遍历，可以提高三维重建速度。由高峰等[3]提出的选取种子体元后根据种子体元衍生出整个器官的等值面，王旭初[4]提出建立接查找子表约束体元搜索路径并根据邻接查找子表特点设计堆栈结构实现搜索算法。王溪[5]提出了表面判断查找表寻找等值面，不需要判断体数据中的每个立方体，但需耗费存储空间存储记录256中每个状态下每个三角形的等值边与面的相对情况，以及计算每条边的情况。本文在此基础上直接通过顶点状态实现判断，而不需要建立查找表及计算等值边。GPU从出现后发展飞速，并且被利用到并行处理中，而GPU强大的计算能力及高效处理效率，可以更快提高MC方法的处理速度。因此，基于GPU的MC方法可以实现更快的重建速度。

1 MarchingCubes算法基本原理

MC算法也就是移动立方体算法，就是利用一个个立方体来获取等值面信息，而等值面则被化为三角面片，MC算法获得的就是三角面片信息。立方体也称体元，是由相邻8个顶点构成，扫描每个体元的8个顶点信息，根据给定的阈值，若包含阈值范围，则体元中必定存在等值面，等值面与体元相交的点成为等值点，相交的边称为等值边。将等值点连接则生成所需的等值面，具体过程如下。

1.1 确定含有等值面的立方体

将三维离散数据场看作立体栅格系统，图像信息读入后，相邻两幅图像构成一层，各自相对应相近的4个点构成一个立方体。由于数据场沿着立方体是线性变化的，所以，若一条边上的2个顶点的灰度值分别大于、小于所规定的阈值，则该边上有且仅有一个与等值面相交的点，即等值点。立方体8个顶点的灰度值分别与规定的阈值比较，以此来判断等值面是否与该立方体相交。根据阈值可以将顶点标记为标记和不标记2种状态：若顶点的阈值大于规定的阈值，则顶点在等值面外部，不标记该顶点；若顶点的阈值小于规定的阈值，则顶点在等值面内部，标记该顶点。每个顶点有2种状态，则1个立方体的8个顶点有256种状态。

1.2 提取等值面

分析立方体的256种状态可以發现，立方体顶点状态具有反对称性和旋转性，即旋转不会改变等值面的拓扑结构，反对称性就是所有顶点状态相反不会改变等值面的连接方式。根据这种情况可以将256种情况化简为15种情况，如图1所示。

可以建立边索引表，将每条边进行编号，边索引中代表这256种状态下每条边的状态，根据8个顶点的状态可以得到一个顶点索引值，根据该索引值可以查看边索引表，得到等值边情况。同时，还可以建立三角形索引记录256种状态下每个三角面片的顶点情况（用顶点编号记录顶点，每3个代表一个三角面片），用得到的边索引值获取三角形索引的三角面片情况。根据8个顶点的状态获取顶点索引值，根据索引值查找边索引，即现在状态下的12条边情况，获取相对应边索引值。再根据2个索引值查找三角索引，获取该立方体的等值面情况。

1.3 计算等值点信息

获取立方体等值面位置后，就可以计算等值点坐标及法向量等信息。等值点计算采取插值计算，常用的有线性插值法和中点选择法。中点选择法是取该等值边的中点，若用p1、p2分别代表2个顶点的坐标，n1、n2代表2个顶点的法向量，则公式可表示为：

MC算法的算法步骤为：①逐步扫描两层数据图像，构建体元；②将体元的所有顶点的灰度值与规定的阈值做比较，建立顶点索引；③若体元含有等值面，利用灰度差分法计算立方体各顶点的梯度；④根据边顶点索引查找在边索引表，获得存在等值点所在的边；⑤根据该相交的边的2个顶点，用线性插值法获取等值点的坐标信息及法向量；⑥根据索引值查找三角面片索引表，确定该立方体内三角面片的顶点连接方式；⑦由三角面片构成等值面；⑧重复以上步骤直到所有的体元遍历结束。

2 基于GPU的改进的MC算法

2.1 改进的MC算法算法原理

一些MC算法中为提高计算速度，通常采用中点法来计算等值点，其计算量比线性插值法少。此外，大多数是建立在减少空立方体的遍历上。在获取图像信息时多数区域信息是无效的，因此会建立许多空立方体，若避免对这些立方体的遍历，基本会减少1/2的时间。本文算法主要思想为确定一个种子体元，以此体元为为起点，遍历相邻体元。

若某一立方体的某一面存在着等值边，又因为模型是相连通的，则该面相邻的立方体的同一面必定存在等值边，除了该立方体正好是边界的情况外，含有等值边的立方体一定存在着等值面。因此，通过判断立方体面上是否有等值边，可判断相邻立方体中一定存在等值面的立方体。逐步遍历这些非空立方体，就可避免对非空立方体的处理。又因为模型是连通的，所以基本不会遗漏主要的特征的三角网格。

在观察立方体的状态时可以发现，若一条边的2个顶点状态相反时必定存在一个等值点，当立方体的一个面存在2个等值点时，必定存在等值边。由此可以发现，当某个面的4个顶点存在2个标记、2个不标记或者3个顶点状态相同、另一顶点状态相反时，该面必定存在等值边，如图2所示。

若存在上述状态即只要一个面的顶点状态不全相同，就可判定该面含有等值边，同时也就可以判定该面相邻立方体为非空立方体，则可以遍历相邻该方向上的立方体。可以对立方体的顶点和边进行编号，以及根据立方体6个面定义6个方向，如图3所示。

若某一面存在等值边，则该方向标记为1，否则为0，并将该方向的体元压入栈中，而下一个遍历的体元就从栈中提取。为了防止同一体元压入栈中或者压入处理过的体元，应对体元设置标记1代表不需要再压入体元，否则为0，在每次压入栈前都需查看该体元状态标记。每个体元在处理过程中不仅需要获取等值面信息，还需获取相邻未处理过的存在等值面的体元信息并将其压入栈中，然后下一个体元不需要再去遍历判断是否是非空体元，直接从栈中读取非空体元，在读取体院同时删除栈中该因素，同时该体元状态更新为1。每遍历一个体元，都可能压入新的非空体元，需要不断从栈中提取体元，直至栈中为空。

该算法的算法步骤为：①逐步扫面图像数据构造体元；②从中间体元开始，按原MC算法找到第一个非空体元；③该体元获取完等值面时，获取6个面的方向标记索引；④若某个方向标记为1，并且该方向的体元的状态0，则将该方向体元坐标压入栈中；⑤下个处理体元从栈中读取，读取完后将该体元信息从栈中删除，并将该体元的状态更改为1；⑥以此读取栈中体元，直到栈为空。

2.2 基于GPU的并行化处理

GPU（Graphics Processing Unit）是图形处理器的简称，这个概念是由NVIDIA公司在发布GeForce256绘图处理芯片时首先提出。GPU与CPU类似，为复杂的数学计算和几何计算而制造，在浮点计算及并行计算等方面有着比CPU高几十倍的性能。因此，对于GPU对于浮点计算的高性能，能对MC算法的插值计算等计算步骤进行更快的计算，并且利用GPU的海量线程进行并行计算，在基于GPU的MC方法上具有更高的效率。将原始的MC方法及本文中改进的MC方法在基于GPU的环境下进行并行。

2.2.1 原MC算法并行化。原始MC算法主要借用GPU的浮点计算能力，其算法流程如下：①并行的获取医学图像数据；②根据获取的医学图像构建体元；③根据顶点状态判断是否含有等值点，若没有返回上一步，继续遍历下一个体元，若有等值面则进行下一步；④获取该体元内的等值面信息；⑤GPU中计算，通过线性插值方法，计算出体元边界与等值而的交点，然后利用中心差分法，求出体元各角点处的法向，再通过线性插值方法，求出三角形各顶点处的法向；⑥遍历体元直至获取完所有体元等值面；⑦利用三角面片信息绘制模型。

2.2.2 基于顶点状态的MC算法并行化。改进的在读入离散数据构造立方体后，首先找出第一个非空立方体后，分别对6个面进行是否有等值边的状态判断，并进行标记，获取完该体元的等值面后根据等值边标记将标记为1的方向的相邻体元压入栈中。计算同样是在GPU中进行的。其算法流程如下：①并行获取两层数据构建立方体；②并行的从立方体中找出第1个非空立方体，从中间开始查找；③对立方体按照正常MC方法进行并行获取等值面信息；④对6个面的状态进行判断，获取相邻含等值面的体元信息并压入栈中；⑤并行从栈中获取待处理的体元，同时将该体元从栈中删除以及将该体元状态改为1；⑥多线程处理体元，重复从步骤③开始处理，在步骤④时，压入栈前先判断体元状态是否为0，若为1，则不压入栈中；⑦并行的重复处理栈中体元直至为空为止。

3 试验分析

在visual studio 2013环境下，用OpenGL实现三维绘制，脊椎骨和头颅阈值为60。由试验结果可以看出，改进的MC算法与原MC算法相比，模型样子基本不变，只是头颅颈椎部分少了一部分边缘琐碎部分，这部分不影响模型本身，而脊椎图像基本没有什么变化（见图4）。

由表1结果可以看出，该算法在原算法的基础上节省了大约50%的时间，虽然可能忽略了一些三角面片，但原模型基本特征没有发生变化。

并行化试验结果见表2，并行化后的算法比原算法快了将近一倍的速度，很明显比串行速度快；而改进后的MC算法时间并行化后比最初时间快了将近3倍，速度有了很大的提高。

4 结语

对于医学图像三维重建，本文在原来的MC算法基础上设置立方体方向来描述相邻立方体，用顶点状态确定含有等值面的立方体方向，将其压入栈中待处理，以此避免空立方體的遍历。本文对其进行实验研究，结果表明该算法比原算法提高了将近一半的速度。同时，还将算法进行并行化处理，可以发现并行化的绘制速度明显提高。基于顶点状态的MC算法是在图像模型具有连通性的情况下，在一个模型中可能存在互相不连通的情况，需要建立多个种子节点。同时，算法需要三维数组记录每个立方体的是否处理情况，由于三角面片数量多，所以三维数组占用空间内存较大。

参考文献：

[1]Willian E Lorense，Harvey E Cline. Marching Cubes：A High Resolution 3d Surface Construction Algorithm[J].Computer Graphics，1987（4）：163-169

[2]王旭初，王赞.基于最近邻Marching Cubes的医学图像三维重建[J].计算机工程与应用，2012（18）：154-158.

[3]高峰，付忠良.基于改进移动立方体的医学图像三维重建算法[J].计算机应用，2013（S1）：201-203，213.

[4]]王旭初，王赞，牛彦敏，等.融合构型查找表与邻接查找子表的改进MC方法[J].重庆大学学报，2012（12）：68-77，83.

[5]王溪，秦新强，党发宁，等.医学图像三维重建的规则移动立方体法[J].2009（4）：477-481.