重力作用下肝脏模型变形及切割模拟仿真

张峰峰， 张石， 黄柯， 于凌涛， 孙立宁

(1.苏州大学 机电工程学院,江苏 苏州 215006； 2.苏州大学 苏州纳米科技协同创新中心,江苏 苏州 215123； 3.哈尔滨工程大学 机电工程学院,黑龙江 哈尔滨 150001)

部分肝脏切除术是目前治疗肝癌的根治性手段，同时也是治疗一般良性肝脏疾病的有效手段[1-2]。随着医学技术的发展，将增强现实或者虚拟现实等辅助手术系统应用于肝脏切割手术过程中，能够解决传统肝脏手术过程中存在的术前和术中信息不同步的问题[3-4]。在实际的肝脏手术过程中，由于肝脏是非刚体，当手术器械对肝脏进行按压或者切割时会产生相应的形变，对于肝脏形变模拟的真实程度在很大程度决定了辅助手术系统的准确性，为准确反映实际肝脏形变特征，建立一个合理的肝脏组织生物形变模型至关重要。根据切割肝脏过程中变形是否具有真实依据可分为物理形变模型和非物理形变模型[5]。

非物理形变模型主要根据用户输入或数学计算直接改变组织模型形态对形变进行模拟，主要包括基于控制点的模型[6]和自由形变模型[7]。物理形变模型依据现实的物理定律构建变形的物理模型[8]，主要包括边界元模型[9]、有限元模型[10-11]和质点弹簧模型[12-14]，其具有准确度高、形变真实的特点，也因计算量大导致较低的刷新率。本文采用基于质点弹簧的肝脏模型，分别在重力作用下变形以及重力作用下切割进行分析，依据肝脏的生物力学特性建立了基于质点弹簧的物理形变模型，进行部分肝脏切除术中肝脏变形仿真。

1 质点弹簧软组织模型

1.1 基于Tetgen的四面体剖分

通常由CT重建获得的模型为三角网格表面模型，使用弹簧阻尼结构将模型表面进行连接。重力作用下不足以支撑其三维结构[15]，因此需要在肝脏模型内部插入新顶点，并将其与表面顶点连接，形成四面体或六面体等结构。

由给定点集生成四面体的过程称为四面体剖分，四面体剖分的方法主要有点集网格剖分、区域网格剖分和限定网格剖分。点集网格剖分要求给定点集，生成的四面体顶点均为点集内的点；区域四面体剖分要求给定一个区域边界，在边界内部生成四面体，内部四面体的顶点没有限制，可在边界内部任意插入顶点；可对点、边或面添加控制条件。本文主要针对由CT重建出的三维表面网格模型进行四面体剖分，模型为肝脏外表面的边界数据，因此使用区域固体剖分将肝脏模型四面体化。

Tetgen为开源四面体剖分工具，其可接受3D边界网格模型作为输入，在其内部插入新顶点，生成四面体单元。本文使用QT结合Tetgen编写了一个通用的四面体剖分程序，程序可输入任意三维网格模型，输出四面体模型。程序主要包括网格质量控制部分和输出文件的选择部分。

网格质量控制部分包括一个保留表面网格的选项和一个质量控制滑块。当选择保留表面网格时，程序在进行四面体剖分时将保留原模型的表面网格，只对模型内部进行操作。否则，程序在进行四面体剖分时将同时对表面网格进行重构。质量控制滑块能够对模型内部四面体单元的质量进行控制。四面体单元的质量由2种方法进行控制：半径-边缘比约束和最小二面角约束。四面体T的半径-边缘比ρ(T)为：

(1)

式中:r为T的外接圆半径；d为四面体最短边长度；θmin是四面体T的最小二面角。因此，半径-边缘比约束等价于最小二面角约束。

一个四面体模型由多个文件组成，程序输出的文件为：1)nodes文件为模型顶点列表的储存文件，每一行对应一个顶点的坐标值；2)face文件为模型三角面片列表的储存文件，每一行对应一个三角面片的3个顶点。face文件包括模型表面的三角面片和内部面片。face文件最后一列为标识位，用来区分表面和内部面片，1为表面面片，0为内部面片；3)ele文件为四面体单元列表的储存文件，每一行对应一个四面体单元4个顶点在nodes文件中的索引；4)edge文件为四面体边列表的储存文件，每一行对应一条边的2个顶点。edge文件同样包括模型表面的边和内部的边，最后一列标识位同face文件；5)neigh文件文件为邻接四面体列表的储存文件。neigh文件每行有5列，第1列为编号，后4列为一个四面体4个顶点对面的四面体，若四面体位于表面则该四面体只有3个邻接四面体，另一个位置表示为-1作为占位符。

使用上述四面体剖分程序对肝脏表面网格模型进行四面体剖分，输出结果如图1所示。

图1 四面体剖分结果Fig.1 Tetrahedral meshing results

1.2 质点弹簧力学模型

应力松弛、蠕变和滞后统称为粘弹性特性，肝脏组织的粘弹性特性是肝脏变形过程中的重要性质。一个物体突然发生形变，之后保持该形变，应力随时间逐渐减小的过程称为应力松弛。若一个物体内突然产生应力并保持应力不变，物体持续发生变形的现象称为蠕变。若一个物体承受循环载荷，但加载时应力应变关系与卸载时的应力应变关系不同的现象称为滞后。目前用于描述物体粘弹性性质的模型主要有Maxwell模型、Voigt模型和Kelvin模型[16]，如图2所示，图中μ和η分别为弹簧的弹性系数和粘性系数;σ为应力；ε为应变。Maxwell模型由一个线性弹簧和一个阻尼器串联而成。当物体突然产生应变时，弹簧将储存弹性势能，阻尼器在弹簧的影响下逐渐发生变形。因此，Maxwell模型更加适用于描述物体的应力松弛现象；Voigt模型由一个线性弹簧和一个阻尼器并联而成。当物体突然产生应力时，由于阻尼器的存在，物体并不会立即产生应变，而是在阻尼器的影响下逐渐变形。因此，Voigt更适合描述物体的蠕变现象；Kelvin模型可看作将一个Maxwell模型与一个线性弹簧并联而成。Kelvin模型能够同时描述物体的应力松弛和蠕变现象。

图2 质点弹簧模型示意Fig.2 Schematic diagram of the mass point spring model

本文主要针对肝脏在重力作用下的变形和切割，其过程更加接近于蠕变。因此采用Voigt模型构建肝脏的物理形变模型。根据四面体剖分软件导出的.edge文件将肝脏模型中的顶点使用Voigt元件连接起来，其中每一个顶点受力方程为：

fext(i)+fvoigt(i)=fsum(i)

(2)

式中：fext(i)为质点i受到的外力；fvoigt(i)为质点i受到Voigt元件产生的力；fsum(i)为质点i受到的合力。

质点i受到的外力由手术器械的交互力fc(i)和重力fg(i)组成：

fext(i)=fc(i)+fg(i)=Δs·pi+mig

(3)

Voigt元件产生的力由阻尼力fdmp和弹簧力fspr组成：

(4)

式中：v为质点的移动速度，r为质点的空间坐标，nij为连接质点i和质点j的Voigt元件的阻尼系数，kij为连接质点i和质点j的Voigt元件的刚度系数。

质点i受到的合外力由牛顿第二定律得出：

fsum(i)=miai

(5)

在肝脏受重力变形和切割变形的过程中将肝脏组织视为各向同性的结构，即mi=m,ηij=η,kij=k,则有：

(6)

计算时对时间进行离散化处理，每隔Δt=0.02 s进行一次计算。计算质点的加速度a，则Δt后质点的速度为：

vi+1=vi+Δv=vi+a·Δt

(7)

质点在Δt内的运动可视为匀速运动，因此，Δt后质点的坐标为：

ri+1=ri+vi·Δt

(8)

在Δt时间内计算所有质点的位置坐标即可得出0.02 s后肝脏模型的整体形态，当所有顶点的受力达到平衡时，肝脏模型不再变形，此时即为肝脏模型在重力作用下的形态。

2 离体猪肝切割及变形实验

本文使用质点弹簧模型模拟肝脏在平面上受重力影响的变形和在重力影响下的切割变形，本节以猪肝为研究对象，分别对以上2种情况进行实验验证。

2.1 肝脏模型在重力下的变形实验

要对肝脏模型在重力影响下的变形进行验证，必须得知肝脏在失重状态下的原始形态并获得其失重模型。因此需要模拟一个失重的环境，将猪肝置于其中，获取猪肝在失重状态下的三维数据。

通常模拟失重的方法有直接法和间接法。直接法以某种方式使物体进入真正的失重状态，一种方式是使物体进行自由落体运动，另一种方式为抛物线飞行。间接法为物体提供一个类失重的环境，主要有浸水法和风洞法，还有卧床法和悬吊法等，但后面2种方法旨在模拟失重状态下其他影响，并不能很好地模拟失重状态下物体的形态。

在本实验中为了保持猪肝在失重状态下的形态并获得其三维数据，采用浸水法模拟失重环境。研究表明猪肝的密度与水接近略大于水，因此可将猪肝浸入水中，然后向水中加盐增加水的密度直至猪肝在盐水中处于悬浮状态。为了获取真实的切割路径，在猪肝上固定2个小球，切割时在2个小球之间沿直线切割。当猪肝在盐水中稳定后对其进行CT扫描，经三维重建获得悬浮状态下的猪肝模型，如图3所示。

图3 悬浮状态下的猪肝模型Fig.3 Pork liver model in suspension

CT扫描完成后将猪肝从盐水中取出放置在平面上再次扫描，重建得到猪肝在重力下的形变模型作为对比标准。将猪肝的悬浮模型导入Unity并赋予程序脚本，新建一个平面，将猪肝模型置于略高于平面的上方，运行程序模拟猪肝在重力下的变形。待猪肝完全落地并稳定之后导出变形后的模型，如图4所示。

图4 Unity导出并实体化的重力形变模型Fig.4 Gravity deformation model exported and solidified by Unity

2.2 肝脏模型在重力下的切割实验

将取出猪肝沿预先设定的切割路径进行切割形成切口，具体切割算法流程已在文献[17]进行了阐述，随后再次进行CT扫描，利用Mimics重建得到猪肝在重力下的切割形变模型，如图5所示。

图5 猪肝在重力下的切割形变模型Fig.5 Cutting deformation model of pork liver under gravity

沿与实际切割路径相同的路径对猪肝进行模拟切割，切割完毕导出Unity下的切割形变模型，如图6所示。

图6 Unity导出并实体化的切割模型Fig.6 Unity exports a materialized cutting model

为了更好地将真实模型与虚拟模型进行比对，将三角面片的壳体模型转换为实体模型。该过程使用Geomagic Studio三维逆向软件的精确曲面功能完成。软件自动计算贴合网格模型的曲面，期间会生成少数劣质曲面，需要手动调整。分别将猪肝在平面上的重力形变网格模型与对应虚拟模型转换为实体模型。

3 离体猪肝实验结果讨论

离体猪肝四面体剖分后顶点数量为6 510个，四面体单元为24 398个，三角单元为53 796个。将猪肝的平面重力形变实体模型和切割实体模型与对应的虚拟模型进行对比。模型对比使用Geomagic Quality的3D比较功能完成。设置真实形变模型为参考对象，虚拟模型为测试对象，分别对猪肝在平面上的重力形变模型和切割模型与相应的虚拟模型进行对比分析，并导出误差报告文件。

图7为猪肝在重力下的变形误差云图，图中以真实形变模型为参考，误差为正表示虚拟模型在真实模型外部，误差为负表示虚拟模型在真实模型内部。可以看出整体误差分布在-4～4 mm，由于底面为平面，误差较小，约为-1 mm。正面误差较大，模型前左半部分误差为正，约2.5 mm，虚拟模型略厚于真实模型。模型底面误差为负，虚拟模型的底面略高于真实模型。另外，模型边缘部分误差较大，即对肝脏边缘的形变模拟较差。表1为肝脏重力形变的误差分布表，图8为误差分布的柱状图，可以看到大部分误差分布于-2～2.5 mm，占所有误差的90%左右。

表1 重力变形误差采样点分布

图7 猪肝在重力下的变形误差Fig.7 Deformation error of pork liver under gravity

图8 重力形变误差分布Fig.8 Error distribution of gravity deformation

图9为肝脏重力切割的误差云图，整体误差与重力下的形变误差相似，误差主要分布于肝脏正面。图中切口处可以看出，真实模型与虚拟模型的切口深度基本相同，切口越接近表面的部分误差越大。切口上半部分误差介于-3～-2 mm，误差整体为负，说明虚拟切口的开口大于真实切口。切口下半部分误差介于0～3 mm，为正说明虚拟切口的开口小于真实切口。表2为肝脏切割形变的误差分布表，图10为误差分布的柱状图。由于切口部分所占比例较小，因此，切口处的误差对整体误差分布影响较小，误差分布与重力形变的误差分布基本相同。

图9 猪肝在重力下切割的形变误差Fig.9 Deformation error plot of pork liver cut under gravity

图10 切割形变误差分布Fig.10 Error distribution of cutting deformation

表2 切割形变误差采样点分布

通常计算机图像的刷新率需要达到30 FPS才能感觉到流畅[18]，为了避免出现卡顿用以保证实时性，本文利用Unity中的Compute Shader分别对肝脏模型重力变形以及重力切割程序进行了加速实验。图11和图12分别为经过加速后肝脏模型重力变形刷新率和重力切割程序刷新率，可以看出，刷新率分别达到了300 FPS和200 FPS，完全能够满足流畅性和实时性的要求。

图11 重力变形程序的刷新率Fig.11 Refresh rate of gravity deformation program

图12 重力切割程序的刷新率Fig.12 Refresh rate of gravity cutting program

4 结论

1)利用该模型获得离体猪肝分别在重力作用下的变形及切割误差分别在2 mm及3 mm以内，同时通过GPU加速，模型刷新率在300 FPS和200 FPS，能够完成增强现实手术导航对软组织模型实时性和准确性的要求。

2)在部分肝脏切除手术过程中肝脏受力状况复杂，并不仅仅受到重力作用，还需考虑与手术器械之间的相互作用力，这也是后续研究需要展开的。