本发明涉及一种虚拟手术的模拟方法，尤其是涉及到一种微创血管介入手术中导丝受血流作用的实时运动仿真方法。
背景技术：
随着科学技术的高速发展，虚拟手术的应用为医疗的现代化提供了越来越多的帮助。虚拟的微创血管手术训练系统可以快速有效地帮助实习医生提高手术技能，让更多的心脑血管疾病患者可以接受微创介入治疗。导丝的实时运动仿真是微创血管介入手术系统最重要的一部分，也是众多学者研究的课题。
在导丝实时运动仿真的研究过程中，早期的学者受限于当时硬件的计算水平采用非物理的几何模型，例如B样条曲线，这种非物理模型从几何角度描述导丝形态，虽然计算速度快，但是完全忽略了导丝的物理特性，也无法描述导丝的受力情况，所以真实感比较差；也有其他学者采用基于物理的质点弹簧模型、弹性杆模型等对导丝的运动仿真进行建模，这一类物理模型考虑了导丝的物理特性，且易于对导丝与血管之间的交互进行物理建模，在硬件水平不断发展之后也能达到实时计算的要求，并且相比于非物理模型，在真实感上有比较大的提升，被广泛采用，例如Basdogan等学者在《Virtual Environments for Medical Training:Graphical and Haptic Simulation of Laparoscopic Common Bile Duct Exploration》一文中采用的质点弹簧模型、Tang等学者在《A realistic elastic rodmodel for real-time simulation of minimally invasive vascular interventions》一文中采用的弹性杆模型等。
然而，虽然手术中导丝是在血流中运动的，大多数学者却都没有考虑导丝在血流中运动时手动的血流作用；微创血管介入手术训练系统要为用户提供一个高度还原的具有沉浸感的虚拟训练环境，这就要求仿真的导丝具有很高的物理可信度，因此在考虑导丝本身的物理属性的基础上，还需要对血管内的血流进行计算，分析导丝在运动中收到的血流作用情况。
要对受血流作用的导丝实时仿真主要的挑战有以下几点：
1)由于流体方程本身的复杂性，流体计算通常计算复杂且需要消耗大量的时间；
2)血管是不规则的三维模型，在不规则的模型中建立流体模型比较困难，并且会增加流体计算的复杂度；
3)导丝的物理仿真要求较高的实时性，整合流体计算必然导致仿真效率的降低，甚至难以达到实时；
技术实现要素：
本发明的目的就是为了克服上述现有技术存在的缺陷，通过简化流体方程、不规则血管的离散化描述、采用硬件加速等手段，提供了一种计算稳定、实时性好、真实感强、应用灵活的微创血管介入手术中导丝受血流作用的实时运动仿真方法。
本发明的目的可以通过以下技术方案来实现：
一种微创血管介入手术中导丝受血流作用的实时运动仿真方法，包括以下几个步骤：
第一步，根据病人的CT扫描图像预先分割出血管部位的数据，采用三角形表面网格对血管进行建模；
第二步，求取血管网格的中心线，按照中心线将连续的血管网格模型离散化表示为一组圆柱体状小血管的组合，并抽象为一个有向图；
第三步，根据泊松叶法则对每一段小血管建立流体模型，用于描述每一段小血管中流量与两端压力差的关系；
第四步，根据小血管之间的位置关系、进出的流量一致和小血管两端的压力差这些条件，对抽象为有向图的血管模型整体建立流体计算矩阵，计算血管中的流量和压力分布；
第五步，使用弹性杆模型对导丝建模，并根据导丝在血管中的相对位置，计算血流对导丝的作用力并施加到导丝的仿真模型上，更新导丝的形态；
第六步，对导丝进行三维渲染。
根据中心线将连续的网格模型离散化表示为一组圆管的集合，并抽象为一个有向图，表示为G(Nn，Ne)，其中包含了Nn个节点和Ne条边，节点为离散点，边为圆管。
所述第三步中根据泊松叶法则对每一段圆柱体状的小血管中的流体建模如下：
其中，Q表示这段血管中的流量，Δp表示血管两端的压力差，L表示该段血管的长度，r表示血管的半径，η表示血液的粘度系数，R表示该段血管的阻力，由L、r、η描述；
所述的第四步中血管整体流体计算过程如下:
41)使用节点-边矩阵描述抽象为有向图后的血管的拓扑结构：
其中，如果第j条边从第i个节点出发，则Ai,j＝1，如果第j条边到第i个节点结束，则Ai,j＝-1；
42)由于每个节点流入流量与流出流量是相等的，因此表示为
AQ＝0
其中Q表示这段血管中的流量，A为描述血管拓扑结构的节点-边矩阵；
43)又因为每一段小血管都中的流体都有3种所示的模型，整个血管的整体流体计算表示为：
CQ＝ΔP
其中，C为对角阵Ri表示第i段血管的阻力，其中i＝1···n，ΔP是一个列向量，包括了每一条边的压力差；
44)使用一个正弦函数表示心脏输出的血压在收缩压和舒张压之间的变化psource(t)
其中，P1和P2分别为收缩压和舒张压，ω用来控制输出压力的周期；
45)整合前面列出的方程，整个血管整体血流计算过程表示为：
其中，是一个稀疏的可逆矩阵，P表示心脏输出的压力，方程的求解方式如下
由于矩阵K的值只与血管的拓扑结构有关，而血管的结构是固定不变的，因此K的逆在开始仿真之前提前求好，以提高仿真的速度。
所述的第五步中计算血流对导丝的作用力并施加到导丝的仿真模型的步骤如下：
51)查找导丝质点在血管中的相对位置，确定其所属的小血管；将离散化的血管描述为V1～n，所有的导丝质点描述为m1～n，为了加快查找的速度，不针对每一个质点搜索所有的血管，而只对第一个质点m1搜索所有的血管并找到对应的血管V1′，而对于后续的任意一个质点mi，只需要搜索上一个质点所属的血管V(i-1)′的一个邻域V(i-1)′±n，即快速找到每个质点所在的血管；使用这种查找方法有效地提高仿真的速度；
52)根据质点在血管中的位置，求取血流在质点上的作用力，表示如下
x表示质点在血管轴向方向的距离，d表示质点在血管径向方向的距离，Q表示这段血管中的流量，ΔP表示血管两端的压力差，r表示血管的半径，η表示血液的粘度系数；
53)将血流的作用整合进导丝的物理仿真中，基于弹性杆模型的导丝每一帧的计算如下
其中，Δx和Δv分别表示各个质点的位置和速度变化，Δt是时间间隔，f表示各个质点的受力情况的雅可比矩阵，包括导丝受到的推拉、旋转以及碰撞的力，M是依赖于导丝各个质点的位置信息的一个对角矩阵，vt,xt分别表示在t时刻质点的速度和位置；
在每个质点上都加入血流的作用如下，得到受血流作用影响的导丝仿真计算模型：
与现有技术相比，本发明将血流的作用力整合进了导丝的运动仿真中，具有计算稳定、实时性好、真实感强、应用灵活等优点。
附图说明
图1为本发明的流程图。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。
实施例
本实施例在一台PC上实现，CPU为英特尔至强E3-1230，内存4G，显卡为英伟达GTX760，实现方式为C++编程语言实现。
本实施例的实施流程如图1所示。
第一步，根据病人的CT扫描图像，使用区域增长加手工调整的方法分割出血管部位的数据，通过Marching Cubes算法采用三角形表面网格对血管进行建模；
第二步，求取血管网格的中心线，按照中心线将连续的血管网格模型离散化表示为一组圆柱体状小血管的组合，本实施例中的血管离散化后总共有2477个圆柱体状小血管，然后使用一个有向图抽象描述；
第三步，根据泊松叶法则对每一段小血管建立流体模型，描述每一段小血管中流量与两端压力差的关系；
第四步，根据小血管之间的位置关系、进出的流量一致和小血管两端的压力差这些条件，对抽象为有向图的血管模型整体建立流体计算矩阵，计算血管中的流量和压力分布；采用英特尔数学核心运算库(Intel MKL)加速血流计算。
第五步，根据导丝在血管中的相对位置，计算血流对导丝的作用力并施加到导丝的仿真模型上；
第六步，使用OpenGL对导丝进行三维渲染。