CN102044086A - 一种软组织形变仿真方法 - Google Patents

Abstract

本发明涉及一种基于光滑粒子流体动力学的软组织形变仿真方法，属于图形处理技术领域，该方法选取光滑粒子流体动力学法，以黏弹性力学模型来反映软组织的生物力学特性，包含以下步骤：依据黏弹性模型，构建软组织形变仿真计算相关的一系列方程；选择合适的支持域搜索策略和光滑核函数，采用粒子近似法对方程组的各相关项进行近似计算，通过显示积分法计算各粒子的密度、位置、速度等随时间的变化值；动态将粒子模型每个时间步长的状态输出到屏幕上，并进行纹理光照的渲染，显示软组织器官受力情况下的实时形变过程。该方法无需繁琐的网格计算，可提高软组织形变仿真的准确性和实时性。

Description

一种软组织形变仿真方法

技术领域

本发明涉及虚拟手术技术领域，具体地说是一种基于光滑粒子流体动力学的软组织形变仿真方法，属于图形处理技术领域的物体形变实时模拟技术领域，其IPC专利分类号为G06T17/00。

背景技术

虚拟手术是虚拟现实技术在现代医学中的应用，是一个多学科交叉的研究领域。它主要由医学数据可视化与建模、人体软组织器官受力形变仿真两部分构成，在视觉上与触觉感官上为使用者提供了手术场景的真实再现。可用于外科医生的培训、手术结果的预测、手术计划的辅助制定、手术导航等，具有重大的理论意义和广阔的应用前景。

软组织形变仿真是虚拟手术中的最重要技术之一，能否逼真、实时的模拟软组织器官在外力作用下的形变是整个系统的关键。

常见的软组织形变计算模型分为两大类：基于几何的形变模型和基于物理的形变模型。其中，第一类模型仅仅考虑了几何形态的变化，而忽略了软组织的实际力学本构方程以及形变过程中物体质量、力或其他物理现象的作用，因此不能真实的反映软组织的形变过程，该模型目前已较少使用。而物理模型则基于软组织的力学本构方程，通过相应的计算模型得出组织受力时的形变，能够更加真实的反映组织的形变，因此目前使用较多。质点-弹簧模型(Mass-Spring Model)是较早出现的物理模型，它由无质量的弹簧和阻尼器网格连接的质点集合组成，该模型建模简单、易于实现，且计算速度相对较快；但很难通过实验来获取模型中成千上万的弹簧、质点和阻尼器参数以反映软组织内部真实的力学特性，因此，使用该模型实现形变仿真的逼真度不高；另外，前期需要耗费很多精力来构造弹簧阻尼器网格，仿真计算对网格结构依赖程度大，并且不能准确描述大变形问题。有限元模型(Finite Element Model，简写为FEM)是另一种常用的有网格物理模型，它将软组织离散为若干子单元，通过这些子单元边界上的节点相互结成为组合体，并利用变分原理建立求解应变量的代数方程组，从而计算形变。用有限元模型实现软组织形变仿真的方法已在近几年得到广泛的研究和应用，该方法具有坚实的弹性力学基础，较之其他有网格模型计算精度高，但它的计算复杂度大，很难实现软组织形变的实时模拟。边界元模型(Boundary Element Model，简写为BEM)可以看作对有限元模型的改进，它只对模型的边界进行离散，以降低问题的维数，简化计算。该方法无需考虑内部节点位移，较有限元法计算简单，但这种方法只能解决具有均质各向同性的线性问题。

形变仿真中软组织所采用的力学模型是形变仿真的基础与关键，决定形变仿真的准确性。当前常见的仿真系统为了能够实现组织的实时形变仿真，大多采用较为简化的力学模型，如：弹性力学模型、线弹性力学模型等，这些类型的模型大多直接采用经典材料力学中的力学模型，与软组织的实际力学特性差异很大，因此尽管可以实现能够满足一定速度需求的仿真，但精度很低，形变效果不如人意。

如前所述，现在常用的有网格物理模型需要耗费很大的精力来构造网格模型，后续的计算过程大都紧密的依赖于这种网络结构，并且计算复杂度往往较大，因此无法满足实时性、健壮性、精确性的要求，更难以准确描述软组织大变形问题；另外，经典的、简化的材料力学模型不能真实反映软组织的生物力学特性，依据这种力学模型对软组织的形变进行仿真计算，逼真度也很低。软组织形变仿真是虚拟手术系统中的关键环节，为此，如何设计一个性能更加良好的方法来实现软组织的形变仿真，以尽可能的满足系统关于实时性、健壮性、精确性的要求，成为虚拟手术系统所面临的首要问题。

发明内容

本发明的目的是解决虚拟手术系统中软组织形变仿真的实时性、逼真性和大变形仿真问题，克服了现有技术的不足，以较能真实反映软组织力学特性的黏弹性模型为力学模型，提供一种基于无网格法的软组织形变仿真方法，并实现其中最关键的基础部分，使其满足虚拟手术系统的需要。

为完成本发明的目的，本发明采取的技术方案是：以黏弹性力学模型反映软组织的生物力学特性，采用无网格连接的粒子模型来表示软组织形变计算模型，并选取光滑粒子流体动力学法(Smoothed Particle Hydrodynamics，简写为SPH)作为无网格模型的计算方法，完成软组织形变的仿真计算。其包含以下几个步骤：

步骤1)：采集软组织的数据信息；

步骤2)：选择黏弹性模型，构建用于软组织形变仿真的方程组：应力-应变本构方程(1)、(2)、(3)、应变-位移几何方程(4)和用于加速度计算的动量方程(5)；

步骤3)：依据步骤1)采集的软组织的数据信息，利用各数据点的位置向量，构建一个没有网格连接的粒子模型，并初始化所述粒子模型中各粒子的位置、质量、速度、加速度、受作用力等信息，构建粒子模型的初始化状态；

步骤4)：定义空间栅格，采用链表搜索法搜索粒子的支持域，构建支持域内的光滑核函数W；核函数W选取三次B-Spline光滑函数(6)，对核函数在各方向求导，得到核函数的一阶导数(7)；

步骤5)：应用光滑核函数W及其导数对参考粒子的支持域内所有粒子函数的加权平均近似的方法来构建步骤2)中各方程的SPH格式：

构建密度方程SPH格式(8)；

将动量方程(5)转换为SPH格式(9)；

构建几何方程的SPH格式(10)；

步骤6)：用显示积分法求解步骤5)的常微分方程，计算出粒子的密度、位置、速度等随时间的变化值；

步骤7)：循环执行上述步骤3)～步骤6)，计算出各粒子的状态；

步骤8)：将粒子模型的当前状态输出之屏幕，经过纹理和光照渲染后，得到软组织的动态形变过程。

进一步地，步骤6)的求解过程如下：

1)对模型的某一质点p施加外力f_outp；

2)对模型中的各粒子p_i(i＝1，2，..，M)循环完成下述3)～7)步骤的运算处理；

3)对当前粒子p_i，以h为光滑半径，搜索支持域内的相邻粒子p_j(j＝1，2，..，N)；

4)使用公式(6)(7)，计算当前粒子p_i与支持域内各邻近粒子p_j(j＝1，2，..，N)间的光滑核函数W_ij及其导数；

5)使用方程(8)，采用密度求和法计算粒子的密度ρ_i；

6)计算粒子的加速度a_i，加速度的计算采用如下方法：

$a_i=-\frac{{\mathrm{Dv}}_i}{\mathrm{Dt}}+\frac{{\mathrm{\σ}}_i}{m_i}+\frac{F_i}{m_i}$

其中，F_i表示粒子p_i所受的重力、外力等，σ_i表示粒子p_i所受的内力，即应力，m_i为粒子p_i的质量，使用公式(9)计算

7)计算粒子的速度和位移：

v_i(t)＝v_i(t)+dt·a_i

x_i(t)＝x_i(t)+dt·v_i(t)

其中，v_i(t)为粒子p_i在t时刻的速度向量，x_i(t)为粒子p_i在t时刻的位置向量，a_i为粒子p_i的加速度，dt为时间增量；

8)计算各粒子与初始状态时的位移情况disp_i(t)＝x_i(t)-x_i(t₀)；

9)使用几何方程(4)和(10)，由上一步算出的位移计算各粒子的应变状态

10)记录各粒子当前的体积应变

和形状畸变

并使用方程(2)，由上一步计算出的应变计算各粒子的新的体积应变

和形状畸变

11)计算

并使用黏弹性应力-应变本构方程(3)，计算各粒子的体积应力

和偏应力

12)使用方程(1)，计算出各粒子的应力状态

黏弹性是指在一定条件下同时具有弹性固体和黏性流体两者特征的特性。著名生物力学家冯元桢院士针对肌肉、器官、血管、骨骼等组织的力学特性进行了大量的试验和分析，并指出，几乎所有的生物固体都是黏弹性体，只不过程度上有所差异。软组织的力学特性主要表现为：黏弹性特性显著，比较容易变形，具有一定的抗拉能力，在拉伸载荷作用下往往显示出大变形特性。因此，形变仿真计算中，采用黏弹性力学模型较之目前广泛采用的弹性力学模型更能真实的反应软组织力学特性，仿真计算的准确性和逼真度更高。

光滑粒子流体动力学是一种无网格数值计算方法，其基本思想是在求解区域上任意设置有限个结点，采用结点权函数(或核函数)来近似表示其支持域内的位移函数和物理场函数，进而形成与结点位移和结点物理场相关的系统刚度方程，进行求解。与有限元等有网格法相比，它免除了定义在求解区域上的网格结构，避免了繁琐的网格剖分计算，不受网格结构束缚，可以方便的在求解域内增加和减少结点，因此有更高的计算精度，能够更精确的表示软组织器官的各种复杂几何形状，并且能够很容易的求解软组织大变形问题。如上所述，软组织是同时具有弹性固体和黏性流体两者特征的黏弹性材料，而光滑粒子流体动力学是既适用于固体，又适用于流体的计算方法。因此，用光滑粒子流体动力学实现软组织形变的仿真计算较传统方法具有更高的实时性和准确性。

附图说明

图1为本发明的总体处理流程图；

图2为本发明所采用的黏弹性力学模型示意图；

图3为肝脏器官的粒子模型初始状态图；

图4为本发明采用光滑粒子流体动力学法来实现软组织形变仿真计算的概念图；

图5为搜索粒子的支持域时所用到的粒子模型空间栅格划分示意图；

图6为链表搜索策略中的粒子索引存储示意图；

图7为肝脏器官软组织模型分别在受到拉力和压力情况下的形变效果图。

具体实施方式

本发明的总体处理流程如图1所示，其关键步骤就是构建基于黏弹性力学的控制方程及利用光滑粒子流体动力学法来实现无网格粒子模型的运动状态计算。下面按照总体流程图的顺序详细描述每一步的实施方式。

1.构建软组织形变仿真的控制方程

几乎所有的生物固体都是黏弹性体，本发明克服传统技术为了单一追求实时性而采用简单弹性力学模型来反映软组织力学特性的缺陷，力图提高形变仿真的准确性，并同时兼顾实时性，采用线性黏弹性力学方程来反映软组织的力学特性。

采用图2所示的Kelvin黏弹性模型，用弹簧和黏壶的并联来表示黏弹性。

构建其黏弹性应力-应变本构方程为：

$\mathrm{\σ}=q_0\mathrm{\ϵ}+q_1\stackrel{\·}{\mathrm{\ϵ}}$

其中，σ是应力向量，可由正应力分量σ_xx，σ_yy，σ_zz和剪应力分量σ_xy，σ_yz，σ_zx6个应力分量来表示，所谓应力是作用于单位截面上的描述物体中一部分材料作用于另一部分材料的内力，其单位是N/m²(帕斯卡)；ε是应变向量，由正应变分量ε_xx，ε_yy，ε_zz和剪应变分量ε_xy，ε_yz，ε_zx6个应变分量表示，所谓应变是物体内与应力有关的变形；q₀＝E，q₁＝η，E为材料的弹性模量，η为黏性系数，为应变率。

上述模型是一维模型，将其推广到三维情况。各向同性材料的应力张量σ可分解成它的球形张量和偏斜张量部分，应变张量ε可分离为体积形变和等体积的形状畸变两部分。即：

${\mathrm{\σ}}_{\mathrm{\α\β}}=S_{\mathrm{\α\β}}+\frac{{\mathrm{\δ}}_{\mathrm{\α\β}}{\mathrm{\σ}}_{\mathrm{kk}}}{3}---\left(1\right)$

${\mathrm{\ϵ}}_{\mathrm{\α\β}}=e_{\mathrm{\α\β}}+\frac{{\mathrm{\δ}}_{\mathrm{\α\β}}{\mathrm{\ϵ}}_{\mathrm{kk}}}{3}---\left(2\right)$

α，β＝x，y，z.式中σ_αβ为应力分量，ε_αβ为应变分量；δ_αβ为Kronecker符号，σ_kk＝σ_xx+σ_yy+σ_zz和ε_kk＝ε_xx+ε_yy+ε_zz分别为体积应力和体积应变；S_αβ和e_αβ分别为偏应力张量和偏应变张量的分量。根据Kelvin模型，偏应力张量和偏应变张量之间、体积应力和体积应变之间的三维黏弹性本构关系可表示为：

$S_{\mathrm{\α\β}}=E\·e_{\mathrm{\α\β}}+\mathrm{\η}\·\frac{{\mathrm{de}}_{\mathrm{\α\β}}}{\mathrm{dt}},\mathrm{\α},\mathrm{\β}=x,y,z$

${\mathrm{\σ}}_{\mathrm{kk}}=E\·{\mathrm{\ϵ}}_{\mathrm{kk}}+\mathrm{\η}\·\frac{{\mathrm{d\ϵ}}_{\mathrm{kk}}}{\mathrm{dt}}---\left(3\right)$

上式中，S_αβ和e_αβ分别为偏应力张量和偏应变张量的分量；σ_kk和ε_kk分别为体积应力和体积应变；E为材料的弹性模量，η为黏性系数；

分别为偏应变分量和体积应变对时间的导数，即应变率，t为时间。

应变-位移方程采用弹性力学中的几何方程：

${\mathrm{\ϵ}}_{\mathrm{xx}}=\frac{\∂u}{\∂x}$ ${\mathrm{\ϵ}}_{\mathrm{xy}}=\frac{\∂u}{\∂y}+\frac{\∂v}{\∂x}$

${\mathrm{\ϵ}}_{\mathrm{yy}}=\frac{\∂v}{\∂y}$ ${\mathrm{\ϵ}}_{\mathrm{yz}}=\frac{\∂v}{\∂z}+\frac{\∂w}{\∂y}---\left(4\right)$

${\mathrm{\ϵ}}_{\mathrm{zz}}=\frac{\∂w}{\∂z}$ ${\mathrm{\ϵ}}_{\mathrm{zx}}=\frac{\∂w}{\∂x}+\frac{\∂u}{\∂z}$

其中u，v，w为位移在三个坐标方向的分量；ε_xx，ε_yy，ε_zz为正应变分量；ε_xy，ε_yz，ε_zx为剪应变分量。

构建用于加速度计算的动量方程，其形式如下：

$\frac{\mathrm{dv}}{\mathrm{dt}}=\frac{1}{\mathrm{\ρ}}\frac{\∂\mathrm{\σ}}{\∂x}---\left(5\right)$

其中，v为速度向量，t为时间，ρ为粒子的密度，σ为粒子的应力，x为坐标向量；

2.粒子模型初始化

依据软组织的数据信息，利用各数据点的位置向量，构建一个没有网格连接的粒子模型。并初始化模型中各粒子的位置、质量、速度、加速度、受作用力等信息。图3是以肝脏器官为例，构建粒子模型的初始化状态。

3.SPH法求解粒子模型的运动信息

图4是本发明采用光滑粒子流体动力学法来实现软组织形变仿真计算的概念图。其核心包含两个过程：函数的光滑近似逼近和粒子近似逼近。函数的光滑近似逼近将描述粒子的密度、速度、位移等的函数表示为积分形式；然后通过粒子近似逼近，用影响半径内邻近粒子的速度、位移等运动特征求和平均近似代替参考粒子的运动信息，即用光滑核函数W影响范围内的邻近粒子p_j的运动信息求和平均代替参考粒子p_i的运动状态。图4中，h为影响半径的光滑长度，W为光滑核函数，p_i为参考粒子，p_j为邻近粒子。

光滑核函数的形式很多，本发明计算过程中选取Monaghan和Lattanzio在三次样条函数基础上提出的B-样条光滑函数。构建其三维格式方程如下：

$W_{\mathrm{ij}}=\frac{1}{{\mathrm{\&pi;h}}^3}\&times;\left\{\begin{array}{cc}1-\frac{3}{2}{s}^2+\frac{3}{4}{s}^3& 0\&le;s<1\\ \frac{1}{4}{(2-s)}^3& 1\&le;s<2\\ 0& s\&GreaterEqual;2\end{array}\right.---\left(6\right)$

式中，W_ij为由邻近粒子p_j估计粒子p_i运动信息的光滑核函数，

为粒子p_i与p_j之间的相对距离，r为粒子p_i与p_j之间的距离，h为光滑长度。核函数的一阶导数为：

$\frac{{\&PartialD;W}_{\mathrm{ij}}}{{\&PartialD;x}^{\mathrm{\&beta;}}}=\frac{x_i^{\mathrm{\&beta;}}-x_j^{\mathrm{\&beta;}}}{r}\&CenterDot;\frac{1}{{\mathrm{\&pi;h}}^4}\&times;\left\{\begin{array}{cc}-3s+\frac{9}{4}{s}^2& 0\&le;s<1\\ -\frac{3}{4}{(2-s)}^2& 1\&le;s<2\\ 0& s\&GreaterEqual;2\end{array}\right.---\left(7\right)$

式中，上标β＝x，y，z表示坐标方向，s，r，h的含义同式(6)，

和

分别表示粒子p_i与p_j的位置坐标向量在各方向的分量。

另外，支持域内邻近粒子的搜索，通过定义空间栅格，采用链表搜索法实现，其基本思想是：将整个粒子空间划分成规则的栅格单元，每个粒子都分布在一个栅格单元Cell中，Cell记录了在栅格中的粒子。在粒子搜索时，将光滑半径设为栅格单元的长度和搜索半径，则只需要搜索邻近的栅格即可确定目标粒子在光滑半径之内的粒子，将这些粒子标记为邻接粒子，并储存粒子索引。这样，搜索的范围仅仅限定在中心栅格单元周围这些栅格单元上，而对于其他粒子则不需要再考虑，大大降低了时间开销，提高了效率。图5、6分别是空间栅格划分及粒子索引存储示意图。

依据粒子近似的思想，通过应用光滑核函数W及其导数对参考粒子的支持域内所有粒子函数的加权平均近似的方法来构建各方程的SPH格式。

构建密度方程SPH格式

${\mathrm{\&rho;}}_i=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{m}_j{W}_{\mathrm{ij}}---\left(8\right)$

其中，ρ_i为粒子p_i的密度，m_j为p_i的支持域内邻近粒子p_j的质量，N为粒子p_i的支持域内粒子总数，W_ij的含义同式(6)。

将动量方程(5)转换成SPH格式

$\frac{{\mathrm{Dv}}_i^x}{\mathrm{Dt}}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}[m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{xx}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{xx}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{xy}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{xy}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{xz}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{xz}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}]$

$\frac{{\mathrm{Dv}}_i^y}{\mathrm{Dt}}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}[m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{yy}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{yy}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{yx}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{yx}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{yz}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{yz}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}]---\left(9\right)$

$\frac{{\mathrm{Dv}}_i^z}{\mathrm{Dt}}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}[m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{zz}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{zz}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{zx}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{zx}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{zy}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{zy}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}]$

其中，

分别为粒子p_i的运动速度在各坐标方向的分量；

为粒子p_i应力向量的正应力分量，

为粒子p_i应力向量的剪应力分量；

为粒子p_j应力向量的正应力分量，

为粒子p_j应力向量的剪应力分量；ρ_i，ρ_j分别为粒子p_i和p_j的密度；m_j，W_ij，N的含义同公式(8)

构建几何方程的SPH格式

${\mathrm{\&epsiv;}}_i^{\mathrm{xx}}=\frac{{\&PartialD;u}_i}{\&PartialD;x}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;u_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}$

${\mathrm{\&epsiv;}}_i^{\mathrm{yy}}=\frac{{\&PartialD;v}_i}{\&PartialD;y}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;v_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}$

${\mathrm{\&epsiv;}}_i^{\mathrm{zz}}=\frac{{\&PartialD;w}_i}{\&PartialD;z}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;w_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}---\left(10\right)$

${\mathrm{\&epsiv;}}_i^{\mathrm{xy}}=\frac{{\&PartialD;u}_i}{\&PartialD;y}+\frac{{\&PartialD;v}_i}{\&PartialD;x}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;u_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}+\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;v_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}$

${\mathrm{\&epsiv;}}_i^{\mathrm{yz}}=\frac{{\&PartialD;v}_i}{\&PartialD;z}+\frac{{\&PartialD;w}_i}{\&PartialD;y}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;v_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}+\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;w_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}$

${\mathrm{\&epsiv;}}_i^{\mathrm{zx}}=\frac{{\&PartialD;w}_i}{\&PartialD;x}+\frac{{\&PartialD;u}_i}{\&PartialD;z}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;w_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}+\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;u_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}$

其中，

为粒子p_i应变向量的正应变分量，为粒子p_i应变向量的剪应变分量；u_i，v_i，w_i为粒子p_i的运动位移在各坐标方向的分量，u_j，v_j，w_j为粒子p_j的运动位移在各坐标方向的分量；ρ_j为粒子p_j的密度，m_j，W_ij，N的含义同式(8)。

根据以上各方程，具体描述基于SPH的粒子运动状态求解过程如下：

1)对模型的某一质点p施加外力f_outp；

2)对模型中的各粒子p_i(i＝1，2，..，M)循环完成下述3)～7)步骤的运算处理；

3)对当前粒子p_i，以h为光滑半径，搜索支持域内的相邻粒子p_j(j＝1，2，..，N)；

4)使用公式(6)(7)，计算当前粒子p_i与支持域内各邻近粒子p_j，(j＝1，2，..，N)间的光滑核函数W_ij及其导数；

5)使用方程(8)，采用密度求和法计算粒子的密度ρ_i；

6)计算粒子的加速度a_i，加速度的计算采用如下方法：

$a_i=-\frac{{\mathrm{Dv}}_i}{\mathrm{Dt}}+\frac{{\mathrm{\&sigma;}}_i}{m_i}+\frac{F_i}{m_i}$

其中，F_i表示粒子p_i所受的重力、外力等，σ_i表示粒子p_i所受的内力，即应力，m_i为粒子p_i的质量，的计算见式(9)。

7)计算粒子的速度和位移：

v_i(t)＝v_i(t)+dt·a_i

x_i(t)＝x_i(t)+dt·v_i(t)

其中，v_i(t)为粒子p_i在t时刻的速度向量，x_i(t)为粒子p_i在t时刻的位置向量，a_i为粒子p_i的加速度，dt为时间增量。

8)计算各粒子与初始状态时的位移情况disp_i(t)＝x_i(t)-x_i(t₀)；

9)使用几何方程(4)(10)，由上一步算出的位移计算各粒子的应变状态

10)记录各粒子当前的体积应变

和形状畸变

并使用方程(2)，由上一步计算出的应变计算各粒子的新的体积应变

和形状畸变

11)计算

并使用黏弹性应力-应变本构方程(3)，计算各粒子的体积应力

和偏应力

12)使用方程(1)，计算出各粒子的应力状态

4.输出并渲染软组织器官的运动状态

采用光滑粒子流体动力学法，可计算得到所有粒子在每个时间步长的密度、位置、速度等运动信息，利用OpenGL技术，将每步长时刻计算出的粒子按位置向量输出至屏幕，可得到粒子模型在外力作用下动态的形变过程；而且，将粒子模型经表面三角化处理后，经过纹理和光照渲染，可得到逼真的软组织器官动态形变效果。图7以肝脏器官为例，给出了分别在受到拉力和压力情况下组织的形变效果，其中的圆柱表示模拟手术器械。

Claims (2)

1.一种软组织形变仿真方法和技术，其特征在于包含以下步骤：

步骤1)：采集软组织的数据信息；

步骤2)：选择黏弹性模型，构建用于软组织形变仿真的方程组：

采用一个弹簧和一个黏壶并联的Kelvin黏弹性模型；

首先，构建三维格式的Kelvin黏弹性应力-应变本构方程：

依据Kelvin模型，各向同性材料的应力张量σ可分解成它的球形张量和偏斜张量部分，应变张量ε可分离为体积形变和等体积的形状畸变两部分：

${\mathrm{\&sigma;}}_{\mathrm{\&alpha;\&beta;}}=S_{\mathrm{\&alpha;\&beta;}}+\frac{{\mathrm{\&delta;}}_{\mathrm{\&alpha;\&beta;}}{\mathrm{\&sigma;}}_{\mathrm{kk}}}{3}---\left(1\right)$

${\mathrm{\&epsiv;}}_{\mathrm{\&alpha;\&beta;}}=e_{\mathrm{\&alpha;\&beta;}}+\frac{{\mathrm{\&delta;}}_{\mathrm{\&alpha;\&beta;}}{\mathrm{\&epsiv;}}_{\mathrm{kk}}}{3}---\left(2\right)$

其中，α，β＝x，y，z.σ_αβ为应力分量，ε_αβ为应变分量；δ_αβ为Kronecker符号，σ_kk＝σ_xx+σ_yy+σ_zz和ε_kk＝ε_xx+ε_yy+ε_zz分别为体积应力和体积应变；S_αβ和e_αβ分别为偏应力张量和偏应变张量的分量；

根据Kelvin模型，偏应力张量和偏应变张量之间、体积应力和体积应变之间的三维黏弹性本构关系可表示为：

$S_{\mathrm{\&alpha;\&beta;}}=E\&CenterDot;e_{\mathrm{\&alpha;\&beta;}}+\mathrm{\&eta;}\&CenterDot;\frac{{\mathrm{de}}_{\mathrm{\&alpha;\&beta;}}}{\mathrm{dt}},\mathrm{\&alpha;},\mathrm{\&beta;}=x,y,z$

${\mathrm{\&sigma;}}_{\mathrm{kk}}=E\&CenterDot;{\mathrm{\&epsiv;}}_{\mathrm{kk}}+\mathrm{\&eta;}\&CenterDot;\frac{{\mathrm{d\&epsiv;}}_{\mathrm{kk}}}{\mathrm{dt}}---\left(3\right)$

上式中，S_αβ和e_αβ分别为偏应力张量和偏应变张量的分量；σ_kk和ε_kk分别为体积应力和体积应变；E为材料的弹性模量，η为黏性系数；

分别为偏应变分量和体积应变对时间的导数，即应变率，t为时间；

其次，构建应变-位移几何方程：

${\mathrm{\&epsiv;}}_{\mathrm{xx}}=\frac{\&PartialD;u}{\&PartialD;x}$ ${\mathrm{\&epsiv;}}_{\mathrm{xy}}=\frac{\&PartialD;u}{\&PartialD;y}+\frac{\&PartialD;v}{\&PartialD;x}$

${\mathrm{\&epsiv;}}_{\mathrm{yy}}=\frac{\&PartialD;v}{\&PartialD;y}$ ${\mathrm{\&epsiv;}}_{\mathrm{yz}}=\frac{\&PartialD;v}{\&PartialD;z}+\frac{\&PartialD;w}{\&PartialD;y}---\left(4\right)$

${\mathrm{\&epsiv;}}_{\mathrm{zz}}=\frac{\&PartialD;w}{\&PartialD;z}$ ${\mathrm{\&epsiv;}}_{\mathrm{zx}}=\frac{\&PartialD;w}{\&PartialD;x}+\frac{\&PartialD;u}{\&PartialD;z}$

其中u，v，w为位移在三个坐标方向的分量；ε_xx，ε_yy，ε_zz为正应变分量；ε_xy，ε_yz，ε_zx为剪应变分量；

然后，构建用于加速度计算的动量方程，其形式如下：

$\frac{\mathrm{dv}}{\mathrm{dt}}=\frac{1}{\mathrm{\&rho;}}\frac{\&PartialD;\mathrm{\&sigma;}}{\&PartialD;x}---\left(5\right)$

其中，v为速度向量，t为时间，ρ为粒子的密度，x为坐标向量；

步骤3)：依据步骤1)采集的软组织的数据信息，利用各数据点的位置向量，构建一个没有网格连接的粒子模型，并初始化所述粒子模型中各粒子的位置、质量、速度、加速度、受作用力等信息，构建粒子模型的初始化状态；

步骤4)：定义空间栅格，采用链表搜索法搜索粒子的支持域，构建支持域内的光滑核函数；

核函数选取三次B-Spline光滑函数：

$W_{\mathrm{ij}}=\frac{1}{{\mathrm{\&pi;h}}^3}\&times;\left\{\begin{array}{cc}1-\frac{3}{2}{s}^2+\frac{3}{4}{s}^3& 0\&le;s<1\\ \frac{1}{4}{(2-s)}^3& 1\&le;s<2\\ 0& s\&GreaterEqual;2\end{array}\right.---\left(6\right)$

式中，W_ij为由邻近粒子p_j估计粒子p_i运动信息的光滑核函数，

为粒子p_i与p_j之间的相对距离，r为粒子p_i与p_j之间的距离，h为光滑长度；

对核函数在各方向求导，得到核函数的一阶导数为：

$\frac{{\&PartialD;W}_{\mathrm{ij}}}{{\&PartialD;x}^{\mathrm{\&beta;}}}=\frac{x_i^{\mathrm{\&beta;}}-x_j^{\mathrm{\&beta;}}}{r}\&CenterDot;\frac{1}{{\mathrm{\&pi;h}}^4}\&times;\left\{\begin{array}{cc}-3s+\frac{9}{4}{s}^2& 0\&le;s<1\\ -\frac{3}{4}{(2-s)}^2& 1\&le;s<2\\ 0& s\&GreaterEqual;2\end{array}\right.---\left(7\right)$

式中，上标β＝x，y，z表示坐标方向，s，r，h的含义同式(6)，和

分别表示粒子p_i与p_j的位置坐标向量在各方向的分量；

步骤5)：应用光滑核函数W及其导数对参考粒子的支持域内所有粒子函数的加权平均近似的方法来构建步骤2)中各方程的SPH格式：

构建密度方程SPH格式

${\mathrm{\&rho;}}_i=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{m}_j{W}_{\mathrm{ij}}---\left(8\right)$

其中，ρ_i为粒子p_i的密度，m_j为p_i的支持域内邻近粒子p_j的质量，N为粒子p_i的支持域内粒子总数，W_ij的含义同式(6)；

将动量方程(5)转换为SPH格式

$\frac{{\mathrm{Dv}}_i^x}{\mathrm{Dt}}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}[m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{xx}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{xx}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{xy}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{xy}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{xz}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{xz}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}]$

$\frac{{\mathrm{Dv}}_i^y}{\mathrm{Dt}}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}[m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{yy}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{yy}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{yx}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{yx}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{yz}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{yz}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}]---\left(9\right)$

$\frac{{\mathrm{Dv}}_i^z}{\mathrm{Dt}}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}[m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{zz}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{zz}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{zx}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{zx}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}+m_j(\frac{{\mathrm{\&sigma;}}_i^{\mathrm{zy}}}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&sigma;}}_j^{\mathrm{zy}}}{{\mathrm{\&rho;}}_j^2})\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}]$

其中，

分别为粒子p_i的运动速度在各坐标方向的分量；为粒子p_i应力向量的正应力分量，

为粒子p_i应力向量的剪应力分量；

为粒子p_j应力向量的正应力分量，为粒子p_j应力向量的剪应力分量；ρ_i，ρ_j分别为粒子p_i和p_j的密度；m_j，W_ij，N的含义同公式(8)；

构建几何方程的SPH格式：

${\mathrm{\&epsiv;}}_i^{\mathrm{xx}}=\frac{{\&PartialD;u}_i}{\&PartialD;x}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;u_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}$

${\mathrm{\&epsiv;}}_i^{\mathrm{yy}}=\frac{{\&PartialD;v}_i}{\&PartialD;y}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;v_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}$

${\mathrm{\&epsiv;}}_i^{\mathrm{zz}}=\frac{{\&PartialD;w}_i}{\&PartialD;z}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;w_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}---\left(10\right)$

${\mathrm{\&epsiv;}}_i^{\mathrm{xy}}=\frac{{\&PartialD;u}_i}{\&PartialD;y}+\frac{{\&PartialD;v}_i}{\&PartialD;x}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;u_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}+\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;v_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}$

${\mathrm{\&epsiv;}}_i^{\mathrm{yz}}=\frac{{\&PartialD;v}_i}{\&PartialD;z}+\frac{{\&PartialD;w}_i}{\&PartialD;y}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;v_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}+\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;w_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;y}$

${\mathrm{\&epsiv;}}_i^{\mathrm{zx}}=\frac{{\&PartialD;w}_i}{\&PartialD;x}+\frac{{\&PartialD;u}_i}{\&PartialD;z}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;w_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;x}+\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{m_j}{{\mathrm{\&rho;}}_j}\&CenterDot;u_j\&CenterDot;\frac{{\&PartialD;W}_{\mathrm{ij}}}{\&PartialD;z}$

其中，

为粒子p_i应变向量的正应变分量，

为粒子p_i应变向量的剪应变分量；u_i，v_i，w_i为粒子p_i的运动位移在各坐标方向的分量，u_j，v_j，w_j为粒子p_j的运动位移在各坐标方向的分量；ρ_j为粒子p_j的密度，m_j，W_ij，N的含义同式(8)；

步骤6)：用显示积分法求解步骤5)的常微分方程，计算出粒子的密度、位置、速度等随时间的变化值；

步骤7)：循环执行上述步骤3)～步骤6)，计算出各粒子的状态；

步骤8)：将粒子模型的当前状态输出之屏幕，经过纹理和光照渲染后，得到软组织的动态形变过程。

2.根据权利要求1的软组织形变仿真方法，其特征在于步骤6)的求解过程如下：

1)对模型的某一质点p施加外力f_outp；

2)对模型中的各粒子p_i(i＝1，2，..，M)循环完成下述3)～7)步骤的运算处理；

3)对当前粒子p_i，以h为光滑半径，搜索支持域内的相邻粒子p_j(j＝1，2，..，N)；

4)使用公式(6)(7)，计算当前粒子p_i与支持域内各邻近粒子p_j(j＝1，2，..，N)间的光滑核函数W_ij及其导数；

5)使用方程(8)，采用密度求和法计算粒子的密度ρ_i；

6)计算粒子的加速度a_i，加速度的计算采用如下方法：

$a_i=-\frac{{\mathrm{Dv}}_i}{\mathrm{Dt}}+\frac{{\mathrm{\&sigma;}}_i}{m_i}+\frac{F_i}{m_i}$

其中，F_i表示粒子p_i所受的重力、外力等，σ_i表示粒子p_i所受的内力，即应力，m_i为粒子p_i的质量，使用公式(9)计算

7)计算粒子的速度和位移：

v_i(t)＝v_i(t)+dt·a_i

x_i(t)＝x_i(t)+dt·v_i(t)

其中，v_i(t)为粒子p_i在t时刻的速度向量，x_i(t)为粒子p_i在t时刻的位置向量，a_i为粒子p_i的加速度，dt为时间增量；

8)计算各粒子与初始状态时的位移情况disp_i(t)＝x_i(t)-x_i(t₀)；

9)使用几何方程(4)和(10)，由上一步算出的位移计算各粒子的应变状态

10)记录各粒子当前的体积应变

和形状畸变并使用方程(2)，由上一步计算出的应变计算各粒子的新的体积应变

和形状畸变

11)计算

并使用黏弹性应力-应变本构方程(3)，计算各粒子的体积应力和偏应力

12)使用方程(1)，计算出各粒子的应力状态

Images