基于位置动力学的软组织切割算法的研究与实现

中图分类号：TP391.9论文编号：10006SY1206205硕 士 学 位 论 文基于位置动力学的软组织切割算法的研究与实现作者姓名 白隽瑄学科专业 计算机应用技术指导教师 郝爱民 教授培养学院 计算机学院Research of video saliency detection based on inter-frame correlationanalysisResearch and Implementation of Cutting Simulation for Soft Tissue using Position Based DynamicsA Dissertation Submitted for the Degree of MasterCandidate：Bai Junxuan Supervisor：Prof. Hao AiminSchool of Computer Science & EngineeringBeihang University, Beijing, China中图分类号：TP391.9论文编号：10006SY1206205硕 士 学 位 论 文基于位置动力学的软组织切割算法的研究与实现作者姓名 白隽瑄 申请学位级别 工学硕士指导教师姓名 郝爱民 职 称 教授学科专业 计算机应用技术 研究方向 虚拟现实与可视化技术学习时间自 年 月 日 起 至 年 月 日 止论文提交日期 年 月 日 论文答辩日期 年 月 日学位授予单位 北京航空航天大学 学位授予日期 年 月 日关于学位论文的独创性声明本人郑重声明：所呈交的论文是本人在指导教师指导下独立进行研究工作所取得的成果，论文中有关资料和数据是实事求是的。尽我所知，除文中已经加以标注和致谢外，本论文不包含其他人已经发表或撰写的研究成果，也不包含本人或他人为获得北京航空航天大学或其它教育机构的学位或学历证书而使用过的材料。与我一同工作的同志对研究所做的任何贡献均已在论文中作出了明确的说明。若有不实之处，本人愿意承担相关法律责任。学位论文作者签名： 日期： 年 月 日学位论文使用授权书本人完全同意北京航空航天大学有权使用本学位论文（包括但不限于其印刷版和电子版），使用方式包括但不限于：保留学位论文，按规定向国家有关部门（机构）送交学位论文，以学术交流为目的赠送和交换学位论文，允许学位论文被查阅、借阅和复印，将学位论文的全部或部分内容编入有关数据库进行检索，采用影印、缩印或其他复制手段保存学位论文。保密学位论文在解密后的使用授权同上。学位论文作者签名： 日期： 年 月 日指导教师签名： 日期： 年 月 日摘 要在计算机数字空间中开展手术规划、预演和训练的虚拟手术，是推动现代医学不断向精准化、个性化、微创化发展的前沿技术，也是提高医学教学、医学研究、临床诊疗水平的有效技术手段。虚拟手术中的软组织变形与切割仿真，需要同时满足视觉逼真、物理逼真和实时交互三方面要求，是虚拟手术技术与系统研究中的难点技术问题。本文以腹腔镜微创手术模拟为研究背景，系统研究了基于位置动力学模型的软组织变形和切割算法。取得的主要成果和技术进展包括：1. 在分析并扩展传统位置动力学模型的基础上，给出了一种物理模型驱动几何模型的软组织变形方法。首先，分析了传统位置动力学模型的数学原理、具体算法和约束条件；然后，基于位置动力学易扩展的特性，使用虚拟弹簧定义了能量保持约束，以表现软组织的材质属性，提高了变形的效果；最后，给出了物理模型的生成方法、几何模型与物理模型的耦合方法，并使用基于四面体网格的物理模型驱动基于三角形表面网格的几何模型变形。2. 给出了夹取、按压和切割三种手术操作的模拟方法。首先，使用拖拽最近点的方法实现夹取操作；然后，通过施加冲量的方法实现按压操作；对于切割操作，本文提出了两种实现算法，一种是基于混合网格的切割算法，考虑到其计算复杂性和实现复杂性，本文又提出了一种基于简化模型的切割算法，在保证变形效果和视觉真实性的前提下，提高了算法的实时性能。3. 设计并实现了手术仿真系统中基于力反馈设备的触觉渲染算法。本文使用力反馈设备控制手术器械，对夹取、按压和切割操作增加了触觉编程，提高了仿真的代入感。上述算法已初步集成于腹腔镜虚拟手术模拟器中，实现了具有力反馈的夹取、按压和切割等虚拟手术核心功能。关键词：虚拟手术，位置动力学，物理变形，切割模拟，触觉渲染AbstractVirtual surgery simulation for surgery planning, preview and training plays an important role in developing advanced surgery techniques which are more accurate, personalized and minimally invasive. Moreover, virtual surgery simulation is an effect method for improving medical teaching, medical research and clinical treatment. There are three kinds of requirements for soft tissue deformation and cutting simulation: the visual output should be as real as possible, the deformation should be physically-accurate, and it should be able to interact with users in real-time. These three requirements are extremely challenging in virtual surgery simulation.Taking laparoscopic surgery simulation as research background, this paper systematically studies the deformation and cutting simulation for soft tissue based on position based dynamics. The main contributions of this paper are summarized as follows:1. After analyzing and extending the original position based dynamics, a soft tissue deformation method is proposed that the geometric model (shape) is driven by the physical model. Firstly, this paper analyzes the mathematical principles, algorithms and constraints functions in position based dynamics. Secondly, energy preserving constraint is introduced into position based dynamics which is computed by virtual springs, this new constraint is used to improve soft tissue deformation by considering material properties. Thirdly, this paper presents the method for creating a physical model and the method for coupling geometric model and physical model. Finally, this paper introduces how the geometric model deforms with the physical model, the geometric model employs a triangular surface mesh and the physical model employs a tetrahedral mesh.2. The methods for dragging, extrusion and cutting operation are proposed. The nearest vertex of the physical model is selected to simulate dragging operation. The extrusion operation is implemented by exerting an impulse which is relevant to the device force. For cutting simulation, two mechanisms are designed. The first one subdivides the triangular surface mesh and tetrahedral mesh simultaneously, but it is very complex and not practical. The second one subdivides a simplified model, it can greatly reduce the complexity without sacrificing visual and deforming performance.3. The haptic device is incorporated into the simulation system and the haptic rendering algorithm is implemented. The aforementioned operations exert a haptic effect through the device that can greatly enhance the simulation environment. Finally, all the mentioned algorithms are incorporated into a simple laparoscopic surgery simulator which provides dragging, extrusion and cutting operations.Key words: Virtual Surgery Simulation, Position Based Dynamics, Physical Deformation, Cutting Simulation, Haptic Rendering目 录第一章绪论11.1选题背景与来源11.2研究意义21.3研究目标与内容41.4论文组织结构6第二章关键技术及研究现状72.1引言72.2虚拟手术仿真中的物理变形方法72.2.1质点弹簧模型72.2.2有限元模型82.2.3无网格模型102.2.4锁子甲模型102.2.5位置动力学模型112.3虚拟手术中的切割方法122.3.1直接删除法132.3.2顶点挤压法142.3.3沿边界分裂法152.3.4单元细分法152.3.5单元复制法172.3.6无网格切割法182.4本章小结18第三章扩展位置动力学模型和软组织变形213.1位置动力学模型（Position-based dynamics）213.1.1模型原理213.1.2算法流程223.1.3传统位置动力学模型中的约束条件233.2能量保持约束253.3物理模型的生成方法263.4几何模型和物理模型的耦合273.5软组织变形的实验283.5.1基于位置动力学模型的软组织变形283.5.2时间效率分析303.5.3位置动力学模型与质点弹簧模型的比较303.5.4能量保持约束的影响313.6本章小结32第四章基于位置动力学模型的手术操作模拟354.1虚拟手术仿真中的夹取操作354.2虚拟手术仿真中的按压操作354.3虚拟手术仿真中的切割操作364.3.1混合网格模型的切割操作374.3.2三角形表面网格与四面体网格的更新处理474.3.3实验结果504.3.4混合网格模型切割操作的问题514.3.5简化模型的切割操作524.4本章小结55第五章触觉环境的建立与交互575.1力反馈设备的功能与参数575.2触觉环境的建立585.3手术操作的交互方法595.3.1夹取操作595.3.2按压操作605.3.3切割操作615.4本章小结62第六章简易虚拟手术模拟器的设计与实现636.1系统整体设计636.1.1开发环境636.1.2手术仿真系统的框架646.2重要的数据结构666.2.1三角形网格相关的数据结构666.2.2四面体网格相关的数据结构676.2.3手术器械类696.2.4数学计算相关的类706.3系统界面716.4本章小结72总结与展望73参考文献75攻读硕士学位期间取得的学术成果81致 谢83图 目图1 Mentice VIST-Lab（引自文献3）3图2 Simbionix LAP Mentor（引自文献4）3图3 VirtaMed HystSim（引自文献5）3图4 质点弹簧模型示意图（引自文献8）7图5 XFEM加入浓缩函数的示意图（引自文献21）10图6 使用锁子甲模型的脑部软组织变形（引自文献30）11图7 使用位置动力学模型的布料仿真（引自文献32）12图8 切割方法分类图（引自文献41）13图9 直接删除法的效果图（引自文献17）14图10 改进的直接删除法的效果图（引自文献42）14图11 顶点挤压法的示意图（引自文献43）15图12 沿边界分裂法的示意图（引自文献45）15图13 预定义的分裂方式（引自文献9）16图14 3种面分裂方式和5种体分裂方式（引自文献10）16图15 对六面体单元进行细分的二维示意图（引自文献46）17图16 虚拟节点实现切割的方法（引自文献48）17图17 无网格切割法的示意图（引自文献28）18图18 拉伸约束示意图23图19 体积保持约束示意图24图20 能量保持约束的示意图26图21 NVidia PhysX Viewer界面27图22 重心坐标的示意图28图23 肝脏模型的几何模型与物理模型29图24 拖拽肝脏模型的变形效果图29图25 脾脏和胆囊的变形效果图30图26 位置动力学模型与质点弹簧模型的变形效果对比31图27 增加能量约束后的变形效果图32图28 不同弹性系数在固定网格一端时的比较32图29 不同弹性系数在下落时的比较32图30 夹取操作的示意图35图31 发生按压时表面的受力分析36图32 切割平面的生成37图33 三角形的分裂方法39图34 切割四面体的编码示例40图35 四面体的细分方式41图36 混合网格的切割算法流程图46图37 利用LinkedSVlist更新表面顶点47图38 三角形网格交点寻找包含四面体的方法48图39 交点不在包含端点四面体中的情况48图40 混合网格的肝脏切割效果图50图41 混合网格的躯干切割效果图50图42 和文献26的切割效果比较51图43 混合网格模型错误分配顶点的示意图52图44 产生锯齿状边缘的效果图52图45 新增的四面体细分方式53图46 切割四面体网格的肝脏模型53图47 切割四面体网格的脾脏模型54图48 切割四面体网格的马模型54图49 切割四面体网格的手臂模型54图50 Geomagic Touch力反馈设备（引自文献57）57图51 HDAPI编程流程图58图52 HLAPI编程流程图59图53 代理位置和设备位置的关系61图54 手术仿真系统框架64图55 混合网格的切割算法框架65图56 简化模型的切割算法框架65图57 虚拟手术仿真系统的显示界面71图58 虚拟手术仿真系统的硬件设备72表 目表1 阻尼计算的伪代码22表2 不同模型变形的时间效率30表3 位置动力学模型与质点弹簧模型的时间效率比较31表4 三角形网格的细分算法38表5 三角形边的编码方式38表6 四面体网格的细分算法39表7 四面体边的编码方式40表8 四面体面的编码方式40表9 所有细分方式的查找表41表10 EdgeInfo结构的成员变量42表11 第1种细分方式的边生成方法42表12 第2种细分方式的边生成方法43表13 第3种细分方式的边生成方法43表14 第4种细分方式的边生成方法44表15 第5种细分方式的边生成方法44表16 表面顶点列表（LinkedSVList）的成员变量45表17 与有限元方法的时间效率比较51表18 夹取操作涉及的HDAPI函数60表19 按压操作涉及的HLAPI函数60表20 切割操作涉及的HLAPI函数61表21 MyObjEdge结构的成员变量66表22 MyObjTriangle结构的成员变量66表23 MyTetModel结构的成员变量68表24 ElementInfo结构的成员变量68表25 CScalpel类的成员变量69表26 CScalpel类的函数69表27 CMyMath类的函数70XIII北京航空航天大学硕士学位论文第一章 绪论1.1 选题背景与来源随着计算机硬件和计算机仿真技术的提高，越来越多的交叉学科运用了虚拟现实技术（Virtual Reality Technology）来替代原本昂贵、不可重用的实验，虚拟手术仿真（Virtual Surgery Simulation）是其中之一。传统训练方法使用尸体或由硅胶制成的内脏模型对专业医疗人员进行手术训练，训练设备只能够一次性使用，因此医院在训练设备上的开销巨大。虚拟手术仿真使用虚拟现实技术来模拟手术过程，制造的模拟器被医院用于训练医疗人员，由于模拟器具有可重复使用的特性，从而大大降低了医院的训练成本。现代医学的快速发展产生了多种新型手术方式，这些新产生的手术方式既能减轻病人在手术过程中的痛苦，又能降低手术操作的风险和创伤，微创手术便是其中之一。“爱护组织，减少创伤”是外科医生追求的理念，内窥镜和相关技术与设施的发展，解决了手术入路微型化但能获得有限的充分暴露，促使微创外科得到迅速发展1。典型的微创手术包括腹腔镜手术、宫腔镜手术等。以腹腔镜手术为例，医生在病人被麻醉后在其腹壁打四个直径为10 mm的小孔，在其中一个小孔中放入微型摄像机，在另外两个小孔中放入手术器械，最后一个小孔用于制造气腹腔或牵拉组织等辅助操作。微型摄像机连接到显示器上，并将病人的脏器环境显示出来，医生则看着显示器对病变的组织进行夹取、切割、缝合等操作。最后，医生将切除的器官组织放入标本袋内并取出2。这种手术方法给病人带来了福音，却给医生提出了更高的要求。在传统手术中，医生能够直接观察病灶部位，但在微创手术中，医生必须面对显示器进行手术，医生的手眼协调能力受到了挑战，因此医生必须经过充足的训练才能胜任手术。使用尸体进行训练能够提供一个真实的解剖环境，但尸体与活体有不同的生物特性，如颜色、刚度等，从而使二者对手术器械的响应不同，并且尸体的来源较少。使用人工合成物的价格低廉，但软组织的生物特性并不能很好的体现。使用动物进行活体实验一方面价格昂贵，另一方面易受到舆论上的指责。除此之外，内窥镜手术与传统手术的操作方式不同，使用传统方法对医生进行训练不能提供真实的手术环境，因此会影响训练效果。相对而言，虚拟手术模拟器借鉴先进的计算机图形学技术能提供一个较真实的手术环境，并且具有可重复使用的特点，已然成为外科手术培训的未来发展趋势。在真实的内窥镜手术中，视觉反馈和触觉反馈是外科医生的主要信息来源，因此虚拟手术仿真也主要针对这两方面进行研究。视觉方面，手术模拟主要涉及计算机图形学、生物力学、弹性力学等学科；触觉方面，手术模拟主要涉及机械学。建立虚拟手术模拟器主要包含四部分：首先需要对人体器官进行几何建模，生成器官组织的形状；其次根据器官的种类和手术类型选择合适的物理变形方法，物理变形方法用于驱动几何模型计算变形；然后根据手术类型设计手术操作算法，如夹取、切割、缝合、灼烧等；最后通过显示器和力反馈设备从视觉和触觉上输出软组织的几何模型。目前，几何建模、视觉渲染和触觉渲染都已经比较成熟，物理变形方法和手术操作算法仍是学术界和工业界的研究热点。本课题来源于国家自然科学基金项目“可交互人体器官数字模型及虚拟手术研究”，是该项目的研究内容之一。本课题将致力于研究软组织变形和切割算法的若干问题。1.2 研究意义为了降低手术训练成本，并提高手术训练效果，基于VR技术的虚拟手术模拟器在近10 年发展十分迅速，工业界已经出现较为成熟的训练模拟器。1999年成立的Mentice公司于2001年设计了全世界第一台血管手术模拟器，产品称为VIST-Lab3，如图1所示。VIST-Lab模拟器能够模拟冠状血管造影术、心律管理、房间隔穿刺等手术。美国Simbionix公司的LAP Mentor模拟器可以模拟双手操作的腹腔镜手术4。LAP Mentor模拟器提供缝合、胆囊切除术、胆管造影术、阑尾切除术、切口疝修复术、肾切除术等的仿真和教学评估。使用者根据模拟器的提示按照真实手术的步骤进行训练，训练完成后模拟器对操作的准确性进行评分。此外，模拟器还提供真实手术的视频，训练人员能通过观看视频进行学习。图2是该模拟器的样图，该设备主要由1台主机、1台显示器、2台力反馈设备组成，和真实手术环境一致。针对不同的手术，模拟器提供了不同的操纵手柄，如医用剪刀、解剖刀、抓钳等。图3为瑞士VirtaMed公司设计的子宫颈手术模拟器HystSim5，HystSim模拟器可用于模拟子宫肌瘤剔除术、子宫内膜切除术、息肉切除术等。图1 Mentice VIST-Lab（引自文献3）图2 Simbionix LAP Mentor（引自文献4）图3 VirtaMed HystSim（引自文献5）目前，发达国家的医院已经使用虚拟手术模拟器进行训练，但我国并未大范围使用模拟器，一个重要原因是我国并未掌握手术模拟器的核心技术，不具备独立制造模拟器的能力，致使外国公司的手术模拟器售价高昂，医院不愿意在购买设备上投入巨资。但随着各种先进手术技术越来越多地在医院中被外科医生使用，同时传统方法不能满足训练的要求，医院必须借助手术模拟器进行训练。研究和发展虚拟手术仿真技术对我国公共卫生事业的发展、提高国民生活质量有着巨大的现实意义和推进作用。在虚拟手术仿真中，物理变形方法和手术操作算法一直是学术界和工业界的研究热点，尤其是在计算机图形学领域。手术模拟器一方面需要保证软组织变形的逼真性，一方面需要保证虚拟现实系统的实时性，逼真性和实时性具有一定的对偶性6。尽管使用某些物理精确的变形方法能保证软组织变形的真实性，但付出的代价就是计算时间的增加。一个实时的虚拟现实系统在视觉方面要求提供不低于24 FPS（Frames Per Second，帧/秒）的刷新速度，也就是两帧间隔的时间不能超过40 ms，计算时间太长会影响视觉的连贯性。同时，触觉渲染要求系统保证不低于1 kHz的刷新率。因此，选择或设计近似保持物理精确又能保证实时性的物理变形方法，在研究和设计虚拟手术模拟器时具有理论价值和实用意义。不同的手术类型需要研究和设计各自的手术操作算法，其主要原因来自于：手术模拟系统使用了不同的几何模型；几何模型和物理变形方法的耦合方式不同；手术过程的具体行为不同。不同的器官组织经常需要不同的几何模型，比如血管往往使用单层表面模型，即只使用一层三角形表面网格，对其进行切割和缝合时，要考虑其单层的特点。肝脏是典型的有内部实体的器官，对其进行仿真时一般采用体网格，在设计切割算法时，需要考虑其几何特征。使用物理变形方法时需要提供驱动几何模型变形的“骨架”，即物理模型，并由物理模型带动几何模型发生形变。物理变形方法有多种方式，而物理模型和几何模型的耦合也有多种方法，设计手术操作算法时有必要考虑这两点。不同的手术过程包含不同的手术操作，因此必须针对每种操作设计算法。考虑到上面提及的几点原因，研究和设计合适的手术操作算法很有必要。1.3 研究目标与内容基于研究背景，本文针对软组织变形和切割算法的若干问题进行研究，研究目标是设计和实现一种实时性和真实性较高的软组织变形方法，尝试表现软组织的材质属性，并提供能应用在虚拟手术模拟器中的夹取、按压和切割操作算法。本文的主要研究内容为：对现有物理变形方法和切割算法进行分析，改进传统的位置动力学（Position based dynamics），并使用扩展的位置动力学驱动物理模型发生形变，选择合适的方法建立软组织几何模型与物理模型的耦合，连接力反馈设备建立触觉反馈环境，使用力反馈设备控制手术器械进行夹取、按压和切割等手术操作。具体工作内容分为以下几部分：1. 分析现有物理变形方法和切割算法为满足虚拟手术仿真实时性和真实性的要求，本文分析了几种常用物理变形方法的原理和优缺点，最终选择了位置动力学作为本文使用的物理变形方法。本文还对切割方法分类别进行了调研和阐述，分析了各种方法的优缺点与适用范围。2. 扩展位置动力学模型，驱动软组织模型变形为了提高软组织变形的生物特性，本文改进了传统的位置动力学模型，增加了能量保持约束。本文使用改进的位置动力学模型作为物理变形方法，物理模型在物理变形方法的作用下产生形变。常用的物理模型包括四面体网格、六面体网格、无网格模型（Meshfree model）、元球模型（Blobby model）、骨架模型（Skeleton model）等，考虑到人体器官的几何形态，本文使用了四面体网格。和六面体网格、无网格模型、元球模型相比，对人体器官模型建立四面体网格产生的单元数目更少，这意味着计算物理变形的时间更少；人体器官往往不存在明显的骨架，因此骨架模型不适用于人体器官。选择四面体网格后，本文通过重心坐标（Barycentric coordinates）建立四面体网格和人体器官几何模型的耦合关系。3. 针对软组织模型设计夹取、按压和切割算法本文主要研究对肝脏、脾脏等内部有实体的器官的夹取和切割操作。本文首先研究和设计了对软组织模型的夹取和按压算法，然后实现了基于混合网格（三角形表面网格+四面体网格）的切割操作，又考虑到虚拟手术模拟器的工程需求，本文还实现了基于四面体网格切割操作。4. 建立触觉环境与系统集成为了提高手术模拟真实性，本文增加了力反馈设备，使用者可以通过力反馈设备操纵手术器械进行软组织的夹取、按压和切割操作。1.4 论文组织结构本文按照研究工作的背景、目标和研究内容，分为6章展开论述，各章具体内容如下：第一章 绪论。简要介绍虚拟手术仿真技术的产生背景及发展现状，分析虚拟手术仿真的关键技术。针对虚拟手术模拟器在我国的普及程度和发展趋势，阐述进行虚拟手术仿真方法研究的意义，并介绍本文的研究目标与主要研究内容等。第二章 关键技术及研究现状。介绍本文相关技术的背景，调研相关技术的国内外研究现状。首先介绍模拟人体器官软组织常用的几种物理变形方法，然后分析每种物理变形方法的优缺点和适用范围，之后介绍切割模拟相关技术的研究进展。第三章 扩展位置动力学模型和软组织变形。介绍位置动力学模型的原理，根据位置动力学模型的特点提出针对软组织模型的改进方法，增加能量保持约束。本章还介绍物理模型的生成方法与软组织几何模型的变形方法。第四章 基于位置动力学模型的手术操作模拟。首先介绍夹取、按压操作的模拟方法，然后介绍切割操作的模拟方法。对于切割操作，首先介绍基于混合网格的切割算法，考虑到虚拟手术模拟器的工程需求，又设计了基于四面体网格的切割算法。第五章 触觉环境的建立与交互。介绍建立触觉环境的步骤和触觉编程的方法，并阐述使用力反馈设备与软组织模型交互的方法。第六章 简易虚拟手术模拟器的设计与实现。介绍一个较为简单的虚拟手术模拟器的开发环境、系统架构，以及实现系统的重要数据结构。总结与展望：总结本文的研究工作，指出本文方法的缺陷和不足，并对虚拟手术仿真的相关技术和模拟器开发提出进一步的展望。第二章 关键技术及研究现状2.1 引言虚拟手术仿真涉及医学、计算机图形学、生物力学、弹性力学、应用数学、机械工程学等多领域的知识。鉴于本文的研究目标，即人体软组织器官的变形仿真与切割模拟，本文将相关技术限定在与人体器官有关的手术仿真中。在绪论中提到，计算软组织变形的物理变形方法和各种手术操作的模拟算法是虚拟手术仿真中的关键技术。本章针对这两项关键技术进行论述。2.2 虚拟手术仿真中的物理变形方法软组织变形常用的物理变形方法包括质点弹簧模型、有限元模型、无网格模型、锁子甲模型和位置动力学模型。2.2.1 质点弹簧模型质点弹簧模型（Mass-spring model）7是最简单和快速的物理变形方法。模型由一组粒子（即质点）构成，每个粒子保存各自的质量、位置和速度，相邻的粒子通过一根弹簧连接，这些弹簧将模型中所有粒子组成一个网络，如图4所示。每根弹簧认为是无重量的，并遵循胡克定律（Hookes law）。仿真时，对每个粒子计算外力和由弹簧产生的弹力/拉力的合力，并通过牛顿第二定律计算粒子的加速度，最终得到粒子的位置改变量和新位置。根据求解位置与速度方法的不同，质点弹簧模型又分为显式欧拉方法（Explicit Euler method）、龙格-库塔法（Runge-Kutta method）、隐式欧拉方法（Implicit Euler method）、韦尔莱积分法（Verlet Integration）。质点弹簧模型的结构简单、计算效率高，经常被用在实时系统中。图4 质点弹簧模型示意图（引自文献8）1999年，苏黎世联邦理工学院（ETH Zurich）的Bielser最早将质点弹簧模型引入虚拟手术中9。Bielser在论文中介绍了一种可用于手术仿真的三维软组织建模与切割的方法，文中使用质点弹簧模型模拟软组织的变形，物理模型为四面体网格，和之前使用张量乘积网格相比，使用四面体网格具有更高的拓扑灵活度。论文使用龙格-库塔法进行求解。系统可用二阶线性微分方程表示： （2.1）其中， 是个对角矩阵，且有， 为第个质点的质量，为第个质点在时刻的位置，是阻尼矩阵，是刚度矩阵，是一个维的向量，表示受到的外力。公式中对位置的二阶导数为质点的加速度，对位置的一阶导数为质点的速度。质点弹簧模型由于其简单、求解速度快的特点，常被用于软组织或柔性对象的变形仿真。2000年，Bielser等在文献10中使用了和文献9相同的物理变形方法。同一年年， Ganovelli等在另一种实现软组织切割模拟的文章中也使用了基于四面体网格的质点弹簧模型11。次年，Ganovelli等在文献12同样使用了质点弹簧模型。随后的几年中，不少文章开始采用有限元模型和无网格模型模拟软组织变形。直到2006年，Paloc等在动态重建软组织模型网格的研究中使用了质点弹簧模型13。同年，Steinemann等在文献14中使用了质点弹簧模型。2009年，Pan等使用质点弹簧模型模拟直肠变形15。此后，随着计算机硬件水平的提升和虚拟手术对变形真实度要求的提高，质点弹簧模型慢慢被淘汰。质点弹簧模型由于物理精度不高、系统稳定性差、变形受网格结构影响大的缺点而受到诟病。此外，质点弹簧模型的刚度系数不易设定，为系统调整到一个合适的状态需要花费大量精力。质点弹簧模型的最大问题是不具有体保持效果，使软组织变形缺乏物理真实性16，大部分研究者转而使用有限元模型模拟软组织变形。2.2.2 有限元模型有限元模型是一种物理变形真实性很高的物理变形方法，模型将软组织看作一块连续介质，并使用有限元方法（Finite element method）对连续介质求解关于运动的偏微分方程组，进而得到软组织的位置。有限元模型将软组织离散化到微小的单元格中，分别对每个单元格求解其应力（stress）和应变（strain），应力度量单元格在各个方向上的力，应变描述单元格中顶点的位置改变，其位置改变的基准是静止状态时的位置。应力和应变可以分别用二阶张量和表示，软组织的材质属性用弹性张量C表示，则应力和应变的关系可写成： （2.2）每个单元格中顶点的位置可以表示成应力、外力、材质密度的偏微分方程： （2.3）其中表示空间微分。公式（2.3）满足牛顿第二定律，公式左边表示质量乘以加速度，公式右边第1项表示单元格内部的应力，第2项表示受到的外力，如重力、碰撞力等。使用有限元模型模拟物理变形具有更高的物理真实性，自2000年以来，越来越多的研究人员使用有限元模型模拟软组织变形。在较早的研究中，软组织被认为是线性弹性的，因为当时的计算机硬件条件不足以实时计算非线性问题。使用线性弹性的文献包括17-19。尽管线性弹性的计算效率较高，但变形的真实性却有所缺失，不能模拟大形变。此外，由于线性弹性不能保持柯西应变的旋转不变性，故无法处理旋转的情况20。2009年，Jebkov等使用扩展有限元模型（Extended Finite Element Method，缩写为XFEM）模拟软组织变形，并解决了大形变的问题21。模型中任意点的位移u(x)可以表示为： （2.4）其中是单元格顶点的位移，是形函数，n是单元格顶点的数目。如图5所示，发生切割时，非连续浓缩函数被加入上式中： （2.5）切割后并发生大形变时，使用扩展有限元模型能真实地模拟物理变形。随着计算机硬件水平的提高，使用图形处理器（GPU）并行计算物理变形得到广泛应用。GPU多核心、多线程、支持浮点计算的特点十分适合处理大计算量的问题。使用GPU进行对有限元方法进行加速的文献包括22-27。图5 XFEM加入浓缩函数的示意图（引自文献21）2.2.3 无网格模型在无网格模型中，各个粒子之间没有固定的邻接关系，粒子的位置由当前状态紧邻的粒子决定。无网格模型能解决拓扑改变和大形变的问题，但无网格的计算量比有限元模型更大。2006年，Steinemann等在文献28中使用移动最小二乘法（Moving Least Squares，缩写为MLS）计算变形。在无网格模型中，每个点的位移可以由临近点的位移确定： （2.6）其中为当前点的位移，为临近点的位移，为物质坐标（material coordinate），为形函数，在文献28中，形函数为MLS方法。每个形函数对应一个权重，是点的支持半径，是临近点的物质坐标，关于权重函数的设定可以参考文献29。对比公式（2.4）和公式（2.6），不难发现，有限元模型和无网格模型的机理是相同的，只是形函数定义、网格结构不同。每个点的位置可以写成： （2.7）无网格模型能处理拓扑改变和大形变的情况，基于这个特性，Yang等提出混合使用有限元模型和无网格模型模拟软组织变形与切割的方法26。这种方法在相同材质的部位使用有限元模型进行模拟，在不同材质连接的部位使用无网格模型。2.2.4 锁子甲模型锁子甲模型（ChainMail Model）是一种遵循几何规则的物理变形方法，模型将当前点的位置作为输入，通过几何规则（如边长固定）计算点的可行区域，通过能量最小的条件确定点的静止位置，得到静止位置后再将改变量传给下一个点，直到网格中的点全部更新完成。2003年，Li等在文献30中提出了三维锁子甲模型，图6为使用锁子甲模型的脑部软组织变形效果。图6 使用锁子甲模型的脑部软组织变形（引自文献30）由于锁子甲模型属于几何方法，计算量比有限元模型和无网格模型的计算量小，但物理真实性比上述两种模型的真实性差，因此锁子甲模型不能胜任物理真实性要求高的应用。2.2.5 位置动力学模型本文使用位置动力学模型模拟软组织变形，并对其进行了改进，本节主要介绍位置动力学模型的起源和入门知识，具体的算法细节在第三章介绍。在质点弹簧模型中，求解牛顿运动方程的数值方法有一种韦尔莱积分法（Verlet Integration），设运动点在时刻t的位置为x(t)，运动时间为t，则通过泰勒公式得到： （2.8） （2.9）两式相加得到： （2.10）从公式（2.10）中可以看出，速度并不需要显式出现，而上一时刻的位置需要保存，由这种想法出发，忽略掉公式（2.10）的高阶无穷小后，可以得到： （2.11） （2.12）由于速度隐含在当前位置和上一次位置中，可以通过两个位置直接计算速度，进而可以节省计算速度的开销，这是位置动力学模型的数学基础。位置动力学模型（Position based dynamics）的原型来自2001年Jakobsen设计的物理引擎Fysix31，2006年，Mller等在文献32中正式提出位置动力学模型，并包括拉伸约束、弯曲约束、自碰撞约束等用于模拟布料的几何约束，这些约束在文献33和34中提出。图7所示为使用位置动力学模型的布料仿真效果图。2007年，Kubiak等使用韦尔莱积分法模拟绳子的变形35，尽管文章并没有提到位置动力学模型，但具体的计算方法和位置动力学模型是一致的。图7 使用位置动力学模型的布料仿真（引自文献32）位置动力学模型解决了质点弹簧模型过度调整（overshooting）的问题，由于其稳定、简单、快速的特点被工业界广泛采用，如NVidia公司的物理引擎PhysX36、Havok公司的布料引擎Havok Cloth37、Autodesk公司的Maya nCloth工具38。在上述引擎中，大部分引擎仅用于布料仿真，但PhysX引擎支持体模型的物理变形，因此有研究人员直接使用PhysX引擎模拟软组织的变形，如文献39，40。PhysX引擎支持二维网格的撕裂操作，但不支持三维网格的拓扑改变，本文的贡献之一是对其拓扑改变的情况进行了扩展。位置动力学模型不能表现变形体的材质属性，本文增加了能量保持约束使软组织具有一定的材质属性。2.3 虚拟手术中的切割方法虚拟手术中的切割方法与几何模型和物理变形方法都紧密相关。发生切割时，一是需要改变几何模型中顶点的拓扑关系、二是需要改变物理模型中顶点的拓扑关系。对于不同的物理变形方法，变形效果与顶点拓扑关系的依赖性不同。在有限元模型中，当生成薄片形状（sliver）的单元时，物理变形方法的稳定性开始降低，最后会产生错误的变形结果；质点弹簧模型对网格拓扑的要求比较低，网格改变对物理变形方法的稳定性影响不大；无网格模型比较适合拓扑改变的情况，因为原始网格中不存在固定的拓扑关系；锁子甲模型不支持拓扑改变的情况；位置动力学模型支持二维变形体的拓扑改变，如布料的撕裂。常用的物理模型包括：四面体网格、六面体网格、无网格结构。本节参考文献41的分类方法，从切割方法的实现思想对其进行分类，如图8所示，可分为：直接删除法、顶点挤压法、沿边界分裂法、单元细分法、单元复制法和无网格切割法。由于无网格切割法存在固定的网格结构，故未在图8出现。图8 切割方法分类图（引自文献41）2.3.1 直接删除法早期的研究使用直接删除单元的方法达到切割效果，如文献17。这种方法删除解剖刀碰到的单元。为了使新的