羅 丹,張 禹,李 佳

(沈陽(yáng)工業(yè)大學(xué) 機(jī)械工程學(xué)院,沈陽(yáng) 110870)

1 引 言

虛擬手術(shù)就是依托虛擬現(xiàn)實(shí)技術(shù)和信息技術(shù),在計(jì)算機(jī)中建立一個(gè)虛擬環(huán)境[1],依靠醫(yī)學(xué)影像數(shù)據(jù)和虛擬環(huán)境信息,可以進(jìn)行醫(yī)學(xué)培訓(xùn)、手術(shù)模擬、術(shù)中引導(dǎo)、術(shù)后康復(fù)和遠(yuǎn)程醫(yī)療.虛擬手術(shù)融合了計(jì)算機(jī)技術(shù)、生物學(xué)、現(xiàn)代醫(yī)學(xué)、機(jī)器人學(xué)等多學(xué)科領(lǐng)域,主要由醫(yī)學(xué)數(shù)據(jù)處理、人體組織器官(主要是軟組織)建模、形變及切割仿真、力反饋4部分構(gòu)成,為用戶提供了視覺(jué)和觸覺(jué)的交互感知[2].軟組織的物理建模是虛擬手術(shù)研究的重點(diǎn),物理模型是指人體軟組織模型在與手術(shù)器械交互過(guò)程能夠表現(xiàn)出來(lái)的視覺(jué)和觸覺(jué)信息(包括形變位移和力反饋).它首先通過(guò)對(duì)醫(yī)學(xué)圖像的處理與分割獲得人體生物組織的幾何模型,然后利用對(duì)人體組織生物力學(xué)特性的研究,根據(jù)軟組織在受力的情況下,產(chǎn)生的形變過(guò)程和計(jì)算的反饋力來(lái)進(jìn)行物理建模[3].合理地建立虛擬手術(shù)仿真模型,既能體現(xiàn)真實(shí)性和沉浸感,又能保證虛擬仿真的實(shí)時(shí)性.

由于人體軟組織具有不均勻性、各向異性、非線性和粘彈性等性質(zhì)[4,5],其形變過(guò)程中對(duì)于力與形變的關(guān)系也表現(xiàn)出相應(yīng)的復(fù)雜性.因此,為了保證形變模型的準(zhǔn)確性,更好的將可視化與觸覺(jué)反饋結(jié)合,在虛擬手術(shù)中軟組織的物理建模中必須考慮這些生物力學(xué)特性.目前,描述軟組織變形的物理模型中最簡(jiǎn)單的是根據(jù)胡克定理的線彈性模型(linear elastic model),但這個(gè)模型一般只適用于位置移動(dòng)較小的情況,并且線性模型不足以表達(dá)復(fù)雜的運(yùn)動(dòng)情況,常用非線性模型來(lái)進(jìn)行表述:例如線彈性模型、粘彈性模型、簡(jiǎn)單的Newton 粘性流公式等[6,7].大多數(shù)的人體組織材料既有彈性又有粘性,因此模型的建立必須基于生物力學(xué)特性的基礎(chǔ)上,兼顧實(shí)時(shí)性.在真實(shí)性和實(shí)時(shí)性之間折衷來(lái)模擬軟組織形變.

軟組織物理建模技術(shù)即軟體的變形模擬技術(shù),國(guó)內(nèi)外學(xué)者做了大量的研究.基于物理的方法中最常用的是質(zhì)點(diǎn)彈簧法(Mass Spring Method,MSM).Terzopoulos最早了提出質(zhì)點(diǎn)彈簧模型,對(duì)皮膚的形變過(guò)程進(jìn)行了仿真[8].質(zhì)點(diǎn)彈簧法是將連續(xù)的組織離散化為一系列有質(zhì)量的點(diǎn),質(zhì)點(diǎn)之間通過(guò)彈簧來(lái)進(jìn)行連接.質(zhì)點(diǎn)受到外力、彈力、阻尼力的作用,通過(guò)牛頓第二定律來(lái)計(jì)算其變形[9,10].Waters K.等人使用質(zhì)點(diǎn)彈簧法對(duì)人臉進(jìn)行了仿真[11].國(guó)內(nèi)鮑春波等人在2006年對(duì)質(zhì)點(diǎn)彈簧法進(jìn)行改進(jìn)[12],按照質(zhì)點(diǎn)受力大小進(jìn)行動(dòng)態(tài)細(xì)分,實(shí)現(xiàn)了仿真精度的提高.這種方法是根據(jù)牛頓第二定律計(jì)算出加速度的內(nèi)部和外部力,然后根據(jù)時(shí)間積分方法來(lái)更新速度,最后更新對(duì)象的位置.MSM結(jié)構(gòu)簡(jiǎn)單,計(jì)算成本低,因此在軟組織建模中應(yīng)用廣泛,但在形變量較大時(shí)振蕩效果明顯,容易產(chǎn)生塌陷等變形缺陷,不易表現(xiàn)真實(shí)的形變特性.

基于物理的質(zhì)點(diǎn)彈簧方法(Mass Spring Method,MSM),在呈現(xiàn)出時(shí)域相關(guān)的非線性及粘彈性的力學(xué)特征下,有著局部變形大且失真的弱點(diǎn);而基于幾何的PBD模型很適合應(yīng)用到以實(shí)時(shí)性和手術(shù)操作為主要目標(biāo)的虛擬手術(shù)模擬器中,雖然可以采用不同的約束函數(shù)來(lái)模擬變形特征,且魯棒性好,卻很難表現(xiàn)出時(shí)域特性;XPBD對(duì)不同材質(zhì)的影響不明顯,精度表現(xiàn)一般.根據(jù)MSM和PBD都是通過(guò)空間離散質(zhì)點(diǎn)的運(yùn)動(dòng)位置來(lái)模擬形變特征的共通點(diǎn),本文提出了一種將MSM與PBD相結(jié)合的模型,即融合生物力學(xué)特性的位置動(dòng)力學(xué)模型(Bio-based PBD,B-PBD).兩種模型的整合思路也來(lái)源于XPBD的推導(dǎo)過(guò)程:將非線性、粘彈性被引入到彈力計(jì)算中,由牛頓第二定律根據(jù)彈力及約束力的合力計(jì)算出質(zhì)點(diǎn)位置.模型分為兩個(gè)部分:粘彈性模型反應(yīng)軟組織的生物力學(xué)特性;PBD處理復(fù)雜的效果,如不可壓縮性,并保持模型的穩(wěn)定性.實(shí)現(xiàn)了快速,穩(wěn)定的模擬及復(fù)雜變形的仿真.

2 融合生物特性的位置動(dòng)力學(xué)建模(Bio-based PBD,B-PBD)

2.1 模型建立

本文的物理模型采用四面體網(wǎng)格,基于位置的動(dòng)力學(xué)算法用到兩個(gè)基本條件:即四面體的頂點(diǎn)和約束.假設(shè)四面體網(wǎng)格中頂點(diǎn)數(shù)為N,四面體每一個(gè)頂點(diǎn)的信息包括:質(zhì)量mi、位置Pi、速度Vi,模型中約束條件數(shù)目為M,第j個(gè)約束條件的基數(shù)為nj,每個(gè)約束可以表示成Cj(xi1,…xinj)=0,則j滿足類型相等的約束.

PBD的核心分為兩步:

1)通過(guò)外力的顯式步進(jìn)求到一個(gè)暫時(shí)的形變位置;

2)定義一系列的約束.通過(guò)應(yīng)用一系列的約束來(lái)模擬不同的形變效果,來(lái)改變這個(gè)暫時(shí)的形變位置得到最終的形變位置.

PBD單元模型2在頂點(diǎn)xi之間設(shè)置拉伸約束、在兩個(gè)三角面片之間設(shè)置彎曲約束、兩個(gè)四面體之間設(shè)置體積約束.四面體頂點(diǎn)的位置改變時(shí),約束函數(shù)將每個(gè)點(diǎn)投影到合適的位置,實(shí)現(xiàn)了由位置計(jì)算變形.B-PBD單元模型如圖1(a)所示.在原有PBD約束的基礎(chǔ)上,在頂點(diǎn)xi之間增加了三參量質(zhì)點(diǎn)彈簧模型約束,表達(dá)模型的生物力學(xué)特性.三參量模型如圖1(b)所示.在其頂點(diǎn)處布置具有質(zhì)量、位置、半徑、轉(zhuǎn)動(dòng)慣量等屬性的填充球體,填充球之間采用三參量粘彈性結(jié)構(gòu)進(jìn)行連接.當(dāng)軟組織器官表面受到力的作用后,受力的填充球會(huì)通過(guò)四面體的邊帶動(dòng)與之相連的填充球運(yùn)動(dòng),模擬軟組織變形.在模型外覆蓋面模型,通過(guò)添加紋理使模型更加逼真,滿足視覺(jué)要求.

圖1 B-PBD模型Fig.1 Model of B-PBD

2.2 約束投影

位置動(dòng)力學(xué)使用約束條件來(lái)控制物體的變形,這些約束條件一般基于幾何規(guī)則,根據(jù)幾何規(guī)定定義的約束函數(shù)[16].約束函數(shù)用于確定動(dòng)態(tài)物體的力學(xué)特性.約束投影即講修正質(zhì)點(diǎn)運(yùn)動(dòng)來(lái)滿足約束函數(shù)的過(guò)程.

計(jì)算過(guò)程如下:

1)對(duì)四面體每個(gè)頂點(diǎn)的位置Pi和速度Vi進(jìn)行初始化;

(1)

式中:F為外力;ωi=1/mi,mi為頂點(diǎn)質(zhì)量;Damp(*)為阻尼函數(shù).

根據(jù)外力預(yù)測(cè)下一時(shí)間間隔內(nèi)頂點(diǎn)的位置:

(2)

(3)

(4)

5)循環(huán)計(jì)算,返回第2)步,更新每個(gè)頂點(diǎn)的速度,預(yù)測(cè)下一時(shí)刻的位置.對(duì)于每個(gè)頂點(diǎn)的位置改變量ΔP:

(5)

式中,C(P1,P2,…,Pn)為頂點(diǎn)位置的約束函數(shù);▽PiC(P1,P2,…,Pn)為函數(shù)關(guān)于頂點(diǎn)Pi的梯度.

在模擬軟組織變形中,常用的約束有拉伸約束、彎曲約束、體積約束、能量保持約束等.拉伸約束是指當(dāng)兩個(gè)填充球在外力的作用下距離發(fā)生變化時(shí),拉伸約束使其回到初始位置.對(duì)于拉伸約束可得:

C(P1,P2)=|P1-P2|-d

(6)

其中,P1和P2為兩個(gè)頂點(diǎn),d是兩個(gè)頂點(diǎn)之間的靜止距離,對(duì)P1和P2分別求梯度,并代入到公式(5)中,可得:

(7)

(8)

體積保持約束用于防止出現(xiàn)“空心”的效果,在外力的作用下四面體的體積不發(fā)生變化.體積約束的函數(shù)表達(dá)式為:

(9)

其中,V0是四面體的體積,P1,P2,P3,P4是四面體的4個(gè)頂點(diǎn),對(duì)4個(gè)頂點(diǎn)求梯度,并代入公式(5)中,得到:

(10)

2.3 算法推導(dǎo)

模型對(duì)象由一組質(zhì)點(diǎn),約束和彈簧、阻尼組成.質(zhì)點(diǎn)的位置由其總力負(fù)荷控制.計(jì)算的目標(biāo)是找出每個(gè)步驟中每個(gè)質(zhì)點(diǎn)的實(shí)際位置.基于牛頓第二定律:

(11)

其中x=[x1x2…xn]T是質(zhì)點(diǎn)位置的列向量,是質(zhì)點(diǎn)的加速度.M為質(zhì)量的對(duì)角矩陣,FC、FM和FE分別是勢(shì)能約束力、彈簧和外力的列向量.

約束能量勢(shì)被定義為:

FC=-▽UT(x)

(12)

其中U(x)為能量勢(shì)能.

將方程(11)離散化,上標(biāo)n表示時(shí)間步長(zhǎng)指數(shù),則有:

(13)

(14)

其中C(x)=[C1C2…Cm]T為列的約束函數(shù),α是對(duì)角柔度矩陣,相當(dāng)于對(duì)應(yīng)的約束逆剛度矩陣.則來(lái)自彈性勢(shì)能的力由U的負(fù)梯度:

FC=-▽xUT=-CTα-1C

(15)

(16)

可得:

(17)

(18)

(19)

將方程(18)、(19)分別表示為g和h,則有:

g(x,λ)=0

(20)

h(x,λ)=0

(21)

線性化方程可得:

(22)

其中K=?g/?x.可以針對(duì)Δx和Δλ求解該系統(tǒng),并相應(yīng)地更新乘子和位置:

λi+1=λi+Δλ

(23)

xi+1=xi+Δx

(24)

(25)

得到位置更新:

Δx=M-1▽C(xi)TΔλ

(26)

拉格朗日乘子的修正量:

(27)

在約束求解期間,首先計(jì)算單個(gè)約束的Δλ,然后使用等式(23)、式(24)更新系數(shù)和乘子.

軟組織的粘彈性通過(guò)三參量質(zhì)點(diǎn)彈簧(MSM)模型來(lái)體現(xiàn).該模型的彈性特性通過(guò)彈簧結(jié)構(gòu)來(lái)表示,粘性特性通過(guò)彈簧和阻尼串聯(lián)結(jié)構(gòu)來(lái)表示,如圖1(b)所示.軟組織的形變量即為系統(tǒng)疊加的位移總和.

由圖1(b)可知:

Fn(t)=KnXn(t)

(28)

(29)

Ff(t)=Fg(t)+Fn(t)

(30)

式中,Kn,Kg是彈簧彈性系數(shù),Cg是粘性系數(shù),Fn(t)、Fg(t)、Ff(t)分別是外力,Xn(t)、Xg1(t)、Xg2(t)是由外力引起的位移.

經(jīng)過(guò)三維重建后的面模型,使用Tetgen算法對(duì)四面體網(wǎng)格進(jìn)行劃分.以正四面體為例,如圖2所示.在正四面體OABC中,O,A,B,C分別代表每個(gè)填充球體的球心,在力F(t)的作用下,若只考慮拉伸力,忽略彎曲和扭轉(zhuǎn)的條件下,則球心O移動(dòng)到了O′點(diǎn),OB粘彈性結(jié)構(gòu)長(zhǎng)度發(fā)生了變化,與之相連接的OA,OC長(zhǎng)度也相應(yīng)發(fā)生改變.由于變形量非常小,忽略球心A,B,C球的移動(dòng)量.設(shè)形變前后OB長(zhǎng)度變化量OO′為X(t),CO′與O′B夾角為θ(60°≤θ≤90°),OC長(zhǎng)度設(shè)為y,形變后O′C長(zhǎng)度為x,則變化值ΔX(t)=y-x.

圖2 小變形受力分析Fig.2 Stress analysis of small deformation

根據(jù)幾何關(guān)系得出:

(31)

(32)

由胡克定律得出:

Fk(t)=2K′ΔX(t)cosθ

(33)

令λ滿足:

(34)

則公式(33)可以寫成:

F(t)=FK(t)+Ff(t)

(35)

其中,K′是三參量粘彈結(jié)構(gòu)的彈性系數(shù),Fk(t)是三參量模型受力分量.

根據(jù)小變形理論,則有:

F(t)=FK(t)+Ff(t)

(36)

根據(jù)上述的小變形理論受力方程(35)及(36),建立混合填充球模型的力位移方程,拉氏變換可得:

(37)

其中τ是松弛時(shí)間,滿足:

(38)

1)力-時(shí)間方程建立

(39)

當(dāng)t→0+時(shí),

F(0+)=(λ+Kn+Kg)X0

(40)

當(dāng)t→∞時(shí),

F(∞)=(λ+Kn)X0

(41)

2)位移-時(shí)間方程建立

(42)

其中:

(43)

當(dāng)t→0+時(shí):

(44)

當(dāng)t→∞時(shí):

(45)

每個(gè)質(zhì)點(diǎn)的彈簧力:

FM=Fd+F

(46)

其中:

(47)

b0和b1為兩個(gè)阻尼約束[17],Fd為直接作用在速度上的阻尼,F為三參量模型的總彈性力.

2.4 算法流程

根據(jù)上面的推導(dǎo),可以總結(jié)出算法的流程如表1所示.

表1 算法流程Table 1 Algorithm flow of B-PBD

基于生物力學(xué)特性的位置動(dòng)力學(xué)方法的概念是將行為模型分成兩個(gè)部分:用XPBD以適應(yīng)拓?fù)湫薷暮突A(chǔ)組織的材料非均質(zhì)性,處理復(fù)雜的效果,并保持模型的穩(wěn)定性,變形效果取決于迭代器計(jì)數(shù);三參量MSM模型反映軟組織的物理特性,非線性和粘彈性,并將擴(kuò)展彈簧力作為外力引入以估算質(zhì)點(diǎn)位置,保持了最終質(zhì)點(diǎn)位置的穩(wěn)定性.

與PBD和XPBD相比,B-PBD在預(yù)測(cè)位置步驟中加入了反應(yīng)軟組織生物力學(xué)特性的改進(jìn)的三參量模型.首先在彈簧彈力計(jì)算過(guò)程中體現(xiàn)了非線性與粘彈性,基于牛頓第二定律,在進(jìn)行下一時(shí)刻位置預(yù)測(cè)時(shí),將彈簧力作為外力,采用顯式歐拉法求解.基于PBD的體積不變約束對(duì)質(zhì)點(diǎn)位置進(jìn)行修正.在約束投影步驟中,按照預(yù)先給定的約束的順序分別求出梯度的數(shù)值和拉格朗日乘子的數(shù)值,并更新位置.采用高斯—賽德?tīng)柕?循環(huán)這個(gè)過(guò)程一直到投影次數(shù)達(dá)到最大值,得到在當(dāng)前時(shí)間步長(zhǎng)內(nèi)質(zhì)點(diǎn)的最終確定位置,并根據(jù)位置隨時(shí)間的變化規(guī)律,以差分的形式來(lái)估計(jì)下一時(shí)刻質(zhì)點(diǎn)速度,預(yù)測(cè)下一時(shí)刻質(zhì)點(diǎn)位置.循環(huán)這一過(guò)程,直至仿真結(jié)束.

3 碰撞檢測(cè)

碰撞檢測(cè)是虛擬手術(shù)軟組織形變仿真的重要組成部分,碰撞檢測(cè)的目的是精確快速的判定虛擬手術(shù)中人體器官模型與虛擬手術(shù)設(shè)備之間是否發(fā)生接觸,因此碰撞檢測(cè)方法對(duì)于避免穿刺和觸發(fā)變形計(jì)算十分重要[18,19].為了保證虛擬手術(shù)仿真系統(tǒng)中生物體軟組織形變仿真的實(shí)時(shí)性,本文采用快速拾取控制點(diǎn)的方法來(lái)進(jìn)行碰撞檢測(cè),如圖3所示.肝臟表面小球是由軟組織物理模型的頂點(diǎn)數(shù)據(jù)生成的控制點(diǎn),P點(diǎn)為手術(shù)器械前端頂點(diǎn),當(dāng)手術(shù)器械即將接觸軟組織時(shí),P點(diǎn)搜索與之距離最近的控制點(diǎn)A.控制點(diǎn)A帶動(dòng)相鄰控制點(diǎn)B、C、D等產(chǎn)生移動(dòng),使得虛擬手術(shù)器械在軟組織上產(chǎn)生提拉、壓縮等控制動(dòng)作.控制點(diǎn)通過(guò)B-PBD算法更新每一時(shí)刻的位置,模擬軟組織形變.

圖3 手術(shù)器械與肝臟控制點(diǎn)的碰撞檢測(cè)Fig.3 Collision detection between surgical instrument and liver control point

4 軟組織變形實(shí)驗(yàn)

4.1 實(shí)驗(yàn)平臺(tái)

虛擬手術(shù)系統(tǒng)平臺(tái)的硬件主要包括計(jì)算機(jī)和交互(視覺(jué)交互、觸覺(jué)交互)設(shè)備兩部分.系統(tǒng)平臺(tái)的觸覺(jué)交互設(shè)備由美國(guó) Sensable Technologies Incorporated(STI)公司的開(kāi)發(fā)生產(chǎn)的Geomagic Touch的觸覺(jué)反饋設(shè)備,該設(shè)備提供了HLAPI,用于觸覺(jué)反饋力的計(jì)算,以實(shí)現(xiàn)觸覺(jué)交互操作.Geomagic Touch力反饋設(shè)備具有6 自由度,160mm×120mm×70mm的三維工作空間,0.009mm的空間分辨率,3.3N的最大施加力.圖形工作站是一類高性能計(jì)算設(shè)備,負(fù)責(zé)整個(gè)系統(tǒng)的算法計(jì)算,執(zhí)行及相關(guān)邏輯處理.為了保證算法的計(jì)算速度和模型的執(zhí)行效率,本課題選用的工作站具有較強(qiáng)的圖形圖像處理能力,具體性能參數(shù)如表2所示,圖4為虛擬手術(shù)仿真平臺(tái)硬件設(shè)備.

表2 圖形工作站性能參數(shù)Table 2 Graphic workstation performance parameters

圖4 虛擬手術(shù)仿真平臺(tái)硬件設(shè)備Fig.4 Hardware equipment of virtual surgery simulation platform

本文開(kāi)發(fā)的虛擬手術(shù)系統(tǒng)是在 Windows 10操作系統(tǒng)平臺(tái)、Visual Stdio 2017的編譯環(huán)境下,應(yīng)用程序使用C#實(shí)現(xiàn),構(gòu)建在Unity3D游戲引擎上.建立的虛擬手術(shù)仿真系統(tǒng)可在普通的 PC 機(jī)上運(yùn)行,便于后續(xù)的推廣運(yùn)用.按照實(shí)現(xiàn)的功能劃分,虛擬手術(shù)仿真系統(tǒng)在軟件層面上由三維體模型重構(gòu)模塊、視覺(jué)渲染模塊、觸覺(jué)渲染模塊、切割仿真模塊和輔助模塊等五大模塊組成.

建模數(shù)據(jù)采用Simens 64層CT掃描到的三組醫(yī)療圖像,通過(guò)設(shè)定CT值將對(duì)應(yīng)的組織提取出來(lái),并利用Mimics軟件進(jìn)行表面三維建模.

4.2 肝臟軟組織形變仿真實(shí)驗(yàn)與結(jié)果分析

本文對(duì)融合生物特性的位置動(dòng)力學(xué)(B-PBD)形變過(guò)程進(jìn)行模擬,采用肝臟組織(原始數(shù)據(jù)為肝臟CT數(shù)據(jù))作為仿真對(duì)象.幾何模型采用高精度的三角形表面網(wǎng)格,物理模型采用Tetgen劃分的德洛內(nèi)四面體網(wǎng)格,如圖5所示.面網(wǎng)格中包含了14128個(gè)三角面片及7469個(gè)頂點(diǎn);體網(wǎng)格中包含3833個(gè)四面體及1164個(gè)頂點(diǎn).

4.2.1 粘彈性模型的真實(shí)性分析

1.選取模型參數(shù)

{Kg,Kn,K,τ}是建立本構(gòu)方程所需參數(shù).根據(jù)文獻(xiàn)[20]測(cè)量的正常人體肝臟彈性系數(shù)0.3262N/mm,任意選取多組數(shù)據(jù)進(jìn)行模擬仿真對(duì)比,確定了本構(gòu)方程中所需參數(shù){0.04N/mm,0.20N/mm,0.3262N/mm,218s}.

由公式(34)得出,在確定的應(yīng)變情況下,λ由四面體邊長(zhǎng)決定.在X0=1.5mm,X0=2mm,以及X0=2.5mm的條件下,通過(guò)設(shè)置不同的四面體邊長(zhǎng),求出λ值,將仿真曲線的應(yīng)力值與實(shí)驗(yàn)數(shù)據(jù)對(duì)比,如表3～表5所示.

表3 位移X0=1.5mm時(shí)的形變力值Table 3 Deformation force at displacement X0=1.5mm

表4 位移X0=2.0mm時(shí)的形變力值Table 4 Deformation force at displacement X0=2.0mm

表5 位移X0=2.5mm時(shí)的形變力值Table 5 Deformation force at displacement X0=2.5mm

由表看出,當(dāng)邊長(zhǎng)取1.4mm,模型應(yīng)力值變化接近于肝臟應(yīng)力測(cè)量值,所以本文正面體邊長(zhǎng)取y=1.4mm劃分四面體網(wǎng)格,如圖5(c)所示.

2.粘彈性分析

通過(guò)實(shí)驗(yàn)來(lái)驗(yàn)證本文建立的本構(gòu)方程的粘彈性,如圖6(a)為力-時(shí)間曲線,圖6(b)為時(shí)間-位移曲線.

圖6 粘彈性仿真曲線Fig.6 Viscoelastic simulation curve

在圖6(a)中,當(dāng)位移x=1.5mm,2mm和2.5mm時(shí),外力何時(shí)間成反比關(guān)系,符合粘彈性材料的松弛規(guī)律.在圖6(b)中,當(dāng)F=0.35N,0.65N,1.0N時(shí),位移和時(shí)間成正比關(guān)系,符合粘彈性材料的蠕變規(guī)律.

結(jié)果表明,本文建立的軟組織三參量小變形混合填充球模型,具有良好的應(yīng)力松弛特性和應(yīng)變?nèi)渥兲匦?符合真實(shí)軟組織的變形規(guī)律.

4.2.2 軟組織形變仿真的視覺(jué)效果分析

設(shè)置彈性模量、泊松比,手術(shù)器械分別對(duì)模型施加拉力和壓力,基于B-PBD的方法,并采用拉伸約束和體積約束來(lái)模擬形變效果進(jìn)行形變模擬,肝臟、膽囊及小黃鴨形變仿真結(jié)果如圖7(a)、圖7(b)、圖7(c)所示.

圖7 B-PBD模型及質(zhì)點(diǎn)彈簧形變仿真Fig.7 Deformation simulation of the B-PBD and MSM model

采用質(zhì)點(diǎn)彈簧模型(MSM)來(lái)模擬粘彈性軟組織的計(jì)算步驟中,接觸點(diǎn)附近的彈簧會(huì)產(chǎn)生較大的超彈性效應(yīng),出現(xiàn)塌陷且無(wú)法復(fù)原的現(xiàn)象,如圖7(d)所示.通過(guò)B-PBD模型可以自然地模擬典型的軟組織變形效應(yīng).從圖中可以看出,使用融合位置動(dòng)力學(xué)的物理模型能獲得較真實(shí)的形變效果.

4.2.3 軟組織形變仿真的時(shí)間效率分析

軟組織形變時(shí)間長(zhǎng)短是衡量系統(tǒng)實(shí)時(shí)性的重要指標(biāo),為進(jìn)一步分析時(shí)間性能,采用Unity3D幀速率記錄儀記錄了軟組織形變信息.表6反應(yīng)了集中不同的軟組織模型變形的時(shí)間效率,計(jì)算變形時(shí)的迭代次數(shù)為3.

本文提出的B-PBD模型與PBD、MSM模型進(jìn)行了比較,MSM模型采用顯式歐拉法求解運(yùn)動(dòng)系統(tǒng);PBD模型及B-PBD模型包括拉伸約束和體積保持約束,表7為不同方法對(duì)上述肝臟模型的對(duì)比實(shí)驗(yàn)結(jié)果.

表7 3種方法(MSM/PBD/B-PBD)的時(shí)間效率Table 7 Time efficiency of the three methods (MSM/PBD/B-PBD)

本文采用控制點(diǎn)法進(jìn)行碰撞檢測(cè),碰撞時(shí)間幾乎可以忽略不計(jì).通過(guò)實(shí)驗(yàn)可知,基于B-PBD方法建立的軟組織模型在渲染時(shí)間和物理形變時(shí)間及幀率上都能滿足實(shí)時(shí)性的要求.

4.2.4 軟組織形變仿真的穩(wěn)定性分析

虛擬手術(shù)應(yīng)用中的軟組織變形方法能否滿足穩(wěn)定性要求至關(guān)重要.穩(wěn)定性的好壞直接影響操作者的體驗(yàn)和沉浸感.數(shù)值積分誤差使得傳統(tǒng)的質(zhì)量彈簧法容易產(chǎn)生抖動(dòng)現(xiàn)象.通過(guò)將PBD約束集成到彈簧力中,采用PBD直接控制了模型中頂點(diǎn)的位置,不會(huì)出現(xiàn)軟組織不停抖動(dòng)的問(wèn)題,同時(shí)在三參量質(zhì)點(diǎn)彈簧模型中控制了虛擬組織的剛度和阻尼,因此提高了模型的穩(wěn)定性.類似高斯迭代的約束投影的過(guò)程使該算法接近于無(wú)條件穩(wěn)定.如圖8所示.將肝臟模型提拉至圖8(a)所示位置,松開(kāi)手術(shù)器械,使其回彈至初始狀態(tài)即圖8(d),通過(guò)圖8(b)、圖8(c)可以看出,形變回復(fù)的過(guò)程中無(wú)塌陷、抖動(dòng),穩(wěn)定性較好.

圖8 肝臟模型提拉后回復(fù)效果圖Fig.8 Liver model recovery after lifting

5 結(jié) 論

為了實(shí)時(shí)模擬生物軟組織變形的典型力學(xué)特征,本文提出了一種將粘彈性質(zhì)點(diǎn)彈簧(MSM)與位置動(dòng)力學(xué)(PBD)集成的新方法.根據(jù)MSM和PBD都是通過(guò)空間離散質(zhì)點(diǎn)的運(yùn)動(dòng)位置來(lái)模擬形變特征的共通點(diǎn),并根據(jù)擴(kuò)展PBD算法的推導(dǎo)過(guò)程,將兩種模型進(jìn)行整合,即基于生物特性的位置動(dòng)力學(xué)(B-PBD).在這種方法中,基于牛頓第二定律采用顯示歐拉法進(jìn)行質(zhì)點(diǎn)下一時(shí)刻的位置預(yù)測(cè),改進(jìn)的三參量質(zhì)點(diǎn)彈簧的彈簧力和外力與PBD約束函數(shù)產(chǎn)生的約束力相結(jié)合,通過(guò)約束修正質(zhì)點(diǎn)位置.

針對(duì)改進(jìn)的三參量模型進(jìn)行數(shù)據(jù)仿真,模擬出合適的四面體邊長(zhǎng),并有針對(duì)性的構(gòu)建了力時(shí)間和位移時(shí)間方程,驗(yàn)證了生物力學(xué)的粘彈性特性.采用膽囊、小黃鴨和肝臟模型進(jìn)行力反饋實(shí)驗(yàn),并進(jìn)行了時(shí)間效率和穩(wěn)定性分析.展示了B-PBD方法在保證視覺(jué)效果的情況下,具有更高的穩(wěn)定性.同時(shí)在單線程核心計(jì)算下,碰撞檢測(cè)、形變計(jì)算及渲染的整體時(shí)間能滿足實(shí)時(shí)性的要求.實(shí)驗(yàn)表明,該方法成功地通過(guò)多種約束條件來(lái)控制軟組織的變形行為,優(yōu)于目前主流實(shí)時(shí)形變模擬方法,可以滿足虛擬手術(shù)的需要.