虛擬手術(shù)中軟組織形變建模及力反饋算法研究

陳衛(wèi)東 趙成龍，3 朱奇光* 關(guān)永貞1(燕山大學(xué)信息科學(xué)與工程學(xué)院，秦皇島 066004)

2(河北省特種光纖與光纖傳感重點(diǎn)實(shí)驗(yàn)室，秦皇島 066004)

3(海軍91821 部隊(duì)，潮州 515729)

引言

虛擬手術(shù)仿真系統(tǒng)，是虛擬現(xiàn)實(shí)(virtual reality，VR)技術(shù)在醫(yī)學(xué)領(lǐng)域的一個(gè)典型應(yīng)用。在虛擬手術(shù)訓(xùn)練中，引入力觸覺(jué)反饋，可以使醫(yī)生在訓(xùn)練時(shí)不僅能夠看到而且還能感覺(jué)到手術(shù)器官;通過(guò)與虛擬模型和環(huán)境進(jìn)行交互，形成對(duì)虛擬模型的一個(gè)完整的認(rèn)識(shí)，如同操作真實(shí)物體一樣，這無(wú)疑可使訓(xùn)練更真實(shí)、準(zhǔn)確、可靠［1］。

人體軟組織形變技術(shù)是虛擬手術(shù)仿真系統(tǒng)的關(guān)鍵技術(shù)，在虛擬環(huán)境力觸覺(jué)交互過(guò)程中，建立基于物理意義的力觸覺(jué)模型是非常重要的。大多數(shù)軟組織的物理形變模型都是基于彈性理論研究的。Majumder 等基于有限元素法(FEM)提出了具有物理意義的軟組織形變模型，雖然精確，但要求網(wǎng)格的節(jié)點(diǎn)數(shù)多、計(jì)算成本高［2］。Peterlik 等提出非線性有限元模型，提高了實(shí)時(shí)性，但仍無(wú)法滿足力觸覺(jué)反饋上千赫茲的頻率要求［3］。Wang 等用邊界元法建立模型，雖對(duì)模型的邊界進(jìn)行了離散，簡(jiǎn)化了計(jì)算，但在穩(wěn)定性方面卻存在一定的困難［4］。Zhong等提出具有物理化學(xué)特性的反應(yīng)擴(kuò)散表面模型，難于進(jìn)行復(fù)雜的手術(shù)操作［5］。Chen 在2006 年提出了基于彈簧-質(zhì)點(diǎn)的三角形表面網(wǎng)格拓?fù)浣Y(jié)構(gòu)模型［6］，模型原理簡(jiǎn)單，但穩(wěn)定性有待提高。徐少平等基于無(wú)網(wǎng)格的軟組織切割模型的研究現(xiàn)狀進(jìn)行了分析［7］，指出其未來(lái)發(fā)展方向，但穩(wěn)定性較差。吳涓等提出了一種基于沿徑向方向分割為呈同心圓分布的彈簧-質(zhì)點(diǎn)模型及實(shí)時(shí)力覺(jué)響應(yīng)算法［8］，但該模型各點(diǎn)計(jì)算不一致，不易擴(kuò)展。彈簧-質(zhì)點(diǎn)模型原理簡(jiǎn)單，易建模，能夠滿足實(shí)時(shí)性交互操作要求，但精度有待提高。因此，在保證變形模型準(zhǔn)確性的同時(shí)，提高計(jì)算的實(shí)時(shí)性，是當(dāng)前虛擬力觸覺(jué)交互系統(tǒng)亟待解決的首要問(wèn)題［9］。

傳統(tǒng)的基于表面網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)的模型，在模擬軟組織形變過(guò)程中穩(wěn)定性較差。本研究對(duì)虛擬人體軟組織表面模型在形變過(guò)程中設(shè)置了虛擬彈簧，使得表面模型形變過(guò)程中有了體的信息，以實(shí)現(xiàn)更加真實(shí)的形變效果，形變過(guò)程中算法簡(jiǎn)單，能夠很好地滿足交互系統(tǒng)的真實(shí)性、穩(wěn)定性和實(shí)時(shí)性的要求。

1 軟組織形變模型

1.1 模型原理及結(jié)構(gòu)

從生物力學(xué)的角度來(lái)看，皮膚、脂肪等軟組織可近似認(rèn)為是各向同性且均質(zhì)的［10］。但由于軟組織在受力狀態(tài)下，受力面的大小和方向都在實(shí)時(shí)變化，這使得基于運(yùn)動(dòng)學(xué)模型的解必然會(huì)產(chǎn)生震蕩。傳統(tǒng)的基于表面網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)的模型，在模擬軟組織形變過(guò)程中穩(wěn)定性較差，不能很好地表現(xiàn)軟組織的生物力學(xué)特性。本研究中的模型是基于表面正六邊形模型的改進(jìn)，為了用表面正六邊形模型近似地仿真它們所表示的軟組織結(jié)構(gòu)，得到一個(gè)比較逼真的含有體行為的仿真效果，在模型中為每個(gè)質(zhì)點(diǎn)設(shè)計(jì)了一個(gè)虛擬體彈簧［11］。虛擬體彈簧與受力方向一致，在仿真開(kāi)始之前它是不存在的，仿真開(kāi)始之后它才動(dòng)態(tài)地產(chǎn)生，主要用于力反饋的計(jì)算和增加系統(tǒng)的穩(wěn)定性。

虛擬體彈簧的變形基本思想是:施加在軟組織表面上的力，通過(guò)與用戶直接交互的那個(gè)表面接觸點(diǎn)，傳遞到周?chē)呐c之相連的其他各個(gè)質(zhì)點(diǎn)，這些質(zhì)點(diǎn)產(chǎn)生的應(yīng)力作用在其他相鄰質(zhì)點(diǎn)，這樣把力向周?chē)鷤鬟f，并同時(shí)發(fā)生變形。在系統(tǒng)中，各質(zhì)點(diǎn)表面彈簧形變量的疊加可等效為物體表面形變，各層質(zhì)點(diǎn)虛擬體彈簧彈性力的合力可等效為物體表面的接觸力。

設(shè)有虛擬體彈簧的正六邊形結(jié)構(gòu)模型結(jié)構(gòu)如圖1 所示。圖中以接觸點(diǎn)為中心，用輻射狀的表面同心圓分層。把整個(gè)形變物體分割成一系列均勻的同心環(huán)層，半徑依次為r，2r，3r，…，雖然每個(gè)同心環(huán)層中質(zhì)點(diǎn)數(shù)和半徑不同，但由于表面正六邊形結(jié)構(gòu)的對(duì)稱性，為了計(jì)算方便，可近似地認(rèn)為其半徑相同，且同一環(huán)層內(nèi)質(zhì)點(diǎn)剛性連接。相鄰層間彈簧彈性系數(shù)相同。這樣第一層有6 個(gè)質(zhì)點(diǎn)，6 個(gè)彈簧。每增加一層質(zhì)點(diǎn)數(shù)增加6 個(gè)，彈簧數(shù)增加12個(gè)。第m 層有質(zhì)點(diǎn)6m 個(gè)，彈簧6(2m -1)個(gè)。

圖1 設(shè)有虛擬體彈簧的正六面體結(jié)構(gòu)模型Fig． 1 Regular hexagon structure model with a virtual spring

1.2 模型形變的力反饋算法

當(dāng)力F 作用于質(zhì)點(diǎn)n0時(shí)，n0在力的作用下離開(kāi)初始位置發(fā)生形變，此時(shí)n0處的虛擬體彈簧發(fā)生作用。同時(shí)向下移動(dòng)的質(zhì)點(diǎn)n0通過(guò)表面彈簧帶動(dòng)相鄰層質(zhì)點(diǎn)發(fā)生位移，發(fā)生位移的各質(zhì)點(diǎn)處虛擬體彈簧也發(fā)生作用。當(dāng)達(dá)到平衡時(shí)處質(zhì)點(diǎn)n0運(yùn)動(dòng)到n'0。這里，任取受力法平面一個(gè)方向進(jìn)行分析，如圖2 所示。

圖2 形變法平面某一方向截面Fig．2 The schematic diagram in normal plane

在受力作用點(diǎn)處，根據(jù)力平衡關(guān)系有

對(duì)于第i 層質(zhì)點(diǎn)，有

式中，fTi表示切向力，fNi表示法向力，KNi表示ni質(zhì)點(diǎn)處的虛擬彈簧彈性系數(shù)。ΔZi表示ni質(zhì)點(diǎn)的位移量，即虛擬彈簧的形變長(zhǎng)度。θi為形變后ni質(zhì)點(diǎn)與ni-1質(zhì)點(diǎn)連線與水平方向夾角。

當(dāng)i = m 時(shí)，由式(1)和式(3)得

即

設(shè)相鄰兩層質(zhì)點(diǎn)的切向和法向形變，如圖3所示。

圖3 相鄰層質(zhì)點(diǎn)形變位置Fig．3 The schematic diagram of surface spring and virtual spring between adjacent layers

li為第i 層質(zhì)點(diǎn)法向虛擬彈簧的形變量，r 為表面彈簧的原始長(zhǎng)度，Δri表示切向相鄰質(zhì)點(diǎn)表面彈簧的形變量。

根據(jù)圖2 ～圖3 的關(guān)系，可知各層質(zhì)點(diǎn)的法向形變量，即發(fā)生形變的穿刺深度ΔZ，有

第i 層網(wǎng)格切向表面彈簧的彈性系數(shù)KTi分布是影響接觸力與變形計(jì)算的關(guān)鍵，這里假設(shè)切向表面彈簧的彈性系數(shù)KTi均相等。由式(2)得

令εi=為第i 層彈簧切向表面彈簧的形變率，αi= 2i - 1，則有

法向上第i 層質(zhì)點(diǎn)法向虛擬彈簧的形變量li為

在計(jì)算軟組織受力變形的模型中，計(jì)算量的大小在很大程度上取決于傳播深度m 的大小。這里m 是一個(gè)邊界參數(shù)，它表示在虛擬物體上以受力點(diǎn)為中心劃分的同心圓層數(shù)。根據(jù)文獻(xiàn)［12］的推導(dǎo)，有

式中:lm/r 代表計(jì)算區(qū)域邊界的形變量，也是邊界參數(shù);ε1則代表了第一層切向表面彈簧的形變率。

計(jì)算中選定m 后，可認(rèn)為sinθm+1≈0 。此時(shí)，式(5)變?yōu)?/p>

以上的公式推導(dǎo)是基于交互點(diǎn)n0得出，根據(jù)模型的對(duì)稱性，可以看出模型中各點(diǎn)的計(jì)算都是一致的。式(6)和式(13)中穿刺深度與傳播深度的大小決定了整個(gè)模型形變區(qū)域大小和反饋力的大小。

2 虛擬手術(shù)力反饋系統(tǒng)仿真

2.1 虛擬軟組織實(shí)時(shí)力反饋交互系統(tǒng)

虛擬軟組織實(shí)時(shí)力反饋交互系統(tǒng)的建立，首先要建立軟組織形變的幾何模型。由醫(yī)學(xué)影像數(shù)據(jù)出發(fā)，基于OpenGL 建立和渲染軟組織的幾何模型，然后在三維重構(gòu)的幾何模型中加入了彈簧-質(zhì)點(diǎn)模型，使軟組織能夠真實(shí)地反映器官的物理形變過(guò)程。在虛擬場(chǎng)景中，建立虛擬手術(shù)器械模型和軟組織模型，導(dǎo)入Matlab 的Simulink 中建立仿真模型。在碰撞檢測(cè)［13］模塊中，虛擬手術(shù)器械的交互定位用力反饋設(shè)備由觸覺(jué)筆尖控制。當(dāng)發(fā)生碰撞時(shí)，程序根據(jù)穿刺深度ΔZ 計(jì)算出反饋力的大小，操縱者能夠?qū)崟r(shí)體驗(yàn)到作用的力反饋。同時(shí)，模型形變模塊依據(jù)穿刺深度更新接觸點(diǎn)的位置，依據(jù)傳播深度m計(jì)算出形變區(qū)域，在虛擬場(chǎng)景中更新模型渲染圖，提供視覺(jué)反饋。虛擬手術(shù)仿真實(shí)現(xiàn)的具體流程如圖4 所示。

圖4 虛擬手術(shù)仿真系統(tǒng)流程Fig．4 The schematic diagram of force feedback interaction system

2.2 模型參數(shù)的選擇

軟組織實(shí)時(shí)力反饋交互系統(tǒng)的變形仿真一般要求比較高的實(shí)時(shí)性，這就對(duì)該系統(tǒng)的計(jì)算速度和刷新頻率有很高的要求。要使視覺(jué)、力覺(jué)再現(xiàn)具有連續(xù)性，視覺(jué)顯示的軟組織形變刷新率應(yīng)該不小于30 Hz，力覺(jué)再現(xiàn)的刷新頻率應(yīng)該在1 kHz 以上［14］。即要求力反饋的計(jì)算具有快速性的特點(diǎn)。

觸覺(jué)交互設(shè)備是用 SensAble 科技公司的PHANTOM 設(shè)備，它能夠準(zhǔn)確跟蹤人手的三維運(yùn)動(dòng)，并將由虛擬模型計(jì)算的虛擬力實(shí)時(shí)地反饋給操作者，在力覺(jué)上提供逼真的沉浸感。PHANTOM 設(shè)備的最大輸出反饋力為3.3 N。

傳播深度m 是控制形變區(qū)域的一個(gè)參數(shù)，它的大小不僅影響著計(jì)算量的大小，而且影響著計(jì)算精度。當(dāng)形變穿刺深度一定時(shí)，傳播深度m 越大，精度越高，反饋力越接近真實(shí)力?？紤]到計(jì)算量的問(wèn)題，實(shí)驗(yàn)中可根據(jù)具體情況確定傳播深度m 的大小。根據(jù)式(14)，取虛擬體彈簧彈性系數(shù)KN0= 50 N/m，r = 0.002 m 進(jìn)行模擬仿真，得到反饋力與傳播深度m 的關(guān)系曲線，如圖5 所示。

圖5 反饋力F 與傳播深度m 的關(guān)系Fig．5 The relation between the force feedback and the propagation depth m

由式(13)可知，m 的大小由lm/r 和ε1共同決定，選取m =25。表面彈簧原始長(zhǎng)度r 為固定長(zhǎng)度。根據(jù)式(6)和式(12)，受力時(shí)軟組織虛擬體彈簧形變量與傳播深度m 的關(guān)系如圖6 所示。

圖6 虛擬體彈簧形變量與傳播深度m 的關(guān)系Fig． 6 The deformation of virtual spring with propagation depth m

當(dāng)接觸點(diǎn)形變率ε1取不同值時(shí)，得到不同的形變曲線。圖6 中表明:在相同力的作用下，不同組織的形變率ε1不同，軟組織的形變也不同。ε1越大，表明組織越柔軟，受力變形越大。當(dāng)然，軟組織的柔軟度還與虛擬體彈簧彈性系數(shù)KN的大小［15］相關(guān)，由二者共同決定。

虛擬體彈簧彈性系數(shù)KN的選取不僅要考慮模型形變的大小，還要考慮力反饋設(shè)備反饋力上限的制約因素。根據(jù)式(14)，令輸出力F 小于設(shè)備最大輸出值。

2.3 系統(tǒng)仿真結(jié)果

使用SensAble 科技公司的PHANTOM 觸覺(jué)交互設(shè)備，基于AMD Athlon(tm)Ⅱ2.80 GHz 雙核處理器，1.75GB 內(nèi)存，NVIDIA GeForce7025 顯卡的臺(tái)式電腦，在 WindowS XP 操作系統(tǒng)下，采用Visua1C++6.0 開(kāi)發(fā)環(huán)境并結(jié)合OpenGL 三維圖形標(biāo)準(zhǔn)，搭建了虛擬軟組織實(shí)時(shí)力反饋交互系統(tǒng)。

基于該實(shí)驗(yàn)系統(tǒng)，實(shí)現(xiàn)了虛擬軟組織網(wǎng)格(皮膚、脂肪等)在虛擬手術(shù)器械(虛擬注射器、虛擬手術(shù)刀等)單點(diǎn)交互下的按壓變形仿真。當(dāng)操作者控制觸覺(jué)交互設(shè)備運(yùn)動(dòng)時(shí)，虛擬手術(shù)器械也隨之運(yùn)動(dòng)。當(dāng)碰撞檢測(cè)算法檢測(cè)到虛擬手術(shù)器械虛擬軟組織發(fā)生碰撞時(shí)，在給定接觸點(diǎn)處運(yùn)用本研究建立的形變模型及力反饋算法，提供給操作者以視覺(jué)與力覺(jué)信息，產(chǎn)生真實(shí)的交互。實(shí)驗(yàn)結(jié)果如圖7 所示，圖中顯示的是軟組織受外力在單點(diǎn)按壓作用下的形變效果。

圖7 實(shí)驗(yàn)結(jié)果。(a)交互圖;(b)形變前;(c)形變后Fig．7 The simulation results． (a)Interaction diagram;(b)Before deformation;(c)After deformation

在實(shí)時(shí)力反饋模塊中，提取出實(shí)時(shí)交互過(guò)程中按壓操作的反饋力，如圖8 所示。

圖8 交互過(guò)程中按壓操作的實(shí)時(shí)反饋力Fig．8 The real-time feedback force under the pressure

實(shí)驗(yàn)結(jié)果顯示，在施加外力作用之后，軟組織形變模型就像一個(gè)實(shí)體一樣很有體積感，在形變過(guò)程中不會(huì)發(fā)生形變失真現(xiàn)象，撤銷(xiāo)外力之后表面變形幾乎完全恢復(fù)。反饋力的計(jì)算簡(jiǎn)單，一次刷新小于1 ms 也滿足刷新頻率的要求。在交互過(guò)程中，力覺(jué)感受平穩(wěn)，模擬效果逼真，富有彈性。

3 結(jié)論

針對(duì)實(shí)時(shí)軟組織力觸覺(jué)交互系統(tǒng)的需要，提出了一種具有體信息的表面正六邊形結(jié)構(gòu)模型。與文獻(xiàn)［12］相比，該模型中每層的各個(gè)質(zhì)點(diǎn)受力與變形情況近似相同，計(jì)算量小，保證了力觸覺(jué)反饋計(jì)算的高刷新頻率的要求，模擬效果逼真，各個(gè)質(zhì)點(diǎn)的相對(duì)位移和彈性力的疊加等效為軟組織的變形和接觸力。在PHANTOM 觸覺(jué)交互設(shè)備上，通過(guò)虛擬手術(shù)工具，對(duì)虛擬軟組織進(jìn)行單點(diǎn)的按壓變形仿真。結(jié)果表明，該模型能夠很好地滿足交互系統(tǒng)的真實(shí)性、穩(wěn)定性和實(shí)時(shí)性的要求。

［1］ 謝叻，張艷，張?zhí)煊?，? 虛擬手術(shù)中的力學(xué)變形和力覺(jué)感知［J］. 醫(yī)用生物力學(xué)，2006，21(3):241 -245.

［2］ Santanu M，Amit R，Subrata P. Simulation of hip fracture in sideways fall using a 3D finite element model of pelvis-femur-soft tissue complex with simplified representation of whole body［J］.Medical Engineering ＆ Physics，2007，29:1167 -1178.

［3］ Igor P，Mert S，Cagatay B，et al. 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，et al. Virtual reality simulation of surgery with haptic feedback based on the boundary element method［J］. Computers and Structures，2007，85:331 -339.

［5］ Zhong Yongmin，Bijan S，Julian S，et al. Soft tissue deformation with reaction-diffusion process for surgery simulation［J］. Journal of Visual Languages ＆ Computing，2012，23:1 -12.

［6］ Chen P，Barner KE，Steiner KV. A displacement driven realtime deformable model for haptic surgery simulation ［C］ //Proceedings of 14th Symposium on Haptics Interfaces for Virtual Environment and Teleoperator Systems. Alexandria: IEEE，2006:499 -505.

［7］ 徐少平，李春泉，江順亮，等. 基于無(wú)網(wǎng)格的軟組織切割模型的研究進(jìn)展［J］. 中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，2012，31(1):123-128.

［8］ 吳涓，宋愛(ài)國(guó)，李建清. 用于虛擬力覺(jué)再現(xiàn)的彈簧- 質(zhì)點(diǎn)力/變形模型研究［J］. 系統(tǒng)仿真學(xué)報(bào)，2006，18(11):3152 -3156.

［9］ 張小瑞，宋愛(ài)國(guó)，孫偉，等. 用于軟組織變形仿真的層狀菱形鏈連接模型［J］. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2010，22(3):473 -479.

［10］ 王彥臻，熊岳山，徐凱，等. 改進(jìn)的采用表面網(wǎng)格的彈簧振子模型［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2007，19(2):168-171

［11］ 鮑春波，王博亮，劉卓，等. 一種用于軟組織變形仿真的動(dòng)態(tài)質(zhì)點(diǎn)彈簧模型［J］. 系統(tǒng)仿真學(xué)報(bào)，2006，18(4):847 -851.

［12］ 吳涓，宋愛(ài)國(guó)，李建清. 一種基于物理意義的快速力反饋形變模型及實(shí)時(shí)力覺(jué)響應(yīng)算法［J］. 傳感技術(shù)學(xué)報(bào)，2005，18(1):90 -94.

［13］ Asma AE，Elsayed EH，Magda BF. Rapid collision detection for deformable objects using inclusion-fields applied to cloth simulation［J］. Journal of Advanced Research，2011，16:25.

［14］ Chang Yeonghwa， Chen Yungte， Chang Chiawen， et al.Development scheme of haptic-based system for interactive deformable simulation［J］. Computer-Aided Design，2010，42:414 -424.

［15］ Suriya N，Cenk C. Determination of elasticity parameters in lumped element (mass-spring)models of deformable objects［J］.Graphical Models，2010，72:61 -73.