基于自適應(yīng)混合肝臟模型的建立

陳衛(wèi)東 趙國志 朱奇光*

1(燕山大學(xué)信息科學(xué)與工程學(xué)院，河北 秦皇島 066004)2(河北省特種光纖與光纖傳感重點實驗室，河北 秦皇島 066004)

基于自適應(yīng)混合肝臟模型的建立

陳衛(wèi)東1,2趙國志1,2朱奇光1,2*

1(燕山大學(xué)信息科學(xué)與工程學(xué)院，河北 秦皇島 066004)2(河北省特種光纖與光纖傳感重點實驗室，河北 秦皇島 066004)

基于黏彈性力學(xué)模型，建立一種自適應(yīng)無網(wǎng)格迦遼金與廣義表面網(wǎng)格混合模型。手術(shù)區(qū)采用自適應(yīng)無網(wǎng)格迦遼金法，非手術(shù)區(qū)采用廣義表面網(wǎng)格法，在手術(shù)區(qū)域和非手術(shù)區(qū)域之間建立過渡區(qū)域，過渡結(jié)點滿足力和位移的平衡，以實現(xiàn)兩個區(qū)域的無縫耦合。基于此模型，對人體肝臟進行形變仿真實驗。實驗表明，該模型在與無網(wǎng)格迦遼金模型形變狀況相同的前提下，單步執(zhí)行時間由17.4 ms降低到12.9 ms，幀頻由57.47幀/s提高到77.52幀/s。該模型繼承了無網(wǎng)格迦遼金模型的優(yōu)點，同時提高了模型的實時性。

虛擬手術(shù); 無網(wǎng)格迦遼金; 廣義表面網(wǎng)格; 形變仿真

引言

在基于虛擬現(xiàn)實的仿真系統(tǒng)中，人體軟組織形變模型的建立是虛擬手術(shù)仿真系統(tǒng)的關(guān)鍵技術(shù)。在虛擬環(huán)境力觸覺交互過程中，建立基于物理意義的力觸覺模型是非常重要的[1]。目前，國內(nèi)外學(xué)者所研究的軟組織模型大多是基于彈性理論的物理形變模型。Tagawa等提出了基于同步旋轉(zhuǎn)有限元的自適應(yīng)矩形四面體軟組織形變模型，該模型精確高，與傳統(tǒng)的有限元模型相比，實時性得到了提高，但要求網(wǎng)格的節(jié)點數(shù)仍很多、計算成本較高[2]。Peterlik等提出非線性有限元模型，雖然提高了實時性，仍無法滿足力觸覺反饋上千赫茲的頻率要求[3]。Wang等通過對模型的邊界進行離散簡化計算，但在穩(wěn)定性方面存在一定的困難[4]。Zhong等提出具有物理化學(xué)特性的反應(yīng)擴散表面模型，缺點是難于進行復(fù)雜的手術(shù)操作[5]。陳衛(wèi)東等基于彈簧質(zhì)點的結(jié)構(gòu)模型，建立一種改進的軟組織實時形變模型，在正六邊形拓撲結(jié)構(gòu)軟組織表面模型中，增加了虛擬體彈簧，以實現(xiàn)更加真實的形變效果，模型原理簡單，但是穩(wěn)定性差[6]。吳等提出了一種基于以受力點為中心、沿徑向方向分割為呈同心圓分布的彈簧-質(zhì)點系統(tǒng)，該方法的分割只涉及發(fā)生變形的局部區(qū)域,不需要對整個物體進行大量的矩陣運算，但是該模型各點計算不一致，不易擴展[7]。針對網(wǎng)格模型的缺點，徐少平等提出基于無網(wǎng)格的軟組織切割模型，打破了網(wǎng)格限制，在觸覺視覺逼真度和實時性方面都具有很大優(yōu)勢，但穩(wěn)定性較差[8]。Lennard-Jones勢能法，可以看作是彈簧質(zhì)點模型的一般化，其實時性較好，但是它的觸覺和視覺逼真度比較差[9]。離散機械力學(xué)法，缺乏描述塑性、黏彈性等的能力、觸覺和視覺效果一般[10]。光滑粒子流法，易于編程實現(xiàn)，但是計算量大，難于滿足實時性[11]。PCMFS方法，是一種利用連續(xù)介質(zhì)理論描述軟組織形變的方法，可以看成是基于質(zhì)點的離散系統(tǒng)和連續(xù)介質(zhì)的混合模型，該模型視覺和觸覺逼真度比較好，但是當節(jié)點數(shù)大于500時難于滿足實時性[12]。Muller法，實時性較好，由于應(yīng)變能函數(shù)中只考慮了線彈性材料特性，所以對軟組織性質(zhì)描述能力不足[13]。劉雪梅等基于黏彈性力學(xué)模型，提出了無網(wǎng)格迦遼金與彈簧質(zhì)點耦合的軟組織模型，為軟組織形變仿真提供了一種新的研究思路，高精度及與拓撲無關(guān)的無網(wǎng)格迦遼金法與高效簡單的彈簧質(zhì)點組合來實現(xiàn)軟組織形變仿真，可以有效地解決精度與效率的平衡問題，但是該模型的真實性和實時性還可以進一步提高[14]。

筆者針對目前軟組織模型存在的問題，基于黏彈性力學(xué)模型，提出了自適應(yīng)耦合無網(wǎng)格迦遼金與廣義表面網(wǎng)格混合模型、在非手術(shù)區(qū)廣義表面網(wǎng)格模型，有效地提高模型的實時性，在手術(shù)區(qū)采用自適應(yīng)無網(wǎng)格迦遼金法，使其不僅適合縫合、切割等復(fù)雜操作，而且進一步有效地提高模型的實時性，采用自適應(yīng)邊界，以實現(xiàn)兩個區(qū)域的無縫耦合。

1 模型構(gòu)建方法

軟組織器官具有松弛特性、不可壓縮特性、蠕變特性、各向異性等特性。由于其復(fù)雜性，軟組織模型完全真實地表現(xiàn)軟組織特性非常困難，難以滿足虛擬手術(shù)的實時性。下面所提出的自適應(yīng)耦合無網(wǎng)格迦遼金與廣義表面網(wǎng)格混合模型不僅適合縫合、切割等復(fù)雜操作情況，而且提高了軟組織模型的實時性。

1.1手術(shù)區(qū)域的自適應(yīng)算法

自適應(yīng)無網(wǎng)格迦遼金法的自適應(yīng)分析包括自適應(yīng)加密和誤差估計。自適應(yīng)加密需要確定稀疏和加密的區(qū)域，由于基于應(yīng)變梯度的方法不需要求解二階導(dǎo)數(shù)，因此本模型采用此方法。

求出節(jié)點的應(yīng)變、應(yīng)力和節(jié)點位移，應(yīng)變能為

(1)

式中，DSE是應(yīng)變能函數(shù)，σ為應(yīng)力，ε為位移。

對于三維情況，分別對x、y、z做偏導(dǎo)，應(yīng)變梯度(GSED)為

(2)

(3)

(4)

DM為網(wǎng)格密度，有

(5)

式中，N為區(qū)域內(nèi)節(jié)點數(shù)，V為區(qū)域體積。

應(yīng)變能變化體現(xiàn)為其梯度變化，在算法中采用背景網(wǎng)格進行誤差分析；對于第k個背景網(wǎng)格單元，首先要找到第k個背景網(wǎng)格內(nèi)的所有節(jié)點Nk，確定網(wǎng)格密度(DM)。

網(wǎng)格強度rd是自適應(yīng)方法的關(guān)鍵，將其定義為梯度對網(wǎng)格密度的比值，即

(6)

第k個單元的平均網(wǎng)格強度為

(7)

遍歷所有背景網(wǎng)格單元，獲得網(wǎng)格強度的最大值、最小值，有

Rmax=max(rd)Rmin=min(rd)

(8)

(9)

1.2非手術(shù)區(qū)域模型算法

利用非手術(shù)區(qū)域遠離接觸區(qū)域，沒有切割、縫合等復(fù)雜操作，根據(jù)添加剛體核彈簧質(zhì)點模型的思想來定義一種廣義的表面網(wǎng)格結(jié)構(gòu)，只在非手術(shù)區(qū)表層的1/5厚度具有質(zhì)點分布，每一個質(zhì)點對應(yīng)一個初始位置質(zhì)點，所有初始位置質(zhì)點不會隨著受力的改變而改變，在實際質(zhì)點與初始位置質(zhì)點之間引入一個初始長度為零的廣義彈簧，在非手術(shù)區(qū)域，實際質(zhì)點不但受到結(jié)構(gòu)彈簧的作用，還受到廣義彈簧的作用，廣義彈簧對質(zhì)點i的力為

fii=ksXii

(10)

式中，Xii是質(zhì)點i初始位置與當前位置的差值。

質(zhì)點i和質(zhì)點j之間的彈簧施加給質(zhì)點i的內(nèi)力為

(11)

式中，kd為彈簧的阻尼系數(shù)，Xij是兩個質(zhì)點的位置向量之差，ks是彈簧的勁度系數(shù)，lij是彈簧的原始長度，vj和vi是節(jié)點j和i的速度。

模型中質(zhì)點的運動都遵循牛頓第二定律，第i個質(zhì)點所受的合力fi等于所有彈簧提供的內(nèi)力和它所受外力之和。對于質(zhì)點i，有

(12)

式中，mi為質(zhì)點i的質(zhì)量，fi為質(zhì)點i所受的合力。

將每個質(zhì)點組合起來，可以得到模型的整體運動方程

X=M-1f(X,v)

(13)

式中，X、v和f分別是模型質(zhì)點的位置向量、速度向量和力向量，M為質(zhì)量矩陣。

1.3過渡區(qū)域的劃分

將手術(shù)區(qū)域和非手術(shù)區(qū)域作為兩個實體，通過過渡區(qū)域來實現(xiàn)手術(shù)區(qū)域和非手術(shù)區(qū)域的耦合，如圖1所示。根據(jù)手術(shù)區(qū)域所承受的載荷劃分手術(shù)區(qū)域和非手術(shù)區(qū)域的邊界，過渡結(jié)點需滿足位移及力的平衡條件。首先在原來網(wǎng)格的基礎(chǔ)上細分背景網(wǎng)格，細分后每個小的背景網(wǎng)格滿足的條件為至多有一個質(zhì)點存在，不包含質(zhì)點的背景網(wǎng)格稱為空網(wǎng)格，以空網(wǎng)格為對象，經(jīng)過順序迭代搜索，將符合無網(wǎng)格條件的轉(zhuǎn)化到手術(shù)區(qū)域，將符合廣義表面網(wǎng)格模型條件的轉(zhuǎn)換到非手術(shù)區(qū)域，然后將剩余的作為過渡單元。

圖1 過渡單元Fig.1 The schematic diagram of transition elements

非手術(shù)區(qū)域內(nèi)部轉(zhuǎn)換條件：如果6個方向都存在單元，有4個及以上的單元為彈簧質(zhì)點單元時，將此單元轉(zhuǎn)換為彈簧質(zhì)點單元；如果6個方向都存在單元且都為空單元，將此單元轉(zhuǎn)換為彈簧質(zhì)點單元；如果6個方向不完全存在，有3個及以上的單元為彈簧質(zhì)點單元時，將此單元轉(zhuǎn)換為彈簧質(zhì)點單元。

手術(shù)區(qū)域與非手術(shù)區(qū)域單元轉(zhuǎn)換條件相同。

剩余空單元的特征為各個方向同時包含彈簧質(zhì)點單元和無網(wǎng)格單元。確定過渡單元后，過渡質(zhì)點需滿足所受合力為0，位移相等。遍歷過渡區(qū)域的所有質(zhì)點，得到手術(shù)區(qū)域和非手術(shù)區(qū)域模型的近似位移函數(shù)，計算出二者之間的函數(shù)關(guān)系作為過渡區(qū)域質(zhì)點的位移函數(shù)，實現(xiàn)兩個區(qū)域的光滑過渡。

2 模型實驗

本研究的仿真以肝臟模型為例，其泊松比設(shè)為0.35，楊氏模量設(shè)為5 000 Pa，平均密度為1.055 g/cm3。成年女性肝臟為800～1 200 g，成年男性肝臟為1 200～1 500 g。仿真基于Intel 酷睿i3 4130處理器，主頻為3.40 GHz，內(nèi)存為4.0 GB，顯卡為NVIDIA GeForce GTX 650，系統(tǒng)為WindowS 7系統(tǒng)，開發(fā)環(huán)境為Visual C++6.0，載荷為5.7 N。

用自適應(yīng)方法、無網(wǎng)格方法和耦合方法[14]分別構(gòu)建肝臟模型，在相同的質(zhì)點施加相同的拉力或壓力，分別計算其位移，將3個模型的位移曲線進行對比。

在實驗中，模型的質(zhì)點數(shù)同為295時，施力步長為0.3 N，分別對自適應(yīng)模型、無網(wǎng)格模型和耦合模型[14]進行平行對比實驗。

3 結(jié)果

在虛擬手術(shù)仿真系統(tǒng)中，視覺逼真度是軟組織仿真系統(tǒng)的一個重要指標，是軟組織仿真系統(tǒng)沉浸感的基礎(chǔ)。視覺逼真度與模型的描述能力有關(guān)，描述能力越強，模型的逼真度越高。模型的視覺逼真度用形變渲染圖來體現(xiàn)。圖2(a)為肝臟受拉力時的局部剖面圖，虛線內(nèi)部為手術(shù)區(qū)域，虛線外部為非手術(shù)區(qū)域，虛線區(qū)域為過渡區(qū)域。圖3為肝臟受拉力時的形變渲染截圖，施力器械為鑷子，施力點為圓點處，細線處為過渡區(qū)域。圖2(b)為肝臟受壓力時的局部剖面圖，虛線內(nèi)部為手術(shù)區(qū)域，虛線外部為非手術(shù)區(qū)域，虛線區(qū)域為過渡區(qū)域。圖4為肝臟受壓力時的形變渲染圖，施力器械為鑷子，施力點為圓點處，細線處為過渡區(qū)域。由圖2～4可以看出，過渡區(qū)域沒有出現(xiàn)尖角和斷裂的痕跡，在形變時是連續(xù)光滑的，直觀證明了筆者所提出的自適應(yīng)模型是有效的。

圖2 力作用下的肝臟局部剖面圖。(a) 拉力；(b)壓力Fig.2 The local liver profile under the action of force (a) Pull；(b) Press

圖3 拉力作用下的肝臟渲染圖Fig.3 The liver rendering under the action of pull force

圖4 壓力作用下的肝臟渲染圖Fig.4 The liver rendering under the action of pressure forc

虛擬手術(shù)系統(tǒng)為受訓(xùn)醫(yī)生在視覺和觸覺上提供各種手術(shù)場景的真實再現(xiàn)，體驗和學(xué)習(xí)各種手術(shù)情況，降低手術(shù)中的風(fēng)險，降低醫(yī)生培訓(xùn)成本，對提高醫(yī)療水平有非常重要的意義。虛擬手術(shù)系統(tǒng)的核心模塊是軟組織模型，實時性和真實性是軟組模型的關(guān)鍵指標。

圖5為壓力作用下的曲線，圖6為拉力作用下的曲線。由圖5 和圖6可以看出，在相同壓力或拉力作用下，自適應(yīng)模型的位移與無網(wǎng)格模型的誤差為±0.01 cm，而耦合模型與無網(wǎng)格模型的誤差為±0.018 cm。從數(shù)值上表明，自適應(yīng)模型的真實性比耦合模型的真實性要高。

圖5 壓力作用下的位移曲線Fig.5 The displacement curve under the action of pressure force

圖6 拉力作用下的位移曲線Fig.6 The displacement curve under the action of pull force

軟組織模型的實時交互能力對虛擬手術(shù)仿真系統(tǒng)具有至關(guān)重要的意義。實時交互能力直接決定了模型的仿真效果，是軟組織模型最為關(guān)鍵的性能指標。模型的實時性與模型的單步執(zhí)行時間以及幀頻有關(guān)。

表1為自適應(yīng)模型與迦遼金無網(wǎng)格模型的效率對比。由表1可以看出，筆者所提出的自適應(yīng)耦合無網(wǎng)格迦遼金與廣義表面網(wǎng)格混合模型的單步執(zhí)行時間比無網(wǎng)格模型少了4.5 ms，幀頻提高了20.048幀，這意味著模型的實時性得到了進一步的提高。

表1自適應(yīng)模型與迦遼金無網(wǎng)格模型的效率對比

Table1Theefficiencycontrastofthedoubleadaptivecouplingmodelandtheelement-freeGalerkin

指標無網(wǎng)格模型耦合模型自適應(yīng)模型單步執(zhí)行時間/ms17.414.812.9形變頻率/(幀/s)57.4767.5777.52

4 討論和結(jié)論

本研究所提出的自適應(yīng)耦合無網(wǎng)格迦遼金與廣義表面網(wǎng)格混合模型在真實性和實時性上對軟組織模型有一定提高，但是相對于軟組織的生物力學(xué)特性還有差距，如在各向異性、非線性等方面，因此還需進一步探求能夠?qū)浗M織模型定量準確描述的模型。手術(shù)的要求非常嚴格，因此對于試驗算法，還需要更為完善的驗證方案，如將力反饋、時間、范圍等也引入到驗證方案中，以達到準確驗證的目地。

本研究針對目前軟組織模型存在的問題，基于自適應(yīng)耦合無網(wǎng)格迦遼金與廣義表面網(wǎng)格方法，建立了肝臟模型，在直觀視覺和數(shù)值兩方面驗證了混合模型的有效性，在保證精度的前提下進一步提高了效率。在手術(shù)區(qū)域采用自適應(yīng)無網(wǎng)格迦遼金法，既保證了手術(shù)區(qū)域的精度，又提高了計算效率，能夠滿足拓撲結(jié)構(gòu)改變的情況，比有限元更加適合大縫合、切割等復(fù)雜形變情況；非手術(shù)區(qū)域一般不會出現(xiàn)切割、縫合等復(fù)雜操作，因此在非手術(shù)區(qū)域采用廣義的表面網(wǎng)格模型，可進一步提高模型的計算效率。實驗結(jié)果表明，基于自適應(yīng)耦合無網(wǎng)格迦遼金與廣義無網(wǎng)格建立的肝臟模型，能夠滿足虛擬手術(shù)實時性、穩(wěn)定性、真實性的要求。

[1] 馮元楨,著. 戴克剛, 鞠烽熾, 譯. 生物力學(xué): 活組織的力學(xué)特征[M]. 長沙: 湖南科學(xué)技術(shù)出版社, 1986: 222-296.

[2] Tagawa K, Yamada T, Tanakal HT. A rectangular tetrahedral adaptive mesh based corotated finite element model for interactive soft tissue simulation[C] //The 35th Annual International Conference of the IEEE EMBS. Osaka: IEEE, 2013:7164-7167.

[3] Igor P, Mert S, Cagatay B,etal. Real-time visio-haptic interaction with static soft tissue models having geometric and material nonlinearity[J]. Computers & Graphics, 2010,34: 43-54.

[4] Wang P, Becker AA, Jones IA,etal. Virtual reality simulation of surgery with haptic feedback based on the boundary element method, Computers and Structures[J]. 2007,85: 331-339.

[5] Zhong Yongmin, Bijan S, Julian S,etal. Soft tissue deformation with reaction-diffusion process for surgery simulation[J]. Journal of Visual Languages & Computing, 2012,23: 1-12.

[6] 陳衛(wèi)東, 趙成龍, 朱奇光, 等. 虛擬手術(shù)中軟組織形變建模及力反饋算法研究[J]. 中國生物醫(yī)學(xué)工程學(xué)報, 2013,32(1): 114-123.

[7] 吳涓, 宋愛國, 李建清. 用于虛擬力覺再現(xiàn)的彈簧-質(zhì)點力/變形模型研究[J].系統(tǒng)仿真學(xué)報, 2006,18(11): 3152-3156.

[8] 徐少平, 李春泉, 江順亮, 等. 基于無網(wǎng)格的軟組織切割模型的研究進展[J]. 中國生物醫(yī)學(xué)工程學(xué)報, 2012,31(1): 123-128.

[9] Nealen A, Mtiller M, Keiser R,etal. Physically based deformable models in computer graphics[J]. Comput. Graph. Forum 2006,25(4): 809-836.

[10] Jansson J, Vergeest JSM. A discrete mechanics model for deformable bodies[J]. Computer-Aided Design, 2002,34(12): 913-928.

[11] Hieber SE, Koumoutsakos P. A lagrangian particle method for the simulation of linear and nonlinear elastic models of soft tissue [J]. Journal of Computational Physics, 2008,227(21): 9195-9215.

[12] Lim Y, Suvranu D. Real time simulation of nonlinear tissue response in virtual surgery using the point collocation-based method of finite spheres[J]. Comput Methods Appl Mech. Engrg, 2007,196(31): 3011-3024.

[13] Mtiller M, Heidelberger B, Teschner M,etal. Meshless deformations based on shape matching[J]. ACM Trans Graph, 2005,24(3): 471-478.

[14] 劉雪梅, 毛磊, 李運華, 等. 耦合無網(wǎng)格迦遼金與質(zhì)點彈簧實現(xiàn)軟組織形 變仿真[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2013,25(1): 1-6.

The Construction of a Liver Model Baseed on Adaptive Hybriding Method

Chen Weidong1,2Zhao Guozhi1,2,Zhu Qiguang1,2*

1(Insitute of Information Science and Engineering, Yanshan University, Qinhuangdao 066004, Hebei, China)2(The Key Laboratory for Special Fiber and Fiber Sensor of Hebei Province, Qinhuangdao 066004, Hebei, China)

virtual surgery; element-free Galerkin; generalized surface mesh; deformation simulation

10.3969/j.issn.0258-8021. 2015. 03.014

2014-10-17， 錄用日期:2014-12-28

國家自然科學(xué)基金(61201112,61172044)；河北省自然科學(xué)基金(F2013203250, F2012203169)

R318

D

0258-8021(2015) 03-0360-05

*通信作者(Corresponding author)， E-mail:zhu7880@ysu.edu.cn