吳德道,劉小平
(1. 南昌大學信息工程學院,江西 南昌 330031;2. 景德鎮陶瓷大學,江西 景德鎮 333000;3. 卡爾頓大學系統與計算機工程系,渥太華 加拿大 K1S 5B6)
在虛擬手術仿真系統中建立自然、準確的軟組織切口模型是該領域研究的重點之一[1]。主要體現在虛擬手術切口生成后切口面幾何和拓撲結構發生了明顯的變化,特別是網格化方法中需要重網格化,這將給模型的仿真帶來很大的問題[2],[3]。而經過一段時間的發展,研究者們提出許多模擬軟組織切口模型[4]。根據建模過程中是否需要對虛擬器官和組織對象進行細分和預處理,可以將軟組織切口模型分為基于網格和無網格的兩種建模架構[5]。虛擬手術切口模型的仿真涉及到使用網格切割技術修改多邊形或三維網格。網格切割是網格編輯研究的一個重要組成部分,包括重映射、變形、合并、簡化和平滑操作。目前,關于網格切割的研究大多集中在虛擬手術和醫學仿真領域,其中較為常見的是三維掃描和三維網格表示。而大多數醫學領域以外的軟件只處理多邊形網格的應用,使得該技術在網格表面操作上更加成熟。曲面網格的使用不僅提高了基于多邊形的渲染引擎的可用性和性能,而且還是實時應用的關鍵因素。
早在1997年,Nielsen等首次提出基于網格的單元去除法[6]。為提高切口的仿真效果,Nielsen等提出改進模型,如邊緣折疊法、頂點復制分裂法、曲面約束法等。這種切口模型只處理虛擬物體的表面網格,使用較少的點和面來表示物體的表面形狀,但不能表達物體內部的信息。為更好地模擬和顯示物體被切割后的各種內部變化,Bielser等人提出一種基于體積模型[7]的體積元細分法。之后,相繼出現虛擬節點法、部分剖分法和立方體單元劃分法等改進模型。所有上述模型需要將預處理網格細分[8,9],在切割過程中進行網格更新和重組,成本很多。2009年,丹尼斯等人結合對象空間色散的無網格方法和基于網格表面重建技術,并成功地應用在子宮鏡檢查和切削仿真。但該方法數據結構復雜,算法復雜度高,計算效率低。而目前混合模型最大的困難是如何管理好兩個不同的建模系統[10]。
根據Bruyns等人[11]的研究,他們把網格切割技術分為重網格化、切割路徑的定義和新原語的創建。前兩種方法是在間隔點之間或在選定點之間構造一個連續平面,而它們的區別是如何設置切割路徑。模板與網格基元的碰撞檢測是第三種方法的要求。此外,Bruyns等還描述了漸進切割方案。為實現刀具的原始碰撞,在刀具連續移動位置之間構造一個掃描面。對原語測試曲面的掃描生成邊和面的交點,其中新節點連接到新創建的邊和面,然后替換原語,該方案適合多層目標的碰撞檢測。而對于非漸進式切割方法,存儲的信息大多是冗余的,而重網格化方案只處理邊-邊和邊-面相交的情況。此外,切割暴露了表面網格的內部,因此該方法只適用于多層網格。
Zhang等人[12]提出一種表面網格模型的漸進式切割技術,該技術將網格元素細分,并沿選擇的切割路徑構造內部結構。在切割線與三角面邊緣相交的地方創建兩個新節點。這兩個節點初始是一樣的,但隨著切割完成而分離。這種分離留下一個切口槽,其深度由刀具尖端的切割深度來定義。在此基礎上,他們還提出一種漸進式切割算法,通過執行臨時細分和動態構造切口槽來跟蹤刀具在三角形中的運動。與前一種方案一樣,該方法沒有一個完整的網格重映射方案,并且沒有網格表示的定義,這使得該方法難以重現。
Van der Stappen和Nienhuys[13]定義一種切割技術,通過限制原始元素不具有短邊或大角度使得生成的新元素較少,從而產生形狀良好的網格結構。切割工具的任意一次移動都會移動一個額外的活動節點。在每次移動之后,對應三角形將被重新映射到本地。通過這種方式,刪除切割工具前面的節點,然后將其插入到工具之后。該方法的特點是切割前后目標網格的分辨率保持不變,但切割路徑的精度不能得到提高,而這是由初始網格的分辨率決定的[14]。
Lim等人[15]提出一種基于節點捕獲的漸進式切割算法。第一次碰撞接觸時,最近的節點被拉到目標網格和工具之間的初始接觸位置。在切割過程中,切割路徑上的每個節點被分成兩個方向相反、垂直于切割線的節點。最后,設置一個切割槽來創建體積切割的錯覺。為解決切割路徑的切割分辨率依賴于周圍網格的細分問題,他們還提出一種局部細分算法[16]。這種方法允許更準確地表示切割路徑,但內存成本高且性能降低。
針對上述建立切口模型方法的優缺點,提出一種實時模擬切口模型的方法。該方法直接與表面網格交互,并在程序運行時生成新的幾何拓撲,為后續切口模型的高真實感渲染提供一個切口內腔結構。本文其余部分組織如下:第2節中給出建模方案中使用的模擬方法和相應算法。第3節通過實驗進一步證實仿真方法的有效性。最后,在第4節對本文的研究成果進行總結,并對未來的研究提出建議。
通常當一個鋒利物體穿過真實軟組織表面撕裂表皮使皮下組織可見時,就在軟組織表面留下一個切口。利用網格編輯技術對切割影響的網格區域進行局部修改,從而可在虛擬環境中重現這一現象。本節中,將研究與切口模擬相關的算法和技術,包括將切割線與目標網格融合并將切割線擴展成切口以及生成切割槽。
在目標網格中插入一條切割線段作為模擬切割邊界,會改變原始網格的許多屬性并破壞原始網格的完整性。因此,為避免這一問題,在插入新的切割線之前,需要對切割線周圍的局部幾何拓撲進行細分。圖1描繪出切割線與目標網格的融合。對于每段切割線段,端點可以與三角形網格的節點重合,也可以是它的一條邊,或者兩者都在三角形網格內等。

圖1 插入的切割線融合到目標網格中
按照切割線與目標網格共點情況劃分為三類:第一類是不共點情況,切割線端點與網格頂點沒有重合也是最常出現的情況如圖2上所示。第二類是共點情況,切割線端點與網格面頂點重合,如圖2下所示。第三類是比較特殊的情況,整個切割線包含在一個三角網格內,由于切割線太短,這種情況的切割沒有任何實際意義,因此在程序中忽略處理。切割線段p0p1與三角網格f之間的每種位置關系需要不同的方式將其融合到f中。根據剖分線的定義和網格的屬性對網格表面進行細分。新插入的端點為p0和p1的邊ec被添加到邊ECs的有序集中。圖2給出各位置關系,其中切割線段用虛線表示,網格細分用點劃線表示。具體操作是:

圖2 切割線與三角網格的位置關系
第一類,如圖2上,從左至右分別對應以下3種:
1)p0在網格內p1在邊上,這種位置關系常發生在切割線的初始位置。三角網格f在p0處用3-split算法細分,在p1處執行2-split操作。
2)p0和p1均在邊上,這種位置關系最常見。切割線段p0p1將網格f分成一個三角形和一個四邊形,而為保持三角網格關系,首先在p1處分出兩個三角網格(點劃線),然后在p0處分割一個新的三角子網格。通過比較p0所在的邊及其邊引用,可以確定p0與兩個三角子網格中的哪一個相關聯,均采用2-split操作。
3)p0在邊上p1在網格內,這種位置關系常發生在切割線的末端。將f在p1處按3-split操作進行細分,然后確定p0在哪個子網格上,在p0處繼續進行2-split操作。
第二類,如圖2下,從左至右分別對應以下3種:
1)p0(p1)和網格共點p1(p0)在網格內,p1(p0)處將f細分為三個子網格并被共享。3-split操作在f內部創建三個新的面和邊,它們都與p1(p0)有關。其中一條新邊與ec重合,并被添加到ECs上。
2)p0(p1)和網格共點p1(p0)在邊上,切割線段將f細分為兩個新的三角網格。使用2-split算法將f在p1(p0)處分割成兩個子網格,并更新f的直接鄰域。
3)p0和p1均與網格共點,即切割線段p0p1與網格共邊,直接添加到ECs中。
算法1結合圖3左側詳細描述了2-split操作。與圖3右側相似的組合可以推導出3-split操作。

圖3 面f分別在2-split (左)和3-split操作(右)時的拓撲
算法1 2-split
輸入: 面f,點p和點p所在的邊es
輸出:分割邊ec,將f從p細分到它的對面節點n1
1)獲取節點{n0,n1,n2},頂點{v1,v2,v3},面f上的邊{e0,e1,e2}
2)在p點處創建節點nm,頂點vm
3)在nm和n1之間生成切割邊ec
4)分別在n0和nm,nm和n2之間生成邊ex0,ex1
5)分別在nm、n0和n1,nm、n1和n2之間生成面fc0,fc1
6)設置邊ec的鄰接面fc0,fc1
7)分別設置邊ex0,ex1的鄰接面fc0,fc1
8)設置面fc0的鄰接邊ex0,e0,ec
9)設置面fc1的鄰接邊ec,e1,ex1
10)分別用fc0,fc1取代邊e0的鄰接面f和邊e1的鄰接面f
11)從網格中刪除面f,邊es
12)返回分裂邊ec
該階段由兩部分組成:切割線擴展和切割槽生成。輸入為前一階段創建的一組有序切割線線段集合,分離切割邊的局部幾何形狀,模擬切口。此操作在目標網格中留下一個凹槽,然后融合這個位于當前網格面以下某個深度的新拓撲。圖4描繪出切割線的分離幾何形狀及其內部結構。

圖4 復制相應的頂點并生成切口邊界與其內腔結構
根據產生切口的要求,切割線上的節點和頂點要相互復制,分別移動。將這些元素與邊緣連接在切口周圍形成一個邊界(圖4中的粗實線),該邊界應與切割線兩側的拓撲結構融合。初始切割線上的節點、頂點和邊不再引用,它們在之后被刪除。根據Lim等人[16]測量的切口張開距離(IOD)與切口長度L、深度d的關系,可以得到一般公式
IOD(L,d)=(0.0111+0.0002(d-0.1/0.02))ln(L)
+0.0415+0.0015(d-0.1/0.02)
(1)
通過可調參數拋物線來調節切割線上任意點的切割位移計算IOD。這條拋物線確保IOD在中心近似為最大值,在開始和結束為零,在中間光滑。由此得到計算方程
IODpi=mIOD(L,d)
(2)
式2中,m可由拋物線方程m=ax2+bx+c確定,其中,a、b、c為可調系數,x=i/|Ec|,i=0,…,|Ec|,|Ec|為切割線邊集合Ec的基數。
定義切割線邊界的新幾何點用于創建網格拓撲,并將切割線劃分為兩個相似的部分。算法2詳細描述了這一過程。第一條邊和最后一條邊的處理方式類似,分別對位于p0和p1中的節點和頂點進行重用。切割深度d的計算公式為
d=max(0.1,min(1.0,l/5))
(3)
算法2 構建切割邊界的幾何拓撲
輸入: 有序的切割邊集Ec
While(對每一條邊ec∈Ec)do
1)獲取與邊ec對應的端點
2)在垂直切割線方向上計算vi ± IODpi/2,得到兩條邊界線的新端點
3)將這4個端點賦值給拓撲結構中的節點
4)生成對應4個頂點
5)由
6)由
7)分別將
8)分別將
9)分別將e+和e-添加到集合E+和E-
10)分別將邊ec的面f0和f1添加到集合F+和F-
對于給定的切割邊ec,計算p0i和p1i在切割槽中的位置如下

(4)


(5)
通過算法3,創建切割槽結構并與模型表面網格融合。
算法3 模擬切口內腔的幾何拓撲
輸入: 切割邊集Ec、E+、E-;節點集N+,NI,N-,內腔頂點集W+、WI、 W-
While(對每一條邊ec ∈Ec)do
1)獲取圖5中所示節點和頂點
2)兩點生成邊:
3)3個節點和對應3個頂點生成面:
4)設置邊面關聯:
5)注冊邊面關聯:
6)注冊相交邊面關聯:
7)注冊內腔面中邊面關聯:

圖5 切口內腔的網格元素
實驗環境配置為便攜式計算機,CPU為Intel Core I5 5200,內存為4GB,顯卡為AMD Radeon R5 M320。仿真軟件采用VC++的開發,圖形渲染方面使用OpenGL 圖形庫,材質語言為CG。實驗案例是通過三維掃描重建真實肝臟的數字模型。
在肝臟三維數字模型中,整個肝臟表面由3968個三角網格面包裹。圖6為三角形網格面相互銜接形成的三維網格模型,以及切口生成后目標網格上形成的拓撲結構,著色模型以及切口部分放大后內腔拓撲結構。在圖6中,切割切口完全融入到目標網格中并生成具有一定深度的內腔拓撲。切口平滑自然,沒有破壞目標網格的完整性。

圖6 數字肝臟原始網格模型以及生成切口后的網格拓撲注:按自上而下,從左到右順序,分別是原始網格模型、切口模型、著色后切口模型以及局部放大切口模型。
圖7中可從多個角度觀察切口的視覺層次和切口的真實深度拓撲。最終模擬的切口結構過渡自然,并可以根據目標場景修改切口和腔體實際對應紋理,易于實現。

圖7 著色渲染后仿真肝臟模型和生成切口后多角度演示的模型
與真實肝臟軟組織的切口對比如圖8所示,其中左圖是固定切割裝置旁邊是切割產生切口的肝臟軟組織[17],右圖是采用本文方法對虛擬肝臟切口的渲染。通過對比,可以注意到切口幾何外形基本相似,所不同的是切口周圍血液顏色的渲染,左圖偏向切口外形實驗,右圖兼顧切口外形和周圍表皮顏色的渲染。

圖8 與真實肝臟切口的對比注:左圖為真實肝臟切口,右圖為仿真切口。
圖9中,左上是DFFD方法[18]模擬的手術切口,右上采用的是基于FEM的DFFD方法[19],中間是本文方法,下圖是本文方法切口局部放大。前兩種方法沒有構造切口內腔結構缺乏立體逼真感,手術切口與表皮之間缺乏過渡元素,不利于后期切口的進一步渲染。

圖9 不同仿真方法在相同部位生成的切口
通常12幀/秒的幀速率是人眼實時交互所能達到的最低幀速率。在數字模擬虛擬肝臟模型的三角網格數為3968的情況下,本文模擬手術切口(600PP)的方法獲得的幀速率約為230幀/秒,其它渲染過程所對應幀速率如表1所示。模擬切口在世界坐標(0.02~0.15)下的總時間為6~10ms,保證了后繼操作的實時性。

表1 虛擬肝臟模型渲染過程及手術切口所對應的幀速率
本文提出一種虛擬手術切口模型的仿真方法,簡化了切割線與目標網格面的對應關系,利用可調參數曲線模擬切口邊界。通過減少對切割線周圍網格面的修改來保證網格結構的完整性和穩定性并以虛擬肝臟為實驗對象,實現了對肝臟切口的仿真。該方法不僅適合虛擬手術切口仿真,而且創建了具有立體感的內腔結構,可實現高真實感虛擬手術切口的實時渲染。本文主要貢獻是詳細描述了模擬網格切口所需的技術,并保持原始網格拓撲結構和參數的完整性。仿真程序能夠通過基本的重網格化操作處理所有的切割線與三角形網格之間的關系,并為模擬目標網格上的切割模型提供完整的數據結構和算法。優點是對網格拓撲的干擾最小,因為只修改與切口邊界有關的網格面。本文采用的是基本操作和數據,能夠支持不同類型的虛擬手術切口。在判斷兩個空間點是否重疊時,網格精細度問題會導致拓撲網格失效。在確定網格節點和頂點時,通過使用閾值來衡量兩點之間的距離是否小到可視為重疊的程度。因此,閾值需根據網格表面節點密度等因素進行調整。最后,無需執行額外的網格細化,只對與切割線相關的面進行重映射。在未來,解決方案可以聚焦在新創建切口周圍的網格變形上。