虚拟手术中基于粒子的多相耦合方法
【专利摘要】一种虚拟手术中基于粒子的多相耦合方法，对于虚拟手术中涉及的软组织和刚体，分别根据三维组织表面三角网格模型生成粒子化三维体模型，而对于液体则采用SPH方法建立粒子化液体模型；建立软组织、液体和刚体的粒子化力学模型，包括将软组织和液体定义为流变体，并对软组织粒子和液体粒子的任意两个同相或两个不同相粒子之间定义流变粘滞系数；进行基于边界粒子和多相耦合核函数的碰撞检测及响应。本发明提出一种了应用于虚拟手术中的通用的、精确的和高效的基于粒子的多相耦合方法，能在虚拟手术仿真中实现精确、高效的多相耦合仿真。
【专利说明】虚拟手术中基于粒子的多相耦合方法
【技术领域】
[0001]本发明属于计算机仿真及虚拟现实【技术领域】，特别是涉及一种应用于虚拟手术中的通用的、精确的和高效的基于粒子的多相耦合方法。
【背景技术】
[0002]随着计算机仿真计算和虚拟现实技术的快速发展，虚拟手术仿真研究及其相关应用已成为国际上发展迅速的一个领域，以Simbionix为代表的虚拟手术仿真器已广泛应用于许多种手术的外科医生训练和评价中。真实手术中往往涉及不同性质物体的多相耦合(软组织、液体和刚体之间的耦合)问题，如在腹腔手术中手术器械与体液、软组织的接触耦合，体液与组织、肿瘤的接触耦合，这些耦合现象的仿真对精确性和效率要求较高。
[0003]近年来，研究人员提出了很多基于物理的计算方法用于多相耦合仿真。基于物理的数值计算方法一般分为两类，一类是基于网格的方法或Eluer方法，另一类是基于粒子的方法或Lagrangian方法。基于网格的方法,如有限元法(Finite Element Method,简称FEM方法)等使用广泛，但在处理自由表面问题、边界运动问题、表面运动问题和大形变及碰撞传播问题时并不适用，并且对于复杂的几何模型，要生成高质量的网格是一个困难耗时的过程，无网格法则能有效地解决这些问题[I]。光滑质点流体动力学(Smoothed ParticleHydrodynamics,简称SPH方法)是近年来无网格法中使用最为广泛的一种。SPH方法本身的特性使其非常适合对高度可变形物体进行仿真[2]。
[0004]M.MUller提出了一个基于粒子的方法，借助SPH方法对弹性体、塑料体和融化物体进行建模[I]，能够有效地仿真物体的形变。之后关于物体形变的研究大多基于M.MUller的工作，并且相关改进多集中于计算的简化[3]以及稳定性[4]的提高。一些方法简单地把流体、固体粒子视为整体来计算邻域[4]，但这将导致流体粒子进入固体的不真实现象。一些研究学者采用虚粒子法[5]在固体、流体周围设置虚粒子用于耦合，但这并不能仿真双向耦合过程，且虚粒子的设置与计算过程相当复杂。一种新颖的、通用的可以同时用于单向双向耦合的流固耦合的边界处理方法在文献[6]中被提出。该模型将固体建模成刚体，因为是刚体，所以该模型只需要与固体最外层粒子进行计算，且不需要粒子法向量，同时计算的边界粒子不需要均匀采样，可以稀疏可以浓密，流固边界粒子的受力是对称的。但是该模型只适用于刚体，不能用于流体与实体的可变形的固体之间的耦合计算。
[0005]虽然近年来研究人员在多相耦合的领域取得了一定的成绩，但是要在虚拟手术仿真中实现精确、高效率的多相耦合仿真，还需要解决现有方法所面临的如下一些问题:
[0006](I)随着虚拟手术仿真场景的复杂度的提高，现有方法处理复杂模型场景的复杂度越来越高，越来越不能适应复杂场景的仿真，因此，亟需一种较高真实性的、适用于任意复杂模型、多种复杂场景仿真的通用多相耦合模型。且真实人体软组织形状大小复杂多变，现有方法无法对任意复杂的人体软组织模型精确地生成内部体数据，为了适应复杂虚拟手术仿真场景及模型的满足物理真实性，需要对任意复杂的人体软组织模型精确地生成内部体数据。[0007](2)现有大部分基于物理的耦合方法通常采用类似于SPH方法的无网格法来对流体建模，采用线弹性模型对软组织建模，这些方法都缺乏对流体和软组织性质的考虑。例如在真实世界中，几乎所有的生物软组织都是由固态和液态两种成分共同构成，生物软组织是最典型的流变体，体液具有松弛行为而软组织具有蠕变行为，因此从虚拟手术仿真的真实性的角度考量，在建立流体和固体模型的时候需要引入流体和固体的性质。
[0008](3)现有方法在处理多物体、多相的耦合过程(包括碰撞检测和响应)比较复杂，在碰撞检测和响应的过程中耗费较多资源，且通常选取与运动粒子最近的粒子作为碰撞点，缺乏精确性，容易产生穿透现象，因此需要一个精确高效的碰撞检测和响应算法来实现多物体、多相的耦合过程。
[0009]相关文献:
[0010][I] Mul I er，Μ.，Keiserj R.，Nealenj A.，Pauly, Μ.，Gross, Μ.，&Alexa, Μ.(2004，August).Point based animation of elastic, plastic and melting objects.1n Proceedings of the2004ACM SIGGRAPH/Eurographics symposium on Computeranimation(pp.141-151).Eurographics Association.[0011][2] Liu，G.G.R.，&Liu，Μ.B.(2003).Smoothed particle hydrodynamics: ameshfree particle method.World Scientific Publishing Company.[0012][3] Gerszewskij D.，Bhattacharyaj H.，&Bargtei I, A.W.(2009，August).A point-based method for animating elastoplastic solids.1n Proceedings ofthe2009ACM SIGGRAPH/Eurographics Symposium on Computer Animation(pp.133-138).ACM.[0013][4] Solenthaler, B.，SoMifli5 J_，&Pajarola, R.(2007) _ A unified particlemodel for fluid - sol id interactions.Computer Animation and VirtualWorlds, 18(1),69-82.[0014][5] Schechter，H.，&Bridson，R.(2012).Ghost sph for animating water.ACMTransactions on Graphics(TOG)，31(4)，61.[0015][6] Akincij N.，Ihmsenj M.，Akincij G.，Solenthalerj B.，&Teschner, M.(2012).Versatile rigid-fluid coupling for incompressible SPH.ACM Transactions onGraphics (TOG) ,31 (4) ,62.
【发明内容】
[0016]为了弥补现有耦合方法在虚拟手术仿真中的不足，并在虚拟手术仿真中实现精确、高效的多相耦合仿真，本发明提出一种虚拟手术中基于粒子的多相耦合方法。
[0017]本发明的技术方案为一种虚拟手术中基于粒子的多相耦合方法，包括以下步骤:
[0018]步骤1，对于虚拟手术中涉及的软组织和刚体，分别根据三维组织表面三角网格模型生成粒子化三维组织体模型，包括以下子步骤，
[0019]①对三维组织表面三角网格模型建立最小长方体包围盒；
[0020]②将步骤①中最小长方体包围盒划分为若干大小相同的立方体的体素，在边界不足处采用立方体的体素补充，并以最小长方体包围盒左下角顶点为起点对所有立方体的体素编号，得到各体素的体素序号；[0021]③在三维组织表面三角网格模型的每一个三角面片内部均匀填充顶点，得到填充顶点后的三维组织表面三角网格模型，并为所有顶点编号得到顶点序号；
[0022]④对步骤③中填充顶点后的三维组织表面三角网格模型，确定每一个顶点所位于的体素，并建立体素-顶点关系表，该关系表的内容是步骤②所得每个体素的体素序号及该体素所包含顶点的顶点序号；
[0023]⑤根据步骤④中的体素-顶点关系表，将所有包含顶点的体素定义为边界体素；
[0024]⑥根据步骤⑤中所确定的边界体素，对最小长方体包围盒沿着O-XYZ坐标系的任一坐标轴方向逐层确定边界体素所包围的内部体素，得到所有的内部体素；
[0025]⑦将步骤③中填充顶点后的三维组织表面三角网格模型上每一个顶点定义为一个粒子，得到边界粒子；并将步骤⑥中得到的每一个内部体素取中心位置定义为一个粒子，得到内部粒子；边界粒子和内部粒子构成粒子化三维组织体模型；
[0026]步骤2，建立软组织、液体和刚体的粒子化力学模型，将软组织和液体定义为流变体，并对软组织粒子和液体粒子的任意两个粒子之间定义流变粘滞系数；
[0027]所述刚体的粒子化力学模型，只采用步骤I所得刚体的粒子化三维组织体模型的边界粒子表示刚体，刚体作为一个整体移动和转动，并计算刚体所受的合外力以及转矩；
[0028]所述软组织的粒子化力学模型，采用步骤I所得软组织的粒子化三维软组织体模型表示软组织，根据弹性力学基本方程和SPH方法求解得到各个粒子的物理量，所述物理量包括弹力、粘滞力和重力，计算粘滞力时采用软组织粒子之间的流变粘滞系数；
[0029]所述液体的粒子化力学模型，采用SPH方法建立，根据描述液体运动的Navier-Stokes方程组和SPH方法求解得到各个粒子的物理量，所述物理量包括粘滞力、压力和重力，计算粘滞力时采用液体粒子之间的流变粘滞系数；
[0030]步骤3，进行基于边界粒子和多相耦合核函数的碰撞检测及响应，对于任意两相的物质A和B，A为软组织、液体或刚体，B为软组织或刚体，且A和B不相同，
[0031]碰撞检测过程如下，
[0032]①初始化A型粒子半径和B型粒子半径，
[0033]②对每一个A型粒子，以A型粒子当前位置为圆心，一个时间步长内移动的距离为半径作球形搜索区域，检测出此球形搜索区域内所有的B型边界粒子；计算得到一个时间步长内A型粒子由当前位置O移动到的位置O'，以O为起点指向O'作一条射线；
[0034]③将步骤②中所得的每个B型边界粒子投影至步骤②中的射线所在的直线上，并求解每个B型边界粒子到相应投影点的投影距离；
[0035]④在步骤③中检测出所有满足投影距离小于等于A型粒子半径与B型边界粒子半径之和的B型边界粒；
[0036]⑤在步骤④检测出的所有的B型边界粒子中满足条件的B型边界粒子为碰撞点，所述条件为，该B型边界粒子的投影点位于步骤②中A型粒子当前运动方向的前方，并其投影点距离步骤②中A型粒子圆心O最近；
[0037]碰撞响应过程如下，
[0038]①根据碰撞检测过程，判断液体粒子、软组织粒子和刚体粒子中任意两相粒子是否发生碰撞，设该两相粒子为A型粒子和B型粒子；
[0039]②若A型粒子和B型粒子发生碰撞，则定义A型粒子和B型粒子作用的耦合核函数；
[0040]③根据步骤②中的耦合核函数计算A型粒子和B型粒子之间的耦合作用力，并根据冲量定理求解下一时刻A型粒子和B型粒子的速度大小和方向；
[0041 ] ④根据步骤③的计算结果，代入A型粒子和B型粒子相应的粒子化力学模型，计算出A型粒子和B型粒子新的位移、速度，而后继续进行碰撞检测过程。
[0042]而且，定义粒子化力学模型中两个粒子i和j之间的流变粘滞系数?其中粒子i和j为软组织粒子或液体粒子，
【权利要求】
1.一种虚拟手术中基于粒子的多相耦合方法，其特征在于，包括以下步骤: 步骤I，对于虚拟手术中涉及的软组织和刚体，分别根据三维组织表面三角网格模型生成粒子化三维组织体模型，包括以下子步骤， ①对三维组织表面三角网格模型建立最小长方体包围盒； ②将步骤①中最小长方体包围盒划分为若干大小相同的立方体的体素，在边界不足处采用立方体的体素补充，并以最小长方体包围盒左下角顶点为起点对所有立方体的体素编号，得到各体素的体素序号； ③在三维组织表面三角网格模型的每一个三角面片内部均匀填充顶点，得到填充顶点后的三维组织表面三角网格模型，并为所有顶点编号得到顶点序号； ④对步骤③中填充顶点后的三维组织表面三角网格模型，确定每一个顶点所位于的体素，并建立体素-顶点关系表，该关系表的内容是步骤②所得每个体素的体素序号及该体素所包含顶点的顶点序号； ⑤根据步骤④中的体素-顶点关系表，将所有包含顶点的体素定义为边界体素； ⑥根据步骤⑤中所确定的边界体素，对最小长方体包围盒沿着O-XYZ坐标系的任一坐标轴方向逐层确定边界体素所包围的内部体素，得到所有的内部体素； ⑦将步骤③中填充顶点后的三维组织表面三角网格模型上每一个顶点定义为一个粒子，得到边界粒子；并将步骤⑥中得到的每一个内部体素取中心位置定义为一个粒子，得到内部粒子；边界粒子和内部粒子构成粒子化三维组织体模型」 步骤2，建立软组织、液体和刚体的粒子化力学模型，将软组织和液体定义为流变体，并对软组织粒子和液体粒子的任意两个粒子之间定义流变粘滞系数； 所述刚体的粒子化力学模型，只采用步骤I所得刚体的粒子化三维组织体模型的边界粒子表示刚体，刚体作为一个整体移动和转动，并计算刚体所受的合外力以及转矩； 所述软组织的粒子化力学模型，采用步骤I所得软组织的粒子化三维软组织体模型表示软组织，根据弹性力学基本方程和SPH方法求解得到各个粒子的物理量，所述物理量包括弹力、粘滞力和重力I计算粘滞力时采用软组织粒子之间的流变粘滞系数； 所述液体的粒子化力学模型，采用SPH方法建立,根据描述液体运动的Navier-Stokes方程组和SPH方法求解得到各个粒子的物理量，所述物理量包括粘滞力、压力和和重力，计算粘滞力时采用液体粒子之间的流变粘滞系数； 步骤3，进行基于边界粒子和多相耦合核函数的碰撞检测及响应，对于任意两相的物质A和B，A为软组织、液体或刚体，B为软组织或刚体，且A和B不相同， 碰撞检测过程如下， ①初始化A型粒子半径和B型粒子半径， ②对每一个A型粒子，以A型粒子当前位置为圆心，一个时间步长内移动的距离为半径作球形搜索区域，检测出此球形搜索区域内所有的B型边界粒子；计算得到一个时间步长内A型粒子由当前位置O移动到的位置O'，以O为起点指向O'作一条射线； ③将步骤②中所得的每个B型边界粒子投影至步骤②中的射线所在的直线上，并求解每个B型边界粒子到相应投影点的投影距离； ④在步骤③中检测出所有满足投影距离小于等于A型粒子半径与B型边界粒子半径之和的B型边界粒；⑤在步骤④检测出的所有的B型边界粒子中满足条件的B型边界粒子为碰撞点，所述条件为，该B型边界粒子的投影点位于步骤②中A型粒子当前运动方向的前方，并其投影点距离步骤②中A型粒子圆心O最近； 碰撞响应过程如下， ①根据碰撞检测过程，判断液体粒子、软组织粒子和刚体粒子中任意两相粒子是否发生碰撞，设该两相粒子为A型粒子和B型粒子； ②若A型粒子和B型粒子发生碰撞，则定义A型粒子和B型粒子作用的耦合核函数； ③根据步骤②中的耦合核函数计算A型粒子和B型粒子之间的耦合作用力，并根据冲量定理求解下一时刻A型粒子和B型粒子的速度大小和方向； ④根据步骤③的计算结果，代入A型粒子和B型粒子相应的粒子化力学模型，计算出A型粒子和B型粒子新的位移、速度，而后继续进行碰撞检测过程。
2.根据权利要求1所述虚拟手术中基于粒子的多相耦合方法，其特征在于:定义粒子化力学模型中两个粒子i和j之间的流变粘滞系数其中粒子i和j为软组织粒子或液体粒子，
3.根据权利要求2所述虚拟手术中基于粒子的多相耦合方法，其特征在于:Α型粒子与B型粒子之间作用的耦合核函数如下，
4.根据权利要求3所述虚拟手术中基于粒子的多相耦合方法，其特征在于:A型粒子和B型粒子的耦合作用力计算如下，
【文档编号】G06T17/30GK103559741SQ201310602225
【公开日】2014年2月5日 申请日期:2013年11月25日 优先权日:2013年11月25日
【发明者】袁志勇, 廖祥云, 郭甲翔, 喻思娇 申请人:武汉大学