交互式血管挤压形变仿真方法研究①

高志坚, 陈国栋, 王 娜

(1.福州大学物理与信息工程学院，福建 福州 350100 ;2.福建师范大学福清分校电子与信息工程学院, 福建 福州 350300 3.无损检测技术福建省高等学校重点实验室(福建师范大学福清分校), 福建 福州 350300)

0 引 言

血管形变是心血管类虚拟手术的重要组成部分。现有人体组织形变模型可分为弹簧质点法、有限元法和无网格法[1]。其中，陈寒青等人[2]借助有限元线性模型建立软组织的物理模型，使用基于静力凝聚法的子空间凝聚法进行软组织形变计算，并利用子空间并行迭代计算提高了仿真的实时性，但当全空间半径较小时，仿真的结果失真变形。

真实感和实时性是虚拟手术中难以同时兼顾的难点。考虑到心血管类疾病操作精密，其仿真更应重视真实感。为此，基于弹性力学的形变能够较为逼真的重现手术对象。

1 形变模型

形变模型包括体现物体外观特征的几何模型和描述形变的物理模型。其中，几何模型是物体形变的前提条件[3]，通常借助拓扑结构构建。物理模型由弹性力学理论推导得出。此外，简化血管的受力情况可以降低表达式复杂程度并提高仿真效率。

1.1 几何模型

从视觉效果和简易性出发，血管的外观可用细长管道进行模拟[4]。图1为来自课题组人体可视化计划VHP重建的几何模型，其包括2940个顶点和735个三角面，如图1所示。

图1 血管的几何模型

1.2 物理模型

血管受力后发生的位移由无限三维连续体动力学[5]的时变位移场确定，其表达式如公式(1)：

(1)

式中u为时变位移场，t为时间，b为外力。参数μ和ν分别为弹性剪切模量和泊松比。根据文献[6]对血管的生物力学参数的测定，后续的实验将参数μ设置为14.24kPa，参数ν设置为0.45。

1.3 形变力

在仿真中，通常用矢量代表力。力矢量作用于物体某个节点，进而导致整个物体位置和形状的改变。然而，单一载荷容易导致动力学方程求解出现病态问题。考虑血管壁和血管钳的几何形状，使用标准正态分布函数将外力进行分布，以此替换单一荷载见公式(2)：

(2)

等式左侧为与受力中心C相距d处的形变压力，为矢量。等式右侧的θ为标量，用于调控外力f的大小，且限制为2≤θ≤8。δ(t)为阶跃函数，表示形变压力施加的时刻。

2 碰撞检测

碰撞检测用于解决手术对象和手术器械之间是否碰撞、何时碰撞、何处碰撞等三大问题，现有主流算法可大致分为空间网格划分法和包围盒法[7]。其中，空间网格划分法将三维空间分为若干单元网格，若含有不同物体的子网格产生交互，则认为发生碰撞。包围盒法是利用简单结构的包围体对三维物体进行拟合。若两个包围体产生交互，即认为物体间已经发生碰撞。包围盒法相较于空间网格划分法较为高效但有一定出错的概率。

2.1 碰撞检测的两大阶段

血管形变仿真中的碰撞检测可以分为如图2所示的两个阶段：

(1)检测血管钳是否与血管外壁发生碰撞。若发生，则获取碰撞点的三维坐标和时间t，并对外力f进行分布。

(2)在血管形变过程中检测其内壁是否发生碰撞。若发生，则表明形变已经到达极限，必须固定形变力。

图2 夹持过程中不同阶段钳尖与血管外壁位置

2.2 碰撞检测算法

第一阶段中，使用轴对齐包围盒(AABB)和包围球(Sphere)分别拟合血管壁和血管钳的钳尖部位。即AABB-Sphere碰撞检测问题(图3)。注意，轴对齐包围盒对圆柱体的拟合必然存在空隙，因此实验中必须将钳尖与血管壁垂直。

AABB-Sphere碰撞检测可细分为如下3个步骤：

(1)三维空间下舍弃其中一维度的坐标，将三维对象转化为二维图形。如此，盒子退化为矩形，并用R表示；球体退化为圆形，并用C表示。

图3 AABB-Sphere碰撞检测

(3)从三个维度分别重复步骤1和步骤2，并统计每一维度下的C与R间碰撞检测结果。设六次检测结果分别为O11,O12,O21,O22,O31,O32，其中第一个下标代表舍弃的第i维度，第二个下标代表二维下的第j维度。若上述六个结果在某一时刻全部为真，才认为包围球和包围盒产生了交互。

碰撞检测的第二阶段基于第一阶段。该阶段依靠拟合钳尖的两个球体进行碰撞检测，即Sphere-Sphere碰撞检测问题(图4)。其中，球的半径r与血管壁厚度一致。

图4 Sphere-Sphere碰撞检测

Sphere-Sphere碰撞检测可细分为如下3个步骤：

(1)读取两个包围球的球心O1及O2和半径r1和r2。

(2)计算两个球心之间相对距离的绝对值，即s=|O1-O2|。

(3)判断距离s与半径r1和r2之间的关系。若有s>r1+r2，则两个球体间未发生碰撞。否则发生碰撞。

3 受力-形变算法

在图形学中，物体的几何模型由若干以特定规则相连的顶点构成。任一顶点pi都有指定的几何坐标和法向量。根据各顶点的位置与坐标轴原点O之间的差可以构造顶点位置矢量场p，即：

pi=pi-O

(3)

顶点位置矢量场p可以视为三维空间中的物体在静止状态下一组离散的点。当作用力施加到物体上后，原先位于空间中某个位置的点移动到了新的位置，所有顶点的移动产生了物体位置和形状的改变。而点的位移则使用位移场物u进行描述。

3.1 位移场表达式推导

根据离散信号处理理论[8]和弹性动力学[5]理论，将式(1)分解为单位脉冲函数δ(t)与单位阶跃函数H(t)的线性组合，即：

(4)

图5 算法流程图

图6 血管在不同大小的力度下形变的过程(侧面)

等式左侧为与位置场p及时间t有关的时变位移场。根据点源和场的叠加原理[9]，借助格林求解法对式(4)中的位移场进行求解，有u(pi,t)=D(pi,t)f，其中D(pi,t)即为弹性动力学的格林函数。

由标准正态分布函数的对称性，借助Helmholtz[10]分解法可以将式(2)简化为如下形式：

(5)

图7血管在不同大小的力度下形变的过程(正面)

联合式(4)和式(5)，可得到动力学(1)对应于分布力(2)的解为：

(6)

其中

公式(6)在位移场和形变外力f之间建立了联系，根据该式，空间中所有顶点可通过外力f及其所处位置得到相应的位移矢量。

3.2 形变计算

当血管钳与血管发生碰撞后，借助血管的空间顶点位置矢量场pi和位移场u可以计算血管在受力后的新坐标。即：

pnew=pold+u(pi,t)

(7)

该式表明血管的任一顶点在形变后的空间位置由该点的位置矢量和位移矢量共同决定。通过求解格林函数D(pi,t)，并将以受力载荷为中心的所有顶点的位置矢量进行迭代计算，即可得到血管在受力之后发生的连续形变位移效果。

3.3 程序流程

血管形变仿真流程如图5所示。

4 实验结果与分析

4.1 实验条件与过程

实验在64位Windows7系统下进行。使用C#编程语言和64位2018.3版本的Unity3D图形引擎。硬件为Intel i7-4700 @2. 40GHz，NVIDIA GTX860M和16G内存。

血管形变仿真的工程实现可分为如下六个步骤：

(1)读取OBJ文件，使用结构体存储三维物体的几何信息。

(2)根据血管的几何坐标和原点坐标计算各顶点的位置矢量。

(3)监听鼠标，交互式控制血管钳夹持的位置和夹持力度的大小。

(4)借助Eigen矩阵库求解位移场表达式。

(5)叠加各顶点的位置矢量场和位移场。

(6)重复步骤(4)和步骤(5)，直到位移场为零，结束迭代。

4.2 实验结果展示与分析

图6从侧面角度展示了血管在血管钳夹持下发生的挤压性形变。图6a为夹持之前，血管钳没有与血管接触；图6b和图6c为血管钳夹持过程中血管的形变状态，可以看出血管模型发生了具有较强真实感的形变；图6d为血管上下壁完全闭合时的状态，此时无法继续增大夹持的力度，符合物理直觉且具有良好的视觉效果。

图7从正面角度展示了血管壁在挤压的过程中的形变。从图7b和图7c可以看出血管壁的厚度和体积在形变的过程中始终得到保持，并未发生塌陷；从图7d可以看出血管的上下壁紧密贴合，并未发生穿透。

仿真中的画面每秒传输帧数(FPS)是用于衡量仿真算法实时性的重要指标。表1展示了血管模型在不同数目的顶点下的每秒传输帧数(平均值)。可见随着血管模型顶点数目增加，帧率显著降低，实时性下降。因此在实际应用的过程中，要充分衡量真实感与实时性二者之间的重要程度。

表1 仿真中不同顶点数对应的FPS

5 结 论

为了得到具有较强真实感的交互式仿真效果。使用了动力学方程和正态分布的形变力进行计算。由于血管和血管钳钳尖特殊的几何形状，使用包围盒法进行碰撞检测。根据位移场求解顶点在受力后的新坐标以实现形变效果。实验结果表明，算法基本能够满足系统对血管交互式形变仿真的真实性和实时性要求。但对于其他手术器械对血管的操作仿真将在今后的工作中进一步探讨。