CN109616211A - 一种基于质点弹簧的混合缝合线模型及其打结方法 - Google Patents

Abstract

一种基于质点弹簧的混合缝合线模型及其打结方法，该模型在质点弹簧模型的基础上加入了扭转弹簧与分段计算缝合线张力的方法，并将图元碰撞检测以及碰撞簇方法结合，还加入了对单位时间步内缝合线顶点的最大移动距离的限制。这种混合缝合线模型及其打结方法可以真实的模拟缝合线内部张力，有效避免了误检、漏检碰撞的现象，并提高了碰撞反馈力计算的准确性。实验表明，与传统的FTL缝合线建模方法相比，该模型及其打结方法在视觉渲染上更具真实性，且在处理“称人结”这类复杂结模拟中更具稳定性。

Description

一种基于质点弹簧的混合缝合线模型及其打结方法

技术领域

本发明涉及虚拟手术仿真***中缝合线模型的构建，用于实现在虚拟手术中对伤口缝合与打结过程的模拟。

背景技术

缝合线是一种质量非常小，对拉伸形变有很强的抵抗力，却又几乎不能抵抗弯曲形变的可形变物体(Deformable Linear Objects,DLO)，通常使用样条曲线来进行模拟，一些模拟***中为了简化计算直接将缝合线渲染成一维的线段。良好的缝合线模型是虚拟手术仿真***的基础。对缝合线的模拟不仅需要合适的物理模型，还需要真实的三维渲染效果；同时，准确的计算缝合线内部的张力也是其与软组织交互的基础。此外，缝合线的打结模拟需要精细的自碰撞检测及反馈方法来防止漏检、误检碰撞。

随着计算机硬件水平的不断发展，三维模型的缝合线越来越多的出现在当今的模拟***中，常见的缝合线模拟方法有：FTL(Follow The Leader)方法，质点弹簧模型(Mass-Spring Model,MSM)和在Cosserat理论基础上建立的“曲线-角度”表示的离散弹性杆模型。在早期，Phillips等通过使用一系列彼此重叠的球体来模拟了一条绳子，当其被拉伸时会添加新的球体来保证绳子的完整性，并在收缩时去掉之前填充的球体。后来，Brown等提出一种基于FTL方法的缝合线模型，它是由有限个首尾相连长度相等(且在实验过程中不会发生变化)的刚性杆构成，杆的连接处可以自由转动。FTL模型在处理缝合线自碰撞过程中，当两个刚性杆相交时，对应的顶点在当前时间步长内固定，下个时间步长推开，来模拟摩擦力的效果。由于该算法是一种基于几何模型的建模方法，它并不能模拟缝合线的内部张力的变化。

不同于FTL方法构建的缝合线，在使用质点弹簧模型构建的缝合线中，每段缝合线的长度是可变的线性弹簧，并通过线性弹簧长度的变化来模拟缝合线内部的张力。为了实现更真实的缝合线模型，Payandeh等人把扭转弹簧加入到缝合线模型，同时还考虑了空气阻力对缝合线的作用力，在处理缝合线与软组织的交互过程中，还提出了使用置换算法来判断某个缝合线顶点位于软组织内部还是外部。在对于缝合线张力的计算上，Lenoir等人还考虑了在缝合线穿过软组织之后进行分段处理，同时在分别计算每段缝合线的张力时还可以加入在每个穿刺点处缝合线与软组织之间的摩擦力，从而获得更准确的力反馈效果。这种分段处理缝合线张力的方法也成为缝合线模拟的常用方法。2013年吕梦雅等人使用质点弹簧模型构建了缝合线，在对缝合线顶点的位置更新中，该模型借鉴了FTL算法的更新方法，从而实现了缝合线与软组织之间的力学交互过程的快速实现。但是并没有提供相应的自碰撞检测处理及打结的实现。

Cosserat弹性杆理论首次被应用到虚拟手术中是由Pai等人进行的缝合线模拟。然而对于弹性杆中心线的隐式法大大增加了其在碰撞处理上的计算复杂度。Spillmann等人在此基础上进行了改进，他们使用显示法处理了弹性杆的中心线，使用隐式法来处理弹性杆的横截面方向。Punak等人通过省略动力学效应以及内部摩擦和其他耗散效应来优化性能。这些计算上的简化不仅保证了模型的相对准确，还提高了计算速度，实现了Cosserat弹性杆模型在精度和计算成本之间良好的折衷，但过高的计算复杂度仍然使得该方法很难满足虚拟手术缝合模拟的实时交互操作，特别是在处理打结模拟时将会显著增加计算复杂度。Bergou等人在Cosserat理论的基础上建立了“曲线-角度”表示的离散弹性杆模型，将需要模拟的模型分解表示为无扭转的Bishop杆和关于弹性杆中心切向的旋转角，该方法既避免了变量冗余，又大大简化了对弯曲-扭转耦合的处理。但这种模型在处理约束的时候会出现一些能量泄漏导致在计算最小能量时不准确。

综上所述，传统的基于FTL(follow the leader)方法虽然简单易实现且可以满足缝合模拟中实时性的要求，但无法模拟缝合线内部的张力变化，而缝合线的内部张力是计算缝合线与软组织之间的摩擦力的基础，仅通过人为的添加一些固定的张力无法满足缝合模拟过程的真实性需求。Cosserat弹性杆模型存在过高的计算复杂度使缝合模拟的实时性降低的问题。相较而言，采用质点弹簧模型可以构建出由若干段长度可变的线性弹簧连接而成的缝合线，通过线性弹簧长度的变化来模拟缝合线内部的张力，准确的张力计算不仅提高了缝合线模型的真实性还有利于模拟缝合线与软组织之间的摩擦力。

发明内容

基于上述背景，本发明提出了一种基于质点弹簧的缝合线混合模型及其打结方法。在模型的构建上，该模型在质点弹簧模型的基础上加入了扭转弹簧与分段计算缝合线张力的方法，在对缝合线打结过程的处理上，将图元碰撞检测以及碰撞簇方法结合，并加入了对单位时间步内缝合线顶点的最大移动距离的限制。

本发明是通过以下技术方案实现的。

本发明所述的一种基于质点弹簧的混合缝合线模型及其打结方法，按以下步骤：

在虚拟手术缝合仿真***中，有一根长度为L、半径为R的缝合线，根据质点弹簧建模原理，将其视为由n段原长为l的线性弹簧首尾相连而成，相邻两弹簧的公共端点为缝合线的顶点，将N个顶点记为P_i，i＝1,2,3……N,N＝n+1。

步骤1：首先根据牛顿运动第二定律构建任意顶点P_i的动力学方程，并使用欧拉有限差分方法对上述方程求解，从而获得缝合线N个顶点的位置矢量rⁱ，在每个时间步长内对所有顶点位置信息更新，获取整条缝合线的运动状态；

步骤2：分段计算缝合线的张力：如附图1，首先从缝合线头部对每段缝合线的长度进行检测，检测到缝合线段p_i+1p_i+2没有发生拉伸，该段缝合线张力大小为 且缝合线顶点p_i+1、p_i+2受到该段缝合线的作用力为0，由此可知缝合线段p_ip_i+1上的张力大小为f_i沿缝合线方向的分力大小，即:则缝合线段p_i-1p_i的张力大小为f_i-1、f_i在该段缝合线方向上的合力大小，即： 通过这种方法一直推导到第一段缝合线p₁p₂的张力，同理向顶点p_i+3方向使用同样的方法计算出最后一段缝合线p_np_n+1的张力；

步骤3:计算扭转弹簧弹力与空气阻力：如附图2，除了缝合线p₁、p_n+1两个顶点外，把其余每个顶点都看作一个夹角的顶点，在这个夹角的两条边下形成一个弯曲弹簧模型，使用这种方法，正交于该弹簧段的力的分量将被用在每个夹角的顶点处。计算角度的公式如下：

再计算扭转弹簧力f_i-1,f_i的方向向量n_i-1,n_i，即:

最后计算点P_i-1、P_i+1受到的扭转弹簧力f_i-1、f_i-1分别为：

其中v_i-1,v_i,v_i+1分别是点P_i-1,P_i,P_i+1处的速度，k_td是扭转弹簧力系数,则点P_i处的扭转弹簧力f_i的计算公式为:f_i＝-(f_i-1+f_i+1)，此外，记空气阻力为f_air，其大小是根据每个缝合线顶点的速度乘上相应的空气阻力系数k_air，其方向与该缝合线顶点的运动方向相反。

步骤4：缝合线的立体渲染：如附图3，以缝合线p_ip_i+1为轴心，i＝1,2,……，N,N＝n+1，构建一个高为L、半径为R的圆柱体。再以每个缝合线顶点p_i为球心构建同等半径的小球连接两段相邻的圆柱体，以此填补缝合线弯折时在两段缝合线连接处视觉上产生的空隙，并使用Catmull-Rom样条曲线插值法平滑线条；

步骤5：缝合线的自碰撞检测及反馈：本发明将两段缝合线的碰撞检测视为两个圆柱体之间的碰撞检测，如附图4所示。在时间步长dt内，对于两个半径为r的柱状缝合线p_ip_i+1,p_jp_j+1，首先需要计算出两段柱状缝合线中心线之间的距离d，如果d<2r，则认为缝合线发生了自碰撞，否则没有发生碰撞。设缝合线段p_ip_i+1,p_jp_j+1发生了碰撞，首先需要计算出p_ip_i+1,p_jp_j+1到彼此最短距离的点c_i,c_j，然后计算出从点c_i到点c_j的单位向量n，对顶点i,i+1,j,j+1的位置进行如下调整：

其中D为位移距离，计算公式为：D＝2*r-d+α，其中α为摩擦力替代系数，是一个可调整的非负数，来消除一些因插值操作而产生的误差，同时α等效于摩擦力的存在，α越大则等效摩擦力越小，若某段缝合线在一个时间步内同时与多段缝合线发生碰撞，则其最终位置为受到多次碰撞反馈位置结果的平均值，即D＝(D₁+D₂+......+D_n)/n；

步骤6：改进的穿透深度算法，针对穿透深度法在处理缝合线这类纤薄物体时，由于自身没有太大的深度，会出现排斥力计算不准确的问题，本发明把单位时间步长内缝合线顶点移动的最大距离D_max设置为小于等于r/4的值能保证在自碰撞反馈中排斥力的准确计算，本发明中取D_max≤0.5mm；

步骤7：碰撞簇方法，在碰撞检测遍历过程中，将已经检测到发生碰撞的两段缝合线所对应的四个顶点归为一个碰撞簇，若一个缝合线属于两个碰撞簇，则将两个碰撞簇合并，每次完成碰撞检测遍历后统计每个碰撞簇的顶点个数，若顶点个数到达一定阈值，则将该碰撞簇作为一个缝合结来看待，且不对这些顶点进行碰撞检测处理。

本发明可以真实的模拟缝合线内部张力，有效避免了误检、漏检碰撞的现象，并提高了碰撞反馈力计算的准确性。实验表明，与传统的FTL缝合线建模方法相比，该模型在视觉渲染上更具真实性，且在处理“称人结”这类复杂结模拟中更具稳定性。

附图说明

图1为缝合线的张力分析图。

图2为扭转弹簧受力分析图。

图3为质点弹簧模型模拟的缝合线立体渲染图。

图4为缝合线自碰撞示意图。

图5为“8字结”真实情况、FTL模型、混合模型打结对比。其中，(a)为“8字结”真实打结图，(b)为基于FTL模型打结图，(c)为混合模型打结图。

图6为“称人结”真实情况、混合模型打结对比。其中，(a)为真实打结图(b)为混合模型打结图。由于FTL模型无法完成“称人结”这类复杂结的处理，故没有其效果图。

具体实施方式

本发明将通过以下实例作进一步说明。

在虚拟手术缝合仿真***中，有一根长度为0.4m、半径为2mm的缝合线，根据质点弹簧建模原理，将其视为由80段原长为5mm的线性弹簧首尾相连而成，设相邻两弹簧的公共端点为缝合线的顶点，将N个顶点记为P_i，i＝1,2,3……N,N＝81。

步骤1：首先根据牛顿运动第二定律构建任意顶点P_i的动力学方程，并使用欧拉有限差分方法对上述方程求解，从而获得缝合线N个顶点的位置矢量rⁱ，在每个时间步长内对所有顶点位置信息更新，获取整条缝合线的运动状态；

步骤2：分段计算缝合线的张力，如附图1，首先从缝合线头部对每段缝合线的长度进行检测，检测到缝合线段p_i+1p_i+2没有发生拉伸，该段缝合线张力大小为 且缝合线顶点p_i+1、p_i+2受到该段缝合线的作用力为0，由此可知缝合线段p_ip_i+1上的张力大小为f_i沿缝合线方向的分力大小，即:则缝合线段p_i-1p_i的张力大小为f_i-1、f_i在该段缝合线方向上的合力大小，即： 通过这种方法一直推导到第一段缝合线p₁p₂的张力，同理向顶点p_i+3方向使用同样的方法计算出最后一段缝合线p_np_n+1的张力。

步骤3:计算扭转弹簧弹力与空气阻力，如附图2，除了缝合线p₁、p_n+1两个顶点外，把其余每个顶点都看作一个夹角的顶点，在这个夹角的两条边下形成一个弯曲弹簧模型，使用这种方法，正交于该弹簧段的力的分量将被用在每个夹角的顶点处。计算角度的公式如下：

再计算扭转弹簧力f_i-1,f_i的方向向量n_i-1,n_i，即:

最后计算点P_i-1、P_i+1受到的扭转弹簧力f_i-1、f_i-1分别为：

其中v_i-1,v_i,v_i+1分别是点P_i-1,P_i,P_i+1处的速度，k_td是扭转弹簧力系数,则点P_i处的扭转弹簧力f_i的计算公式为:f_i＝-(f_i-1+f_i+1)，此外，记空气阻力为f_air，其大小是根据每个缝合线顶点的速度乘上相应的空气阻力系数k_air，其方向与该缝合线顶点的运动方向相反,本发明中k_air取值为0.2Ns/m。

步骤4：缝合线的立体渲染，如附图3，以缝合线p_ip_i+1为轴心，i＝1,2,……，N,N＝n+1，构建一个高为L、半径为R的圆柱体。再以每个缝合线顶点p_i为球心构建同等半径的小球连接两段相邻的圆柱体，以此填补缝合线弯折时在两段缝合线连接处视觉上产生的空隙，并使用Catmull-Rom样条曲线插值法平滑线条。

步骤5：缝合线的自碰撞检测及反馈，本发明将两段缝合线的碰撞检测视为两个圆柱体之间的碰撞检测，如附图4所示。在时间步长dt内，对于两个半径为r的柱状缝合线p_ip_i+1,p_jp_j+1，首先需要计算出两段柱状缝合线中心线之间的距离d，如果d<2r，则认为缝合线发生了自碰撞，否则没有发生碰撞。设缝合线段p_ip_i+1,p_jp_j+1发生了碰撞，首先需要计算出p_ip_i+1,p_jp_j+1到彼此最短距离的点c_i,c_j，然后计算出从点c_i到点c_j的单位向量n，对顶点i,i+1,j,j+1的位置进行如下调整：

其中D为位移距离，计算公式为：D＝2*r-d+α，其中α为摩擦力替代系数，是一个可调整的非负数，来消除一些因插值操作而产生的误差，同时α等效于摩擦力的存在，α越大则等效摩擦力越小，若某段缝合线在一个时间步内同时与多段缝合线发生碰撞，则其最终位置为受到多次碰撞反馈位置结果的平均值，即D＝(D₁+D₂+......+D_n)/n。

步骤6：改进的穿透深度算法，针对穿透深度法在处理缝合线这类纤薄物体时，由于自身没有太大的深度，会出现排斥力计算不准确的问题，本发明把单位时间步长内缝合线顶点移动的最大距离D_max设置为小于等于r/4的值将保证在自碰撞反馈中排斥力的准确计算，本发明中取D_max≤0.5mm。

步骤7：碰撞簇方法，在碰撞检测遍历过程中，将已经检测到发生碰撞的两段缝合线所对应的四个顶点归为一个碰撞簇，若一个缝合线属于两个碰撞簇，则将两个碰撞簇合并，每次完成碰撞检测遍历后统计每个碰撞簇的顶点个数，若顶点个数到达一定阈值，则将该碰撞簇作为一个缝合结来看待，且不对这些顶点进行碰撞检测处理。

Claims (1)

1.一种基于质点弹簧的混合缝合线模型及其打结方法，其特征是按以下步骤：

在虚拟手术缝合仿真***中，一根长度为L、半径为R的缝合线，根据质点弹簧建模原理，将其视为由n段原长为l的线性弹簧首尾相连而成，相邻两弹簧的公共端点为缝合线的顶点，将N个顶点记为P_i，i＝1,2,3……N,N＝n+1；

步骤1：首先根据牛顿运动第二定律构建任意顶点P_i的动力学方程，并使用欧拉有限差分方法对上述方程求解，从而获得缝合线N个顶点的位置矢量rⁱ，在每个时间步长内对所有顶点位置信息更新，获取整条缝合线的运动状态；

步骤2：分段计算缝合线的张力：首先从缝合线头部对每段缝合线的长度进行检测，检测到缝合线段p_i+1p_i+2没有发生拉伸，该段缝合线张力大小为且缝合线顶点p_i+1、p_i+2受到该段缝合线的作用力为0，由此可知缝合线段p_ip_i+1上的张力大小为f_i沿缝合线方向的分力大小，即：则缝合线段p_i-1p_i的张力大小为f_i-1、f_i在该段缝合线方向上的合力大小，即：通过这种方法一直推导到第一段缝合线p₁p₂的张力，同理向顶点p_i+3方向使用同样的方法计算出最后一段缝合线p_np_n+1的张力；

步骤3：计算扭转弹簧弹力与空气阻力：除了缝合线p₁、p_n+1两个顶点外，把其余每个顶点都看作一个夹角的顶点，在这个夹角的两条边下形成一个弯曲弹簧模型，正交于该弹簧段的力的分量将被用在每个夹角的顶点处，计算角度的公式如下：

再计算扭转弹簧力f_i-1,f_i的方向向量n_i-1,n_i，即：

最后计算点P_i-1、P_i+1受到的扭转弹簧力f_i-1、f_i-1分别为：

其中v_i-1,v_i,v_i+1分别是点P_i-1,P_i,P_i+1处的速度，k_td是扭转弹簧力系数，则点P_i处的扭转弹簧力f_i的计算公式为:f_i＝-(f_i-1+f_i+1)，此外，记空气阻力为f_air，其大小是根据每个缝合线顶点的速度乘上相应的空气阻力系数k_air，其方向与该缝合线顶点的运动方向相反；

步骤4：缝合线的立体渲染：以缝合线p_ip_i+1为轴心，i＝1,2,……，N,N＝n+1，构建一个高为L、半径为R的圆柱体，再以每个缝合线顶点p_i为球心构建同等半径的小球连接两段相邻的圆柱体，以此填补缝合线弯折时在两段缝合线连接处视觉上产生的空隙，并使用Catmull-Rom样条曲线插值法平滑线条；

步骤5：缝合线的自碰撞检测及反馈：将两段缝合线的碰撞检测视为两个圆柱体之间的碰撞检测在时间步长dt内，对于两个半径为r的柱状缝合线p_ip_i+1，p_jp_j+1，首先需要计算出两段柱状缝合线中心线之间的距离d，如果d<2r，则认为缝合线发生了自碰撞，否则没有发生碰撞，设缝合线段p_ip_i+1,p_jp_j+1发生了碰撞，首先需要计算出p_ip_i+1，p_jp_j+1到彼此最短距离的点c_i，c_j，然后计算出从点c_i到点c_j的单位向量n，对顶点i，i+1，j，j+1的位置进行如下调整：

其中D为位移距离，计算公式为：D＝2*r-d+α，其中α为摩擦力替代系数，是一个可调整的非负数，来消除一些因插值操作而产生的误差，同时α等效于摩擦力的存在，α越大则等效摩擦力越小，若某段缝合线在一个时间步内同时与多段缝合线发生碰撞，则其最终位置为受到多次碰撞反馈位置结果的平均值，即D＝(D₁+D₂+......+D_n)/n；

步骤6：改进的穿透深度算法：针对穿透深度法在处理缝合线这类纤薄物体时，由于自身没有太大的深度，会出现排斥力计算不准确的问题，把单位时间步长内缝合线顶点移动的最大距离D_max设置为小于等于r/4的值能保证在自碰撞反馈中排斥力的准确计算；

步骤7：碰撞簇方法：在碰撞检测遍历过程中，将已经检测到发生碰撞的两段缝合线所对应的四个顶点归为一个碰撞簇，若一个缝合线属于两个碰撞簇，则将两个碰撞簇合并，每次完成碰撞检测遍历后统计每个碰撞簇的顶点个数，若顶点个数到达一定阈值，则将该碰撞簇作为一个缝合结来看待，且不再对这些顶点进行碰撞检测处理。