基于sph的虚拟血管造影手术造影剂扩散过程模拟方法
【技术领域】
[0001 ] 本发明涉及一种基于SPH(Smoothed Particle Hydrodynamics)的虚拟血管造影 手术造影剂扩散过程模拟方法，属于计算机仿真建模和虚拟手术领域。
【背景技术】
[0002] 心血管疾病目前是世界上死亡人数第一的疾病，微创介入手术是现代医学手术的 一个巨大突破，其能够有效的治疗各类心血管疾病并且具有出血少、创伤低、恢复快、并发 症少等优点。其中血管造影术是微创介入手术中一个非常重要的部分。在进行微创介入手 术时医生通过观看X-ray产生的医学图像来进行手术，而X-ray会直接穿透血管使得最后在 医学成像设备上无法清晰呈现血管的结构与状态。由于造影剂原子量高、比重大在X-ray下 能够清晰的成像，所以医生通过外部注射器注射造影剂血管进行显影来对血管病变处的血 管结构和病变位置进行诊断。目前有许多研究者对造影剂在血管中的扩散提出了许多方 法：
[0003] Cotin S等在"ICTS，an interventional cardiology training system" (Studies in healthtechnology and informatics ,vol. 70 ,pp. 59-65,2000)一文中米用 层流对血流进行建模，采用对流-扩散模型对造影剂混合与扩撒进行模拟。
[0004] Duriez C等在 "New approaches to catheter navigation for interventional radiology simulation，(Computer Aided Surgery,vol. 11 ,pp.300-308,2006)一文中米 用泊肃叶定律对血流进行建模并采用对流-扩散方程来模拟造影剂在血流中的扩散。
[0005] ffu X 等在 "Real-time modeling of vascular flow for angiography simulation"（Proc.Medicallmage Computing and Computer-Assisted Intervention (MICCAI 07)，pp. 557-565，2007)-文中同样采用泊肃叶定律对血流进行建模，其次采用对 流-扩散方程来模拟造影剂在血流中的扩散。
[0006] F.Wang等在"A computer-based real-time simulation of interventional radiology"（Proceedings of the 29th Annual International Conference of the IEEE EMBS，pp. 108-115,2007)-文中使用血管中心线并采用恒定的速度对造影剂在血管 中进行传播扩散。
[0007] 以上这些方法都采用了简单化的Navier-stokes(NS)方程即泊肃叶定律或者层流 的方式对血流在血管中的流动进行建模，并对造影剂在血管中的扩撒采用对流扩散模型进 行扩散仿真模拟。以上方法忽略了"血液-造影剂"二相混合流体间的动力学交互对造影结 果的影响。由于造影剂具有良好水溶性，当造影注入到血管内后与血液迅速在血液中扩散 开来并与血液混合形成"血液-造影剂"二相混合流体，所以我们采用基于SPH的多相混合流 模型对"血液-造影剂"二相混合流体进行建模并模拟造影剂在血液中扩散造影的过程。
【发明内容】
[0008] 鉴于以上方法所述现有技术存在的不足，本发明的目的在于提供一种基于SPH的 虚拟血管造影手术造影剂扩散过程模拟方法，能够实时准确的模拟造影剂在血液中的造影 扩散过程。
[0009] 为达到上述目的，本发明采用以下技术方案：
[0010] 一种基于SPH的虚拟血管造影手术造影剂扩散过程模拟方法，包括模拟初始化和 模拟循环两个步骤，具体如下：
[0011] (1)模拟初始化:对血管边界条件进行构建，设置血管入口与出口边界条件并对 "血液-造影剂"二相混合流粒子位置进行初始化；
[0012] (2)模拟循环:通过"血液-造影剂"二相混合流体模型并采用边界粒子修正边界处 的混合流密度和应用无滑移边界条件进行循环模拟的步骤依次为：更新空间网格邻居信 息，更新边界粒子属性，更新混合流体粒子每一相的体积分数并对粒子每一相体积分数及 压强值进行修正，计算混合流体粒子的加速度，进行碰撞检测与粒子速度、位置的更新，造 影渲染。
[0013] ⑴模拟初始化
[0014] 首先，对边界粒子的初始化步骤如下：
[0015] 1.以血管模型的每一个三角形面片的边作为单位采样向量，即^和为采 样间隔。因此对每个三角形采样点的位置为：
[0016]
[0017] c\ =/7/,(,(., < 1
[0018] < l
[0019] 其中血和肥为大于〇的正整数。
[0020] 2.对步骤1中采样得到的点沿着采样点对应三角形向外的法向量从1到η倍进行外 推，其中每一次外推的单位长度为采样间隔|，η的最大值满足和为SPH 模型中的光滑核半径。外推得到的点的位置即为边界粒子的位置，外推点对应三角形的法 向量即为外推向量，为距离最近三角形的距离。
[0021] 3.优化通过步骤2得到的边界粒子使其均匀覆盖在均匀覆盖在血管外壁表面。首 先，去除通过外推后位置在血管内部的边界粒子;其次，通过两个相邻粒子的外推向量的内 积?)判断相邻粒子所处的位置，?./ζ 20说明粒子处于相邻三角形间，反之处于血管 的分叉处，再分别对处于分叉和相邻三角形间密集的边界粒子分别采用在半径阈值1和丫2 内进行优化，其中Τ, <€，Τ2 <￡ β通过中间位置的粒子来替换需要进行优化的两个相邻的 粒子，设置新粒子的外推向量为前两个相邻粒子外推向量和的规范化的值，并设置新粒子 的离最近三角形的距离的值为前两个相邻粒子的均值。
[0022] 其次，设置血管入口与出口边界条件血管入口采用入口压强Pir^P血流入口初速 度，出口采用出口压强P?t;初始化"血液-造影剂"二相混合流粒子的位置使其充满血管。 [0023] (2)模拟循环
[0024] 1.更新空间网格邻居信息
[0025]采用空间稀疏哈希的方法对空间模拟粒子的邻居信息进行更新。首先对每个模拟 空间内的粒子根据其位置哈希映射得到其对应的空间网格索引。将该空间网格索引作为 key，模拟的粒子索引作为value进行快速计数排序，排序后在同一网格下的粒子索引即为 连续存储，最后再得到排序后的同一网格下粒子索引的起始位置和结束位置。这样当在进 行邻居查找时，对当前网格和其周围的26个网格，通过起始和结束位置访问排序后的结果 来得到同一网格下的所有粒子索引，即可快速完成邻居查找。
[0026] 2.更新边界粒子属性
[0027] 边界粒子的位置和质量是固定的，边界粒子的密度和压强值为距离最近的混合流 体粒子对应的属性值，边界粒子的速度为：
[0028]
[0029]
[0030] 其中设置β为约束ζ的大小，db为边界粒子距离最近三角形的距离，pb为边界粒子 的位置，Pf为流体粒子的位置，?为边界粒子的外推向量。
[0031] 3.更新混合流体粒子每一相的体积分数并对粒子每一相体积分数及压强值进行 修正
[0032] (a)混合流体中的每一相k的每一个粒子i的连续方程求解体积分数的变化，该连 续方程为：
[0033]
[0034] 其中，ak为粒子中相k的体积分数，而且需要满足1% = ，um为混合流体 t 速度，Ui为相k漂移速度，漂移速度Umk表示为：
[0035]
[0036]其中，τ为控制惯性和压力大小的参数，τ的值越大惯性和压力效果越大，〇为控制 ahpk 扩散速度的参数，σ的值越大其扩散速度越快，％ = 为粒子i相k的质量分数，加速 Pmi = g-((U", ·ν)ιι," + 粒子i相k压强梯度VA与体积分数梯度^%分别表示为： dt ，
[0039] pki = akipmi[0040] !>m = ^P： - P,m)
[0037]
[0038]
[0041]
[0042] 其中，Pkl为粒子i相k的压强，pmi为粒子i混合流体压强，光滑核函数，κ为气态 常数，为粒子j的质量，Α为粒子i的对所有粒子j通过光滑核函数插值得到的插值密 度，其中粒子j包括边界粒子。
[0043] (b)通过(a)中求得体积分数的变化，其后通过前一时刻的体积分数即可求得当前 时刻粒子每一相的体积分数。由于当前时刻粒子i每一相a k不一定满足 k 所以需要对体积分数进行调整。首先对于粒子每一相体积分数ak<0,设置ak = 0,其后再缩 放该粒子中的所有相体积分来满足Σ%=1°其中对于ak2 1或者Xkak2 1的情况，在体积缩 k 放后即可完成调整。由于体积分数参与压力项的计算，所以需要对压力项进行调整，对于修 正后变化的Δ Σ-，所以修正的压强为氣^ = & +Δ/\。. k
[0044] 4.计算混合流体粒子的加速度
[0045] 混合流体中的每一相k的每一个粒子i的动量方程求解流体粒子的加速度，该动量 方程为：
[0046]
[0047] 其中，Α? =_1〃^_为混合流体密度，Pk为相k的静止密度，g为重力加速度，Tm为混 k 合流体的黏滞应力张量，TDm为相间对流的动量转移量，Fe为血管壁对混合流粒子的外力。其 中边界粒子不参与流体体积分数的传递，边界粒子对流体加速度的贡献通过压强的梯度 ▽p,"和粘滞应力张量▽ ·τΜ来满足无滑移边界条件。其中动量方程右边各项分别表示为：
[0055]其中，μL= Skakyk为粒子i中各相粘度yk的聚合粘度，ri为粒子i的位置。f r为流体 粒子与血管内表面碰撞时发生的摩擦力，μ为血管内表面的静摩擦系数，为粒子i碰撞时 的切线速度。Fn为混合粒子碰撞血管壁的法向压力，ri为接触刚度，δ为粒子碰撞时的穿透深 度。fs为血管内表面对流体粒子的粘性力，当它们间的距离小于阈值ds时对粒子才产生相应 的粘性力，△ t为模拟时间步长，vk为血管内表面对各相的粘性系数，cU为粒子i距离表面最 近的距离，η为血管内表面的法向量。fb为血压对流体粒子产生的推力，通过设置入口和出 口压强以及入口血液初速度来描述。
[0056] 5.碰撞检测与粒子速度、位置的更新
[0057] (a)碰撞检测：以模拟的流体粒子的当前位置p?r为原点，以速度方向￥_为射线方 向，以v A t为射线长度，对血管做AABB快速射线检测，对发送碰撞的粒子进行标记。
[0058] (b)位置与速度更新:粒子与血管壁发生碰撞的粒子设置其位置移动到血管壁表 面上，未与血管壁发生碰撞的粒子则通过速度更新粒子至新的位置pnext = pcur+VAt。其中， 为当前时刻的位置，为当前时刻粒子速度，A t为模拟时间间隔。再将所有的粒子的 速度