交互式血管擠壓形變仿真方法研究①

高志堅, 陳國棟, 王 娜

(1.福州大學(xué)物理與信息工程學(xué)院，福建 福州 350100 ;2.福建師范大學(xué)福清分校電子與信息工程學(xué)院, 福建 福州 350300 3.無損檢測技術(shù)福建省高等學(xué)校重點(diǎn)實驗室(福建師范大學(xué)福清分校), 福建 福州 350300)

0 引 言

血管形變是心血管類虛擬手術(shù)的重要組成部分?，F(xiàn)有人體組織形變模型可分為彈簧質(zhì)點(diǎn)法、有限元法和無網(wǎng)格法[1]。其中，陳寒青等人[2]借助有限元線性模型建立軟組織的物理模型，使用基于靜力凝聚法的子空間凝聚法進(jìn)行軟組織形變計算，并利用子空間并行迭代計算提高了仿真的實時性，但當(dāng)全空間半徑較小時，仿真的結(jié)果失真變形。

真實感和實時性是虛擬手術(shù)中難以同時兼顧的難點(diǎn)?？紤]到心血管類疾病操作精密，其仿真更應(yīng)重視真實感。為此，基于彈性力學(xué)的形變能夠較為逼真的重現(xiàn)手術(shù)對象。

1 形變模型

形變模型包括體現(xiàn)物體外觀特征的幾何模型和描述形變的物理模型。其中，幾何模型是物體形變的前提條件[3]，通常借助拓?fù)浣Y(jié)構(gòu)構(gòu)建。物理模型由彈性力學(xué)理論推導(dǎo)得出。此外，簡化血管的受力情況可以降低表達(dá)式復(fù)雜程度并提高仿真效率。

1.1 幾何模型

從視覺效果和簡易性出發(fā)，血管的外觀可用細(xì)長管道進(jìn)行模擬[4]。圖1為來自課題組人體可視化計劃VHP重建的幾何模型，其包括2940個頂點(diǎn)和735個三角面，如圖1所示。

圖1 血管的幾何模型

1.2 物理模型

血管受力后發(fā)生的位移由無限三維連續(xù)體動力學(xué)[5]的時變位移場確定，其表達(dá)式如公式(1)：

(1)

式中u為時變位移場，t為時間，b為外力。參數(shù)μ和ν分別為彈性剪切模量和泊松比。根據(jù)文獻(xiàn)[6]對血管的生物力學(xué)參數(shù)的測定，后續(xù)的實驗將參數(shù)μ設(shè)置為14.24kPa，參數(shù)ν設(shè)置為0.45。

1.3 形變力

在仿真中，通常用矢量代表力。力矢量作用于物體某個節(jié)點(diǎn)，進(jìn)而導(dǎo)致整個物體位置和形狀的改變。然而，單一載荷容易導(dǎo)致動力學(xué)方程求解出現(xiàn)病態(tài)問題。考慮血管壁和血管鉗的幾何形狀，使用標(biāo)準(zhǔn)正態(tài)分布函數(shù)將外力進(jìn)行分布，以此替換單一荷載見公式(2)：

(2)

等式左側(cè)為與受力中心C相距d處的形變壓力，為矢量。等式右側(cè)的θ為標(biāo)量，用于調(diào)控外力f的大小，且限制為2≤θ≤8。δ(t)為階躍函數(shù)，表示形變壓力施加的時刻。

2 碰撞檢測

碰撞檢測用于解決手術(shù)對象和手術(shù)器械之間是否碰撞、何時碰撞、何處碰撞等三大問題，現(xiàn)有主流算法可大致分為空間網(wǎng)格劃分法和包圍盒法[7]。其中，空間網(wǎng)格劃分法將三維空間分為若干單元網(wǎng)格，若含有不同物體的子網(wǎng)格產(chǎn)生交互，則認(rèn)為發(fā)生碰撞。包圍盒法是利用簡單結(jié)構(gòu)的包圍體對三維物體進(jìn)行擬合。若兩個包圍體產(chǎn)生交互，即認(rèn)為物體間已經(jīng)發(fā)生碰撞。包圍盒法相較于空間網(wǎng)格劃分法較為高效但有一定出錯的概率。

2.1 碰撞檢測的兩大階段

血管形變仿真中的碰撞檢測可以分為如圖2所示的兩個階段：

(1)檢測血管鉗是否與血管外壁發(fā)生碰撞。若發(fā)生，則獲取碰撞點(diǎn)的三維坐標(biāo)和時間t，并對外力f進(jìn)行分布。

(2)在血管形變過程中檢測其內(nèi)壁是否發(fā)生碰撞。若發(fā)生，則表明形變已經(jīng)到達(dá)極限，必須固定形變力。

圖2 夾持過程中不同階段鉗尖與血管外壁位置

2.2 碰撞檢測算法

第一階段中，使用軸對齊包圍盒(AABB)和包圍球(Sphere)分別擬合血管壁和血管鉗的鉗尖部位。即AABB-Sphere碰撞檢測問題(圖3)。注意，軸對齊包圍盒對圓柱體的擬合必然存在空隙，因此實驗中必須將鉗尖與血管壁垂直。

AABB-Sphere碰撞檢測可細(xì)分為如下3個步驟：

(1)三維空間下舍棄其中一維度的坐標(biāo)，將三維對象轉(zhuǎn)化為二維圖形。如此，盒子退化為矩形，并用R表示；球體退化為圓形，并用C表示。

圖3 AABB-Sphere碰撞檢測

(3)從三個維度分別重復(fù)步驟1和步驟2，并統(tǒng)計每一維度下的C與R間碰撞檢測結(jié)果。設(shè)六次檢測結(jié)果分別為O11,O12,O21,O22,O31,O32，其中第一個下標(biāo)代表舍棄的第i維度，第二個下標(biāo)代表二維下的第j維度。若上述六個結(jié)果在某一時刻全部為真，才認(rèn)為包圍球和包圍盒產(chǎn)生了交互。

碰撞檢測的第二階段基于第一階段。該階段依靠擬合鉗尖的兩個球體進(jìn)行碰撞檢測，即Sphere-Sphere碰撞檢測問題(圖4)。其中，球的半徑r與血管壁厚度一致。

圖4 Sphere-Sphere碰撞檢測

Sphere-Sphere碰撞檢測可細(xì)分為如下3個步驟：

(1)讀取兩個包圍球的球心O1及O2和半徑r1和r2。

(2)計算兩個球心之間相對距離的絕對值，即s=|O1-O2|。

(3)判斷距離s與半徑r1和r2之間的關(guān)系。若有s>r1+r2，則兩個球體間未發(fā)生碰撞。否則發(fā)生碰撞。

3 受力-形變算法

在圖形學(xué)中，物體的幾何模型由若干以特定規(guī)則相連的頂點(diǎn)構(gòu)成。任一頂點(diǎn)pi都有指定的幾何坐標(biāo)和法向量。根據(jù)各頂點(diǎn)的位置與坐標(biāo)軸原點(diǎn)O之間的差可以構(gòu)造頂點(diǎn)位置矢量場p，即：

pi=pi-O

(3)

頂點(diǎn)位置矢量場p可以視為三維空間中的物體在靜止?fàn)顟B(tài)下一組離散的點(diǎn)。當(dāng)作用力施加到物體上后，原先位于空間中某個位置的點(diǎn)移動到了新的位置，所有頂點(diǎn)的移動產(chǎn)生了物體位置和形狀的改變。而點(diǎn)的位移則使用位移場物u進(jìn)行描述。

3.1 位移場表達(dá)式推導(dǎo)

根據(jù)離散信號處理理論[8]和彈性動力學(xué)[5]理論，將式(1)分解為單位脈沖函數(shù)δ(t)與單位階躍函數(shù)H(t)的線性組合，即：

(4)

圖5 算法流程圖

圖6 血管在不同大小的力度下形變的過程(側(cè)面)

等式左側(cè)為與位置場p及時間t有關(guān)的時變位移場。根據(jù)點(diǎn)源和場的疊加原理[9]，借助格林求解法對式(4)中的位移場進(jìn)行求解，有u(pi,t)=D(pi,t)f，其中D(pi,t)即為彈性動力學(xué)的格林函數(shù)。

由標(biāo)準(zhǔn)正態(tài)分布函數(shù)的對稱性，借助Helmholtz[10]分解法可以將式(2)簡化為如下形式：

(5)

圖7血管在不同大小的力度下形變的過程(正面)

聯(lián)合式(4)和式(5)，可得到動力學(xué)(1)對應(yīng)于分布力(2)的解為：

(6)

其中

公式(6)在位移場和形變外力f之間建立了聯(lián)系，根據(jù)該式，空間中所有頂點(diǎn)可通過外力f及其所處位置得到相應(yīng)的位移矢量。

3.2 形變計算

當(dāng)血管鉗與血管發(fā)生碰撞后，借助血管的空間頂點(diǎn)位置矢量場pi和位移場u可以計算血管在受力后的新坐標(biāo)。即：

pnew=pold+u(pi,t)

(7)

該式表明血管的任一頂點(diǎn)在形變后的空間位置由該點(diǎn)的位置矢量和位移矢量共同決定。通過求解格林函數(shù)D(pi,t)，并將以受力載荷為中心的所有頂點(diǎn)的位置矢量進(jìn)行迭代計算，即可得到血管在受力之后發(fā)生的連續(xù)形變位移效果。

3.3 程序流程

血管形變仿真流程如圖5所示。

4 實驗結(jié)果與分析

4.1 實驗條件與過程

實驗在64位Windows7系統(tǒng)下進(jìn)行。使用C#編程語言和64位2018.3版本的Unity3D圖形引擎。硬件為Intel i7-4700 @2. 40GHz，NVIDIA GTX860M和16G內(nèi)存。

血管形變仿真的工程實現(xiàn)可分為如下六個步驟：

(1)讀取OBJ文件，使用結(jié)構(gòu)體存儲三維物體的幾何信息。

(2)根據(jù)血管的幾何坐標(biāo)和原點(diǎn)坐標(biāo)計算各頂點(diǎn)的位置矢量。

(3)監(jiān)聽鼠標(biāo)，交互式控制血管鉗夾持的位置和夾持力度的大小。

(4)借助Eigen矩陣庫求解位移場表達(dá)式。

(5)疊加各頂點(diǎn)的位置矢量場和位移場。

(6)重復(fù)步驟(4)和步驟(5)，直到位移場為零，結(jié)束迭代。

4.2 實驗結(jié)果展示與分析

圖6從側(cè)面角度展示了血管在血管鉗夾持下發(fā)生的擠壓性形變。圖6a為夾持之前，血管鉗沒有與血管接觸；圖6b和圖6c為血管鉗夾持過程中血管的形變狀態(tài)，可以看出血管模型發(fā)生了具有較強(qiáng)真實感的形變；圖6d為血管上下壁完全閉合時的狀態(tài)，此時無法繼續(xù)增大夾持的力度，符合物理直覺且具有良好的視覺效果。

圖7從正面角度展示了血管壁在擠壓的過程中的形變。從圖7b和圖7c可以看出血管壁的厚度和體積在形變的過程中始終得到保持，并未發(fā)生塌陷；從圖7d可以看出血管的上下壁緊密貼合，并未發(fā)生穿透。

仿真中的畫面每秒傳輸幀數(shù)(FPS)是用于衡量仿真算法實時性的重要指標(biāo)。表1展示了血管模型在不同數(shù)目的頂點(diǎn)下的每秒傳輸幀數(shù)(平均值)?？梢婋S著血管模型頂點(diǎn)數(shù)目增加，幀率顯著降低，實時性下降。因此在實際應(yīng)用的過程中，要充分衡量真實感與實時性二者之間的重要程度。

表1 仿真中不同頂點(diǎn)數(shù)對應(yīng)的FPS

5 結(jié) 論

為了得到具有較強(qiáng)真實感的交互式仿真效果。使用了動力學(xué)方程和正態(tài)分布的形變力進(jìn)行計算。由于血管和血管鉗鉗尖特殊的幾何形狀，使用包圍盒法進(jìn)行碰撞檢測。根據(jù)位移場求解頂點(diǎn)在受力后的新坐標(biāo)以實現(xiàn)形變效果。實驗結(jié)果表明，算法基本能夠滿足系統(tǒng)對血管交互式形變仿真的真實性和實時性要求。但對于其他手術(shù)器械對血管的操作仿真將在今后的工作中進(jìn)一步探討。