一種基于三角網(wǎng)格結(jié)構(gòu)的醫(yī)學(xué)虛擬切割算法

劉青，姚莉秀

0 引言

隨著計(jì)算機(jī)圖形學(xué)的發(fā)展，其在虛擬醫(yī)學(xué)導(dǎo)航技術(shù)中的應(yīng)用成為仿真系統(tǒng)研發(fā)中的熱點(diǎn)。由于外科手術(shù)大多需要切割操作，三維模型的虛擬切割交互算法也得到了人們的廣泛關(guān)注。

被切割物體的建模方法主要有面模型和體模型。面模型僅包含器官的表面形狀信息，計(jì)算數(shù)據(jù)量小，但不能表示器官內(nèi)部結(jié)構(gòu)。體模型能夠表達(dá)物體內(nèi)部特征，數(shù)據(jù)量特別大，很難在實(shí)時(shí)性要求比較高的場(chǎng)合應(yīng)用。本文中切割系統(tǒng)對(duì)計(jì)算精度要求不高，對(duì)實(shí)時(shí)性要求較高，所以采用面模型建模。

在表面網(wǎng)格模型中，虛擬切割算法的仿真主要包括物體與切割工具的碰撞檢測(cè)，物體網(wǎng)格的切割算法和切割后網(wǎng)格間的分離3部分。文中采用OBB包圍盒方法進(jìn)行碰撞檢測(cè)計(jì)算，切割部分采用基于頂點(diǎn)移動(dòng)的算法實(shí)現(xiàn)，并對(duì)切割后的封閉區(qū)間采用廣度遍歷算法進(jìn)行網(wǎng)格分離。

1 碰撞的檢測(cè)

碰撞檢測(cè)就是在虛擬場(chǎng)景下，判斷物體之間是否發(fā)生碰撞，也即物體間的求交檢測(cè)。目前，空間幾何模型間的碰撞檢測(cè)算法大致可分為兩類，空間分解法和層次包圍盒法[1]。

空間分割法是將整個(gè)空間劃分成體積相等的小單元格，對(duì)占據(jù)了相同單元格或相鄰單元格的幾何模型進(jìn)行相交測(cè)試。典型的空間分解法有八叉樹（Octree）法和二叉空間剖分（Binary Space Partition）法。

層次包圍盒法的基本思想是用體積稍大且特性簡(jiǎn)單的幾何體（包圍盒）來近似代替復(fù)雜幾何對(duì)象進(jìn)行相交測(cè)試，并通過構(gòu)造樹狀層次結(jié)構(gòu)逐漸逼近對(duì)象的幾何特性。進(jìn)行重疊測(cè)試時(shí)只有在包圍盒重疊時(shí)才對(duì)其重疊部分進(jìn)行進(jìn)一步的相交測(cè)試，大大減少了參與相交測(cè)試的包圍盒的數(shù)目，提高了檢測(cè)的效率。本文中主要介紹層次包圍盒法。

1.1 包圍盒

目前用于碰撞檢測(cè)的包圍盒主要有包圍球，沿坐標(biāo)軸的包圍盒（AABB axis-aligned bounding boxes）、方向包圍盒（OBB oriented bounding box）等，如圖1所示。

圖1 包圍盒示意圖

給定物體的包圍球是指包含該對(duì)象的最小球體。計(jì)算包圍球的方法非常簡(jiǎn)單。包圍球間的相交測(cè)試也比較簡(jiǎn)單，若兩包圍球球心之間的距離小于兩半徑之和則兩包圍球相交，否則不相交。物體包圍球的緊致性比較差。

給定物體的AABB包圍盒被定義為包含該對(duì)象且各邊平行于坐標(biāo)軸的最小的六面體。AABB包圍盒的相交測(cè)試也比較簡(jiǎn)單，兩個(gè)AABB相交當(dāng)且僅當(dāng)他們?cè)?個(gè)坐標(biāo)軸上的投影區(qū)間均重疊，只要存在一個(gè)方向上的投影布重疊它們就不相交。

給定物體的OBB包圍盒被定義為包含該對(duì)象且相對(duì)于坐標(biāo)軸方向任意的最小的正六面體。設(shè)物體由三角形網(wǎng)格構(gòu)成，構(gòu)造物體的OBB包圍盒需要根據(jù)網(wǎng)格的頂點(diǎn)數(shù)據(jù)計(jì)算均值和協(xié)方差矩陣，然后以協(xié)方差矩陣的特征向量作為包圍盒局部坐標(biāo)的軸方向，通過計(jì)算物體在軸方向上的最大最小值確定包圍盒的大小。

OBB包圍盒常用的一種相交測(cè)試方法叫做“軸投影”。該方法將包圍盒投影在空間的某個(gè)軸上，如果兩個(gè)包圍盒在該軸上的投影不相交，這兩個(gè)包圍盒是分離的，這條軸稱作這兩個(gè)包圍盒的“分離軸”，否則這兩個(gè)包圍盒的關(guān)系不確定，需要繼續(xù)判斷。

兩個(gè)空間多面體是分離的當(dāng)且僅當(dāng)存在一條分離軸，該分離軸垂直于兩個(gè)多面體中的一個(gè)面或者同時(shí)垂直于兩個(gè)多面體中的一條邊。由于每個(gè)空間包圍盒有3個(gè)不同的面方向，3個(gè)不同的邊方向，所以共需要15次分離軸判斷。如果兩個(gè)包圍盒是分離的，則上述的15個(gè)軸方向中必有一個(gè)方向?yàn)榉蛛x軸，如果重疊則不存在分離軸。

1.2 包圍盒樹

采用層次包圍盒法的基礎(chǔ)是構(gòu)造包圍盒樹。通常構(gòu)造包圍盒樹的方法主要有兩種，自頂向下法和自底向上法。這里主要介紹自頂向下法。

自頂向下法以全集S作為根節(jié)點(diǎn)，利用全局信息遞歸地對(duì)節(jié)點(diǎn)進(jìn)行劃分以形成子節(jié)點(diǎn)，直到達(dá)到葉節(jié)點(diǎn)，葉節(jié)點(diǎn)中僅包含基本幾何元素。該方法的核心是如何把一個(gè)集合劃分為若干不相交的子集。對(duì)集合的劃分通常采用分裂平面法，即選定一個(gè)空間平面，根據(jù)集合元素相對(duì)于平面的位置進(jìn)行劃分。一個(gè)基本元素或者屬于平面的左半空間，或者屬于平面的右半空間，再或者與平面相交，跨越這兩個(gè)半空間。對(duì)于前兩種情況可把元素分到相應(yīng)的子集中，而后一種情況需要依據(jù)元素中心點(diǎn)相對(duì)平面的位置來確定，但若中心點(diǎn)恰好位于平面上，則將元素分給所含元素較少的子集。圖2中三角面片1屬于S1子集，三角面片2和3屬于S2子集。

圖2 包圍盒元素子集劃分

分裂平面法的關(guān)鍵是分裂面的選取。分裂軸的確定通常選擇最長(zhǎng)軸法，即選擇分裂軸方向?yàn)榘鼑虚L(zhǎng)度最長(zhǎng)的方向。分裂點(diǎn)的選擇有平均值方法和中值方法，平均值法即選擇集合中元素在分裂軸上投影的平均值，該法通?？墒棺庸?jié)點(diǎn)包圍盒體積更小，碰撞檢測(cè)速度提高地更多。中值方法是計(jì)算幾何元素在分裂軸上投影的中值，該法簡(jiǎn)單快速，且可以得到兩個(gè)大小相等的子集，從而得到一棵基本平衡的包圍盒樹，尤其當(dāng)集合中的基本元素大小相等不多時(shí)，該法更有效。圖3為二維情況下OBB層次包圍盒樹的構(gòu)造過程。

圖3 OBB層次包圍盒樹構(gòu)造過程

2 基于頂點(diǎn)移動(dòng)的表面網(wǎng)格切割算法

針對(duì)表面網(wǎng)格結(jié)構(gòu)目前常用的切割算法，主要有面片剖切法和頂點(diǎn)移動(dòng)法。面片剖切法依照切割路徑對(duì)三角面片進(jìn)行剖分，該法可以保存原始網(wǎng)格的拓?fù)浣Y(jié)構(gòu)，但容易產(chǎn)生畸形三角面片而對(duì)后續(xù)操作產(chǎn)生巨大的影響。

頂點(diǎn)移動(dòng)方法[2]的基本思想是在切割過程中，計(jì)算出切割工具與每個(gè)三角面片相交的交點(diǎn)，然后計(jì)算需要移動(dòng)的頂點(diǎn)位置，移動(dòng)頂點(diǎn)后沿切割邊分離切割網(wǎng)格，切割過程如圖4所示。本文中采用移動(dòng)頂點(diǎn)移動(dòng)的方法對(duì)三角網(wǎng)格進(jìn)行剖分。

圖4 基于頂點(diǎn)移動(dòng)的表面網(wǎng)格切割過程

2.1 切割路徑

通常切割工具的仿真有三種。一種將切割工具視為與網(wǎng)格表面的一系列碰撞點(diǎn)，一種將切割工具視為簡(jiǎn)單的平面或直線，還有一種由多個(gè)面片構(gòu)成切割工具進(jìn)行渲染，但在具體計(jì)算中使用其中軸線進(jìn)行計(jì)算[3]。

本文中采用第三種仿真方法，如圖5所示。切割刀片采用單一的三角面片，刀柄由三角網(wǎng)格構(gòu)成。在碰撞檢測(cè)中僅使用切割刀片刀刃所在的直線進(jìn)行相交計(jì)算。切割工具沿刀片所在的平面移動(dòng)，切割過程中切割工具可以進(jìn)行旋轉(zhuǎn)。

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

切割路徑定義為切割工具移動(dòng)過程中與表面網(wǎng)格的一系列交點(diǎn)。在實(shí)現(xiàn)過程中，計(jì)算刀刃進(jìn)入三角面片時(shí)與入邊的交點(diǎn)和移出三角面片時(shí)與出邊的交點(diǎn)，并依據(jù)兩個(gè)交點(diǎn)與頂點(diǎn)之間的關(guān)系進(jìn)行頂點(diǎn)移動(dòng)。

2.2 移動(dòng)頂點(diǎn)[2]

在計(jì)算得到入邊和出邊交點(diǎn)后進(jìn)行頂點(diǎn)移動(dòng)。頂點(diǎn)移動(dòng)遵循的原則是就近原則。即根據(jù)交點(diǎn)坐標(biāo)位置判斷與交點(diǎn)相鄰的兩個(gè)頂點(diǎn)到該交點(diǎn)的距離大小，將距離交點(diǎn)較小的頂點(diǎn)移動(dòng)到交點(diǎn)的位置處。交點(diǎn)移動(dòng)過程中有兩種特殊情況需要注意：

情況一：兩出入交點(diǎn)距離三角形某個(gè)頂點(diǎn)非常接近，且該頂點(diǎn)到兩個(gè)交點(diǎn)的距離大小相等。在這種情況下將頂點(diǎn)移動(dòng)到交點(diǎn)所在邊的邊長(zhǎng)較大的交點(diǎn)處，另一個(gè)交點(diǎn)也移動(dòng)到該處。

情況二：兩出入交點(diǎn)分別位于出入邊的中點(diǎn)處，則將除兩邊公共頂點(diǎn)之外的兩個(gè)頂點(diǎn)移動(dòng)到相應(yīng)得交點(diǎn)位置。

圖6 頂點(diǎn)移動(dòng)時(shí)的特殊情況

可見，經(jīng)過頂點(diǎn)移動(dòng)后，不再出現(xiàn)狹長(zhǎng)三角形或與原模型中的三角形不同數(shù)量級(jí)的小三角形，有利于在網(wǎng)格上進(jìn)行后續(xù)操作，如繼續(xù)切割或變形計(jì)算等。

2.3 切割路徑分離

頂點(diǎn)移動(dòng)后，通過將位于切割路徑上的頂點(diǎn)分裂為兩個(gè)來分離網(wǎng)格，顯示切割路徑。兩個(gè)分裂點(diǎn)間的連線垂直于切割平面，且二者與切割平面間的距離相等，如圖7所示。

圖7 切割路徑分離方法

3 切割搜索

切割過程結(jié)束后如果被切下的網(wǎng)格部分與原始網(wǎng)格不連通則需要將二者分離，這就需要采用網(wǎng)格搜索算法。網(wǎng)格搜索算法中包含網(wǎng)格拓?fù)浣Y(jié)構(gòu)表示和在拓?fù)浣Y(jié)構(gòu)下的搜索算法。

3.1 AIF數(shù)據(jù)結(jié)構(gòu)

AIF（adjacency and incidence framework）數(shù)據(jù)結(jié)構(gòu)[5]通過構(gòu)建多邊形間的鄰接入射關(guān)系來表示多邊形網(wǎng)格的拓?fù)浣Y(jié)構(gòu)。該法可以高效地實(shí)現(xiàn)創(chuàng)建和搜索，且在切割過程中可以方便地維護(hù)。

AIF通過描述網(wǎng)格結(jié)構(gòu)的基本元素來表達(dá)網(wǎng)格的拓?fù)浣Y(jié)構(gòu)。當(dāng)一種元素x是另一種元素y的邊界之一，且其幾何維數(shù)低于y，稱x鄰接于y，表示為x＜y，或y入射于x，表示為y＞x。

當(dāng)網(wǎng)格由三角面片組成時(shí)，基本元素包括點(diǎn)、線和面。利用網(wǎng)格結(jié)構(gòu)的鄰接和入射關(guān)系，AIF數(shù)據(jù)結(jié)構(gòu)構(gòu)造包含點(diǎn)集、邊集和面集及表征它們之間關(guān)系的數(shù)據(jù)集合。各集合中每個(gè)元素都包含與之有鄰接和入射關(guān)系的其他元素的子集合，即每個(gè)頂點(diǎn)數(shù)據(jù)中包含所有通過該頂點(diǎn)的邊數(shù)據(jù)的集合，每條邊數(shù)據(jù)中包含所有經(jīng)過該邊的所有面的數(shù)據(jù)集合，每個(gè)面數(shù)據(jù)中包含所有構(gòu)成該面的所有邊的數(shù)據(jù)集合。通過AIF數(shù)據(jù)結(jié)構(gòu)，可以利用一種元素得到所有與之有鄰接關(guān)系的數(shù)據(jù)元素。數(shù)據(jù)結(jié)構(gòu)實(shí)現(xiàn)如下：

Class Vertex{

float x，y，z；//點(diǎn)的空間坐標(biāo)

vector＜Edge*＞edgeSet；//經(jīng)過該點(diǎn)的所有邊的集合

}

Class Edge{

Vertex* p1，*p2；//組成該邊的兩個(gè)頂點(diǎn)的指針

vector＜Face*＞faceSet；//經(jīng)過該邊的所有面的集合

}

Class Face{

vector＜Edge*＞edgeSet；//構(gòu)成該面的所有邊的集合

}

3.2 基于AIF結(jié)構(gòu)的搜索算法

利用AIF數(shù)據(jù)結(jié)構(gòu)，可以從一種網(wǎng)格基本元素得到與之相鄰的其他網(wǎng)格基本元素?？刹捎脠D的廣度優(yōu)先算法來進(jìn)行網(wǎng)格數(shù)據(jù)搜索。

圖的廣度優(yōu)先遍歷算法的基本思想是首先從圖中某個(gè)頂點(diǎn)V0出發(fā)，訪問該頂點(diǎn)后訪問各個(gè)未曾訪問過的鄰接點(diǎn)W1，W2，…，Wk，然后依次從W1，W2，…，Wk出發(fā)訪問各自未被訪問的鄰接點(diǎn)。如此反復(fù)直到所有點(diǎn)都被訪問過為止。

為了實(shí)現(xiàn)上述算法，本文在AIF的點(diǎn)和面結(jié)構(gòu)中加入了reach屬性，表示點(diǎn)和面是否已經(jīng)被訪問過。同時(shí)還引入隊(duì)列結(jié)構(gòu)輔助算法實(shí)現(xiàn)。首先將初始頂點(diǎn)的reach屬性設(shè)置為true并將該點(diǎn)加入隊(duì)列中，網(wǎng)格結(jié)構(gòu)中其他元素的reach屬性設(shè)置為false。若隊(duì)列非空，取出隊(duì)列中的首元素，通過AIF數(shù)據(jù)結(jié)構(gòu)得到所有通過該點(diǎn)的網(wǎng)格面，若網(wǎng)格面的reach屬性為false，設(shè)置其屬性為true并利用AIF數(shù)據(jù)結(jié)構(gòu)得到該網(wǎng)格面的所有頂點(diǎn)，設(shè)置reach屬性為false的網(wǎng)格點(diǎn)為true并將這些網(wǎng)格點(diǎn)加入隊(duì)列尾部，如此反復(fù)直到隊(duì)列為空。

網(wǎng)格搜索算法比較簡(jiǎn)單，但搜索過程中初始點(diǎn)的選擇非常重要。在網(wǎng)格切割過程中，通常被切下的網(wǎng)格部分相比原來的網(wǎng)格結(jié)構(gòu)是非常小的，即包含數(shù)量較少的點(diǎn)、線和面數(shù)據(jù)。此時(shí)如果選取被切下網(wǎng)格中的頂點(diǎn)作為搜索起點(diǎn)，則搜索過程中需要遍歷的網(wǎng)格節(jié)點(diǎn)較少，計(jì)算量較小，而如果選擇原始網(wǎng)格上的網(wǎng)格頂點(diǎn)作為搜索起始點(diǎn)，則搜索算法的計(jì)算量很大。

圖8中原始網(wǎng)格被剖分為兩個(gè)網(wǎng)格G1和G2，G2為一封閉區(qū)域且不與G1連通。切割的起始點(diǎn)分別為S1和S2?？梢娙魪腉1中的頂點(diǎn)（如S1）開始進(jìn)行網(wǎng)格搜索，最終將遍歷G1，并將其從兩個(gè)網(wǎng)格結(jié)構(gòu)中分離，但G1包含的數(shù)據(jù)量遠(yuǎn)遠(yuǎn)大于G2，從切割的實(shí)時(shí)性考慮應(yīng)從G2中的頂點(diǎn)（如S2）開始搜索并將G2分離出來。

圖8 網(wǎng)格分離時(shí)起始點(diǎn)的選擇

通常選擇G2中的S2作為起始點(diǎn)比較方便。為了選擇起始搜索點(diǎn)S2，在切割開始時(shí)保存切割起始點(diǎn)S1和S2，并在每次切割后計(jì)算包含S1的切割路徑path1的長(zhǎng)度和包含起始切割點(diǎn)S2的切割路徑path2的長(zhǎng)度。由上圖可見，由于閉合的切割路徑path2在path1內(nèi)部，所以path2的長(zhǎng)度小于path1，所以切割結(jié)束后選擇切割路徑長(zhǎng)度小的起始切割點(diǎn)作為網(wǎng)格搜索的起始點(diǎn)。

4 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)中將上述算法應(yīng)用到了CT顱面數(shù)據(jù)中。首先采用marching cube算法提取等值面網(wǎng)格，網(wǎng)格結(jié)構(gòu)由三角面片組成，然后移動(dòng)切割工具并進(jìn)行碰撞檢測(cè)和網(wǎng)格切割，切割結(jié)束后找到被切下網(wǎng)格上的頂點(diǎn)進(jìn)行網(wǎng)格分離。

在切割交互中，移動(dòng)切割工具并進(jìn)行其與物體網(wǎng)格的碰撞檢測(cè)。碰撞檢測(cè)采用1中所述的OBB包圍盒樹算法。雖然OBB包圍盒在物體拓?fù)浣Y(jié)構(gòu)發(fā)生改變后需要更新，且更新算法比較復(fù)雜，但考慮到首次碰撞后切割工具連續(xù)移動(dòng)，所以僅在初次碰撞后采用碰撞檢測(cè)得到的三角面片與切割三角面片進(jìn)行相交計(jì)算，后續(xù)的相交計(jì)算則依據(jù)切割工具的移動(dòng)方向和AIF數(shù)據(jù)結(jié)構(gòu)進(jìn)行。因切割工具沿切割刀片所在的平面移動(dòng)，所以在計(jì)算得到當(dāng)前三角面片的出邊信息后由AIF結(jié)構(gòu)找到下一個(gè)三角面片，如此反復(fù)計(jì)算直到得到某個(gè)三角面片與當(dāng)前切割三角面片相交，以該三角面片為當(dāng)前碰撞面并進(jìn)行相交計(jì)算。

圖9為本文算法在醫(yī)學(xué)顱面整復(fù)手術(shù)中的切割應(yīng)用。圖（a）為由marching cube算法提取的三維等值面網(wǎng)格結(jié)構(gòu)，其包含兩層等值面結(jié)構(gòu)。外層為顱面表皮，內(nèi)層為骨組織表面。圖（b）和圖（c）為圖（a）網(wǎng)格數(shù)據(jù)的填充渲染模式。圖（d）為顱面表皮被切割后的切割路徑顯示。圖（e）為分離被切下網(wǎng)格結(jié)構(gòu)與原始網(wǎng)格結(jié)構(gòu)后顯示出內(nèi)層骨組織。

圖9 實(shí)驗(yàn)結(jié)果

5 結(jié)語

文中針對(duì)醫(yī)學(xué)導(dǎo)航系統(tǒng)中的虛擬切割算法進(jìn)行了系統(tǒng)仿真。對(duì)基于表面重建的三角網(wǎng)格數(shù)據(jù)建立AIF數(shù)據(jù)結(jié)構(gòu)以維護(hù)其空間拓?fù)浣Y(jié)構(gòu)，然后在移動(dòng)切割工具的過程中采用OBB包圍盒樹的方法進(jìn)行碰撞檢測(cè)，檢測(cè)到碰撞發(fā)生后采用基于頂點(diǎn)移動(dòng)的網(wǎng)格切割方法進(jìn)行網(wǎng)格切割，切割結(jié)束后若切割路徑為封閉曲線，則對(duì)由切割路徑包圍的網(wǎng)格結(jié)構(gòu)進(jìn)行搜索遍歷并與原網(wǎng)格分離。實(shí)驗(yàn)結(jié)果表明文中的系統(tǒng)算法能夠較好的實(shí)現(xiàn)切割要求。通常人體組織為軟組織，在切割過程中會(huì)有形變的產(chǎn)生，本文中沒有涉及變形計(jì)算，后續(xù)研究將在三角網(wǎng)格結(jié)構(gòu)的基礎(chǔ)上考慮加入切割過程中的網(wǎng)格變形計(jì)算，使切割仿真更加真實(shí)。

[1]潘振寬，崔樹娟，張繼萍等.基于層次包圍盒的碰撞檢測(cè)方法[J].青島大學(xué)學(xué)報(bào):自然科學(xué)版，2005，18:71 76.

[2]王鈺，于素平.針對(duì)面模型的頂點(diǎn)移動(dòng)切割算法研究[J].系統(tǒng)仿真技術(shù)，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的三角形網(wǎng)格切割方法[J].上海交通大學(xué)學(xué)報(bào)，2008，42（004）: 564-568.