一种改进的虚拟手术中人体软组织形变方法

何巍，王魏平，师为礼，苗语，何飞，杨华民

（长春理工大学　计算机科学与技术学院，长春　130022）

一种改进的虚拟手术中人体软组织形变方法

何巍，王魏平，师为礼，苗语，何飞，杨华民

（长春理工大学计算机科学与技术学院，长春130022）

虚拟手术是将计算机虚拟现实技术应用于现代医学，利用医学图像数据，重构出虚拟人体软组织模型，利用触觉交互设备，在视觉和触觉上为使用者提供手术场景的模拟和仿真。其中，软组织形变的仿真和建模是虚拟手术系统核心的研究内容之一。针对虚拟手术中软组织形变仿真，提出了一种改进的弹簧—质点模型。在正四边形网格拓扑结构表面模型中，增加了结构弹簧、剪切弹簧和弯曲弹簧，使表面模型具有体模型的特征。该模型继承了传统弹簧—质点模型的建模快、原理简单、仿真速度快等优点，同时还具有控制形变区域的能力，从而提高了仿真模型的精确度，满足交互系统的真实性和稳定性需求。实验表明该方法有效增强了形变仿真的体积特征，适合于软组织形变较大时的实时仿真，满足了虚拟手术系统实时性的要求。

虚拟现实；弹簧—质点模型；软组织形变；碰撞检测

随着计算机科学技术的不断发展，虚拟现实技术（Virtual Reality Technology VRT）已被广泛应用于现代科学研究的各个领域［1，2］，其中利用虚拟现实技术进行生物医学仿真是当前的一个研究热点，虚拟手术系统是其中的典型应用。

虚拟手术系统是对医学数据进行可视化并实现实时交互，模拟手术过程中组织器官变形、感官反馈等可能遇到的各种现象的虚拟现实应用系统。其中在软组织器官形变方面，生物力学领域的研究提出了很多复杂的数学模型，能够比较准确地描述软组织的形变，软组织大多数物理模型都是建立在弹性力学理论上的。Santanu等人的基于有限元方法（Finite Element Model，FEM）提出了具有生物特征的物理形变模型，虽然精确度很高，但是网格数量很多，导致模型实时性差［3］；Zhong等人提出具有物理化学特征的表面模型，但很难进行复杂的虚拟手术操作［4］；陈卫东等改进了基于表面正六边形模型，虽然计算量小，实时性高，但是仿真精度有待提高［5］；王瑞艺提出的基于无网格SPH的物理模型，虽然能够动态地反映软组织的实时形变，但交互方面需要进一步改善［6］；李艳东，朱玲等人结合ALMSDM的特点提出了顶点法向量局部更新与预计算策略，提高了实时性，但在大范围变形的情况下，稳定性不足［7］；Wang等用边界元法建立模型，虽对模型的边界进行了离散，简化了计算，但稳定性较差［8］。因此，在保证模型形变精确度的同时，提高交互的实时性，是我们虚拟手术操作中需要解决的首要问题。传统的基于表面网格拓扑结构的物理模型，存在两方面问题，没有体模型的特征或者形变真实性不足。

本论文在保证实时性的同时，在虚拟的人体软组织表面模型中设置了结构弹簧、剪切弹簧和弯曲弹簧，使模型具备了体模型的特征，从而提高了仿真精度，达到了更好的仿真效果。同时保证形变实时性高，很好地满足虚拟手术中交互过程的真实性、实时性和稳定性需求。

1　软组织形变物理模型

弹簧-质点模型（Mass Spring Method，MSM）是一种经典的物理建模方法，其主要思想是把仿真的软组织用质点离散，质点与质点之间用满足胡克定律的弹簧连接。该模型通常基于一定的几何拓扑结构，形变过程是当一个质点受外力作用时，产生的应力作用在其它相邻的质点上，同时通过弹簧传递力，带动其相邻的质点运动。模型的建立比较简单，形变的计算量较低，能够很好满足实时性的要求。

1.1软组织的生物力学特征

对人体软组织器官（如血管、肝脏、皮肤、气管等）形变特性的研究属于生物力学范畴［9］。随着生物力学研究的发展，软组织的生物力学特征通常表现为非线性、粘弹性、不均匀性和各向异性等［10］。因此，虚拟手术仿真系统中的软组织建模必须考虑上述生物力学特性，并根据实际情况做相应的简化处理。通常采用弹性模量、阻尼系数等物理量来表示软组织的各种生物力学特征。

1.2弹簧—质点模型建立

假设弹簧—质点模型由n个质点组成，在任意时刻t每个质点Ni的运动方程都满足牛顿第二定律，因此每个质点Ni都具有如下的动力学微分方程（1）：

其中，mi表示第i个质点质量；xi表示质点Ni的形变位移；ci表示阻尼系数；kij表示弹簧的弹性系数；lij表示质点Ni和Nj间弹簧的初始长度，lij0表示形变后的长度；Fext表示质点Ni所受的外力的综合。

在弹簧—质点模型模拟人体软组织形变时，力的作用导致质点发生位移。力的构成，除了外力之外，还有弹簧的变形产生的内力。不同的弹簧布局在仿真过程中会产生不同的力，质点也会产生不同的位移，使模型产生不同的变形效果。当采用弹簧—质点模型进行软组织变形仿真实验时，核心之一是建立合适的几何拓扑结构。

本文对人体的皮肤进行仿真实验，采用面模型。因为人体的皮肤质量分布相对均匀，其网格可均匀划分。网格划分的密度应根据仿真需求的精度和交互的实时性来确定。网格的拓扑结构从不同的精度考虑采用图1（b）、（c）、（d）三种方式［11］。

图1　人体软组织面模型拓扑结构

为提高仿真精度，我们选用图1（d）中拓扑结构作为物理模型。弹簧可分为3类：结构弹簧、剪切弹簧和弯曲弹簧，如图2所示。

图2　弹簧—质点模型的三种弹簧类型

（1）结构弹簧：质点Ni，j和Ni+1，j之间和质点Ni，j和Ni，j+1之间的弹簧，是防止质点之间过度的拉压形变。

（2）剪切弹簧：质点Ni，j和Ni+1，j+1之间和质点Ni，j+1和 Ni+1，j之间的弹簧，是阻止结构的斜向过度变形，给结构一个剪切刚度。

（3）弯曲弹簧：质点Ni，j和Ni+2，j之间和质点Ni，j和Ni，j+2之间的弹簧，是为了防止结构过度的不自然弯曲。

2　软组织模型算法分析

2.1动力学方程

基于上述弹簧—质点模型拓扑结构，结合受外力情况下模型的运动规律，对质点建立动力学方程，通过求方程的解来实现软组织的形变。

力学模型的建立基于以下两点：（1）模型能够达到实时的仿真效果；（2）在现有的硬件条件下能够实现相关算法的复杂度要求。现有的基于物理模型的仿真算法一般采用两种基本方法：力法和能量法。为了满足仿真的实时交互需求，本文采用力法。以下是对弹簧—质点模型的动力学分析［12，13］：

由牛顿第二定律可得，模型中每个质点Ni的运动微分方程由它所受到的合力（内力和外力组成）决定，对于质点N 有：

其中，Fext（i）是质点Ni所受的外力，Fint（i）是质点Ni所受的内力（弹性力、阻尼力等）。根据胡克定律，质点Ni所受弹簧的拉（压）力如公式（3）所示：

其中，Kijs是质点Ni和Nj间的弹性系数，Iij是该时刻弹簧的矢量，Xi和Xj是质点Ni和Nj的位置。

质点Ni和Nj之间的弹簧阻尼力如下：

其中，Kijd是质点Ni和Nj之间的弹簧阻尼系数。

根据公式（3）和（4），可以得到质点Ni所受的内力：

其中，K是与质点Ni连接的质点总数。

一个由N个质点组成的模型，所受的内力由结构力，弯曲力和剪切力组成；可以表述为：

其中，Fstruct、Fshearing、Fbending分别是与质点Ni连接的弹簧产生的结构力、剪切力和弯曲力，k1、k2、k3分别表示与质点Ni相连的结构弹簧、剪切弹簧和弯曲弹簧的数量。

一般情况下，质点Ni所受到的外力包括重力和拉力，则受到的外力可以表述为：

对于上述的动力学微分方程，模型x并不是现实给出的，为了方便求解，需要对二阶微分方程进行处理，将其转化为一阶常微分方程进行求解。则，每个质点Ni的速度和位移都可以用以下公式表示：

2.2模型的数值分析

首先考虑一阶常微分方程的初值求解问题：

其中f（ ）x，y为已知函数；y0为初值。

数值积分的原理就是，在区间［a，b］上取n+1个节 a=x0＜x1＜x2＜…＜xn=b。 hi=xi-xi-1（1，2，…） 表示相邻两个节点的间距，称为步长。这些hi可以不相等，我们取相等考虑。则将这些节点离散化，就可以将初值问题转化成离散变量的相关问题。把相应问题的解yn作为y（xn）的近似值。这样解得yn就是上述初值问题在节点xn的数值解，不同的离散方法产生不同的结果。

一阶常微分方程初值问题的几种常用的数值积分算法包括欧拉方法、中点方法以及四阶龙格—库塔方法等。对于质点数一定的软组织模型，形变的实时性和稳定性取决于算法。主要从以下几个方面来考虑：

（1）步长：表示为了达到一定数值稳定性所需要的时间离散化程度，表明方法的数值稳定性；

（2）迭代次数：反映了采用积分方法的计算量，影响总的计算速度；

（3）误差：数值求解精度受误差影响，包括截断误差、舍入误差等；当确定步长后，微分方程的阶次越高，截断误差越小；表1列出了以上三种数值算法的截断误差。

表1　三种数值算法的截断误差

从计算精度上来看，步长h确定的情况下，四阶龙格—库塔法的计算精度最高，但是执行一步所需的计算量也最大，影响了仿真形变的实时性。只从单步积分来看，步长设置越小，模型离散化程度就越高，局部截断误差也就越小。但是离散化程度越大，同样意味着在一定范围内需要的步数越多，不仅增加了计算量，还有可能引起舍入误差的积累过大，导致精度变低。

3　碰撞检测

碰撞检测称为接触检测或者干涉检测，基本目的是确定两个或两个以上物体彼此之间是否发生接触或穿透并计算出碰撞发生的位置。在虚拟手术中，需要碰撞检测的场景主要有手术器械与人体器官之间、人体器官的自碰撞和手术器械的自碰撞。这些是整个虚拟手术系统中最重要的部分之一，决定着下一步的软组织变形、切割、缝合等的实现。

如图3所示，质点t0时刻在软组织外部，而在下一个时间步长t0+Δt时刻已经在软组织内部，则质点和软组织在t0+Δt和t0+Δt之间的时刻发生了碰撞。

图3　碰撞对象

目前常见的碰撞检测算法主要有：层次包围盒法［14］和空间分解法［15］。

层次包围盒的核心思想是用体积略大而几何特征简单的包围盒近似表述复杂的几何对象，只对包围盒相交时包裹的对象进行相交实验。此外，通过引入树状层次结构快速删除不发生碰撞的部分，减少了不必要的相交实验，提高检测效率。层次包围盒法应用较为广泛，适用于复杂环境下的碰撞检测。根据所采用的包围体类型的不同可对层次包围盒进行区分，比较典型的有沿坐标轴的包围盒AABB、方向包围盒OBB和包围球等。

空间分解法的核心思想是将全部空间分解成体积相等的小单元格，只对同一单元格或相邻单元格的几何对象进行相交测试。比较典型的有八叉树、BSP树等。

虚拟手术仿真系统需要很好地满足交互实时性的要求，因此，我们使用AABB层次包围盒的改进算法作为本文虚拟手术仿真系统的碰撞检测算法。

AABB方法［16］沿坐标轴的包围盒是一种使用广泛的碰撞检测算法，包含几何对象且边平行于坐标轴的最小六面体。因此只需要六个标量就能描述一个AABB方法。

AABB树是基于AABB方法自底至上动态更新生成的层次结构二叉树。AABB树算法的核心是通过遍历两个碰撞物体的AABB树来判断该位置是否发生碰撞。步骤如下：

步骤1：构造碰撞刚体和软组织的AABB树；

步骤2：检测最大包围盒是否相交；

步骤3：递归处理刚体和软组织的AABB树，如果发现子树没有发生相交，停止并得出结论没有发生碰撞。如果发现子树相交，则进一步处理它的子树直到到达叶子节点，并最终得出结论。

软组织形变主要是软组织被提拉或者挤压，所以软组织形变的更新主要是AABB树的更新。发生形变时，只有接触点周围的一些质点发生了位移形变。相对于整个组织，变形只发生在局部区域，所以模型的拓扑结构变化较小，更新时只对发生了形变的质点进行更新。由于各个节点形变程度不同，不能用统一函数来描述。因此，先更新发生形变处节点的AABB包围盒，再对整个包围盒树进行更新。父结点的包围盒更新可通过两个更新了的子节点表示。对于整个物体包围盒的更新，如果检测到上一层包围盒没有相交，其后的包围盒就无需再更新［16］。

4　虚拟手术仿真结果

使用PC INTEL E3 3.10GHZ，8GBRAM，NVIDIA NVS 300显卡计算机，以OPENGL图形库为基础，使用VS2010和C++开发环境实现软组织仿真的局部动态弹性形变。

由于弹簧—质点模型是模拟皮肤表面的形变现象，只是一种近似仿真，而且模型的参数也无法选取真实软组织的生物力学参数，所以通过大量实验，对质点质量、弹簧的弹性系数和阻尼系数等进行对比后选取如下：时间步长0.01，质点质量0.01，重力加速度9.8，弹簧的结构弹簧、剪切弹簧和弯曲弹簧的弹性系数分别为2，2，1.2；阻尼系数分别为0.4，0.4，0.24。虚拟体弹簧的弹性系数为表面弹簧弹性系数的两倍。

以平面模型为例进行仿真，如图4所示。分别选取了225个质点，450个面片（a）和625个质点1250个面片（b）的面模型进行试验。

图4　不同质点的面模型

以图4（a）为模型，对固定形变区域和动态形变区域进行对比。当外力较小时，若形变区域范围很大，则计算形变花费的时间多，所以我们采用固定形变区域，如图5（a）所示；当外力大时，若形变区域范围很小，则形变效果较差，所以我们采用动态形变区域，如图5（b）所示。

图5　形变效果图

以100个时间步长为例，虚拟一个较小的外力，选用图5（b）所示模型。在进行形变计算时，选取一个合适的临界值，顶点的形变大于该值时，激活下一个领域的顶点。若在某个时间步长内，碰撞点在四个方向的形变量均小于该临界值，则退出循环。在相同的外力作用下，不同的临界值不仅影响形变区域的大小，同时还影响时间步长。临界值设置越小形变越精细，相应的形变计算步骤越多，时间越长。本文选取临界值为0.001，虚拟外力为5N的情况下，分别用欧拉法，中点法，四阶龙格—库塔法进行仿真实验，计算时间如表2所示，仿真效果如图6所示。

表2　网格模型三种积分方法计算时间比较

实验结果表明网格质点越多，实时性越差。从3种不同的积分效果去考虑，可以发现，实验的真实性越来越好，但实时性越来越差。在网格精度相对较高的情况下，欧拉法和中点法能够达到很好的仿真效果且计算时间能够达到实时性的需求，但四阶龙格—库塔法的计算时间明显提高。综合比较后，在满足实时性的同时选择中点法，从而达到较好的仿真精度和效果。

图6　实验结果图

5　结论

本文主要针对虚拟手术仿真过程中实时性和真实性两方面的需求，采用较能反映软组织形变“体特征”的质点—弹簧体模型实现了软组织形变，并在此基础上提出了更加精确的拓扑结构。为了增强形变过程中的稳定性，设置了结构弹簧、剪切弹簧和弯曲弹簧；基于模型形变的动力学方程采用欧拉方法，中点法和四阶龙格-库塔法作为单个质点的运动算法，并通过对虚拟软组织器官的提拉和按压手术实验。结果表明，本文提出的方法能够很好地满足模型的精度和仿真实时性要求，增强了形变的仿真体积感，适合于软组织形变较大时的实时仿真。

［1］ 孙连荣，姜元章，任玲.虚拟现实技术及其在高校教学中的应用［J］.石油教育，2003（5）：41-44.

［2］ 戴雯惠.虚拟现实技术的应用现状及发展［J］.信息技术，2006，35（6）：170-174.

［3］Santanu M，Amit R.Simulation of hip fracture in sideways fall using a 3D finite element model of pelvis-femur-softtissuecomplexwithsimplified representation of whole body［J］.Medical Engineering&Physics，2007，29：1167-1178.

［4］ Zhong Yongmin，Bjian S，et al.Virtual reality simulation of surgery with haptic feedback based on the boundary element method［J］.Computers and Structures，2007，85：331-339.

［5］ 陈卫东，赵成龙.虚拟手术中软组织形变建模及力反馈算法研究［J］.中国生物医学工程学报，2013，32（1）：113-118.

［6］ 王瑞艺.虚拟手术中力反馈的研究和实现［D］.郑州：华北水利大学，2014.

［7］Li Yandong，Zhu Ling.Modeling and simulation of softtissuedeformationbasedonlocaldynamic model［J］.Computer Science，2013，40（10）：283-288.

［8］ Wang P，Becker A A，Jones I A，et al.Virtual reality simulation of surgery with haptic feedback based on the boundary element method［J］.Computers and Structures，2007，85：331-339.

［9］冯元祯.生物力学［M］.北京：科学出版社，1983：124-158.

［10］Meredith M，Maddock S.Real-time inverse kinematics：the return of the jacobian［J］.Department of ComputerScience，TheUniversityofSheffield. 2004：51-60.

［11］ 乔冰.基于质点—弹簧模型的人体软组织形变技术的研究［D］.哈尔滨：哈尔滨工程大学，2008.

［12］ 褚莲娣.基于质点一弹簧模型的布模拟方法［J］.嘉兴学院学报，2002，14（3）：57-60.

［13］ 潘振宽，高波.手术方针中基于质点一弹簧模型的人体组织变形仿真［J］.青岛大学学报，2003，18（3）：9-14.

［14］ 朱元峰，孟军，谢光华，等.基于复合层次包围盒的实时碰撞检测研究［J］.系统仿真学报，2008，20（2）：372-377.

［15］ 邹益胜，丁国富，许明恒，等.实时碰撞检测算法综述［J］.计算机应用研究，2008，25（1）：8-12.

［16］ 赵成龙.虚拟手术中软组织建模与力反馈算法研究［D］.秦皇岛：燕山大学，2013.

An Improved Method on Soft Tissue Deformation in Virtual Surgery

HE Wei，WANG Weiping，SHI Weili，MIAO Yu，HE Fei，YANG Huamin
（School of Computer Science and Technology，Changchun University of Science and Technology，Changchun 130022）

Virtual surgery provides the surgery simulation for users on vision and force.It applies computer virtual reality in modern medicine，uses medical image data and reconstructs virtual human soft tissue models.And the simulation and modeling of soft tissue deformation has become one of the most important research contents in virtual surgery system.An improved mass spring model is proposed in this paper for the soft tissue deformation simulation.Bending spring，shear spring and structure spring are added onto the surface model of quadrilateral mesh topological structural for soft tissues，which will make the surface model having volume feature.This model not only has the advantages of traditional spring mass model such as rapid modeling，simple principle and quick simulation，but also has the ability to control the deformation regions，which improves the accuracy of the simulation model，and meets the requirements of the realness and stability of interactive system.Experimental results show that the method is suitable for the real-time simulation of soft tissues with large deformation.It can effectively enhance the volume of deformation simulation and meet the requirements of real-time interactive systems.

virtual reality technology；mass spring model；soft tissue deformation；collision detection

TP317.4

A

1672-9870（2015）06-0118-05

2015-10-30

何巍（1978-），女，博士研究生，讲师，E-mail：hw@cust.edu.cn

杨华民（1963-），男，教授，博士生导师，E-mail：yhm@cust.edu.cn