基于彈簧-質點模型的基礎上提出一種新的可變菱形拓撲結構模型,通過改變菱形邊彈簧的長度及彈簧系數、初始夾角,可實現該模型對不同器官的形變模擬。此外在模型中加入虛擬彈簧提供約束力,避免外力過大出現超彈性現象,并用于形變過程中反饋力的計算。在力反饋形變算法中,篩選有效的移動質點,并把這些移動質點形成區域面積的改變與作用力聯系起來,計算量小,從而簡化了傳統的彈簧-質點模型的力反饋算法。基于此模型采用PHANTOM觸覺交互設備進行受力形變模擬實驗。實驗結果表明,該模型結構簡單易于實現,采用的力反饋形變算法,降低了形變計算的次數,提高了形變的實時性,具有更逼真的形變效果。
引用本文： 陳衛東, 陳攀攀, 朱奇光. 基于改進彈簧-質點模型的形變建模及力反饋算法研究. 生物醫學工程學雜志, 2015, 32(5): 989-996. doi: 10.7507/1001-5515.20150176 復制
引言
虛擬手術仿真系統是虛擬現實技術在醫學領域的一種應用[1]，在手術前對手術過程涉及的人體軟組織進行虛擬手術仿真具有重大的學術、經濟和社會意義，采用快速準確的形變模型和形變算法是整個手術仿真系統的關鍵，也是當前虛擬現實技術研究的難點和熱點。
目前虛擬建模已有多種方法[2]，其中主要以彈簧-質點模型、有限元模型[3]、邊界元模型等為主。各種建模方法都有其特有的優點和缺陷，相比之下，雖然彈簧-質點模型在大形變時同樣會出現失真等問題，但是彈簧-質點模型計算量較少，能夠很好地滿足仿真的實時性要求。而在一個虛擬手術系統中，實時性是衡量系統優劣最為關鍵的指標之一。因此，多數模型是在基于彈簧-質點模型[4-6]的基礎上進行改進的。吳涓等[7]的基于同心圓分割的彈簧-質點模型及實時力覺相應算法雖然計算速度快，滿足力反饋的實時性要求，但各點計算不一致，不能進一步擴展。張小瑞等[8]的層狀菱形鏈連接模型雖然力反饋和變形計算具有較高精度，但是菱形邊的長度在受力時無法改變，影響形變效果的真實性。混合模型[9]的提出，雖然滿足力反饋計算的實時性，但是混合模型的過渡區不易確定。因此如何建立穩定可靠，同時又計算簡便、實時性好的物理模型，仍然是虛擬現實領域亟待解決的一個問題[10]。
為了達到更好的模擬效果，在基于彈簧-質點模型的基礎上進行改進，提出了一種基于可變菱形的彈簧-質點模型，能實現用一種模型對不同的器官進行形變模擬的目的，提高形變效果的真實性，同時通過在模型中加入虛擬彈簧，用于形變過程中力反饋的計算，降低了形變計算的次數，得到一個更加逼真有效的形變效果，更好地滿足了交互系統的真實性、穩定性和實時性的要求。加入虛擬彈簧的同時在兩質點間也加入一個無質量的不可拉伸的輕繩，在物體受外力作用時起限定作用，當質點受力大于輕繩能承受的力的上限時，會出現輕繩的斷裂，兩質點脫離，從而模擬了現實中柔性物體的破裂現象。
1 軟組織形變模型原理及結構
1.1 軟組織形變模型原理
由于軟組織在受力狀態下，受力面的位置、施加力的大小和方向都在實時變化，給系統的穩定性和受力分析帶來了困難，而傳統的基于表面網格拓撲結構的模型，在模擬軟組織形變過程中穩定性較差，不能很好地表現軟組織的生物力學特性[11]。針對該問題，在基于層狀菱形鏈連接模型的基礎上，提出了一種改進的基于可變菱形拓撲結構的彈簧-質點模型。
由于人體軟組織器官的多樣性和復雜性，軟組織器官通常表現為不均勻性、各向異性、準不可壓縮性等材料性質，雖然目前已有很多算法模型并在實際應用中取得了良好的效果，但是在穩定性和可靠性方面還存在不少問題，建立完全真實的模型很難實現。基于可變菱形的拓撲結構，針對不同的組織器官及其本身特性，考慮到軟組織的不均勻性，例如肺癌早期，肺表面會出現比正常肺組織更硬的硬塊，在相同外力作用下肺硬塊比正常肺組織對應的形變量要小一些，因此可以設置模型中的參數，如肺硬塊菱形邊的彈簧長度長一些、彈簧系數大一點、彈簧間的夾角大一些，那么在相同作用力下，產生的形變量就會小一些，以更好地表現形變的細節，利于辨別不同的組織，為手術的成功提供一定的保障。
1.2 形變模型結構
設有虛擬彈簧的菱形結構的模型如圖 1所示。該模型由多層菱形結構單元組成，每個菱形結構單元由4個質點組成，各質點之間通過彈簧相連，施加在軟組織表面上的力，通過與其直接交互的結構單元傳遞到周圍相連的各個結構單元，由于結構單元各質點之間由彈簧連接，若彈簧一端質點發生移動，另一端質點受到彈簧力的作用位置也隨著改變。
圖1
設有虛擬彈簧的菱形結構模型
Figure1.
Diamond structure model with the virtual spring
在對整體模型進行受力分析時，物體的表面形變就是各質點的彈簧形變的疊加，而物體表面的接觸力就等效為各質點對應的虛擬彈簧彈力的合力，因此在進行受力分析時不需要分析整個虛擬軟組織，只涉及發生形變的的局部區域內的局部質點。在虛擬器械與受力柔性體接觸過程中，按照質點的連接順序劃分層次，以受力網格作為第一層，對柔性體表面每一個質點以接觸點為中心按先內后外的順序分層，這樣第一層就有4個質點，4根彈簧。每增加一層質點數增加8個，彈簧數增加16個，以此類推，第m層質點的個數為4+8(m-1)，彈簧的個數為4+16(m-1)，形成類似m層同心環狀的質點分布結構。
為了驗證本文模型拓撲結構的優越性，在相同質點數目下，將本文模型的拓撲結構與常用模型的拓撲結構進行比較，常用模型的拓撲結構如圖 2所示。多組相同質點數目下彈簧數量的比較結果如表 1所示。
圖2
常用模型的拓撲結構
Figure2.
Topology of the commonly used model
表1
彈簧數量對比
Table1.
Comparison of spring number
從表 1中數據可知，改進菱形結構分別與傳統三角形結構、改進正六邊形結構相比，在質點數量均為144時，彈簧數量各減少了32％、6%左右，在質點數量均為324時，彈簧數量各減少了33％、4%左右，在質點數量均為676時，彈簧數量各減少了33％、4%左右，因而建模復雜度低，實時性較好；而與層狀菱形鏈結構相比，在質點數量均為144時，彈簧數量增加了37.5％左右，在質點數量均為324時，彈簧數量增加了42%左右，在質點數量均為676時，彈簧數量增加了44%左右，因而穩定性較好，計算精度較高。因此，改進的拓撲結構在計算精度和實時性兩方面進行了折中，能較好地描述軟組織形變的物理特性。
2 模型形變的力反饋計算
2.1 模型形變的力反饋算法
假設質點均是理想質點，連接質點的彈簧也是理想彈簧，彈簧的阻尼力忽略不計，彈簧初始長度的變化量為零，彈簧彈力與形變量滿足胡克定律，在外力撤消后還能恢復到原狀，模型中的每一個質點的運動均滿足牛頓第二定律。在變形計算中，按照先整體后局部的思想進行。當外力F作用于任一質點時，物體的表面形變就是各質點的彈簧形變的疊加，而物體表面的接觸力就等效為各質點對應的虛擬彈簧彈力的合力，當形變層數i=m時，計算公式為
| $ F = 4\sum\limits_{i = 1}^{i = m} {{f_{{N_i}}}} + 8\sum\limits_{i = 1}^{i = m} {\left( {n - 1} \right)f{'_{{N_i}}}} , $ |
式中fNi為位于菱形頂點處質點的虛擬彈簧的彈力，fNi′為位于菱形邊長上質點的虛擬彈簧的彈力。
這里假設受力點為n0，當達到平衡狀態時質點n0運動到n0′。任取形變平面的一層進行受力分析，其它層的受力和變形情況與此類似，受力分析如圖 3所示。
圖3
形變平面的任一層截面及受力分析圖
Figure3.
Schematic diagram of force analysis in deformed plane
當受力分析點位于菱形的頂點時，其受力分析如下：
切線方向:
| $ {f_{{T_i}}}\cos {\theta _i} = {f_{{T_{i + 1}}}}\cos {\theta _{i + 1}}, $ |
法線方向:
| $ {f_{{N_i}}} = 2{f_{{T_i}}}\sin {\theta _i} - 2{f_{{T_{i + 1}}}}sin{\theta _{i + 1}}, $ |
當受力分析點位于菱形的邊長上時，其受力分析如下：
切線方向:
| $ 3{{f'}_T}_{_i}\cos {\theta _i} = {{f'}_{{T_{i + 1}}}}\cos {\theta _{i + 1}} $ |
法線方向:
| $ {{f'}_{{N_i}}} = 3f{'_T}_i\sin {\theta _i} - f{'_{{T_{i + 1}}}}\sin {\theta _{i + 1}}, $ |
式中fTi為第i層彈簧的彈力，fTi+1為第i+1層彈簧的彈力，θi為形變后相鄰層質點連線與水平方向的夾角。
因為各層菱形中只有位于邊長上的個別質點與受力點直接連接，在受到外力作用時，可以看作是同一菱形層上邊長上的質點帶動頂點處的質點移動，這樣就不需要對整個模型進行網格劃分，計算只涉及柔性物體表面的局部區域內的局部質點。其中需要計算的質點分布如圖 4所示。
圖4
模型受力質點分布圖
Figure4.
Stressed particle distribution of the model
由圖 4所示，式(1)可簡化為
| $ F = 4f{'_{{N_i}}} + \sum\limits_{i = 1}^{i = m} {{f_{{N_i}}}} $ |
由圖 4所示的受力質點可以看出，每一層菱形上的受力質點均可以看作是位于一系列同心圓上的質點，任一層區域面積的改變即可看作是同心圓之間環形面積的改變。形變達到平衡狀態時，任一層的形變面積可用下式計算,即
| $ \begin{array}{l} \Delta {S_i} = \\ \pi \left\{ {{{\left[{\left( {r + \Delta {r_i}} \right)\cos {\theta _i} + \left( {r + \Delta {r_{i + 1}}} \right)\cos {\theta _{i + 1}}} \right]}^2} - {{\left[{\left( {r + \Delta {r_{i + 1}}} \right)\cos {\theta _{i + 1}}} \right]}^2}} \right\} \end{array} $ |
其中r為表面彈簧的原始長度，Δri為表面彈簧的形變量。這里認為作用力和形變量是緩慢均勻的。當形變層i=m時，角度形變量很小，cosθi=1，則
| $ \Delta {S_i} = \pi {r^2}\left[{{{\left( {1 + \frac{{\Delta {r_i}}}{r}} \right)}^2} + 2\left( {1 + \frac{{\Delta {r_i}}}{r}} \right)\left( {1 + \frac{{\Delta {r_{i + 1}}}}{r}} \right)} \right] $ |
令，它表示第i層表面彈簧的形變率，則
| $ \Delta {S_i} = \pi {r^2}{\left( {1 + {\varepsilon _i}} \right)^2}\left[{{{\left( {1 + {\varepsilon _i}} \right)}^2} + 2\left( {1 + {\varepsilon _{i + 1}}} \right)} \right] $ |
假設任一層的虛擬彈力與形變的環形面積滿足一定的比例關系，即
| $ f{'_{{N_i}}} = k\Delta {S_i} $ |
則對于接觸力的計算，只要已知各層的形變率，就可以求出其形變區域的環形面積，同時該層的虛擬彈力也可以求出，從而算出物體表面的接觸力，這種算法大大降低了形變計算量的復雜度及次數，并保證了形變的實時性及穩定性。
2.2 模型參數的選擇
由式(10)可知，各層形變量的大小和比例系數k是影響接觸力和形變計算的關鍵，這里假設表面各層彈簧的彈性系數關系如式(11)所示：
| $ \frac{{{k_{{T_i}}}}}{{{k_{{T_1}}}}} = {e^{\beta \left( {i - 1} \right)}} $ |
則由式(4)可得
| $ 3{k_{{T_i}}}\Delta {r_i}\frac{r}{{r + \Delta {r_i}}} = {k_{{T_{i + 1}}}}\Delta {r_{i + 1}}\frac{r}{{r + \Delta {r_{i + 1}}}} $ |
| $ \frac{{\Delta {r_{i + 1}}}}{{\Delta {r_i}}} = 3\frac{{{k_{{T_i}}}}}{{{k_{{T_{i + 1}}}}}}\frac{{r + \Delta {r_{i + 1}}}}{{r + \Delta {r_i}}} $ |
| $ \frac{{\Delta {r_i}}}{{\Delta {r_1}}} = {3^{i - 1}}\frac{{{k_{{T_{i + 1}}}}}}{{{k_{{T_i}}}}}\frac{{r + \Delta {r_i}}}{{r + \Delta {r_1}}} $ |
| $ \frac{{{\varepsilon _i}}}{{{\varepsilon _1}}} = {\left( {\frac{3}{{{e^\beta }}}} \right)^{i - 1}}\frac{{1 + {\varepsilon _i}}}{{1 + {\varepsilon _1}}} $ |
為計算方便，令，則
| $ {\varepsilon _i} = \frac{{{\varepsilon _1}}}{{{a_i} + \left( {{a_i} - 1} \right){\varepsilon _1}}} $ |
以上公式說明，軟組織接觸力的計算與形變的計算是密切相關的，在整個形變計算過程中主要的計算集中于式(9)、(10)。
經典的彈簧-質點模型在進行受力分析時，系統中所有的質點在任一時刻都遵循牛頓第二定律，動力學方程為
| $ {m_i}\frac{{{\partial ^2}{x_i}}}{{\partial {t^2}}} + f_{{\rm{in}}}^i = f_{{\rm{ext}}}^i, $ |
其中mi、xi分別為質點i的質量和位置矢量，fiin為與質點i相連的所有彈簧施加的內力，fexti為質點i所受外力(摩擦力、重力、用戶施加的力等)的總和。由公式可以看出在形變計算過程中，計算所有發生移動的質點能保證反饋力的計算精度，但是計算量很大，計算速度比較慢，而在本文中，將形變區域的改變與作用力聯系起來，用虛擬彈簧產生的彈性力來計算反饋力，在一定程度上降低了計算量的復雜度，節省了計算時間，并保證了形變的真實效果。
3 仿真實驗及數據分析
3.1 虛擬手術實時仿真系統流程
建立虛擬軟組織的實時力反饋交互系統[12]，首先是形變軟組織的模型構建。從醫學圖像中提取有效數據[13]，使用3DS MAX軟件與OpenGL圖形繪制功能相結合，采用彈簧-質點模型，能夠真實模擬軟組織的形變過程。在虛擬環境中，分別建立手術器械和軟組織的物理模型，在碰撞檢測模塊[14]，虛擬手術器械與軟組織的交互定位由力反饋設備的觸覺筆尖控制。當檢測到碰撞時，根據形變區域的面積可以計算出反饋力的大小，同時操作者也能感受到反饋力的作用，形變模塊根據發生的形變量實時更新接觸點的位置和形變區域，其仿真實現流程圖如圖 5所示。
圖5
軟組織形變仿真流程圖
Figure5.
low chart of soft tissue deformation simulation
3.2 系統仿真結果
使用SensAble科技公司的PHANTOM觸覺交互設備，基于Inter(R) Pentium(R) 4.32 GHz雙核處理器，0.99 GB內存，Inter(R) 82945 G顯卡的臺式電腦，在windows XP操作系統下，采用Visual C++ 6.0開發環境并結合OpenGL三維圖形標準，搭建了虛擬軟組織實時力反饋交互系統。基于該實驗系統，實現了改進的菱形拓撲結構模型的構建，并采用一個由3 470個頂點、6 080個面組成的球體表面模型進行了實驗，實現了虛擬軟組織網格在虛擬手術器械單點和多點交互下的變形仿真，實驗結果如圖 6所示。
圖6
實驗結果
Figure6.
Experimental results
由式(6)、(9)和(10)可知模型的形變區域及接觸力的計算復雜度與形變率εi和比例系數k的選取有關，在取值時要綜合考慮精度和計算量兩方面的因素。形變深度是控制形變區域的一個參數，在實驗中取m=15，選取不同的形變率和k值，組織的形變不同，得到的反饋力也不同。本文的實驗基于matlab軟件進行仿真，分別顯示了形變率與比例系數的不同對反饋力的計算產生的影響。圖 7選擇比例系數為定值，比較了形變率對反饋力計算的影響；圖 8選擇形變率為定值，比較了比例系數對反饋力計算的影響；圖 9則比較了不同的形變率和比例系數對反饋力計算的影響。
圖7
比例系數一定，整體形變及反饋力仿真圖
Figure7.
Simulation diagram of integral deformation and force feedback with spring coefficient as a certain value
圖8
形變率一定，整體形變及反饋力仿真圖
Figure8.
Simulation diagram of integral deformation and force feedback with deformation rate as a certain value
圖9
比例系數、形變率均不同時得到的整體形變與反饋力仿真圖
Figure9.
Simulation diagram of integral deformation and force feedback with spring coefficient and deformation rate as different values
由以上三個實驗結果可知，物體與外力直接接觸的局部區域形變明顯，隨著力向周圍傳遞并逐漸減弱，形變程度也逐漸減弱，形變區域不會無限擴大，當形變邊界點的形變率達到一定值時，認為在此邊界以外外力產生的形變量為零，它的形變主要是由物體的膨脹內力產生的。模型的形變區域及接觸力的計算復雜度與形變率εi和比例常數k的選取有關，其中形變率主要影響形變區域的大小及反饋力的計算，形變率越小，計算精度越高，反之形變率越大，形變區域越大，反饋力計算量也就越大；而比例系數主要影響反饋力的計算，比例系數越大，反饋力越大，計算也就越復雜，因此在取值時要綜合考慮精度和計算量兩方面的因素。
為了驗證本文算法的可行性，采用matlab仿真軟件將本文算法與吳涓等[7]提出的同心圓分割算法進行對比，其仿真圖如圖 10所示。
圖10
相同條件下，不同形變算法的對比
Figure10.
Comparison of different deformation algorithms under the same condition
由仿真結果可以看出本文的算法更符合形變的模擬，在一定深度范圍內外力是隨著形變深度穩定增加，超過這個形變深度后力基本保持不變且趨于穩定，而同心圓的算法力是隨著形變深度遞增的，而且增長的速度越來越快，不符合形變的真實性，由此可知本文的算法具有真實性和可行性。
4 結論
針對虛擬手術中的軟組織形變仿真，提出了一種簡單可靠的彈簧-質點優化模型及形變算法。模型通過篩選有效的移動質點，計算各層形變區域的面積，計算量小，從而簡化了傳統的彈簧-質點模型的力反饋算法，形變仿真的實時性得以提高，虛擬彈簧的加入提高了形變效果的逼真度。通過對球體的形變仿真實驗，結果達到了滿意的效果，實驗結果同時表明形變率的取值對形變區域及反饋力的計算都有影響，而比例系數的取值則主要影響反饋力的計算。總之，在進行人體軟組織形變模擬時應在計算精度與計算實時性之間進行折中考慮。
引言
虛擬手術仿真系統是虛擬現實技術在醫學領域的一種應用[1]，在手術前對手術過程涉及的人體軟組織進行虛擬手術仿真具有重大的學術、經濟和社會意義，采用快速準確的形變模型和形變算法是整個手術仿真系統的關鍵，也是當前虛擬現實技術研究的難點和熱點。
目前虛擬建模已有多種方法[2]，其中主要以彈簧-質點模型、有限元模型[3]、邊界元模型等為主。各種建模方法都有其特有的優點和缺陷，相比之下，雖然彈簧-質點模型在大形變時同樣會出現失真等問題，但是彈簧-質點模型計算量較少，能夠很好地滿足仿真的實時性要求。而在一個虛擬手術系統中，實時性是衡量系統優劣最為關鍵的指標之一。因此，多數模型是在基于彈簧-質點模型[4-6]的基礎上進行改進的。吳涓等[7]的基于同心圓分割的彈簧-質點模型及實時力覺相應算法雖然計算速度快，滿足力反饋的實時性要求，但各點計算不一致，不能進一步擴展。張小瑞等[8]的層狀菱形鏈連接模型雖然力反饋和變形計算具有較高精度，但是菱形邊的長度在受力時無法改變，影響形變效果的真實性。混合模型[9]的提出，雖然滿足力反饋計算的實時性，但是混合模型的過渡區不易確定。因此如何建立穩定可靠，同時又計算簡便、實時性好的物理模型，仍然是虛擬現實領域亟待解決的一個問題[10]。
為了達到更好的模擬效果，在基于彈簧-質點模型的基礎上進行改進，提出了一種基于可變菱形的彈簧-質點模型，能實現用一種模型對不同的器官進行形變模擬的目的，提高形變效果的真實性，同時通過在模型中加入虛擬彈簧，用于形變過程中力反饋的計算，降低了形變計算的次數，得到一個更加逼真有效的形變效果，更好地滿足了交互系統的真實性、穩定性和實時性的要求。加入虛擬彈簧的同時在兩質點間也加入一個無質量的不可拉伸的輕繩，在物體受外力作用時起限定作用，當質點受力大于輕繩能承受的力的上限時，會出現輕繩的斷裂，兩質點脫離，從而模擬了現實中柔性物體的破裂現象。
1 軟組織形變模型原理及結構
1.1 軟組織形變模型原理
由于軟組織在受力狀態下，受力面的位置、施加力的大小和方向都在實時變化，給系統的穩定性和受力分析帶來了困難，而傳統的基于表面網格拓撲結構的模型，在模擬軟組織形變過程中穩定性較差，不能很好地表現軟組織的生物力學特性[11]。針對該問題，在基于層狀菱形鏈連接模型的基礎上，提出了一種改進的基于可變菱形拓撲結構的彈簧-質點模型。
由于人體軟組織器官的多樣性和復雜性，軟組織器官通常表現為不均勻性、各向異性、準不可壓縮性等材料性質，雖然目前已有很多算法模型并在實際應用中取得了良好的效果，但是在穩定性和可靠性方面還存在不少問題，建立完全真實的模型很難實現。基于可變菱形的拓撲結構，針對不同的組織器官及其本身特性，考慮到軟組織的不均勻性，例如肺癌早期，肺表面會出現比正常肺組織更硬的硬塊，在相同外力作用下肺硬塊比正常肺組織對應的形變量要小一些，因此可以設置模型中的參數，如肺硬塊菱形邊的彈簧長度長一些、彈簧系數大一點、彈簧間的夾角大一些，那么在相同作用力下，產生的形變量就會小一些，以更好地表現形變的細節，利于辨別不同的組織，為手術的成功提供一定的保障。
1.2 形變模型結構
設有虛擬彈簧的菱形結構的模型如圖 1所示。該模型由多層菱形結構單元組成，每個菱形結構單元由4個質點組成，各質點之間通過彈簧相連，施加在軟組織表面上的力，通過與其直接交互的結構單元傳遞到周圍相連的各個結構單元，由于結構單元各質點之間由彈簧連接，若彈簧一端質點發生移動，另一端質點受到彈簧力的作用位置也隨著改變。
圖1
設有虛擬彈簧的菱形結構模型
Figure1.
Diamond structure model with the virtual spring
在對整體模型進行受力分析時，物體的表面形變就是各質點的彈簧形變的疊加，而物體表面的接觸力就等效為各質點對應的虛擬彈簧彈力的合力，因此在進行受力分析時不需要分析整個虛擬軟組織，只涉及發生形變的的局部區域內的局部質點。在虛擬器械與受力柔性體接觸過程中，按照質點的連接順序劃分層次，以受力網格作為第一層，對柔性體表面每一個質點以接觸點為中心按先內后外的順序分層，這樣第一層就有4個質點，4根彈簧。每增加一層質點數增加8個，彈簧數增加16個，以此類推，第m層質點的個數為4+8(m-1)，彈簧的個數為4+16(m-1)，形成類似m層同心環狀的質點分布結構。
為了驗證本文模型拓撲結構的優越性，在相同質點數目下，將本文模型的拓撲結構與常用模型的拓撲結構進行比較，常用模型的拓撲結構如圖 2所示。多組相同質點數目下彈簧數量的比較結果如表 1所示。
圖2
常用模型的拓撲結構
Figure2.
Topology of the commonly used model
表1
彈簧數量對比
Table1.
Comparison of spring number
從表 1中數據可知，改進菱形結構分別與傳統三角形結構、改進正六邊形結構相比，在質點數量均為144時，彈簧數量各減少了32％、6%左右，在質點數量均為324時，彈簧數量各減少了33％、4%左右，在質點數量均為676時，彈簧數量各減少了33％、4%左右，因而建模復雜度低，實時性較好；而與層狀菱形鏈結構相比，在質點數量均為144時，彈簧數量增加了37.5％左右，在質點數量均為324時，彈簧數量增加了42%左右，在質點數量均為676時，彈簧數量增加了44%左右，因而穩定性較好，計算精度較高。因此，改進的拓撲結構在計算精度和實時性兩方面進行了折中，能較好地描述軟組織形變的物理特性。
2 模型形變的力反饋計算
2.1 模型形變的力反饋算法
假設質點均是理想質點，連接質點的彈簧也是理想彈簧，彈簧的阻尼力忽略不計，彈簧初始長度的變化量為零，彈簧彈力與形變量滿足胡克定律，在外力撤消后還能恢復到原狀，模型中的每一個質點的運動均滿足牛頓第二定律。在變形計算中，按照先整體后局部的思想進行。當外力F作用于任一質點時，物體的表面形變就是各質點的彈簧形變的疊加，而物體表面的接觸力就等效為各質點對應的虛擬彈簧彈力的合力，當形變層數i=m時，計算公式為
| $ F = 4\sum\limits_{i = 1}^{i = m} {{f_{{N_i}}}} + 8\sum\limits_{i = 1}^{i = m} {\left( {n - 1} \right)f{'_{{N_i}}}} , $ |
式中fNi為位于菱形頂點處質點的虛擬彈簧的彈力，fNi′為位于菱形邊長上質點的虛擬彈簧的彈力。
這里假設受力點為n0，當達到平衡狀態時質點n0運動到n0′。任取形變平面的一層進行受力分析，其它層的受力和變形情況與此類似，受力分析如圖 3所示。
圖3
形變平面的任一層截面及受力分析圖
Figure3.
Schematic diagram of force analysis in deformed plane
當受力分析點位于菱形的頂點時，其受力分析如下：
切線方向:
| $ {f_{{T_i}}}\cos {\theta _i} = {f_{{T_{i + 1}}}}\cos {\theta _{i + 1}}, $ |
法線方向:
| $ {f_{{N_i}}} = 2{f_{{T_i}}}\sin {\theta _i} - 2{f_{{T_{i + 1}}}}sin{\theta _{i + 1}}, $ |
當受力分析點位于菱形的邊長上時，其受力分析如下：
切線方向:
| $ 3{{f'}_T}_{_i}\cos {\theta _i} = {{f'}_{{T_{i + 1}}}}\cos {\theta _{i + 1}} $ |
法線方向:
| $ {{f'}_{{N_i}}} = 3f{'_T}_i\sin {\theta _i} - f{'_{{T_{i + 1}}}}\sin {\theta _{i + 1}}, $ |
式中fTi為第i層彈簧的彈力，fTi+1為第i+1層彈簧的彈力，θi為形變后相鄰層質點連線與水平方向的夾角。
因為各層菱形中只有位于邊長上的個別質點與受力點直接連接，在受到外力作用時，可以看作是同一菱形層上邊長上的質點帶動頂點處的質點移動，這樣就不需要對整個模型進行網格劃分，計算只涉及柔性物體表面的局部區域內的局部質點。其中需要計算的質點分布如圖 4所示。
圖4
模型受力質點分布圖
Figure4.
Stressed particle distribution of the model
由圖 4所示，式(1)可簡化為
| $ F = 4f{'_{{N_i}}} + \sum\limits_{i = 1}^{i = m} {{f_{{N_i}}}} $ |
由圖 4所示的受力質點可以看出，每一層菱形上的受力質點均可以看作是位于一系列同心圓上的質點，任一層區域面積的改變即可看作是同心圓之間環形面積的改變。形變達到平衡狀態時，任一層的形變面積可用下式計算,即
| $ \begin{array}{l} \Delta {S_i} = \\ \pi \left\{ {{{\left[{\left( {r + \Delta {r_i}} \right)\cos {\theta _i} + \left( {r + \Delta {r_{i + 1}}} \right)\cos {\theta _{i + 1}}} \right]}^2} - {{\left[{\left( {r + \Delta {r_{i + 1}}} \right)\cos {\theta _{i + 1}}} \right]}^2}} \right\} \end{array} $ |
其中r為表面彈簧的原始長度，Δri為表面彈簧的形變量。這里認為作用力和形變量是緩慢均勻的。當形變層i=m時，角度形變量很小，cosθi=1，則
| $ \Delta {S_i} = \pi {r^2}\left[{{{\left( {1 + \frac{{\Delta {r_i}}}{r}} \right)}^2} + 2\left( {1 + \frac{{\Delta {r_i}}}{r}} \right)\left( {1 + \frac{{\Delta {r_{i + 1}}}}{r}} \right)} \right] $ |
令，它表示第i層表面彈簧的形變率，則
| $ \Delta {S_i} = \pi {r^2}{\left( {1 + {\varepsilon _i}} \right)^2}\left[{{{\left( {1 + {\varepsilon _i}} \right)}^2} + 2\left( {1 + {\varepsilon _{i + 1}}} \right)} \right] $ |
假設任一層的虛擬彈力與形變的環形面積滿足一定的比例關系，即
| $ f{'_{{N_i}}} = k\Delta {S_i} $ |
則對于接觸力的計算，只要已知各層的形變率，就可以求出其形變區域的環形面積，同時該層的虛擬彈力也可以求出，從而算出物體表面的接觸力，這種算法大大降低了形變計算量的復雜度及次數，并保證了形變的實時性及穩定性。
2.2 模型參數的選擇
由式(10)可知，各層形變量的大小和比例系數k是影響接觸力和形變計算的關鍵，這里假設表面各層彈簧的彈性系數關系如式(11)所示：
| $ \frac{{{k_{{T_i}}}}}{{{k_{{T_1}}}}} = {e^{\beta \left( {i - 1} \right)}} $ |
則由式(4)可得
| $ 3{k_{{T_i}}}\Delta {r_i}\frac{r}{{r + \Delta {r_i}}} = {k_{{T_{i + 1}}}}\Delta {r_{i + 1}}\frac{r}{{r + \Delta {r_{i + 1}}}} $ |
| $ \frac{{\Delta {r_{i + 1}}}}{{\Delta {r_i}}} = 3\frac{{{k_{{T_i}}}}}{{{k_{{T_{i + 1}}}}}}\frac{{r + \Delta {r_{i + 1}}}}{{r + \Delta {r_i}}} $ |
| $ \frac{{\Delta {r_i}}}{{\Delta {r_1}}} = {3^{i - 1}}\frac{{{k_{{T_{i + 1}}}}}}{{{k_{{T_i}}}}}\frac{{r + \Delta {r_i}}}{{r + \Delta {r_1}}} $ |
| $ \frac{{{\varepsilon _i}}}{{{\varepsilon _1}}} = {\left( {\frac{3}{{{e^\beta }}}} \right)^{i - 1}}\frac{{1 + {\varepsilon _i}}}{{1 + {\varepsilon _1}}} $ |
為計算方便，令，則
| $ {\varepsilon _i} = \frac{{{\varepsilon _1}}}{{{a_i} + \left( {{a_i} - 1} \right){\varepsilon _1}}} $ |
以上公式說明，軟組織接觸力的計算與形變的計算是密切相關的，在整個形變計算過程中主要的計算集中于式(9)、(10)。
經典的彈簧-質點模型在進行受力分析時，系統中所有的質點在任一時刻都遵循牛頓第二定律，動力學方程為
| $ {m_i}\frac{{{\partial ^2}{x_i}}}{{\partial {t^2}}} + f_{{\rm{in}}}^i = f_{{\rm{ext}}}^i, $ |
其中mi、xi分別為質點i的質量和位置矢量，fiin為與質點i相連的所有彈簧施加的內力，fexti為質點i所受外力(摩擦力、重力、用戶施加的力等)的總和。由公式可以看出在形變計算過程中，計算所有發生移動的質點能保證反饋力的計算精度，但是計算量很大，計算速度比較慢，而在本文中，將形變區域的改變與作用力聯系起來，用虛擬彈簧產生的彈性力來計算反饋力，在一定程度上降低了計算量的復雜度，節省了計算時間，并保證了形變的真實效果。
3 仿真實驗及數據分析
3.1 虛擬手術實時仿真系統流程
建立虛擬軟組織的實時力反饋交互系統[12]，首先是形變軟組織的模型構建。從醫學圖像中提取有效數據[13]，使用3DS MAX軟件與OpenGL圖形繪制功能相結合，采用彈簧-質點模型，能夠真實模擬軟組織的形變過程。在虛擬環境中，分別建立手術器械和軟組織的物理模型，在碰撞檢測模塊[14]，虛擬手術器械與軟組織的交互定位由力反饋設備的觸覺筆尖控制。當檢測到碰撞時，根據形變區域的面積可以計算出反饋力的大小，同時操作者也能感受到反饋力的作用，形變模塊根據發生的形變量實時更新接觸點的位置和形變區域，其仿真實現流程圖如圖 5所示。
圖5
軟組織形變仿真流程圖
Figure5.
low chart of soft tissue deformation simulation
3.2 系統仿真結果
使用SensAble科技公司的PHANTOM觸覺交互設備，基于Inter(R) Pentium(R) 4.32 GHz雙核處理器，0.99 GB內存，Inter(R) 82945 G顯卡的臺式電腦，在windows XP操作系統下，采用Visual C++ 6.0開發環境并結合OpenGL三維圖形標準，搭建了虛擬軟組織實時力反饋交互系統。基于該實驗系統，實現了改進的菱形拓撲結構模型的構建，并采用一個由3 470個頂點、6 080個面組成的球體表面模型進行了實驗，實現了虛擬軟組織網格在虛擬手術器械單點和多點交互下的變形仿真，實驗結果如圖 6所示。
圖6
實驗結果
Figure6.
Experimental results
由式(6)、(9)和(10)可知模型的形變區域及接觸力的計算復雜度與形變率εi和比例系數k的選取有關，在取值時要綜合考慮精度和計算量兩方面的因素。形變深度是控制形變區域的一個參數，在實驗中取m=15，選取不同的形變率和k值，組織的形變不同，得到的反饋力也不同。本文的實驗基于matlab軟件進行仿真，分別顯示了形變率與比例系數的不同對反饋力的計算產生的影響。圖 7選擇比例系數為定值，比較了形變率對反饋力計算的影響；圖 8選擇形變率為定值，比較了比例系數對反饋力計算的影響；圖 9則比較了不同的形變率和比例系數對反饋力計算的影響。
圖7
比例系數一定，整體形變及反饋力仿真圖
Figure7.
Simulation diagram of integral deformation and force feedback with spring coefficient as a certain value
圖8
形變率一定，整體形變及反饋力仿真圖
Figure8.
Simulation diagram of integral deformation and force feedback with deformation rate as a certain value
圖9
比例系數、形變率均不同時得到的整體形變與反饋力仿真圖
Figure9.
Simulation diagram of integral deformation and force feedback with spring coefficient and deformation rate as different values
由以上三個實驗結果可知，物體與外力直接接觸的局部區域形變明顯，隨著力向周圍傳遞并逐漸減弱，形變程度也逐漸減弱，形變區域不會無限擴大，當形變邊界點的形變率達到一定值時，認為在此邊界以外外力產生的形變量為零，它的形變主要是由物體的膨脹內力產生的。模型的形變區域及接觸力的計算復雜度與形變率εi和比例常數k的選取有關，其中形變率主要影響形變區域的大小及反饋力的計算，形變率越小，計算精度越高，反之形變率越大，形變區域越大，反饋力計算量也就越大；而比例系數主要影響反饋力的計算，比例系數越大，反饋力越大，計算也就越復雜，因此在取值時要綜合考慮精度和計算量兩方面的因素。
為了驗證本文算法的可行性，采用matlab仿真軟件將本文算法與吳涓等[7]提出的同心圓分割算法進行對比，其仿真圖如圖 10所示。
圖10
相同條件下，不同形變算法的對比
Figure10.
Comparison of different deformation algorithms under the same condition
由仿真結果可以看出本文的算法更符合形變的模擬，在一定深度范圍內外力是隨著形變深度穩定增加，超過這個形變深度后力基本保持不變且趨于穩定，而同心圓的算法力是隨著形變深度遞增的，而且增長的速度越來越快，不符合形變的真實性，由此可知本文的算法具有真實性和可行性。
4 結論
針對虛擬手術中的軟組織形變仿真，提出了一種簡單可靠的彈簧-質點優化模型及形變算法。模型通過篩選有效的移動質點，計算各層形變區域的面積，計算量小，從而簡化了傳統的彈簧-質點模型的力反饋算法，形變仿真的實時性得以提高，虛擬彈簧的加入提高了形變效果的逼真度。通過對球體的形變仿真實驗，結果達到了滿意的效果，實驗結果同時表明形變率的取值對形變區域及反饋力的計算都有影響，而比例系數的取值則主要影響反饋力的計算。總之，在進行人體軟組織形變模擬時應在計算精度與計算實時性之間進行折中考慮。
