CN110289104B - 软组织按压和形变恢复的模拟方法 - Google Patents

Abstract

本发明公开了一种软组织按压和形变恢复的模拟方法，该方法用于模拟虚拟手术中对肝脏进行的按压排查过程，使用添加弯曲弹簧的质点弹簧模型对肝脏表面构建面模型，并采用粒子群算法对模型参数进行优化，然后使用基于层状的内部质点域约束对肝脏内部构建体模型，并在层状质点域之间采用约束弹簧元件进行连接。本发明方法不仅提高了质点弹簧模型在产生大形变时的形变恢复能力和计算效率，增强了其形变准确性和真实性，而且解决了传统质点弹簧模型无法模拟内部体积力的难点，并模拟了软组织形变过程的粘弹特性。

Description

软组织按压和形变恢复的模拟方法

技术领域

本发明涉及一种软组织力触觉模拟方法，尤其涉及一种用于模拟软组织按压和形变恢复的模拟方法。

背景技术

在对人体软组织，例如肝脏，实施手术的过程中，首先需要对肝脏进行按压感知以确定病变区域。现有技术中，提到了通过虚拟手术训练系统来为外科医生提供手术训练学习的条件与平台，实现虚拟手术系统的一个重要技术难点就是对虚拟肝脏或者其他软组织的高精度模拟。

在现有的虚拟软组织模拟方法中，有限元模型将问题求解域划分为若干个互不重叠的单元，通过单元形函数和节点插值函数建立刚度矩阵进而构建描述软组织材料力学特性的近似函数模型，由函数模型获得各单元内节点位移就可用来表现软组织的形变，因此其形变仿真精度最高，但在形变过程中网格拓扑结构会不断重组，从而导致计算复杂度太高和计算量太大，无法应用在实时性较强的虚拟手术中。为了降低计算的复杂度提高计算效率，提出了质点弹簧模型来模拟软组织的形变。因为该模型结构简单、计算复杂度低、易于实现并且能满足仿真实时性要求，但是传统的质点弹簧模型稳定性差，形变精确度低，且由于模型参数与软组织的物理性质没有直接联系，因此很难在保持模型真实性的情况下确定参数的大小，而且此模型在撤销外力后的形变恢复过程与实际软组织的恢复情况无法很好的吻合，恢复过程的模拟精度较压缩过程更差，甚至在产生大形变时恢复速度过于缓慢，与实际的软组织的恢复能力大相径庭。在一般情况下，通常使用质点弹簧面模型进行形变仿真，因而其无法模拟软组织内部体特性。可以看到，现有的简单质点-弹簧模型的模拟效果并不理想，尤其在形变恢复的模拟过程中达不到精度的要求。

发明内容

发明目的：针对以上问题，本发明提出一种用于虚拟手术系统的软组织模拟方法，模拟软组织在外力按压作用下的形变以及撤销外力后的形变恢复情况，该方法既能满足仿真实时性要求，又能同时提高对软组织形变以及形变恢复能力的模拟精确度，解决了撤销外力后软组织形变恢复的模拟问题，也解决了使用面模型模拟形变的情况下无法描述软组织体积特性的问题。

技术方案：本发明所采用的技术方案是一种软组织按压和形变恢复的模拟方法，包括以下步骤：

(1)采集软组织的医学影像数据，使用基于三角拓扑结构的质点弹簧模型构建软组织面模型，所述质点弹簧模型是将每个三角拓扑结构的三个顶点作为质点，质点之间由弹簧组件连接，所述弹簧组件由结构弹簧、弯曲弹簧和阻尼器并联组成；其中，弯曲弹簧对软组织表面内的某个质点x_i的弯曲弹簧力

表达式如下：

其中

表示连接质点x_i和x_j的弯曲弹簧的弹性系数，p⁰(j)-p⁰(i)表示质点x_i和x_j间的初始位移向量，p(j)-p(i)表示质点x_i和x_j间的当前位移向量，θ_ij表示初始位移向量和当前位移向量之间的夹角，d_ij为弯曲弹簧力的方向向量。

(2)使用粒子群算法对质点弹簧模型中的结构弹簧和弯曲弹簧的弹性系数、阻尼器的阻尼系数进行优化；所述的粒子群算法中定义适应度函数为：

其中Θ表示两个模型中所有对应质点位置差距的总和，p_i表示第i个质点在质点弹簧模型中所处的位置，q_i表示第i个质点在有限元模型中所处的位置，k^s表示结构弹簧的弹性系数，k^f表示弯曲弹簧的弹性系数，d表示弹簧的阻尼系数，E表示杨氏模量，v表示泊松比。

所述的粒子群算法中的更新公式为：

v_i＝w·v_i+c₁·rand()·(pbest_i-x_i)+c₂·rand()·(gbest-x_i)

x_i＝x_i+v_i

其中v_i，x_i分别表示第i个粒子的速度和位置，w表示惯性因子，一般取[0，1]区间的数，这里假设其值为0.5，c₁，c₂表示学习因子，通常取固定值2.0，rand()表示介于(0，1)之间的随机数，pbest_i表示第i个粒子的历史最优解，gbest表示粒子群体的全局最优解。第二个式子表示循环中对x_i的更新。

(3)建立软组织内部体模型，所述内部体模型是将软组织内部划分为一系列与软组织表面相平行的平面，每个平面上存在有数目相同的有限个质点，每个质点在其相邻的层间都有其对应的质点，且相邻平面之间通过约束弹簧相连，所述约束弹簧由阻尼器串接第二弹簧后与第一弹簧并联组成；其中，约束弹簧的本构方程为：

上式中σ₁表示应力，σ₂表示应力的导数，ε₁表示应变，ε₂表示应变的导数，E₁，E₂分别表示两个弹簧的弹性系数，η表示阻尼器的阻尼系数。

(4)将步骤(2)中获得的弹性系数和阻尼系数的全局最优解代入建立的软组织表面质点的动力学方程，计算得到软组织中表面质点的形变值，并通过软组织内部层间约束模型的应变分布函数，得到软组织内部质点域的形变分布情况。其中，软组织表面质点的动力学方程为：

其中

式中m_i，x_i分别表示第i个质点的质量和位置，

分别表示第i个质点所受的阻尼力、结构弹簧力、弯曲弹簧力和外力，n表示与第i个质点相连的质点数。

软组织内部层间约束模型的应变分布函数为：

式中y_i表示与软组织表面受外力作用的质点P_i对应的内部质点域中的质点P_i′的形变量，y_i表示内部质点域中质点P′_j的形变量，d_ij表示内部质点域中质点P_i′和P′_j间的距离。

有益效果：相对于现有技术，本发明具有以下优点：(1)使用改进的质点弹簧模型构建软组织面模型，添加了弯曲弹簧来恢复撤销外力作用后软组织表面发生的弯曲形变，提高了软组织模拟的形变精确度和形变恢复能力；(2)使用粒子群算法对质点弹簧模型参数进行优化，建立了其参数与描述软组织物理性质的杨氏模量、泊松比之间的函数关系，提高了软组织模拟的真实性；(3)使用基于层状的内部质点域约束构建软组织体模型，并在层状平面之间采用约束弹簧进行连接，以此模拟软组织内部体积力和形变过程的粘弹特性。

附图说明

图1是本发明的流程图；

图2是本发明所述的质点弹簧模型中弹簧组件的模型结构图；

图3是本发明所述的软组织内部质点约束模型示意图；

图4是本发明所述的内部约束模型中约束弹簧的模型结构图。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

本发明所述的软组织按压和形变恢复的模拟方法，流程如图1所示，包括四个部分：对软组织表面进行添加弯曲弹簧元件的质点弹簧面模型构建、使用粒子群算法优化模型参数以及基于约束弹簧的内部质点域约束对软组织内部建模。具体实施步骤如下：

(1)根据医学影像数据，使用基于三角拓扑结构的质点弹簧模型构建软组织面模型，并添加弯曲弹簧元件提高软组织形变恢复的模拟精度。使用OpenGL将CT扫描采集的软组织医学图像数据进行三维几何重现，并采用三角网格化法处理图像数据得到软组织表面的三角网格模型，然后使用质点弹簧模型对软组织表面建模。三角网格化法是一种成熟算法，将图形表面使用三角拓扑元件来构建。质点弹簧模型是将每个三角拓扑元件的三个顶点作为质点，质点之间由弹簧组件连接，在此我们添加了弯曲弹簧元件，即每个弹簧组件分别由结构弹簧、弯曲弹簧和阻尼器并联组成，如图2所示。将外力f_ext作用于质点弹簧模型的某一质点x_i，根据胡克定律，结构弹簧力

如下式计算：

其中，l_ij＝p(j)-p(i)

上式中

表示连接质点x_i和x_j的结构弹簧的弹性系数，l_ij ⁰表示质点x_i和x_j的初始距离，l_ij表示当前质点x_i和x_j间的位移向量，由质点i和质点j的位置p(i)、p(j)获得，||l_ij||为向量l_ij的模。

同样地，作用于质点x_i的阻尼力

如下式计算：

其中，l_ij＝p(j)-p(i)

上式中

表示连接质点x_i和x_j的阻尼器的阻尼系数，表示

表示应变的导数。

当撤销作用于软组织表面的外力f_ext后，质点x_i回到其初始位置的恢复速度缓慢甚至在大形变下无法精准地回到其初始位置，因此引入弯曲弹簧力

确保软组织恢复至初始形状，如下式计算：

上式中

表示连接质点x_i和x_j的弯曲弹簧的弹性系数，p⁰(j)-p⁰(i)表示质点x_i和x_j间的初始位移向量，p(j)-p(i)表示质点x_i和x_j间的当前位移向量，θ_ij表示初始位移向量和当前位移向量之间的夹角，d_ij为弯曲弹簧力的方向向量。

(2)使用粒子群算法对质点弹簧模型的弹性系数和阻尼系数进行优化。我们假设将基于质点弹簧模型和基于有限元模型的软组织表面离散为相同的N个质点，并认为所有结构弹簧都具有相同大小的弹性系数，所有弯曲弹簧都具有相同大小的弹性系数，所有阻尼器都具有相同大小的阻尼系数，然后在模型表面作用相同的外力f_ext，设置优化问题的目标为最小化质点弹簧模型和有限元模型中所有同一质点位置间的差距，采用该目标是为了保证模拟的精确性，因为有限元模型的模拟结果是目前模型中精度最高的。在计算过程中，有限元模型的模拟结果作为先验数据被直接调用，不需要重新模拟计算。在此采用了粒子群算法，定义适应度函数：

上式中Θ表示两个模型中所有对应质点位置差距的总和，p_i表示第i个质点在质点弹簧模型中所处的位置，q_i表示第i个质点在有限元模型中所处的位置，k^s表示结构弹簧的弹性系数，k^f表示弯曲弹簧的弹性系数，d表示弹簧的阻尼系数，E表示杨氏模量，v表示泊松比。

然后根据粒子群算法的更新公式更新质点的速度和位置。更新公式如下所示：

v_i＝w·v_i+c₁·rand()·(pbest_i-x_i)+c₂·rand()·(gbest-x_i)

x_i＝x_i+v_i

上式中v_i，x_i分别表示第i个粒子的速度和位置，w表示惯性因子，一般取[0，1]区间的数，这里假设其值为0.5，c₁，c₂表示学习因子，通常取固定值2.0，rand()表示介于(0，1)之间的随机数，pbest_i表示第i个粒子的历史最优解，gbest表示粒子群体的全局最优解。第二个式子表示循环中对x_i的更新。

最后通过迭代更新公式最小化适应度函数来获取全局最优解，得到的全局最优解即为弹性系数和阻尼系数的最优解，使得质点弹簧模型的形变仿真效果最大程度的接近于有限元模型的形变仿真效果。

(3)使用基于层状的内部质点域约束对软组织内部构建体模型，在层状平面之间采用约束弹簧进行连接，以此模拟内部体积力和形变过程的粘弹性特性。我们将软组织内部划分为一系列与软组织表面相平行的平面，平面上的质点域是由表面S上的质点投影到平面T上所得质点组成的集合，且每个平面之间通过约束弹簧相连，如图3所示。

图3中，当模型表面S的某一质点P_i受到外力f_ext的作用时，产生的应力会通过约束弹簧传递给P_i在平面T上相应的质点P_i′，并一层一层的往下传递来模拟形变过程的粘弹性特性，最后当作用于P_i的外力撤销时，内部质点将通过约束弹簧慢慢恢复至其初始位置。

约束弹簧模型如图4所示，由阻尼器串接第二弹簧后与第一弹簧并联组成，其应力和应变的本构方程为：

上式中σ₁表示应力，σ₂表示应力的导数，ε₁表示应变，ε₂表示应变的导数，E₁，E₂分别表示两个弹簧的弹性系数，η表示阻尼器的阻尼系数。

由约束弹簧模型的本构方程推导出约束模型的应力松弛方程：

上式中c_b，c₁表示与软组织本身有关的材料系数，τ₁表示时间常量，ε表示应变，t表示时间。

(4)当软组织表面质点P_i受到外力作用后，通过步骤(2)获得弹性系数和阻尼系数的全局最优解，建立软组织表面质点的动力学方程，并把得到的最优解代入动力学方程中计算获得软组织模型中表面质点的形变值。同时，通过约束弹簧向内部质点域中其对应的点P_i′传递应力，从而使点P_i′产生应变，并通过软组织内部层间约束模型的应变分布函数，将点P_i′的应变传递给质点域中的其他质点。当撤销作用于表面的外力时，内部质点域上的质点不再计算应变分布，软组织表面的质点在弯曲弹簧的约束下恢复到其初始位置，内部质点域的质点在约束弹簧的作用下恢复到其初始位置，两个过程满足外力为0的软组织表面质点的动力学方程。

软组织表面质点的动力学方程为：

其中

上式中m_i，x_i分别表示第i个质点的质量和位置，

分别表示第i个质点所受的阻尼力、结构弹簧力、弯曲弹簧力和外力，n表示与第i个质点相连的质点数。

软组织内部层间约束模型的应变分布函数为：

上式中y_i表示与软组织表面受外力作用的质点P_i对应的内部质点域中的质点P_i′的形变量，y_j表示内部质点域中质点P′_j的形变量，d_ij表示内部质点域中质点P_i′和P′_j间的距离。

本发明算法可以保证计算简单的同时提高模拟肝脏形变的精确度、计算效率以及形变恢复过程的模拟精度，让操作者在人机交互的过程中可以感受到软组织逼真的体积特性。

Claims (5)

1.一种软组织按压和形变恢复的模拟方法，其特征在于，包括以下步骤：

(1)采集软组织的医学影像数据，使用基于三角拓扑结构的质点弹簧模型构建软组织面模型，所述质点弹簧模型是将每个三角拓扑结构的三个顶点作为质点，质点之间由弹簧组件连接，所述弹簧组件由结构弹簧、弯曲弹簧和阻尼器并联组成；

(2)使用粒子群算法对质点弹簧模型中的结构弹簧和弯曲弹簧的弹性系数、阻尼器的阻尼系数进行优化；

(3)建立软组织内部体模型，所述内部体模型是将软组织内部划分为一系列与软组织表面相平行的平面，每个平面上存在有数目相同的有限个质点，每个质点在其相邻的层间都有其对应的质点，且相邻平面之间通过约束弹簧相连，所述约束弹簧由阻尼器串接第二弹簧后与第一弹簧并联组成；

(4)将步骤(2)中获得的弹性系数和阻尼系数的全局最优解代入建立的软组织表面质点的动力学方程，计算得到软组织中表面质点的形变值，并通过软组织内部层间约束模型的应变分布函数，得到软组织内部质点域的形变分布情况。

2.根据权利要求1所述的软组织按压和形变恢复的模拟方法，其特征在于：步骤(1)中所述的弯曲弹簧对软组织表面内的某个质点x_i的弯曲弹簧力

表达式如下：

其中

表示连接质点x_i和x_j的弯曲弹簧的弹性系数，p⁰(j)-p⁰(i)表示质点x_i和x_j间的初始位移向量，p(j)-p(i)表示质点x_i和x_j间的当前位移向量，θ_ij表示初始位移向量和当前位移向量之间的夹角，d_ij为弯曲弹簧力的方向向量。

3.根据权利要求1所述的软组织按压和形变恢复的模拟方法，其特征在于：步骤(2)中所述的粒子群算法中定义适应度函数为：

其中Θ表示两个模型中所有对应质点位置差距的总和，p_i表示第i个质点在质点弹簧模型中所处的位置，q_i表示第i个质点在有限元模型中所处的位置，k^s表示结构弹簧的弹性系数，k^f表示弯曲弹簧的弹性系数，d表示弹簧的阻尼系数，E表示杨氏模量，v表示泊松比；

所述的粒子群算法中的更新公式为：

v_i＝w·v_i+c₁·rand()·(pbest_i-x_i)+c₂·rand()·(gbesf-x_i)

x_i＝x_i+v_i

其中v_i，x_i分别表示第i个粒子的速度和位置，w表示惯性因子，一般取[0，1]区间的数，这里假设其值为0.5，c₁，c₂表示学习因子，通常取固定值2.0，rand()表示介于(0，1)之间的随机数，pbest_i表示第i个粒子的历史最优解，gbest表示粒子群体的全局最优解。

4.根据权利要求1所述的软组织按压和形变恢复的模拟方法，其特征在于：步骤(3)中所述的约束弹簧的本构方程为：

式中σ₁表示应力，σ₂表示应力的导数，ε₁表示应变，ε₂表示应变的导数，E₁，E₂分别表示第一弹簧和第二弹簧的弹性系数，η表示阻尼器的阻尼系数。

5.根据权利要求1所述的软组织按压和形变恢复的模拟方法，其特征在于：步骤(4)中所述的软组织表面质点的动力学方程为：

其中

式中m_i，x_i分别表示第i个质点的质量和位置，f_i ^d，f_i ^s，f_i ^f，

分别表示第i个质点所受的阻尼力、结构弹簧力、弯曲弹簧力和外力，n表示与第i个质点相连的质点数；

软组织内部层间约束模型的应变分布函数为：

式中y_i表示与软组织表面受外力作用的质点P_i对应的内部质点域中的质点P_i′的形变量，y_j表示内部质点域中质点P′_j的形变量，d_ij表示内部质点域中质点P_i′和P′_j间的距离。

Images