吳 昊，張峰峰，2*，詹 蔚，王榮淼

(1.蘇州大學(xué)機(jī)電工程學(xué)院，江蘇 蘇州 215006；2.蘇州大學(xué)蘇州納米科技協(xié)同創(chuàng)新中心，江蘇 蘇州 215123；3.蘇州大學(xué)附屬第一醫(yī)院，江蘇 蘇州 215006)

1 引言

虛擬手術(shù)是虛擬現(xiàn)實(shí)技術(shù)在現(xiàn)代醫(yī)學(xué)領(lǐng)域的一個(gè)典型應(yīng)用，是一種可以替代傳統(tǒng)醫(yī)學(xué)培訓(xùn)方法并且具有巨大應(yīng)用價(jià)值的新興學(xué)科[1]。目前，虛擬手術(shù)的研究主要集中在軟組織切割變形及縫合模擬，其中縫合線的打結(jié)作為手術(shù)縫合模擬中的關(guān)鍵環(huán)節(jié)也成為了一項(xiàng)重要研究?jī)?nèi)容[2]。手術(shù)縫合線作為一種柔性體，具有很強(qiáng)的抗拉伸變形能力并且能被隨意彎曲、折疊。由于縫合線獨(dú)特的形變特性，對(duì)于縫合線建模、動(dòng)力學(xué)方程計(jì)算以及自碰撞檢測(cè)算法等提出了更高的要求。

目前，針對(duì)虛擬手術(shù)縫合線等其它柔性體的仿真研究國(guó)內(nèi)外已有眾多研究成果，大致可以分成基于幾何特征的仿真研究和基于物理特征的仿真研究[2]。Lenoir[4]等人使用基于樣條曲線的幾何仿真來(lái)模擬縫合線。根據(jù)能量方程，將縫線的變形最終轉(zhuǎn)化為樣條之間控制點(diǎn)的運(yùn)動(dòng)方程，并采用隱式歐拉積分法進(jìn)行計(jì)算，但該方法不能很好的描述物理特性和反饋力。Joel Brown[5]等人提出了另一種幾何仿真方法，用一系列的圓柱體通過(guò)首尾相連的方式表示線模型，利用FTL(Follow The Leader)算法，通過(guò)施加相應(yīng)的物理約束保證線的特性。該方法的優(yōu)點(diǎn)為簡(jiǎn)化大量計(jì)算，但是運(yùn)動(dòng)簡(jiǎn)單，穩(wěn)定性較差。賈世宇[6]等人基于Joel Brown中的剛性桿模型和FTL(Follow The Leader)形變算法，并通過(guò)剛性桿結(jié)構(gòu)和AABB層次包圍結(jié)構(gòu)實(shí)現(xiàn)繩子的自碰撞響應(yīng)以及繩子與其它剛體之間的碰撞，并設(shè)置臨時(shí)接觸約束實(shí)現(xiàn)等效摩擦力，完成實(shí)時(shí)繩子打結(jié)仿真并帶有力反饋；但是未考慮縫合線的物理特性，缺少動(dòng)力學(xué)方程。Lang[7]提出了基于位置動(dòng)力學(xué)(PBD)的Cosserat桿理論的手術(shù)線的實(shí)時(shí)仿真方法。該方法實(shí)現(xiàn)了手術(shù)線的穩(wěn)定纏繞和打結(jié)，但是該算法基于物理模型，需要求解多個(gè)微分方程，計(jì)算量較大。梁民倉(cāng)[9]等人對(duì)柔性繩索進(jìn)行了研究，基于質(zhì)點(diǎn)彈簧模型的結(jié)構(gòu)彎曲彈簧模型，通過(guò)添加彎曲彈簧以限制繩子質(zhì)點(diǎn)處的不正常彎曲現(xiàn)象，保證柔性繩索的穩(wěn)定性和精度。王崴[10]等人提出了一種蜂窩狀彈簧-質(zhì)點(diǎn)模型應(yīng)用于繩索形變仿真，質(zhì)點(diǎn)間設(shè)置4種彈簧模擬繩索內(nèi)部彎曲、扭轉(zhuǎn)、拉伸等形變特性，采用非共線過(guò)濾器的剔除算法實(shí)現(xiàn)自碰撞檢測(cè)，滿足繩索模擬時(shí)的實(shí)時(shí)性、快速性，但模型結(jié)構(gòu)過(guò)于復(fù)雜。

與上述文獻(xiàn)[4]～[10]的研究相比，本文針對(duì)縫合線的物理特性，提出了MSS的物理建模方法并結(jié)合物體非線性及彎曲力學(xué)特性對(duì)其進(jìn)行動(dòng)力學(xué)解算。其次通過(guò)跟蹤控制點(diǎn)實(shí)現(xiàn)了觸覺設(shè)備對(duì)縫合線的運(yùn)動(dòng)控制，有效解決了縫合線運(yùn)動(dòng)仿真時(shí)效果失真等缺點(diǎn)。在傳統(tǒng)的層次包圍盒碰撞檢測(cè)方法的基礎(chǔ)上，提出了一種模擬受力方法實(shí)現(xiàn)縫合線打結(jié)過(guò)程中的自碰撞響應(yīng)，增強(qiáng)了縫合線打結(jié)過(guò)程中的實(shí)時(shí)性和逼真度，最后完成了縫合線實(shí)時(shí)打結(jié)仿真。

2 縫合線模型描述

2.1 模型建立

根據(jù)手術(shù)縫合線特性，采用質(zhì)點(diǎn)彈簧系統(tǒng)(MSS)建立縫合線物理模型，如圖1所示。MSS方法是將縫合線離散成一系列有質(zhì)量的粒子(質(zhì)點(diǎn))，其位置信息為Pi(i=0，1，2…n)。質(zhì)點(diǎn)間通過(guò)n個(gè)無(wú)質(zhì)量的彈簧連接，設(shè)置質(zhì)點(diǎn)的質(zhì)量為m以及彈簧的初始長(zhǎng)度為l0。縫合線的物理特性由其內(nèi)部的質(zhì)點(diǎn)彈簧模型決定，只需調(diào)整質(zhì)點(diǎn)彈簧的物理參數(shù)就能得到不同性能的縫合線?？p合線的變形由質(zhì)點(diǎn)運(yùn)動(dòng)產(chǎn)生，并受到連接彈簧的控制。

圖1 縫合線模型

為了增加縫合線的光滑度，避免鋸齒現(xiàn)象，將長(zhǎng)為L(zhǎng)、半徑為r的一個(gè)完整的圓柱體作為縫合線的幾何模型。將物理模型放置在圓柱體的幾何中心線上，并構(gòu)建模型表面網(wǎng)格頂點(diǎn)與內(nèi)部質(zhì)點(diǎn)間的映射關(guān)系，使得內(nèi)部質(zhì)點(diǎn)彈簧模型的形變效果轉(zhuǎn)化為外表面幾何模型的變形效果，實(shí)現(xiàn)縫合線的變形與打結(jié)的圖像渲染效果[11]。

2.2 模型受力分析

手術(shù)縫合線為一種典型的形變線性體，其軸向長(zhǎng)度比其它兩個(gè)維度大得多。因此，縫合線在變形時(shí)表現(xiàn)出明顯的非線性特性，其彎曲應(yīng)變較大，且具有很強(qiáng)的抗拉伸能力。針對(duì)縫合線物理特性，需要對(duì)MSS模型進(jìn)行合理的受力分析[12]。

在理想的彈簧模型中，根據(jù)霍克定律，彈簧拉伸力與彈簧的伸長(zhǎng)量成正比。然而，現(xiàn)實(shí)中的縫合線具有很強(qiáng)的抗拉伸的特性，當(dāng)變形加大時(shí)，其彈力會(huì)急劇增加，防止其形變過(guò)大。計(jì)算彈簧拉伸力需要綜合考慮縫合線作為非線性物體的特性?？p合線的非線性拉伸力計(jì)算公式為

(1)

為防止彈簧發(fā)生過(guò)度形變，維護(hù)系統(tǒng)穩(wěn)定，質(zhì)點(diǎn)運(yùn)動(dòng)時(shí)應(yīng)受到彈簧的阻礙。根據(jù)粘滯阻尼理論，彈簧阻尼力與相連質(zhì)點(diǎn)速度有關(guān)，其公式如下

(2)

其中kd為彈簧阻尼系數(shù)，vi和vi+1分別為質(zhì)點(diǎn)Pi和Pi+1的速度，可由式(8)求得。

縫合線模型中由于未能限制質(zhì)點(diǎn)之間彈簧的彎曲度，造成線模型因受力發(fā)生彎曲或折疊形變過(guò)度，不能很好的模擬縫合線的幾何效果。因此，模擬縫合線的彎曲特性是保證線模型平滑運(yùn)動(dòng)的基礎(chǔ)?？p合線的彎曲特性是通過(guò)添加質(zhì)點(diǎn)上的彎曲力來(lái)實(shí)現(xiàn)的。彎曲力是指柔性體由彎曲狀態(tài)向伸直狀態(tài)恢復(fù)時(shí)所產(chǎn)生的力。質(zhì)點(diǎn)所受的彎曲力是由它前后兩個(gè)質(zhì)點(diǎn)間的彈簧彎曲程度決定，基本思想是：將兩個(gè)連接的彈簧設(shè)置為三角形(兩個(gè)彈簧不在同一直線上)，假設(shè)彎曲力是將端點(diǎn)推到完全展開的位置。兩個(gè)連接段的長(zhǎng)度保持不變，只有垂直于彈簧線段的力分量才用于端點(diǎn)。如圖2為彎曲彈簧模型。

圖2 彎曲彈簧模型

彎曲彈簧計(jì)算公式如下

(3)

(4)

(5)

3 縫合線的運(yùn)動(dòng)模擬

3.1 運(yùn)動(dòng)模型求解

縫合線的質(zhì)量質(zhì)點(diǎn)受到內(nèi)力和外力的共同作用，根據(jù)牛頓第二定律，其動(dòng)力學(xué)方程可以表示為

(6)

顯式Euler法計(jì)算簡(jiǎn)單，但精度較低，系統(tǒng)運(yùn)行不穩(wěn)定。隱式Euler法在大步長(zhǎng)仿真時(shí)系統(tǒng)穩(wěn)定性較高，但計(jì)算量較大、速度慢，降低了系統(tǒng)仿真的實(shí)時(shí)性。Verlet法為居中計(jì)算，精度比基于前向計(jì)算的顯式歐拉方法高，穩(wěn)定性高，計(jì)算復(fù)雜度也較低。本文采用Verlet法對(duì)質(zhì)點(diǎn)彈簧模型的動(dòng)力學(xué)方程進(jìn)行求解。針對(duì)式(6)使用verlet法求解質(zhì)點(diǎn)位置和速度方程如下

(7)

(8)

其中Pi(t)、Fi(t)、vi(t)分別是在t時(shí)刻下質(zhì)點(diǎn)的位置、所受的總力以及速度，Pi(t-Δt)為上一時(shí)刻質(zhì)點(diǎn)的位置，Δt為時(shí)間步長(zhǎng)。由上式可知，質(zhì)點(diǎn)在t時(shí)刻下的位置Pi(t)和所受的總力Fi(t)及t-Δt時(shí)刻下的位置Pi(t-Δt)，并根據(jù)間隔步長(zhǎng)為Δt，可得到在t時(shí)刻的質(zhì)點(diǎn)速度以及t+Δt時(shí)刻下質(zhì)點(diǎn)位置。只要設(shè)置的時(shí)間步長(zhǎng)Δt足夠短，就能保持積分算法穩(wěn)定，滿足實(shí)時(shí)交互的需求。

3.2 基于控制點(diǎn)的運(yùn)動(dòng)軌跡追蹤

基于MSS方法建立縫合線模型，由于縫合線形變較大、難以控制，會(huì)降低了線模型運(yùn)動(dòng)模擬的真實(shí)性，所以采用跟蹤控制點(diǎn)的方法對(duì)線模型運(yùn)動(dòng)軌跡追蹤。根據(jù)虛擬耦合技術(shù)[13]將觸覺反饋設(shè)備的位置信息轉(zhuǎn)化為縫合線上的任意質(zhì)點(diǎn)的位置信息，該質(zhì)點(diǎn)即為控制點(diǎn)。在每一個(gè)時(shí)間步長(zhǎng)下，觸覺工具對(duì)控制點(diǎn)施加外力，控制點(diǎn)進(jìn)行位置更新。與此同時(shí)，求得與控制點(diǎn)相連的質(zhì)點(diǎn)上的總力，并通過(guò)動(dòng)力學(xué)方程計(jì)算該質(zhì)點(diǎn)的位移與速度，即與控制點(diǎn)相連的質(zhì)點(diǎn)發(fā)生位移變化，再由這些質(zhì)點(diǎn)繼續(xù)帶動(dòng)周圍的質(zhì)點(diǎn)移動(dòng)。由此可知，任何一個(gè)質(zhì)點(diǎn)都是依據(jù)與它相鄰的質(zhì)點(diǎn)的狀態(tài)來(lái)計(jì)算的，縫合線模型為一個(gè)動(dòng)態(tài)的形變模型。

圖3顯示了跟蹤控制點(diǎn)的縫合線運(yùn)動(dòng)模擬。質(zhì)點(diǎn)4為控制點(diǎn)，且在外力的作用下發(fā)生位移變化，進(jìn)而引起與之相連的質(zhì)點(diǎn)3和質(zhì)點(diǎn)5發(fā)生位移變化，然后質(zhì)點(diǎn)3和質(zhì)點(diǎn)5又帶動(dòng)與之相連的其它質(zhì)點(diǎn)運(yùn)動(dòng)，從而使線模型產(chǎn)生跟隨控制點(diǎn)的平滑運(yùn)動(dòng)效果。

圖3 跟蹤控制點(diǎn)

4 縫合線的自碰撞檢測(cè)與響應(yīng)

縫合線的自碰撞是實(shí)現(xiàn)打結(jié)過(guò)程真實(shí)性和實(shí)時(shí)性的主要因素，也是打結(jié)仿真能否成功的關(guān)鍵。由于縫合線本身是一個(gè)細(xì)長(zhǎng)的柔性體，沒有太大的深度，并且隨時(shí)間線模型形狀不斷發(fā)生變化。因此，快速精確地碰撞檢測(cè)對(duì)縫合線打結(jié)模擬的逼真程度有著重要影響[14]。

根據(jù)縫合線分段式物理模型，使用傳統(tǒng)的基于層次包圍體結(jié)構(gòu)(Bounding Volume Hierarchy，BVH)的動(dòng)態(tài)碰撞檢測(cè)算法，包圍體為球體。根據(jù)縫合線物理模型形狀自底而上構(gòu)建層次包圍球體樹。每一個(gè)線段被一個(gè)最底層葉子包圍球包圍。相鄰兩個(gè)包圍球形成一個(gè)父節(jié)點(diǎn)包圍球，以此類推根節(jié)點(diǎn)為整個(gè)縫合線的包圍球?？p合線的自碰撞檢測(cè)算法是將層次包圍球的自身與自身比較，根據(jù)自上而下搜索原則，搜索所有兩個(gè)父節(jié)點(diǎn)包圍球，如果兩個(gè)父節(jié)點(diǎn)球體沒有相交則一定不包含碰撞片段，不需要搜索其葉子球體。如果父節(jié)點(diǎn)球體相交則繼續(xù)搜索其葉子球體。若葉子球體相交，需檢測(cè)這兩個(gè)球體中心的距離是否小于其半徑之和，若球心距離小于半徑之和表示發(fā)生碰撞。在縫合線模型上任意兩個(gè)相鄰的葉子球體都不發(fā)生碰撞。由于縫合線是可形變物體，層次包圍球需要根據(jù)模型形狀在每個(gè)時(shí)間步長(zhǎng)下更新。

如圖4所示PiPi+1和PjPj+1為縫合線任意兩個(gè)不相鄰的線段彈簧，縫合線碰撞檢測(cè)是通過(guò)計(jì)算兩線段包圍球是否相交，即球心距離Oab(即為碰撞距離)是否小于包圍球半徑之和R=Ra+Rb。如果球心距離小于包圍球半徑，則認(rèn)為線段彈簧發(fā)生碰撞。對(duì)于發(fā)生碰撞的兩線段彈簧，碰撞響應(yīng)處理算法應(yīng)是將二者推開一定距離使其不再碰撞。本文采用一種模擬受力的方法，對(duì)彈簧兩端質(zhì)點(diǎn)添加排斥力避免縫合線模型發(fā)生穿透現(xiàn)象，該排斥力通過(guò)穿透深度法計(jì)算得出，其計(jì)算公式為

(9)

圖4 線模型碰撞

5 縫合線打結(jié)仿真

仿真程序運(yùn)行在惠普Z(yǔ)820工作站上。其配置為雙Intel(R)Xeon(R)3.3GHz CPU。操作系統(tǒng)為Windows 7 Professional 。觸覺反饋設(shè)備是瑞士Force Dimension公司制造的Omega 7.0。仿真程序是在 Microsoft Visual studio 2010平臺(tái)下，使用標(biāo)準(zhǔn)C++編寫，結(jié)合CHAI 3D中封裝的OpenGL搭建縫合線仿真系統(tǒng)。

圖5 系統(tǒng)仿真平臺(tái)

基于物理模型的縫合線仿真中參數(shù)的選擇尤為重要，模型的參數(shù)決定了整個(gè)仿真系統(tǒng)的穩(wěn)定性以及效果的真實(shí)性。經(jīng)過(guò)多次試驗(yàn)，確定了表1中縫合線的參數(shù)。

表1 縫合線參數(shù)設(shè)定

縫合線運(yùn)動(dòng)仿真如圖6所示，圖中的黑色小球代表觸覺反饋設(shè)備的作用點(diǎn)，也是縫合線的運(yùn)動(dòng)控制點(diǎn)，用戶通過(guò)觸覺反饋設(shè)備拖動(dòng)該控制點(diǎn)引導(dǎo)縫合線在虛擬空間運(yùn)動(dòng)和變形。圖中可見通過(guò)跟蹤控制點(diǎn)方法可以實(shí)現(xiàn)縫合線的運(yùn)動(dòng)控制，且縫合線的運(yùn)動(dòng)線條及大小拐彎處的表現(xiàn)非常自然流暢。

圖6 縫合線運(yùn)動(dòng)控制

縫合線打結(jié)仿真如圖7所示。由于只有一臺(tái)觸覺反饋設(shè)備，需將縫合線模型的一端固定，并利用觸覺設(shè)備移動(dòng)縫合線的另一端完成縫合線打結(jié)仿真。圖7a為縫合線線運(yùn)動(dòng)初始狀態(tài)。在每一個(gè)時(shí)間步長(zhǎng)下，層次包圍球算法將不斷地對(duì)縫合線模型進(jìn)行自碰撞檢測(cè)。若檢測(cè)到線段彈簧發(fā)生自碰撞，則在相關(guān)彈簧線段兩端質(zhì)點(diǎn)處強(qiáng)加排斥力，對(duì)線段間的運(yùn)動(dòng)產(chǎn)生“阻撓”，阻止縫合線穿透現(xiàn)象的發(fā)生，如圖7b。圖7c為縫合線的結(jié)處于松散狀態(tài)。當(dāng)用戶通過(guò)觸覺設(shè)備拉動(dòng)縫合線的一端時(shí)，縫合線逐漸收緊，結(jié)以自然的方式結(jié)合得更緊，并隨著距離的減小排斥力呈指數(shù)級(jí)增長(zhǎng)，最終結(jié)將處于平衡狀態(tài)，如圖7d所示形成了一個(gè)緊密的結(jié)。

圖7 縫合線打結(jié)過(guò)程

圖8 排斥力

圖8顯示了縫合線模型在2000個(gè)時(shí)間步長(zhǎng)內(nèi)(包括一個(gè)打結(jié)操作和一系列線運(yùn)動(dòng)的過(guò)程)排斥力的變化過(guò)程。從圖中可以看出，一些明顯的峰值出現(xiàn)在靠近打結(jié)處的質(zhì)點(diǎn)上，并隨著打結(jié)完成排斥力將逐漸趨于0。在線運(yùn)動(dòng)的過(guò)程中，排斥力繼續(xù)作用。

經(jīng)上述實(shí)驗(yàn)可知，模擬縫合線模型內(nèi)部非線性拉伸力、阻尼力以及添加彈簧彎曲特性，能有效地呈現(xiàn)出縫合線的形變特性，實(shí)現(xiàn)縫合線模型的平滑穩(wěn)定運(yùn)動(dòng)。通過(guò)模擬受力的方法避免了縫合線自碰撞部分穿透現(xiàn)象的發(fā)生，這種模擬受力的方法不僅能夠很好地實(shí)現(xiàn)縫合線的打結(jié)模擬，還能夠?qū)崿F(xiàn)力的交互，大大增加了虛擬仿真的沉浸感和真實(shí)性。在整個(gè)縫合線打結(jié)仿真過(guò)程中，圖像刷新頻率均滿足30Hz以上的要求，該仿真過(guò)程實(shí)時(shí)性較好，圖像反饋真實(shí)、流暢。

6 結(jié)束語(yǔ)

本文針對(duì)虛擬縫合手術(shù)中縫合線打結(jié)進(jìn)行仿真。根據(jù)縫合線形變特性，在傳統(tǒng)的MSS物理模型中引入了非線性拉伸、阻尼以及彎曲特性，避免了線模型形變過(guò)大且拐點(diǎn)不真實(shí)的現(xiàn)象。為了保證模型穩(wěn)定性及精度要求，采用Verlet法求解動(dòng)力學(xué)方程提高計(jì)算速度。通過(guò)跟蹤控制點(diǎn)方法實(shí)現(xiàn)力反饋設(shè)備對(duì)縫合線的運(yùn)動(dòng)控制?；趥鹘y(tǒng)的層次包圍球碰撞檢測(cè)算法，提出了一種模擬受力方法解決打結(jié)過(guò)程中縫合線模型自穿透問題。最后完成了縫合線的打結(jié)仿真并帶有良好的觸覺交互。通過(guò)仿真及實(shí)驗(yàn)，驗(yàn)證了本文算法的有效性，可以逼真地模擬縫合打結(jié)操作過(guò)程，且實(shí)時(shí)性較好，穩(wěn)定性高。下一步考慮結(jié)合真實(shí)的手術(shù)縫合線打結(jié)，將縫合線與手術(shù)縫合針連接實(shí)現(xiàn)剛?cè)峤佑|的動(dòng)力學(xué)打結(jié)仿真。