虚拟手术中表面网格模型的力反馈算法与仿真

潘家辉　朱玲利

1(华南师范大学软件学院　广东 南海 528225)2(洛阳师范学院信息技术学院　河南 洛阳 471002)



虚拟手术中表面网格模型的力反馈算法与仿真

潘家辉1朱玲利2

1(华南师范大学软件学院广东 南海 528225)2(洛阳师范学院信息技术学院河南 洛阳 471002)

摘要提出一个带力反馈的虚拟腹部外科手术系统，并重点研究及实现了基于表面网格模型的力反馈算法。为了提高系统的稳定性，使用基于线段与表面三角网格的碰撞检测方法；基于四阶Runge-Kutta方法的表面网格mass-spring模型进行变形仿真；基于力的广度优先传播来解决变形仿真的局部性。通过仿真实验验证了该原型系统可以实现与可变形的三维模型的实时力觉交互。最后，把该方法应用到虚拟腹部外科手术系统中，并实现带力反馈的手术仿真操作。

关键词虚拟手术表面网格质点—弹簧模型四阶Runge-Kutta算法广度优先的力传播模型

0引言

随着虚拟现实领域研究的进展和外科手术方式的复杂化和精细化，虚拟手术渐渐成为热门的研究课题。传统的外科手术面临如下困难：一方面，手术过程越来越复杂，需要详细的手术规划和术前预演。通过预计在实际手术中可能出现的复杂和险要情况，制定合理的个体化手术方案，提高手术成功率[2]。另一方面，医学院医生在手术室中的训练计划越来越少，需要通过计算机辅助技术或虚拟手术器械进行手术过程的仿真和模拟，来延长医生的学习曲线。此外，虚拟手术还可以将熟练医生的手术操作进行保存，用于解剖学和临床手术的教学，这也有利于医生重复训练重要的手术操作[3]。总之，虚拟手术具有非常巨大的社会价值和经济价值。本课题将建立一个用于手术训练的虚拟腹部外科手术系统。

腹部外科手术一直被认为是难度高、风险大的手术，主要因为腹部器官内部结构的复杂性和变异性[4]。虚拟手术系统在计算机中建立一个虚拟3D环境,通过虚拟手术器械, 提供精确的视觉仿真和细腻的力觉反馈，仿真模拟手术过程。人体对触觉相较视觉而言更加敏感，25Hz以上的视觉图像看起来就是连续的，而力觉信息的频率要达到250Hz以上才会平滑。也就是说，力觉绘制线程每4ms要运行1次以上，将耗费大量的CPU时间。在实现物理真实的形变模拟的同时，还要进行真实的力觉绘制，所以系统需要解决模拟的真实性和系统效率的问题。

对虚拟手术系统而言，力反馈具有特殊的价值，因为触觉手感是外科医生手术技巧的基础。目前，虚拟手术的研究和开发己经有了很大的进展。法国的INRIA研究小组一直致力于肝脏切割模拟器的研究，并在此基础上实现了对肝脏等软组织的切割和撕裂等仿真操作[5]。德国Karlsruhe研究中心成功开发了用于妇科手术的内窥镜虚拟手术系统，并能实现对软组织进行切割、烧灼和捏等手术操作的仿真。与国外研究情况相比，我国的研究才刚刚起步。上海交通大学谢叻等开发了一种多功能的虚拟手术器械[7]，可实现手术刀、手术剪和手术钳等操作，但还不具备力反馈功能。国防科技大学与301医院联合开发了带有力反馈的虚拟膝关节镜手术系统[8]和虚拟心脏介入手术系统[9]。为了提高虚拟切割的性能，青岛大学贾世宇等提出了一种先细分后分裂的新式交互切割方法[10]。针对虚拟软组织形变多为大变形的特点，哈尔滨工程大学朱玲等提出了一种改进的质点积分无网格的力反馈技术[11]。

目前使用的基于物理的变形方法主要有mass-spring方法和FEM[12]。相较FEM模型，mass-spring模型的计算复杂度低，拓扑结构修改方便，所以目前我们使用mass-spring模型来模拟变形。但mass-spring模型基于连续物体离散化的理论，精度系数和阻尼系数不易确定，在模拟的真实性方面不如FEM模型。文献[13]提出了使用模拟退火法参照FEM模型来调整mass-spring的劲度系数和阻尼系数的方法。用mass-spring模拟整体形变效果也是一个挑战，文献[14]实现了使用广义弹簧保持几何形状和体积、使用形状匹配实现整体形变、使用逆向动力学修正“超弹性”的mass-spring模型。而在实时力觉视觉交互的系统之中，即使易于求解的mass-spring模型的效率也有待提高。

1虚拟手术系统的设计

整个系统分为预处理和实时计算两大部分，如图1所示。

图1　虚拟手术系统的系统框架

目前我们系统的基本架构已经基本完成[15]。系统可以读入.obj格式的三维模型，使用线段来表示手术器械的形状，并进行线段与表面三角网格的碰撞检测；进行弹簧—振子的触觉绘制；使用广度优先的力传播模型来实现局部变形，降低了系统的计算复杂度；使用基于四阶Runge-Kutta方法的表面网格mass-spring模型进行变形仿真，从而在使用更大的仿真步长的同时还能得到更好的仿真效果，减少了单位时间的计算量；可以实现实时的视觉和力觉反馈。

2力反馈算法的设计

2.1手术器械表示

图2　手术刀的形状表示

在众多的手术器械当中，手术刀和电刀无疑是应用最为广泛的。它们是在腹部外科手术操作过程中，在组织或器官上造成切口的必备工具。本文分别采用线段来表示其几何特征，其他笔杆状器械也可采用该形式表示。如图2所示，灰色区域代表外形；点P是手术器械的质点，对应着笔杆状器械的尖点；点Q是器械的手柄代表点，对应着笔杆状器械的末端；方向n是器械表面的法线方向；点R是手术器械的受力代理点。本文通过这样的方式把笔杆状的手术器械表示为线段PQ，进行有关碰撞检测、力传播模型、触觉绘制方面的交互计算。

2.2碰撞检测

图3　线段与三角形的相交测试

由于本文用一条线段来表示笔杆状的手术器械，用三角网格表示腹部外科组织表面，所以手术器械和脏器组织之间的碰撞检测问题就可以简化为线段与三角网格的相交检测问题。如图3所示，线段的两个端点为P(Px,Py,Pz)和Q(Qx,Qy,Qz)，三角形的三个顶点分别为P0(x0,y0,z0)、P1(x1,y1,z1)和P2(x2,y2,z2)。相交检测的任务则是判断线段有否与模型的三角片相交，若相交，则求出交点R(Rx,Ry,Rz)。关于线段和三角形的相交检测方法分两步：

第一步是判断线段与包含三角形的平面是否相交，若相交，则求交点。

点R在线段PQ上，同时点R在三角形P0P1P2上，则有：

R=(1-t)P+tQR·n=d

(1)

其中t是比例系数，n为三角形的法向量，d为常量。可以求出：

(2)

若比例系数t满足0≤t≤1，则线段与该平面相交，代入式(1)，求出交点R。

第二步是判断点R是否在三角形P0P1P2内。

图4　三角形投影到某一平面

为了减少计算量，本文通过三角形投影，把3D问题转化到2D中。通过降维操作，将三角形投影到某一个基本水平面上。为了避免垂直或者接近垂直的情况，一般挑选投影面积最大的投影面。这可以通过检查三角形平面的法向量做到，去掉其绝对值最大的分量对应的坐标。如图4所示，这里假设投影面为XOY平面，去掉z坐标。

此外，本文通过计算交点的重心坐标，来判断它是否在三角形中。其中，重心坐标(α，β，γ)即是三角形所在平面的任意点都能表示为顶点的加权平均值。若满足0<α<1，0<β<1，0<γ<1，则点R在三角形P0P1P2中。其具体计算公式如下：

Rx=αP0x+βP1x+γP2xRy=αP0y+βP1y+γP2yα+β+γ=1

(3)

解方程组可得 ：

(4)

2.3触觉绘制

在检测出手术器械与物理模型碰撞之后，系统将构造反馈力并在力反馈设备Phantom上进行触觉绘制。对于这类笔杆状的手术器械，本文假定切割中的反馈力作用于尖端P上(如图2所示)。软组织对器械的约束力fh通过定义手术器械与接触点的弹性约束kconstraint来实现触觉绘制。kconstraint描述的是一个虚拟弹簧的约束力刚度，如图5所示。虚拟弹簧的两个端点分别是手术器械上的接触点和物体上的接触点。刚开始两个接触点重合，虚拟弹簧初始长度drelease= 0，当手术器械进入物体，手术器械上的接触点移动，弹簧的长度变为d，那么反馈力是：

fh=-kconstraint(d-drelease)

(5)

图5　力反馈模型示意图

2.4广度优先的力传播模型

很多手术操作并不会引起组织的全局变形，例如手术刀在胆囊上切割一道入口、止血钳钳住血管等。我们假设以线段表示的手术器械和物体只有一个接触点R，当这个接触点在外力的作用下发生运动时，此接触点就成为受激发的节点，在该中心受力节点上产生的应力通过与其相连接的弹簧作用在其他相邻质点上，从而把力向周围传递，带动相邻的节点运动。这样，物体的变形就由于节点的运动而产生了。广度优先的力传播模型如图6所示。

图6　广度优先的力传播模型

图6中，接触点R是广度优先搜索的第一层变形顶点，它的相邻点形成第二层变形顶点，第二层变形顶点的没有被访问过的相邻点形成第三层变形顶点。

本文使用广度优先力传播模型来实现变形的局部性[16,17]，其基本思想是外力先作用于该接触点，然后以广度优先遍历的顺序，先传播到接触点的最近邻点，然后是最近邻点的未被访问的相邻点，这样传播直到最大的遍历深度为止。如图6所示，在最大深度为2的情况下，按照广度优先的方法，接触点的反馈力在第一层传播中更新与该节点相连接的6个邻接节点的位置信息，在第二层传播中则依次更新这6个节点的邻接节点(12个邻接节点)的位置信息。具体步骤描述如下：(1) 确定遍历深度n，一般遍历深度取值越大，形变仿真的逼真程度则越高，这需要在计算复杂度和虚拟仿真度之间做出平衡，本文设定n=3。(2) 建立一个空集L0，将接触点集V0中的初始接触点加入L0，其中，在初始状态下，V0只有一个接触点R(手术器械的受力代理点)。(3) 依次创建Li={Li| i = 1，2，…，n}，并将接触点集V={Vi| i = 1，2，…，n}中所有的点加入L0，对每一个在Li中的顶点，访问其相邻点，如果该相邻点已经被访问过，不做任何工作，否则置该相邻点的状态为被访问并将该相邻点加入Li+1。(4) 当达到最大的遍历深度n时，输出变形顶点集合{L0，L1，…，Ln}。

对于手术工具和物体有多个接触点的情况，则独立地以各个接触点为树根做广度优先的遍历。若遍历到的顶点没有被访问过(可能是来自任何接触点的遍历的访问)，则将顶点状态置为已访问并将之加入变形顶点集合，否则不做任何工作。

2.5基于四阶Runge-Kutta方法的表面网格mass-spring模型

本文采用四阶Runge-Kutta方法来求解mass-spring网络，以容忍较大的仿真步长，从而减轻系统单位时间内的计算负载。Euler和Runge-Kutta方法都是求解离散差分问题的经典方法。其中，Euler方法计算简单，但精确性和稳定性不够强，而且要求较小的时间步长；而四阶Runge-Kutta方法在提高精确度的同时增加了计算量，但仍不失为一个平衡了精确度和计算量之后较好的选择[18]。

在系统中，我们定义tn表示第n+1次仿真迭代开始时的系统时间，那么自变量：

xn=tn

(6)

定义yn表示mass-spring系统的状态集合：

(7)

其中Xn是质点空间坐标的集合，Vn是质点速度的集合。

那么，对于每一质点，作用力是：

(8)

(9)

式中，ks表示弹簧的劲度系数。Δdij表示弹簧的弹性形变：

Δdij=d-drelease

(10)

式中，d表示弹簧当前长度，drelease表示零张量长度。

由式(1)，得到质点的加速度，即速度的差分：

ΔVn=f/m

(11)

式中，m是质点的质量。

而位置的差分：

ΔXn=V

(12)

式中，V是质点的当前速度。

故f函数计算差分的表达式是：

(13)

使用四阶Runge-Kutta方法求解，系统在h=40ms的步长下，得到了较为真实的模拟效果。

3实验与结果

在虚拟手术框架(如图1所示)下，本文利用力反馈设备PHANToM和图形工作站(CPU:Xeon2.80GHz，内存:4GB，显卡：NVIDIAQuadroFX5500)的平台，使用VC++和OpenHaptic开发包开发了带力反馈的虚拟手术系统MIPS，实现了实时力觉和视觉交互的物理变形。为了验证本文方法的有效性，本文进行了力反馈的物理变形仿真、带力反馈的虚拟手术两方面的实验。

3.1力反馈的物理变形仿真

由于力反馈在精度和实时性方面的要求较高，数据计算量比较大。因此，本实验目的在于测试本文力反馈算法在不同模型规模下的性能。如图7所示，系统可以实时模拟外力作用下的三维形变。如表1所示，对不同规模三角面片数的obj文件，使用广度优先的力传播模型后，依次求解mass-spring网络的时间大致相同。此外，图像绘制和触觉绘制的帧频与模型结点数是呈线性关系的。可以推算当表面网格的结点数为2500以内的模型，其图像绘制的频率可达30FPS，力反馈频率可达300FPS。此算法的运算速度基本能满足虚拟手术的需求，具有满足交互的实时性。

图7　力的渗透深度为2时的变形效果图

模型结点数面片数力的渗透深度弹簧质点网格求解耗时(ms)图像绘制帧频(FPS)触觉绘制帧频(FPS)689135531638.17325.902703496831627.65246.7871201433931612.11125.46

3.2带力反馈的虚拟手术

为验证所开发的虚拟手术系统的可操作性和变形效果，设计了左半肝切除的实验任务，由一名外科医生进行力反馈的手术仿真操作。由于PHANToM(虚拟手术刀)操纵杆也是笔杆状，因此手持操纵杆的方式可以模仿真实手术，采用指压式、握持式或执笔式。医生首先把肝脏模型导入系统，并把组成表面网格的三角面皮可视化，如图8所示。然后，通过操纵PHANToM手柄来操作虚拟手术刀，进行带有力反馈的切割操作，效果如图9所示。

图8　基于表面网格的肝脏　　图9　手术刀切割肝脏的切口

医生操纵虚拟手术刀进行左半肝切除的一系列操作。在操作过程中，医生可明显感受到来自PHANToM的作用力。图10把操作过程中虚拟手术刀的三维力反馈信息显示出来。其中，x轴为仿真手术的时间，y轴为力反馈的大小。从图10中可以看到，在左半肝切除的200秒仿真手术过程中，虚拟手术刀产生的作用力在三个维度的分量(x方向、y方向和z方向)的变化呈现一定的规律。例如，在78秒和176秒，虚拟手术刀在x、y、z方向上均呈现出一个反方向的反馈力。此外，反馈力在x、y、z方向的大小在200秒虚拟手术过程中均保持在[-1.5，1.5]的范围内。这表明反馈力分量的变化基本符合医生在左半肝切除手术过程中的运动期望，表现出稳定而连续的力觉交互效果。对于本仿真过程，医生认为虽然只是粗略模拟了仿真手术的过程，但他始终能感受到连续变化的反馈力，并在这种牵引力的引导下操纵虚拟手术刀完成模拟操作。虚拟手术中脏器的形变、手术刀的力反馈、以及整个交互操作等效果已达到仿真训练的最简单要求。

图10　手术刀切割时的三个维度的反馈力变化

4结语

本文从虚拟手术仿真系统的设计与实现出发，围绕着表面网格的力反馈算法，对手术器械表示、碰撞检测、触觉绘制、力反馈传播模型以及mass-spring模型的求解等方面进行深入的研究。下一步的工作是构造手术钳、剪刀等手术器械的几何模型和物理模型，并实现相应的手术操作。在触觉绘制方面，构造虚平面来产生更细腻、更真实、更复杂的反馈力。

参考文献

[1]SatavaRM.MedicalVirtualReality:TheCurrentstatusofthefuture[J].InProc.4thConf.MedicineMeetVirtualReality(MMVRIV),SanDiego,CA,1996,29(1):100-106.

[2]SatavaRM.Historicalreviewofsurgicalsimulation—apersonalperspective[J].Plasticandreconstructivesurgery,2004,32(2):141-148.

[3]AiKadiAS,DonnonT,PaolucciEO,etal.Theeffectofsimulationinimprovingstudents’performanceinlaparoscopicsurgery:ameta-analysis[J].Surgicalendoscopy,2012,26(11):3215-3224.

[4] 李忠廉,崔乃强,苗彬,等.腹部外科再手术原因分析:附828例报告[J].中国普通外科杂志,2007,16(2):148-150.

[5]CourtecuisseH,AllardJ,KerfridenP,etal.Real-timesimulationofcontactandcuttingofheterogeneoussoft-tissues[J].Medicalimageanalysis,2014,18(2):394-410.

[6]ColesTR,MeglanD,JohnN.Theroleofhapticsinmedicaltrainingsimulators:asurveyofthestateoftheart[J].Haptics,IEEETransactionson,2011,4(1):51-66.

[7] 谢叻,张艳,张天宇,等.虚拟手术中的力学变形和力觉感知[J].医用生物力学,2006,21(3):241-245.

[8] 徐凯.基于力反馈的虚拟膝关节镜手术系统的研究与实现[D].长沙:国防科技大学,2005.

[9] 任巨.虚拟心脏介入手术中的血流模拟和特效场景处理[D].长沙:国防科学技术大学,2005.

[10] 贾世宇,潘振宽.先细分后分裂的新式四面体网格交互切割方法[J].系统仿真学报,2012,23(12):2704-2708.

[11] 朱玲.虚拟手术中软组织形变与切割技术研究[D].哈尔滨工程大学,2012.

[12] 刘秀玲,陈栋,董亚龙,等.真实软组织特性的肝脏组织物理建模及受力分析优化[J].计算机辅助设计与图形学学报,2013,25(7):1029-1035.

[13] 蔡及时,孙汉秋,王平安.应用于网上虚拟手术的自适应变形模型[J].软件学报,2002,13(9):1893-1898.

[14] 王彦臻,熊岳山,徐凯,等.改进的采用表面网格的弹簧振子模型[J].计算机辅助设计与图形学学报,2007,19(2):168-171.

[15] 潘家辉,鲍苏苏.一种腹部外科虚拟手术仿真系统[J].计算机应用研究,2010,27(2):566-568.

[16]ChoiKS,SunH,HengPA.Anefficientandscalabledeformablemodelforvirtualreality-basedmedicalapplications[J].Artificialintelligenceinmedicine,2004,32(1):51-69.

[17]ZhongY,ShirinzadehB,SmithJ,etal.Softtissuedeformationwithreaction-diffusionprocessforsurgerysimulation[J].JournalofVisualLanguages&Computing,2012,23(1):1-12.

[18] 刘利斌,王伟.四阶Runge-Kutta算法的优化分析[J].成都大学学报:自然科学版,2007,26(1):19-21.

SURFACE GRID MODEL-BASED FORCE FEEDBACK ALGORITHM ANDSIMULATIONINVIRTUALSURGERY

Pan Jiahui1Zhu lingli2

1(School of Software,South China Normal University,Nanhai 528225,Guangdong,China)2(School of Information Technology,Luoyang Normal University,Luoyang 471002,Henan,China)

AbstractIn this paper we introduce a virtual abdominal surgery system with force feedback, and emphatically study and realise the surface grid model-based force feedback algorithm. To improve the stability of the system, we used the collision detection which is based on line and surface triangular grid, the surface grid mass-spring model which is based on fourth-order Runge-Kutta algorithm to carry out deformation simulation, and the force-based breadth-first propagation to solve the locality of deformation simulation. Through simulation experiment we verified that the prototype system was able to achieve real-time force interaction with deformable 3D model. Finally, we applied the method to virtual abdominal surgery system, and realised the surgery simulation operation with force feedback.

KeywordsVirtual surgerySurface gridMass-spring modelFourth-order Runge-Kutta algorithmBreadth-first force propagation model

收稿日期：2015-01-31。国家高技术研究发展计划项目(2012AA 021105)；广东省中科院产学研合作研究项目(2010A090100032)。潘家辉，讲师，主研领域：模式识别与智能系统，数字医学，脑机接口。朱玲利，讲师。

中图分类号TP3

文献标识码A

DOI:10.3969/j.issn.1000-386x.2016.06.062