一種改進的虛擬手術中人體軟組織形變方法

何巍，王魏平，師為禮，苗語，何飛，楊華民

（長春理工大學　計算機科學與技術學院，長春　130022）

一種改進的虛擬手術中人體軟組織形變方法

何巍，王魏平，師為禮，苗語，何飛，楊華民

（長春理工大學計算機科學與技術學院，長春130022）

虛擬手術是將計算機虛擬現(xiàn)實技術應用于現(xiàn)代醫(yī)學，利用醫(yī)學圖像數(shù)據(jù)，重構(gòu)出虛擬人體軟組織模型，利用觸覺交互設備，在視覺和觸覺上為使用者提供手術場景的模擬和仿真。其中，軟組織形變的仿真和建模是虛擬手術系統(tǒng)核心的研究內(nèi)容之一。針對虛擬手術中軟組織形變仿真，提出了一種改進的彈簧—質(zhì)點模型。在正四邊形網(wǎng)格拓撲結(jié)構(gòu)表面模型中，增加了結(jié)構(gòu)彈簧、剪切彈簧和彎曲彈簧，使表面模型具有體模型的特征。該模型繼承了傳統(tǒng)彈簧—質(zhì)點模型的建?？?、原理簡單、仿真速度快等優(yōu)點，同時還具有控制形變區(qū)域的能力，從而提高了仿真模型的精確度，滿足交互系統(tǒng)的真實性和穩(wěn)定性需求。實驗表明該方法有效增強了形變仿真的體積特征，適合于軟組織形變較大時的實時仿真，滿足了虛擬手術系統(tǒng)實時性的要求。

虛擬現(xiàn)實；彈簧—質(zhì)點模型；軟組織形變；碰撞檢測

隨著計算機科學技術的不斷發(fā)展，虛擬現(xiàn)實技術（Virtual Reality Technology VRT）已被廣泛應用于現(xiàn)代科學研究的各個領域［1，2］，其中利用虛擬現(xiàn)實技術進行生物醫(yī)學仿真是當前的一個研究熱點，虛擬手術系統(tǒng)是其中的典型應用。

虛擬手術系統(tǒng)是對醫(yī)學數(shù)據(jù)進行可視化并實現(xiàn)實時交互，模擬手術過程中組織器官變形、感官反饋等可能遇到的各種現(xiàn)象的虛擬現(xiàn)實應用系統(tǒng)。其中在軟組織器官形變方面，生物力學領域的研究提出了很多復雜的數(shù)學模型，能夠比較準確地描述軟組織的形變，軟組織大多數(shù)物理模型都是建立在彈性力學理論上的。Santanu等人的基于有限元方法（Finite Element Model，F(xiàn)EM）提出了具有生物特征的物理形變模型，雖然精確度很高，但是網(wǎng)格數(shù)量很多，導致模型實時性差［3］；Zhong等人提出具有物理化學特征的表面模型，但很難進行復雜的虛擬手術操作［4］；陳衛(wèi)東等改進了基于表面正六邊形模型，雖然計算量小，實時性高，但是仿真精度有待提高［5］；王瑞藝提出的基于無網(wǎng)格SPH的物理模型，雖然能夠動態(tài)地反映軟組織的實時形變，但交互方面需要進一步改善［6］；李艷東，朱玲等人結(jié)合ALMSDM的特點提出了頂點法向量局部更新與預計算策略，提高了實時性，但在大范圍變形的情況下，穩(wěn)定性不足［7］；Wang等用邊界元法建立模型，雖對模型的邊界進行了離散，簡化了計算，但穩(wěn)定性較差［8］。因此，在保證模型形變精確度的同時，提高交互的實時性，是我們虛擬手術操作中需要解決的首要問題。傳統(tǒng)的基于表面網(wǎng)格拓撲結(jié)構(gòu)的物理模型，存在兩方面問題，沒有體模型的特征或者形變真實性不足。

本論文在保證實時性的同時，在虛擬的人體軟組織表面模型中設置了結(jié)構(gòu)彈簧、剪切彈簧和彎曲彈簧，使模型具備了體模型的特征，從而提高了仿真精度，達到了更好的仿真效果。同時保證形變實時性高，很好地滿足虛擬手術中交互過程的真實性、實時性和穩(wěn)定性需求。

1　軟組織形變物理模型

彈簧-質(zhì)點模型（Mass Spring Method，MSM）是一種經(jīng)典的物理建模方法，其主要思想是把仿真的軟組織用質(zhì)點離散，質(zhì)點與質(zhì)點之間用滿足胡克定律的彈簧連接。該模型通?；谝欢ǖ膸缀瓮負浣Y(jié)構(gòu)，形變過程是當一個質(zhì)點受外力作用時，產(chǎn)生的應力作用在其它相鄰的質(zhì)點上，同時通過彈簧傳遞力，帶動其相鄰的質(zhì)點運動。模型的建立比較簡單，形變的計算量較低，能夠很好滿足實時性的要求。

1.1軟組織的生物力學特征

對人體軟組織器官（如血管、肝臟、皮膚、氣管等）形變特性的研究屬于生物力學范疇［9］。隨著生物力學研究的發(fā)展，軟組織的生物力學特征通常表現(xiàn)為非線性、粘彈性、不均勻性和各向異性等［10］。因此，虛擬手術仿真系統(tǒng)中的軟組織建模必須考慮上述生物力學特性，并根據(jù)實際情況做相應的簡化處理。通常采用彈性模量、阻尼系數(shù)等物理量來表示軟組織的各種生物力學特征。

1.2彈簧—質(zhì)點模型建立

假設彈簧—質(zhì)點模型由n個質(zhì)點組成，在任意時刻t每個質(zhì)點Ni的運動方程都滿足牛頓第二定律，因此每個質(zhì)點Ni都具有如下的動力學微分方程（1）：

其中，mi表示第i個質(zhì)點質(zhì)量；xi表示質(zhì)點Ni的形變位移；ci表示阻尼系數(shù)；kij表示彈簧的彈性系數(shù)；lij表示質(zhì)點Ni和Nj間彈簧的初始長度，lij0表示形變后的長度；Fext表示質(zhì)點Ni所受的外力的綜合。

在彈簧—質(zhì)點模型模擬人體軟組織形變時，力的作用導致質(zhì)點發(fā)生位移。力的構(gòu)成，除了外力之外，還有彈簧的變形產(chǎn)生的內(nèi)力。不同的彈簧布局在仿真過程中會產(chǎn)生不同的力，質(zhì)點也會產(chǎn)生不同的位移，使模型產(chǎn)生不同的變形效果。當采用彈簧—質(zhì)點模型進行軟組織變形仿真實驗時，核心之一是建立合適的幾何拓撲結(jié)構(gòu)。

本文對人體的皮膚進行仿真實驗，采用面模型。因為人體的皮膚質(zhì)量分布相對均勻，其網(wǎng)格可均勻劃分。網(wǎng)格劃分的密度應根據(jù)仿真需求的精度和交互的實時性來確定。網(wǎng)格的拓撲結(jié)構(gòu)從不同的精度考慮采用圖1（b）、（c）、（d）三種方式［11］。

圖1　人體軟組織面模型拓撲結(jié)構(gòu)

為提高仿真精度，我們選用圖1（d）中拓撲結(jié)構(gòu)作為物理模型。彈簧可分為3類：結(jié)構(gòu)彈簧、剪切彈簧和彎曲彈簧，如圖2所示。

圖2　彈簧—質(zhì)點模型的三種彈簧類型

（1）結(jié)構(gòu)彈簧：質(zhì)點Ni，j和Ni+1，j之間和質(zhì)點Ni，j和Ni，j+1之間的彈簧，是防止質(zhì)點之間過度的拉壓形變。

（2）剪切彈簧：質(zhì)點Ni，j和Ni+1，j+1之間和質(zhì)點Ni，j+1和 Ni+1，j之間的彈簧，是阻止結(jié)構(gòu)的斜向過度變形，給結(jié)構(gòu)一個剪切剛度。

（3）彎曲彈簧：質(zhì)點Ni，j和Ni+2，j之間和質(zhì)點Ni，j和Ni，j+2之間的彈簧，是為了防止結(jié)構(gòu)過度的不自然彎曲。

2　軟組織模型算法分析

2.1動力學方程

基于上述彈簧—質(zhì)點模型拓撲結(jié)構(gòu)，結(jié)合受外力情況下模型的運動規(guī)律，對質(zhì)點建立動力學方程，通過求方程的解來實現(xiàn)軟組織的形變。

力學模型的建立基于以下兩點：（1）模型能夠達到實時的仿真效果；（2）在現(xiàn)有的硬件條件下能夠?qū)崿F(xiàn)相關算法的復雜度要求?，F(xiàn)有的基于物理模型的仿真算法一般采用兩種基本方法：力法和能量法。為了滿足仿真的實時交互需求，本文采用力法。以下是對彈簧—質(zhì)點模型的動力學分析［12，13］：

由牛頓第二定律可得，模型中每個質(zhì)點Ni的運動微分方程由它所受到的合力（內(nèi)力和外力組成）決定，對于質(zhì)點N 有：

其中，F(xiàn)ext（i）是質(zhì)點Ni所受的外力，F(xiàn)int（i）是質(zhì)點Ni所受的內(nèi)力（彈性力、阻尼力等）。根據(jù)胡克定律，質(zhì)點Ni所受彈簧的拉（壓）力如公式（3）所示：

其中，Kijs是質(zhì)點Ni和Nj間的彈性系數(shù)，Iij是該時刻彈簧的矢量，Xi和Xj是質(zhì)點Ni和Nj的位置。

質(zhì)點Ni和Nj之間的彈簧阻尼力如下：

其中，Kijd是質(zhì)點Ni和Nj之間的彈簧阻尼系數(shù)。

根據(jù)公式（3）和（4），可以得到質(zhì)點Ni所受的內(nèi)力：

其中，K是與質(zhì)點Ni連接的質(zhì)點總數(shù)。

一個由N個質(zhì)點組成的模型，所受的內(nèi)力由結(jié)構(gòu)力，彎曲力和剪切力組成；可以表述為：

其中，F(xiàn)struct、Fshearing、Fbending分別是與質(zhì)點Ni連接的彈簧產(chǎn)生的結(jié)構(gòu)力、剪切力和彎曲力，k1、k2、k3分別表示與質(zhì)點Ni相連的結(jié)構(gòu)彈簧、剪切彈簧和彎曲彈簧的數(shù)量。

一般情況下，質(zhì)點Ni所受到的外力包括重力和拉力，則受到的外力可以表述為：

對于上述的動力學微分方程，模型x并不是現(xiàn)實給出的，為了方便求解，需要對二階微分方程進行處理，將其轉(zhuǎn)化為一階常微分方程進行求解。則，每個質(zhì)點Ni的速度和位移都可以用以下公式表示：

2.2模型的數(shù)值分析

首先考慮一階常微分方程的初值求解問題：

其中f（ ）x，y為已知函數(shù)；y0為初值。

數(shù)值積分的原理就是，在區(qū)間［a，b］上取n+1個節(jié) a=x0＜x1＜x2＜…＜xn=b。 hi=xi-xi-1（1，2，…） 表示相鄰兩個節(jié)點的間距，稱為步長。這些hi可以不相等，我們?nèi)∠嗟瓤紤]。則將這些節(jié)點離散化，就可以將初值問題轉(zhuǎn)化成離散變量的相關問題。把相應問題的解yn作為y（xn）的近似值。這樣解得yn就是上述初值問題在節(jié)點xn的數(shù)值解，不同的離散方法產(chǎn)生不同的結(jié)果。

一階常微分方程初值問題的幾種常用的數(shù)值積分算法包括歐拉方法、中點方法以及四階龍格—庫塔方法等。對于質(zhì)點數(shù)一定的軟組織模型，形變的實時性和穩(wěn)定性取決于算法。主要從以下幾個方面來考慮：

（1）步長：表示為了達到一定數(shù)值穩(wěn)定性所需要的時間離散化程度，表明方法的數(shù)值穩(wěn)定性；

（2）迭代次數(shù)：反映了采用積分方法的計算量，影響總的計算速度；

（3）誤差：數(shù)值求解精度受誤差影響，包括截斷誤差、舍入誤差等；當確定步長后，微分方程的階次越高，截斷誤差越小；表1列出了以上三種數(shù)值算法的截斷誤差。

表1　三種數(shù)值算法的截斷誤差

從計算精度上來看，步長h確定的情況下，四階龍格—庫塔法的計算精度最高，但是執(zhí)行一步所需的計算量也最大，影響了仿真形變的實時性。只從單步積分來看，步長設置越小，模型離散化程度就越高，局部截斷誤差也就越小。但是離散化程度越大，同樣意味著在一定范圍內(nèi)需要的步數(shù)越多，不僅增加了計算量，還有可能引起舍入誤差的積累過大，導致精度變低。

3　碰撞檢測

碰撞檢測稱為接觸檢測或者干涉檢測，基本目的是確定兩個或兩個以上物體彼此之間是否發(fā)生接觸或穿透并計算出碰撞發(fā)生的位置。在虛擬手術中，需要碰撞檢測的場景主要有手術器械與人體器官之間、人體器官的自碰撞和手術器械的自碰撞。這些是整個虛擬手術系統(tǒng)中最重要的部分之一，決定著下一步的軟組織變形、切割、縫合等的實現(xiàn)。

如圖3所示，質(zhì)點t0時刻在軟組織外部，而在下一個時間步長t0+Δt時刻已經(jīng)在軟組織內(nèi)部，則質(zhì)點和軟組織在t0+Δt和t0+Δt之間的時刻發(fā)生了碰撞。

圖3　碰撞對象

目前常見的碰撞檢測算法主要有：層次包圍盒法［14］和空間分解法［15］。

層次包圍盒的核心思想是用體積略大而幾何特征簡單的包圍盒近似表述復雜的幾何對象，只對包圍盒相交時包裹的對象進行相交實驗。此外，通過引入樹狀層次結(jié)構(gòu)快速刪除不發(fā)生碰撞的部分，減少了不必要的相交實驗，提高檢測效率。層次包圍盒法應用較為廣泛，適用于復雜環(huán)境下的碰撞檢測。根據(jù)所采用的包圍體類型的不同可對層次包圍盒進行區(qū)分，比較典型的有沿坐標軸的包圍盒AABB、方向包圍盒OBB和包圍球等。

空間分解法的核心思想是將全部空間分解成體積相等的小單元格，只對同一單元格或相鄰單元格的幾何對象進行相交測試。比較典型的有八叉樹、BSP樹等。

虛擬手術仿真系統(tǒng)需要很好地滿足交互實時性的要求，因此，我們使用AABB層次包圍盒的改進算法作為本文虛擬手術仿真系統(tǒng)的碰撞檢測算法。

AABB方法［16］沿坐標軸的包圍盒是一種使用廣泛的碰撞檢測算法，包含幾何對象且邊平行于坐標軸的最小六面體。因此只需要六個標量就能描述一個AABB方法。

AABB樹是基于AABB方法自底至上動態(tài)更新生成的層次結(jié)構(gòu)二叉樹。AABB樹算法的核心是通過遍歷兩個碰撞物體的AABB樹來判斷該位置是否發(fā)生碰撞。步驟如下：

步驟1：構(gòu)造碰撞剛體和軟組織的AABB樹；

步驟2：檢測最大包圍盒是否相交；

步驟3：遞歸處理剛體和軟組織的AABB樹，如果發(fā)現(xiàn)子樹沒有發(fā)生相交，停止并得出結(jié)論沒有發(fā)生碰撞。如果發(fā)現(xiàn)子樹相交，則進一步處理它的子樹直到到達葉子節(jié)點，并最終得出結(jié)論。

軟組織形變主要是軟組織被提拉或者擠壓，所以軟組織形變的更新主要是AABB樹的更新。發(fā)生形變時，只有接觸點周圍的一些質(zhì)點發(fā)生了位移形變。相對于整個組織，變形只發(fā)生在局部區(qū)域，所以模型的拓撲結(jié)構(gòu)變化較小，更新時只對發(fā)生了形變的質(zhì)點進行更新。由于各個節(jié)點形變程度不同，不能用統(tǒng)一函數(shù)來描述。因此，先更新發(fā)生形變處節(jié)點的AABB包圍盒，再對整個包圍盒樹進行更新。父結(jié)點的包圍盒更新可通過兩個更新了的子節(jié)點表示。對于整個物體包圍盒的更新，如果檢測到上一層包圍盒沒有相交，其后的包圍盒就無需再更新［16］。

4　虛擬手術仿真結(jié)果

使用PC INTEL E3 3.10GHZ，8GBRAM，NVIDIA NVS 300顯卡計算機，以OPENGL圖形庫為基礎，使用VS2010和C++開發(fā)環(huán)境實現(xiàn)軟組織仿真的局部動態(tài)彈性形變。

由于彈簧—質(zhì)點模型是模擬皮膚表面的形變現(xiàn)象，只是一種近似仿真，而且模型的參數(shù)也無法選取真實軟組織的生物力學參數(shù)，所以通過大量實驗，對質(zhì)點質(zhì)量、彈簧的彈性系數(shù)和阻尼系數(shù)等進行對比后選取如下：時間步長0.01，質(zhì)點質(zhì)量0.01，重力加速度9.8，彈簧的結(jié)構(gòu)彈簧、剪切彈簧和彎曲彈簧的彈性系數(shù)分別為2，2，1.2；阻尼系數(shù)分別為0.4，0.4，0.24。虛擬體彈簧的彈性系數(shù)為表面彈簧彈性系數(shù)的兩倍。

以平面模型為例進行仿真，如圖4所示。分別選取了225個質(zhì)點，450個面片（a）和625個質(zhì)點1250個面片（b）的面模型進行試驗。

圖4　不同質(zhì)點的面模型

以圖4（a）為模型，對固定形變區(qū)域和動態(tài)形變區(qū)域進行對比。當外力較小時，若形變區(qū)域范圍很大，則計算形變花費的時間多，所以我們采用固定形變區(qū)域，如圖5（a）所示；當外力大時，若形變區(qū)域范圍很小，則形變效果較差，所以我們采用動態(tài)形變區(qū)域，如圖5（b）所示。

圖5　形變效果圖

以100個時間步長為例，虛擬一個較小的外力，選用圖5（b）所示模型。在進行形變計算時，選取一個合適的臨界值，頂點的形變大于該值時，激活下一個領域的頂點。若在某個時間步長內(nèi)，碰撞點在四個方向的形變量均小于該臨界值，則退出循環(huán)。在相同的外力作用下，不同的臨界值不僅影響形變區(qū)域的大小，同時還影響時間步長。臨界值設置越小形變越精細，相應的形變計算步驟越多，時間越長。本文選取臨界值為0.001，虛擬外力為5N的情況下，分別用歐拉法，中點法，四階龍格—庫塔法進行仿真實驗，計算時間如表2所示，仿真效果如圖6所示。

表2　網(wǎng)格模型三種積分方法計算時間比較

實驗結(jié)果表明網(wǎng)格質(zhì)點越多，實時性越差。從3種不同的積分效果去考慮，可以發(fā)現(xiàn)，實驗的真實性越來越好，但實時性越來越差。在網(wǎng)格精度相對較高的情況下，歐拉法和中點法能夠達到很好的仿真效果且計算時間能夠達到實時性的需求，但四階龍格—庫塔法的計算時間明顯提高。綜合比較后，在滿足實時性的同時選擇中點法，從而達到較好的仿真精度和效果。

圖6　實驗結(jié)果圖

5　結(jié)論

本文主要針對虛擬手術仿真過程中實時性和真實性兩方面的需求，采用較能反映軟組織形變“體特征”的質(zhì)點—彈簧體模型實現(xiàn)了軟組織形變，并在此基礎上提出了更加精確的拓撲結(jié)構(gòu)。為了增強形變過程中的穩(wěn)定性，設置了結(jié)構(gòu)彈簧、剪切彈簧和彎曲彈簧；基于模型形變的動力學方程采用歐拉方法，中點法和四階龍格-庫塔法作為單個質(zhì)點的運動算法，并通過對虛擬軟組織器官的提拉和按壓手術實驗。結(jié)果表明，本文提出的方法能夠很好地滿足模型的精度和仿真實時性要求，增強了形變的仿真體積感，適合于軟組織形變較大時的實時仿真。

［1］ 孫連榮，姜元章，任玲.虛擬現(xiàn)實技術及其在高校教學中的應用［J］.石油教育，2003（5）：41-44.

［2］ 戴雯惠.虛擬現(xiàn)實技術的應用現(xiàn)狀及發(fā)展［J］.信息技術，2006，35（6）：170-174.

［3］Santanu M，Amit R.Simulation of hip fracture in sideways fall using a 3D finite element model of pelvis-femur-softtissuecomplexwithsimplified representation of whole body［J］.Medical Engineering&Physics，2007，29：1167-1178.

［4］ Zhong Yongmin，Bjian S，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］ 陳衛(wèi)東，趙成龍.虛擬手術中軟組織形變建模及力反饋算法研究［J］.中國生物醫(yī)學工程學報，2013，32（1）：113-118.

［6］ 王瑞藝.虛擬手術中力反饋的研究和實現(xiàn)［D］.鄭州：華北水利大學，2014.

［7］Li Yandong，Zhu Ling.Modeling and simulation of softtissuedeformationbasedonlocaldynamic model［J］.Computer Science，2013，40（10）：283-288.

［8］ Wang P，Becker A A，Jones I A，et al.Virtual reality simulation of surgery with haptic feedback based on the boundary element method［J］.Computers and Structures，2007，85：331-339.

［9］馮元禎.生物力學［M］.北京：科學出版社，1983：124-158.

［10］Meredith M，Maddock S.Real-time inverse kinematics：the return of the jacobian［J］.Department of ComputerScience，TheUniversityofSheffield. 2004：51-60.

［11］ 喬冰.基于質(zhì)點—彈簧模型的人體軟組織形變技術的研究［D］.哈爾濱：哈爾濱工程大學，2008.

［12］ 褚蓮娣.基于質(zhì)點一彈簧模型的布模擬方法［J］.嘉興學院學報，2002，14（3）：57-60.

［13］ 潘振寬，高波.手術方針中基于質(zhì)點一彈簧模型的人體組織變形仿真［J］.青島大學學報，2003，18（3）：9-14.

［14］ 朱元峰，孟軍，謝光華，等.基于復合層次包圍盒的實時碰撞檢測研究［J］.系統(tǒng)仿真學報，2008，20（2）：372-377.

［15］ 鄒益勝，丁國富，許明恒，等.實時碰撞檢測算法綜述［J］.計算機應用研究，2008，25（1）：8-12.

［16］ 趙成龍.虛擬手術中軟組織建模與力反饋算法研究［D］.秦皇島：燕山大學，2013.

An Improved Method on Soft Tissue Deformation in Virtual Surgery

HE Wei，WANG Weiping，SHI Weili，MIAO Yu，HE Fei，YANG Huamin
（School of Computer Science and Technology，Changchun University of Science and Technology，Changchun 130022）

Virtual surgery provides the surgery simulation for users on vision and force.It applies computer virtual reality in modern medicine，uses medical image data and reconstructs virtual human soft tissue models.And the simulation and modeling of soft tissue deformation has become one of the most important research contents in virtual surgery system.An improved mass spring model is proposed in this paper for the soft tissue deformation simulation.Bending spring，shear spring and structure spring are added onto the surface model of quadrilateral mesh topological structural for soft tissues，which will make the surface model having volume feature.This model not only has the advantages of traditional spring mass model such as rapid modeling，simple principle and quick simulation，but also has the ability to control the deformation regions，which improves the accuracy of the simulation model，and meets the requirements of the realness and stability of interactive system.Experimental results show that the method is suitable for the real-time simulation of soft tissues with large deformation.It can effectively enhance the volume of deformation simulation and meet the requirements of real-time interactive systems.

virtual reality technology；mass spring model；soft tissue deformation；collision detection

TP317.4

A

1672-9870（2015）06-0118-05

2015-10-30

何?。?978-），女，博士研究生，講師，E-mail：hw@cust.edu.cn

楊華民（1963-），男，教授，博士生導師，E-mail：yhm@cust.edu.cn