一种基于位置约束和非线性弹簧的软组织建模方法
技术领域
本发明属于软组织形变模型建模
技术领域
，具体涉及一种基于位置约束和非线性弹簧的软组织建模方法，用于解决虚拟手术实验中软组织建模问题。背景技术
虚拟手术是一种可以在低成本模式下逼真地模拟手术环境的重要的外科手术辅助训练系统。虚拟手术不仅能为使用者提供可视化需求，还能使医生熟悉手术器械的使用，积累手术技能。目前，一个完善的虚拟手术系统应该同时具备三个需求：实时性，精确性，稳定性。软组织建模是虚拟手术中的核心模块之一，模型的精细程度对系统的可靠性起着决定性的作用。
软组织建模方法主要包括几何建模法和物理建模法。几何建模法通过对医学图像进行三维重建重绘软组织的几何模型，而物理建模法则是根据软组织的生物力学和运动学特性建立物理模型。几何建模法计算简单，但没有考虑物体材料的本构关系，因此无法保证模型的精确性(脑外科虚拟手术中软组织形变及撕裂模型研究，南昌大学，2016)。目前，常见的物理建模方法包括有限元法(Finite Element Method，FEM)和质点弹簧法(Mass-spring Method，MSM)(Soft Tissue Deformation and Optimized Data Structures forMass Spring Methods，IEEE International Conference on Bioinformatics&Bioengineering.IEEE，2009：22-24)。
FEM是一种将连续域问题离散化为若干有限元的基于连续介质理论的建模方法，模型的形变是根据生物力学作为基础求解而得，基于FEM的软组织模型具有较高的精确性(An Introduction to The Finite Element Method，Topics in EngineeringMathematics，1992(81)：37-60)。因此，有限元法在软组织建模中得到了广泛的应用。虽然FEM能精确地模拟力学特性，然而人体软组织的生物特性通常较为复杂，内部结构精细，这将导致FEM的求解消耗大量资源，无法满足虚拟手术在实时交互上的需求(A SurfaceMass-Spring Model With New Flexion Springs and Collision Detection AlgorithmsBased on Volume Structure for Real-Time Soft-Tissue Deformation Interaction，IEEE Access，2018，6：75572-75597)。
与连续物理模型不同，MSM将软组织模型离散为多个质点，质点间的相互作用由连接着质点的弹簧实现(Calibration of Mass Spring Models for Organ Simulations，IEEE/RSJ International Conference on Intelligent Robots&Systems.IEEE，2007：370-375)。与FEM相比，MSM具有计算效率高，非常适用于拓扑结构频繁变化的模型。在计算机图形领域，基于MSM的软组织建模已经得到了广泛应用。目前，已经有研究者将连续碰撞检测算法用于虚拟手术中。
然而，传统的MSM依然存在缺陷(一种用于软组织变形仿真的支撑球弹簧模型，计算机应用与软件，2013(01)：109-112)：1、模型的精度与弹簧系数密切相关，而系数大小通常仅凭经验决定；2、传统MSM的力学原理是基于胡克定律的，因此无法精确地模拟生物软组织的生物特性，如非线性、体积保持性等。3、微分方程的计算结果带来不可避免的误差，影响了模型的精确性和稳定性。研究者们为了改善传统的MSM做了许多贡献。
Qin等人提出了一种多层MSM模型框架，通过引入双弹簧系统来模拟人体皮肤和骨骼肌的生物力学特性(A Novel Modeling Framework for Multilayered Soft TissueDeformation in Virtual Orthopedic Surgery，Journal ofMedical Systems，2010，34(3)：261-271)。Duan等人提出一种适用于肝脏胆囊的MSM模型，他们采用隐式积分并对弹簧的长度施加约束以保持模型的体积不变性(Synchronous Simulation for Deformationof Liver and Gallbladder with Stretch and Compression Compensation，Conferenceproceedings：.Annual International Conference of the IEEE Engineering inMedicine and Biology Society.IEEE Engineering in Medicine and BiologySociety.Conference，2013，2013(2013)：4941-4944)。Li等人提出了弯曲弹簧提升了传统MSM的大形变下的精确性(A Surface Mass-Spring Model With New Flexion Springsand Collision Detection Algorithms Based on Volume Structure for Real-TimeSoft-Tissue Deformation Interaction，IEEE Access，2018，6：75572-75597)。然而，他们的方案均没有完整地解决上述传统模型的缺陷。
综上所述，由于人体软组织结构的复杂性，软组织模型的建模方法一直是虚拟手术系统中关键且困难的一部分，且当前鲜有能满足虚拟手术现实需求的软组织建模方法。
发明内容
针对现有技术中的不足与难题，本发明旨在提供一种基于位置约束和非线性弹簧的软组织建模方法，本发明提适用于各种形状和尺寸的组织器官，能使模型较好的满足虚拟手术的现实需求。
本发明通过以下技术方案予以实现：
一种基于位置约束和非线性弹簧的软组织建模方法，该方法包括以下步骤：
步骤1、利用三维建模技术为采集的医学图像信息构建软组织模型，所述模型由非线性弹簧和虚拟体弹簧形成的四面体单元构成，四面体模型结构用于获取表面节点信息以及模型内部体积信息；
步骤2、模拟软组织模型的形变行为，所述模型在形变过程中，非线性弹簧使模型的应力与应变呈现线性和非线性关系，虚拟体弹簧能抑制模型体积的变化，它们将共同作用来模拟软组织的生物力学特性；
步骤3、通过动力学方程预估每个节点的近似位置，然后判断节点约束方程的状态是否改变，节点约束包括弹簧长度约束和弹簧弯曲度约束，在节点的位置引入了上述两种新的约束，以增加模型的稳定性，新的约束不仅补偿了节点之间的极限距离，而且根据弹簧的弯曲角度修正了节点的相对位置；
若方程的初始状态改变，将修正节点到合理范围内的正确位置；若方程的初始状态未改变，那么非线性弹簧和虚拟体弹簧将模拟软组织的形变效果；
步骤4、计算形变的单步迭代结束，进入下一轮循环。
进一步地，在步骤3的形变计算过程中，根据牛顿运动学，计算节点受到的合力并初步预估节点位置，再依据初步计算的节点位置判断其是否满足弹簧长度和弯曲角度的约束方程，若满足则非线性弹簧和虚拟体弹簧将模拟软组织的形变效果；否则节点的位置和速度将被修正以满足约束方程，这使得模型具有较好的可控性。
步骤3的具体的步骤为：
3.1对于模型中任意节点i，根据力学微分方程预估p的近似位置：
这里，Mi代表节点i的质量，表示节点i的加速度，F(ui)表示节点i所受内力和，Fext_i表示它的外力。F(ui)由弹簧力和弹簧阻尼力构成：
上式表示表明节点i受到的内力等于所有与其相邻节点的弹簧力与阻尼力、体积力之和。其中，kij表示节点i和节点j之间的弹性系数，cij表示节点i和节点j之间的弹簧阻尼系数，表示对位移的一阶导数即节点i的速度。lij＝ui-uj表示节点i和节点j之间弹簧的距离矢量，表示节点i和节点j之间弹簧的初始距离。
为了实现生物软组织的非线性形变行为，非线性弹簧力及其相应系数的表达式应该如下：
式中，kij为连接节点i和节点j的弹簧弹性系数；k1和k2是常数。Δlij表示弹簧长度的变化。可见，在位移较小时，弹簧力和长度变化的关系为三次多项式，而较大位移时为线性关系。
体积力的公式如下：
这里Vj表示包含节点i的四面体单元的体积，k3为虚拟体弹簧系数。和分别表示四面体单元的当前体积和初始体积。
弹簧的阻尼力可以被用来模拟软组织的粘弹性。为了实现这种粘弹性特征，阻尼力的公式如下：
为实现软组织的体积保持性，不仅要考虑软组织模型表面节点间的受力，还要考虑内部结构对表面节点的影响。
以上，可以导出节点i的微分方程组：
这样，可得节点i的位置为：
ut+Δt＝ut+vi·Δt
3.2判断节点i的位置是否满足约束方程，若节点位置满足约束方程，则保留3.1中预估的节点位置；若不满足则依据约束方程修正节点位置。具体地，本发明提出了弹簧的长度约束和弯曲度约束。
弹簧长度约束适用于相邻节点之间的距离。若节点间弹簧长度形变率大于预先设定的形变率，则直接对节点位置进行修正以满足临界形变率。弹簧长度的形变率如下：
τ＝|lcur-linit|/linit
这里，τ表示弹簧的形变率linit表示弹簧的初始长度，lcur表示弹簧的当前长度。过形变率可以判断弹簧压缩或拉伸的程度，进而使用约束方程为节点产生位置修正。对于每个节点，定义弹簧长度约束方程为：
若弹簧长度超过阈值而使约束方程的状态发生变化，则弹簧两端的节点pi与pj做出相应修正：
与弹簧长度约束类似，定义δ为弹簧弯曲度，具体如下：
其中，p′i和p′j为弹簧两端节点的初始位置，pi和pj分别为两端节点的当前位置。若δ0表示弹簧弯曲度的临界值，则弹簧弯曲约束方程为：
若弹簧形变的弯曲度大于临界值，则弹簧两端的节点pi与pj位置被重新修正以满足约束方程：
与现有技术相比，本发明采用新的非线性弹簧力系统和节点位置约束的机制，使模型能有效地反馈软组织的生物力学特性此外，在计算节点位移时，不同于直接采用力学微分方程的求解结果，而是预先判断节点位置是否满足提出的约束方程——包括新的弹簧长度约束及弹簧弯曲约束。新的弹簧约束能限制弹簧形变强度的临界值，如极限拉伸长度或弯曲角度。经过位置约束后的系统是无条件稳定的，弥补了传统MSM模型不稳定的缺点。
附图说明
图1为本发明MSM模型形变计算的流程图。
图2为本发明构成软组织模型的基础四面体单元结构示意图。
图3为本发明四面体单元中的虚拟体弹簧示意图。
图4为本发明弹簧的长度约束示意图。
图5为本发明弹簧的弯曲度约束示意图。
具体实施方式
下面结合附图，对本发明作进一步地说明。
基于位置约束和非线性弹簧的软组织建模方法，如图1流程图所示，该方法步骤包括：
S1.利用三维建模技术为采集的医学图像信息构建软组织模型，该模型如图2和如图3所示由非线性弹簧和虚拟体弹簧形成的四面体单元构成；
S2.模拟软组织模型的形变行为，模型在形变过程中，非线性弹簧使模型的应力与应变呈现线性和非线性关系，虚拟体弹簧能抑制模型体积的变化，它们将共同作用来模拟软组织的生物力学特性；
S3.通过动力学方程预估每个节点的近似位置，然后判断节点约束方程的状态是否改变，节点约束包括如图4所示的弹簧长度约束和如图5所示的弹簧弯曲度约束；若方程的初始状态改变，将修正节点到合理范围内的正确位置；若方程的初始状态未改变，那么非线性弹簧和虚拟体弹簧将模拟软组织的形变效果；
3.1对于模型中任意节点i，根据力学微分方程预估p的近似位置：
这里，Mi代表节点i的质量，表示节点i的加速度，F(ui)表示节点i所受内力和，Fext_i表示它的外力；
F(ui)由弹簧力和弹簧阻尼力构成：
上式表示表明节点i受到的内力等于所有与其相邻节点的弹簧力与阻尼力、体积力之和。其中，kij表示节点i和节点j之间的弹性系数，cij表示节点i和节点j之间的弹簧阻尼系数，表示对位移的一阶导数即节点i的速度。lij＝ui-uj表示节点i和节点j之间弹簧的距离矢量，表示节点i和节点j之间弹簧的初始距离。
为了实现生物软组织的非线性形变行为，非线性弹簧力及其相应系数的表达式应该如下：
式中，kij为连接节点i和节点j的弹簧弹性系数；k1和k2是常数。Δlij表示弹簧长度的变化；可见，在位移较小时，弹簧力和长度变化的关系为三次多项式，而较大位移时为线性关系；
体积力的公式如下：
这里Vj表示包含节点i的四面体单元的体积，k3为虚拟体弹簧系数。和分别表示四面体单元的当前体积和初始体积；
弹簧的阻尼力可以被用来模拟软组织的粘弹性，为了实现这种粘弹性特征，阻尼力的公式如下：
为实现软组织的体积保持性，不仅要考虑软组织模型表面节点间的受力，还要考虑内部结构对表面节点的影响。
以上，可以导出节点i的微分方程组：
这样，可得节点i的位置为：ut+Δt＝ut+vi·Δt
3.2判断节点i的位置是否满足约束方程，若节点位置满足约束方程，则保留3.1中预估的节点位置；若不满足则依据约束方程修正节点位置，。
弹簧长度约束适用于相邻节点之间的距离。若节点间弹簧长度形变率大于预先设定的形变率，则直接对节点位置进行修正以满足临界形变率。弹簧长度的形变率如下：τ＝|lcur-linit|/linit
这里，τ表示弹簧的形变率，linit表示弹簧的初始长度，lcur表示弹簧的当前长度。过形变率可以判断弹簧压缩或拉伸的程度，进而使用约束方程为节点产生位置修正。对于每个节点，定义弹簧长度约束方程为：
若弹簧长度超过阈值而使约束方程的状态发生变化，则弹簧两端的节点pi与pj做出相应修正：
与弹簧长度约束类似，定义δ为弹簧弯曲度，具体如下：
其中，p′i和p′j为弹簧两端节点的初始位置，pi和pj分别为两端节点的当前位置。若δ0表示弹簧弯曲度的临界值，则弹簧弯曲约束方程为：
若弹簧形变的弯曲度大于临界值，则弹簧两端的节点pi与pj位置被重新修正以满足约束方程：
S4.计算形变的单步迭代结束，进入下一轮循环。
本实施例提出的MSM软组织模型由2032块四面体单元、722块三角面、515个节点构成。模型中节点的质量均为1g，非线性弹簧的参数k1为4，k2为2，Δlc为1.5mm，A为3.5，B为4。阻尼力系数b0为0.5，b1为0.1，虚拟体弹簧的参数k3为1。弹簧长度约束的临界形变率τ0为20％，弹簧弯曲约束的弯曲临界值δ0为30°；对于相邻节点p1、p2，若p1所在的四面体单元体积变化为5mm3，其四面体中心的坐标为(1，0，0)；p1的速度为4mm/s，方向为初始坐标为(2，0，0)，当前坐标为(3，0，0)；p2速度为0mm/s，初始坐标和当前坐标均为(0，0，0)，时间步长为0.1s，计算p1在t时刻的位置。在实际应用中，模型内所有节点的计算过程均与本例中相同。
第一步，根据力学微分方程预估p1的近似位置：
(1)弹簧的形变量为1mm，小于Δlc，则节点p1的非线性弹簧力为4*1+2*13＝6N，方向为
(2)p1所受阻尼力为(0.5+4*0.1)*1＝0.9N，方向为
(3)p1所受虚拟体弹簧力为1*5＝5N，方向为
(4)p1的速度变化为：
(5)p1的新坐标变化为：
第二步，判断节点p1的位置是否满足约束方程：
(1)判断节点p1的位置是否满足弹簧长度约束。弹簧的长度变化量为21％，不满足弹簧的长度约束，则节点p1的修正Δp1为0.01mm，方向为节点p2的修正Δp2为0.01mm，方向为
(2)判断节点p1的位置是否满足弹簧弯曲约束。弹簧的弯曲度为0°，未达到临界值，因此满足弯曲约束方程。
最终，计算出节点p1及节点p2的最终位置为(2.20，0，0)，(0.01，0，0)。
以上所述仅表达了本发明的优选实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形、改进及替代，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。
- 上一篇：石墨接头机器人自动装卡簧、装栓机
- 下一篇：三维重建方法、装置、电子设备及存储介质
