重力作用下肝臟模型變形及切割模擬仿真

張峰峰， 張石， 黃柯， 于凌濤， 孫立寧

(1.蘇州大學 機電工程學院,江蘇 蘇州 215006； 2.蘇州大學 蘇州納米科技協同創新中心,江蘇 蘇州 215123； 3.哈爾濱工程大學 機電工程學院,黑龍江 哈爾濱 150001)

部分肝臟切除術是目前治療肝癌的根治性手段，同時也是治療一般良性肝臟疾病的有效手段[1-2]。隨著醫學技術的發展，將增強現實或者虛擬現實等輔助手術系統應用于肝臟切割手術過程中，能夠解決傳統肝臟手術過程中存在的術前和術中信息不同步的問題[3-4]。在實際的肝臟手術過程中，由于肝臟是非剛體，當手術器械對肝臟進行按壓或者切割時會產生相應的形變，對于肝臟形變模擬的真實程度在很大程度決定了輔助手術系統的準確性，為準確反映實際肝臟形變特征，建立一個合理的肝臟組織生物形變模型至關重要。根據切割肝臟過程中變形是否具有真實依據可分為物理形變模型和非物理形變模型[5]。

非物理形變模型主要根據用戶輸入或數學計算直接改變組織模型形態對形變進行模擬，主要包括基于控制點的模型[6]和自由形變模型[7]。物理形變模型依據現實的物理定律構建變形的物理模型[8]，主要包括邊界元模型[9]、有限元模型[10-11]和質點彈簧模型[12-14]，其具有準確度高、形變真實的特點，也因計算量大導致較低的刷新率。本文采用基于質點彈簧的肝臟模型，分別在重力作用下變形以及重力作用下切割進行分析，依據肝臟的生物力學特性建立了基于質點彈簧的物理形變模型，進行部分肝臟切除術中肝臟變形仿真。

1 質點彈簧軟組織模型

1.1 基于Tetgen的四面體剖分

通常由CT重建獲得的模型為三角網格表面模型，使用彈簧阻尼結構將模型表面進行連接。重力作用下不足以支撐其三維結構[15]，因此需要在肝臟模型內部插入新頂點，并將其與表面頂點連接，形成四面體或六面體等結構。

由給定點集生成四面體的過程稱為四面體剖分，四面體剖分的方法主要有點集網格剖分、區域網格剖分和限定網格剖分。點集網格剖分要求給定點集，生成的四面體頂點均為點集內的點；區域四面體剖分要求給定一個區域邊界，在邊界內部生成四面體，內部四面體的頂點沒有限制，可在邊界內部任意插入頂點；可對點、邊或面添加控制條件。本文主要針對由CT重建出的三維表面網格模型進行四面體剖分，模型為肝臟外表面的邊界數據，因此使用區域固體剖分將肝臟模型四面體化。

Tetgen為開源四面體剖分工具，其可接受3D邊界網格模型作為輸入，在其內部插入新頂點，生成四面體單元。本文使用QT結合Tetgen編寫了一個通用的四面體剖分程序，程序可輸入任意三維網格模型，輸出四面體模型。程序主要包括網格質量控制部分和輸出文件的選擇部分。

網格質量控制部分包括一個保留表面網格的選項和一個質量控制滑塊。當選擇保留表面網格時，程序在進行四面體剖分時將保留原模型的表面網格，只對模型內部進行操作。否則，程序在進行四面體剖分時將同時對表面網格進行重構。質量控制滑塊能夠對模型內部四面體單元的質量進行控制。四面體單元的質量由2種方法進行控制：半徑-邊緣比約束和最小二面角約束。四面體T的半徑-邊緣比ρ(T)為：

(1)

式中:r為T的外接圓半徑；d為四面體最短邊長度；θmin是四面體T的最小二面角。因此，半徑-邊緣比約束等價于最小二面角約束。

一個四面體模型由多個文件組成，程序輸出的文件為：1)nodes文件為模型頂點列表的儲存文件，每一行對應一個頂點的坐標值；2)face文件為模型三角面片列表的儲存文件，每一行對應一個三角面片的3個頂點。face文件包括模型表面的三角面片和內部面片。face文件最后一列為標識位，用來區分表面和內部面片，1為表面面片，0為內部面片；3)ele文件為四面體單元列表的儲存文件，每一行對應一個四面體單元4個頂點在nodes文件中的索引；4)edge文件為四面體邊列表的儲存文件，每一行對應一條邊的2個頂點。edge文件同樣包括模型表面的邊和內部的邊，最后一列標識位同face文件；5)neigh文件文件為鄰接四面體列表的儲存文件。neigh文件每行有5列，第1列為編號，后4列為一個四面體4個頂點對面的四面體，若四面體位于表面則該四面體只有3個鄰接四面體，另一個位置表示為-1作為占位符。

使用上述四面體剖分程序對肝臟表面網格模型進行四面體剖分，輸出結果如圖1所示。

圖1 四面體剖分結果Fig.1 Tetrahedral meshing results

1.2 質點彈簧力學模型

應力松弛、蠕變和滯后統稱為粘彈性特性，肝臟組織的粘彈性特性是肝臟變形過程中的重要性質。一個物體突然發生形變，之后保持該形變，應力隨時間逐漸減小的過程稱為應力松弛。若一個物體內突然產生應力并保持應力不變，物體持續發生變形的現象稱為蠕變。若一個物體承受循環載荷，但加載時應力應變關系與卸載時的應力應變關系不同的現象稱為滯后。目前用于描述物體粘彈性性質的模型主要有Maxwell模型、Voigt模型和Kelvin模型[16]，如圖2所示，圖中μ和η分別為彈簧的彈性系數和粘性系數;σ為應力；ε為應變。Maxwell模型由一個線性彈簧和一個阻尼器串聯而成。當物體突然產生應變時，彈簧將儲存彈性勢能，阻尼器在彈簧的影響下逐漸發生變形。因此，Maxwell模型更加適用于描述物體的應力松弛現象；Voigt模型由一個線性彈簧和一個阻尼器并聯而成。當物體突然產生應力時，由于阻尼器的存在，物體并不會立即產生應變，而是在阻尼器的影響下逐漸變形。因此，Voigt更適合描述物體的蠕變現象；Kelvin模型可看作將一個Maxwell模型與一個線性彈簧并聯而成。Kelvin模型能夠同時描述物體的應力松弛和蠕變現象。

圖2 質點彈簧模型示意Fig.2 Schematic diagram of the mass point spring model

本文主要針對肝臟在重力作用下的變形和切割，其過程更加接近于蠕變。因此采用Voigt模型構建肝臟的物理形變模型。根據四面體剖分軟件導出的.edge文件將肝臟模型中的頂點使用Voigt元件連接起來，其中每一個頂點受力方程為：

fext(i)+fvoigt(i)=fsum(i)

(2)

式中：fext(i)為質點i受到的外力；fvoigt(i)為質點i受到Voigt元件產生的力；fsum(i)為質點i受到的合力。

質點i受到的外力由手術器械的交互力fc(i)和重力fg(i)組成：

fext(i)=fc(i)+fg(i)=Δs·pi+mig

(3)

Voigt元件產生的力由阻尼力fdmp和彈簧力fspr組成：

(4)

式中：v為質點的移動速度，r為質點的空間坐標，nij為連接質點i和質點j的Voigt元件的阻尼系數，kij為連接質點i和質點j的Voigt元件的剛度系數。

質點i受到的合外力由牛頓第二定律得出：

fsum(i)=miai

(5)

在肝臟受重力變形和切割變形的過程中將肝臟組織視為各向同性的結構，即mi=m,ηij=η,kij=k,則有：

(6)

計算時對時間進行離散化處理，每隔Δt=0.02 s進行一次計算。計算質點的加速度a，則Δt后質點的速度為：

vi+1=vi+Δv=vi+a·Δt

(7)

質點在Δt內的運動可視為勻速運動，因此，Δt后質點的坐標為：

ri+1=ri+vi·Δt

(8)

在Δt時間內計算所有質點的位置坐標即可得出0.02 s后肝臟模型的整體形態，當所有頂點的受力達到平衡時，肝臟模型不再變形，此時即為肝臟模型在重力作用下的形態。

2 離體豬肝切割及變形實驗

本文使用質點彈簧模型模擬肝臟在平面上受重力影響的變形和在重力影響下的切割變形，本節以豬肝為研究對象，分別對以上2種情況進行實驗驗證。

2.1 肝臟模型在重力下的變形實驗

要對肝臟模型在重力影響下的變形進行驗證，必須得知肝臟在失重狀態下的原始形態并獲得其失重模型。因此需要模擬一個失重的環境，將豬肝置于其中，獲取豬肝在失重狀態下的三維數據。

通常模擬失重的方法有直接法和間接法。直接法以某種方式使物體進入真正的失重狀態，一種方式是使物體進行自由落體運動，另一種方式為拋物線飛行。間接法為物體提供一個類失重的環境，主要有浸水法和風洞法，還有臥床法和懸吊法等，但后面2種方法旨在模擬失重狀態下其他影響，并不能很好地模擬失重狀態下物體的形態。

在本實驗中為了保持豬肝在失重狀態下的形態并獲得其三維數據，采用浸水法模擬失重環境。研究表明豬肝的密度與水接近略大于水，因此可將豬肝浸入水中，然后向水中加鹽增加水的密度直至豬肝在鹽水中處于懸浮狀態。為了獲取真實的切割路徑，在豬肝上固定2個小球，切割時在2個小球之間沿直線切割。當豬肝在鹽水中穩定后對其進行CT掃描，經三維重建獲得懸浮狀態下的豬肝模型，如圖3所示。

圖3 懸浮狀態下的豬肝模型Fig.3 Pork liver model in suspension

CT掃描完成后將豬肝從鹽水中取出放置在平面上再次掃描，重建得到豬肝在重力下的形變模型作為對比標準。將豬肝的懸浮模型導入Unity并賦予程序腳本，新建一個平面，將豬肝模型置于略高于平面的上方，運行程序模擬豬肝在重力下的變形。待豬肝完全落地并穩定之后導出變形后的模型，如圖4所示。

圖4 Unity導出并實體化的重力形變模型Fig.4 Gravity deformation model exported and solidified by Unity

2.2 肝臟模型在重力下的切割實驗

將取出豬肝沿預先設定的切割路徑進行切割形成切口，具體切割算法流程已在文獻[17]進行了闡述，隨后再次進行CT掃描，利用Mimics重建得到豬肝在重力下的切割形變模型，如圖5所示。

圖5 豬肝在重力下的切割形變模型Fig.5 Cutting deformation model of pork liver under gravity

沿與實際切割路徑相同的路徑對豬肝進行模擬切割，切割完畢導出Unity下的切割形變模型，如圖6所示。

圖6 Unity導出并實體化的切割模型Fig.6 Unity exports a materialized cutting model

為了更好地將真實模型與虛擬模型進行比對，將三角面片的殼體模型轉換為實體模型。該過程使用Geomagic Studio三維逆向軟件的精確曲面功能完成。軟件自動計算貼合網格模型的曲面，期間會生成少數劣質曲面，需要手動調整。分別將豬肝在平面上的重力形變網格模型與對應虛擬模型轉換為實體模型。

3 離體豬肝實驗結果討論

離體豬肝四面體剖分后頂點數量為6 510個，四面體單元為24 398個，三角單元為53 796個。將豬肝的平面重力形變實體模型和切割實體模型與對應的虛擬模型進行對比。模型對比使用Geomagic Quality的3D比較功能完成。設置真實形變模型為參考對象，虛擬模型為測試對象，分別對豬肝在平面上的重力形變模型和切割模型與相應的虛擬模型進行對比分析，并導出誤差報告文件。

圖7為豬肝在重力下的變形誤差云圖，圖中以真實形變模型為參考，誤差為正表示虛擬模型在真實模型外部，誤差為負表示虛擬模型在真實模型內部。可以看出整體誤差分布在-4～4 mm，由于底面為平面，誤差較小，約為-1 mm。正面誤差較大，模型前左半部分誤差為正，約2.5 mm，虛擬模型略厚于真實模型。模型底面誤差為負，虛擬模型的底面略高于真實模型。另外，模型邊緣部分誤差較大，即對肝臟邊緣的形變模擬較差。表1為肝臟重力形變的誤差分布表，圖8為誤差分布的柱狀圖，可以看到大部分誤差分布于-2～2.5 mm，占所有誤差的90%左右。

表1 重力變形誤差采樣點分布

圖7 豬肝在重力下的變形誤差Fig.7 Deformation error of pork liver under gravity

圖8 重力形變誤差分布Fig.8 Error distribution of gravity deformation

圖9為肝臟重力切割的誤差云圖，整體誤差與重力下的形變誤差相似，誤差主要分布于肝臟正面。圖中切口處可以看出，真實模型與虛擬模型的切口深度基本相同，切口越接近表面的部分誤差越大。切口上半部分誤差介于-3～-2 mm，誤差整體為負，說明虛擬切口的開口大于真實切口。切口下半部分誤差介于0～3 mm，為正說明虛擬切口的開口小于真實切口。表2為肝臟切割形變的誤差分布表，圖10為誤差分布的柱狀圖。由于切口部分所占比例較小，因此，切口處的誤差對整體誤差分布影響較小，誤差分布與重力形變的誤差分布基本相同。

圖9 豬肝在重力下切割的形變誤差Fig.9 Deformation error plot of pork liver cut under gravity

圖10 切割形變誤差分布Fig.10 Error distribution of cutting deformation

表2 切割形變誤差采樣點分布

通常計算機圖像的刷新率需要達到30 FPS才能感覺到流暢[18]，為了避免出現卡頓用以保證實時性，本文利用Unity中的Compute Shader分別對肝臟模型重力變形以及重力切割程序進行了加速實驗。圖11和圖12分別為經過加速后肝臟模型重力變形刷新率和重力切割程序刷新率，可以看出，刷新率分別達到了300 FPS和200 FPS，完全能夠滿足流暢性和實時性的要求。

圖11 重力變形程序的刷新率Fig.11 Refresh rate of gravity deformation program

圖12 重力切割程序的刷新率Fig.12 Refresh rate of gravity cutting program

4 結論

1)利用該模型獲得離體豬肝分別在重力作用下的變形及切割誤差分別在2 mm及3 mm以內，同時通過GPU加速，模型刷新率在300 FPS和200 FPS，能夠完成增強現實手術導航對軟組織模型實時性和準確性的要求。

2)在部分肝臟切除手術過程中肝臟受力狀況復雜，并不僅僅受到重力作用，還需考慮與手術器械之間的相互作用力，這也是后續研究需要展開的。