CN107330972A - 模拟生物力学特性的实时软组织形变方法和*** - Google Patents

Abstract

一种模拟生物力学特性的实时软组织形变方法，其包括如下步骤：S1、基于软组织器官三维可视化数据生成软组织物理模型数据；S2、基于软组织物理模型数据生成碰撞检测模型数据；S3、通过游戏引擎加载物理模型数据和碰撞检测模型数据进行碰撞检测；S4、通过约束的最优化求解方法来计算软组织变形后的状态。

Description

模拟生物力学特性的实时软组织形变方法和***

技术领域

本发明涉及计算机图形学技术领域，特别涉及一种模拟生物力学特性的实时软组织形变方法和***。

背景技术

手术模拟器能够无损耗的帮助训练者熟悉真实手术流程和操作。它通过软件上的可视化渲染和硬件机械手上的交互反馈实现了科学和可重复的训练。相较于基于生物器官的手术训练，它的学习成本低，而且训练有效性也得到保证。软组织变形的模拟涉及到手术模拟器的三个关键部分:手术环境的渲染，器械器官的交互以及力反馈。因此它是手术模拟器开发中最重要的部分。通常一个优秀的形变算法使得模拟器更接近临床表现。同时准确的形变，撕裂，流血及力反馈让训练更有效。为了达到流畅和真实的操作。软组织变形模拟算法实现时需要达到真实性和实时性的兼容。然而软组织作为粘弹性体，表现出更复杂的非线性形变行为，比较于完全弹性和完全黏性体。包括蠕变，应力松弛或者不可压缩等几何特征。因此开发一个稳定而快速的实时模拟软组织变形的非线性行为的方法很重要。实时形变算法的发展伴随着硬件计算能力的提升。在早期工作中，由于较低的计算表现，几何技术被用于模拟形变。在这些情况中，物理的准确性被牺牲给计算效率，***也不了解对象的材料特性。随着计算能力的进步，基于物理的方法发展了起来。在这个领域也有着两个不同的方向。一个方向是求解基于连续物理模型的偏微分方程。比较典型的有基于网格划分的有限元方法，该方法使用同样的策略比如有限体积法和有限差分法。同时另一类基于邻近支持域离散的无网格方法最近也得到应用，如SPH，RPIM等。基于连续力学的物理模型的优势是准确性足够到能够反映真实对象的物理属性。他们被广泛应用于结构分析和电磁辐射模拟计算。但是在手术模拟的应用中的目标是达到能够反映软组织的重要物理特征，比如粘弹性不可压缩性等。因此，在这些实时的领域，选择另一类所谓看起来像物理的方法是合适的，比如弹簧振子模型，链式模型和基于位置的动态模型。通过抽象出典型的特征到这类模型中，他们能够达到很高的计算效率。

发明内容

有鉴于此，本发明提出一种模拟生物力学特性的实时软组织形变方法和***。

一种模拟生物力学特性的实时软组织形变方法，其包括如下步骤：

S1、基于软组织器官三维可视化数据生成软组织物理模型数据；

S2、基于软组织物理模型数据生成碰撞检测模型数据；

S3、通过游戏引擎加载物理模型数据和碰撞检测模型数据进行碰撞检测；

S4、通过约束的最优化求解方法来计算软组织变形后的状态。

在本发明所述的模拟生物力学特性的实时软组织形变方法中，

所述步骤S1包括：

分析软组织的拓扑结构，引入相应的约束和质点分布实现表征器官的变形变现与拓扑结构相关的效果。

在本发明所述的模拟生物力学特性的实时软组织形变方法中，

所述步骤S2包括：

通过三角形和球形的基本单元来构成软组织物理模型；

通过改进的空间hash算法对软组织物理模型中基本单元进行AABB包围盒的计算；通过空间hash算法将AABB包围盒索引存储于一个hash数组中；遍历所有的hash数组对数组内的元素进行一一碰撞检测；AABB包围盒用于投影基本单元对象到坐标轴的三个不同方向，同时维持一个二维的数组来存储每一个基元与基元之间的碰撞检测状态来避免重复检测碰撞；

在碰撞检测结束后，进行碰撞响应的计算实现达到移动基本单元到不碰撞的状态。

在本发明所述的模拟生物力学特性的实时软组织形变方法中，

所述步骤S3包括：

通过Unity3D游戏引擎进行虚拟手术场景的渲染构建；

渲染构建包括场景的碰撞检测和物理计算，场景的渲染和机械手设备数据采集处理功能被置于Unity3D运行时环境中。

在本发明所述的模拟生物力学特性的实时软组织形变方法中，

所述步骤S4包括：

将需要模拟的对象以质点的集合，约束的集合以及弹簧的集合来进行离散，在每一个固定时间步长中，需要通过计算得到的数据为质点集合的位置；

对于生物力学相关的非线性和粘弹性通过弹簧振子模型来描述，对于拓扑结构相关的特性通过基于位置的动态模拟来描述。

本发明还提供一种模拟生物力学特性的实时软组织形变***，其包括如下单元：

软组织物理模型生成单元，用于基于软组织器官三维可视化数据生成软组织物理模型数据；

碰撞检测模型生成单元，用于基于软组织物理模型数据生成碰撞检测模型数据；

碰撞检测单元，用于通过游戏引擎加载物理模型数据和碰撞检测模型数据进行碰撞检测；

软组织变形计算单元，用于通过约束的最优化求解方法来计算软组织变形后的状态。

在本发明所述的模拟生物力学特性的实时软组织形变***中，

所述软组织物理模型生成单元包括：

分析软组织的拓扑结构，引入相应的约束和质点分布实现表征器官的变形变现与拓扑结构相关的效果。

在本发明所述的模拟生物力学特性的实时软组织形变***中，

所述碰撞检测模型生成单元包括：

通过三角形和球形的基本单元来构成软组织物理模型；

通过改进的空间hash算法对软组织物理模型中基本单元进行AABB包围盒的计算；通过空间hash算法将AABB包围盒索引存储于一个hash数组中；遍历所有的hash数组对数组内的元素进行一一碰撞检测；AABB包围盒用于投影基本单元对象到坐标轴的三个不同方向，同时维持一个二维的数组来存储每一个基元与基元之间的碰撞检测状态来避免重复检测碰撞；

在碰撞检测结束后，进行碰撞响应的计算实现达到移动基本单元到不碰撞的状态。

在本发明所述的模拟生物力学特性的实时软组织形变***中，

所述碰撞检测单元包括：

通过Unity3D游戏引擎进行虚拟手术场景的渲染构建；

渲染构建包括场景的碰撞检测和物理计算，场景的渲染和机械手设备数据采集处理功能被置于Unity3D运行时环境中。

在本发明所述的模拟生物力学特性的实时软组织形变***中，

所述软组织变形计算单元包括：

将需要模拟的对象以质点的集合，约束的集合以及弹簧的集合来进行离散，在每一个固定时间步长中，需要通过计算得到的数据为质点集合的位置；

对于生物力学相关的非线性和粘弹性通过弹簧振子模型来描述，对于拓扑结构相关的特性通过基于位置的动态模拟来描述。

实施本发明提供的模拟生物力学特性的实时软组织形变方法和***与现有技术相比具有以下有益效果：首先针对于实体模型和空腔模型采用了不同的特征约束来表现拓扑结构上的不同，而对于非线性，粘弹性等生物力学特征通过粘弹簧模型来表现。在最后质点位置求解时，通过基于高斯赛德尔迭代的位置动态模拟实现无条件稳定的求解。此外本发明的算法简单，速度快，可有效地用于虚拟手术***的软组织形变仿真中。

附图说明

图1是本发明实施例的模拟生物力学特性的实时软组织形变方法流程图；

图2是实体肝脏物理模型以及空腔胆囊物理模型构建过程示意图；

图3是构建腹腔镜胆囊切除场景图；

图4是a部分为原始状态，b部分为器械提拉肝脏的效果，c部分为器械同时夹持肝脏和胆囊的效果；

图5是点与点碰撞约束示意图；

图6是点与三角形碰撞约束示意图；

图7是变形测试图；

图8是体约束效果图；

图9是全局约束效果图；

图10是方法稳定性测试图。

具体实施方式

如图1所示，本发明实施例提供一种模拟生物力学特性的实时软组织形变方法，包括如下步骤:

S1、基于软组织器官三维可视化数据生成软组织物理模型数据；

S2、基于软组织物理模型数据生成碰撞检测模型数据；

S3、通过游戏引擎加载物理模型数据和碰撞检测模型数据进行碰撞检测；

S4、通过约束的最优化求解方法来计算软组织变形后的状态。

在本发明一实施例中，所述步骤S1的具体实现方式为:首先分析软组织的拓扑结构，因为器官的变形变现与组织结构特征具有重要联系。因此表征这些与拓扑结构相关的效果，通过引入相应的约束和质点分布能够实现。举例实现的两类典型器官，肝脏和胆囊。在组织结构特点上，肝脏内部由紧密排列的组织细胞构成，而胆囊是由表皮细胞成囊状排列构成。对于肝脏的物理模型，使用四面体的拓扑结构来描述。通过四面体的边点面体信息来生成约束。而对于胆囊，使用三角面片表面约束来描述。通过三角面片边点面来生成约束。仔细来说，在模型构建上，首先是基于中国数字人的彩色切片数据进行三维重建得到器官的表面几何模型。而对于肝脏使用Tetgen进行有效的四面体网格划分，分别得到点边面体信息。除了基本生物组织模拟结构部分。对每一个四面体施加体积不变约束，如下，来表示实体模型的体积不变效果。而对于胆囊模型，直接使用表面几何特征，点表面信息。除了基本生物组织模拟部分。对于封闭表面施加了一个全局的体积约束，如下，来表示囊状模型的体积不变效果。

在本发明一实施例中，所述步骤S2的具体实现方式为:为了进行准确而快速的碰撞检测，使用了改进的空间hash算法。物理模型的基本单元为三角形和球形来构成。然后对基元进行AABB包围盒的计算。通过空间hash算法将AABB包围盒索引存储于一个hash数组中。遍历所有的hash数组对数组内的元素进行一一碰撞检测。AABB包围盒只是被用来投影基元对象到坐标轴的三个不同方向。同时维持了一个二维的数组来存储每一个基元与基元之间的碰撞检测状态来避免重复检测碰撞。在碰撞检测结束后，需要进行碰撞响应的计算来达到移动基元到不碰撞的状态。碰撞响应的过程，借助于传统的PBD，通过三种类型的距离约束来进行碰撞的移动。

在本发明一实施例中，所述步骤S3的具体实现方式为:通过Unity3D游戏引擎进行虚拟手术场景的渲染构建。包括场景的碰撞检测和物理计算，场景的渲染和机械手设备数据采集处理等功能被至于Unity3D运行时环境中。

在本发明一实施例中，所述步骤S4的具体实现方式为:将需要模拟的对象以质点的集合，约束的集合以及弹簧的集合来进行离散，在每一个固定时间步长中，需要通过计算得到的数据为质点集合的位置。质点的位置由作用于其上的力载荷所控制.力载荷分为约束力f_constraint,弹簧弹力f_spring and外力f_out.因此，基于牛顿第二定律有.

在上面的式子中,x＝[x₁ x₂ … x_n]^T质点位置列向量,M是质点质量对角矩阵,f_constraint＝[f_c1 f_c2 … f_cm]^T约束力列向量。f_spring和f_external是弹簧的弹力和外力列向量。为了导出f_constraint,约束力势能被定义为

C(x)＝[C₁ C₂ … C_m]^T是约束函数列向量,α是与约束强度相关的顺应程度.对应于重力势能和弹性势能等保守力势能.力方向是平行于能量变化最大的方向.因此约束力被定义为势能的梯度方向

为了求解质点的位置，将等式(1)以时间步长Δt进行离散,n+1代表需要迭代到的时间步长位置。

通过引入拉格朗日乘子来分解力为标量和向量部分。通过定义来简化后面的描述

通过将等式(5)(3)代入等式(4)中然后重新分布得到

等式右边的部分可以化简为这是在约束投影前需要进行计算的预测位置。为了求解时间步长为n+1的质点位置。可以得到接下来的约束最优化求解方程组。

分析上面的方程，的目标是获取λ和x在时间步长为n+1的时候.上述方程组的变量个数是等于方程的个数，因此该方程组具有唯一解.为了求解该方程组，通过牛顿拉夫逊方法来线性化然后通过迭代的数值解法来求解.将迭代次数简化为i然后省略掉n+1.得到

的目标是在每一帧的情况下获取Δx和Δλ.当迭代完成时通过λ_i+1＝λ_i+Δλ和x_i+1＝x_i+Δx来更新点的位置.应用了两个近似的方法来简化求解的过程.首先,考虑积分迭代次数为0的时候i＝0而g(x_i，λ_i)＝0,当迭代进行的过程中g(x_i，λ_i)保持一个很小的值,因此省略掉g(x_i，λ_i).第二,省略掉第一部分并不会对结果的准确性有什么影响.因此让经过这两个步骤之后，求解线性化后的方程(9)(10)为

计算等式(11)获取将其代入到等式(12)获得

当Δλ被解出时,Δx也被求解出。

现在回过头来导出弹力f_spring,非线性弹力方程如下，对于质点i有

衰减力f_damper为

弹簧弹力为

k_ij与弹簧的伸长量直接相关.因此在伸长量较小时为伸长量的三阶相关，当伸长量较长时为线性相关。

从上面的过程可以看出，的模型是通过两部分来分别反应软组织力学特性。对于生物力学相关的非线性和粘弹性就通过弹簧振子模型来描述，对于拓扑结构相关的特性就通过基于位置的动态模拟来描述。

本部分结合附图对发明的具体实施步骤流程进行解释。

本发明提供了一种模拟生物力学特性的实时软组织形变方法，主要步骤介绍如下:

1、特征器官模型处理

首先依据需要模拟的器官组织结构，将器官分为实体模型和空腔囊状模型。对于这两类模型，首先根据形状特征构建三维表面面片数据。该数据的定义格式为fbx文件格式。基本的数据单元为，顶点集合，三角形集合和贴图材质等渲染数据。如图5，对于实体模型采集三维表面面片数据中的顶点信息和三角形信息利用TetGen四面体生成算法生成空间四面体信息。对于表面模型采集顶点信息和三角形信息判断表面模型的封闭情况，消除重复顶点和边，实现完整封闭的表面网格模型。

在预处理步骤结束后，对于实体模型遍历每一个空间四面体，生成对应的四面体体积不变约束。体积不变约束对应于顶点坐标为P1,P2,P3,P4的四面体对应的约束函数如下:

对于空腔模型直接对所有的表面三角形，生成对应的全局体积不变约束。对应的全局体积不变约束函数如下:

在特征模型约束生成完毕后，对模型的通用约束进行生成。包括特定的距离约束，过拉伸约束。为了表现出生物力学特性中的非线性和粘弹性，需要生成对应的弹簧振子模型，综合模型的控制参数如下:

表1.物理控制参数描述

依据以上的控制参数能够调试模型的模拟效果。

2、碰撞检测数据生成

在物理模型数据生成完毕后，为了实现器械和器官变形之间的交互功能，需要对模型进行碰撞检测数据的生成，如图2，图2中其中从a到c为实体肝脏物理模型构建过程，从d到f为空腔胆囊物理模型构建过程。如图3所示，首先对于器械使用空间球结构来进行模拟，对于器官组织使用空间三角形结构来进行模拟。通过遍历所有的碰撞检测基元单元(分别包括，空间三角形和空间球)直接生成AABB碰撞盒集合。AABB碰撞盒集合作为欧拉空间坐标系下的轴向六面体，通过判断空间基元在x轴y轴和z轴三个方向的最大最小值来构建。构建完成所有的AABB碰撞盒集合后。通过单元长度的大小，将碰撞检测空间离散为等大小的立方体集合。该集合为空间分布的三维数组，分别具有三个方向的参数。使用一个三维到一维的空间hash映射函数能够实现将立方体单元映射为哈希数组的一个元素。因此在该步骤中，能够获取一个一维分布的空间哈希数组。哈希数组里面存储的结构为空间立方体的索引。而在同一个索引中的碰撞基元才可能发生碰撞，因此只检测同一索引中的基元的碰撞情况。图4构建腹腔镜胆囊切除场景图。a部分为原始状态。b部分为器械提拉肝脏的效果。c部分为器械同时夹持肝脏和胆囊的效果

而对于基元的碰撞情况分为两类进行求解。一类碰撞情况是空间球和空间球之间的碰撞响应，使用点与点之间的距离约束求解，如图4，图4中a部分为原始状态。b部分为器械提拉肝脏的效果。c部分为器械同时夹持肝脏和胆囊的效果。二类碰撞情况是空间球和空间三角形之间的碰撞响应，使用点与三角形之间的距离约束求解，如图5。

3、场景渲染和物理计算

在物理模型数据生成和碰撞检测数据生成结束后，为了渲染出用户操作可见的虚拟现实场景。直接将场景数据导入到Unity3D游戏引擎中，该游戏引擎能够实现次时代的三维场景渲染。而的模拟代码直接运行在Unity3D支持的Mono运行环境时，代码以C#语言结构写成.整个方法的流程被总结如下。首先非线性和粘弹性等力学特性在弹簧受力计算步骤被表现。然后通过扩展后的弹力作为外力来估计质点位置。在这一步骤结束之后，约束投影步骤开始来修正质点的位置来满足约束函数的方程。尽管使用扩展后的约束动态来避免迭代次数和形变效果相关。但是，得益于高斯塞德尔方法，最终求得的质点位置会保持稳定。

在每帧迭代求解完成之后，顶点的信息数据被输出到Unity3D的Mesh面片结构中，用于修改场景的模型几何，达到场景中的器官变形效果。在应用本发明实现的腹腔镜胆囊切除手术模拟场景中，如图6，能够看到实际的效果，图6为以垂直于组织上表面的外力拉伸组织变形的效果。

图7至图10为生物力学特性和组织结构特性测试效果图。

图7表示以垂直于组织上表面的外力拉伸组织变形的效果。

图8表示在重力作用下下垂过程中，悬梁臂模型保持局部体积不变效果。

图9表示在平板压缩过程中，气球模型保持全局体积不变效果。

图10表示从a到c显示了本发明的方法从模型压缩后恢复原状的过程，从d到f显示了传统的弹簧振子模型从模型压缩后无法恢复原状。

本发明的技术关键点和欲保护点在于：

(1)技术关键点一：将生物力学特性和拓扑结构特性分别使用非线性粘弹簧和基于位置的动态方法来模拟，实现了快速稳定的求解。

(2)技术关键点二：对于不同组织结构的器官，采用具有拓扑结构特征的不同约束来表现，实现了组织特性的针对性表现。

可以理解的是，对于本领域的普通技术人员来说，可以根据本发明的技术构思做出其它各种相应的改变与变形，而所有这些改变与变形都应属于本发明权利要求的保护范围。

Claims (10)

1.一种模拟生物力学特性的实时软组织形变方法，其特征在于，其包括如下步骤：

S1、基于软组织器官三维可视化数据生成软组织物理模型数据；

S2、基于软组织物理模型数据生成碰撞检测模型数据；

S3、通过游戏引擎加载物理模型数据和碰撞检测模型数据进行碰撞检测；

S4、通过约束的最优化求解方法来计算软组织变形后的状态。

2.如权利要求1所述的模拟生物力学特性的实时软组织形变方法，其特征在于，

所述步骤S1包括：

分析软组织的拓扑结构，引入相应的约束和质点分布实现表征器官的变形变现与拓扑结构相关的效果。

3.如权利要求2所述的模拟生物力学特性的实时软组织形变方法，其特征在于，

所述步骤S2包括：

通过三角形和球形的基本单元来构成软组织物理模型；

通过改进的空间hash算法对软组织物理模型中基本单元进行AABB包围盒的计算；通过空间hash算法将AABB包围盒索引存储于一个hash数组中；遍历所有的hash数组对数组内的元素进行一一碰撞检测；AABB包围盒用于投影基本单元对象到坐标轴的三个不同方向，同时维持一个二维的数组来存储每一个基元与基元之间的碰撞检测状态来避免重复检测碰撞；

在碰撞检测结束后，进行碰撞响应的计算实现达到移动基本单元到不碰撞的状态。

4.如权利要求3所述的模拟生物力学特性的实时软组织形变方法，其特征在于，

所述步骤S3包括：

通过Unity3D游戏引擎进行虚拟手术场景的渲染构建；

渲染构建包括场景的碰撞检测和物理计算，场景的渲染和机械手设备数据采集处理功能被置于Unity3D运行时环境中。

5.如权利要求4所述的模拟生物力学特性的实时软组织形变方法，其特征在于，

所述步骤S4包括：

将需要模拟的对象以质点的集合，约束的集合以及弹簧的集合来进行离散，在每一个固定时间步长中，需要通过计算得到的数据为质点集合的位置；

对于生物力学相关的非线性和粘弹性通过弹簧振子模型来描述，对于拓扑结构相关的特性通过基于位置的动态模拟来描述。

6.一种模拟生物力学特性的实时软组织形变***，其特征在于，其包括如下单元：

软组织物理模型生成单元，用于基于软组织器官三维可视化数据生成软组织物理模型数据；

碰撞检测模型生成单元，用于基于软组织物理模型数据生成碰撞检测模型数据；

碰撞检测单元，用于通过游戏引擎加载物理模型数据和碰撞检测模型数据进行碰撞检测；

软组织变形计算单元，用于通过约束的最优化求解方法来计算软组织变形后的状态。

7.如权利要求6所述的模拟生物力学特性的实时软组织形变***，其特征在于，

所述软组织物理模型生成单元包括：

分析软组织的拓扑结构，引入相应的约束和质点分布实现表征器官的变形变现与拓扑结构相关的效果。

8.如权利要求7所述的模拟生物力学特性的实时软组织形变方法，其特征在于，

所述碰撞检测模型生成单元包括：

通过三角形和球形的基本单元来构成软组织物理模型；

通过改进的空间hash算法对软组织物理模型中基本单元进行AABB包围盒的计算；通过空间hash算法将AABB包围盒索引存储于一个hash数组中；遍历所有的hash数组对数组内的元素进行一一碰撞检测；AABB包围盒用于投影基本单元对象到坐标轴的三个不同方向，同时维持一个二维的数组来存储每一个基元与基元之间的碰撞检测状态来避免重复检测碰撞；

在碰撞检测结束后，进行碰撞响应的计算实现达到移动基本单元到不碰撞的状态。

9.如权利要求8所述的模拟生物力学特性的实时软组织形变***，其特征在于，

所述碰撞检测单元包括：

通过Unity3D游戏引擎进行虚拟手术场景的渲染构建；

渲染构建包括场景的碰撞检测和物理计算，场景的渲染和机械手设备数据采集处理功能被置于Unity3D运行时环境中。

10.如权利要求9所述的模拟生物力学特性的实时软组织形变***，其特征在于，

所述软组织变形计算单元包括：

将需要模拟的对象以质点的集合，约束的集合以及弹簧的集合来进行离散，在每一个固定时间步长中，需要通过计算得到的数据为质点集合的位置；

对于生物力学相关的非线性和粘弹性通过弹簧振子模型来描述，对于拓扑结构相关的特性通过基于位置的动态模拟来描述。