CN104318056B - 基于位置动力学的软组织变形和切割模拟方法 - Google Patents

Abstract

本发明提供一种基于位置动力学的软组织变形和切割模拟方法，包括了四个步骤：网格预处理阶段，根据原始三角形网格和四面体网格，计算每个顶点的纹理坐标和计算网格中边的集合，以及设置动力学模拟的参数，生成改进的四面体网格；软组织对象的变形计算阶段，根据网格顶点的当前位置和受力情况，按照位置动力学计算下一时刻网格顶点的位置；网格拓扑改变时的更新阶段，根据切割平面与网格交点的相交情况剖分原始四面体，更新几何模型和位置动力学模型的参数；视觉渲染和触觉渲染阶段，用于显示软组织对象和产生力反馈。本发明可真实模拟虚拟手术中切割软组织的过程，并具有较高的实时性。

Description

基于位置动力学的软组织变形和切割模拟方法

技术领域

本发明涉及一种基于位置动力学的软组织变形和切割模拟方法。

背景技术

随着计算机硬件处理性能的提升，基于虚拟现实技术的手术模拟器得到了广泛研究。而变形和切割模拟是虚拟手术中的关键技术。模拟软组织变形和切割主要包含三方面内容：第一是使用合适的物理模型进行变形计算；第二是处理切割时物理模型与几何模型的耦合问题；第三是系统的碰撞检测、力反馈和控制等机制。在设计变形和切割算法时要同时考虑以上三个问题，保证系统具有较好的实时性和真实性。

切割模拟在虚拟手术中非常关键，尽管前人进行了深入的研究，但切割模拟到目前为止仍然是一个很有挑战性的问题。从物理计算模型角度可以将切割模拟划分为三类：第一类，基于有限元的方法，这种方法满足软组织的生物力学属性，易于并行处理，但计算量大，很难满足实时性要求；第二类，基于几何的方法，使用几何原理近似模拟物理属性，这种方法计算速度快，但变形效果不够真实，并有可能出现变形不稳定的情况；第三类，无网格方法，这种方法适用于大形变的情况，但由于采样点比较密集，计算负担比较重。

本发明提出了一种基于位置动力学模型(position based dynamics)的软组织变形与切割模拟的方法。该方法是一种基于几何的方法，和其他几何方法不同，该方法具有较高的稳定性，并能满足实时性要求。

发明内容

本发明解决的技术问题是：克服了几何方法变形缺乏真实性和不稳定的问题，提供了一种基于位置动力学的软组织变形和切割模拟方法。该方法满足了虚拟手术系统实时性和真实性的要求。

本发明采用的技术方案为：一种基于位置动力学的软组织变形和切割模拟方法，包括以下四个步骤：

步骤(1)、网格预处理：根据输入的三角形网格和四面体网格计算四面体网格的凸包，根据最近点原则设置四面体网格顶点的纹理坐标，遍历四面体网格，计算网格中边的集合，使其中的每条边都是唯一的，得到包含纹理信息的四面体网格后，对四面体网格设置位置动力学参数。

步骤(2)、软组织对象的变形计算：根据网格顶点的当前位置和受力情况，按照位置动力学计算下一时刻网格顶点的位置，位置动力学模型中包含拉力约束和体积保持约束，采用预测-修正的方法计算位置改变量。

步骤(3)、网格拓扑改变时的更新：根据相邻时间间隔解剖刀的位置生成切割平面，后根据切割平面与网格交点的相交情况剖分原始四面体，向网格中增加新点，并修改位置动力学模型的参数设置。

步骤(4)、视觉渲染和触觉渲染：完成步骤(2)、步骤(3)后，显示软组织对象，并根据新的网格重新生成力反馈对象，产生力反馈。

本发明的原理在于：

(1)原始的位置动力学不支持体模型拓扑改变的情况。为了适应位置动力学和支持切割操作，本发明提出了一种支持拓扑改变的四面体数据组织方式。

(2)为了提高软组织变形模拟的真实性和稳定性，本发明的物理计算模型采用了较为新颖的位置动力学模型，该模型能提供较真实和稳定的变形效果。

(3)为了模拟切割过程中软组织对象因拓扑改变而发生变形的现象，本发明提出了一种基于插值思想的位置动力学参数修改的方法，用于支持位置动力学在拓扑改变时的情况。

(4)为了增加系统真实性，本发明使用了力反馈设备操作手术器械。

本发明与现有技术相比的优点在于：

1、本发明采用的物理计算模型较其他模型具有实时性高、变形较真实的特点，并易于对软组织对象进行拖拽等操作。

2、本发明提出的四面体数据组织方式，扩展了位置动力学的适用范围，并使其支持拓扑改变的情况。

附图说明

图1：基于位置动力学的软组织变形和切割模拟方法的处理流程图；

图2：用于模拟的四面体数据结构；

图3：两个点的拉力约束示意图；

图4：四面体的体积保持约束示意图；

图5：切割算法流程图；

图6：切割平面的示意图；

图7：切割平面与四面体相交的示意图；

图8：普通的分裂方法；

图9：特殊的分裂方法；

图10：第三种特殊分裂的示意图；

图11：重心坐标示意图；

图12：系统的设备和环境；

图13：系统架构；

图14-(a)：拖拽脾脏实验图1；

图14-(b)：拖拽脾脏实验图2；

图15-(a)：开始切割脾脏的实验图；

图15-(b)：结束切割的实验图；

图15-(c)：拖拽切割后脾脏的实验图1；

图15-(d)：拖拽切割后脾脏的实验图2；

图16：切割肝脏的实验图；

图17：切断肝脏的实验图；

图18：切割手臂的实验图；

图19：拖拽被切割马模型的实验图。

具体实施方式

图1给出了基于位置动力学的软组织变形和切割模拟方法的处理流程，下面结合其他附图及具体实施方式进一步说明本发明。

本发明提供一种基于位置动力学的软组织变形和切割模拟方法，主要步骤介绍如下：

1、网格预处理

该方法根据输入的三角形网格、四面体网格和动力学参数生成用于模拟变形和切割的模拟网格。最初的输入为包含纹理信息的三角形网格，按照obj格式存储，对三角形网格进行Polygon Triangulation得到四面体网格，该方法使得四面体网格上的顶点包含三角形网格的顶点。生成的四面体网格包括顶点坐标和每个四面体的顶点索引。

遍历三角形网格的每个顶点，对每个顶点寻找最近的四面体网格的顶点，将这个四面体顶点的纹理坐标设置为该三角形顶点的纹理坐标。经过这次遍历后，四面体网格包含了软组织对象的纹理信息。

遍历四面体网格的每个单元，向网格数据结构中添加边的集合与面的集合，其中边和面的集合中的对象都是唯一的。在生成面的集合时判断该面被共享的次数，如果共享次数等于1，则标记该面为外表面，若共享次数小于1，则标记该面为内部面。该过程结束后，网格数据结构包含了边集和面集。

遍历网格结构边集的每条边，计算每条边两个顶点的距离，设置拉力约束的拉力系数和静止距离，产生该边的拉力约束。遍历结束后，再遍历网格结构的每个四面体，计算四面体的体积，设置体积保持约束的约束系数和静止体积，产生该四面体单元的体保持约束。

经过上述过程后，网格预处理结束，得到用于模拟变形和切割的四面体网格数据结构。数据结构的主要内容如图2所示。

为了提高数据在内存中的连续性，TetModel结构存储全部顶点、边、面、四面体，以及约束条件信息。TetVertex结构存储网格顶点的位置、速度、质量和纹理坐标，速度和质量用于变形计算。TetEdge结构存储网格的边信息，并且网格中的每条边是唯一的。TetFace结构存储网格的面信息，三个布尔变量标识该面是否已分裂、是否为表面、是否为切口。TetElement结构存储四面体单元的信息，包括顶点索引、边索引、面索引、删除标记等。StretchConstraint结构存储拉力约束，包括顶点索引、拉力系数和静止距离。VolumeConstraint结构存储体积保持约束，包括顶点索引、体保持系数和静止体积。

2、软组织对象的变形计算

软组织的变形计算采用位置动力学模型作为物理计算模型。该模型包括如下五步：

第一步，初始化四面体每个顶点的位置P_i和速度V_i。

第二步，在每个时间间隔Δt中，按照公式(1)更新顶点的速度。

V_i'＝V_i+FΔtw_i+Damp(V_i) (1)

其中，V_i'是更新后的速度，F是外力，w_i＝1/m_i，m_i是该顶点的质量，Damp(*)是阻尼函数。

然后预测顶点位置，如下所示：

P_i'＝P_i+V_i'Δt (2)

第三步，以预测位置P_i'为参数，遍历所有约束条件，并通过迭代方法求解动态系统，得到更新后的位置信息P_i″。

第四步，按照公式(3)和(4)更新最终的位置和速度

第五步，开始下一时间间隔的计算，返回第二步，更新顶点的速度并预测位置。

位置改变量ΔP_i的计算公式为：

其中，C(P₁,…,P_n)是以所有顶点位置为输入的约束函数，是函数关于顶点P_i的梯度。

模拟软体变形时，我们考虑了两种约束条件：一种是拉力约束，另一种是体积保持约束。每一时刻都要根据两种约束进行变形计算。

(1)拉力约束

拉力约束的示意图如图3所示。针对每个拉力约束，受约束两点的位置改变量为公式(6)和公式(7)。

其中，d是P₁和P₂的静止距离。

(2)体积保持约束

体积保持约束的示意图如图4所示。体积保持的约束函数如下所示：

根据体积保持约束函数，推导得到四面体单元中每个点的梯度为：

每个顶点的位置改变量为：

每一时刻，先设置每个点位置和速度的初始值，计算每个点受到的外力，根据外力产生该点的速度。然后计算阻力，降低速度。得到速度后，进入位置更新迭代。在每次迭代中按照上面的方法所有拉力约束和所有体积保持约束。如果想提高仿真的精度，可增加迭代次数。

完成该时刻的所有计算后，四面体顶点的位置得到更新。

3、网格拓扑改变时的更新

(1)四面体分裂

四面体的切割流程按照图5进行。分裂四面体网格时，解剖刀被抽象为一条只有长度没有厚度的直线段。每一时刻的刀尖位置和刀柄位置均被记录下来，相邻时刻的对应端点连接形成切割平面P_tP_t+1Q_t+1Q_t，如图6所示。使用切割平面计算解剖刀和四面体的交点，当解剖刀和四面体接触时，开始计算四面体网格的边和面与切割平面的交点。

计算交点时，需要进行两种测试：一种是相邻时刻刀尖连线和面的相交测试，另一种是边和切割平面的相交测试。第一种相交测试计算出四面体的面交点，第二种相交测试计算出四面体的边交点，如图7所示，FV为面交点，EV为边交点。算法将新产生的交点保存在TetModel中保存下来，并在TetElement结构记录新交点的索引。

当切割平面与四面体有一个交点后，TetElement结构的bIsEnterd变量标记为true，当解剖到离开四面体时，bIsLeft标记为true。如果这两个变量都为true，则将该四面体加入到缓存列表中开始分裂。

四面体的边分裂状态和面分裂状态组成一张查找表，当解剖刀离开当前四面体时，按照查找表进行分裂。将四面体的每条边和每个面按照2的幂分配一个编码，可以保证该切割状态是唯一的，边的分配方式如表1所示，面的编码方式如表2所示。

表1四面体边的编码方式

边	AB	AC	AD	BC	BD	CD
							编码	32	16	8	4	2	1

表2四面体面的编码方式

边	ABC	ABD	ACD	BCD
					编码	8	4	2	1

切割平面和单个四面体相交有5种情况，按照图8的方式分裂。以第一种分裂方式为例，若切割平面与AB、AC、AD相交，则四面体的边状态编码为56，后按照56对应的分裂方式进行剖分。

如果边和切割平面有交点，则将TetEdge的bIsSplit变量标记为true，并在TetElement的eintrindices数组相应位置记录新增交点的索引，需要保存两个交点，因为两个交点属于不同的子四面体。如果面和切割平面有交点，则将TetFace的bIsSplit变量标记为true，并在TetElement的fintrindices数组相应位置记录新增面交点的索引。如果面是切割产生的断口，即该面不是四面体的外部面，则将TetFace的bIsIncision变量标记为true。

还有几种特殊情况，如果切割平面和四面体的顶点相交，四面体将按照图9所示的方式进行分裂，从而避免出现裂缝。情况(a)时，解剖刀和四面体的两个顶点和一条边相交。情况(b)时，解剖刀和四面体的一个顶点和两条边相交。情况(c)时，四面体的三条边V₁V₂，V₁V₃和V₃V₄和解剖刀相交于点EV₁，EV₂和EV₃，深色区域是形成的切口，分裂后的结果如图10所示。当分裂结束后，内部网格沿着切割平面分离。由于四面体单元的体积比较大，在一次切割模拟时不会产生退化四面体，因此可以避免小四面体的合并。

(2)拓扑改变时位置动力学的参数修改

切割发生时会生成子四面体，因此需要对原有的约束进行更新，一是删除被切割四面体的拉力约束和体积保持约束，二是产生子四面体拉力约束和体积保持约束。

子四面体的拉力约束需要获得子边的静止距离。边交点产生子边的静止距离如公式(14)所示。

式中，d_AP是顶点A到交点P的静止距离，l_AP是顶点A到交点P的当前距离，l_PB是交点P到顶点B的当前距离，d是线段AB的静止距离。

对于面交点P，计算子边AP的静止距离d_AP为例。首先计算P在三角形ABC内的重心坐标(λ₁,λ₂,λ₃)。然后计算静止时的情况，如图11所示，射线AP与边相交与点D，可得到线段AD的静止长度：

式中，λ₂和λ₃为P在三角形ABC内重心坐标(λ₁,λ₂,λ₃)的两个分量，d₁、d₂和d₃为AB、BC和AC的静止长度。线段AP的静止距离为：

d_AP＝(λ₂+λ₃)d_AD (16)

子四面体的体积保持约束需要计算子四面体的体积，我们根据四个顶点的位置直接计算获得。尽管子四面体的体积可能不是静止体积，但在变形程度不大的情况下，这种近似方法也可以模拟切割操作。

4、视觉渲染和触觉渲染

完成拓扑更新和参数更新和变形计算后，开始进行视觉渲染和触觉渲染。视觉渲染使用OpenGL。触觉渲染使用一台Geomagic Tough力反馈设备，该设备可以提供碰撞检测和解剖刀的位置、方向。当解剖刀或抓钳挤压软组织时，设备会根据虎克定律产生反弹力，如果解剖器械在软组织表面移动时，设备会产生摩擦力。解剖刀和抓钳的位置可以通过设备获得。系统的设备和环境如图12所示。

使用OpenHaptics时，通过如下几步建立触觉渲染：

第一步，初始化力反馈设备并创建上下文。

第二步，定义触觉帧，并将软组织对象定义在触觉帧内。

第三步，指定对象的形状并渲染对象。

在渲染时需要注意，该接口认为按顺时针定义的面是正向面。

上述几部分的关系如图13所示。

实验使用的设备为NVIDIA GeForce GTX 580、Intel(R)Xeon(R)CPU(2.53GHz，8核)和12GB RAM，运行在Windows 764位系统上。

进行了四组实验：脾脏、肝脏、手臂、马。实验的时间性能如表3所示。

表3不同模型的时间效率

表3中，拓扑更新的时间包括相交计算。从表中可以看出，物理变形时间与四面体的数量有关，降低四面体的数量可以减小物理变形的时间。拓扑更新的时间与切割区域的单元密度有关，因此会出现单元数量少的马模型比手臂模型更耗时的情况。实验数据表明本方法可以满足实时性要求，即每秒30帧。

图14至图19为变形和切割的效果图。

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (3)

1.一种基于位置动力学的软组织变形和切割模拟方法，其特征在于包括以下四个步骤：

步骤(1)、网格预处理：根据输入的三角形网格和四面体网格计算四面体网格的凸包，根据最近点原则设置四面体网格顶点的纹理坐标，遍历四面体网格，计算网格中边的集合，得到包含纹理信息的四面体网格后，对四面体网格设置位置动力学参数；

步骤(2)、软组织对象的变形计算：根据网格顶点的当前位置和受力情况，按照位置动力学计算下一时刻网格顶点的位置，位置动力学模型中包含拉力约束和体积保持约束，采用预测-修正的方法计算位置改变量；

步骤(3)、网格拓扑改变时的更新：根据相邻时间间隔解剖刀的位置生成切割平面，后根据切割平面与网格交点的相交情况剖分原始四面体，向网格中增加新点，并修改位置动力学模型的参数设置；

步骤(4)、视觉渲染和触觉渲染：完成步骤(2)、步骤(3)后，显示软组织对象，并根据新的网格重新生成力反馈对象，产生力反馈；

所述步骤(1)中所述的网格预处理方法，首先遍历三角形网格的每个顶点，对每个顶点寻找最近的四面体网格的顶点，将这个四面体顶点的纹理坐标设置为该三角形顶点的纹理坐标；经过这次遍历后，四面体网格包含了软组织对象的纹理信息；

遍历四面体网格的每个单元，向网格数据结构中添加边的集合与面的集合，其中边和面的集合中的对象都是唯一的；在生成面的集合时判断该面被共享的次数，如果共享次数等于1，则标记该面为外表面，若共享次数小于1，则标记该面为内部面；该判断该面被共享的次数过程结束后，网格数据结构包含了边集和面集；

遍历网格结构边集的每条边，计算每条边两个顶点的距离，设置拉力约束的拉力系数和静止距离，产生该边的拉力约束；遍历结束后，再遍历网格结构的每个四面体，计算该四面体的体积，设置体积保持约束的约束系数和静止体积，产生该四面体的体积保持约束；

为了提高数据在内存中的连续性，提出TetModel的数据结构，用来存储全部顶点、边、面、四面体，以及约束条件信息；TetVertex结构存储网格顶点的位置、速度、质量和纹理坐标，速度和质量用于变形计算；TetEdge结构存储网格的边信息，并且网格中的每条边是唯一的；TetFace结构存储网格的面信息，三个布尔变量标识该面是否已分裂、是否为表面、是否为切口；TetElement结构存储四面体单元的信息，包括顶点索引、边索引、面索引、删除标记；StretchConstraint结构存储拉力约束，包括顶点索引、拉力系数和静止距离；VolumeConstraint结构存储体积保持约束，包括顶点索引、体保持系数和静止体积。

2.根据权利要求1所述的基于位置动力学的软组织变形和切割模拟方法，其特征在于：步骤(2)中软组织的变形计算采用位置动力学模型作为物理计算模型,并考虑了两种约束条件：一种是拉力约束，另一种是体积保持约束；每一时刻都要根据两种约束进行变形计算；基于位置动力学模型的计算步骤如下：

第一步，初始化四面体每个顶点的位置P_i和速度V_i；

第二步，在每个时间间隔Δt中，按照公式(1)更新顶点的速度；

V_i'＝V_i+FΔtw_i+Damp(V_i) (1)

其中，V_i'是更新后的速度，F是外力，w_i＝1/m_i，m_i是该顶点的质量，Damp(*)是阻尼函数；

然后预测顶点位置，如下所示：

P_i'＝P_i+V_i'Δt (2)

第三步，以预测位置P_i'为参数，遍历所有约束条件，并通过迭代方法求解动态系统，得到更新后的位置信息P_i”；

第四步，按照公式(3)和(4)更新最终的位置和速度；

第五步，开始下一时间间隔的计算，返回第二步，更新顶点的速度并预测位置；

位置改变量ΔP_i的计算公式为：

其中，C(P₁,...,P_n)是以所有顶点位置为输入的约束函数，是函数关于顶点P_i的梯度；

另外，两种约束的表达式如下：

(1)拉力约束

针对每个拉力约束，受约束两点的位置改变量为公式(6)和公式(7)；

其中，d是P₁和P₂的静止距离；

(2)体积保持约束

体积保持的约束函数如下所示：

根据体积保持约束函数，推导得到四面体单元中每个点的梯度为：

每个顶点的位置改变量为：

3.根据权利要求1所述的基于位置动力学的软组织变形和切割模拟方法,其特征在于：步骤(3)中所述的四面体剖分和拓扑改变时位置动力学的参数修改；

四面体的边分裂状态和面分裂状态组成一张查找表，当解剖刀离开当前四面体时，按照查找表进行分裂；将四面体的每条边和每个面按照2的幂分配一个编码，保证该切割状态是唯一的；

拓扑改变时位置动力学参数的更新包括两种：一是删除被切割四面体的拉力约束和体积保持约束，二是产生子四面体拉力约束和体积保持约束；

子四面体的拉力约束需要获得子边的静止距离；边交点产生子边的静止距离如公式(14)所示；

式中，d_AP是顶点A到交点P的静止距离，l_AP是顶点A到交点P的当前距离，l_PB是交点P到顶点B的当前距离，d是线段AB的静止距离；

对于面交点P，以计算子边AP的静止距离d_AP为例；首先计算P在三角形ABC内的重心坐标(λ₁,λ₂,λ₃)；然后计算静止时的情况，射线AP与边相交与点D，可得到线段AD的静止长度：

式中，λ₂和λ₃为P在三角形ABC内重心坐标(λ₁,λ₂,λ₃)的两个分量，d₁、d₂和d₃为AB、BC和AC的静止长度；线段AP的静止距离为：

d_AP＝(λ₂+λ₃)d_AD (16)

子四面体的体积保持约束需要计算子四面体的体积，根据四个顶点的位置直接计算获得。