一种基于三角网格结构的医学虚拟切割算法

刘青，姚莉秀

0 引言

随着计算机图形学的发展，其在虚拟医学导航技术中的应用成为仿真系统研发中的热点。由于外科手术大多需要切割操作，三维模型的虚拟切割交互算法也得到了人们的广泛关注。

被切割物体的建模方法主要有面模型和体模型。面模型仅包含器官的表面形状信息，计算数据量小，但不能表示器官内部结构。体模型能够表达物体内部特征，数据量特别大，很难在实时性要求比较高的场合应用。本文中切割系统对计算精度要求不高，对实时性要求较高，所以采用面模型建模。

在表面网格模型中，虚拟切割算法的仿真主要包括物体与切割工具的碰撞检测，物体网格的切割算法和切割后网格间的分离3部分。文中采用OBB包围盒方法进行碰撞检测计算，切割部分采用基于顶点移动的算法实现，并对切割后的封闭区间采用广度遍历算法进行网格分离。

1 碰撞的检测

碰撞检测就是在虚拟场景下，判断物体之间是否发生碰撞，也即物体间的求交检测。目前，空间几何模型间的碰撞检测算法大致可分为两类，空间分解法和层次包围盒法[1]。

空间分割法是将整个空间划分成体积相等的小单元格，对占据了相同单元格或相邻单元格的几何模型进行相交测试。典型的空间分解法有八叉树（Octree）法和二叉空间剖分（Binary Space Partition）法。

层次包围盒法的基本思想是用体积稍大且特性简单的几何体（包围盒）来近似代替复杂几何对象进行相交测试，并通过构造树状层次结构逐渐逼近对象的几何特性。进行重叠测试时只有在包围盒重叠时才对其重叠部分进行进一步的相交测试，大大减少了参与相交测试的包围盒的数目，提高了检测的效率。本文中主要介绍层次包围盒法。

1.1 包围盒

目前用于碰撞检测的包围盒主要有包围球，沿坐标轴的包围盒（AABB axis-aligned bounding boxes）、方向包围盒（OBB oriented bounding box）等，如图1所示。

图1 包围盒示意图

给定物体的包围球是指包含该对象的最小球体。计算包围球的方法非常简单。包围球间的相交测试也比较简单，若两包围球球心之间的距离小于两半径之和则两包围球相交，否则不相交。物体包围球的紧致性比较差。

给定物体的AABB包围盒被定义为包含该对象且各边平行于坐标轴的最小的六面体。AABB包围盒的相交测试也比较简单，两个AABB相交当且仅当他们在3个坐标轴上的投影区间均重叠，只要存在一个方向上的投影布重叠它们就不相交。

给定物体的OBB包围盒被定义为包含该对象且相对于坐标轴方向任意的最小的正六面体。设物体由三角形网格构成，构造物体的OBB包围盒需要根据网格的顶点数据计算均值和协方差矩阵，然后以协方差矩阵的特征向量作为包围盒局部坐标的轴方向，通过计算物体在轴方向上的最大最小值确定包围盒的大小。

OBB包围盒常用的一种相交测试方法叫做“轴投影”。该方法将包围盒投影在空间的某个轴上，如果两个包围盒在该轴上的投影不相交，这两个包围盒是分离的，这条轴称作这两个包围盒的“分离轴”，否则这两个包围盒的关系不确定，需要继续判断。

两个空间多面体是分离的当且仅当存在一条分离轴，该分离轴垂直于两个多面体中的一个面或者同时垂直于两个多面体中的一条边。由于每个空间包围盒有3个不同的面方向，3个不同的边方向，所以共需要15次分离轴判断。如果两个包围盒是分离的，则上述的15个轴方向中必有一个方向为分离轴，如果重叠则不存在分离轴。

1.2 包围盒树

采用层次包围盒法的基础是构造包围盒树。通常构造包围盒树的方法主要有两种，自顶向下法和自底向上法。这里主要介绍自顶向下法。

自顶向下法以全集S作为根节点，利用全局信息递归地对节点进行划分以形成子节点，直到达到叶节点，叶节点中仅包含基本几何元素。该方法的核心是如何把一个集合划分为若干不相交的子集。对集合的划分通常采用分裂平面法，即选定一个空间平面，根据集合元素相对于平面的位置进行划分。一个基本元素或者属于平面的左半空间，或者属于平面的右半空间，再或者与平面相交，跨越这两个半空间。对于前两种情况可把元素分到相应的子集中，而后一种情况需要依据元素中心点相对平面的位置来确定，但若中心点恰好位于平面上，则将元素分给所含元素较少的子集。图2中三角面片1属于S1子集，三角面片2和3属于S2子集。

图2 包围盒元素子集划分

分裂平面法的关键是分裂面的选取。分裂轴的确定通常选择最长轴法，即选择分裂轴方向为包围盒长度最长的方向。分裂点的选择有平均值方法和中值方法，平均值法即选择集合中元素在分裂轴上投影的平均值，该法通常可使子节点包围盒体积更小，碰撞检测速度提高地更多。中值方法是计算几何元素在分裂轴上投影的中值，该法简单快速，且可以得到两个大小相等的子集，从而得到一棵基本平衡的包围盒树，尤其当集合中的基本元素大小相等不多时，该法更有效。图3为二维情况下OBB层次包围盒树的构造过程。

图3 OBB层次包围盒树构造过程

2 基于顶点移动的表面网格切割算法

针对表面网格结构目前常用的切割算法，主要有面片剖切法和顶点移动法。面片剖切法依照切割路径对三角面片进行剖分，该法可以保存原始网格的拓扑结构，但容易产生畸形三角面片而对后续操作产生巨大的影响。

顶点移动方法[2]的基本思想是在切割过程中，计算出切割工具与每个三角面片相交的交点，然后计算需要移动的顶点位置，移动顶点后沿切割边分离切割网格，切割过程如图4所示。本文中采用移动顶点移动的方法对三角网格进行剖分。

图4 基于顶点移动的表面网格切割过程

2.1 切割路径

通常切割工具的仿真有三种。一种将切割工具视为与网格表面的一系列碰撞点，一种将切割工具视为简单的平面或直线，还有一种由多个面片构成切割工具进行渲染，但在具体计算中使用其中轴线进行计算[3]。

本文中采用第三种仿真方法，如图5所示。切割刀片采用单一的三角面片，刀柄由三角网格构成。在碰撞检测中仅使用切割刀片刀刃所在的直线进行相交计算。切割工具沿刀片所在的平面移动，切割过程中切割工具可以进行旋转。

图5 本文中切割工具仿真方法

切割路径定义为切割工具移动过程中与表面网格的一系列交点。在实现过程中，计算刀刃进入三角面片时与入边的交点和移出三角面片时与出边的交点，并依据两个交点与顶点之间的关系进行顶点移动。

2.2 移动顶点[2]

在计算得到入边和出边交点后进行顶点移动。顶点移动遵循的原则是就近原则。即根据交点坐标位置判断与交点相邻的两个顶点到该交点的距离大小，将距离交点较小的顶点移动到交点的位置处。交点移动过程中有两种特殊情况需要注意：

情况一：两出入交点距离三角形某个顶点非常接近，且该顶点到两个交点的距离大小相等。在这种情况下将顶点移动到交点所在边的边长较大的交点处，另一个交点也移动到该处。

情况二：两出入交点分别位于出入边的中点处，则将除两边公共顶点之外的两个顶点移动到相应得交点位置。

图6 顶点移动时的特殊情况

可见，经过顶点移动后，不再出现狭长三角形或与原模型中的三角形不同数量级的小三角形，有利于在网格上进行后续操作，如继续切割或变形计算等。

2.3 切割路径分离

顶点移动后，通过将位于切割路径上的顶点分裂为两个来分离网格，显示切割路径。两个分裂点间的连线垂直于切割平面，且二者与切割平面间的距离相等，如图7所示。

图7 切割路径分离方法

3 切割搜索

切割过程结束后如果被切下的网格部分与原始网格不连通则需要将二者分离，这就需要采用网格搜索算法。网格搜索算法中包含网格拓扑结构表示和在拓扑结构下的搜索算法。

3.1 AIF数据结构

AIF（adjacency and incidence framework）数据结构[5]通过构建多边形间的邻接入射关系来表示多边形网格的拓扑结构。该法可以高效地实现创建和搜索，且在切割过程中可以方便地维护。

AIF通过描述网格结构的基本元素来表达网格的拓扑结构。当一种元素x是另一种元素y的边界之一，且其几何维数低于y，称x邻接于y，表示为x＜y，或y入射于x，表示为y＞x。

当网格由三角面片组成时，基本元素包括点、线和面。利用网格结构的邻接和入射关系，AIF数据结构构造包含点集、边集和面集及表征它们之间关系的数据集合。各集合中每个元素都包含与之有邻接和入射关系的其他元素的子集合，即每个顶点数据中包含所有通过该顶点的边数据的集合，每条边数据中包含所有经过该边的所有面的数据集合，每个面数据中包含所有构成该面的所有边的数据集合。通过AIF数据结构，可以利用一种元素得到所有与之有邻接关系的数据元素。数据结构实现如下：

Class Vertex{

float x，y，z；//点的空间坐标

vector＜Edge*＞edgeSet；//经过该点的所有边的集合

}

Class Edge{

Vertex* p1，*p2；//组成该边的两个顶点的指针

vector＜Face*＞faceSet；//经过该边的所有面的集合

}

Class Face{

vector＜Edge*＞edgeSet；//构成该面的所有边的集合

}

3.2 基于AIF结构的搜索算法

利用AIF数据结构，可以从一种网格基本元素得到与之相邻的其他网格基本元素。可采用图的广度优先算法来进行网格数据搜索。

图的广度优先遍历算法的基本思想是首先从图中某个顶点V0出发，访问该顶点后访问各个未曾访问过的邻接点W1，W2，…，Wk，然后依次从W1，W2，…，Wk出发访问各自未被访问的邻接点。如此反复直到所有点都被访问过为止。

为了实现上述算法，本文在AIF的点和面结构中加入了reach属性，表示点和面是否已经被访问过。同时还引入队列结构辅助算法实现。首先将初始顶点的reach属性设置为true并将该点加入队列中，网格结构中其他元素的reach属性设置为false。若队列非空，取出队列中的首元素，通过AIF数据结构得到所有通过该点的网格面，若网格面的reach属性为false，设置其属性为true并利用AIF数据结构得到该网格面的所有顶点，设置reach属性为false的网格点为true并将这些网格点加入队列尾部，如此反复直到队列为空。

网格搜索算法比较简单，但搜索过程中初始点的选择非常重要。在网格切割过程中，通常被切下的网格部分相比原来的网格结构是非常小的，即包含数量较少的点、线和面数据。此时如果选取被切下网格中的顶点作为搜索起点，则搜索过程中需要遍历的网格节点较少，计算量较小，而如果选择原始网格上的网格顶点作为搜索起始点，则搜索算法的计算量很大。

图8中原始网格被剖分为两个网格G1和G2，G2为一封闭区域且不与G1连通。切割的起始点分别为S1和S2。可见若从G1中的顶点（如S1）开始进行网格搜索，最终将遍历G1，并将其从两个网格结构中分离，但G1包含的数据量远远大于G2，从切割的实时性考虑应从G2中的顶点（如S2）开始搜索并将G2分离出来。

图8 网格分离时起始点的选择

通常选择G2中的S2作为起始点比较方便。为了选择起始搜索点S2，在切割开始时保存切割起始点S1和S2，并在每次切割后计算包含S1的切割路径path1的长度和包含起始切割点S2的切割路径path2的长度。由上图可见，由于闭合的切割路径path2在path1内部，所以path2的长度小于path1，所以切割结束后选择切割路径长度小的起始切割点作为网格搜索的起始点。

4 实验结果

实验中将上述算法应用到了CT颅面数据中。首先采用marching cube算法提取等值面网格，网格结构由三角面片组成，然后移动切割工具并进行碰撞检测和网格切割，切割结束后找到被切下网格上的顶点进行网格分离。

在切割交互中，移动切割工具并进行其与物体网格的碰撞检测。碰撞检测采用1中所述的OBB包围盒树算法。虽然OBB包围盒在物体拓扑结构发生改变后需要更新，且更新算法比较复杂，但考虑到首次碰撞后切割工具连续移动，所以仅在初次碰撞后采用碰撞检测得到的三角面片与切割三角面片进行相交计算，后续的相交计算则依据切割工具的移动方向和AIF数据结构进行。因切割工具沿切割刀片所在的平面移动，所以在计算得到当前三角面片的出边信息后由AIF结构找到下一个三角面片，如此反复计算直到得到某个三角面片与当前切割三角面片相交，以该三角面片为当前碰撞面并进行相交计算。

图9为本文算法在医学颅面整复手术中的切割应用。图（a）为由marching cube算法提取的三维等值面网格结构，其包含两层等值面结构。外层为颅面表皮，内层为骨组织表面。图（b）和图（c）为图（a）网格数据的填充渲染模式。图（d）为颅面表皮被切割后的切割路径显示。图（e）为分离被切下网格结构与原始网格结构后显示出内层骨组织。

图9 实验结果

5 结语

文中针对医学导航系统中的虚拟切割算法进行了系统仿真。对基于表面重建的三角网格数据建立AIF数据结构以维护其空间拓扑结构，然后在移动切割工具的过程中采用OBB包围盒树的方法进行碰撞检测，检测到碰撞发生后采用基于顶点移动的网格切割方法进行网格切割，切割结束后若切割路径为封闭曲线，则对由切割路径包围的网格结构进行搜索遍历并与原网格分离。实验结果表明文中的系统算法能够较好的实现切割要求。通常人体组织为软组织，在切割过程中会有形变的产生，本文中没有涉及变形计算，后续研究将在三角网格结构的基础上考虑加入切割过程中的网格变形计算，使切割仿真更加真实。

[1]潘振宽，崔树娟，张继萍等.基于层次包围盒的碰撞检测方法[J].青岛大学学报:自然科学版，2005，18:71 76.

[2]王钰，于素平.针对面模型的顶点移动切割算法研究[J].系统仿真技术，2008，4（2）:117 121.

[3]Bruyns C，Senger S，Menon A，Montgomery K，Wildermuth S，Boyle R.A Survey of Interactive Mesh Cutting Techniques and A New Method for Implementing Generalized Interactive Mesh Cutting Using Virtual Tools[J].The Journal of Visualization and Computer Animation，2002，13（1）:21-42.

[4]Yi-Je L，Hu J，Chu-Yin C.et al.Soft Tissue Deformation and Cutting Simulation for the Multimodal Surgery Training[C]// Proceedings of the 19th IEEE Symposium on Computer-Based Medical Systems.New York: IEEE Press，2006:635-640.

[5]黄洁，杨杰.基于 AIF的三角形网格切割方法[J].上海交通大学学报，2008，42（004）: 564-568.