質點彈簧模型是虛擬手術中常用模型之一, 但由于其模型參數無明確物理意義, 參數設定上存在諸多不便。由此提出基于遺傳算法確定質點彈簧模型參數的方法。采用計算機輔助斷層掃描數據(CAT data)確定質點質量值; 彈簧彈性系數及阻尼系數通過遺傳算法計算獲得; 并利用參考形變同模擬形變間差異大小作為適應度函數求得模型參數近似最優解。實驗結果證明該方法能在較低計算成本的前提下獲得彈簧參數的近似最優解, 能夠使虛擬模型較準確地再現實際模型形變效果。
引用本文： 陳穎, 胡序一, 朱奇光. 虛擬手術中基于遺傳算法的質點彈簧模型參數確定. 生物醫學工程學雜志, 2015, 32(6): 1202-1206. doi: 10.7507/1001-5515.20150213 復制
0 引言
虛擬手術中最重要的一個問題為真實性表現問題[1]。形變模型物理表現和真實性決定了一個模型質量的優劣，而有限元模型和質點彈簧模型為現今主要使用的兩種模型[2]。
對于有限元模型，其實際為離散區域上的表示函數，很容易依據待模擬材料的特性方程獲得表示函數的參數。有限元模型能夠很好地表示各種交互類型，但涉及到模型拓撲結構改變時，有限元往往具有極大的局限性[3]。
近幾年研究人員設計了一些對拓撲結構改變適應性較強的模型，用于虛擬手術中對切割縫合等涉及模型拓撲結構改變的手術模擬操作[4]。以質點彈簧模型為例，雖然在模型準確性以及參數確定方面存在劣勢，但在拓撲結構改變時優良的適應性以及在圖形渲染方面極小的資源要求使其在虛擬手術建模中仍具有很大的應用價值。質點彈簧模型中，采用普適值的模型間由于網格結構不同，并不能保證有相似的形變結果[5]。不同于有限元模型，其他實例中參數設定經驗不一定適用于當前模型，因而針對某個模型對其進行參數確定是有必要的。
為保持模型真實性而對其參數進行確定一般有兩種方法。第一種方法試圖建立從真實材料物理性質到模型參數的轉換標準，基于楊氏模量、拉梅系數、泊松比等組織材料可測得物理量，并通過理論推導獲得模型質量值、倔強系數和阻尼系數等參數。Deussen等[6]通過采用Voronoi單元法，利用積分將值賦給每一個質點對象。Gelder[7]提出了一種通過比較相鄰兩個三角單元面積及楊氏模量來確定彈簧彈性的方法。但以上兩種方法對于相同模型、不同拓撲結構會產生不同仿真結果。Mollemans等[8]通過從核磁共振圖像中提取生成包含體積守恒的拓撲結構，并采用平均密度和楊氏模量對四面體彈簧參數進行設定，這種設定下的模型不能很好地復現預期形變效果。Cotin等[9]采用類似方法進行改進，獲得包含各向異性的標準彈簧質點模型，但同樣存在復現性較差的問題。第二種方法通過使用最小化程序，尋找與實際材料形變行為最接近的模型參數。Nürnberger等[10]采用神經網絡來實現對質點彈簧模型的參數設定，但網絡參數選定較為復雜。Bianchi等[11]和Louchet等[12]采用遺傳算法對質點彈簧模型參數設定進行了實驗，前者試圖重建網格拓撲來確定彈簧參數，并采用有限元模型形變作為參考形變，但參考形變同實際軟組織物理實驗得到的參數并沒有聯系，同時模型需要非均勻參數分布才能實現對線彈性材料的形變模擬；后者沒有對模型質量值進行設定，其最佳適應度函數的選取中沒有考慮位置信息同速度信息間的權重關系，且該方法只適用于面結構物體的形變模擬，對于本文中體組織的模擬效果不佳。
本文提出一種新方法來設定四面體網格質點彈簧模型參數，模型可以是任意形狀的拓撲結構。使用計算機輔助斷層掃描獲取醫學數據圖像，依據圖像信息獲取組織各個部分的位移及密度信息。以實際組織中若干點的實際形變為參考，對虛擬模型參數進行優化求解，利用遺傳算法獲得模型參數的近似最優解。
1 質點彈簧模型參數確定
1.1 質點彈簧模型數學表示
質點彈簧模型可以看成一系列由彈簧連接而成的點，這些點組成模型的網格，模型質量分配給各個離散點，模型的彈性、體積保持等物理特性通過連接各點間的彈簧體現。通常用四面體對這些點進行細小劃分，為了更好地表達組織實際特性，有些模型中附加了與彈簧并聯的阻尼器。一個標準質點彈簧模型中一般含有三類參數，即質點質量值m、質點間彈簧倔強系數k以及阻尼器阻尼系數μ。
對于質點彈簧模型中的兩個點i和j，坐標為ri和rj，連接兩點間彈簧的倔強系數和阻尼系數分別為kij和μij。彈簧初始長度用lij0表示，由此兩質點間相互作用力可表示為
| ${F_{ij}}=\left({{k_{ij}}\left({l_{^{ij}}^0-\left\| {{r_i}-{r_j}} \right\|} \right)-{\mu _{ij}}\frac{{{\text{d}}\left\| {{r_i}-{r_j}} \right\|}}{{{\text{d}}t}}} \right)\frac{{{r_i}-{r_j}}}{{\left\| {{r_i}-{r_j}} \right\|}}$ |
對于模型網格中某一個質點i，求出所有與其相連彈簧的作用力，求和得到該點受到的內力FiINT。
為體現模型體行為，四面體網格中質點i同時受到四面體t給予的體保持力Fit，其表達式為
| ${F_{ij}}=\left({\sum\limits_{j=1}^4 {\left\| {{r_i}-{r_b}} \right\|}-\sum\limits_{j=1}^4 {\left\| {r_j^0-r_b^0} \right\|} } \right)\frac{{{r_i}-{r_b}}}{{\left\| {{r_i}-{r_b}} \right\|}}$ |
式中rj(j從1~4)表示i質點所在四面體t的4個頂點坐標。rb表示四面體中心點坐標。將i點所在的所有四面體對其的體保持力累加，得到i點受到的所有體積力FiVOL。
如果i點體保持力、內力以及受到的外力FiEXT和i點質量已知，質點i的運動學方程可以表示為
| ${m_i}\frac{{{{\text{d}}^2}{r_i}}}{{{\text{d}}{t^2}}}{\text{=}}F_i^{{\text{INT}}} + F_i^{{\text{EXT}}} + F_i^{{\text{VOL}}}$ |
實際操作中外力FiEXT一般給定或由力反饋設備獲得，將式(1)、(2)代入式(3)中，整個方程包含m、k、μ以及r，前三個量需要通過遺傳算法確定。
1.2 質點質量值的確定
由于目標為保持虛擬手術真實性，選擇從計算機輔助斷層掃描數據(CAT data)獲得想要模擬器官的局部密度信息。為了解決醫學數據同模型網格的匹配問題，通過Amira軟件處理CAT數據[13]。Amira是由澳大利亞公司開發的一款醫學圖像處理軟件。通過軟件，可以得到CAT掃描的器官四面體結構。由CAT數據提取出組織器官的密度信息步驟在此不再贅述，具體過程參考文獻[14]。
由此獲得的密度信息為離散單個體素密度。通過三線性插值得到組織內任意點的密度信息，從而使模型網格中每個點的密度均可得到。針對模型中每一個四面體，計算其體積，根據其四個頂點所在位置密度信息求平均密度，將平均密度乘以四面體體積即可得到該四面體質量。得到四面體單元質量后，則質點i的質量mi可以表示為
| ${m_i}=\frac{1}{4}\sum\limits_{j=1}^n {m_j^t} $ |
每個質點被n個四面體單元復用，mjt為第j個包含質點i的四面體單元質量。為減少計算時間，采用式(4)將四面體質量分配至其4個頂點。依此做相關實驗處理，實驗結果呈現在2.1節中。結果表明，對于一般肝臟模型，模型網格尺寸在6 000以上時，文中提出的質量分配方法均可取得較好的實驗效果。
1.3 彈簧參數設定
通過遺傳算法獲得彈簧參數最優解。遺傳算法通常適用于解決最小化問題，對于本文，待求解模型的真實性問題轉化為模擬形變(即虛擬模型在力反饋設備操作下發生的形變)與參考形變(即實際生物組織在力作用下的形變)二者間誤差最小化的問題。需要通過遺傳算法找到一組k和μ的值，使以此為彈簧參數的模型形變同實際參考形變間具有最小的形變誤差。
無論是彈簧的k值還是μ值，二者都與質點彈簧模型局部楊氏模量有關，采用遺傳算法搜索k和μ值的最優解問題轉換為搜索局部楊氏模量的最優分配問題。k和μ同局部楊氏模量間關系為[15]
| ${k_1}=\frac{{{E_1} + {E_2}}}{2}\frac{{\sum\limits_j {{V_j}} }}{{l_i^2}}, $ |
| ${\mu _i}=2\frac{{\sqrt {{k_i}\left({{m_1} + {m_2}} \right)} }}{{{l_i}}}, $ |
其中E1和E2以及m1和m2分別為彈簧兩端質點的局部楊氏模量和質量，j為包含該段彈簧的所有四面體，Vj為四面體體積，lj為該段彈簧初始長度。
若模型拓撲結構和各個質點質量值已知，未知量僅為各質點局部楊氏模量，因此遺傳算法的目標為找到模型各質點的局部楊氏模量值。
遺傳算法中需確定的問題主要包括選擇、交叉、變異等行為的概率和操作運算[16]。在算法每一次迭代中替換掉50%的個體，對個體間采用中間重組的方式以0.9的概率進行交叉。之后對群中每一個體以變異概率p判斷是否進行變異操作。變異概率p定義如下：
| $p=\frac{1}{{n + 2}} + 0.002$ |
對于確定變異個體，按照50%概率選擇變異位進行變異。變異后個體值需滿足以變異前個體值為均值、以std(n)為偏差的正態分布。
| $std\left(n \right)=\left({{E_{\max }}-{E_{\min }}} \right)\left({\frac{1}{{n + 1}} + 0.002} \right)$ |
式(7)、(8)中的n均為當前迭代次數，Emax和Emin為事先設定局部楊氏模量的上下限，交叉概率設定考慮了隨著迭代次數上升個體間差異將下降，交叉操作對種群影響也會下降等因素。綜合考慮并結合其他類似問題的概率值，選擇0.9這個比較大的交叉概率。若交叉概率較小，將會大大增加算法收斂時間。選定交叉概率后，依據文獻[17]所述概率選擇經驗結論，設定式(7)所示的變異概率。
對個體適應度函數的確定同樣為重要步驟，適應度函數要求單值連續且具有合理一致性，故采用目標形變同當前形變間的差異評價個體優劣。遺傳算法的本質是尋找最高適應度的個體，因此采用模擬形變中關鍵點位置和速度與參考形變的差值加和后取反，作為個體適應度函數，用式(9)表示，即
| ${f_i}=-\sum\limits_{t=1}^n {\mathop \Sigma \limits^p } \left({r_p^i\left(t \right)-r_p^s\left(t \right) + \alpha \left({\frac{{{\text{d}}r_p^i\left(t \right)}}{{{\text{d}}t}}-\frac{{{\text{d}}r_p^s\left(t \right)}}{{{\text{d}}t}}} \right)} \right)$ |
其中n為虛擬模型形變步數，p為所有與參考模型比較的關鍵點，rpi(t)和rps(t)分別為時間步數t時刻p點在虛擬形變和參考形變中的坐標。α為調整位置誤差和速度誤差間的權重參數，一般根據形變類型設為0.002 5。由此通過遺傳算法獲得最佳適應度個體，分配模型中各點的局部楊氏模量。在這種利用參考形變作為遺傳算法期望模型的方法中，獲得的最優虛擬模型能夠很好地模擬與實際效果相似的虛擬形變。
2 仿真實驗研究
2.1 質點質量分配效果
為了更好地體現質點質量分配效果，采用Amira軟件對肝臟組織模型進行質量分配，不同精細程度模型網格均源于同一醫學數據(CAT data)。具體數據如表 1所示。
表1
質點質量分配表
Table1.
Mass distribution table
可以看出，隨著肝臟模型網格精細程度增加，分配后模型總體質量有小幅提高，符合高精度模型更多地占滿組織器官體積空間的實際情況。總體質量相互之間相差不大于3.5%，同時平均密度沒有太大變化，說明此方法質量分配效果良好，針對不同精細程度的模型網格均能較好適應。
2.2 彈簧參數優化效果
采用1.3節所述步驟，構造一個具有120個點、377個四面體、590個彈簧的立方體網格模型，分別對其頂視、左視、前視三個面中心施加0.001 N指向立方體中心的外力。在C++環境下采用GAlib類庫設計實現遺傳算法。整個仿真系統如圖 1所示，采用力反饋設備實現人機交互，計算虛擬形變與參考形變間誤差的大小。
圖1
仿真系統圖
Figure1.
Simulation system
算法執行中個體最大適應度與迭代次數之間的關系如圖 2所示。
圖2
個體最大適應度同迭代次數的關系
Figure2.
Relationship between largest individual fitness and iterations
由圖 2可以得出，在迭代次數達到1 000次左右時，個體最大適應度近似收斂到最優解，此時適應度達到-0.05左右，整個優化時間約為1 h，證明該方法在時間上可行。
將迭代后最優解作為質點彈簧模型參數進行運算，分析其在三個方向上與參考形變間誤差大小與模型運算步長數之間的關系，結果如圖 3所示。圖中從左至右分別表示前視、左視、頂視三個方向上受力時模擬形變與參考形變間的誤差。
圖3
虛擬形變與參考形變間的誤差
Figure3.
Deviations between virtual deformations and reference
圖 3中橫軸表示彈簧質點模型虛擬形變時所花費時間步長的個數，縱軸表示模型在各個方向上與參考形變各點間的誤差距離。在虛擬形變的整個時間范圍內，隨著步數增加誤差有所上升，但均小于5 mm。誤差增大是由于步數增加導致每一步誤差累積所致，符合式(9)中誤差公式的表述。左視面受力時其z軸誤差較其他兩軸異號，這是由于實驗所用立方體的四面體網格建模時四面體底面處于左視面上，反映了建模方式的差異同樣會影響誤差分布。實驗所用參考形變采用從頂視面施加壓力，因而模擬形變在頂視方向上誤差最小，在其余方向誤差較頂視面大。但整體誤差均處于較小范圍，模型總體的形變表現良好。
3 結論
針對質點彈簧模型在參數設定上的諸多不確定性，提出以計算機輔助斷層掃描等醫學數據為基礎的質點彈簧模型各質點質量分配方法，完成了對不同精細程度網格模型質量參數的確定，效果良好，對不同網格均具有較好的適應性。同時采用遺傳算法，以參考形變為基礎，獲得彈簧阻尼系數和倔強系數的最佳設定。優化后的質點彈簧模型在各方向均能較好地復現參考模型的形變效果，方法真實可靠，在質點彈簧模型的基礎上大大提高了其形變的真實程度。
0 引言
虛擬手術中最重要的一個問題為真實性表現問題[1]。形變模型物理表現和真實性決定了一個模型質量的優劣，而有限元模型和質點彈簧模型為現今主要使用的兩種模型[2]。
對于有限元模型，其實際為離散區域上的表示函數，很容易依據待模擬材料的特性方程獲得表示函數的參數。有限元模型能夠很好地表示各種交互類型，但涉及到模型拓撲結構改變時，有限元往往具有極大的局限性[3]。
近幾年研究人員設計了一些對拓撲結構改變適應性較強的模型，用于虛擬手術中對切割縫合等涉及模型拓撲結構改變的手術模擬操作[4]。以質點彈簧模型為例，雖然在模型準確性以及參數確定方面存在劣勢，但在拓撲結構改變時優良的適應性以及在圖形渲染方面極小的資源要求使其在虛擬手術建模中仍具有很大的應用價值。質點彈簧模型中，采用普適值的模型間由于網格結構不同，并不能保證有相似的形變結果[5]。不同于有限元模型，其他實例中參數設定經驗不一定適用于當前模型，因而針對某個模型對其進行參數確定是有必要的。
為保持模型真實性而對其參數進行確定一般有兩種方法。第一種方法試圖建立從真實材料物理性質到模型參數的轉換標準，基于楊氏模量、拉梅系數、泊松比等組織材料可測得物理量，并通過理論推導獲得模型質量值、倔強系數和阻尼系數等參數。Deussen等[6]通過采用Voronoi單元法，利用積分將值賦給每一個質點對象。Gelder[7]提出了一種通過比較相鄰兩個三角單元面積及楊氏模量來確定彈簧彈性的方法。但以上兩種方法對于相同模型、不同拓撲結構會產生不同仿真結果。Mollemans等[8]通過從核磁共振圖像中提取生成包含體積守恒的拓撲結構，并采用平均密度和楊氏模量對四面體彈簧參數進行設定，這種設定下的模型不能很好地復現預期形變效果。Cotin等[9]采用類似方法進行改進，獲得包含各向異性的標準彈簧質點模型，但同樣存在復現性較差的問題。第二種方法通過使用最小化程序，尋找與實際材料形變行為最接近的模型參數。Nürnberger等[10]采用神經網絡來實現對質點彈簧模型的參數設定，但網絡參數選定較為復雜。Bianchi等[11]和Louchet等[12]采用遺傳算法對質點彈簧模型參數設定進行了實驗，前者試圖重建網格拓撲來確定彈簧參數，并采用有限元模型形變作為參考形變，但參考形變同實際軟組織物理實驗得到的參數并沒有聯系，同時模型需要非均勻參數分布才能實現對線彈性材料的形變模擬；后者沒有對模型質量值進行設定，其最佳適應度函數的選取中沒有考慮位置信息同速度信息間的權重關系，且該方法只適用于面結構物體的形變模擬，對于本文中體組織的模擬效果不佳。
本文提出一種新方法來設定四面體網格質點彈簧模型參數，模型可以是任意形狀的拓撲結構。使用計算機輔助斷層掃描獲取醫學數據圖像，依據圖像信息獲取組織各個部分的位移及密度信息。以實際組織中若干點的實際形變為參考，對虛擬模型參數進行優化求解，利用遺傳算法獲得模型參數的近似最優解。
1 質點彈簧模型參數確定
1.1 質點彈簧模型數學表示
質點彈簧模型可以看成一系列由彈簧連接而成的點，這些點組成模型的網格，模型質量分配給各個離散點，模型的彈性、體積保持等物理特性通過連接各點間的彈簧體現。通常用四面體對這些點進行細小劃分，為了更好地表達組織實際特性，有些模型中附加了與彈簧并聯的阻尼器。一個標準質點彈簧模型中一般含有三類參數，即質點質量值m、質點間彈簧倔強系數k以及阻尼器阻尼系數μ。
對于質點彈簧模型中的兩個點i和j，坐標為ri和rj，連接兩點間彈簧的倔強系數和阻尼系數分別為kij和μij。彈簧初始長度用lij0表示，由此兩質點間相互作用力可表示為
| ${F_{ij}}=\left({{k_{ij}}\left({l_{^{ij}}^0-\left\| {{r_i}-{r_j}} \right\|} \right)-{\mu _{ij}}\frac{{{\text{d}}\left\| {{r_i}-{r_j}} \right\|}}{{{\text{d}}t}}} \right)\frac{{{r_i}-{r_j}}}{{\left\| {{r_i}-{r_j}} \right\|}}$ |
對于模型網格中某一個質點i，求出所有與其相連彈簧的作用力，求和得到該點受到的內力FiINT。
為體現模型體行為，四面體網格中質點i同時受到四面體t給予的體保持力Fit，其表達式為
| ${F_{ij}}=\left({\sum\limits_{j=1}^4 {\left\| {{r_i}-{r_b}} \right\|}-\sum\limits_{j=1}^4 {\left\| {r_j^0-r_b^0} \right\|} } \right)\frac{{{r_i}-{r_b}}}{{\left\| {{r_i}-{r_b}} \right\|}}$ |
式中rj(j從1~4)表示i質點所在四面體t的4個頂點坐標。rb表示四面體中心點坐標。將i點所在的所有四面體對其的體保持力累加，得到i點受到的所有體積力FiVOL。
如果i點體保持力、內力以及受到的外力FiEXT和i點質量已知，質點i的運動學方程可以表示為
| ${m_i}\frac{{{{\text{d}}^2}{r_i}}}{{{\text{d}}{t^2}}}{\text{=}}F_i^{{\text{INT}}} + F_i^{{\text{EXT}}} + F_i^{{\text{VOL}}}$ |
實際操作中外力FiEXT一般給定或由力反饋設備獲得，將式(1)、(2)代入式(3)中，整個方程包含m、k、μ以及r，前三個量需要通過遺傳算法確定。
1.2 質點質量值的確定
由于目標為保持虛擬手術真實性，選擇從計算機輔助斷層掃描數據(CAT data)獲得想要模擬器官的局部密度信息。為了解決醫學數據同模型網格的匹配問題，通過Amira軟件處理CAT數據[13]。Amira是由澳大利亞公司開發的一款醫學圖像處理軟件。通過軟件，可以得到CAT掃描的器官四面體結構。由CAT數據提取出組織器官的密度信息步驟在此不再贅述，具體過程參考文獻[14]。
由此獲得的密度信息為離散單個體素密度。通過三線性插值得到組織內任意點的密度信息，從而使模型網格中每個點的密度均可得到。針對模型中每一個四面體，計算其體積，根據其四個頂點所在位置密度信息求平均密度，將平均密度乘以四面體體積即可得到該四面體質量。得到四面體單元質量后，則質點i的質量mi可以表示為
| ${m_i}=\frac{1}{4}\sum\limits_{j=1}^n {m_j^t} $ |
每個質點被n個四面體單元復用，mjt為第j個包含質點i的四面體單元質量。為減少計算時間，采用式(4)將四面體質量分配至其4個頂點。依此做相關實驗處理，實驗結果呈現在2.1節中。結果表明，對于一般肝臟模型，模型網格尺寸在6 000以上時，文中提出的質量分配方法均可取得較好的實驗效果。
1.3 彈簧參數設定
通過遺傳算法獲得彈簧參數最優解。遺傳算法通常適用于解決最小化問題，對于本文，待求解模型的真實性問題轉化為模擬形變(即虛擬模型在力反饋設備操作下發生的形變)與參考形變(即實際生物組織在力作用下的形變)二者間誤差最小化的問題。需要通過遺傳算法找到一組k和μ的值，使以此為彈簧參數的模型形變同實際參考形變間具有最小的形變誤差。
無論是彈簧的k值還是μ值，二者都與質點彈簧模型局部楊氏模量有關，采用遺傳算法搜索k和μ值的最優解問題轉換為搜索局部楊氏模量的最優分配問題。k和μ同局部楊氏模量間關系為[15]
| ${k_1}=\frac{{{E_1} + {E_2}}}{2}\frac{{\sum\limits_j {{V_j}} }}{{l_i^2}}, $ |
| ${\mu _i}=2\frac{{\sqrt {{k_i}\left({{m_1} + {m_2}} \right)} }}{{{l_i}}}, $ |
其中E1和E2以及m1和m2分別為彈簧兩端質點的局部楊氏模量和質量，j為包含該段彈簧的所有四面體，Vj為四面體體積，lj為該段彈簧初始長度。
若模型拓撲結構和各個質點質量值已知，未知量僅為各質點局部楊氏模量，因此遺傳算法的目標為找到模型各質點的局部楊氏模量值。
遺傳算法中需確定的問題主要包括選擇、交叉、變異等行為的概率和操作運算[16]。在算法每一次迭代中替換掉50%的個體，對個體間采用中間重組的方式以0.9的概率進行交叉。之后對群中每一個體以變異概率p判斷是否進行變異操作。變異概率p定義如下：
| $p=\frac{1}{{n + 2}} + 0.002$ |
對于確定變異個體，按照50%概率選擇變異位進行變異。變異后個體值需滿足以變異前個體值為均值、以std(n)為偏差的正態分布。
| $std\left(n \right)=\left({{E_{\max }}-{E_{\min }}} \right)\left({\frac{1}{{n + 1}} + 0.002} \right)$ |
式(7)、(8)中的n均為當前迭代次數，Emax和Emin為事先設定局部楊氏模量的上下限，交叉概率設定考慮了隨著迭代次數上升個體間差異將下降，交叉操作對種群影響也會下降等因素。綜合考慮并結合其他類似問題的概率值，選擇0.9這個比較大的交叉概率。若交叉概率較小，將會大大增加算法收斂時間。選定交叉概率后，依據文獻[17]所述概率選擇經驗結論，設定式(7)所示的變異概率。
對個體適應度函數的確定同樣為重要步驟，適應度函數要求單值連續且具有合理一致性，故采用目標形變同當前形變間的差異評價個體優劣。遺傳算法的本質是尋找最高適應度的個體，因此采用模擬形變中關鍵點位置和速度與參考形變的差值加和后取反，作為個體適應度函數，用式(9)表示，即
| ${f_i}=-\sum\limits_{t=1}^n {\mathop \Sigma \limits^p } \left({r_p^i\left(t \right)-r_p^s\left(t \right) + \alpha \left({\frac{{{\text{d}}r_p^i\left(t \right)}}{{{\text{d}}t}}-\frac{{{\text{d}}r_p^s\left(t \right)}}{{{\text{d}}t}}} \right)} \right)$ |
其中n為虛擬模型形變步數，p為所有與參考模型比較的關鍵點，rpi(t)和rps(t)分別為時間步數t時刻p點在虛擬形變和參考形變中的坐標。α為調整位置誤差和速度誤差間的權重參數，一般根據形變類型設為0.002 5。由此通過遺傳算法獲得最佳適應度個體，分配模型中各點的局部楊氏模量。在這種利用參考形變作為遺傳算法期望模型的方法中，獲得的最優虛擬模型能夠很好地模擬與實際效果相似的虛擬形變。
2 仿真實驗研究
2.1 質點質量分配效果
為了更好地體現質點質量分配效果，采用Amira軟件對肝臟組織模型進行質量分配，不同精細程度模型網格均源于同一醫學數據(CAT data)。具體數據如表 1所示。
表1
質點質量分配表
Table1.
Mass distribution table
可以看出，隨著肝臟模型網格精細程度增加，分配后模型總體質量有小幅提高，符合高精度模型更多地占滿組織器官體積空間的實際情況。總體質量相互之間相差不大于3.5%，同時平均密度沒有太大變化，說明此方法質量分配效果良好，針對不同精細程度的模型網格均能較好適應。
2.2 彈簧參數優化效果
采用1.3節所述步驟，構造一個具有120個點、377個四面體、590個彈簧的立方體網格模型，分別對其頂視、左視、前視三個面中心施加0.001 N指向立方體中心的外力。在C++環境下采用GAlib類庫設計實現遺傳算法。整個仿真系統如圖 1所示，采用力反饋設備實現人機交互，計算虛擬形變與參考形變間誤差的大小。
圖1
仿真系統圖
Figure1.
Simulation system
算法執行中個體最大適應度與迭代次數之間的關系如圖 2所示。
圖2
個體最大適應度同迭代次數的關系
Figure2.
Relationship between largest individual fitness and iterations
由圖 2可以得出，在迭代次數達到1 000次左右時，個體最大適應度近似收斂到最優解，此時適應度達到-0.05左右，整個優化時間約為1 h，證明該方法在時間上可行。
將迭代后最優解作為質點彈簧模型參數進行運算，分析其在三個方向上與參考形變間誤差大小與模型運算步長數之間的關系，結果如圖 3所示。圖中從左至右分別表示前視、左視、頂視三個方向上受力時模擬形變與參考形變間的誤差。
圖3
虛擬形變與參考形變間的誤差
Figure3.
Deviations between virtual deformations and reference
圖 3中橫軸表示彈簧質點模型虛擬形變時所花費時間步長的個數，縱軸表示模型在各個方向上與參考形變各點間的誤差距離。在虛擬形變的整個時間范圍內，隨著步數增加誤差有所上升，但均小于5 mm。誤差增大是由于步數增加導致每一步誤差累積所致，符合式(9)中誤差公式的表述。左視面受力時其z軸誤差較其他兩軸異號，這是由于實驗所用立方體的四面體網格建模時四面體底面處于左視面上，反映了建模方式的差異同樣會影響誤差分布。實驗所用參考形變采用從頂視面施加壓力，因而模擬形變在頂視方向上誤差最小，在其余方向誤差較頂視面大。但整體誤差均處于較小范圍，模型總體的形變表現良好。
3 結論
針對質點彈簧模型在參數設定上的諸多不確定性，提出以計算機輔助斷層掃描等醫學數據為基礎的質點彈簧模型各質點質量分配方法，完成了對不同精細程度網格模型質量參數的確定，效果良好，對不同網格均具有較好的適應性。同時采用遺傳算法，以參考形變為基礎，獲得彈簧阻尼系數和倔強系數的最佳設定。優化后的質點彈簧模型在各方向均能較好地復現參考模型的形變效果，方法真實可靠，在質點彈簧模型的基礎上大大提高了其形變的真實程度。
