一種基于三角網格結構的醫學虛擬切割算法

劉青，姚莉秀

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.