CN108710735A - 一种实时交互的无网格软组织形变模拟方法 - Google Patents

Abstract

本发明公开了一种实时交互的无网格软组织形变模拟方法，包括以下步骤：(1)将软组织的CT图像转换成OBJ文件以便获取顶点信息；(2)使用移动最小二乘法MLS来构建无网格形函数；(3)将Kelvin粘弹性纳入组织模型；(4)根据步骤(3)中计算的给定应力相对应的所有节点的位移，使用Marquardt算法来预先拟合软组织上的力和每个节点的位移之间的数学关系；(5)引入富集函数来模拟因切割等操作导致的组织表面不连续性裂缝。本发明在人体组织模型中纳入了Kelvin粘弹性力学模型，可更好表现软组织的粘弹性，提高模拟真实度，使用Marquardt算法预先拟合软组织上的力和每个节点的位移之间的数学关系，提高模拟实时性，同时引入权重函数富集来实现软组织的交互式模拟。

Description

一种实时交互的无网格软组织形变模拟方法

技术领域

本发明涉及软组织变形模拟方法，尤其涉及一种实时交互的无网格软组织形变模拟方法。

背景技术

拜师学艺是中华儿女从古至今学习技能的重要方式，医术学习也不例外，新手通过向有经验技术的老医者学习，提高自身的医学技术水平。尽管传统的收徒授艺为医学教育发展奠定了基础，但是这种培训方式周期长，代价大。随着计算机科学技术的发展，基于计算机的虚拟手术模拟逐渐解决了上述缺陷。

在虚拟手术模拟中人体软组织建模是其中的一项重要研究领域，人体软组织可以通过模型准确获得在外力作用下发生的形变，这有助于实现人与模型之间的实时交互。目前已经提出了许多模型来模拟人体软组织，常见的有：(1)弹簧质点模型，由于其结构简单和计算速度快常常被用于手术模拟，但该模型受到许多约束，迭代计算不稳定；(2)有限元模型，该模型精度高、适应性强，但是计算量大、复杂性高，并且很难实现实时模拟；(3)无网格模型，该模型只需要一组离散节点的信息，不需要预先处理网格数据进而不会出现网格失真或者纠缠等问题，并且相对于有限元模型而言更能够实现软组织实时模拟变形，但是有时模拟不够真实。因此开发一个具有良好实时性和真实性的人体软组织模型对于虚拟手术系统来说至关重要。

发明内容

发明目的：本发明目的是着力解决软组织形变模拟中实时性不强，真实性不足的问题，提出了一种实时交互的无网格软组织形变模拟方法。

技术方案：本发明包括以下步骤：

(1)将软组织的CT图像导入到Mimics软件中，并导出三维模型的STL文件，使用MeshLab软件将STL文件转换成OBJ文件以便获取顶点信息；

(2)使用移动最小二乘法MLS来构建无网格形函数；

(3)将Kelvin粘弹性纳入组织模型；

(4)根据步骤(3)中计算的给定应力相对应的所有节点的位移使用Marquardt算法来预先拟合软组织上的力和每个节点的位移之间的数学关系；

(5)引入富集函数来模拟因切割等操作导致的组织表面不连续性裂缝。

所述步骤(2)中的无网格形函数可以表示为：

Φ^T(x)＝(Φ₁(x),Φ₂(x),......Φ_n(x))＝P^T(x)A^-1(x)B(x) (1)

其中，Φ^T(x)表示由MLS构建的无网格形函数，Φ_i(x)(i＝1,2,…,n)表示点x_i的形函数，P^T(x)是多项式基函数矩阵P(x)的转置，A^-1(x)是第一加权瞬时矩阵A(x)的逆形式，B(x)是第二加权瞬时矩阵。

所述步骤(3)中的Kelvin模型可以表示为：

其中，σ、η、K₂、K₁、ε和ε_d分别表示应力、阻尼器、模型中第二个弹簧的刚度、应力的时间导数、模型中第一个弹簧的刚度、应变、应变的时间导数。

所述步骤(4)中的力在x,y,z三个方向的分力和每个节点的位移之间的数学关系为：

其中是对应三个分量的系统参数，λ_x，λ_y，λ_z是由交叉验证确定的系统自变量。

所述步骤(5)中的富集函数可以表示为：

其中，h(x,y)表示从0到1平滑过渡的富集函数，表示二维距离函数d₂(x,y)对法线坐标s的偏导数。

有益效果：本发明具有良好的实时性和真实性，在人体组织模型中纳入了Kelvin粘弹性力学模型，可更好表现软组织的粘弹性，提高模拟真实度，使用Marquardt算法预先拟合软组织上的力和每个节点的位移之间的数学关系，提高模拟实时性，同时引入权重函数富集来模拟因切割等操作导致的组织表面不连续性裂缝，实现软组织的交互式模拟。

附图说明

图1为本发明的流程示意图；

图2为本发明的Kelvin粘弹性模型示意图；

图3为本发明的具有不连续性裂缝的软组织分析域示意图。

具体实施方式

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

如图1所示，本发明包括以下步骤：

(1)将软组织的CT图像导入到软件Mimics中并导出三维模型的STL文件，使用MeshLab软件将STL文件转换成OBJ文件以便获取顶点信息。

(2)由于移动最小二乘法(MLS)收敛速度快，故使用MLS来构建无网格形函数，具体如下：

Φ^T(x)＝(Φ₁(x),Φ₂(x),......Φ_n(x))＝P^T(x)A^-1(x)B(x) (1)

其中，Φ^T(x)表示由MLS构建的无网格形函数，Φ_i(x)(i＝1,2,…,n)表示点x_i的形函数，P^T(x)是多项式基函数矩阵P(x)的转置，A^-1(x)是第一加权瞬时矩阵A(x)的逆形式，B(x)是第二加权瞬时矩阵，具体如下：

A(x)＝P^T(x)W(x)P(x) (2)

B(x)＝P^T(x)W(x) (3)

w_i(x)＝w(x-x_i) (5)

其中，W(x)是加权函数矩阵，w_i(x)(i＝1,2,…,n)表示第i个节点的权重函数，

其中，P(x)是多项式基函数矩阵,p_j(x_i)是多项式基函数，j＝1,2,…,m，i＝1,2,…,n，m为P(x)列数，n为P(x)行数，f^h(x)是场函数，f_i是在节点x_i处的节点场值。

(3)为了提高模拟的真实度，将Kelvin粘弹性纳入组织模型，Kelvin模型是一个标准的线性模型，模型中的弹簧表示软组织的线弹性特点，阻尼器表示软组织结构变化时的阻尼特性，具体如下：

Kelvin模型的本构方程如下：

其中，σ、η、K₂、K₁、ε和ε_d分别表示应力、阻尼器、模型中第二个弹簧的刚度、应力的时间导数、模型中第一个弹簧的刚度、应变、应变的时间导数。在Kelvin粘弹性模型中，软组织在应力加载过程中的松弛本构关系可表示为：

其中，表示材料的松弛响应，t表示时间变量，E(t-τ)表示松弛模量，τ表示阶跃应力加载时间。将变形模拟时间T划分为n个时间片t₁,t₂,…,t_n,Δt＝T/n称为一个时间增量，从t_n到t_n+1的时间段，位移，应力，应变的增量分别为Δu_n，Δσ_n，Δε_n。Kelvin模型在时间t_n到t_n+1的应力增量可表示如下：

Δσ_n＝E_kΔε_n+σ_0,n (10)

其中，E_k是非线性松弛系数，σ_0，n是初始应力。

在t_n+1时，位移，应力，应变分别增变为：

u_n+1＝u_n+Δu_n (11)

σ_n+1＝σ_n+Δσ_n (12)

ε_n+1＝ε_n+Δε_n (13)

同时，应变与应力的关系可表示为：

ε_n+1＝Lu_n+1+LΔu_n (14)

其中，L表示偏微分算子。由于MLS不显示Kroneckerδ函数的性质，所以引用了EFG的弱形式：

其中，Ω表示分析域，δ表示Kronecker函数值，u表示位移，D表示弹性常数矩阵，u^T表示位移矢量，b表示物理矢量，Γ_t表示自然边界条件，是对应的给予表面的力，Γ_u表示本质边界条件，表示对应于本质边界条件的位移，α是惩罚因子。

将公式(10)-(14)代入公式(15)可得到粘弹性无网格求解方程的增量形式：

其中，K_n表示粘弹性刚度矩阵，表示由EFG中形函数以及形函数导数确定的惩罚刚度矩阵,ΔR_n表示不平衡力矢量。给定相应的材料参数和时间步长，可以根据公式(16)计算出Δu_n，最后计算每个点新的位移，应力和应变。

(4)为了提高模拟虚拟手术过程中变形的实时性，根据步骤(3)中计算的给定应力相对应所有节点的位移，使用Marquardt算法来预先拟合软组织上的力和每个节点的位移之间的数学关系，当软组织受到压力时，通过调用该拟合关系生成拟合曲面可以立即显示它们的变形，具体如下：

为了表示施加力与变形表面之间的关系，建立了施加力σ与受力节点ο之间的拟合关系，假设有n个力σ被施加到节点上，将σ在三个方向上分解为σ_x、σ_y和σ_z，并且在这三个方向上受力节点ο处相应的位移分量为Δu_x，Δu_y，Δu_z，各个分力相应的节点位移函数为：

其中是对应三个分量的系统参数，λ_x，λ_y，λ_z是由交叉验证确定的系统自变量。

当一个节点受到力的作用时，其他点也会改变位移，导致形变发生，因此，建立受力节点位移与其他节点位移之间的关系来表示变形表面，在分力σ_z的作用下，为了简化，假设表面节点具有相同的z坐标，诱导表面函数为：

其中，Δx，Δy，Δz表示在分力σ_z作用下其他点在三个方向上的诱导位移， 表示三个表面变形函数参数，a₁，a₂，a₃，b₁，b₂，b₃是通过交叉验证确定的表面变形顺序，同理可以求得分力σ_x和σ_y的诱导表面函数。

(5)人体软组织用于例如切割等交互式手术模拟是必不可少的，故引入权重函数富集来模拟因切割等操作导致的组织表面不连续性裂缝，其关键思想是用富集函数乘以相应节点的权重函数来改变其形函数。为了降低计算成本，只对节点距离小于它们的膨胀参数ρ_i的节点重新计算形函数，其他节点形函数保持不变。具体如下：

其中，h(x,y)是从0到1平滑过渡的富集函数，表示二维距离函数d₂(x,y)对法线坐标s的偏导数，表示有符号距离函数的正部分，ζ是切割局部坐标系的局部坐标，d_s(ζ)是具有端点(ζ₁,ζ₂)的一维带符号距离函数。

Claims (5)

1.一种实时交互的无网格软组织形变模拟方法，其特征在于，包括以下步骤：

(1)将软组织的CT图像导入到Mimics软件中，并导出三维模型的STL文件，使用MeshLab软件将STL文件转换成OBJ文件以便获取顶点信息；

(2)使用移动最小二乘法MLS来构建无网格形函数；

(3)将Kelvin粘弹性纳入组织模型；

(4)根据步骤(3)中计算的给定应力相对应的所有节点的位移，使用Marquardt算法来预先拟合软组织上的力和每个节点的位移之间的数学关系；

(5)引入富集函数来模拟因切割等操作导致的组织表面不连续性裂缝。

2.根据权利要求1所述的一种实时交互的无网格软组织形变模拟方法，其特征在于，所述步骤(2)中的无网格形函数可以表示为：

Φ^T(x)＝(Φ₁(x),Φ₂(x),......Φ_n(x))＝P^T(x)A^-1(x)B(x) (1)

其中，Φ^T(x)表示由MLS构建的无网格形函数，Φ_i(x)(i＝1,2,…,n)表示点x_i的形函数，P^T(x)是多项式基函数矩阵P(x)的转置，A^-1(x)是第一加权瞬时矩阵A(x)的逆形式，B(x)是第二加权瞬时矩阵。

3.根据权利要求1所述的一种实时交互的无网格软组织形变模拟方法，其特征在于，所述步骤(3)中的Kelvin模型可以表示为：

其中，σ、η、K₂、K₁、ε和ε_d分别表示应力、阻尼器、模型中第二个弹簧的刚度、应力的时间导数、模型中第一个弹簧的刚度、应变、应变的时间导数。

4.根据权利要求1所述的一种实时交互的无网格软组织形变模拟方法，其特征在于，所述步骤(4)中的力在x,y,z三个方向的分力和每个节点的位移之间的数学关系为：

其中是对应三个分量的系统参数，λ_x，λ_y，λ_z是由交叉验证确定的系统自变量。

5.根据权利要求1所述的一种实时交互的无网格软组织形变模拟方法，其特征在于，所述步骤(5)中的富集函数可以表示为：

其中，h(x,y)表示从0到1平滑过渡的富集函数，表示二维距离函数d₂(x,y)对法线坐标s的偏导数。