CN103559741A - 虚拟手术中基于粒子的多相耦合方法 - Google Patents

Abstract

一种虚拟手术中基于粒子的多相耦合方法，对于虚拟手术中涉及的软组织和刚体，分别根据三维组织表面三角网格模型生成粒子化三维体模型，而对于液体则采用SPH方法建立粒子化液体模型；建立软组织、液体和刚体的粒子化力学模型，包括将软组织和液体定义为流变体，并对软组织粒子和液体粒子的任意两个同相或两个不同相粒子之间定义流变粘滞系数；进行基于边界粒子和多相耦合核函数的碰撞检测及响应。本发明提出一种了应用于虚拟手术中的通用的、精确的和高效的基于粒子的多相耦合方法，能在虚拟手术仿真中实现精确、高效的多相耦合仿真。

Description

虚拟手术中基于粒子的多相耦合方法

技术领域

本发明属于计算机仿真及虚拟现实技术领域，特别是涉及一种应用于虚拟手术中的通用的、精确的和高效的基于粒子的多相耦合方法。

背景技术

随着计算机仿真计算和虚拟现实技术的快速发展，虚拟手术仿真研究及其相关应用已成为国际上发展迅速的一个领域，以Simbionix为代表的虚拟手术仿真器已广泛应用于许多种手术的外科医生训练和评价中。真实手术中往往涉及不同性质物体的多相耦合（软组织、液体和刚体之间的耦合）问题，如在腹腔手术中手术器械与体液、软组织的接触耦合，体液与组织、肿瘤的接触耦合，这些耦合现象的仿真对精确性和效率要求较高。

近年来，研究人员提出了很多基于物理的计算方法用于多相耦合仿真。基于物理的数值计算方法一般分为两类，一类是基于网格的方法或Eluer方法，另一类是基于粒子的方法或Lagrangian方法。基于网格的方法，如有限元法（Finite Element Method，简称FEM方法）等使用广泛，但在处理自由表面问题、边界运动问题、表面运动问题和大形变及碰撞传播问题时并不适用，并且对于复杂的几何模型，要生成高质量的网格是一个困难耗时的过程，无网格法则能有效地解决这些问题[1]。光滑质点流体动力学（Smoothed Particle Hydrodynamics，简称SPH方法）是近年来无网格法中使用最为广泛的一种。SPH方法本身的特性使其非常适合对高度可变形物体进行仿真[2]。

M.Müller提出了一个基于粒子的方法，借助SPH方法对弹性体、塑料体和融化物体进行建模[1]，能够有效地仿真物体的形变。之后关于物体形变的研究大多基于M.Müller的工作，并且相关改进多集中于计算的简化[3]以及稳定性[4]的提高。一些方法简单地把流体、固体粒子视为整体来计算邻域[4]，但这将导致流体粒子进入固体的不真实现象。一些研究学者采用虚粒子法[5]在固体、流体周围设置虚粒子用于耦合，但这并不能仿真双向耦合过程，且虚粒子的设置与计算过程相当复杂。一种新颖的、通用的可以同时用于单向双向耦合的流固耦合的边界处理方法在文献[6]中被提出。该模型将固体建模成刚体，因为是刚体，所以该模型只需要与固体最外层粒子进行计算，且不需要粒子法向量，同时计算的边界粒子不需要均匀采样，可以稀疏可以浓密，流固边界粒子的受力是对称的。但是该模型只适用于刚体，不能用于流体与实体的可变形的固体之间的耦合计算。

虽然近年来研究人员在多相耦合的领域取得了一定的成绩，但是要在虚拟手术仿真中实现精确、高效率的多相耦合仿真，还需要解决现有方法所面临的如下一些问题：

（1）随着虚拟手术仿真场景的复杂度的提高，现有方法处理复杂模型场景的复杂度越来越高，越来越不能适应复杂场景的仿真，因此，亟需一种较高真实性的、适用于任意复杂模型、多种复杂场景仿真的通用多相耦合模型。且真实人体软组织形状大小复杂多变，现有方法无法对任意复杂的人体软组织模型精确地生成内部体数据，为了适应复杂虚拟手术仿真场景及模型的满足物理真实性，需要对任意复杂的人体软组织模型精确地生成内部体数据。

（2）现有大部分基于物理的耦合方法通常采用类似于SPH方法的无网格法来对流体建模，采用线弹性模型对软组织建模，这些方法都缺乏对流体和软组织性质的考虑。例如在真实世界中，几乎所有的生物软组织都是由固态和液态两种成分共同构成，生物软组织是最典型的流变体，体液具有松弛行为而软组织具有蠕变行为，因此从虚拟手术仿真的真实性的角度考量，在建立流体和固体模型的时候需要引入流体和固体的性质。

（3）现有方法在处理多物体、多相的耦合过程（包括碰撞检测和响应）比较复杂，在碰撞检测和响应的过程中耗费较多资源，且通常选取与运动粒子最近的粒子作为碰撞点，缺乏精确性，容易产生穿透现象，因此需要一个精确高效的碰撞检测和响应算法来实现多物体、多相的耦合过程。

相关文献：

[1]Müller,M.,Keiser,R.,Nealen,A.,Pauly,M.,Gross,M.,&Alexa,M.(2004,August).Pointbased animation of elastic,plastic and melting objects.In Proceedings of the2004ACMSIGGRAPH/Eurographics symposium on Computer animation(pp.141-151).EurographicsAssociation.

[2]Liu,G.G.R.,&Liu,M.B.(2003).Smoothed particle hydrodynamics:a meshfree particlemethod.World Scientific Publishing Company.

[3]Gerszewski,D.,Bhattacharya,H.,&Bargteil,A.W.(2009,August).A point-based method foranimating elastoplastic solids.In Proceedings of the2009ACM SIGGRAPH/EurographicsSymposium on Computer Animation(pp.133-138).ACM.

[4]Solenthaler,B.,J.,&Pajarola,R.(2007).A unified particle model for fluid–solidinteractions.Computer Animation and Virtual Worlds,18(1),69-82.

[5]Schechter,H.,&Bridson,R.(2012).Ghost sph for animating water.ACM Transactions onGraphics(TOG),31(4),61.

[6]Akinci,N.,Ihmsen,M.,Akinci,G.,Solenthaler,B.,&Teschner,M.(2012).Versatile rigid-fluidcoupling for incompressible SPH.ACM Transactions on Graphics(TOG),31(4),62.

发明内容

为了弥补现有耦合方法在虚拟手术仿真中的不足，并在虚拟手术仿真中实现精确、高效的多相耦合仿真，本发明提出一种虚拟手术中基于粒子的多相耦合方法。

本发明的技术方案为一种虚拟手术中基于粒子的多相耦合方法，包括以下步骤：

步骤1，对于虚拟手术中涉及的软组织和刚体，分别根据三维组织表面三角网格模型生成粒子化三维组织体模型，包括以下子步骤，

①对三维组织表面三角网格模型建立最小长方体包围盒；

②将步骤①中最小长方体包围盒划分为若干大小相同的立方体的体素，在边界不足处采用立方体的体素补充，并以最小长方体包围盒左下角顶点为起点对所有立方体的体素编号，得到各体素的体素序号；

③在三维组织表面三角网格模型的每一个三角面片内部均匀填充顶点，得到填充顶点后的三维组织表面三角网格模型，并为所有顶点编号得到顶点序号；

④对步骤③中填充顶点后的三维组织表面三角网格模型，确定每一个顶点所位于的体素，并建立体素-顶点关系表，该关系表的内容是步骤②所得每个体素的体素序号及该体素所包含顶点的顶点序号；

⑤根据步骤④中的体素-顶点关系表，将所有包含顶点的体素定义为边界体素；

⑥根据步骤⑤中所确定的边界体素，对最小长方体包围盒沿着O-XYZ坐标系的任一坐标轴方向逐层确定边界体素所包围的内部体素，得到所有的内部体素；

⑦将步骤③中填充顶点后的三维组织表面三角网格模型上每一个顶点定义为一个粒子，得到边界粒子；并将步骤⑥中得到的每一个内部体素取中心位置定义为一个粒子，得到内部粒子；边界粒子和内部粒子构成粒子化三维组织体模型；

步骤2，建立软组织、液体和刚体的粒子化力学模型，将软组织和液体定义为流变体，并对软组织粒子和液体粒子的任意两个粒子之间定义流变粘滞系数；

所述刚体的粒子化力学模型，只采用步骤1所得刚体的粒子化三维组织体模型的边界粒子表示刚体，刚体作为一个整体移动和转动，并计算刚体所受的合外力以及转矩；

所述软组织的粒子化力学模型，采用步骤1所得软组织的粒子化三维软组织体模型表示软组织，根据弹性力学基本方程和SPH方法求解得到各个粒子的物理量，所述物理量包括弹力、粘滞力和重力，计算粘滞力时采用软组织粒子之间的流变粘滞系数；

所述液体的粒子化力学模型，采用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型粒子新的位移、速度，而后继续进行碰撞检测过程。

而且，定义粒子化力学模型中两个粒子i和j之间的流变粘滞系数

其中粒子i和j为软组织粒子或液体粒子，

${\mathrm{\μ}}_{j\↔i}=K{\left(\frac{\left|v_{\mathrm{\τ}}\right|}{|x_i-x_j|}\right)}^{\mathrm{\β}}---\left(1\right)$

其中，K、β为与软组织和液体相关的系数，设当两个粒子为液体粒子时，K、β的取值为K1、β1，当两个粒子为软组织粒子时，K、β的取值为K2、β2，当两个粒子一个是液体粒子而一个是软组织粒子时，则K、β分别取（K1+K2）/2、(β1+β2）/2，K1<K2,β1<β2；

v_τ为粒子i相对粒子j运动切向方向的相对速度，x_i和x_j、v_i和v_j分别是粒子i和粒子j的位置和速度，

$\left|v_{\mathrm{\&tau;}}\right|=|v_i-v_j|\sin \mathrm{\&theta;}=\sqrt{1-{\cos }^2\mathrm{\&theta;}}|v_i-v_j|---\left(2\right)$

其中，

$\cos \mathrm{\&theta;}=\frac{(v_i-v_j)(x_j-x_i)}{|v_i-v_j||x_j-x_i|}---\left(3\right)$

角度θ根据上式得到。

而且，A型粒子与B型粒子之间作用的耦合核函数如下，

$W_{A\&LeftRightArrow;B}={\mathrm{\&delta;}}_A^{A\&LeftRightArrow;B}{W}_A+{\mathrm{\&delta;}}_B^{A\&LeftRightArrow;B}{W}_B---\left(13\right)$

其中，

分别表示在B型粒子的作用域内A型粒子与B型粒子数目所占的比例，且n_A,n_B分别为B型粒子作用域内A型粒子和B型粒子的数目，W_A表示A型粒子间作用的核函数，W_B表示B型粒子间作用的核函数。

而且，A型粒子和B型粒子的耦合作用力

计算如下，

$f_{A\&LeftRightArrow;B}=-{\mathrm{\&mu;}}_{A\&LeftRightArrow;B}{m}_A{m}_A\&dtri;W_{A\&LeftRightArrow;B}---\left(14\right)$

其中，m_A，m_B分别表示A型粒子与B型粒子的质量，

为A型粒子与B型粒子之间的流变粘滞系数。

本发明的优点在于，提出一种了应用于虚拟手术中的通用的、精确的和高效的基于粒子的多相耦合方法，该方法能高效精确地处理软组织、液体和刚体粒子之间的碰撞检测及响应，并根据所建立的粒子化力学模型精确地模拟软组织、液体和刚体在多相耦合作用过程中的变化情况，能在虚拟手术仿真中实现精确、高效的多相耦合仿真，能够应用于复杂环境，并且节约了系统资源。

附图说明

图1是本发明实施例的流程图。

图2是本发明实施例的粒子化力学模型中软组织粒子和液体粒子中任意两个粒子之间的流变粘滞系数定义示意图。

图3是本发明实施例的以液体粒子和软组织粒子为例的碰撞检测方法示意图。

图4是本发明实施例的多相耦合核函数的示意图。

图5是本发明实施例的寻找内部体素示意图。

具体实施方式

为了在虚拟手术中实现精确、高效的多相耦合仿真，本发明提出一种基于粒子的多相耦合方法，具体内容包括三个部分：一种对任意三维软组织表面三角网格模型生成粒子化三维软组织体模型的方法、一种基于流变力学的软组织、液体和刚体的粒子化力学模型和一种基于边界粒子和多相耦合核函数的碰撞检测及响应方法。该方法可采用计算机软件技术实现自动运行流程，实施例的处理流程如图1所示，每个步骤的详细说明如下：

步骤1，对于虚拟手术中涉及的软组织和刚体，分别根据三维组织表面三角网格模型生成粒子化三维组织体模型。

在虚拟手术仿真中，通常使用表面三角网格数据进行仿真。本发明采用粒子法来表示软组织，为了满足物理真实性，需将软组织建模成连续介质模型，不仅需要表面粒子，更需要内部粒子以完成各物理量的传递。然而真实人体软组织形状大小复杂多变，为了适应复杂手术仿真场景的需要，需对任意复杂的人体软组织模型精确地生成内部体数据。本发明提出了一种对任意三维软组织表面三角网格模型生成粒子化三维软组织体模型的技术方案，该技术方案所生成的粒子化三维软组织体模型可适用于虚拟手术中任意复杂场景下的多相耦合仿真，对刚体也可以采用同样的方式建模。

实施例根据软组织或刚体的三维组织表面三角网格模型生成粒子化三维组织体模型，具体步骤如下：

1.设三维组织表面三角网格模型置于O-XYZ坐标系中，对三维组织表面三角网格模型建立最小长方体包围盒，该包围盒的一条空间对角线的两个顶点坐标分别是此三维组织表面三角网格模型的所有顶点的x,y,z坐标值的最小和最大值构成，这两个顶点的坐标分别记为(x_min,y_min,z_min)和(x_max,y_max,z_max)。

2.将步骤1中最小长方体包围盒划分为若干大小相同的立方体体素，在边界不足处采用立方体体素补充，并以顶点(x_min,y_min,z_min)为起始点对所有立方体体素编号，具体实

施时，本领域技术人员可自行设定立方体体素的尺寸。

3.在三维组织表面三角网格模型的每一个三角面片内部均匀填充顶点，增加三维组织表面模型顶点密度，得到填充顶点后的三维组织表面网格模型，这样做的目的是增加三维组织表面三角网格模型的顶点密度，使得在后续步骤5中得到的边界体素完全包裹住三维组织表面三角网格。具体填充方式可为，对每一个三角面片而言，在任一条边上作等分点，其等分数量应不小于该三角形最长边除以立方体体素边长的倍数的最小整数。选取三角形任一条边上的所有等分点，作三角形另外两条边的平行线，与另外两条边的相交点恰好是此两边上的等分点，这些平行线在三角形中相交形成的顶点即为填充的顶点，可对填充顶点后的三维组织表面三角网格模型上每一个顶点分配顶点序号。

4.确定步骤3中填充顶点后的三维组织表面三角网格模型上每一个顶点所位于的体素，并建立体素-顶点关系表，该关系表的内容可以是每个体素序号及该体素所包含的三维组织表面三角网格顶点序号。可根据三维组织表面三角网格模型中任一顶点与步骤2中的起始点的位置关系，判断该顶点所在的体素位置。

5.体素可分为外部体素、边界体素和内部体素。三维组织表面三角网格所位于的体素称为边界体素，位于三维组织表面三角网格的内部的体素称为内部体素，除了边界体素和内部体素，剩下的为外部体素。实施例根据步骤4中的体素-顶点关系表，将所有包含三维组织表面三角网格模型顶点的体素定义为边界体素。由于在步骤3中增加了三维组织表面三角网格顶点密度，因此本步骤所得到的边界体素完全包裹住三维组织表面三角网格，即三维组织表面网格上的任意部分都位于边界体素内。

6.由于在步骤5中得到的边界体素全包裹住三维组织表面三角网格，因此边界体素也将完全包裹所有内部体素。根据步骤5中所确定的边界体素，可对最小长方体包围盒沿着O-XYZ坐标系的X、Y、Z任一坐标轴方向逐层确定边界体素所包围的内部体素。例如沿着X轴从小到大每次增加一个体素边长，对每一层平行于OYZ平面的体素寻找其内部体素，每一层中的边界体素所包围的体素即为内部体素，参见图5。

7.将步骤3中填充顶点后的三维组织表面三角网格模型上每一个顶点定义为一个粒子，得到边界粒子；将步骤6中得到的每一个内部体素取其体素中心位置定义为一个粒子，得到内部粒子，这样便最终得到粒子化三维组织体模型。

步骤2，建立软组织、液体和刚体的粒子化力学模型，将软组织和液体定义为流变体，并对软组织粒子和液体粒子的任意两个粒子之间定义流变粘滞系数。

在真实世界中，几乎所有的生物软组织都是由固态和液态两种成分共同构成的，因此生物软组织是最典型的流变体，如生物粘液（体液）、肌肉、器官等。为了更真实的对手术现象进行仿真，本发明将软组织和液体定义为流变体，结合流变力学来建立软组织和液体的粒子化力学模型。对软组织和液体各自的任意两个粒子，引入软组织和液体各自粒子间的流变粘滞系数，该系数在常温常压下与软组织粒子间或者液体粒子间的切变速度有幂律关系。

液体的粒子化力学模型的采用Navier-Stokes方程组来描述其运动状态，采用光滑粒子流体动力学SPH(Smoothed Particle Hydrodynamics)方法求解Navier-Stokes方程，并通过引入流变粘滞系数计算液体粒子间的粘滞力为体液赋予流变性质。同样，对于软组织的粒子化力学模型，采用基本弹性力学方程求解软组织粒子各物理量的方法，采用SPH方法计算粘滞力，并代入流变粘滞系数为软组织粒子赋予流变性质。对于刚体粒子模型，因为刚体只有平移和转动两种状态，因此只采用表面一层粒子表示刚体。

因此，实施例的该粒子化力学模型的具体内容如下：

1.粒子模型：软组织、液体和刚体都表示成粒子模型；

2.流变性质：将软组织和液体定义为流变体，并对软组织粒子和液体粒子的任意两个粒子之间定义流变粘滞系数，如图2所示，本发明采用如下方法定义粒子化力学模型中两个粒子i和j之间的流变粘滞系数其中粒子i和j为软组织粒子或液体粒子：

${\mathrm{\&mu;}}_{j\&LeftRightArrow;i}=K{\left(\frac{\left|v_{\mathrm{\&tau;}}\right|}{|x_i-x_j|}\right)}^{\mathrm{\&beta;}}---\left(1\right)$

其中，K、β为与软组织和液体相关的系数，可根据实验预设取值，当两个粒子为液体粒子时，这两个量较小；当两个粒子为软组织粒子时，这两个量较大；当两个粒子一个是液体粒子而一个是软组织粒子时，则K、β分别取两个粒子都是软组织粒子和都是液体粒子时对应系数分别的平均值。即设当两个粒子为液体粒子时，K、β的取值为K1、β1，当两个粒子为软组织粒子时，K、β的取值为K2、β2，当两个粒子一个是液体粒子而一个是软组织粒子时，则K、β分别取（K1+K2）/2、(β1+β2）/2，K1<K2,β1<β2。一般，K可取1000≤K≤2000，β可取2.0≤β≤6.0，K和β的取值随着物质的粘性增大而取较大值。v_τ为粒子i相对粒子j运动切向方向的相对速度，x_i和x_j、v_i和v_j分别是粒子i和粒子j的位置和速度，则有：

$\left|v_{\mathrm{\&tau;}}\right|=|v_i-v_j|\sin \mathrm{\&theta;}=\sqrt{1-{\cos }^2\mathrm{\&theta;}}|v_i-v_j|---\left(2\right)$

其中，角度θ根据下式得到：

$\cos \mathrm{\&theta;}=\frac{(v_i-v_j)(x_j-x_i)}{|v_i-v_j||x_j-x_i|}---\left(3\right)$

3.刚体粒子模型力学建模：只有平移和转动，无形变，因此对刚体按上述基于三角网格模型的模型建立方法生成粒子化三维软组织体模型后，只采用边界粒子表示刚体，即只采用表面粒子模型，所以刚体粒子即为刚体边界粒子。由于刚体粒子之间保持相对静止，所以在刚体作为一个整体移动和转动。

为便于实施参考起见，提供刚体粒子模型力学建模的相关物理量计算公式如下：刚体受到的合外力设为

即为每个刚体边界粒子所受到的外力之和，即

$F_{\mathrm{rigid}}^{\mathrm{total}}=\underset{i=1}{\overset{n_r}{\mathrm{\&Sigma;}}}{f}_i---\left(4\right)$

其中n_r是刚体边界粒子的数量，f_i是第i个刚体粒子所受到的外力，此外力即为刚体粒子和其他类型粒子之间的耦合作用力和重力之和，每个刚体粒子的重力即为刚体粒子的质量乘以重力加速度。

刚体所受的外力除了对刚体有平动作用，还对刚体具有转动作用，外力作用在刚体上产生的转矩为：

${\mathrm{\&tau;}}_{\mathrm{rigid}}^{\mathrm{total}}=\underset{i=1}{\overset{n_r}{\mathrm{\&Sigma;}}}(x_i-x_r)\&times;f_i---\left(5\right)$

其中x_i和f_i分别是第i个刚体粒子的位置以及所受到的外力，x_r是所有刚体粒子的中心位置。计算出刚体所受的合外力以及转矩之后，可根据经典牛顿力学求解刚体在每个时刻的速度、加速度和位移。

4.软组织粒子模型力学建模：采用按上述基于三角网格模型的模型建立方法生成的粒子化三维软组织体模型表示软组织，采用现有经典弹性力学基本方程计算，并采用SPH方法求解各个粒子的物理量，例如弹力、粘滞力、重力、位移、速度和加速度，因三维软组织具有流变特性，因此采用经典SPH方法计算粘滞力，但是粘滞系数采用公式(1)定义的流变粘滞系数。

为便于实施参考起见，提供软组织粒子模型力学建模的相关物理量计算公式如下：

第i个软组织粒子的弹力f_i-elastic和粘滞力f_i-viscosity的计算公式如下：

$f_{i-\mathrm{elastic}}=-2\frac{m_i}{{\mathrm{\&rho;}}_i}{\mathrm{\&sigma;}}_i(I+\&dtri;u_i^T){\left(\underset{j=1}{\overset{n_i^s}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ij}}{x}_{\mathrm{ij}}^T{W}_s\right)}^{-1}(-\underset{j=1}{\overset{n_i^s}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ij}}{W}_s)---\left(6\right)$

$f_{i-\mathrm{vis}\cos \mathrm{ity}}=\underset{j=1}{\overset{n_i^s}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{j\&LeftRightArrow;i}(v_j-v_i)\frac{m_j}{{\mathrm{\&rho;}}_j}{\&dtri;}^2{W}_s---\left(7\right)$

其中C是弹性矩阵，u_i是第i个软组织粒子的位移，▽表示梯度算子，I表示单位矩阵，

表示第i个软组织粒子的作用域内软组织粒子的数量，x_ij表示第i个软组织粒子和它的作用域内的第j个软组织粒子的相对位置，ρ_i和ρ_j为软组织粒子密度，v_i和v_j为第i个软组织粒子和它的作用域内的第j个软组织粒子分别的速度，m_i和m_j为第i个软组织粒子和它的作用域内的第j个软组织粒子分别的质量，W_s表示软组织粒子间作用的核函数，

为第i个软组织粒子和它的作用域内的第j个软组织粒子之间的流变粘滞系数。每个软组织粒子的重力即为软组织粒子的质量乘以重力加速度。在计算完每个软组织粒子所受的弹力、粘滞力和重力之后，可根据经典牛顿力学求解软组织粒子在每个时刻的速度、加速度和位移。

5.液体粒子模型力学建模：液体采用现有的SPH方法建模得到粒子模型，而不是上述基于三角网格模型的模型建立方法，因此液体粒子不区分内部粒子、边界粒子，所有的液体粒子都是同一性质的。液体粒子的计算使用SPH方法求解描述液体运动的Navier-Stokes方程组，可得到各个粒子的物理量，例如压力、粘滞力、重力、位移、速度和加速度。因液体具有流变特性，因此具有流变体性质，采用经典SPH方法计算粘滞力，在求解粘滞力时采用公式(1)定义的流变粘滞系数。

为便于实施参考起见，提供液体粒子模型力学建模后相关物理量计算公式如下：

第i个液体粒子的粘滞力

和压力的计算公式如下：

$f_i^{\mathrm{vis}\cos \mathrm{ity}}=\underset{j=1}{\overset{n_i^f}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_{j\&LeftRightArrow;i}(v_j-v_i)\frac{m_j}{{\mathrm{\&rho;}}_j}{\&dtri;}^2{W}_f----\left(8\right)$

$f_i^{\mathrm{pressure}}=\underset{j=1}{\overset{n_i^f}{\mathrm{\&Sigma;}}}\frac{m_j({\mathrm{\&rho;}}_i+{\mathrm{\&rho;}}_j)}{2{\mathrm{\&rho;}}_j}\&dtri;W_f----\left(9\right)$

其中

表示第i个液体粒子的作用域内液体粒子的数量，ρ_i和ρ_j为第i个液体粒子和它的作用域内的第j个液体粒子分别的密度，v_i和v_j为第i个液体粒子和它的作用域内的第j个液体粒子分别的速度，m_j为第i个液体粒子的作用域内的第j个液体粒子质量，W_f表示液体粒子间作用的核函数，

为第i个液体粒子和它的作用域内的第j个液体粒子之间的流变粘滞系数。

每个液体粒子的重力即为液体粒子的质量乘以重力加速度。在计算完每个液体粒子所受的压力、粘滞力和重力之后，可根据经典牛顿力学求解液体粒子在每个时刻的速度、加速度和位移。

SPH方法、经典弹性力学基本方程、Navier-Stokes方程组为现有技术，本发明不予赘述。

步骤3，进行基于边界粒子和多相耦合核函数的碰撞检测及响应。

在虚拟手术仿真中，多相耦合过程包含多相物质的碰撞检测以及响应过程，在这个过程中，液体、软组织以及刚体粒子之间具有相互作用，它们之间的碰撞检测过程通过检测液体粒子、软组织粒子和刚体粒子模型中的表面粒子之间的碰撞可大大的提高碰撞检测的效率，碰撞响应则需要考虑到多相粒子直接的相互影响。

现有算法在进行碰撞检测时，一般选取离目标粒子最近的粒子作为碰撞点，但是当所选取的最近粒子当前时刻不在目标粒子运动轨迹速度方向上甚至是速度相反方向上时，很明显的是下一时刻目标粒子并不会与所选取的最近粒子相碰撞，这样则会产生不真实的碰撞现象。为了防止这种不真实的现象，实现与客观世界更接近的碰撞检测，本发明所提出的碰撞检测方法考虑了粒子的运动方向及轨迹对碰撞检测的影响。

虚拟手术仿真场景中通常会遇到多相耦合仿真，需要处理多相物质之间的碰撞，本发明首先对每相物质的粒子进行标记，以区分软组织、液体和刚体粒子的类型。对于软组织粒子和刚体粒子，可在粒子模型生成过程中标记边界粒子，然后只采用边界粒子来进行多相粒子间碰撞检测，对其中任何一相粒子，需要将其与其他另外两相粒子一一进行碰撞检测。找到碰撞点后，再根据碰撞响应施加在粒子上的作用力进行单相粒子模型的计算。

对于任意两相的物质A和B，A为软组织、液体或刚体，B为软组织或刚体，且A和B不相同，碰撞检测及响应方法的具体步骤如下：

碰撞检测过程为：

⑥初始化A型粒子半径和B型粒子半径，以A型粒子和B型粒子的坐标为各自圆心，以A、B两型粒子分别的半径作小球，如果液体粒子和软组织粒子圆心之间的聚类小于A、B两型粒子分别的半径之和，则发生碰撞；

⑦对每一个A型粒子，以A型粒子当前位置为圆心，一个时间步长内移动的距离为半径作球形搜索区域，检测出此球形区域内所有的B型边界粒子。时间步长可由本领域技术人员预先设定。计算可得一个时间步长内A型粒子由当前位置O应该移动到位置O′，以O为起点指向O′作一条射线；

⑧将步骤②中所得的每个B型边界粒子投影至步骤②中的射线所在的直线上，并求解每个B型边界粒子到相应投影点的投影距离；

⑨在步骤③中检测出所有满足以下条件的B型边界粒：其投影距离小于等于A型粒子半径与B型边界粒子半径之和；

⑩在步骤④检测出的所有的B型边界粒子中满足以下条件的B型边界粒子为即为碰撞点：该B型边界粒子的投影点位于步骤②中A型粒子当前运动方向的前方，并其投影点距离步骤②中A型粒子圆心O最近。

以A型粒子为软组织粒子和B型粒子为液体粒子为例，如图3所示，碰撞检测具体步骤为：

1.初始化液体粒子半径和软组织粒子半径，以液体粒子和软组织粒子的坐标为各自圆心，以两种粒子分别的半径作小球，可以认为，如果液体粒子和软组织粒子圆心之间的聚类小于其半径之和，则发生碰撞。p_f为液体粒子，其圆心为O，半径为r_f；p₁,p₂,…,p₇是软组织边界粒子，其圆心为S₁,S₂,…,S₇，软组织边界粒子半径为r_s；

2.考察液体粒子p_f，设p_f在一个时间步长内移动的距离大小为r₀，以r₀的大小为半径作球形搜索区域，检测出此球形搜索区域内所有的软组织边界粒子，当前液体粒子的速度为v_cur，一个时间步长内液体粒子移动到了O′，则在该球形搜索区域内的软组织边界粒子为p₁,p₂,p₃,p₄，l为以O为起点指向O′的一条射线；

3.将步骤2中所得的软组织边界粒子p₁,p₂,p₃,p₄投影至射线l所在的直线上，求解得每个软组织边界粒子的投影点为S₁′,S₂′,S₃′,S₄′，相应的到投影点的距离为S₁S₁′,S₂S₂′,S₃S₃′,S₄S₄′；

4.根据步骤3所得S₁S₁′,S₂S₂′,S₃S₃′,S₄S₄′检测出所有满足以下条件的软组织边界粒子：其投影距离小于等于液体粒子半径r_f与软组织边界粒子半径r_s之和，满足此条件的软组织边界粒子为p₁,p₂,p₃；

5.在步骤4检测出的所有软组织边界粒子p₁,p₂,p₃中，以满足以下条件的软组织边界粒子为碰撞点：位于步骤2中粒子p_f当前运动方向的前方，并其投影点距离步骤2中当前液体粒子圆心O最近。如图3所示，满足条件的碰撞点是p₁，因为p₁的投影点S₁′位于射线l上并且距离点O最近。

在单相粒子的条件下，同相粒子之间的作用包括压力、粘性力等，皆由此相粒子的性质决定。而多相耦合中，多相粒子之间相互接触过程中会受到彼此间相互作用力的影响，这种力阻碍对方粒子的运动，将其定义为粒子间的耦合作用力。同样，为了描述多相粒子间发生碰撞后的响应作用，本发明定义了多相耦合核函数，以定义多相粒子间的作用力。其设计原理是：单相粒子条件下，粒子间的作用力由单相粒子的核函数决定，例如软组织粒子间作用的核函数为W_s，液体粒子间作用的核函数为W_f，由于刚体粒子间相对位置保持不变，因此刚体粒子间没有核函数，为了方便计算实施例定义一个刚体粒子对外作用的初始核函数W_r，核函数W_s、W_f和W_r可采用SPH方法中典型的核函数。例如，软组织粒子间作用的核函数W_s可以采用式(10)：

$W_s(r_s,h_s)=\left\{\begin{array}{cc}\frac{1}{\mathrm{\&pi;}{h_s}^3}(1-\frac{3{l_s}^2}{2{h_s}^2}+\frac{3{l_s}^3}{4{h_s}^3}),& 0\&le;l_s<h_s\\ \frac{1}{4\mathrm{\&pi;}{h_s}^3}(2-\frac{l_s}{h_s}),& h_s\&le;l_s<{2h}_s\\ 0,& l_s\&GreaterEqual;{2h}_s\end{array}\right.---\left(10\right)$

流体粒子间作用的核函数W_f可以采用式(11)：

$W_f(r_f,h_f)=\left\{\begin{array}{cc}\frac{15}{\mathrm{\&pi;}{h_f}^6}{(h_f-l_f)}^3,& 0\&le;l_f<h_f\\ 0,& l_f\&GreaterEqual;h_f\end{array}\right.----\left(11\right)$

刚体粒子对外作用的初始核函数W_r可以采用式(12)：

$W_r(r_{\mathrm{rigid}},h_r)=\left\{\begin{array}{cc}\frac{4}{\mathrm{\&pi;}{h_r}^8}{({h_r}^2-l_{\mathrm{rigid}}^2)}^3,& 0\&le;l_{\mathrm{rigid}}<h_r\\ 0,& l_{\mathrm{rigid}}\&GreaterEqual;h_r\end{array}\right.---\left(12\right)$

其中l_s＝|r_s|，l_f＝|r_f|，l_rigid＝|r_rigid|，r_s,r_f,r_rigid分别是软组织、液体和刚体的中心粒子指向其他粒子的向量，h_s、h_f、h_r是软组织粒子、液体粒子和刚体粒子的光滑核半径。

在多相粒子条件下，则无法采用单相粒子的核函数来计算多相粒子间的作用力，这种情况下，本发明提出多相耦合核函数，即多相粒子间作用时，其作用核函数由各相粒子的核函数乘以各自权值线性组合而成，而此权值则为碰撞点作用域内各相粒子的数量所占比值。如图4所示，白色粒子为液体粒子，右斜线粒子为软组织粒子，左斜线粒子为刚体粒子（黑色实线包围区域内），指向两相物质的双箭头符号表示该两相物质间具有耦合作用。

碰撞响应过程为：

⑤根据碰撞检测过程，判断液体粒子、软组织粒子和刚体粒子中任意两相粒子是否发生碰撞，设该两相粒子为A型粒子和B型粒子，A为软组织、液体或刚体，B为软组织或刚体，且A和B不相同；即碰撞检测过程对某A型粒子找到了满足条件的B型边界粒子为碰撞点时，此A型粒子和B型边界粒子会发生碰撞。

⑥若A型粒子和B型粒子发生碰撞，则定义A型粒子和B型粒子耦合作用的核函数；

⑦根据步骤②中的耦合核函数计算A型粒子和B型粒子之间的耦合作用力，并根据冲量定理求解下一时刻A型粒子和B型粒子的速度大小和方向；

⑧根据步骤③的计算结果，代入步骤2中建立的A型粒子和B型粒子相应的粒子化力学模型，计算出各相粒子新的位移、速度。得到响应结果后，继续返回进行碰撞检测过程，根据两个粒子新的位移、速度继续寻找接下来的碰撞点，并根据碰撞响应过程计算新的位移、速度，直到结束循环。具体实施时，可由用户决定循环是否结束或者预设循环结束条件。

当A为软组织、液体或刚体，B为软组织或刚体，且A和B不相同，B型粒子的作用域半径是r_B，在确定所碰撞的B型粒子之后，定义A型粒子与B型粒子之间作用的耦合核函数为：

$W_{A\&LeftRightArrow;B}={\mathrm{\&delta;}}_A^{A\&LeftRightArrow;B}{W}_A+{\mathrm{\&delta;}}_B^{A\&LeftRightArrow;B}{W}_B---\left(13\right)$

其中

分别表示在B型粒子的作用域内A型粒子与B型粒子数目所占的比例，且n_A,n_B分别为B型粒子作用域内A型粒子和B型粒子的数目。W_A表示A型粒子间作用的核函数，W_B表示B型粒子间作用的核函数。

步骤③中，A型粒子和B型粒子的耦合作用力

采用公式(14)计算：

$f_{A\&LeftRightArrow;B}=-{\mathrm{\&mu;}}_{A\&LeftRightArrow;B}{m}_A{m}_B\&dtri;W_{A\&LeftRightArrow;B}---\left(14\right)$

其中m_A，m_B分别表示A型粒子与B型粒子的质量

为A型粒子与B型粒子之间作用的耦合核函数，▽表示梯度算子，

为根据公式(1)定义的A型粒子与B型粒子之间的流变粘滞系数，由于A型粒子和B型粒子是类型不相同的粒子，如果A型粒子与B型粒子一个类型是软组织粒子而一个类型是液体粒子时，则

的值为液体粒子和软组织粒子之间的流变粘滞系数，即粒子i和j分别为A型粒子与B型粒子，根据公式(1)得到的

作为

两个粒子i和j之间的流变粘滞系数

定义中，粒子i和j为软组织粒子或液体粒子，因此如果A型粒子与B型粒子其中有一个类型是刚体粒子时，那么

的值为另一个类型粒子之间的流变粘滞系数，即当A型粒子为刚体粒子时，粒子i和j分别为B型粒子与B型粒子，当B型粒子为刚体粒子时，粒子i和j分别为A型粒子与A型粒子，根据公式(1)得到作为

采用上述方法求解每相粒子所受的耦合作用力后，在一个步长的时间内，采用冲量定理进行碰撞响应的计算，求得每相粒子碰撞后的速度大小和方向。冲量定理为现有技术，本发明不予赘述。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (4)

1.一种虚拟手术中基于粒子的多相耦合方法，其特征在于，包括以下步骤：

步骤1，对于虚拟手术中涉及的软组织和刚体，分别根据三维组织表面三角网格模型生成粒子化三维组织体模型，包括以下子步骤，

①对三维组织表面三角网格模型建立最小长方体包围盒；

②将步骤①中最小长方体包围盒划分为若干大小相同的立方体的体素，在边界不足处采用立方体的体素补充，并以最小长方体包围盒左下角顶点为起点对所有立方体的体素编号，得到各体素的体素序号；

③在三维组织表面三角网格模型的每一个三角面片内部均匀填充顶点，得到填充顶点后的三维组织表面三角网格模型，并为所有顶点编号得到顶点序号；

④对步骤③中填充顶点后的三维组织表面三角网格模型，确定每一个顶点所位于的体素，并建立体素-顶点关系表，该关系表的内容是步骤②所得每个体素的体素序号及该体素所包含顶点的顶点序号；

⑤根据步骤④中的体素-顶点关系表，将所有包含顶点的体素定义为边界体素；

⑥根据步骤⑤中所确定的边界体素，对最小长方体包围盒沿着O-XYZ坐标系的任一坐标轴方向逐层确定边界体素所包围的内部体素，得到所有的内部体素；

⑦将步骤③中填充顶点后的三维组织表面三角网格模型上每一个顶点定义为一个粒子，得到边界粒子；并将步骤⑥中得到的每一个内部体素取中心位置定义为一个粒子，得到内部粒子；边界粒子和内部粒子构成粒子化三维组织体模型；

步骤2，建立软组织、液体和刚体的粒子化力学模型，将软组织和液体定义为流变体，并对软组织粒子和液体粒子的任意两个粒子之间定义流变粘滞系数；

所述刚体的粒子化力学模型，只采用步骤1所得刚体的粒子化三维组织体模型的边界粒子表示刚体，刚体作为一个整体移动和转动，并计算刚体所受的合外力以及转矩；

所述软组织的粒子化力学模型，采用步骤1所得软组织的粒子化三维软组织体模型表示软组织，根据弹性力学基本方程和SPH方法求解得到各个粒子的物理量，所述物理量包括弹力、粘滞力和重力，计算粘滞力时采用软组织粒子之间的流变粘滞系数；

所述液体的粒子化力学模型，采用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为软组织粒子或液体粒子，

${\mathrm{\&mu;}}_{j\&LeftRightArrow;i}=K{\left(\frac{\left|v_{\mathrm{\&tau;}}\right|}{|x_i-x_j|}\right)}^{\mathrm{\&beta;}}---\left(1\right)$ 其中，K、β为与软组织和液体相关的系数，设当两个粒子为液体粒子时，K、β的取值为K1、β1，当两个粒子为软组织粒子时，K、β的取值为K2、β2，当两个粒子一个是液体粒子而一个是软组织粒子时，则K、β分别取（K1+K2）/2、(β1+β2）/2，K1<K2,β1<β2；

v_τ为粒子i相对粒子j运动切向方向的相对速度，x_i和x_j、v_i和v_j分别是粒子i和粒子j的位置和速度，

$\left|v_{\mathrm{\&tau;}}\right|=|v_i-v_j|\sin \mathrm{\&theta;}=\sqrt{1-{\cos }^2\mathrm{\&theta;}}|v_i-v_j|---\left(2\right)$

其中，

$\cos \mathrm{\&theta;}=\frac{(v_i-v_j)(x_j-x_i)}{|v_i-v_j||x_j-x_i|}---\left(3\right)$ 角度θ根据上式得到。

3.根据权利要求2所述虚拟手术中基于粒子的多相耦合方法，其特征在于：A型粒子与B型粒子之间作用的耦合核函数如下，

$W_{A\&LeftRightArrow;B}={\mathrm{\&delta;}}_A^{A\&LeftRightArrow;B}{W}_A+{\mathrm{\&delta;}}_B^{A\&LeftRightArrow;B}{W}_B---\left(13\right)$

其中

分别表示在B型粒子的作用域内A型粒子与B型粒子数目所占的比例，且n_A,n_B分别为B型粒子作用域内A型粒子和B型粒子的数目，W_A表示A型粒子间作用的核函数，W_B表示B型粒子间作用的核函数。

4.根据权利要求3所述虚拟手术中基于粒子的多相耦合方法，其特征在于：A型粒子和B型粒子的耦合作用力

计算如下，

$f_{A\&LeftRightArrow;B}=-{\mathrm{\&mu;}}_{A\&LeftRightArrow;B}{m}_A{m}_A\&dtri;W_{A\&LeftRightArrow;B}---\left(14\right)$

其中m_A，m_B分别表示A型粒子与B型粒子的质量，

为A型粒子与B型粒子之间的流变粘滞系数。

Images