CN103400023B - 软组织形变仿真方法 - Google Patents

Abstract

本发明涉及一种软组织形变仿真方法，包括以下步骤：建立软组织的生物力学模型、并对其中的各个质点进行初始化；力反馈设备对软组织施加作用力，进行碰撞检测；用改进的欧拉算法计算运动状态信息；将模型每个时间步长的状态输出到显示屏上，动态显示软组织形变过程；计算反馈力并输出触觉反馈。通过采用上述步骤，能够有效解决虚拟手术仿真中的实时性、精确性以及反馈力光滑性问题，提高了软组织形变仿真的精度和实时性，从而满足虚拟手术仿真的需要。

Description

软组织形变仿真方法

技术领域

本发明涉及一种软组织形变仿真方法。

背景技术

虚拟手术仿真是虚拟现实技术的重要应用之一，其目标是利用各种医学影像数据以及虚拟现实技术在计算机中建立逼真的虚拟手术环境，使得医生或者受训者可以借助该虚拟环境进行手术训练。

在虚拟手术仿真中，软组织器官的形变模型决定了虚拟手术仿真的视觉效果和力反馈精度。常见的软组织形变计算模型分为两大类：基于几何的形变模型和基于物理的形变模型。基于几何的形变模型仅仅考虑了几何形态的变化，而忽略了软组织的实际力学本构方程以及形变过程中物体质量、力或其他物理现象的作用，因此不能真实的反映软组织的形变过程，该模型目前已较少使用；基于物理的形变模型则基于软组织的力学本构方程，通过相应的计算模型得出组织受力时的形变，能够更加真实的反映组织的形变，因此目前使用较多。目前，基于物理的形变模型主要有有限元模型和质点-弹簧模型。有限元模型的优点是参数连续化，在大形变情况下也能模拟组织形变，具有较高的精度和适应性，但其建模复杂，计算量大，需要大量的预处理工作，尤其在切割过程中难以实现快速的网格重构，因此，实时性较差。与有限元模型相比，质点-弹簧模型不需要将参数连续化，容易实现，运算速度较快。而质点-弹簧模型一般分为面模型和体模型两种，面模型的网格密度可以根据仿真的要求进行划分，计算效率高，可以满足仿真的实时性，但是，只能反映组织表面的形变特征，不能刻画三维实体的内部结构，是以牺牲精度为代价来保证实时性，体模型是对具有一定体积的实体组织进行建模，能够同时表现表面和内部结构，具有较高的精度，但实时性较差。

对于上述模型，目前求解其二阶微分方程的主要方法有欧拉法和龙格-库塔法。显式欧拉法计算量小，但是一般只一阶收敛，精度不高，为实现逼真的模拟效果，时间步长需要设置很小，从而导致整个形变过程延长，并且由于存在进退性冲击波，不适合用于质点-弹簧模型。隐式欧拉法精度高，稳定性好，但是计算量却较大。龙格-库塔法也具有较高的精度和稳定性，但是，在一步计算中需要计算四次函数的值，计算量大，实时性较差。在虚拟手术仿真中，光滑的触觉反馈对实时性的要求更高。因此，如何设计出一种支持光滑触觉反馈的快速且逼真的软组织形变仿真方法，尽可能的满足虚拟手术对实时性以及精确性的要求，已成为虚拟手术面临的首要问题。

发明内容

为克服以上现有技术的不足，本发明要解决的技术问题是提供一种软组织形变仿真方法，能够有效构建软组织生物力学模型，并对其中质点的速度、位移以及受力信息进行高效求解，从而有效解决虚拟手术仿真中的实时性、精确性以及反馈力光滑性问题，进而满足虚拟手术仿真的需要。

本发明的技术方案是：

一种软组织形变仿真方法，包含以下步骤：

步骤1）：采集软组织的数据信息，并采用基于四面体的质点-弹簧体模型建立软组织的生物力学模型，该生物力学模型由n个质点组成，并且对于其中任意一个质点i，满足以下方程：

$m_i\frac{{\∂}^2{x}_i}{{\∂t}^2}+F_{\mathrm{in}}^i=F_{\mathrm{ext}}^i$

$F_{\mathrm{in}}^i=\underset{j=j1}{\overset{\mathrm{jq}}{\mathrm{\Σ}}}\mathrm{\μ}\left(\right|x_j-x_i|-{|x_j-x_i|}^0)\frac{(x_j-x_i)}{|x_j-x_i|}+\underset{j=j1}{\overset{\mathrm{jq}}{\mathrm{\Σ}}}\mathrm{\η}(v_j-v_i)$

$a_i=\frac{F_{\mathrm{ext}}^i-F_{\mathrm{in}}^i}{m_i}$

其中，m_i为质点i的质量，x_i表示质点i的位置矢量，表示与质点i相连的所有质点j（包括j1至jq）施加于质点i的内力，x_j表示质点j的位置矢量，表示质点i所受外力，|x_j-x_i|表示两个质点间位置矢量之差，|x_j-x_i|⁰是发生形变前弹簧的长度，v_i为质点i的速度，v_j为质点j的速度，a_i为质点i的加速度，μ为弹簧的弹性系数，η为阻尼器的阻尼系数；

步骤2）：依据步骤1）所建的生物力学模型，对其中的各个质点进行初始化，构建模型的初始状态，并且计算该模型中的每个质点与相连质点间的弹簧初始长度；

步骤3）：外接的力反馈设备对软组织施加作用力，进行碰撞检测，确定软组织上发生碰撞的质点和受力发生形变的区域，以及软组织被按压或者拉伸的长度；

步骤4）：用改进的欧拉算法计算发生形变区域内的各个质点的运动状态信息；

步骤5）：依据步骤3）和步骤4）计算反馈力，并将该反馈力输出至力反馈设备；

步骤6）：循环执行步骤4）至步骤5），计算每个时间步长中各质点的运动状态，并在显示器上动态显示软组织的形变过程。

上述软组织形变仿真方法，其中步骤2）中的初始化包括初始化各个质点的位置、质量、速度、加速度和受力信息。

上述软组织形变仿真方法，其中步骤4）中的运动状态信息包括发生形变区域内的各个质点的位置、速度和受力随时间变化的信息。

上述软组织形变仿真方法，其中步骤4）的计算过程具体为：

步骤4.1）：经力反馈设备对软组织施加外作用力、并经碰撞检测，返回发生碰撞的质点的序列号r及所受外力并且当时刻k=0时，

步骤4.2）：计算质点r的加速度：

$a_r\left(k\right)=\frac{F_{\mathrm{ext}}^r-F_{\mathrm{in}}^r\left(k\right)}{m_r}$

步骤4.3）：用改进的欧拉算法求解质点r的速度和位移，其中用显式欧拉法对速度v_r进行迭代求解，用隐式欧拉法对位置矢量x_r进行求解：

$v_r^{k+1}=v_r^k+\mathrm{\Δt}\·a_r\left(k\right)$

$x_r^{k+1}=x_r^k+\frac{\mathrm{\Δt}}{2}\·(v_r^k+v_r^{k+1})$

其中，为质点r在时刻k的速度向量，为质点r在时刻k+1的速度向量，为质点r在时刻k的位置向量，为质点r在时刻k+1的位置向量；

步骤4.4）：计算质点r开始运动后与其相连质点间弹簧的长度及弹簧的形变量，并计算质点r所受内力

步骤4.5）：循环执行步骤4.1）至4.4），计算受力区域内其他质点的位置、速度及受力信息。

本发明的有益效果是：本发明通过有效构建软组织生物力学模型，利用改进的欧拉算法求解软组织形变过程，对其中质点的速度、位移以及受力信息进行高效求解，并实现了光滑反馈力的输出，从而有效解决虚拟手术仿真中的实时性、精确性以及反馈力光滑性问题，提高了软组织形变仿真的精度和实时性，进而满足了虚拟手术仿真的需要。

附图说明

下面结合附图对本发明的具体实施方式作进一步详细的说明。

图1是软组织形变仿真方法的总体流程图示意图；

图2是软组织形变仿真方法中的质点-弹簧体模型的拓扑结构示意图。

具体实施方式

如图1和图2所示，一种软组织形变仿真方法，包含以下步骤：

步骤1）：采集软组织的数据信息，并采用基于四面体的质点-弹簧体模型建立软组织的生物力学模型，该生物力学模型由n个质点组成，并且对于其中任意一个质点i，满足以下该动力学方程：

$m_i\frac{{\∂}^2{x}_i}{{\∂t}^2}+F_{\mathrm{in}}^i=F_{\mathrm{ext}}^i$

本发明采用Kelvin黏弹性模型来反映软组织的力学特性，即质点与质点间由弹簧和阻尼器并联连接。依据Kelvin模型，当力F作用到弹簧上时，弹簧服从胡克定律，伸长量u与F成正比，即F=μu；当力F作用到阻尼器上时，阻尼器服从牛顿黏性定律，形变速度u′与力F成正比，即F=μu′。因此，当软组织受到外力作用发生形变时，两个相连质点之间的内力由弹簧产生的弹性力和阻尼器产生的粘滞力组成，即f_in=μu+ηu′。因此，施加于某质点的内力的计算方法为：

$F_{\mathrm{in}}^i=\underset{j=j1}{\overset{\mathrm{jq}}{\mathrm{\Σ}}}{f}_{\mathrm{in}}^i=\underset{j=j1}{\overset{\mathrm{jq}}{\mathrm{\Σ}}}\mathrm{\μ}\left(\right|x_j-x_i|-{|x_j-x_i|}^0)\frac{(x_j-x_i)}{|x_j-x_i|}+\underset{j=j1}{\overset{\mathrm{jq}}{\mathrm{\Σ}}}\mathrm{\η}(v_j-v_i)$

并且，质点的加速度计算方法为：

$a_i=\frac{F_{\mathrm{ext}}^i-F_{\mathrm{in}}^i}{m_i}$

上述公式中，m_i为质点i的质量，x_i表示质点i的位置矢量，表示与质点i相连的所有质点j（包括j1至jq）施加于质点i的内力，x_j表示质点j的位置矢量，表示质点i所受外力，|x_j-x_i|表示两个质点间位置矢量之差，|x_j-x_i|⁰是发生形变前弹簧的长度，v_i为质点i的速度，v_j为质点j的速度，a_i为质点i的加速度，μ为弹簧的弹性系数，η为阻尼器的阻尼系数。

步骤2）：依据步骤1）所建的生物力学模型，对其中的各个质点进行初始化，可包括初始化各个质点的位置、质量、速度、加速度和受力信息，并构建模型的初始状态，计算该模型中的每个质点与相连质点间的弹簧初始长度。上述初始化信息均存入邻接链表中。

步骤3）：外接的力反馈设备对软组织施加作用力，进行碰撞检测，确定软组织上发生碰撞的质点和受力发生形变的区域，以及软组织被按压或者拉伸的长度。

步骤4）：用改进的欧拉算法计算发生形变区域内的各个质点的运动状态信息，其中运动状态信息可包括发生形变区域内的各个质点的位置、速度和受力随时间变化的信息。

步骤5）：依据步骤3），当外接的力反馈设备对软组织施加作用力时，计算设备按压或者拉伸软组织的长度，并结合依据步骤4）得出的各个质点的运动状态信息，计算反馈力，并将该反馈力输出至力反馈设备。此时，用户可以通过操纵杆感受到力的反馈；

步骤6）：循环执行步骤4）至步骤5），计算每个时间步长中各质点的运动状态，利用OpenGL技术，将每个时间步长内计算出的质点按照位移向量输出至显示器，即可得到软组织模型在外力作用下的形变过程。

对于上述软组织形变仿真方法，其中步骤4）的计算过程可具体包括为：

步骤4.1）：经力反馈设备对软组织施加外作用力、并经碰撞检测，返回发生碰撞的质点的序列号r及所受外力并且当时刻k=0时，

步骤4.2）：计算质点r的加速度：

$a_r\left(k\right)=\frac{F_{\mathrm{ext}}^r-F_{\mathrm{in}}^r\left(k\right)}{m_r}$

步骤4.3）：用改进的欧拉算法求解质点r的速度和位移，其中用显式欧拉法对速度v_r进行迭代求解，用隐式欧拉法对位置矢量x_r进行求解：

$v_r^{k+1}=v_r^k+\mathrm{\Δt}\·a_r\left(k\right)$

$x_r^{k+1}=x_r^k+\frac{\mathrm{\Δt}}{2}\·(v_r^k+v_r^{k+1})$

其中，为质点r在时刻k的速度向量，为质点r在时刻k+1的速度向量，本质点r在时刻k的位置向量，为质点r在时刻k+1的位置向量；

步骤4.4）：计算质点r开始运动后与其相连质点间弹簧的长度及弹簧的形变量，并计算质点r所受内力

步骤4.5）：循环执行步骤4.1）至4.4），计算受力区域内其他质点的位置、速度及受力信息。

在上述步骤和方法中，通过有效构建软组织生物力学模型，利用改进的欧拉算法求解软组织形变过程，对其中质点的速度、位移以及受力信息进行高效求解，并实现了光滑反馈力的输出，从而有效解决虚拟手术仿真中的实时性、精确性以及反馈力光滑性问题，提高了软组织形变仿真的精度和实时性，进而满足了虚拟手术仿真的需要。

上面结合附图对本发明优选的具体实施方式和实施例作了详细说明，但是本发明并不限于上述实施方式和实施例，在本领域技术人员所具备的知识范围内，还可以在不脱离本发明构思的前提下作出各种变化。

Claims (1)

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

步骤1)：采集软组织的数据信息，并采用基于四面体的质点-弹簧体模型建立软组织的生物力学模型，该生物力学模型由n个质点组成，并且对于其中任意一个质点i，满足以下方程：

$m_i\frac{{\∂}^2{x}_i}{\∂t^2}+F_{in}^i=F_{ext}^i$

$F_{in}^i={\mathrm{\Σ}}_{j=j1}^{jq}\mathrm{\μ}\left(\right|x_j-x_i|-|x_j-x_i{|}^0)\frac{(x_j-x_i)}{|x_j-x_i|}+{\mathrm{\Σ}}_{j=j1}^{jq}\mathrm{\η}(v_j-v_i)$

$a_i=\frac{F_{ext}^i-F_{in}^i}{m_i}$

其中，m_i为质点i的质量，x_i表示质点i的位置矢量，表示与质点i相连的所有质点j施加于质点i的内力，与质点i相连的所有质点j包括j1至jq，x_j表示质点j的位置矢量，表示质点i所受外力，|x_j-x_i|表示两个质点间位置矢量之差，|x_j-x_i|⁰是发生形变前弹簧的长度，v_i为质点i的速度，v_j为质点j的速度，a_i为质点i的加速度，μ为弹簧的弹性系数，η为阻尼器的阻尼系数；

步骤2)：依据步骤1)所建的生物力学模型，对其中的各个质点进行初始化，包括初始化各个质点的位置、质量、速度、加速度和受力信息，构建模型的初始状态，并且计算该模型中的每个质点与相连质点间的弹簧初始长度；

步骤3)：外接的力反馈设备对软组织施加作用力，进行碰撞检测，确定软组织上发生碰撞的质点和受力发生形变的区域，以及软组织被按压或者拉伸的长度；

步骤4)：用改进的欧拉算法计算发生形变区域内的各个质点的运动状态信息，其中运动状态信息包括发生形变区域内的各个质点的位置、速度和受力随时间变化的信息；计算过程具体为：

步骤4.1)：经力反馈设备对软组织施加外作用力、并经碰撞检测，返回发生碰撞的质点的序列号r及所受外力并且当时刻k＝0时，

步骤4.2)：计算质点r的加速度：

$a_r\left(k\right)=\frac{F_{ext}^r-F_{in}^r\left(k\right)}{m_r}$

步骤4.3)：用改进的欧拉算法求解质点r的速度和位移，其中用显示欧拉法对速度v_r进行迭代求解，用隐式欧拉法对位置矢量x_r进行求解：

$v_r^{k+1}=v_r^k+\mathrm{\Δ}t\·a_r\left(k\right)$

$x_r^{k+1}=x_r^k+\frac{\mathrm{\Δ}t}{2}\·(v_r^k+v_r^{k+1})$

其中，为质点r在时刻k的速度向量，为质点r在时刻k+1的速度向量，为质点r在时刻k的位置向量，为质点r在时刻k+1的位置向量；

步骤4.4)：计算质点r开始运动后与其相连质点间弹簧的长度及弹簧的形变量，并计算质点r所受内力

步骤4.5)：循环执行步骤4.1)至4.4)，计算受力区域内其他质点的位置、速度及受力信息；

步骤5)：依据步骤3)和步骤4)计算反馈力，并将该反馈力输出至力反馈设备；

步骤6)：循环执行步骤4)至步骤5)，计算每个时间步长中各质点的运动状态，并在显示器上动态显示软组织的形变过程。