本发明涉及虚拟切割，尤其是基于纳入开尔文粘弹性模型的网格模型的虚拟切割方法。
背景技术
近年来，随着虚拟现实技术的不断发展，通过虚拟现实进行模拟手术得到可能。虚拟手术给予医护人员极大的便利，他们可以使用该技术进行反复的模拟练习，提升自己的技术。目前的虚拟手术通过建模、渲染、计算将模拟手术的过程呈现在平台上，使用了有限元模型、无网格模型等。然而，有限元模型、无网格模型都有一些缺点。对于有限元模型，它深深依赖于网格，而扭曲的或者低质量的网格会造成很大的误差。在重新啮合过程中，产生的失真元素甚至会导致模型不稳定。相比于经典的基于网格模型的基础结构，它不适合模拟切割网格结构和连续性的切割。而无网格模型的目的是克服与有限元模型相关的问题。与有限元模型相比，无网格模型在离散和分离点元素的基础上重建虚拟软组织，并且每个点元素之间的关系不与网格关联。因此，点元素是随机的，且不受网格约束，适用于不连续的场景。虽然上述无网格模型在模拟软组织的切割过程中是有前景的，但是虚拟手术器械与软组织间的相互作用成为了一个问题。为了简化仿真过程，大多数方法都考虑只要软组织被虚拟手术刀扫过，组织便被分开。然而，结果并不像预期的那么简单。在切割过程中，可以观察到在软组织切开之前可能出现明显的变形。正如我们所能理解的，不同的器械与软组织相互作用通常会导致不同形式的变形。此外，有限元模型、无网格模型所对应的方法的计算量均较大，且较为复杂。
技术实现要素：
发明目的：针对上述现有技术存在的缺陷，本发明旨在提供一种基于纳入开尔文粘弹性模型的网格模型的虚拟切割方法，利用纳入开尔文粘弹性模型的网格提高虚拟变形以及切割的效率。
技术方案：一种基于纳入开尔文粘弹性模型的网格模型的虚拟切割方法，构建开尔文粘弹性模型，所述开尔文粘弹性模型的一般离散控制方程为：
其中，kn是全局刚度矩阵，
施加滞纳函数的全局滞纳刚度矩阵：
其中b，φ，δ∈n中的元素均为常量；δ是虚拟切割材料的参数，b是常数，
由虚拟切割材料的参数和切割时间求解[tn，tn+1]时间中的位移增量
进一步的，所述计算每个节点的新位移、应变和应力具体为依据[tn，tn+1]时间中的位移增量
在tn+1时，位移、应力、应变的增量依次为：
σn+1＝σn+δσn；
εn+1＝εn+δεn；
其中，
其中l为常数；
其中，开尔文粘弹性模型在[tn，tn+1]中的应力增量为：δσn＝δεnεk+σ0，n
；其中，应力与应变的关系为：
其中，εk是线性松弛系数，表示在时间间隔[tn，tn+1]中由单位阶跃应变增量引起的应力变化，c0，c1，τ1为材料参数；
tn+1时的初始应力为：
进一步的，所述开尔文粘弹性模型的本构方程是：
其中，σ1表示应力，η是阻尼器的阻尼系数，σ2表示应力的时间导数，ε2和ε1分别表示两个弹簧的刚性，ε1表示应变，ε2表示应变的时间导数。
进一步的，所述开尔文粘弹性模型中应变和应力的本构关系是：
其中，σ表示应力，ε表示弹性系数(杨氏模型)，ε表示应变，c0和c1是材料参数，t表示时间，τ1是时间常数。
进一步的，所述构建开尔文粘弹性模型之前还包括受力判断步骤：记碰撞区域为a，设定阈值f1，d1和d2，若力小于f1，或a的宽度大于等于d2，则该网格模型仅发生形变；若力大于等于f1，且a的宽度小于d1且大于0，则属于切割情况一；若力大于等于f1，且a的宽度大于等于d1且小于d2，则属于切割情况二；
所述切割情况一具体为：将a视为一条直线l1，读取端点位置，将两个端点设为刚性核，固定不动，复制l1，得到l1、l2，将复制前连在l1上的网格线与a的a上的左侧交点连在l1上，把右侧交点连在l2上，l1、l2因为受力不平衡根据胡克定律f＝-kδx产生弹性形变，其中k为弹簧劲度系数，仅考虑水平方向受力；
所述切割情况二具体为：将a视为矩形区域，固定宽不动，删除a内的网格线，两条长边因受力不平衡根据胡克定律f＝-kδx产生弹性形变，仅考虑水平方向受力。
有益效果：本发明为网格纳入了开尔文粘弹性模型，粘弹性是生物软组织的基本特性，使软组织具有粘弹性进行模拟切割更加符合实际，且更有利于后续力的反馈的设定。本发明方法立足高质量的网格模型，以受力变形代替近似计算，解决了有限元模型中的网格失真、连续切割的问题与无网格模型中有关不同变形形式的问题。且本发明提出的网格模型以力的平衡取代部分点的位置的计算，大大减少了计算量，优化了性能。通过此模型，针对变形，只需计算施力时点的位置，而撤力后，由于力的平衡被打破，网格自动复原，减少了计算量；针对切割，剔除开口间的具有粘弹性的网格线，打破力的平衡状态，使网格自动变形，相比于传统的近似开口形状、大小，更符合实际，且操作简单。由于变形与切割更加符合实际，后期渲染得到的模型效果卓越。
附图说明
图1是本发明的流程图；
图2是开尔文粘弹性模型；
图3是网格模型；
图4是切割情况一的示意图；
图5是切割情况二的示意图。
具体实施方式
下面通过一个最佳实施例并结合附图对本技术方案进行详细说明。
如图1所示，一种基于纳入开尔文粘弹性模型的网格模型的虚拟切割方法，其主要步骤包括为网格纳入开尔文粘弹性模型、粘弹性的结合、受力并判断相关条件、切割，具体如下：
步骤1：为网格纳入开尔文粘弹性模型。
现实生活中的软组织有滞后，松弛和蠕变的性质，这些性质统称为粘弹性。该粘弹性机制模型可以描述软组织的粘弹性，这在生物学特性中非常重要。可以通过修改在体内实验中可以获得的相关参数来描述不同软组织的生物力学特征。本文使用的开尔文粘弹性模型是一种标准线性模型，其结构如图2所示；
模型中的弹簧表示软组织的线性弹性特征，阻尼器表示软组织变化时的阻尼特性。粘弹性力学的几何与运动方程与弹性力学中的几何与运动方程相同。粘弹性边界值的解可以通过求解运动方程，几何方程，边界条件和初始条件的本构方程得到；
该开尔文粘弹性模型的本构方程是：
其中，σ1表示应力，η是阻尼器的阻尼系数，σ2表示应力的时间导数，ε2和ε1分别表示两个弹簧的刚性，ε1表示应变，ε2表示应变的时间导数；
应变和应力的本构关系是：
其中，σ表示应力，ε表示弹性系数(杨氏模型)，ε表示应变，c0和c1是材料参数，t表示时间，τ1是时间常数；
步骤2：粘弹性的结合。
为了将粘弹性纳入到变形模型中，使用增量形式的粘弹性模型。首先，将变形模拟时间t分为n个时间片t1，t2，...，tn。每个时间间隔
在开尔文粘弹性模型中，软组织处于外力下的情况的松弛度可以通过松弛本构关系表示：
当处于tn与tn+1，且δt→0时，应力增量为：
其中，ε表示弹簧的刚性，τ表示时间常数；
开尔文粘弹性模型的应力与应变关系由(2)表示，结合(4)得到：
其中εk是线性松弛系数，表示在时间间隔[tn，tn+1]中由单位阶跃应变增量引起的应力变化。c0，c1，τ1为材料参数。通过下式得到tn+1时的应力σ0，n：
通过下式获得开尔文粘弹性模型在[tn，tn+1]中的应力增量：
δσn＝δεnεk+σ0，n(7)
在tn+1时，位移，应力，应变的增量分别是：
σn+1＝σn+δσn(9)
εn+1＝εn+δεn(10)
其中，
其中l为常数。在软组织变形过程中，应力、体力和外力满足平衡条件，即根据变形的虚拟工作原理，总虚拟作业为零，可得到粘弹性模型的一般离散控制方程：
其中，kn是全局刚度矩阵，
其中，材料是均质的，b，φ，δ∈n中的元素均为常量。δ是材料参数，b是常数，
步骤3：受力并判断相关条件。
网格模型如图3，记碰撞区域为a，设定阈值f1，d1和d2，若力小于f1，或a的宽度大于等于d2，则该网格模型仅发生形变；若力大于等于f1，且a的宽度小于d1且大于0，则属于切割情况一(如图4)；若力大于等于f1，且a的宽度大于等于d1且小于d2，则属于切割情况二(如图5)；
步骤4：变形。
根据受力情况，完成粘弹性结合部分的算法，计算位移；
步骤5：切割。
当情况属于情况1或情况2时，该网格模型将产生切口；
5-1：切割情况一，将a视为一条直线l1，读取端点位置，将两个端点设为刚性核，固定不动，复制l1，得到l1、l2，将复制前连在l1上的网格线与a的a上的左侧交点连在l1上，把右侧交点连在l2上，这样，l1、l2因为受力不平衡根据胡克定律f＝-kδx(其中k为弹簧劲度系数)，产生弹性形变，仅考虑水平方向受力，最终达到受力平衡状态，从而产生切口；
5-2：切割情况二，将a视为矩形区域，固定宽不动，删除a内的网格线，从而两条长因受力不平衡根据胡克定律产生弹性形变，仅考虑水平方向受力，最终达到受力平衡状态，从而产生切口。