CN102968811B - 微创血管介入手术中导丝的实时运动仿真方法 - Google Patents

Abstract

本发明涉及一种微创血管介入手术中导丝的实时运动仿真方法，包括以下步骤：1)采用三角形表面网格对血管进行建模；2)采用基尔霍夫弹性杆模型，对导丝进行建模，并利用拉格朗日方程计算导丝在仿真过程产生的形变力；3)进行导丝的渲染；4)采用K-DOPS树进行导丝与血管之间碰撞检测；5)采用非迭代的约束方法计算导丝与血管之间的接触力并进行碰撞响应；6)使用Verlet积分公式来迭代更新导丝的位置与速度；7)利用拉格朗日乘数和快速投影方法来实现导丝的不可伸缩特性；8)使用力反馈设备Phantom Omni来渲染导丝给用户的作用力。与现有技术相比，本发明具有稳定可靠、实现方便、真实高效、工程应用性强等优点。

Description

微创血管介入手术中导丝的实时运动仿真方法

技术领域

本发明涉及一种虚拟手术的模拟方法，尤其是涉及一种微创血管介入手术中导丝的实时运动仿真方法。 

背景技术

随着科学技术的高速发展，虚拟手术的应用为医疗的现代化提供了越来越多的帮助。虚拟的微创血管介入手术系统能够帮助实习医生快速提高手术技能，让更多的心脑血管疾病患者能够进行微创介入治疗。导丝的实时运动仿真作为微创血管介入手术系统中不可缺少的一部分，成为许多学者研究的课题。在导丝实时运动仿真的研究过程中，早期的一部分学者采用几何模型对导丝进行建模，这种方法完全忽略了导丝的物理性质，然而如果完全忽略了导丝的物理性质，与其他模型之间的交互将难以处理。也有一部分学者采用质点弹簧这样简单的物理模型来对导丝进行建模，这种方法虽然计算效率很高，但是由于模型过于简单，根本无法真实地模拟导丝在血管中的运动。本文实现的微创血管介入手术中导丝的实时运动仿真对导丝的模拟有很高的要求，不仅要求导丝要有高的真实性，还对导丝与血管之间交互的模拟要有可靠的物理模型，并且由于要应用于虚拟手术，这就要求模型系统达到高的更新速度，并且能达到实时。但是微创血管介入手术中的导丝的实时运动仿真对血管的模拟要求较低。一个完整的微创血管介入手术中的导丝的实时运动仿真包括了血管模拟、导丝模拟、碰撞检测及响应、导丝的不可伸缩和力反馈渲染。 

对于血管的模拟，由于血管在仿真的过程中几乎不产生形变，只需要采用三角形表面网格这样简单的表面模型对血管建模就足够了。不需要利用四面体网格这样复杂的立体模型来模拟血管。这样有效地降低了计算成本和存储空间。 

导丝的模拟，是微创血管介入手术中导丝的实时运动仿真中的重点部分，算法设计的好坏关系到整个仿真系统的真实性、准确性、高效性。目前，大部分的方法是多刚体运动模型，弹簧模型，有限元模型，B-Spline等理论，这其中弹簧模型和 B-Spline被常用于导丝的模拟，由于其效率高能达到实时。然而弹簧模型过于简单，忽略了导丝很多重要的物理性质而导致真实感很差。B-Spline的缺点是缺乏物理支持，只是从几何的角度来仿真导丝在血管中的形态，无法知道导丝的受力情况。因此无法为后面的力反馈设备提供真实的力反馈。有限元模型能够提供很强的真实性。但是有限元方法有一个致命的弱点：实时性差。不能满足在仿真过程中与用户进行实时交互的要求。为了同时达到真实性和实时性的目标，本文采用了基尔霍夫弹性杆对导丝进行物理建模。 

主流的碰撞检测方法之一是，采用层次结构(hierarchy)为基础核心，例如Octree，k-d树，BSP树等，利用包围体的hierarchy结构也被广泛地应用于各种碰撞检测算法(同时这种方法也被应用于其他领域，如：光线追踪)，包围体的选择，学者们也提出了许多方法，包括，球包围体，有轴向的长方体(AABBs)另外一些碰撞检测方法，包括采用space-time bounds以及采用四维几何体来包围物体，近些年，学者们也提出基于Voronoi图的碰撞检测方法。 

导丝与血管之间的碰撞响应，目前的两种主流算法是惩罚法和约束法。惩罚法的真实性和稳定性差，然而约束法需要进行多次迭代才能确保导丝不刺穿血管壁，计算成本很高。结合惩罚法与约束法的优点，本文采用一种非迭代的约束方法来计算导丝与血管之间的接触力并进行碰撞响应。 

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种稳定可靠、实现方便、真实高效的微创血管介入手术中导丝的实时运动仿真方法。 

本发明的目的可以通过以下技术方案来实现： 

一种微创血管介入手术中导丝的实时运动仿真方法，其特征在于，该方法包括以下步骤： 

1)由于血管在仿真的过程中不产生形变，采用三角形表面网格对血管进行建模； 

2)采用基尔霍夫弹性杆模型，对导丝进行建模，并利用拉格朗日方程计算导丝在仿真过程产生的形变力； 

3)进行导丝的渲染，相邻的两导丝质点之间绘制一个圆柱体，同时在相邻的两圆柱体之间绘制一个圆球进行连接； 

4)采用K-DOPS树进行导丝与血管之间碰撞检测； 

5)采用非迭代的约束方法计算导丝与血管之间的接触力并进行碰撞响应； 

6)使用Verlet积分公式来迭代更新导丝的位置与速度； 

7)利用拉格朗日乘数和快速投影方法来实现导丝的不可伸缩特性； 

8)使用力反馈设备Phantom Omni来渲染导丝给用户的作用力。 

所述的采用基尔霍夫弹性杆模型，对导丝进行建模，并利用拉格朗日方程计算导丝在仿真过程产生的形变力的具体实现过程为： 

21)将导丝描述成一条被坐标化的曲线Г＝{t，m₁，m₂}，其中t表示沿着曲线的单位切线，m₁，m₂表示曲线的两条单位法向量，该三个正交向量构成一个局部移动坐标系{t(s)，m₁(s)，m₂(s)}，其中s表示曲线的弧长； 

22)利用连续的基尔霍夫弹性杆模型计算被坐标化的导丝的连续弹性能量E(Г)，其中E(Г)＝E_弯曲(Г)+E_扭曲(Г)； 

23)利用离散的基尔霍夫弹性杆模型对连续的弹性能量进行离散化处理，推导出导丝的离散弹性能量的公式； 

24)利用拉格朗日方程对导丝的离散弹性能量进行导丝位置的求导，计算出导丝的形变力。 

所述的连续弹性能量E(Г)＝E_弯曲(Г)+E_扭曲(Г)的具体公式为： 

其中二维向量ω＝(ω₁，ω₂)^T表示导丝在局部移动坐标系中的曲率向量，m＝m₁′m₂表示导丝在局部移动坐标系中的扭曲率；其中α表示导丝的弯曲模数，β表示导丝的扭曲系数；ω₁，ω₂分别表示导丝在局部移动坐标系m₁，m₂上的曲率，ds表示对导丝的弧长进行积分。 

所述的离散弹性能量E(Г)＝E_弯曲(Г)+E_扭曲(Г)的具体公式为： 

其中i＝1，2，3，…n表示将导丝进行离散分段后的段号，表示第i段导丝的长度， θ_i表示导丝的第i个顶点的旋转角度，ω_i表示第i段导丝在局部移动坐标系中的曲率向量，m_i表示第i段导丝在局部移动坐标系中的扭曲率。 

所述的采用非迭代的约束方法的具体实现过程为： 

首先计算导丝在血管中的穿刺向量d；然后根据穿刺向量计算接触力

最后利用步骤6)提到的Verlet积分公式来进行导丝的位置更新，确保导丝仅在血管内部运动，不穿透血管壁。 

所述的快速投影方法的具体实现过程为： 

71)构造导丝的长度约束函数C(p₁，p₂)＝|p₁，p₂|-d，其中p₁，p₂为相邻两导丝质点的位置坐标，d为该相邻两质点的原始长度； 

72)如果长度约束函数的值超过了最低阈值，则对质点p₁，p₂进行位置修正并且循环执行该过程直至长度约束函数的值小于最低阈值，其中修正公式为： 因此p₁，p₂的修正量Δp₁，Δp₂为： 

${\mathrm{\Δp}}_1=-\frac{1}{2}\left(\right|p_1-p_2|-d)\frac{p_1-p_2}{|p_1-p_2|}$

${\mathrm{\Δp}}_1=+\frac{1}{2}\left(\right|p_1-p_2|-d)\frac{p_1-p_2}{|p_1-p_2|}$

所述的力反馈渲染的具体实现方法为： 

导丝的仿真频率为50Hz，而力反馈设备的更新频率要达到1000Hz才能提供平滑的力触觉，为了实现两者的同步，采用线性插值的方法来计算力反馈设备的力输入参数。 

与现有技术相比，本发明具有以下优点： 

(1)采用基尔霍夫弹性杆模型，作为导丝的模拟基础，并且采用Verlet更新方法，有效地提高了数值的稳定性和精确性，实现了真实高效的更新； 

(2)采用非迭代的约束方法来计算导丝与血管之间的接触力并进行碰撞响应，该方法避免了传统的约束方法需多次迭代的问题。 

(3)采用了K-DOPS树，提高了碰撞检测的效率，提高了微创血管介入手术中导丝的实时运动仿真的真实性和有效性。 

(4)利用拉格朗日乘数和快速投影方法来实现导丝的不可伸缩特性，该方法能够快速、有效地实现导丝的不可拉伸。 

附图说明

图1为本发明的流程示意图； 

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。 

实施例 

本实施例在一台拥有Pentium(R)Dual-Core CPU、一个NVIDIA GeForce GT 

240和2G内存的个人计算机中实现，整个虚拟手术训练系统采用C++语言编写，采用QT编写界面。系统在实现的过程中使用的第三方库有OpenGL、Boost和CGAL。 

本实施例的实施流程如图1所示，一种微创血管介入手术中导丝的实时运动仿真方法，具体包括以下步骤： 

第一步，采用三角形表面网格，对血管进行建模。所述的血管三角形表面网格，首先是从人体CT数据中分割出STL格式的血管数据文件，然后利用Netgen生成三角形表面网格，最后对三角形表面网格进行平滑处理生成最终系统可用的三角形表面网格。 

第二步，采用基尔霍夫弹性杆模型，对导丝进行建模，并利用拉格朗日方程计算导丝在仿真过程产生的形变力。所述的采用基尔霍夫弹性杆模型，对导丝进行建模，并利用拉格朗日方程计算导丝在仿真过程产生的形变力的具体实现过程为： 

将导丝描述成被一条被坐标化的曲线Г＝{t，m₁，m₂}，其中t表示沿着曲线的单位切线，m₁，m₂表示曲线的两条单位法向量，该三个正交向量构成一个局部移动坐标系{t(s)，m₁(s)，m₂(s)}，其中s表示曲线的弧长。利用连续的基尔霍夫弹性杆模型计算被坐标化的导丝的连续弹性能量E(Г)，其中E(Г)＝E_弯曲(Г)+E_扭曲(Г)。然后再利用离散的基尔霍夫弹性杆模型对连续的弹性能量进行离散化处理，推导出导丝的离散弹性能量的公式。最后利用拉格朗日方程对导丝的离散弹性能量进行导丝位置的求导，推导出导丝的形变受力公式。 

所述的连续弹性能量E(Г)＝E_弯曲(Г)+E_扭曲(Г)的具体公式为： 

其中二维向量ω＝(ω₁，ω₂)^T表示导丝在局部移动坐标系中的曲率向量，m＝m′₁m₂表示导丝在局部移动坐标系中的扭曲率；其中α表示导丝的弯曲模数，β表示导丝的扭曲系数。ω₁，ω₂分别表示导丝在局部移动坐标系m₁，m₂上的曲率，ds表示对导丝的弧长进行积分。 

所述的离散弹性能量E(Г)＝E_弯曲(Г)+E_弯曲(Г)的具体公式为： 

其中i＝1，2，3，…n表示将导丝进行离散分段后的段号，表示第i段导丝的长度， 

θ_i表不导丝的第i顶点的旋转角度。 

第三步，进行导丝的渲染，相邻的两导丝质点之间绘制一个圆柱体，同时在相邻的两圆柱体之间绘制一个小圆球进行连接。 

第四步，采用K-DOPS树进行导丝与血管之间的碰撞检测。碰撞检测中的基元是三角面片与直线段(三角面片用于模拟血管，直线段用于模拟导丝)：对三角面片和直线段，在进行碰撞检测和响应时，视为三角棱柱和两头有相同半径的半球状的圆柱体，当线段进入圆柱体，即判定二者相互碰撞。 

第五步，采用非迭代的约束方法计算导丝与血管之间的接触力并进行碰撞响应。所述的采用非迭代的约束方法的具体实现过程为：首先计算导丝在血管中的穿刺向量d，然后根据穿刺向量计算接触力最后利用Verlet积分公式来进行导丝的位置更新，确保导丝仅在血管内部运动，而不穿透血管壁。 

第六步，使用Verlet积分公式来迭代更新导丝的位置与速度。本文采用的是速度Verlet积分公式，具体的速度Verlet积分公式如下： 

$x(t+\mathrm{\Δt})=x\left(t\right)+v\left(t\right)\mathrm{\Δt}+\frac{1}{2}a\left(t\right){\mathrm{\Δt}}^2,$

$v(t+\mathrm{\Δt})=v\left(t\right)+\frac{a\left(t\right)+a(t+\mathrm{\Δt})}{2}\mathrm{\Δt}.$

第七步，利用拉格朗日乘数和一种快速投影方法来实现导丝的不可伸缩特性。所述的快速投影方法的具体实现过程为： 

701)构造导丝的长度约束函数C(p₁，p₂)＝|p₁，p₂|-d，其中p₁，p₂为导丝中相邻两质点的位置坐标，d为该相邻两质点的原始长度。 

702)如果长度约束函数的值超过了最低阈值，则对质点p₁，p₂进行位置修正并且循环执行该过程直至长度约束函数的值小于最低阈值，其中修正公式为： 因此p₁，p₂的修正量ΔP₁，ΔP₂为： 

${\mathrm{\Δp}}_1=-\frac{1}{2}\left(\right|p_1-p_2|-d)\frac{p_1-p_2}{|p_1-p_2|},$

${\mathrm{\Δp}}_1=+\frac{1}{2}\left(\right|p_1-p_2|-d)\frac{p_1-p_2}{|p_1-p_2|}$

第八步，力反馈渲染，使用力反馈设备Phantom Omni渲染导丝给用户的作用力。所述的力反馈渲染的具体实现方法为：导丝的仿真频率只有50Hz，而力反馈设备的更新频率只要要达到1000Hz才能提供平滑的力触觉，为了实现两者的同步，我们采用线性插值的方法来计算力反馈设备的力输入参数。即对系统计算的相邻时间间隔(20ms)的力f₁，f₂进行线性插值(即20等分)，则力反馈设备的力输入参数为f_i＝f₁+(f₂-f₁)i/20。 

Claims (5)

1.一种微创血管介入手术中导丝的实时运动仿真方法，其特征在于，该方法包括以下步骤：

1)由于血管在仿真的过程中不产生形变，采用三角形表面网格对血管进行建模；

2)采用基尔霍夫弹性杆模型，对导丝进行建模，并利用拉格朗日方程计算导丝在仿真过程产生的形变力；

3)进行导丝的渲染，相邻的两导丝质点之间绘制一个圆柱体，同时在相邻的两圆柱体之间绘制一个圆球进行连接；

4)采用K-DOPS树进行导丝与血管之间的碰撞检测；

5)采用非迭代的约束方法计算导丝与血管之间的接触力并进行碰撞响应；

6)使用Verlet积分公式来迭代更新导丝的位置与速度；

7)利用拉格朗日乘数和快速投影方法来实现导丝的不可伸缩特性；

8)使用力反馈设备Phantom Omni来渲染导丝给用户的作用力；

采用非迭代的约束方法的具体实现过程为：

首先计算导丝在血管中的穿刺向量d；

然后根据穿刺向量计算接触力M为导丝的质量。

2.根据权利要求1所述的一种微创血管介入手术中导丝的实时运动仿真方法，其特征在于，所述的采用基尔霍夫弹性杆模型，对导丝进行建模，并利用拉格朗日方程计算导丝在仿真过程产生的形变力的具体实现过程为：

21)将导丝描述成一条被坐标化的曲线Γ＝{t,m₁,m₂}，其中t表示沿着曲线的单位切线，m₁,m₂表示曲线的两条单位法向量，该三个正交向量构成一个局部移动坐标系{t(s),m₁(s),m₂(s)}，其中s表示曲线的弧长；

22)利用连续的基尔霍夫弹性杆模型计算被坐标化的导丝的连续弹性能量E(Γ)，其中E(Γ)＝E_弯曲(Γ)+E_扭曲(Γ)；

23)利用离散的基尔霍夫弹性杆模型对连续的弹性能量进行离散化处理，推导出导丝的离散弹性能量的公式；

24)利用拉格朗日方程对导丝的离散弹性能量进行导丝位置的求导，计算出导丝的形变力。

3.根据权利要求2所述的一种微创血管介入手术中导丝的实时运动仿真方法，其特征在于，所述的连续弹性能量E(Γ)＝E_弯曲(Γ)+E_扭曲(Γ)的具体公式为：

其中二维向量ω＝(ω₁,ω₂)^T表示导丝在局部移动坐标系中的曲率向量，m＝m′₁m₂表示导丝在局部移动坐标系中的扭曲率；其中α表示导丝的弯曲模数，β表示导丝的扭曲系数；ω₁,ω₂分别表示导丝在局部移动坐标系m₁,m₂上的曲率，ds表示对导丝的弧长进行积分。

4.根据权利要求3所述的一种微创血管介入手术中导丝的实时运动仿真方法，其特征在于，所述的离散弹性能量E(Γ)^/＝E_弯曲(Γ)^/+E_扭曲(Γ)^/的具体公式为：

其中i＝1,2,3,…n表示将导丝进行离散分段后的段号，表示第i段导丝的长度，θ_i表示导丝的第i个顶点的旋转角度，ω_i表示第i段导丝在局部移动坐标系中的曲率向量，m_i表示第i段导丝在局部移动坐标系中的扭曲率。

5.根据权利要求1所述的一种微创血管介入手术中导丝的实时运动仿真方法，其特征在于，所述的力反馈渲染的具体实现方法为：

导丝的仿真频率为50Hz，而力反馈设备的更新频率要达到1000Hz才能提供平滑的力触觉，为了实现两者的同步，采用线性插值的方法来计算力反馈设备的力输入参数。