一种基于硅胶愈合模型的虚拟切割算法的制作方法

本发明属于医疗模拟技术领域，尤其涉及一种基于硅胶愈合模型的虚拟切割算法。

背景技术

近年来，随着虚拟现实技术的不断发展，通过虚拟现实平台模拟手术训练得到可能。虚拟手术给予医护人员极大的便利，他们可以使用该技术进行反复的模拟练习，提升自己的技术。目前的虚拟手术通过建模、渲染、计算将模拟手术的过程呈现在平台上，使用了有限元模型、无网格模型等。然而，目前虚拟手术技术考虑的方面仍不全面，真实性仍旧不够。在进行虚拟切割的过程中，生物软组织对手术刀的摩擦力影响到了手术医师把握施力程度，且身体不同部位，切割难易程度往往不同，这就需要手术医师把握力度从而把握切割深度；生物软组织的自愈能力也是非常重要的一环，根据切口大小、形状、深度的不同，将切口类型大致分为两种：需缝合切口与无需缝合切口，前者需要对切口进行缝补才可愈合，后者可自愈。

技术实现要素：

本发明目的在于克服现有技术的不足，本发明提出一种基于硅胶愈合模型的虚拟切割算法，该方法能够在不降低变形精度的情况下，通过改进后的网格模型更好地实现力的反馈并能使手术医生把控手术力度，所得到的模型更加具有真实性，具体由以下技术方案实现：

所述基于硅胶愈合模型的虚拟切割算法，具体包括如下步骤：

步骤1)建立有限元控制方程；

步骤2)检测切割条件产生切口，并实时产生纳入硅胶愈合模型的网格模型；

步骤3)实时状态检测；

步骤4)实时渲染：对切口进行实时渲染，展现愈合效果。

所述基于硅胶愈合模型的虚拟切割算法的进一步设计在于，所述步骤1)包括如下步骤：

步骤1-1)根据式(1)得到有限元控制方程的矢量化形式：

其中，m是是元素的质量矩阵，u是节点的位移向量，f^ext，f^int，f^damp分别表示单个结点上的外力，内力和阻尼力；

步骤1-2)采用三角形网格作为虚拟血管模型的有限元，引入单个三角单元的局部刚度矩阵如式(2)：

[k^e]＝∫[b]^t[e][b]ds(2)

其中，k^e是单个单元的全局刚度矩阵，b是位移矩阵的应变，e是组织特性的弹性模量矩阵；

步骤1-3)设定三角形单元中每个结点有三个自由度，设定每个结点的位置为xi(i＝1,2,3,4)，有限元的位置矩阵如(3)：

步骤1-4)根据质量集中机制把三角形的质量分配给每个结点，根据式(4)从方程中得到一个三角形的质量：

其中，mi表示第i个三角形的质量，e由与第i个结点相邻的所有三角形有限元组成，ρ是组织材料的密度，se三角形有限元e表示的区域；

步骤1-5)给出系统t0的已知位置x(t0)与速度根据系统的动态方程计算得到系统t0+h的位置x(t0+h)与速度

所述基于硅胶愈合模型的虚拟切割算法的进一步设计在于，所述步骤1-2)中设定血管表面是等距的，进而组织特性的弹性模量矩阵表示为：

其中，λ和μ是lamé常数。

所述基于硅胶愈合模型的虚拟切割算法的进一步设计在于，步骤1-4)zhong1根据式(6)使每个结点上的动态方程(1)离散化：

其中，μi(t)＝[μi(t),υi(t),ωi(t)]^t是第i个结点的位移矢量；

根据式(7)得到第i个结点的内部力为：

根据式(8)得到第i个结点上的阻尼力为：

其中，表示第i个结点的速度；

所述基于硅胶愈合模型的虚拟切割算法的进一步设计在于，步骤1-5)中，首先通过时变的偏微分方程表示模型的推进，再将离散化后的偏微分方程作为一个常微分方程，如式(9)：

其中，向量x表示软模型的位置，对角矩阵mx表示软模型的质量分布，e表示模型的内能，f表示作用于模型的其他力；

再根据牛顿定律，得到系统的动力学方程为：

接着通过将系统的速度υ定义为转化为一阶微分方程：

隐式前向欧拉法将离散动态方程定义为：

式中，δυ表示相邻位置的速度差，δx分别相邻位置的间隔距离；

根据式(13)将一个泰勒级数展开式应用于f并做出一阶近似：

根据式(14)重塑式(13)；

最后计算出δx＝h(υ0+δυ)，从而计算出下个位置x(t0+h)＝x0+δx和速度υ(t0+h)＝υ0+δυ。

所述基于硅胶愈合模型的虚拟切割算法的进一步设计在于，所述步骤2)中，把切割过程分成三个阶段，采用了包含非线性粘弹性的变形模型，

根据式(15)构建平面方程：

ax+by+cz+d＝0(15)

当c≠0时，有：

定义有：

z＝a0x+a1y+a2(17)

设定n个碰撞点pi(xi,yi,zi)，构成的切割面：

当时，式(18)达到最小值，因此得到式(19)与式(20)

求解式(19)与式(20)，得到系数a0,a1,a2，进而得到切割面。

所述基于硅胶愈合模型的虚拟切割算法的进一步设计在于，所述步骤2)中随着虚拟手术刀切割，及时产生加入硅胶愈合模型的网格模型：在网格与手术刀接触面实时产生网格模型，并在网格模型表面加入一层硅胶。

所述基于硅胶愈合模型的虚拟切割算法的进一步设计在于，所述步骤3)中通过硅胶与刀之间的摩擦力模拟生物软组织对手术刀的摩擦力，通过硅胶的粘性模拟生生物软组织的自发愈合现象，具体为：

定义切割过程中产生的摩擦力均为滑动摩擦力，定义：

ff＝μfn1(21)

其中，ff为当前硅胶对手术刀的摩擦力，μ为摩擦因数，由实验得出，表示正压力与产生的摩擦力之间的线性关系，fn1表示当前切割力对网格面的正压力。根据式(22)、式(23)构建力模型：

其中，f1为切割力，为切口顶角的一半，f2为切割力平行于切割平面的分力。设定切口顶端的重叠部分与缝合图中缝合点与切口顶端的重叠部分代表硅胶接触，如式(24)：

α＝f(t)(24)

其中，α表示愈合程度(用百分数表示)，t表示时间，f表示一个分段函数，可表示为：

其中，k1,k2为两个常数，t1为常数，根据所述控制方程，得到相应的位移与速度，进而得到下一时刻的模型。

所述基于硅胶愈合模型的虚拟切割算法的进一步设计在于，所述步骤3)中随着时间增加，愈合程度增加，在切口图中，固定顶角不变，顶点随愈合程度上移，更新模型；在缝合图中，固定缝合点不动，保持缝合点处的顶角与切口顶角不变，切口处顶点随愈合程度上移，更新模型。

本发明的优点如下：

(1)本文方法为现有网格纳入了硅胶愈合模型，其中生物软组织对虚拟手术刀的摩擦力是生物软组织的基本特性，存在摩擦力进行模拟切割更加符合实际，更有利于手术医师把握施力程度。

(2)该方法以硅胶的粘性模拟切口随着时间的自动愈合情况，更加符合实际，适用于缝合等手术。

(3)该方法提出的硅胶愈合模型具有较强的鲁棒性，适用于绝大多数虚拟手术，有利于手术医师在虚拟现实系统中磨练手术技能。

附图说明

图1是切割示意图。

图2是硅胶愈合网格模型图。

图3是图2所示硅胶愈合网格模型图的侧视图。

图4是力模型示意图。

图5是切口示意图。

图6是缝合示意图。

图中，1-网格，2-硅胶，3-缝合点，4-待愈合切口，5-接触面，6-虚拟手术刀。

具体实施方式

结合具体实施例与附图对本发明的技术方案进一步说明。

如图1，本实施例法的基于硅胶愈合模型的虚拟切割算法，其主要步骤包括建立有限元控制方程、检测切割条件产生切口、实时状态检测、实时渲染，具体如下：

步骤1：建立有限元控制方程。

物理模拟是手术模拟器的虚拟世界中的关键部分，基于物理的容积模型为用户提供更真实和交互的性能。然而，由于计算有限元的刚度矩阵或者整个模型降低了模拟系统的流畅性和效率，组织模型的大量数据导致相对复杂的计算问题。msm可以很容易地克服这个问题，并以更有效的方式对体积虚拟模型进行建模，但msm的建模参数并不是真实的组织属性。张量-质量模型在计算复杂度上与msm模型一样有效，但其参数是基于物理的，并且在介入手术模拟器中应用于我们的血管模型。本文基于三角张量-质量模型，结合硅胶愈合模型，提出一种虚拟切割算法。

有限元控制方程的矢量化形式可以表示为：

其中，m是是元素的质量矩阵，u是节点的位移向量，f^ext，f^int，f^damp分别表示单个结点上的外力，内力和阻尼力。f^ext和f^damp以下述形式派生：

其中，k是全局刚度矩阵，c是阻尼力矩阵，根据rayleigh阻尼方程，两者的关系被表示为：

{c}＝α[m]+β[k](3)

其中，α，β是rayleigh系数。

tmm的关键思想是简单地结合局部单元刚度矩阵形成全局刚度矩阵从而打破传统的约束。tmm方法将刚度矩阵分解为节点分量和边缘分量，并将两个分量分别累加到节点和边上。本文采用三角形网格作为虚拟血管模型的有限元，参见图2、图3，引入单个单元的局部刚度矩阵：

[k^e]＝∫[b]^t[e][b]ds(4)

其中，k^e是单个单元的全局刚度矩阵，b是位移矩阵的应变，e是组织特性的弹性模量矩阵(杨氏模量)。假设血管表面是等距的，则软组织的性质在任何方向上都具有相同的性能。有限元的杨氏模量矩阵表示为：

其中λ和μ是lamé常数(拉梅常数)，是应力-应变关系中的两个材料相关量。

在一个三角形单元中有三个结点，每个结点有三个自由度(dof)。设定每个结点的位置为xi(i＝1,2,3,4)，有限元的位置矩阵可以表示为：

将应变和位移矩阵b表示为下述数组：

[b]＝[b1b2b3](7)

将(7)带入(4)，得到3×3的局部单元刚度矩阵的阵列，表示为：

其中，是tmm模型中的张量，表示第i个结点，表示边ij，i和j是一个三角形有限元中结点的局部指数，在tmm力学模型中，局部三角形有限元的张量组合成一个全局张量，表示为kij，i,j是全局定义的。

另外一个机制是质量集中，它把三角形的质量分配给每个结点。从方程中得到一个三角形的质量：

其中，mi表示第i个三角形的质量，e由与第i个结点相邻的所有三角形有限元组成，ρ是组织材料的密度，se三角形有限元e表示的区域。用质量集中法得到的质量矩阵m是关于对角线对称的。使每个结点上的动态方程(1)离散化：

其中，μi(t)＝[μi(t),υi(t),ωi(t)]^t是第i个结点的位移矢量。第i个结点的内部力为：

第i个结点上的阻尼力为：

其中，表示第i个结点的速度。

大部分将模拟器向前推进的方法都有一个共同点，可以表示为一个时变的偏微分方程，它在离散化后作为一个常微分方程被求解：

其中，向量x和对角矩阵m表示软模型的位置和质量分布，e表示模型的内能，f表示作用于模型的其他力。

给出系统t0的已知位置x(t0)与速度目的是得到系统t0+h的位置x(t0+h)与速度为了使用隐式技术来计算新的状态，需要求解系统的动态方程。根据牛顿定律，系统的动力学方程可表示为：

通过将系统的速度υ定义为可以将其转化为一阶微分方程：

隐式前向欧拉法将离散动态方程定义为：

将一个泰勒级数展开式应用于f并做出一阶近似：

用(19)重塑(18)得到：

计算出δx＝h(υ0+δυ)，从而计算出下个位置x(t0+h)＝x0+δx和速度υ(t0+h)＝υ0+δυ。

步骤2)检测切割条件产生切口(如图1)，并实时产生纳入硅胶愈合模型的网格模型(如图2、图3)。

在真正的切割中，在手术刀穿透软组织的过程中发生以下两个重要的相互作用。首先，当切割力增加时，软组织的变形不断发生。直到施加到软组织上的切割力超过阈值水平才会发生分离。其次，当施加在软组织上的应力达到断裂应力时，软组织将被切开，并随着切割力的急剧下降而开始快速膨胀。一旦软组织被切开，软组织将以与切割平面中心相反的方向变形。因此，可以把切割过程分成三个阶段。采用了包含非线性粘弹性的变形模型，更切合实际地表示了变形。

手术刀通常只是沿着刀片的方向切割软组织，所以在计算软组织变形之前，必须确定切割面。由于切割面是虚构的，使用碰撞检测得到的点来适应切割面。在几何中，一个平面的方程可以表示为：

ax+by+cz+d＝0(20)

当c≠0时，有：

定义有：

z＝a0x+a1y+a2(22)

假设有n个碰撞点pi(xi,yi,zi)，这些点构成的切割面符合要求：

当时，上式达到最小值。因此：

求解上式，可以得到系数a0,a1,a2，进而得到切割面。

随着虚拟手术刀切割，及时产生加入硅胶愈合模型的网格模型：在网格-手术刀接触面实时产生网格模型，并在网格模型表面加入一层硅胶。

步骤3：实时状态检测。

生物软组织在手术刀切割时，会对手术刀产生一个较小的摩擦力；而切割产生的开口因生物软组织的粘弹性等会自发局部愈合。在虚拟手术系统中加入这两点，能更好的模拟真实，使手术医师能更好的磨练技术。

本文通过硅胶与刀之间的摩擦力，硅胶的粘性分别模拟生物软组织对手术刀的摩擦力与生物软组织的自发愈合现象。

考虑切割过程中产生的摩擦力均为滑动摩擦力，不考虑无法切割而产生静摩擦力的情况，定义：

ff＝μfn1(26)

其中，ff为当前硅胶对手术刀的摩擦力，μ为摩擦因数，由实验得出，表示正压力与产生的摩擦力之间的线性关系，fn1表示当前切割力对网格面的正压力。

如图4的力模型图，其中：

其中，f1为切割力，为切口顶角的一半，f2为切割力平行于切割平面的分力。

切口图(参见图5)切口顶端的重叠部分与缝合图(参见图6)中缝合点与切口顶端的重叠部分代表硅胶接触，定义：

α＝f(t)(29)

其中，α表示愈合程度(用百分数表示)，t表示时间，f表示一个分段函数，可表示为：

其中，k1,k2为两个常数，由实验得出，t1为常数，此时愈合程度为100％。结合步骤一建立的控制方程，可以得到相应的位移与速度，进而得到下一时刻的模型。

随着时间增加，愈合程度增加，在切口图中，固定顶角不变，顶点随愈合程度上移，更新模型；在缝合图中，固定缝合点不动，保持缝合点处的顶角与切口顶角不变，切口处顶点随愈合程度上移，更新模型。

步骤4：实时渲染，对切口进行实时渲染，展现愈合效果。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。