CN108536936A - 一种多重优化的无网格软组织形变模拟方法 - Google Patents

Abstract

本发明公开了一种多重优化的无网格软组织形变模拟方法，包括：(1)获取软组织CT图像的顶点信息，得到点云数据；(2)采用基于八叉树编码的统一流线法将所述点云数据进行简化；(3)采用移动最小二乘法对简化后的点云数据构建无网格形函数；(4)采用Kelvin粘弹性模型构造非线性粘弹性模型；(5)根据所述非线性粘弹性模型和无网格形函数采用EFG法构建非线性粘弹性模型的无网格方程；(6)基于高分辨率近似算法在每个时间片采用所述非线性粘弹性模型的无网格方程计算每个点的新位置，实现软组织形变模拟。本发明成本低，模拟效果好。

Description

一种多重优化的无网格软组织形变模拟方法

技术领域

本发明涉及软组织模拟方法，尤其涉及一种多重优化的无网格软组织形变模拟方法。

背景技术

随着科学技术的不断进步，虚拟手术模拟系统得到了前所未有的发展，这些系统对于训练医院外科手术实习生有着十分重要的意义，一方面能够迅速提升实习生的手术技能，并降低培训成本，另一方面能够提高手术成功率，也是对病人安全方面的负责。在医学领域中真实的软组织在形变的过程中会表现出复杂的生物力学特性，这给软组织建模带来了巨大困难。

因此，近年来关于软组织形变的模拟成为虚拟物体仿真研究的焦点。目前已经提出了许多变形模型，常见的有：(1)弹簧-质量模型由于其简单的结构和快速的计算常常被用于手术模拟，但难以设置合适的参数值，模型不稳定；(2)有限元模型精度高但是计算复杂度高；(3)无网格模型不需要维护数据点之间的拓扑信息，避免了网格模型中固有复杂拓扑结构(如：病态网格，网格重构等)的相关问题，具有较强的适应性，但是现有与无网格模型结合的模型还有很多缺点，如：线性粘弹性无网格模型，虽然线性弹性模型可以减少计算问题但是同时也使许多组织生物力学特性丧失，其他与无网格相关的一些模型虽然能够较为真实的模拟软组织，但是没有对对象的数据点进行很好的优化和存储，计算成本较高。因此开发一个计算成本低，模拟效果好的模型至关重要。

发明内容

发明目的：本发明针对现有技术存在的问题，提供一种多重优化的无网格软组织形变模拟方法，成本低，模拟效果好。

技术方案：本发明所述的一种多重优化的无网格软组织形变模拟方法包括：

(1)获取软组织CT图像的顶点信息，得到点云数据；

(2)采用基于八叉树编码的统一流线法将所述点云数据进行简化；

(3)采用移动最小二乘法对简化后的点云数据构建无网格形函数；

(4)采用Kelvin粘弹性模型构造非线性粘弹性模型；

(5)根据所述非线性粘弹性模型和无网格形函数采用EFG法构建非线性粘弹性模型的无网格方程；

(6)基于高分辨率近似算法在每个时间片采用所述非线性粘弹性模型的无网格方程计算每个点的新位置，实现软组织形变模拟。

进一步的，步骤(1)具体包括：

获取软组织的CT图像，并导入到软件Mimics中，导出三维模型的STL文件，使用MeshLab软件将STL文件转换成OBJ文件，从而获取顶点信息，得到点云数据。

进一步的，步骤(2)具体包括：

(2-1)将点云数据的空间划分为固定边长的多个立方体，保留每个立方体中心最近点；

(2-2)计算每个点的八叉树编码，按编码值从小到大恢复点，将具有相同编码值的点存储在同一链表中，同一链表中的多个点只保留离中心最近的点，从而实现简化。

进一步的，步骤(3)中移动最小二乘法形函数具体为：

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

其中，Φ_i(x)表示点x_i的形函数，

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

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

P^T(x)是多项式基函数矩阵P(x)的转置，A^-1(x)是加权瞬时矩阵A(x)的逆形式，B(x)是加权瞬时矩阵，W(x)是加权函数矩阵，p_j(x_i)是多项式基函数，j＝1,2…,m，i＝1,2…,n，m为P(x)列数，n为P(x)行数，x_i表示Θ中第i个点，Θ为从简化后的点云数据中选取的点数量大于预设值的支持域，为矩形或者圆形。

进一步的，步骤(4)中的非线性粘弹性模型具体包括：

a.本构方程：

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

b.应力与应变的本构关系：

式中，σ(t)表示材料的松弛响应，E(t)表示松弛模量，ε表示应变，c₀和c₁表示材料参数，t表示时间变量，τ₁表示时间常数；

c.软组织在应力加载过程中的松弛本构关系：

式中，调整函数N(ε)＝ε+kε²，τ表示阶跃应力加载时间，k是常系数，是应变参数；

d.应力增量：Δσ_n＝E_kΔε_n+σ_0,n；

式中，Δσ_n表示从t_n时刻到t_n+1时刻的应力增量，E_k为非线性松弛系数，Δε_n表示应变增量，初始应力

e.应变与应力关系：ε_n+1＝Lu_n+1+LΔu_n；

式中，ε_n+1、u_n+1分别表示t_n+1时刻的应变、位移，L表示偏微分算子，Δu_n表示从t_n时刻到t_n+1时刻的位移增量。

进一步的，步骤(5)中所述非线性粘弹性模型的无网格方程具体为：

式中，K_n表示粘弹性刚度矩阵，表示由EFG中形函数以及形函数导数确定的惩罚刚度矩阵，α是惩罚因子，ΔR_n表示不平衡力矢量，Δu_n表示位移增量，求解该方程得到软组织的非线性和粘弹性。

进一步的，步骤(6)具体包括：

(6-1)对简化后施加力的区域的点进行采样；

(6-2)计算采样点中每个点的l个表面最近邻点，得到所有非采样表面点中d_i的l个最近邻点的集合N_u(d_i)以及包括表面和内部采样点的所有体积数据中d_i的l个最近邻居的集合N_s(d_i)，其中，d_i∈D,D表示施加力区域的点集合；

(6-3)在每个时间步，利用所述非线性粘弹性模型的无网格方程计算每个点的新位置，然后根据N_s(d_i)中每个点的初始位置和当前位置得到d_i的仿射变换矩阵；

(6-4)根据仿射变换矩阵计算出属于N_u(d_i)的每个点的新位置并进行渲染，实现软组织形变模拟。

有益效果：本发明与现有技术相比，其显著优点是：经过多重优化的无网格软组织形变模型具有计算快，模拟效果逼真的优势，具体如下：(1)采用基于八叉树编码的统一流线法简化点，缩小了存储空间，并用八叉树对点进行存储，大大缩短了点的搜索时间；(2)在组织模型中纳入了Kelvin粘弹性力学模型，可更好表现软组织的粘弹性，提高模拟真实度；(3)在粘弹性无网格模型中加入调整函数，展现应力与应变之间非线性关系，可有效的模拟软组织的真实属性；(4)运用高分辨近似算法进行优化，渲染效果好且降低了计算成本。

附图说明

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

图2是Kelvin粘弹性模型图；

图3是八叉树存储结构定义图；

图4是表面点的邻点集合示意图。

具体实施方式

本实施例提供了一种多重优化的无网格软组织形变模拟方法，如图1包括：

(1)获取软组织CT图像的顶点信息，得到点云数据。具体包括：

获取软组织的CT图像，并导入到软件Mimics中，导出三维模型的STL文件，使用MeshLab软件将STL文件转换成OBJ文件，从而获取顶点信息，得到点云数据。

(2)采用基于八叉树编码的统一流线法将所述点云数据进行简化。

由MeshLab导出的几何模型包含大量的数据，故采用基于八叉树编码的统一流线法来简化点。具体包括：

(2-1)将点云数据的空间划分为固定边长的多个立方体，保留每个立方体中心最近点；

(2-2)计算每个点的八叉树编码，按编码值从小到大恢复点，将具有相同编码值的点存储在同一链表中，同一链表中的多个点只保留离中心最近的点，从而实现简化。

(3)采用移动最小二乘法对简化后的点云数据构建无网格形函数。

移动最小二乘法收敛速度快，故使用该方法来构建无网格形函数，同时在无网格模型中计算点插值的准确度依赖于支持域中点的数量，故需要选择合适的近似支持域Θ，一般为矩形或者圆形。移动最小二乘法形函数具体为：

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

其中，Φ_i(x)表示点x_i的形函数，

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

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

P^T(x)是多项式基函数矩阵P(x)的转置，A^-1(x)是加权瞬时矩阵A(x)的逆形式，B(x)是加权瞬时矩阵，W(x)是加权函数矩阵，p_j(x_i)是多项式基函数，j＝1,2…,m，i＝1,2…,n，m为P(x)列数，n为P(x)行数，是场函数，u_i是在点x_i处的节点场值，x_i表示Θ中第i个点，Θ为从简化后的点云数据中选取的点数量大于预设值的支持域。

(4)采用Kelvin粘弹性模型构造非线性粘弹性模型。

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

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

应力与应变的本构关系为：

式中，σ(t)表示材料的松弛响应，E(t)表示松弛模量，ε表示应变，c₀和c₁表示材料参数，t表示时间变量，τ₁表示时间常数。

在Kelvin粘弹性模型中软组织在应力加载过程中的松弛本构关系可表示为：

为了显示应力和应变之间的非线性系，在粘弹性松弛本构关系上增加调整函数N(ε)：

式中，N(ε)＝ε+kε²，τ表示阶跃应力加载时间，k是常系数，是应变参数。

将变形模拟时间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 (9)

式中，σ_n表示在t_n时刻的应力，σ_n+1表示在t_n+1时刻的应力，σ_0,n是初始应力，E_k为非线性松弛系数。

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

u_n+1＝u_n+Δu_n (10)

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

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

应变与应力关系为：

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

式中，ε_n+1、u_n+1分别表示t_n+1时刻的应变、位移，L表示偏微分算子。

(5)根据所述非线性粘弹性模型采用EFG法构建非线性粘弹性模型的无网格方程。

由于移动最小二乘法不显示Kroneckerδ函数的性质，所以引用了EFG的弱形式：

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

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

式中，K_n表示粘弹性刚度矩阵，表示由EFG中形函数以及形函数导数确定的惩罚刚度矩阵，Δu_n表示位移增量，ΔR_n表示不平衡力矢量，求解该方程得到软组织的非线性和粘弹性。

(6)基于高分辨率近似算法在每个时间片采用所述非线性粘弹性模型的无网格方程计算每个点的新位置，实现软组织形变模拟。

采用高分辨近似算法进行优化，该方法是基于最佳逼近单个时间步长的采样点以及相邻点的运动的原理。具体包括：

(6-1)对简化后施加力的区域的点进行采样；

(6-2)计算采样点中每个点的l个表面最近邻点，得到所有非采样表面点中d_i的l个最近邻点的集合N_u(d_i)以及包括表面和内部采样点的所有体积数据中d_i的l个最近邻居的集合N_s(d_i)，其中，d_i∈D,D表示施加力区域的点集合，如图4所示；

(6-3)在每个时间步，利用所述非线性粘弹性模型的无网格方程计算每个点的新位置，然后根据N_s(d_i)中每个点的初始位置和当前位置得到d_i的仿射变换矩阵；

(6-4)根据仿射变换矩阵计算出属于N_u(d_i)的每个点的新位置并进行渲染，实现软组织形变模拟。

以上所揭露的仅为本发明一种较佳实施例而已，不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。

Claims (7)

1.一种多重优化的无网格软组织形变模拟方法，其特征在于该方法包括：

(1)获取软组织CT图像的顶点信息，得到点云数据；

(2)采用基于八叉树编码的统一流线法将所述点云数据进行简化；

(3)采用移动最小二乘法对简化后的点云数据构建软无网格形函数；

(4)采用Kelvin粘弹性模型构造非线性粘弹性模型；

(5)根据所述非线性粘弹性模型和无网格形函数采用EFG法构建非线性粘弹性模型的无网格方程；

(6)基于高分辨率近似算法在每个时间片采用所述非线性粘弹性模型的无网格方程计算每个点的新位置，实现软组织形变模拟。

2.根据权利要求1所述的一种多重优化的无网格软组织形变模拟方法，其特征在于：步骤(1)具体包括：

获取软组织的CT图像，并导入到软件Mimics中，导出三维模型的STL文件，使用MeshLab软件将STL文件转换成OBJ文件，从而获取顶点信息，得到点云数据。

3.根据权利要求1所述的一种多重优化的无网格软组织形变模拟方法，其特征在于：步骤(2)具体包括：

(2-1)将点云数据的空间划分为固定边长的多个立方体，保留每个立方体中心最近点；

(2-2)计算每个点的八叉树编码，按编码值从小到大恢复点，将具有相同编码值的点存储在同一链表中，同一链表中的多个点只保留离中心最近的点，从而实现简化。

4.根据权利要求1所述的一种多重优化的无网格软组织形变模拟方法，其特征在于：步骤(3)中的移动最小二乘法形函数具体为：

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

其中，Φ_i(x)表示点x_i的形函数，

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

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

P^T(x)是多项式基函数矩阵P(x)的转置，A^-1(x)是加权瞬时矩阵A(x)的逆形式，B(x)是加权瞬时矩阵，W(x)是加权函数矩阵，p_j(x_i)是多项式基函数，j＝1,2…,m，i＝1,2…,n，m为P(x)列数，n为P(x)行数，x_i表示Θ中第i个点，Θ为从简化后的点云数据中选取的点数量大于预设值的支持域，为矩形或者圆形。

5.根据权利要求1所述的一种多重优化的无网格软组织形变模拟方法，其特征在于：步骤(4)中的非线性粘弹性模型具体包括：

a.本构方程：

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

b.应力与应变的本构关系：

式中，σ(t)表示材料的松弛响应，E(t)表示松弛模量，ε表示应变，c₀和c₁表示材料参数，t表示时间变量，τ₁表示时间常数；

c.软组织在应力加载过程中的松弛本构关系：

式中，调整函数N(ε)＝ε+kε²，τ表示阶跃应力加载时间，k是常系数，是应变参数；

d.应力增量：Δσ_n＝E_kΔε_n+σ_0,n；

式中，Δσ_n表示从t_n时刻到t_n+1时刻的应力增量，E_k为非线性松弛系数，Δε_n表示应变增量，初始应力

e.应变与应力关系：ε_n+1＝Lu_n+1+LΔu_n；

式中，ε_n+1、u_n+1分别表示t_n+1时刻的应变、位移，L表示偏微分算子，Δu_n表示从t_n时刻到t_n+1时刻的位移增量。

6.根据权利要求1所述的一种多重优化的无网格软组织形变模拟方法，其特征在于：步骤(5)中所述非线性粘弹性模型的无网格方程具体为：

式中，K_n表示粘弹性刚度矩阵，表示由EFG中形函数以及形函数导数确定的惩罚刚度矩阵，α是惩罚因子，Δu_n表示位移增量，ΔR_n表示不平衡力矢量，求解该方程得到软组织的非线性和粘弹性。

7.根据权利要求1所述的一种多重优化的无网格软组织形变模拟方法，其特征在于：步骤(6)具体包括：

(6-1)对简化后施加力的区域的点进行采样；

(6-2)计算采样点中每个点的l个表面最近邻点，得到所有非采样表面点中d_i的l个最近邻点的集合N_u(d_i)以及包括表面和内部采样点的所有体积数据中d_i的l个最近邻居的集合N_s(d_i)，其中，d_i∈D,D表示施加力区域的点集合；

(6-3)在每个时间步，利用所述非线性粘弹性模型的无网格方程计算每个点的新位置，然后根据N_s(d_i)中每个点的初始位置和当前位置得到d_i的仿射变换矩阵；

(6-4)根据仿射变换矩阵计算出属于N_u(d_i)的每个点的新位置并进行渲染，实现软组织形变模拟。