CN105912859A - 一种基于质点弹簧和流体力学的组织形变方法 - Google Patents

Abstract

一种基于质点弹簧和流体力学的组织形变方法，包括如下步骤：第一步，由医学数据对所需研究的组织建立均匀的三角面片网格；第二步，在构建的组织网格模型上，建立质点弹簧物理模型；第三步，建立定刚体核及广义弹簧；第四步，在构建完的质点弹簧网格模型内填充SPH模型流体；第五步，设定一个占据一定空间的刚性体，用于与流体一起传递负载；质点弹簧模型模拟组织表面的皮肤形变，SPH模型流体和刚性体则模拟组织内部体液的运行机制。本发明提供一种在具有良好实时性的同时提高仿真的真实性的基于质点弹簧和流体力学的组织形变方法。

Description

一种基于质点弹簧和流体力学的组织形变方法

技术领域

本发明涉及计算机图形学生物组织力学领域，具体涉及一种基于质点弹簧和流体力学的组织形变方法。

背景技术

计算机图形学的快速发展，为虚拟手术训练系统研发提供了坚实的基础与保障。外科医生可以通过在虚拟手术训练系统上模拟手术的过程，为其在实施真实的手术之前提供宝贵的经验，使外科手术变得更加安全、可靠、精确。

完整的虚拟手术系统，需要跨越很多的技术领域，涉及了很多学科，包括医学数据可视化、计算机图形学、牛顿动力学、生物力学、连续介质力学、流体力学、数值计算、物理建模与几何算法等领域，对计算机编程能力要求极高，是一个综合性极强的研究方向，一般具有以下两个特点：所要虚拟的场景及其复杂、需要高强度的实时人机交互。在虚拟手术中存在复杂的手术交互操作：碰撞、挤压、抓取、切割、吸除、剪断、针穿刺、灼烧、缝合等；同时为了达到高度逼真性，使医生具有沉浸其中的感觉，也需要模拟手术过程中的流血、灼烧烟雾等特殊的渲染效果。

软组织形变是手术交互操作中最基本核心的模块。碰撞、挤压等操作都需要基于弹性形变加以呈现，但不会改变组织的固有结构，因此对于碰撞、挤压主要需要研究的就是组织的形变特性，是手术仿真中最基本的交互方式。而切割、吸除等属于非线性形变，这种操作使软组织的几何结构发生改变，涉及到拓扑结构重组、断裂力学理论等，仿真难度较大，但基础理论还是软组织形变。

在虚拟手术系统中软组织行为的计算是研究的核心问题，国内外很多学者也提出了许多软组织建模的方法。这些方法中需要面临的最关键的问题是实时性与逼真度之间的协调关系。逼真性和实时性这对矛盾是软组织建模中最为主要的问题。逼真度越高，就需要所建立模型越精确复杂，而复杂的模型通常会导致计算时间的大幅增加，从而导致了实时性下降。

发明内容

为了克服已有虚拟手术系统中软组织行为的无法兼顾仿真者真实性和实时性的不足，本发明提供一种在具有良好实时性的同时提高仿真的真实性的基于质点弹簧和流体力学的组织形变方法。

本发明解决其技术问题所采用的技术方案是：

一种基于质点弹簧和流体力学的组织形变方法，所述组织形变方法包括如下步骤：

第一步，由医学数据对所需研究的组织建立均匀的三角面片网格；

第二步，在构建的组织网格模型上，建立质点弹簧物理模型；

第三步，建立定刚体核及广义弹簧；

第四步，在构建完的质点弹簧网格模型内填充SPH模型流体；

第五步，设定一个占据一定空间的刚性体，用于与流体一起传递负载；质点弹簧模型模拟组织表面的皮肤形变，SPH模型流体和刚性体则模拟组织内部体液的运行机制。

进一步，所述第一步中，所述的三角面片网格的建立，由医学所得CT、核磁共振数据进过图像处理的算法获得相应的三维网格模型。

再进一步，所述第二步中，所述的质点弹簧物理模型，是在三角面片网格基础上，以网格节点作为质点，各节点之间相连线段视为弹簧阻尼器，从而构建质点弹簧物理模型。

更进一步，所述的质点弹簧物理模型，其每个质点的动力学微分方程如下：

$m_i{\stackrel{\·\·}{x}}_i+c_i{\stackrel{\·}{x}}_i+k_{ij}(l_{ij}-l_{ij}^0)=f_{ext}$

其中，m_i为第i个质点(N_i)的质量；x_i为质点N_i的形变位移；c_i为阻尼系数；k_ij为弹簧弹性系数；为质点N_i与N_j间连接弹簧的初始长度；f_ext为质点N_i所受外力的总和。

所述第三步中，刚体核记录组织质点初始位置的不变点，广义弹簧是组织质点与刚体核之间的虚拟弹簧，初始长度为0。

所述第四步中，基于流体不可压缩假设，采用Navier-Stokes方程描述流体的运动：

$f_i=-\▿P+\mathrm{\μ}{\▿}^{2}V+\mathrm{\ρ}g$

$f_i^p=-\▿P=-\Σm_j\frac{P_i+P_j}{2{\mathrm{\ρ}}_j}\▿W(x_i-x_j,h)$

$f_i^v=\mathrm{\μ}{\▿}^{2}V=\mathrm{\μ}\Σm_j\frac{V_i-V_j}{{\mathrm{\ρ}}_j}{\▿}^{2}W(x_i-x_j,h)$

$f_i^g={\mathrm{\ρ}}_{i}g$

在上式中，W(x_i-x_j,h)表示光滑核函数，与分别表示光滑核函数的一阶和二阶导数；分别表示粒子的压力、粘滞力和重力，μ为粘滞系数，ρ为粒子的密度，P为粒子的压强，可由公式：P＝k(ρ-ρ₀)进行计算(其中，k为气体常数，ρ₀为初始密度)，V为粒子速度。

所述第五步中，设定一个将组织模型缩小1/2的刚性体，填充至组织内，并设定该刚性体的重心与组织模型的重心初始位置重合，并且在传递负载过程中刚性体的重心位置只在一定范围中运动，确保刚性体能够不与组织表面直接接触。

所述的占据一定空间的刚性体，其只能做刚性运动，将刚性体整体看作一个质点进行力的传递，而刚性体的表面则限制了流体的表面，从而达到刚性体与流体共同传递负载。

本发明的技术构思为：先在组织的表面网格的基础上建立质点弹簧模型，并添加刚体核及广义弹簧，来模拟组织表面的皮肤形变，再向模型体内填充SPH模型流体及一块刚性体，来模拟组织内部的体液运行机制，并用于传递负载。

如图1所示是在组织模型表面的三角网格面片的基础上建立质点弹簧模型，该模型将组织模型表面的几何坐标点看做有质量的节点,将各个节点相连的拓扑线段看成连接质点的弹簧，当质点受到外力的作用，试图产生位移时，就会受到与其相连的弹簧附加的收缩力，而外力负载也在弹簧的作用下传递到整个模型中去。经过一段的时间，整个弹簧系统达到动态平衡，或者说软组织形变的过程就是一种向平衡过度的过程，直到被外界的外力打破这种平衡，继而再次寻找平衡的过程。

质点弹簧模型的简单数学描述如下：

将软组织的质点弹簧模型抽象成由n个离散化后的质点和质点间连接的虚拟弹簧阻尼器构成。其中,用P_i(i＝1,2,......n)表示离散后的抽象质点；用l_ij(i,j∈[1,n],i≠j)表示质点间连接的弹簧阻尼器。整个质点弹簧模型在力学的控制下通过更新每个质点的位移从而描述了虚拟人体软组织的形变过程，而每个质点都是软组织的特性点,也称为控制点。

在任一时刻,质点弹簧模型模拟的组织形变可由n个微分方程来描述的,对于每一个抽象的质点的运动都满足如下的动力学微分方程:

$m_i{\stackrel{\·\·}{x}}_i+c_i{\stackrel{\·}{x}}_i+k_{ij}(l_{ij}-l_{ij}^0)=f_{ext}$

其中,m_i为第i个质点(P_i)的质量；x_i为质点P_i的形变位移；c_i为阻尼系数；k_ij为弹簧的弹性系数；为质点P_i与P_j间连接弹簧的初始长度；f_ext为质点P_i所受的外力的总和。

为保证组织模型在形变后不受外力作用下能够有恢复到原来的趋势，本发明中添加了刚体核及广义弹簧。如图2，刚体核是记录了组织质点初始位置的不变点，广义弹簧是组织质点与刚体核之间的虚拟弹簧，初始长度为0。当组织质点受力远离刚体核时，广义弹簧就被拉伸，产生使质点趋向刚体核的力。

在完成质点弹簧建模后，要向组织模型中填充入SPH流体，用于仿真生物组织体内的体液的运行机制，从而传递负载。

SPH的本质是一种基于粒子的插值方式，SPH方程的构造有两个关键的步骤，第一步为积分表示法，即“核函数逼近”；第二步为“粒子逼近”。整个流场变成了一系列粒子的表达，每个粒子又由其支持域内的粒子表达，从而使整个流场的力学量都由这些粒子来负载。如图3表示了粒子i的一个支持域示意图，从而将连续的流场离散成支持域中的有限个粒子的和。采用Navier-Stokes方程描述流体的运动，N-S方程如下：

$f_i=-\▿P+\mathrm{\μ}{\▿}^{2}V+\mathrm{\ρ}g$

$f_i^p=-\▿P=-\Σm_j\frac{P_i+P_j}{2{\mathrm{\ρ}}_j}\▿W(x_i-x_j,h)$

$f_i^v=\mathrm{\μ}{\▿}^{2}V=\mathrm{\μ}\Σm_j\frac{V_i-V_j}{{\mathrm{\ρ}}_j}{\▿}^{2}W(x_i-x_j,h)$

$f_i^g={\mathrm{\ρ}}_{i}g$

在上式中，W(x_i-x_j,h)表示光滑核函数，与分别表示光滑核函数的一阶和二阶导数；分别表示粒子的压力、粘滞力和重力，μ为粘滞系数，ρ为粒子的密度，P为粒子的压强，可由公式：P＝k(ρ-ρ₀)进行计算(其中，k为刚体系数，ρ₀为初始密度)，V为粒子速度。

光滑核函数的选取：

除了粘滞力和压力外的所有物理量，如密度、质量等，我们使用w_poly6(r,h)核函数，该函数详细描述为：

$W_{poly6}(r,h)=\frac{315}{64{\mathrm{\πh}}^9}\left\{\begin{array}{cc}{(h^2-r^2)}^3& 0\≤r\≤h\\ 0& other\end{array}\right.$

我们使用w_spiky(r,h)核函数来计算压力，该函数详细描述为：

$W_{spiky}(r,h)=\frac{15}{{\mathrm{\πh}}^6}\left\{\begin{array}{cc}{(h-r)}^3& 0\≤r\≤h\\ 0& other\end{array}\right.$

我们使用w_viscosity(r,h)核函数计算粒子间的粘滞力，该函数详细描述为：

$W_{vis\cos ity}(r,h)=\frac{15}{2{\mathrm{\πh}}^3}\left\{\begin{array}{cc}-\frac{r^3}{2h^3}+\frac{r^2}{h^2}+\frac{h}{2r}-1& 0\≤r\≤h\\ 0& other\end{array}\right.$

本发明为了避免因流体的粒子数量过多而导致流体运算时间增加，影响实时性，设计了一个不能形变的刚性体，与流体一起填充入组织模型中，在模型中占据一定空间，并与流体一同传递负载，从而降低了流体粒子的数量，节省了运算时间。

本发明的有益效果主要表现在：在质点弹簧模型的基础上，用流体填充到组织模型内，以模拟生物体组织中液体的运作形式。该方法既保留了质点弹簧模型实时性好的优点，又改进了质点弹簧模型是由局部到整体来传递力的问题，并提高了仿真的真实性。

附图说明

图1为本发明的组织质点弹簧模型示意图。

图2为本发明的刚体核及广义弹簧示意图。

图3为本发明的流体粒子支持域示意图。

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图3，一种基于质点弹簧和流体力学的组织形变方法，所述组织形变方法包括如下步骤：

第一步，由医学数据对所需研究的组织建立均匀的三角面片网格；

第二步，在构建的组织网格模型上，建立质点弹簧物理模型；

第三步，建立定刚体核及广义弹簧；

第四步，在构建完的质点弹簧网格模型内填充SPH模型流体；

第五步，设定一个占据一定空间的刚性体，用于与流体一起传递负载；质点弹簧模型模拟组织表面的皮肤形变，SPH模型流体和刚性体则模拟组织内部体液的运行机制。

所述的质点弹簧模型模拟组织表面的皮肤形变，流体和刚性体模拟组织内部体液的运行机制，既保留了质点弹簧模型实时性好的优点，又改进了质点弹簧模型是由局部到整体来传递力的问题，并提高了仿真的真实性。

本实施例是针对虚拟手术系统中快速逼真的组织形变的方法，其先在组织的表面网格的基础上建立质点弹簧模型，并添加刚体核及广义弹簧，来模拟组织表面的皮肤形变，再向模型体内填充SPH模型流体及一块刚性体，来模拟组织内部的体液运行机制，并用于传递负载，最终在质点弹簧模型和流体的物理力学作用下达到一个动态平衡状态。该方法的具体步骤：

第一步：构建组织表面的三角网格。

第二步：如图1所示是在组织模型表面的三角网格面片的基础上建立质点弹簧模型，该模型将组织模型表面的几何坐标点看作有质量的节点,将各个节点相连的拓扑线段看成连接质点的弹簧，当质点受到外力的作用，试图产生位移时，就会受到与其相连的弹簧附加的收缩力，而外力负载也在弹簧的作用下传递到整个模型中去。经过一段的时间，整个弹簧系统达到动态平衡，或者说软组织形变的过程就是一种向平衡过度的过程，直到被外界的外力打破这种平衡，继而再次寻找平衡的过程。

质点弹簧模型的简单数学描述如下：

将软组织的质点弹簧模型抽象成由n个离散化后的质点和质点间连接的虚拟弹簧阻尼器构成。其中,用P_i(i＝1,2,......n)表示离散后的抽象质点；用l_ij(i,j∈[1,n],i≠j)表示质点间连接的弹簧阻尼器。整个质点弹簧模型在力学的控制下通过更新每个质点的位移从而描述了虚拟人体软组织的形变过程，而每个质点都是软组织的特性点,也称为控制点。

在任一时刻,质点弹簧模型模拟的组织形变可由n个微分方程来描述的,对于每一个抽象的质点的运动都满足如下的动力学微分方程:

$m_i{\stackrel{\·\·}{x}}_i+c_i{\stackrel{\·}{x}}_i+k_{ij}(l_{ij}-l_{ij}^0)=f_{ext}$

其中,m_i为第i个质点(P_i)的质量；x_i为质点P_i的形变位移；c_i为阻尼系数；k_ij为弹簧的弹性系数；为质点P_i与P_j间连接弹簧的初始长度；f_ext为质点P_i所受的外力的总和。

第三步：建立刚体核及广义弹簧。为保证组织模型在形变后不受外力作用下能够有恢复到原来的趋势，本发明中添加了刚体核及广义弹簧。如图2，刚体核是记录了组织质点初始位置的不变点，广义弹簧是组织质点与刚体核之间的虚拟弹簧，初始长度为0。当组织质点受力远离刚体核时，广义弹簧就被拉伸，产生使质点趋向刚体核的力。

第四步：在完成质点弹簧建模后，要向组织模型中填充入SPH流体，用于仿真生物组织体内的体液的运行机制，从而传递负载。

SPH的本质是一种基于粒子的插值方式，SPH方程的构造有两个关键的步骤，第一步为积分表示法，即“核函数逼近”；第二步为“粒子逼近”。整个流场变成了一系列粒子的表达，每个粒子又由其支持域内的粒子表达，从而使整个流场的力学量都由这些粒子来负载。如图3表示了粒子i的一个支持域示意图，从而将连续的流场离散成支持域中的有限个粒子的和。采用Navier-Stokes方程描述流体的运动，N-S方程如下：

$f_i=-\▿P+\mathrm{\μ}{\▿}^{2}V+\mathrm{\ρ}g$

$f_i^p=-\▿P=-\Σm_j\frac{P_i+P_j}{2{\mathrm{\ρ}}_j}\▿W(x_i-x_j,h)$

$f_i^v=\mathrm{\μ}{\▿}^{2}V=\mathrm{\μ}\Σm_j\frac{V_i-V_j}{{\mathrm{\ρ}}_j}{\▿}^{2}W(x_i-x_j,h)$

$f_i^g={\mathrm{\ρ}}_{i}g$

在上式中，W(x_i-x_j,h)表示光滑核函数，与分别表示光滑核函数的一阶和二阶导数；分别表示粒子的压力、粘滞力和重力，μ为粘滞系数，ρ为粒子的密度，P为粒子的压强，可由公式：P＝k(ρ-ρ₀)进行计算(其中，k为刚体系数为(10²～10³m²/s²)，ρ₀为初始密度，在本发明中取(1.05～1.06×10³kg/m³))，V为粒子速度。

第五步：为了避免因流体的粒子数量过多而导致流体运算时间增加，影响实时性，设计了一个不能形变的刚性体，与流体一起填充入组织模型中，在模型中占据一定空间，并与流体一同传递负载，从而降低了流体粒子的数量，节省了运算时间。设定一个将组织模型缩小1/2的刚性体，填充至组织内，并设定该刚性体的重心与组织模型的重心初始位置重合，并且在传递负载过程中限制刚性体的重心位置只在一定范围中运动，确保刚性体能够不与组织表面直接接触。

以上描述了本发明的基本原理和主要特征以及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (8)

1.一种基于质点弹簧和流体力学的组织形变方法，其特征在于：所述组织形变方法包括如下步骤：

第一步，由医学数据对所需研究的组织建立均匀的三角面片网格；

第二步，在构建的组织网格模型上，建立质点弹簧物理模型；

第三步，建立定刚体核及广义弹簧；

第四步，在构建完的质点弹簧网格模型内填充SPH模型流体；

第五步，设定一个占据一定空间的刚性体，用于与流体一起传递负载；质点弹簧模型模拟组织表面的皮肤形变，SPH模型流体和刚性体则模拟组织内部体液的运行机制。

2.如权利要求1所述的一种基于质点弹簧和流体力学的组织形变方法，其特征在于：所述第一步中，所述的三角面片网格的建立，由医学所得CT、核磁共振数据进过图像处理的算法获得相应的三维网格模型。

3.如权利要求1或2所述的一种基于质点弹簧和流体力学的组织形变方法，其特征在于：所述第二步中，所述的质点弹簧物理模型，是在三角面片网格基础上，以网格节点作为质点，各节点之间相连线段视为弹簧阻尼器，从而构建质点弹簧物理模型。

4.如权利要求3所述的一种基于质点弹簧和流体力学的组织形变方法，其特征在于：所述的质点弹簧物理模型，其每个质点的动力学微分方程如下：

$m_i{\stackrel{\·\·}{x}}_i+c_i{\stackrel{\·}{x}}_i+k_{ij}(l_{ij}-l_{ij}^0)=f_{ext}$

其中，m_i为第i个质点(N_i)的质量；x_i为质点N_i的形变位移；c_i为阻尼系数；k_ij为弹簧弹性系数；为质点N_i与N_j间连接弹簧的初始长度；f_ext为质点N_i所受外力的总和。

5.如权利要求1或2所述的一种基于质点弹簧和流体力学的组织形变方法，其特征在于：所述第三步中，刚体核记录组织质点初始位置的不变点，广义弹簧是组织质点与刚体核之间的虚拟弹簧，初始长度为0。

6.根据权利要求1或2所述的一种基于质点弹簧和流体力学的组织形变方法，其特征在于：所述第四步中，基于流体不可压缩假设，采用Navier-Stokes方程描述流体的运动：

$f_i=-\▿P+\mathrm{\μ}{\▿}^{2}V+\mathrm{\ρ}g$

$f_i^p=-\▿P=-{\mathrm{\Σm}}_j\frac{P_i+P_j}{2{\mathrm{\ρ}}_j}\▿W(x_i-x_j,h)$

$f_i^v=\mathrm{\μ}{\▿}^{2}V={\mathrm{\μ\Σm}}_j\frac{V_i-V_j}{{\mathrm{\ρ}}_j}{\▿}^{2}W(x_i-x_j,h)$

f_i ^g＝ρ_ig

在上式中，W(x_i-x_j,h)表示光滑核函数，与分别表示光滑核函数的一阶和二阶导数；f_i ^p、f_i ^v、f_i ^g分别表示粒子的压力、粘滞力和重力，μ为粘滞系数，ρ为粒子的密度，P为粒子的压强，可由公式：P＝k(ρ-ρ₀)进行计算，k为气体常数，ρ₀为初始密度，V为粒子速度。

7.如权利要求1或2所述的一种基于质点弹簧和流体力学的组织形变方法，其特征在于：所述第五步中，设定一个将组织模型缩小1/2的刚性体，填充至组织内，并设定该刚性体的重心与组织模型的重心初始位置重合，并且在传递负载过程中刚性体的重心位置只在一定范围中运动，确保刚性体能够不与组织表面直接接触。

8.如权利要求7所述的一种基于质点弹簧和流体力学的组织形变方法，其特征在于：所述的占据一定空间的刚性体，其只能做刚性运动，将刚性体整体看作一个质点进行力的传递，而刚性体的表面则限制了流体的表面，从而达到刚性体与流体共同传递负载。