基于MLS-MWS混合模型的軟組織建模

陳衛(wèi)東 劉 波 朱奇光*

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

基于MLS-MWS混合模型的軟組織建模

陳衛(wèi)東1,2劉 波1朱奇光1,2*

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

針對無網(wǎng)格模型計算量大和邊界處理繁瑣的問題，提出一種基于改進最小二乘近似函數(shù)的無網(wǎng)格強弱勢混合模型(MLS-MWS)。首先，基于場理論將模型劃分為內部域和外部域，在內部域中，不與自然邊界相連接的節(jié)點通過無網(wǎng)格強勢法來構建平衡方程，該方法無數(shù)值積分，從而達到計算量的減少；在外部域中，通過無網(wǎng)格局部弱勢法來弱化邊界問題。其次，在強勢方法中，需利用最小二乘近似式表示位移的二階導數(shù)，針對其計算量大和效率低的問題，提出改進移動最小二乘近似函數(shù)和導數(shù)，通過將節(jié)點值進行兩次MLS導數(shù)插值，求出該節(jié)點的二階導數(shù)，進一步簡化計算；最后采用PHANTOM觸覺交互設備，對肝臟模型進行拉伸與擠壓的形變模擬實驗。 實驗結果表明，MLS-MWS混合模型相比無網(wǎng)格伽遼金模型，單步執(zhí)行時間從17.3 ms減少到13.8 ms，頻率從56.80幀/s提高到72.46幀/s。該混合模型結合無網(wǎng)格強弱勢法和雙重無網(wǎng)格配點法，達到優(yōu)勢互補，有效提高模型的實時性。

軟組織建模；無網(wǎng)格法；雙重配點無網(wǎng)格法；MLS-MWS模型

引言

虛擬手術仿真系統(tǒng)在醫(yī)學教學與實踐中占據(jù)重要地位，通過該仿真系統(tǒng)，醫(yī)生可以進行大量實踐來豐富自己的臨床經(jīng)驗[1]。軟組織形變模型的構建是虛擬手術仿真系統(tǒng)的核心步驟，軟組織模型運用計算機圖形學與處理、可視化技術和生物醫(yī)學等相關學科知識，獲得可視化圖像和力學反饋數(shù)據(jù)。通過對所獲得數(shù)據(jù)的處理，可以模擬出軟組織模型的組織形態(tài)特性與力學形態(tài)特性，因此軟組織建模研究具有重要的理論依據(jù)與實際意義。

軟組織模型通常分為彈簧質點模型和無網(wǎng)格模型。Vollinger等提出了改進的質點彈簧模型[2]，模型參數(shù)的選擇是該模型的難題，在某些特殊的仿真環(huán)境中，系統(tǒng)會因參數(shù)選擇不當而引起震蕩的不穩(wěn)定現(xiàn)象。閆雒恒等提出了一種改進四面體網(wǎng)格的彈簧質點模型[3]，通過從新定義彈性力學方程來描述形變過程中應力的變化，提高了仿真的逼真度和實時性，但該模型還存在穩(wěn)定性差的問題。陳衛(wèi)東等提出了一種簡單的彈簧質點混合模型[4]，通過選擇有效的節(jié)點，計算各層形變區(qū)域，相比傳統(tǒng)的彈簧質點模型提高了逼真度，但穩(wěn)定性相對較差。宋愛國等提出了一種高效的層狀菱形鏈接連接模型[5]，通過搭建雙通道力/觸覺交互的虛擬肺仿真系統(tǒng)，提高了逼真度和計算速度，但是該模型的穩(wěn)定性相對較差。Francois等提出了運用超彈性質點連接來快速計算軟組織的形變與實時仿真[6]，其仿真的逼真度比較高，但模型的穩(wěn)定性相對較差，容易變形。Zhang等采用全拉格朗日近似的無網(wǎng)格伽遼金方法建立軟組織模型[7]，通過使用大量節(jié)點來建立形函數(shù)，在大形變中可以顯著延遲網(wǎng)格畸變，提高了仿真的逼真度，但該方法計算過程復雜。Steinemann等提出了一種基于無網(wǎng)格和有限元方法的混合模型[8]，但該方法計算復雜，計算效率低。劉雪梅等提出了無網(wǎng)格弱勢法與質點彈簧模型[9]，將形變區(qū)域分為手術區(qū)域和其他區(qū)域，手術區(qū)域采用無網(wǎng)格弱勢法，具有較好的精度與穩(wěn)定性，但該模型的缺點是計算量大。許少劍等提出了一種最小二乘法的無網(wǎng)格方法對模型進行研究[10]，減少了畸形網(wǎng)格生成對虛擬手術仿真系統(tǒng)所產(chǎn)生的不穩(wěn)定現(xiàn)象，無需網(wǎng)格，減少了計算量，但是如何處理邊界問題沒有得到解決。

本研究針對無網(wǎng)格模型計算量大和邊界處理繁瑣的問題，提出了一種基于改進最小二乘近似函數(shù)的無網(wǎng)格強弱勢混合模型。首先，基于場理論將模型劃分為內部域和外部域，在內部域中所有節(jié)點通過無網(wǎng)格強勢方法建立平衡方程，減小計算量，在外部域中利用無網(wǎng)格弱勢方法弱化邊界問題；其次，在強勢方法中，利用雙重無網(wǎng)格配點法改進近似函數(shù)二階導數(shù)的求導過程，進一步簡化計算。該混合模型繼承了無網(wǎng)格強弱勢法和雙重無網(wǎng)格配點法的優(yōu)點，從而有效提高了模型的實時性。

1 模型構建方法

在建模過程中，無網(wǎng)格算法只需對場內節(jié)點建立離散方程，從而構建出模型的形狀，但是傳統(tǒng)的無網(wǎng)格算法存在計算量大和邊界問題，因此本研究提出了一種改進最小二乘近似函數(shù)的無網(wǎng)格強弱勢法，包括內部域與外部域的劃分和利用改進的最小二乘函數(shù)建立無網(wǎng)格強弱勢離散方程。MLS-MWS方法的核心思想是：基于場理論知識，將場內任意節(jié)點分別通過無網(wǎng)格強勢和無網(wǎng)格局部弱勢來建立離散方程。

改進最小二乘近似函數(shù)的無網(wǎng)格強弱勢方法具有相同的離散形式，因此在整體離散域Ω上，MLS-MWS的離散方程為

(1)

式中，n為求解域的節(jié)點數(shù)，u為形函數(shù)的函數(shù)值，K為系統(tǒng)的剛度矩陣，F(xiàn)載荷向量。

如圖1所示，求解域Ω和邊界Γ由一些任意分布的節(jié)點組成，Ωq是積分域的任意節(jié)點。如果Ωq不與自然邊界相連接，則采用無網(wǎng)格強勢方法；反之，采用無網(wǎng)格局部弱勢方法[11]。將求解域劃分為內部求解域Ωin和邊界域ΩB，有

(2)

式中,CI表示以場點為場域圓心、半徑為κdIS的圓，dIS為節(jié)點支持域的半徑，κ為改變圓形場域的控制參數(shù)。

利用調整控制參數(shù)κ，通過改變內部域Ωin與邊界域ΩB的節(jié)點數(shù)，可以有效改變內部域與外部域的大小。

圖1 MLS-MWS方法示意圖Fig.1 A schematic diagram of MLS-MWS method

在軟組織建模時，時間的快慢將嚴重影響虛擬仿真系統(tǒng)整體性能。Γu、Γt分別表示位移邊界和面力邊界，考慮強勢邊界真值問題時，其控制方程表示如下：

在Ω中，平衡方程為

(3)

在Γu和Γt上，邊界條件分別為

(4)

(5)

1.1 無網(wǎng)格內部域的算法

在內部域Ωin中，數(shù)值積分和邊界問題不做考慮，只需要滿足平衡方程。平衡方程之間采用強勢方法離散，使得內部域Ωin中的節(jié)點精確滿足平衡方程，位移μ的二階導數(shù)利用最小二乘法近似式表示，即可得到Ωin中節(jié)點的離散方程，有

(6)

在z=0時，KIJ,in、uJ和FI,in分別為

(8)

式中，KIJ,in、FI，in分別為系統(tǒng)剛度矩陣[12]與載荷向量[12]，E、μ分別表示彈性模量和泊松比，φJ,ij是通過式(33)得出的。

1.2 無網(wǎng)格邊界域的算法

盡管無網(wǎng)格局部弱勢法已經(jīng)向理想無網(wǎng)格法的發(fā)展邁出了重要一步，尤其在復雜形狀的邊界問題上，數(shù)值積分仍然是一項繁瑣的任務。因此，應該盡量減少使用數(shù)值積分。利用邊界區(qū)域場點少，沒有切割和縫合等復雜問題。根據(jù)無網(wǎng)格局部弱勢法，可以弱化邊界問題。在邊界域中，根據(jù)彈性力學中的本構關系、幾何方程以及應力-面力關系，可以求得應力δij和面力ti的移動最小二乘近似式，可以得到ΩB中節(jié)點的離散方程，有

(9)

在z=0時，KIJ,B和FI,B[12]分別為

(10)

(11)

其他相關矩陣表達式如下：

(12)

(13)

(14)

(15)

在邊界域ΩB中，子域通過加權殘值法，將本質邊界條件通過罰函數(shù)引入到殘值方程，得到

(16)

式中：α為罰因子；Ωq為局部積分域；Γqu為局部積分域邊界?Ωq與位移邊界Γu相交的部分；γi為權函數(shù)，任意形狀的局部積分域均可使用，本研究的局部積分域采用圓形域。

在局部積分域中，邊界Γq可以包括三部分：內部邊界Γqi、局部積分域邊界?Ωq與位移邊界Γu相交的部分Γqu、自然邊界Γqt。文獻[13]對式(16)利用散度定理進行處理，可以得到邊界域ΩB中控制方程的局部弱形式，有

(17)

1.3 移動最小二乘近似函數(shù)和導數(shù)

假設待求函數(shù)u(x)在求解域Ω中N個節(jié)點XI(I=1,2,3，…，n)的函數(shù)值已知，即uI=u(xI)，則u(x)局部近似函數(shù)為

(18)

二維空間單項式基函數(shù)中的線性基和二次基分別為

(19)

(20)

三維空間單項式基函數(shù)中的線性基和二次基分別為

(21)

(22)

(23)

令J取最小值，即

(24)

通過式(24)得

(25)

(26)

式中

(27)

由式(26),可得待定系數(shù)向量α(x),即

(28)

將式(28)代入式(18)，得

(29)

其中，形函數(shù)

(30)

在配點無網(wǎng)格法中，利用近似函數(shù)的二階導數(shù)來表示位移的二階導數(shù)，由于其計算兩次形函數(shù)導數(shù)，使得計算量大、效率低。本研究利用雙重無網(wǎng)格配點法[14]來進行計算。首先引入兩組點，即節(jié)點xI和計算點xke。

(31)

(32)

式中

(33)

利用改進的最小二乘近似函數(shù)求二階導數(shù)，降低了計算的復雜度，避免出現(xiàn)病態(tài)方程，使方程的穩(wěn)定性增強，且通過對節(jié)點進行兩次MLS導數(shù)插值，使結果具有光滑的連續(xù)性，在軟組織形變過程中更能逼真地進行形變。

2 模型實驗

以肝臟模型為例，其泊松比為0.35。采用SensAble科技公司的PHANTOM觸覺交互設備，基于酷睿i5處理器，4GB內存，NVIDIAGeForceGTX960M顯卡臺式電腦。在Windows7操作系統(tǒng)下，采用VisualC++ 6.0 開發(fā)環(huán)境，并結合OpenGL三維圖形標準，搭建了虛擬軟組織實時力反饋交互系統(tǒng)。

在實驗中，分別用無網(wǎng)格伽遼金方法和MLS-MWS方法構建肝臟模型，并在相同場點施加相同載荷的拉力和壓力，然后計算出各自的位移大小，進行對比實驗。同時，通過設置參數(shù)κ(無明確的物理意義)的大小來控制內部域與邊界域大小。由于MLS-MWS法的計算效率隨著邊界域節(jié)點數(shù)占總節(jié)點數(shù)的比例值減小而提高，參數(shù)κ通過幾組數(shù)據(jù)測試來進行設置，當0.5≤κ≤2.0時，使得MLS-MWS法減少了計算復雜度，提高了實時性。

在實驗中，無網(wǎng)格伽遼金模型和MLS-MWS模型在質點個數(shù)為293時，施加力F=0.3N，分別對無網(wǎng)格伽遼金模型和MLS-MWS模型進行平行對比實驗。

3 形變仿真結果

在虛擬手術仿真系統(tǒng)中，逼真度與實時性是決定仿真成功的指標。視覺逼真度是描述模型仿真的一個能力，描述能力越強，模型的逼真度越高。圖2為肝臟模型的拉伸仿真，圖3是肝臟模型壓迫仿真圖。

圖2 拉力作用下的受力Fig.2 The liver rendering under action of pull force

圖3 壓力作用下的受力Fig.3 The liver rendering under action of pressure force

在實驗中，分別運用無網(wǎng)格伽遼金法[9]和所提出的MLS-MWS方法構造肝臟模型，并在相同場點施加相同載荷的拉力和壓力，然后計算出各自的位移大小，進行對比實驗。圖4、5分別是拉力與壓力作用下的位移曲線。通過仿真實驗結果可知，在施加相同拉力或壓力下，MLS-MWS混合模型與無網(wǎng)格模型的位移誤差為±0.12。從數(shù)值上可以看出，該混合模型的精度與無網(wǎng)格模型的精度基本相同，可以說明MLS-MWS混合模型有效。

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

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

為了驗證MLS-MWS法的計算效率，在實驗中，當κ=0.8、1.0、1.5、2.0時，分別對無網(wǎng)格伽遼金方法和MLS-MWS方法進行了平行對比實驗，只有當κ=1.0時效果最好。在實驗中，當κ=1.0時，內部域節(jié)點數(shù)為232，邊界域節(jié)點數(shù)為61，分別施加不同載荷。由表1的實驗數(shù)據(jù)可以得出，MLS-MWS混合模型的單步執(zhí)行時間比無網(wǎng)格模型的單步執(zhí)行時間少了3.8ms，頻率提高了15.64幀。結果說明，MLS-MWS混合模型的計算速度快，減少了計算時間，實時性優(yōu)越。同時，在不影響精度的前提下，MLS-MWS混合模型單步執(zhí)行時間比無網(wǎng)格模型的單步執(zhí)行時間約提高了20%，頻率提高近15.64幀，進一步說明該混合模型實時性和逼真度優(yōu)于無網(wǎng)格模型。

表1 無網(wǎng)格伽遼金法與MLS-MWS法計算效率對比

4 討論和結論

虛擬手術是虛擬現(xiàn)實技術在醫(yī)學方面上的一個重要應用，它利用醫(yī)學圖像數(shù)據(jù)和計算機圖形學構建出軟組織模型，模擬出虛擬的手術環(huán)境，并利用觸覺交互設備與之進行交互。軟組織形變模型的構建是虛擬手術仿真系統(tǒng)的核心步驟，由于軟組織模型在生物黏彈性、各向異性等方面與軟組織器官有一定差距，因此構建一個實時性與逼真度高的軟組織模型至關重要。隨著虛擬手術系統(tǒng)被引入臨床教學中，醫(yī)生通過大量的訓練并結合臨床經(jīng)驗，減少了手術的風險系數(shù)。因此，通過虛擬手術仿真系統(tǒng)可以減少手術的風險，降低培養(yǎng)一個優(yōu)秀外科醫(yī)生的費用。同時，通過典型病例并利用該仿真系統(tǒng)對醫(yī)生進行培訓，對提高醫(yī)生的操作能力、質量和醫(yī)學水平有著深遠的意義。

本研究提出MLS-MWS混合模型，通過調節(jié)參數(shù)κ的大小來控制內部域與外部域的范圍。在內部域中，采用無數(shù)值積分的無網(wǎng)格強勢方法，減少了計算的復雜程度；在外部域中，采用無網(wǎng)格局部弱勢方法，使形狀函數(shù)具有Kroneckerdelta函數(shù)性質，可以有效處理邊界問題，使計算結果準確穩(wěn)定。同時，在強勢方法中，利用雙重無網(wǎng)格配點法改進近似函數(shù)二階導數(shù)的求導過程，進一步簡化了計算。實驗結果表明，MLS-MWS混合模型單步執(zhí)行時間比無網(wǎng)格模型的單步執(zhí)行時間約提高了20%，頻率提高了近15.642幀，極大地提升了醫(yī)療手術過程中軟組織的實時形態(tài)及其逼真程度，降低了手術風險，提高了手術成功幾率。

本研究提出了一種MLS-MWS混合模型，減少了建模過程中的計算量，弱化了邊界問題，大幅地提高了模型的實時性和逼真程度。但是，在無網(wǎng)格軟組織建模中，邊界域節(jié)點數(shù)占總節(jié)點數(shù)的比重值越小，軟組織模型的計算效率越高。本研究提出MLS-MWS混合模型，通過控制參數(shù)κ的大小來控制邊界域的節(jié)點個數(shù)，但是κ值的選取需進行大量實驗，從而導致該方法選取的κ值存在一定誤差，因此如何有效控制內部域與外部域的范圍是下一步研究的重點工作。

[1] 吳崢, 謝叻, 馬浩博. 虛擬手術實時物體碰撞檢測和軟組織變形研究[J]. 計算機仿真, 2010, 27(2):255-259.

[2] V?llinger U, Stier H, Priesnitz J, et al. Evolutionary optimization of mass-spring models[J]. Cirp Journal of Manufacturing Science & Technology, 2009, 1(3):137-141.

[3] 閆雒恒. 基于改進彈簧振子模型的軟組織形變仿真[J]. 計算機仿真, 2012(7):330-333.

[4] 陳衛(wèi)東, 陳攀攀, 朱奇光. 基于改進彈簧-質點模型的形變建模及力反饋算法研究[J]. 生物醫(yī)學工程學雜志, 2015，32(5)：989-996.

[5] 張小瑞, 孫偉, 宋愛國,等. 雙通道力/觸覺交互的虛擬肺手術仿真系統(tǒng)[J]. 儀器儀表學報, 2012, 33(2):421-428.

[6] Goulette F, Chen Zhuowei. Fast computation of soft tissue deformations in real-time simulation with Hyper-Elastic Mass Links[J]. Computer Methods in Applied Mechanics & Engineering, 2015, 295:18-38.

[7] Zhang GY, Wittek A, Joldes GR, et al. A three-dimensional nonlinear meshfree algorithm for simulating mechanical responses of soft tissue[J]. Engineering Analysis with Boundary Elements, 2014, 42(42):60-66.

[8] Steinemann D, Otaduy MA, Gross M. Splitting meshless deforming objects with explicit surface tracking[J]. Graphical Models, 2009, 71(6):209-220.

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

[10] 許少劍, 陳國棟. 基于移動最小二乘法的無網(wǎng)格肝臟形變研究[J]. 計算機仿真, 2015, 32(4):408-413.

[11] 肖毅華, 胡德安, 韓旭. 彈性靜力問題的無網(wǎng)格弱-強形式結合法[J]. 計算力學學報, 2010, 27(5):764-769.

[12] Gu YTGR. A meshfree weak-strong (MWS) form method for time dependent problems[J]. Computational Mechanics, 2005, 35(2):134-145.

[13] 龍述堯, 姜琛, 鄭娟. 三維彈性靜力問題的無網(wǎng)格局部Petrov-Galerkin法[J]. 湖南大學學報(自然科學版),2013, 40(12):55-61.

[14] 張雄, 劉巖, 馬上. 無網(wǎng)格法的理論及應用[J]. 力學進展, 2009, 39(1):1-36.

Soft Tissue Modeling Based on MLS-MWS Mixing Method

Chen Weidong1,2Liu Bo1Zhu Qiguang1,2*

1(InsituteofInformationScienceandEngineering,YanshanUniversity,Qinhuangdao066004,China)2(TheKeyLaboratoryforSpecialFiberandFiberSensorofHebeiProvince,Qinhuangdao066004,China)

Targeting to problem of the large computation amount of mesh-free model and tedious boundary processing, MLS-MWS based on the improved least square approximation was proposed in this work. At first, the model was divided into the internal domain and external domain based on the field theory. In the internal domain, nodes far from the natural boundary constructed the equation of equilibrium with the mesh-free method. The method has no numerical integration, which may reduce the computation amount. In the external domain, the boundary problem was weakened through the mesh-free partial weak method. Secondly in the strong method, second-order derivative of the displacement should be indicated with the least square approximation equation. Targeting tothe large computation amount of mesh-free model and low efficiency, it was proposed to improve the moving least square approximation function and derivate, and MLS derivative interpolation was carried out for the node value, for working out the second-order derivative and further simplifying the computation. Eventually, deformation simulation experiment of stretching and squeezing the liver model was carried out with the PHANTOM touch interactive equipment. According to the results, the single-step execution time of MLS-MWS combined model was reduced to 13.8 ms from the original 17.3 ms of the mesh-free model, while the frequency was improved from 56.80 frames/s to 72.46 frames/s. The mixed model integrated advantages of the mesh-free weak-strong method and dual mesh-free point-allotting method, which improved the instantaneity of the model effectively.

soft tissue modeling; meshless method; double match point meshless method; MLS-MWS hybrid model

10.3969/j.issn.0258-8021. 2016. 06.009

2016-01-12， 錄用日期:2016-05-19

國家自然科學基金(61201112)；河北省自然科學基金(F2016203245)；河北省普通高等學校青年拔尖人才計劃(BJ2014056)

R318

A

0258-8021(2016) 06-0699-06

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