CN107411819B - 球囊血管成形手术过程实时模拟方法 - Google Patents

球囊血管成形手术过程实时模拟方法 Download PDF

Info

Publication number
CN107411819B
CN107411819B CN201710234745.8A CN201710234745A CN107411819B CN 107411819 B CN107411819 B CN 107411819B CN 201710234745 A CN201710234745 A CN 201710234745A CN 107411819 B CN107411819 B CN 107411819B
Authority
CN
China
Prior art keywords
collision
balloon
tetrahedral
mesh
simulation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710234745.8A
Other languages
English (en)
Other versions
CN107411819A (zh
Inventor
唐鹏斌
王音
宫若彬
黄东晋
谢志峰
丁友东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Shanghai for Science and Technology
Original Assignee
University of Shanghai for Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Shanghai for Science and Technology filed Critical University of Shanghai for Science and Technology
Priority to CN201710234745.8A priority Critical patent/CN107411819B/zh
Publication of CN107411819A publication Critical patent/CN107411819A/zh
Application granted granted Critical
Publication of CN107411819B publication Critical patent/CN107411819B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B34/00—Computer-aided surgery; Manipulators or robots specially adapted for use in surgery
    • A61B34/10—Computer-aided planning, simulation or modelling of surgical operations
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B34/00—Computer-aided surgery; Manipulators or robots specially adapted for use in surgery
    • A61B34/10—Computer-aided planning, simulation or modelling of surgical operations
    • A61B2034/101—Computer-aided simulation of surgical operations
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B34/00—Computer-aided surgery; Manipulators or robots specially adapted for use in surgery
    • A61B34/10—Computer-aided planning, simulation or modelling of surgical operations
    • A61B2034/101—Computer-aided simulation of surgical operations
    • A61B2034/102—Modelling of surgical devices, implants or prosthesis
    • A61B2034/104—Modelling the effect of the tool, e.g. the effect of an implanted prosthesis or for predicting the effect of ablation or burring

Landscapes

  • Health & Medical Sciences (AREA)
  • Surgery (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biomedical Technology (AREA)
  • Robotics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提出了一种球囊血管成形手术过程实时模拟方法。该方法包括球囊模拟和球囊与血管碰撞交互模拟两部分。球囊模拟部分，采用基于位置的方法并结合连续介质材料模型对球囊材质进行模拟；球囊与血管碰撞交互模拟部分，通过构建自适应空间碰撞网格来加速球囊和血管之间的碰撞检测。本方法能够实时准确的模拟球囊与血管的交互扩张过程。

Description

球囊血管成形手术过程实时模拟方法
技术领域
本发明涉及一种球囊血管成形手术过程实时模拟方法，属于计算机仿真建模和虚拟手术 领域。
背景技术
心血管疾病目前是世界上死亡人数第一的疾病，其中动脉粥样硬化是心血管疾病中最常 见且危害较大的一种疾病。这种心血管疾病大都是胆固醇和胆固醇吸附细胞产生的废物、钙 质和纤维蛋白等物质堆积在血管内部所形成的。这些斑块的堆积造成血管阻塞使得血流量的 减少，从而导致心脏死亡或者脑死亡等急性病发死亡。球囊成形术是一种治疗这种动脉阻塞 常见且实用的微创血管介入手术方法。这种治疗手段通过股动脉经皮穿刺在X-Ray下将导丝 推送到狭窄血管的部位，随后将球囊导管顺着导丝到达病处后通过外部的加压注射器使球囊 膨胀，膨胀的球囊扩大被粥样斑块阻塞的动脉，通过这样的治疗手段即可使扩张后的动脉血 管的血流量恢复到一个可接受的范围。这种实际的操作往往都是在x-ray所形成的影像下完成 的具有很高的难度，因此这种手术在实际对人进行手术前的技能训练是很有必要的。现在随 着虚拟现实技术的不断发展，采用VR技术对手术的训练越来越有效而且便捷。在进行虚拟 手术的训练中，虚拟手术的真实性是影响训练效果的关键，因此需要采用真实的物理模型来 对整个过程进行模拟。在整个球囊成形手术过程中，球囊和血管的扩张的结果很大程度上依 赖于整个过程中的各个部分的几何结构和物理属性，对气囊的建模和其对血管模型的碰撞交 互是球囊成形术中至关重要的一部分。目前有许多研究者提出了许多相关的方法：
Luboz V等在“Real-time stent and balloon simulation for stenosistreatment”(The Visual Computer,2014,30(3):341-349)一文中描述了支架以及球囊的仿真，在其系统中球囊表示为 一系列重叠的粒子依附在导丝上，使用这些粒子的扩张来表示球囊的扩展，当粒子扩张并与 血管发生碰撞后则施加额外的外力来约束球囊在血管内部。这种采用粒子来模拟球囊的方法 不能精确的模拟球囊扩张以及缩紧的动力学过程，并且粒子也无法详细表示气囊与血管的碰 撞交互过程。另一方面，其系统中的球囊与血管的动力学交互碰撞模块采用多个不同预处理 后的血管作为球囊膨胀后的碰撞响应，即当球囊粒子在血管固定的狭窄处扩张到某个阈值后 血管采用下一个层次的血管结构作为碰撞交互的结果。这种采用不同层次血管的碰撞响应的 方法虽然可以加速系统的速度，但是生成这些不同层次的血管比较繁琐而且只能模拟生成的 这个狭窄血管的固定位置，并且预生成的这几种不同状态的血管也不一定符合球囊当前的膨 胀状态，同样这种碰撞响应的方式忽略了球囊和血管相互碰撞间的动力学交互的过程。
Wang Q等在“Cute balloons with thickness”(International Conference onImage and Graphics.Springer International Publishing,2015:75-89)一文中采用三棱柱元素构建球囊的表 面并基于FEM方法求解气囊膨胀的大变形过程，这种基于FEM的方法虽然能够精确的模拟 薄膜的非线性形变但是往往很难达到实时。
Skouras M等在“Computational design of rubber balloons”(ComputerGraphics Forum. Blackwell Publishing Ltd,2012,31(2pt4):835-844)一文中橡胶材质的气球采用连续材料介质力 学的方法对气球的力学性质进行建模，并通过实验数据拟合出所需要的材料模拟对橡胶材质 的气球进行模拟。
Skouras M等在“Designing inflatable structures”(ACM Transactions onGraphics(TOG), 2014,33(4):63)一文中通过张力场理论来模拟充气薄膜起皱的状态。
Bridson等在“Robust treatment of collisions,contact and friction forcloth animation”(ACM Transactions on Graphics(TOG),2002,21(3):594-603)一文中采用连续碰撞检测(CCD)提出 了稳定的碰撞处理的方法来模拟布料网格间的碰撞，随后既有很多不同的方法被提出来对该 方法进行优化。在某些情况下可变形物体是不可避免的与其他物体发生碰撞并紧密接触，例 如球囊和血管网格的相互碰撞和挤压，这种情况很容易导致这些物体的网格边界进行纠缠状 态，这种CCD方法不能很好的解决种问题。
Selle等在“Robust high-resolution cloth using parallelism,history-based collisions,and accurate friction”(IEEE Transactions on Visualizationand Computer Graphics,2009,15(2): 339-350)一文中采用基于历史的方法来求解施加在碰撞对的力，但是这种方法在每一帧仿真 前需要重新初始化网格状态。
Müller等在“Air meshes for robust collision handling”(ACM Transactionson Graphics (TOG),2015,34(4):133)一文中通过初始化生成模拟的空间网格之间的四面体碰撞网格，并 通过约束四面体碰撞网格来加速碰撞检测。
Si H在“TetGen,a Delaunay-based quality tetrahedral mesh generator”(ACM Transactions on Mathematical Software(TOMS),2015,41(2):11)一文中提出一个新的四面体网格生成工具 TetGen。
Shewchuk J R在“Two discrete optimization algorithms for thetopological improvement of tetrahedral meshes”(Unpublished manuscript,2002,65.)一文中采用边和面的增删的方法来对 四面体网格进行质量结构优化。
Müller M等在“Position based dynamics”(Journal of Visual Communicationand Image Representation,2007,18(2):109-118)一文中提出了基于位置的动力学约束的方法来实时模拟 各种形变物体。
Bender J等在“Position-based simulation of continuous materials”(Computers&Graphics, 2014,44:1-10)一文中提出基于位置的能量递减方法来模拟连续材料介质的固体以及布料。
发明内容
鉴于以上已有的球囊扩张方法的不足，以及目前相关球囊物理建模和碰撞检测方法的深 入研究，本发明的目的在于提供一种球囊血管成形手术过程实时模拟方法，能够实时准确的 模拟球囊与血管的交互扩张过程。
为达到上述目的，本发明采用以下技术方案：
一种球囊血管成形手术过程实时模拟方法，包括球囊模拟和球囊与血管的碰撞交互模拟， 具体如下：
(一)球囊模拟：采用基于位置的方法并结合连续介质材料模型对球囊材质进行模拟；
(二)球囊与血管碰撞交互模拟：通过构建自适应空间碰撞网格来加速球囊和血管之间 的碰撞检测。
(一)球囊模拟
(1-1)基于位置方法是一种采用约束函数进行求解动力学的一种方法其包括模拟初始化 和模拟循环两部分。
(1-1-1)模拟初始化：初始化仿真模拟所有顶点的位置、速度、质量以及相关初始化约束条件。
(1-1-2)模拟循环：循环模拟初始化设定后的所有顶点，其中模拟循环分为如下若干步骤：
(a)通过时间积分预求每一个顶点的新速度
Figure RE-GDA0001440107500000031
和位置xt＝xt-1+Δtvt，其中vt和vt+1分别表示当前模拟时 刻和上一模拟时刻顶点的速度，Δt为模拟的时间步长，M为质量对角矩 阵，
Figure RE-GDA0001440107500000032
为当前时刻的外力，xt和xt-1分别表示当前模拟时刻和上一 模拟时刻顶点的位置。
(b)约束函数C(x)来约束求得每一个顶点的修正位移
Figure BDA0001267611290000034
使其修 正后的顶点位置xt＝xt+Δx满足约束函数，其中
Figure BDA0001267611290000035
为每个顶点质量 的倒数，
Figure BDA0001267611290000041
为拉格朗日乘子，
Figure BDA0001267611290000042
为约束函数的梯 度。其后通过并行的jacobi迭代约束其收敛稳定。
(c)通过当前时刻和上一时刻的位置更新每个顶点的速度为
Figure BDA0001267611290000043
(d)碰撞检测与响应。
减弱每个顶点的速度使系统稳定。
(1-2)通过线性拉格朗日形状函数并采用闭合的三角形网格对气囊模型进行离散求 解。
对于一个封闭的球囊三角形网格模型，
Figure BDA0001267611290000044
表示第i个未变形的顶点的位置,xi表示第i 个变形后的顶点的位置，对应未变形状态下的边矢量
Figure BDA0001267611290000045
和变形状态下的边矢量 eij＝xj-xi。对一个初始状态的三角形t其三个顶点分别为
Figure BDA0001267611290000046
和
Figure BDA0001267611290000047
我们首先构建该三 角形的材料结构矩阵为
Figure BDA0001267611290000048
其中
Figure BDA0001267611290000049
和
Figure BDA00012676112900000410
分别为：
Figure BDA00012676112900000411
Figure BDA00012676112900000412
对一个三角形体元其变形梯度
Figure BDA00012676112900000413
表示为:
Figure BDA00012676112900000414
其中，Ds为变形状态下的矩阵，Dm为未变形状态下的参考矩阵。这两个矩阵分别表示为：
Ds＝[e31 e32]
Figure BDA00012676112900000415
因此2×2的右柯西格林张量即可定义为：
C＝FTF
第一和第三柯西格林不变量张量分别为：
I1＝tr(C)
I3＝det(C)
这里我们假设模拟的材料不存在横向剪切应力而且材料满足不可压条件，因此对应当 前状态下厚度表示为：
Figure BDA0001267611290000051
其中，H为初始未变形状态下的厚度，I3的物理意义表示为体元的体积变化。
由于不同材料其具有不用的性质，并且其材料的能力密度函数也不相同。球囊是一种 具有不可压性且超弹性的薄膜体，目前气囊大都使用聚氨酯材料，且目前大部分基于有 限元的数值模拟的方法采用Mooney-Rivlin模型来模拟气囊的这种聚氨酯橡胶材料。该Mooney-Rivlin材料的应变能量密度函数为：
Figure BDA0001267611290000052
Figure BDA0001267611290000053
Figure BDA0001267611290000054
其中，
Figure BDA0001267611290000055
和
Figure BDA0001267611290000056
为第一个第二拉伸不变性,参数C1和C2为聚氨酯橡胶材料的参数，λ1和λ2为”面内”主要拉伸量，该拉伸量为右柯西格林张量C的特征值，并分别表示为：
Figure BDA0001267611290000057
Figure BDA0001267611290000058
λ3为厚度拉伸量可通过不可压条件求得：
λ1λ2λ3＝1
Figure BDA00012676112900000510
变形梯度，我们即可求得一个三角形元的应变能量为：
E(x)＝Ψ(C)·hA
其中，h为当前状态下三角形元的厚度，A为其当前状态下的面积。
我们采用Bender J一文中的基于位置的能量减弱方法进行求解，因此该能量约束方程 为C(x)＝E(x)＝0，通过求解能量梯度来得到三角形元每个顶点的移动值，该能量梯度描 述为物体的内力场的，该梯度表达为：
Figure BDA0001267611290000059
其中第一皮奥拉-柯克霍夫应力张量P(F)定义为：
Figure BDA0001267611290000061
其中，对于I1有：
Figure BDA0001267611290000062
Figure BDA0001267611290000063
Figure BDA0001267611290000064
Figure BDA0001267611290000065
Figure BDA0001267611290000066
对于I3有：
Figure BDA0001267611290000067
Figure BDA0001267611290000068
Figure BDA0001267611290000069
Figure BDA00012676112900000610
Figure BDA00012676112900000611
因此，三角形元上的三个点的应变能量梯度为：
Figure BDA00012676112900000612
Figure BDA00012676112900000613
(1-3)通过构建内部四面体碰撞网格进行球囊的膨胀与紧缩。
首先，求球囊三角形网格与最近的导丝节并通过该三角形网格与导丝节点形成一个四面 体网格。然后，对于每一个四面体采用气态方程PV＝nRT对四面体进行膨胀与紧缩。其中， P为压强，V为体积，n为空气的量，R为常量，T为温度。因此我们设定在一个四面体i内其压强为
Figure BDA0001267611290000071
其中k为常量，Q为空气的量，V为体积。球囊在血管内膨胀和收缩，我 们设定冠状动脉的平均压强为Pblood＝10kPa。最后，通过给予位置的方法来求解约束Cair_tet＝Pi-Pblood。
(二)球囊与血管碰撞交互模拟
(2-1)构建自适应空间碰撞网格
当球囊导管上的球囊到达病处后，球囊即开始膨胀随后即与发生狭窄的血管壁发生碰撞， 我们根据球囊膨胀并与血管发生碰撞区域的固定的特点提出一种自适应的四面体碰撞网格生 成方法，该方法步骤如下：
(2-1-1)我们根据血管管状的特点，采用截取的方式确定血管与球囊可能的碰撞区 域。(a)当球囊在直管处时。我们根据球囊的两端A和B分别查找求得离其最近的两个 中心线控制点并求这两个中心线控制点半径的和作为圆截面的半径。距离A点最近的两 个最近的中心线控制点为xi-1和xi，对应中心线控制点的半径分别为ri-1和ri，因此端点A 的半径rA＝ri-1+ri，同样对于端点B有rB＝ri+ri+1。再根据端点A和B处导管切线方向即 可确定以A和B为中心的两个圆截面，这样也确定了中间可能碰撞的区域。(b)当球囊 在分叉处时。同(a)中相同首先确定A和B处的两个圆截面，其次将分叉处中心线在另 一个子分叉的节点x'i+1设置为该分叉圆截面的中心点C，该圆截面的半径为rC＝2r'i+1。
(2-1-2)通过获得的圆截面得到其间的三角形血管网格，如图2所示，为了方便我们 以直管的位置为例，分叉位置与直管位置处理方法相同。这里我们首先通过两个圆截面分别求血管与其交叉所形成的两个由三角形面片组成的圆环。其次对于其中一个圆环上的三角形面片，通过点A和B形成的矢量来确定三角形上的点到截面的有符号距离，这 样通过该点的有符号的距离即可得到该点是在两个截面的内部还是在外部。最后，我们 即可从上一步获得的截面内的点作为起始点通过网格点与点之间的连接性进行传播迭代 来求得整个内部网格。
(2-1-3)通过步骤(2)得到的血管网格和已知的球囊网格生成其间的四面体碰撞网 格。这里直接采用Si H一文中的四面体网格生成工具TetGen对两种网格生成其间的四面体 碰撞网格。
(2-2)四面体碰撞网格碰撞处理方法
通过(2-1)中初始化生成四面体碰撞网格，我们采用如下三个步骤处理球囊和血管网格 间的碰撞：
(2-2-1)四面体碰撞网格约束：通过对球囊和血管网格间的每一个四面体碰撞网格 添加单侧体积约束来处理其间的相互碰撞。该碰撞检测与处理采用基于位置的动力学约 束的方法，为了防止这些网格出现交叉的情况，对每一个网格添加一个单侧的体积约束， 其约束方程如下：
C(x)＝det[p2-p1 p3-p1 p4-p1]≥0
其中p1,p2,p3和p4分别是四面体碰撞网格的四个顶点。随着球囊与血管网格相互逐渐 靠近，对应四面体碰撞网格的体积逐渐减小，当四面体碰撞网格的体积为负时(C(x)＜0)， 则判定四面球囊和血管网格发生了碰撞。
(2-2-2)四面体碰撞网格优化：由于在物体进行较大的旋转或位移时可能会产生四 面体碰撞网格与实际球囊和血管网格间位置不符而导致四面体碰撞网格体积小于0，这样 会导致碰撞检测发生错误使得物体无法正常的旋转或是位移。因此采用四面体网格优化 方法来优化四面体碰撞网格的结构和质量，优化方法采用Shewchuk J R一文中的边和面增 删的方法来优化四面体，其中在四面体网格的质量的计算方法为：
Figure BDA0001267611290000081
Figure BDA0001267611290000082
其中，lrms是四面体的六条边长度的均方根，V为四面体的体积。
(2-2-3)接触摩擦处理：通过四面体碰撞网格的约束碰撞后，需要施加摩擦力给球囊与血管壁。我们处理步骤如下：
(a')采用步骤(2)中的四面体碰撞网格的体积正负作为是否添加摩擦力的触 发条件，其中四面体碰撞网格的四个点中一个点为球囊顶点，其他三个顶 点为血管网格顶点所组成的三角形面片。当这种四面体碰撞体积为负时， 即发生碰撞时我们计算球囊上的点相对于对应三角形面片上的切向位移 x⊥。
(b')通过步骤(a')中得到的切向位移x⊥，则摩擦力对球囊顶点产生的静摩擦 和动摩擦位移修正分别为：
Figure BDA0001267611290000083
其中，
Figure BDA0001267611290000084
为四面体碰撞网格中一个球囊顶点质量的倒数，
Figure BDA0001267611290000091
为四面体碰撞中剩余三个血管顶点质量的倒数和，μs为静摩擦系数，μk为动摩擦系数，d为设定的距离阈值。
(c')摩擦力作用在球囊上，同样其反作用力作用在血管三个顶点上修正的位移分别有：
Figure BDA0001267611290000092
Figure BDA0001267611290000093
Figure BDA0001267611290000094
(2-2-4)接触粘滞处理：由于血管表面和球囊表面都是湿的表面，在球囊碰撞并与血管表面接触时由于其间潮湿的表面会产生较大的粘滞力，因此其间的粘滞力是不可忽略的。在碰撞检测过程中，四面体网格的优化使得球囊和血管网格的四面体始终维持较 好的结构，即球囊网格上的顶点和血管三角形网格片面上的顶点处于面对面的状态，因 此我们对摩擦力中所述的四面体碰撞网格添加粘滞力约束。我们对四面体碰撞网格的体 积分为三类进行处理：
(a”)当四面体碰撞网格的体积Vair≤Vth1时，则触发该四面体碰撞网格的粘滞项，但 此时并不产生粘滞力。
(b”)当四面体碰撞网格的体积满足Vth1＜Vair＜Vth2，并且已经触发四面体的粘滞项 后，我们采用体积约束Cair＝Vair≤Vth1来表示粘滞力，其中采用三角形的正反 法线方向来表示各部分的粘滞力的方向。
(c”)当四面体碰撞网格的体积Vair≥Vth2时,则取消四面体碰撞网格的粘滞项。
本发明与现有技术相比较，具有如下显而易见的突出实质性特点和显著性技术进步：
(1)提出了基于位置的方法的气囊物理模型。本发明首次引入连续材料介质对球囊成 形术中的球囊模型进行物理建模，与现有的气囊物理模型相比，本发明提出的气 囊物理模型不仅考虑到了真实气囊的材质的力学性质，在平面网格上加入了气囊 厚度属性。
(2)针对复杂的血管系统和网格数量大难以实时完成球囊与血管碰撞检测的问题，本 发明考虑碰撞球囊的局部性，提出一个碰撞网格预处理方法自适应生成碰撞网格来加速球囊和血管的碰撞检测，由此来实时和稳定的完成球囊和血管系统的碰撞 交互。
附图说明
图1是一种球囊血管成形手术过程实时模拟方法程序框图。
图2是球囊在直管处碰撞区域示例图。
图3是球囊在分叉处碰撞区域示例图。
图4是碰撞区域内部血管碰撞三角形网格图。
图5是球囊血管成形术模拟过程图。
图6是球囊血管成形术模拟过程摩擦位移增量随时间变化图。
图7是球囊血管成形术模拟过程粘滞位移增量随时间变化图。
具体实施方式
本发明的优选实施例结合附图说明如下：
实施例一：
参见图1，一种球囊血管成形手术过程实时模拟方法，包括模拟初始化和模拟循环两个 步骤，具体如下：
(1)模拟初始化：初始化球囊的顶点的位置，速度，质量等信息，以及球囊三角形网格 的材料结构矩阵；初始化构建血管与球囊间的四面体碰撞网格。
(2)模拟循环：首先通过显示时间积分预求顶点新的速度以及位置；其次迭代多次来求 解球囊三角形网格的拉格朗日乘子和能量梯度，并更新每个顶点的位置；通过约束四面体碰 撞网格体积大于0来完成球囊和血管网格间的碰撞检测，并对四面体碰撞网格进行优化；最 后求解球囊与血管间的摩擦力和粘滞力，并减弱顶点的速度使系统稳定。
实施例二：
参见图1～4，本实施例与实施例一基本相同，特别之处如下：
(1)模拟初始化：
(1-1)初始化球囊所有顶点的位置、速度、质量以及相关初始化约束条件
(1-2)初始化构建自适应空间四面体碰撞网格，该自适应的空间四面体碰撞网格生 成方法步骤如下：
(1-2-1)确定血管与球囊可能的碰撞区域。
(a)当球囊在直血管处时，参见图2。根据球囊的两端A和B分别查找求得离 其最近的两个中心线控制点并求这两个中心线控制点半径的和作为圆截 面的半径；距离A点最近的两个最近的中心线控制点为xi-1和xi，对应中 心线控制点的半径分别为ri-1和ri，因此端点A的半径rA＝ri-1+ri，同样对 于端点B有rB＝ri+ri+1；再根据端点A和B处导管切线方向即可确定以A 和B为中心的两个圆截面，这样也确定了中间可能碰撞的区域；
(b)当球囊在分叉血管处时，参见图3。同(a)中相同首先确定A和B处的 两个圆截面，其次将分叉处中心线在另一个子分叉的节点x'i+1设置为该分 叉圆截面的中心点C，该圆截面的半径为rC＝2r'i+1；
(1-2-2)参见图4，通过获得的圆截面得到其间的三角形血管网格。
(a)通过圆截面分别求血管与其交叉所形成的由三角形面片组成的圆环，
(b)确定三角形面片所组成的圆环上的顶点是在圆截面的内部还是在外部，
(c)从上一步获得的圆截面内的点作为起始点通过网格点与点之间的连接性 进行传播迭代来求得整个内部网格；
(1-2-3)通过步骤(1-2-2)得到的血管网格和球囊网格生成其间的四面体碰撞网格；
(2)模拟循环：
(2-1)球囊模拟
(1-2-1)通过时间积分预求每一个顶点的新速度
Figure BDA0001267611290000111
和位置 xt＝xt-1+Δtvt，其中vt和vt+1分别表示当前模拟时刻和上一模拟时刻顶点的速度，Δt为模拟 的时间步长，M为质量对角矩阵，
Figure BDA0001267611290000112
为当前时刻的外力，xt和xt-1分别表示当前模拟时刻和 上一模拟时刻顶点的位置；
(1-2-2)约束函数C(x)来约束求得每一个顶点的修正位移
Figure BDA0001267611290000113
使其修正后的顶点位置xt＝xt+Δx满足约束函数，其中
Figure BDA0001267611290000114
为每个顶点质量的倒数，
Figure BDA0001267611290000115
为拉格朗日乘子，
Figure BDA0001267611290000116
为约束函数的梯度。其后通过并行的jacobi 迭代约束其收敛稳定；
(2-2)通过当前时刻和上一时刻的位置更新每个顶点的速度为
Figure BDA0001267611290000117
(2-3)四面体碰撞网格碰撞处理
(2-3-1)四面体碰撞网格约束：对球囊和血管网格间的每一个四面体碰撞网格添加单侧体积约束来处理其间的相互碰撞，该碰撞检测与处理采用基于位置的动力学约束的方法，为了防止这些网格出现交叉的情况，对每一个网格添加一个单侧的体积约束， 其约束方程如下：
C(x)＝det[p2-p1 p3-p1 p4-p1]≥0
其中p1,p2,p3和p4分别是四面体碰撞网格的四个顶点；
(2-3-2)四面体碰撞网格优化：采用边和面增删的方法来优化四面体碰撞网格的结构和质量，在优化过程中四面体网格的质量的计算方法为：
Figure BDA0001267611290000121
Figure BDA0001267611290000122
其中，lrms是四面体的六条边长度的均方根，V为四面体的体积；
(2-3-3)接触摩擦处理：通过四面体碰撞网格的约束碰撞后，需要施加摩擦力给球囊与血管壁，接触摩擦处理步骤如下：
(a')采用步骤(2-3-2)中的四面体碰撞网格的体积正负作为是否添加摩擦 力的触发条件，其中四面体碰撞网格的四个点中一个点为球囊顶点， 其他三个顶点为血管网格顶点所组成的三角形面片；当这种四面体碰 撞体积为负时，即发生碰撞时，计算球囊上的点相对于对应三角形面 片上的切向位移x⊥；
(b')通过步骤(a')中得到的切向位移x⊥，则摩擦力对球囊顶点产生的静 摩擦和动摩擦位移修正分别为：
Figure BDA0001267611290000123
其中，
Figure BDA0001267611290000124
为四面体碰撞网格中一个球囊顶点质量的倒数，
Figure BDA0001267611290000125
为四面体碰撞中剩余三个血管顶点质量的倒数和， μs＝0.5为静摩擦系数，μk＝0.3为动摩擦系数，d＝0.5为设定的距离阈值。
(c')球囊上反作用力作用在血管三个顶点上修正的位移分别有：
Figure BDA0001267611290000131
Figure BDA0001267611290000132
Figure BDA0001267611290000133
(2-2-4)接触粘滞处理：通过对四面体碰撞网格的体积分为三类分别进行粘滞力处理：
(a”)当四面体碰撞网格的体积Vair≤Vth1时，则触发该四面体碰撞网格的粘 滞项，但此时并不产生粘滞力；
(b”)当四面体碰撞网格的体积满足Vth1＜Vair＜Vth2，并且已经触发四面体的 粘滞项后，采用体积约束Cair＝Vair≤Vth1来表示粘滞力，其中采用三角 形的正反法线方向来表示各部分的粘滞力的方向；
(c”)当四面体碰撞网格的体积Vair≥Vth2时,则取消四面体碰撞网格的粘滞 项，其中Vth1＝0.0001，Vth2＝0.0003。
(2-4)减弱每个顶点的速度使系统稳定。
实施例三：
参见图2～4，为实施例一所述模拟初始化步骤构建血管与球囊间的四面体碰撞网格示意 图。
(1)采用截取的方式确定血管与球囊可能的碰撞区域。(a)当球囊在直管处时,参见图2。 根据球囊的两端A和B分别查找求得离其最近的两个中心线控制点并求这两个中心线控制点 半径的和作为圆截面的半径。距离A点最近的两个最近的中心线控制点为xi-1和xi，对应中心 线控制点的半径分别为ri-1和ri，因此端点A的半径rA＝ri-1+ri，同样对于端点B有rB＝ri+ri+1。 再根据端点A和B处导管切线方向即可确定以A和B为中心的两个圆截面，这样也确定了中 间可能碰撞的区域。(b)当球囊在分叉处时,参见图3。同(a)中相同首先确定A和B处的 两个圆截面，其次将分叉处中心线在另一个子分叉的节点x'i+1设置为该分叉圆截面的中心点 C，该圆截面的半径为rC＝2r'i+1。
(2)参见图4通过获得的圆截面得到其间的血管碰撞网格，为了方便我们以直管的位置 为例，分叉位置与直管位置处理方法相同。这里我们首先通过两个圆截面分别求血管与其交 叉所形成的两个由三角形面片组成的圆环。其次对于其中一个圆环上的三角形面片，通过点A 和B形成的矢量来确定三角形上的点到截面的有符号距离，这样通过该点的有符号的距离即 可得到该点是在两个截面的内部还是在外部。最后，我们即可从上一步获得的截面内的点作 为起始点通过网格点与点之间的连接性进行传播迭代来求得整个内部网格。
实施例四：
实施例一所述模拟循环中三角形网格球囊的材料模型采用Mooney-Rivlin模型来模拟气 囊的这种聚氨酯橡胶材料，该材料模型的参数系数参见表1。
表1Mooney-Rivlin模型材料的参数
Figure BDA0001267611290000141
实施例五：
参见图5为球囊血管成形手术过程实时模拟，所述模拟循环对球囊进行充气，球囊接触 血管并挤压血管，当球囊完全膨胀即完成该血管成形术，最后收缩并收回球囊。
实施例六：
参见图6～7分别为实施例一所述模拟循环中在实施例四求解球囊与血管间的摩擦和粘滞 的位移增量。图6为模拟产生的摩擦位置增量变化情况，其中摩擦项中μs＝0.5为静摩擦系数， μk＝0.3为动摩擦系数，d＝0.5为设定的距离阈值；图7为模拟产生的粘滞位置增量变化情 况，粘滞项中Vth1＝0.0001，Vth2＝0.0003。

Claims (6)

1.一种球囊血管成形手术过程实时模拟方法，其特征在于：包括球囊模拟和球囊与血管的碰撞交互模拟两部分，具体如下：
(一)球囊模拟：采用基于位置的方法并结合连续介质材料模型对球囊材质进行模拟；
(二)球囊与血管的碰撞交互模拟：通过构建自适应空间四面体碰撞网格来加速球囊和血管之间的碰撞检测，并通过四面体碰撞网格添加摩擦力和粘滞力；
所述步骤(二)中球囊与血管的碰撞交互模拟的具体步骤如下：
(2-1)初始化构建自适应空间四面体碰撞网格，该自适应的空间四面体碰撞网格生成方法步骤如下：
(2-1-1)确定血管与球囊可能的碰撞区域：
(a)当球囊在直血管处时：根据球囊的两端A和B分别查找求得离其最近的两个中心线控制点并求这两个中心线控制点半径的和作为圆截面的半径；距离A点最近的两个最近的中心线控制点为xi-1和xi，对应中心线控制点的半径分别为ri-1和ri，因此端点A的半径rA＝ri-1+ri，同样对于端点B有rB＝ri+ri+1；再根据端点A和B处导管切线方向即可确定以A和B为中心的两个圆截面，这样也确定了中间可能碰撞的区域；
(b)当球囊在分叉血管处时：同(a)中相同首先确定A和B处的两个圆截面，其次将分叉处中心线在另一个子分叉的节点x'i+1设置为该分叉圆截面的中心点C，该圆截面的半径为rC＝2r'i+1；
(2-1-2)通过获得的圆截面得到其间的三角形血管网格：
(a)通过圆截面分别求血管与其交叉所形成的由三角形面片组成两个由三角形面片组成的圆环；
(b)对于其中一个圆环上的三角形面片，通过点A和B形成的矢量来确定三角形上的点到截面的有符号距离，确定三角形面片所组成的圆环上的顶点是在圆截面的内部还是在外部，
(c)从上一步获得的圆截面内的点作为起始点通过网格点与点之间的连接性进行传播迭代来求得整个内部网格；
(2-1-3)通过步骤(2-1-2)得到的血管网格和球囊网格生成其间的四面体碰撞网格；
(2-2)四面体碰撞网格碰撞处理，对(2-1)中初始化生成的四面体碰撞网格，球囊和血管网格间的碰撞处理步骤如下：
(2-2-1)四面体碰撞网格约束：对球囊和血管网格间的每一个四面体碰撞网格添加单侧体积约束来处理其间的相互碰撞，该碰撞检测与处理采用基于位置的动力学约束的方法，为了防止这些网格出现交叉的情况，对每一个网格添加一个单侧的体积约束，其约束方程如下：
C(x)＝det[p2-p1 p3-p1 p4-p1]≥0
其中p1,p2,p3和p4分别是四面体碰撞网格的四个顶点；
(2-2-2)四面体碰撞网格优化：采用边和面增删的方法来优化四面体碰撞网格的结构和质量，在优化过程中四面体网格的质量的计算方法为：
Figure FDA0003006375700000021
Figure FDA0003006375700000022
其中，lrms是四面体的六条边长度的均方根，V为四面体的体积；
(2-2-3)接触摩擦处理：通过四面体碰撞网格的约束碰撞后，需要施加摩擦力给球囊与血管壁，接触摩擦处理步骤如下：
(a')采用步骤(2-2-2)中的四面体碰撞网格的体积正负作为是否添加摩擦力的触发条件，其中四面体碰撞网格的四个点中一个点为球囊顶点，其他三个顶点为血管网格顶点所组成的三角形面片；当这种四面体碰撞体积为负时，即发生碰撞时，计算球囊上的点相对于对应三角形面片上的切向位移x⊥；
(b')通过步骤(a')中得到的切向位移x⊥，则摩擦力对球囊顶点产生的静摩擦和动摩擦位移修正分别为：
Figure FDA0003006375700000023
其中，
Figure FDA0003006375700000024
为四面体碰撞网格中一个球囊顶点质量的倒数，
Figure FDA0003006375700000025
为四面体碰撞中剩余三个血管顶点质量的倒数和，μs为静摩擦系数，μk为动摩擦系数，d为设定的距离阈值；
(c')球囊上反作用力作用在血管三个顶点上修正的位移分别有：
Figure FDA0003006375700000026
Figure FDA0003006375700000031
Figure FDA0003006375700000032
(2-2-4)接触粘滞处理：在碰撞检测过程中，四面体网格的优化使得球囊和血管网格的四面体始终维持较好的结构，即球囊网格上的顶点和血管三角形网格片面上的顶点处于面对面的状态；对摩擦力中所述的四面体碰撞网格添加粘滞力约束，通过对四面体碰撞网格的体积分为三类分别进行粘滞力处理：
(a”)当四面体碰撞网格的体积Vair≤Vth1时，则触发该四面体碰撞网格的粘滞项，但此时并不产生粘滞力；
(b”)当四面体碰撞网格的体积满足Vth1＜Vair＜Vth2，并且已经触发四面体的粘滞项后，采用体积约束Cair＝Vair≤Vth1来表示粘滞力，其中采用三角形的正反法线方向来表示各部分的粘滞力的方向；
(c”)当四面体碰撞网格的体积Vair≥Vth2时,则取消四面体碰撞网格的粘滞项，其中Vth1＝0.0001，Vth2＝0.0003；
(2-3)减弱每个顶点的速度使系统稳定。
2.根据权利要求1所述的球囊血管成形手术过程实时模拟方法，其特征在于：所述步骤(一)中基于位置的方法为通过约束函数来约束顶点的位置，并通过并行的雅克比迭代法来收敛其约束的顶点。
3.根据权利要求1所述的球囊血管成形手术过程实时模拟方法，其特征在于：所述步骤(一)中连续介质材料模型采用三角形网格对球囊进行离散求解；在变形状态下，通过三角形元的3×2的变形梯度F求解2×2的右柯西格林张量C＝FTF，其次通过应变能密度函数求得三角形元的变形能量E＝∫ΩΨ(C)dx，最后求解应变能量梯度
Figure FDA0003006375700000033
求得三角形每个顶点的修正位移Δx。
4.根据权利要求1所述的球囊血管成形手术过程实时模拟方法，其特征在于：所述步骤(一)中球囊模拟的具体步骤如下：
(1-1)模拟初始化：初始化球囊所有顶点的位置、速度、质量以及相关初始化约束条件；
(1-2)模拟循环：循环模拟初始化设定后的所有顶点，其中模拟循环分为如下若干步骤：
(1-2-1)通过时间积分预求每一个顶点的新速度
Figure FDA0003006375700000041
和位置xt＝xt-1+Δtvt，其中vt和vt+1分别表示当前模拟时刻和上一模拟时刻顶点的速度，Δt为模拟的时间步长，M为质量对角矩阵，
Figure FDA0003006375700000042
为当前时刻的外力，xt和xt-1分别表示当前模拟时刻和上一模拟时刻顶点的位置；
(1-2-2)约束函数C(x)来约束求得每一个顶点的修正位移
Figure FDA0003006375700000043
使其修正后的顶点位置xt＝xt+Δx满足约束函数，其中
Figure FDA0003006375700000044
为每个顶点质量的倒数，
Figure FDA0003006375700000045
为拉格朗日乘子，
Figure FDA0003006375700000046
为约束函数的梯度，其后通过并行的jacobi迭代约束其收敛稳定；
(1-2-3)通过当前时刻和上一时刻的位置更新每个顶点的速度为
Figure FDA0003006375700000047
(1-2-4)碰撞检测与响应。
5.根据权利要求1所述的球囊血管成形手术过程实时模拟方法，其特征在于，所述步骤(二)中自适应空间四面体碰撞网格根据球囊和与球囊可能发生碰撞的局部血管网格来得到，在仿真模拟过程中并通过四面体网格优化对空间四面体碰撞网格进行四面体质量优化。
6.根据权利要求1或5所述的球囊血管成形手术过程实时模拟方法，其特征在于，所述步骤(二)中自适应空间四面体碰撞网格的碰撞检测，为通过约束血管网格和球囊网格间的空间四面体体积大于0。
CN201710234745.8A 2017-04-12 2017-04-12 球囊血管成形手术过程实时模拟方法 Active CN107411819B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710234745.8A CN107411819B (zh) 2017-04-12 2017-04-12 球囊血管成形手术过程实时模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710234745.8A CN107411819B (zh) 2017-04-12 2017-04-12 球囊血管成形手术过程实时模拟方法

Publications (2)

Publication Number Publication Date
CN107411819A CN107411819A (zh) 2017-12-01
CN107411819B true CN107411819B (zh) 2021-08-10

Family

ID=60423976

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710234745.8A Active CN107411819B (zh) 2017-04-12 2017-04-12 球囊血管成形手术过程实时模拟方法

Country Status (1)

Country Link
CN (1) CN107411819B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6113395A (en) * 1998-08-18 2000-09-05 Hon; David C. Selectable instruments with homing devices for haptic virtual reality medical simulation
US7970719B2 (en) * 2008-06-06 2011-06-28 Siemens Aktiengesellschaft Method and simulation device for structurally individualized simulation of the introduction of a wall support element into a section of a tubular structure
WO2012063266A2 (en) * 2010-11-10 2012-05-18 Perfint Healthcare Pvt. Ltd Systems and methods for planning image-guided interventional procedures
CN104081401A (zh) * 2011-11-10 2014-10-01 西门子公司 用于对冠脉循环进行多尺度的解剖学和功能建模的方法和系统
CN106067269A (zh) * 2016-05-13 2016-11-02 中国科学院自动化研究所 虚拟心血管介入手术培训系统中反馈力的确定方法及系统

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002070980A1 (en) * 2001-03-06 2002-09-12 The Johns Hopkins University School Of Medicine Simulation system for image-guided medical procedures
WO2002097735A1 (en) * 2001-05-31 2002-12-05 Kent Ridge Digital Labs System and method of anatomical modeling
EP1805744A2 (en) * 2004-08-10 2007-07-11 The General Hospital Corporation Methods and apparatus for simulation of endovascular and endoluminal procedures

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6113395A (en) * 1998-08-18 2000-09-05 Hon; David C. Selectable instruments with homing devices for haptic virtual reality medical simulation
US7970719B2 (en) * 2008-06-06 2011-06-28 Siemens Aktiengesellschaft Method and simulation device for structurally individualized simulation of the introduction of a wall support element into a section of a tubular structure
WO2012063266A2 (en) * 2010-11-10 2012-05-18 Perfint Healthcare Pvt. Ltd Systems and methods for planning image-guided interventional procedures
CN104081401A (zh) * 2011-11-10 2014-10-01 西门子公司 用于对冠脉循环进行多尺度的解剖学和功能建模的方法和系统
CN106067269A (zh) * 2016-05-13 2016-11-02 中国科学院自动化研究所 虚拟心血管介入手术培训系统中反馈力的确定方法及系统

Also Published As

Publication number Publication date
CN107411819A (zh) 2017-12-01

Similar Documents

Publication Publication Date Title
Wolters et al. A patient-specific computational model of fluid–structure interaction in abdominal aortic aneurysms
Cani-Gascuel et al. Animation of deformable models using implicit surfaces
Bridson et al. Robust treatment of collisions, contact and friction for cloth animation
Tezduyar et al. Modelling of fluid–structure interactions with the space–time finite elements: solution techniques
Borazjani Fluid–structure interaction, immersed boundary-finite element method simulations of bio-prosthetic heart valves
Desbrun et al. Animating soft substances with implicit surfaces
Teran et al. Adaptive physics based tetrahedral mesh generation using level sets
Tang et al. Simulating cyclic artery compression using a 3D unsteady model with fluid–structure interactions
CN107274414A (zh) 基于改进局部信息的cv模型的图像分割方法
Kim et al. Anisotropic elasticity for inversion-safety and element rehabilitation
Vigmostad et al. Fluid–structure interaction methods in biological flows with special emphasis on heart valve dynamics
Xiong et al. Virtual interventions for image-based blood flow computation
Diaby et al. Buckling and wrinkling of prestressed membranes
Zhong et al. A physically based method for triangulated surface flattening
CN107411819B (zh) 球囊血管成形手术过程实时模拟方法
Vassilev et al. A mass-spring model for real time deformable solids
Jaramillo et al. Deformable part inspection using a spring–mass system
CN101276483A (zh) 一种实现平移敏感的Laplacian网格编辑方法
Thornton et al. Three-dimensional stress analysis of polypropylene leaflets for prosthetic heart valves
Sarayi et al. A Nonlinear Mechanics-based Virtual Coiling Method For Intracranial Aneurysm
CN109003319B (zh) 角色动画中带有动力学约束的元球驱动蒙皮方法
Petit et al. Tracking fractures of deformable objects in real-time with an RGB-D sensor
Gostaf et al. Pressure boundary conditions for blood flows
Egger et al. Fast self-collision detection and simulation of bifurcated stents to treat abdominal aortic aneurysms (AAA)
Koh et al. Hierarchical shape representation using locally adaptive finite elements

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant