CN110289103A - 软组织局部压缩形变模拟方法 - Google Patents

Abstract

本发明公开了一种软组织局部压缩形变模拟方法，包括：(1)采集软组织的医学图像数据，通过四面体网格生成方法将软组织的医学图像数据处理为三维四面体网格模型；(2)构造四面体质点‑弹簧模型并建立与质点‑弹簧模型对应的运动学方程；(3)利用模拟退火算法优化质点‑弹簧模型中的弹簧刚度和阻尼系数；(4)将软组织压缩变形划分为局部变形区和非变形区；(5)利用Verlet显示积分法对局部变形区的点进行运动学方程求解，获得在压力作用下软组织的压缩形变；(6)设定过度压缩约束条件，对不满足过度压缩约束条件的质点进行校正。本发明在保证计算效率的同时提高了模拟的精度，更快速、真实的反映了软组织的形变情况。

Description

软组织局部压缩形变模拟方法

技术领域

本发明涉及一种软组织形变模拟方法，尤其涉及一种软组织局部压缩形变模拟方法。

背景技术

近年来虚拟现实技术蓬勃发展，虚拟手术模拟系统就是其重要应用之一。在医学手术过程中要经常对软组织进行按压操作，医生或医学生可通过虚拟手术模拟系统进行按压手术训练，积累经验提高手术技能，同时降低培训成本。在虚拟手术模拟中，建立真实的软组织形变模型至关重要。

使软组织形变模型满足真实性和实时性是最具挑战性问题，现今主要使用的模型包括有限元模型和质点-弹簧模型。有限元模型的优点是参数连续化，精度高，适应性强，可模拟大的软组织形变，但其建模复杂，需进行大量复杂的数值计算，实时性较差；质点-弹簧模型参数不需要连续化，计算快、效率高，但其精度较差。因此建立一个具有良好变形效果且计算效率高的模型对虚拟手术模拟系统来说至关重要。

发明内容

发明目的：针对以上问题，本发明提出一种软组织局部压缩形变模拟方法，在保证计算效率的同时提高模拟精度。

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

(1)采集软组织的医学图像数据，通过四面体网格生成方法将软组织的医学图像数据处理为三维四面体网格模型；

(2)构造四面体质点-弹簧模型并建立与质点-弹簧模型对应的运动学方程；其中，四面体质点-弹簧模型是将三维四面体网格中四面体的四个顶点作为质点，质点之间由弹簧连接。与该质点-弹簧模型对应的一种运动学方程表达式为：

其中m_i表示质点i的质量，x_i表示质点的位移，t表示时间，F_i ^int表示质点i所受内力，F_i ^vol表示质点i所受体积力，F_i ^ext表示质点i所受外力；

在该表达式中，质点i所受内力F_i ^int被设置为与弹簧刚度系数和阻尼系数相关的量，具体的表达式为：

其中，任意质点i及其任意一个邻点j的位置分别为p_i和p_j，连接两质点间的弹簧刚度系数和阻尼系数分别为k_ij和d_ij，表示弹簧初始长度，S_i表示为所有连接到质点i的弹簧集合，v_j和v_i分别表示质点j和i的速度；

为保证体积守恒，在每个点上定义了体积力F_i ^vol为

其中，p_wb表示四面体w的重心，Q表示所有包含质点i的四面体集合，和分别表示四面体w当前的体积和初始体积。

(3)利用模拟退火算法优化质点-弹簧模型中的弹簧刚度和阻尼系数；一种可执行的模拟退火算法优化过程包括以下步骤：

(3.1)给定初始温度、终止温度和退火系数，随机初始化一组参数S＝(k_ijd_ij)进行循环计算；

(3.2)在温度T_n下，在当前参数S的领域内随机产生新的一组参数S′，计算目标函数值增量Δf＝f(S′)-f(S)，若Δf＜0则令S＝S′，反之根据Metropolis准测，以min{1，exp(-Δf/T_n)}＞random[0，1]的概率接受S′作为当前较优的一组参数S，其中目标函数定义为用于图像匹配的Hausdorff距离函数：

f(S)＝f(k_ij，d_ij)＝H(N^MSM，N^FEM)＝max{h(N^MSM，N^FEM)，h(N^FEM，N^MSM)}

其中N^MSM、N^FEM、p_m、p_f、d(p_m，p_f)分别表示质点-弹簧模型关键点集合、有限元模型关键点集合、质点-弹簧模型关键点位置、有限元模型关键点位置、点p_m和p_f之间的绝对距离。

(3.3)判断在温度T_n下是否达到规定的迭代次数，若满足条件则转到步骤(3.4)，否则返回步骤(3.2)；

(3.4)进行降温操作，令n＝n+1，T_n+1＝λT_n，其中n为迭代次数，λ为退火系数；

(3.5)判断是否达到设定的终止温度，若满足终止条件则输出一组最优参数S＝(k_ij，d_ij)，否则返回步骤(3.2)。

(4)将软组织压缩变形划分为局部变形区和非变形区；具体的，对两个区域的划分包括以下过程：

(4.1)当一个压头在均匀的压力下作用于软组织，垂直应力F_⊥与接触压力的关系为：

其中F_⊥表示垂直应力，F_c表示施加在软组织上的接触压力，h表示负载中心点在接触表面以下的垂直深度，r表示压头半径。

其中，接触压力F_c表示为：

其中表示为单向外部负载，E表示杨氏模量，μ表示泊松比，x_z表示在单向外部负载作用下软组织变形位移。

(4.2)根据步骤(4.1)中的公式计算出最小垂直应力分量F_⊥min对应的负载中心点在接触表面以下的垂直深度h_max，h_max作为局部变形范围的半径，半径以内区域为局部变形区，半径以外区域为非变形区。

(5)利用Verlet显示积分法对局部变形区的点进行运动学方程求解，获得在压力作用下软组织的压缩形变；

(6)设定过度压缩约束条件，对不满足过度压缩约束条件的质点进行校正。校正过程包括：

(6.1)定义过度压缩约束：

其中，τ_c表示临界压缩比，p_i和p_j分别为任意质点i及其任意一个邻点j的位置，表示质点间的初始长度；

(6.2)若不满足约束，对点进行校正p_i′＝p_i+Δp_i，p_j′＝p_j+Δp_j，

其中，

有益效果：相对于现有技术，本发明具有以下优点：(1)利用模拟退火算法确定质点-弹簧模型的最佳弹簧刚度和阻尼系数，提高了模拟的精度；(2)设定过度压缩约束，防止局部压力过大导致模型出现“超弹性”效应，从而更真实模拟软组织变形；(3)根据最小垂直应力计算对应负载中心点在接触表面以下的垂直深度，以此为半径将软组织划分为局部变形区域和非变形区域，只对局部变形区域点利用Verlet显示积分法进行计算，降低计算成本，计算效率高。

附图说明

图1是本发明的流程图；

图2是本发明所述的四面体网格模型的拓扑结构图；

图3是本发明所述的局部变形区域图；

图4是本发明所述的压头和软组织在单向外部负载作用下相互作用图。

具体实施方式

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

本发明所述的一种软组织局部压缩形变模拟方法应用在虚拟手术仿真系统中，利用虚拟现实技术进行医学手术的模拟训练。虚拟手术是借由医学图像数据出发，应用计算机图形学重构出虚拟人体软组织模型，模拟出虚拟的医学环境，并利用触觉交互设备与训练者进行交互的手术系统。虚拟手术仿真系统为医生提供了一个虚拟的3D环境以及可交互操作平台，逼真的模拟临床手术的全过程。本发明所述的软组织局部压缩形变模拟方法，属于该系统中虚拟医学环境这一部分，其具体的流程如图1所示，包括以下步骤：

(1)对软组织进行医学图像数据采集，例如电子计算机断层扫描(CT)，然后利用Amira软件对软组织的CT数据进行处理导出三维四面体网格模型。如图2所示，为四面体网格模型的拓扑结构图。

(2)将三维四面体网格模型中四面体的四个顶点作为质点，质点之间认为由弹簧连接，由此构造出四面体质点-弹簧模型。建立与质点-弹簧模型对应的运动学方程，具体如下：

设质点-弹簧模型中的任意质点i及其任意一个邻点j的位置分别为p_i和p_j，连接两质点间的弹簧刚度系数和阻尼系数分别为k_ij和d_ij，弹簧初始长度用表示，则质点i所受内力F_i ^int可以表示为

其中S_i表示为所有连接到质点i的弹簧集合，v_j和v_i分别表示质点j和i的速度。

为保证体积守恒，在每个点上定义了体积力，设四面体w的重心为p_wb，则质点i所受体积力F_i ^vol为

其中Q_i表示所有包含质点i的四面体集合，和分别表示四面体w当前的体积和初始体积。

计算内力和体积力后，带入质点i的运动学方程，表示如下

其中m_i表示质点i的质量，F_i ^ext表示质点i所受外力。

(3)利用模拟退火算法优化质点-弹簧模型参数，该步骤是为提高软组织模型模拟的真实性(精度)，具体如下：

(3.1)给定初始温度T_n＝100(n＝0)，终止温度T_f＝0.001，退火系数λ＝0.98，内循环次数L＝200，随机初始化一组参数S＝(k_ijd_ij)；

(3.2)在温度T_n下，在当前的一组参数S的领域内随机产生新的一组参数S′，计算目标函数值增量Δf＝f(S′)-f(S)，若Δf＜0则令S＝S′，反之根据Metropolis准测，以min{1，exp(-Δf/T_n)}＞random[0，1]的概率接受S′作为当前较优的一组参数S，即S＝S′，其中目标函数定义为用于图像匹配的Hausdorff距离函数：

f(S)＝f(k_ij，d_ij)＝H(N^MSM，N^FEM)＝max{h(N^MSM，N^FEM)，h(N^FEM，N^MSM)} (4)

其中N^MSM、N^FEM、p_m、p_f、d(p_m，p_f)分别表示质点-弹簧模型关键点集合、有限元模型关键点集合、质点-弹簧模型关键点位置、有限元模型关键点位置、点p_m和p_f之间的绝对距离。

(3.3)判断在温度T_n下是否达到规定的迭代次数，若满足条件则转到步骤(3.4)，否则返回步骤(3.2)；

(3.4)进行降温操作，令n＝n+1，T_n+1＝λT_n；

(3.5)判断是否达到设定的终止温度，若满足终止条件则输出一组最优参数S＝(k_ijd_ij)，否则返回步骤(3.2)。

(4)如图3所示，软组织变形一般在局部发生，局部外变形影响很小。为降低计算成本，将软组织压缩变形划分为局部变形区和非变形区，同时根据研究表明垂直应力分量在材料失效中起主要作用，垂直应力F是软组织模型在受到单向外部负载时内部各部分相互作用的内力。考虑垂直应力分量在软组织压缩中的作用，具体如下：

如图4所示是一个压头和软组织在单向外部负载作用下的相互作用图。假设压头在均匀的压力下作用于软组织，则垂直应力F_⊥为：

其中F_⊥、F_c、h、r分别表示垂直应力、施加在软组织上的接触压力、负载中心点在接触表面以下的垂直深度、压头半径。

接触压力F_c可以表示为：

其中E、μ、x_z分别表示为单向外部负载、杨氏模量、泊松比、在单向外部负载作用下软组织变形位移。单向外部负载即压头对模型施加的外力。

将最小垂直应力分量F_⊥min对应的负载中心点在接触表面以下的垂直深度h_max作为局部变形范围的半径，半径以内区域为局部变形区，半径以外区域为非变形区，F_⊥min可通过实验计算确定，然后利用公式(7)和(8)计算出局部变形半径h_max，在每次迭代中根据接触点外部负载和位移动态确定局部变形范围，非变形区的点将不进行变形计算。

(5)对局部变形区的点进行变形计算，由于Verlet显示积分法计算快、精度高，故用其求解运动学方程，具体如下：

其中x_i(t+Δt)、x_i(t)、v_i(t)、Δt、F_i(t)、v_i(t+Δt)、v_i(()、F_i(t+Δt)分别表示质点i在t+Δt时刻、t时刻的位移、质点i在t时刻的速度、一个时间步长、质点i在t时刻所受合力、质点i在在t+Δt时刻、t时刻的速度、质点i在t+Δt时刻所受合力。弹簧刚度和阻尼系数与质点所受内力的关系见公式(1)，其影响了质点所受内力，进而影响了运动学方程中的合力F_i(t)。

(6)当一个较大的力集中作用于软组织局部时，质点-弹簧模型会陷入局部过度压缩的问题，即“超弹性”效应，故设定了过度压缩约束。当弹簧长度小于时，利用过度压缩约束将弹簧移回到过度压缩约束定义如下：

其中rc表示临界压缩比，若不满足不等式(11)，就对点进行校正p_i′＝p_i+Δp_i，p_j′＝p_j+Δp_j，校正相的计算式如下：

至此，我们通过本发明所述的软组织局部压缩形变模拟方法，获得了更精确的形变模拟结果，并且在提高精度的前提下没有对计算效率做出牺牲，保证了在虚拟手术中的实时性。

Claims (6)

1.一种软组织局部压缩形变模拟方法，其特征在于，包括以下步骤：

(1)采集软组织的医学图像数据，通过四面体网格生成方法将软组织的医学图像数据处理为三维四面体网格模型；

(2)构造四面体质点-弹簧模型并建立与质点-弹簧模型对应的运动学方程；

(3)利用模拟退火算法优化质点-弹簧模型中的弹簧刚度和阻尼系数；

(4)将软组织压缩变形划分为局部变形区和非变形区；

(5)利用Verlet显示积分法对局部变形区的点进行运动学方程求解，获得在压力作用下软组织的压缩形变；

(6)设定过度压缩约束条件，对不满足过度压缩约束条件的质点进行校正。

2.根据权利要求1所述的软组织局部压缩形变模拟方法，其特征在于，步骤(2)中所述的四面体质点-弹簧模型，是将三维四面体网格中四面体的四个顶点作为质点，质点之间由弹簧连接。

3.根据权利要求2所述的软组织局部压缩形变模拟方法，其特征在于，步骤(2)中所述的与质点-弹簧模型对应的运动学方程为：

其中m_i表示质点i的质量，x_i表示质点的位移，t表示时间，F_i ^int表示质点i所受内力，F_i ^vol表示质点i所受体积力，F_i ^ext表示质点i所受外力；

质点i所受内力F_i ^int为：

其中，任意质点i及其任意一个邻点j的位置分别为p_i和p_j，连接两质点间的弹簧刚度系数和阻尼系数分别为k_ij和d_ij，表示弹簧初始长度，S_i表示为所有连接到质点i的弹簧集合，v_j和v_i分别表示质点j和i的速度；

质点i所受体积力F_i ^vol为

其中，p_wb表示四面体w的重心，Ω_i表示所有包含质点i的四面体集合，和分别表示四面体w当前的体积和初始体积。

4.根据权利要求1所述的软组织局部压缩形变模拟方法，其特征在于，所述步骤(3)包括以下过程：

(3.1)给定初始温度、终止温度和退火系数，随机初始化一组参数S＝(k_ij，d_ij)进行循环计算，其中k_ij和d_ij分别为弹簧刚度系数和阻尼系数；

(3.2)在温度T_n下，在当前参数S的领域内随机产生新的一组参数S′，计算目标函数值增量Δf＝f(S′)-f(S)，若Δf＜0则令S＝S′，反之根据Metropolis准测，以min{1，exp(-Δf/T_n)}＞random[0，1]的概率接受S′作为当前较优的一组参数S，其中目标函数定义为用于图像匹配的Hausdorff距离函数：

f(S)＝f(k_ij，d_ij)＝H(N^MSM，N^FEM)＝max{h(N^MSM，N^FEM)，h(N^FEM，N^MSM)}

其中N^MSM、N^FEM、p_m、p_f、d(p_m，p_f)分别表示质点-弹簧模型关键点集合、有限元模型关键点集合、质点-弹簧模型关键点位置、有限元模型关键点位置、点p_m和p_f之间的绝对距离。

(3.3)判断在温度T_n下是否达到规定的迭代次数，若满足条件则转到步骤(3.4)，否则返回步骤(3.2)；

(3.4)进行降温操作，对温度迭代，T_n+1＝λT_n，其中n为迭代次数，λ为退火系数；

(3.5)判断是否达到设定的终止温度，若满足终止条件则输出一组最优参数S＝(k_ij，d_ij)，否则返回步骤(3.2)。

5.根据权利要求1所述的软组织局部压缩形变模拟方法，其特征在于，所述步骤(4)包括以下过程：

(4.1)当一个压头在均匀的压力下作用于软组织，垂直应力F_⊥与施加在软组织上的接触压力F_c的关系为：

其中，h表示负载中心点在接触表面以下的垂直深度，r表示压头半径，接触压力F_c表达式为：

其中表示为单向外部负载，E表示杨氏模量，μ表示泊松比，x_z表示在单向外部负载作用下软组织变形位移。

(4.2)根据步骤(4.1)中的公式计算出最小垂直应力分量F_⊥min对应的负载中心点在接触表面以下的垂直深度h_max，h_max作为局部变形范围的半径，半径以内区域为局部变形区，半径以外区域为非变形区。

6.根据权利要求1所述的软组织局部压缩形变模拟方法，其特征在于，所述步骤(6)中包括以下过程：

(6.1)定义过度压缩约束：

其中，τ_c表示临界压缩比，p_i和p_j分别为任意质点i及其任意一个邻点j的位置，表示弹簧的初始长度；

(6.2)若不满足约束，对点进行校正p_i′＝p_i+Δp_i，p_j′＝p_j+Δp_j，

其中，