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

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

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

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

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

1 質(zhì)點(diǎn)彈簧軟組織模型

1.1 基于Tetgen的四面體剖分

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

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

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

網(wǎng)格質(zhì)量控制部分包括一個(gè)保留表面網(wǎng)格的選項(xiàng)和一個(gè)質(zhì)量控制滑塊。當(dāng)選擇保留表面網(wǎng)格時(shí)，程序在進(jìn)行四面體剖分時(shí)將保留原模型的表面網(wǎng)格，只對(duì)模型內(nèi)部進(jìn)行操作。否則，程序在進(jìn)行四面體剖分時(shí)將同時(shí)對(duì)表面網(wǎng)格進(jìn)行重構(gòu)。質(zhì)量控制滑塊能夠?qū)δＰ蛢?nèi)部四面體單元的質(zhì)量進(jìn)行控制。四面體單元的質(zhì)量由2種方法進(jìn)行控制：半徑-邊緣比約束和最小二面角約束。四面體T的半徑-邊緣比ρ(T)為：

(1)

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

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

使用上述四面體剖分程序?qū)Ω闻K表面網(wǎng)格模型進(jìn)行四面體剖分，輸出結(jié)果如圖1所示。

圖1 四面體剖分結(jié)果Fig.1 Tetrahedral meshing results

1.2 質(zhì)點(diǎn)彈簧力學(xué)模型

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

圖2 質(zhì)點(diǎn)彈簧模型示意Fig.2 Schematic diagram of the mass point spring model

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

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

(2)

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

質(zhì)點(diǎn)i受到的外力由手術(shù)器械的交互力fc(i)和重力fg(i)組成：

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

(3)

Voigt元件產(chǎn)生的力由阻尼力fdmp和彈簧力fspr組成：

(4)

式中：v為質(zhì)點(diǎn)的移動(dòng)速度，r為質(zhì)點(diǎn)的空間坐標(biāo)，nij為連接質(zhì)點(diǎn)i和質(zhì)點(diǎn)j的Voigt元件的阻尼系數(shù)，kij為連接質(zhì)點(diǎn)i和質(zhì)點(diǎn)j的Voigt元件的剛度系數(shù)。

質(zhì)點(diǎn)i受到的合外力由牛頓第二定律得出：

fsum(i)=miai

(5)

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

(6)

計(jì)算時(shí)對(duì)時(shí)間進(jìn)行離散化處理，每隔Δt=0.02 s進(jìn)行一次計(jì)算。計(jì)算質(zhì)點(diǎn)的加速度a，則Δt后質(zhì)點(diǎn)的速度為：

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

(7)

質(zhì)點(diǎn)在Δt內(nèi)的運(yùn)動(dòng)可視為勻速運(yùn)動(dòng)，因此，Δt后質(zhì)點(diǎn)的坐標(biāo)為：

ri+1=ri+vi·Δt

(8)

在Δt時(shí)間內(nèi)計(jì)算所有質(zhì)點(diǎn)的位置坐標(biāo)即可得出0.02 s后肝臟模型的整體形態(tài)，當(dāng)所有頂點(diǎn)的受力達(dá)到平衡時(shí)，肝臟模型不再變形，此時(shí)即為肝臟模型在重力作用下的形態(tài)。

2 離體豬肝切割及變形實(shí)驗(yàn)

本文使用質(zhì)點(diǎn)彈簧模型模擬肝臟在平面上受重力影響的變形和在重力影響下的切割變形，本節(jié)以豬肝為研究對(duì)象，分別對(duì)以上2種情況進(jìn)行實(shí)驗(yàn)驗(yàn)證。

2.1 肝臟模型在重力下的變形實(shí)驗(yàn)

要對(duì)肝臟模型在重力影響下的變形進(jìn)行驗(yàn)證，必須得知肝臟在失重狀態(tài)下的原始形態(tài)并獲得其失重模型。因此需要模擬一個(gè)失重的環(huán)境，將豬肝置于其中，獲取豬肝在失重狀態(tài)下的三維數(shù)據(jù)。

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

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

圖3 懸浮狀態(tài)下的豬肝模型Fig.3 Pork liver model in suspension

CT掃描完成后將豬肝從鹽水中取出放置在平面上再次掃描，重建得到豬肝在重力下的形變模型作為對(duì)比標(biāo)準(zhǔn)。將豬肝的懸浮模型導(dǎo)入U(xiǎn)nity并賦予程序腳本，新建一個(gè)平面，將豬肝模型置于略高于平面的上方，運(yùn)行程序模擬豬肝在重力下的變形。待豬肝完全落地并穩(wěn)定之后導(dǎo)出變形后的模型，如圖4所示。

圖4 Unity導(dǎo)出并實(shí)體化的重力形變模型Fig.4 Gravity deformation model exported and solidified by Unity

2.2 肝臟模型在重力下的切割實(shí)驗(yàn)

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

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

沿與實(shí)際切割路徑相同的路徑對(duì)豬肝進(jìn)行模擬切割，切割完畢導(dǎo)出Unity下的切割形變模型，如圖6所示。

圖6 Unity導(dǎo)出并實(shí)體化的切割模型Fig.6 Unity exports a materialized cutting model

為了更好地將真實(shí)模型與虛擬模型進(jìn)行比對(duì)，將三角面片的殼體模型轉(zhuǎn)換為實(shí)體模型。該過程使用Geomagic Studio三維逆向軟件的精確曲面功能完成。軟件自動(dòng)計(jì)算貼合網(wǎng)格模型的曲面，期間會(huì)生成少數(shù)劣質(zhì)曲面，需要手動(dòng)調(diào)整。分別將豬肝在平面上的重力形變網(wǎng)格模型與對(duì)應(yīng)虛擬模型轉(zhuǎn)換為實(shí)體模型。

3 離體豬肝實(shí)驗(yàn)結(jié)果討論

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

圖7為豬肝在重力下的變形誤差云圖，圖中以真實(shí)形變模型為參考，誤差為正表示虛擬模型在真實(shí)模型外部，誤差為負(fù)表示虛擬模型在真實(shí)模型內(nèi)部?？梢钥闯稣w誤差分布在-4～4 mm，由于底面為平面，誤差較小，約為-1 mm。正面誤差較大，模型前左半部分誤差為正，約2.5 mm，虛擬模型略厚于真實(shí)模型。模型底面誤差為負(fù)，虛擬模型的底面略高于真實(shí)模型。另外，模型邊緣部分誤差較大，即對(duì)肝臟邊緣的形變模擬較差。表1為肝臟重力形變的誤差分布表，圖8為誤差分布的柱狀圖，可以看到大部分誤差分布于-2～2.5 mm，占所有誤差的90%左右。

表1 重力變形誤差采樣點(diǎn)分布

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

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

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

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

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

表2 切割形變誤差采樣點(diǎn)分布

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

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

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

4 結(jié)論

1)利用該模型獲得離體豬肝分別在重力作用下的變形及切割誤差分別在2 mm及3 mm以內(nèi)，同時(shí)通過GPU加速，模型刷新率在300 FPS和200 FPS，能夠完成增強(qiáng)現(xiàn)實(shí)手術(shù)導(dǎo)航對(duì)軟組織模型實(shí)時(shí)性和準(zhǔn)確性的要求。

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