虚拟手术系统中形变物体的实时模拟系统

摘要

一个虚拟手术系统中的实时模拟软组织物体形变系统，属于图形处理技术领域。本发明中：边界元模块首先为形变物体建立物理模型，划分物体表面网格单元，对边界积分方程进行离散化并求解出每个网格单元的形变向量。形状匹配模块根据划分好的网格单元，对每一个网格单元的形变位置和初始位置建立一一映射关系，根据网格单元的位置和物体材料属性建立运动学方程，计算物体的形变回复形状。有限状态机在软组织物体形变过程中，对物体的受力状态进行分析，并通过基本状态之间的切换来控制当前使用的形变运算模块。本发明实现了在交互式系统中实时并且精确对软组织形变进行模拟，特别是虚拟手术系统中手术器械与软组织交互变得更精确，更快速。

说明书

技术领域

本发明涉及的是一种图形处理技术领域的实时模拟物体形变系统，特别是一种虚拟手术系统中形变物体的实时模拟系统。

背景技术

随着科学技术的高速发展，高科技医疗设备的不断涌现为医疗的现代化提供了越来越多的帮助。虚拟手术平台的建立在医生的培训、手术导航等方面起到重要作用。在虚拟手术训练中，医生需要通过控制虚拟环境中的手术器械来对虚拟的人体器官、血管、软组织等进行操作。在这种交互过程中，模拟虚拟器官、组织在手术器械操作下的形变情况是虚拟手术系统中重要的一部分。一个精确，实时的模拟软组织器官形变系统能够很大程度上提升手术系统的真实性和实用性。

尽管对软组织形变模拟的研究在虚拟现实和计算机图形学领域有着很长的历史，这些研究成果在虚拟手术领域的应用并不广泛。这其中的主要原因是创建一个形变模拟系统框架需要同时具有处理形变精确性模块和处理形变实时性模块。一方面，在计算机上进行手术训练不同于其他诸如视频游戏之类的人机交互应用，训练者是以掌握一种真实的操作技术为目的的，因此他们需要训练系统能够提供和真实手术中非常类似甚至完全一样的反馈信息，这其中包括视觉信息(软组织的移动和变形)和触觉信息(通过器械感受到的反馈力)。形变模拟系统大多数只具有其中之一的模块，因此很难完整的模拟手术中的器官形变。

经过对现有技术的文献检索发现，Doug James等在《Computer GraphicsProceedings》(计算机图形学学报)，Annual Conference Series，ACMSIGGRAPH99(ACM SIGGRAPH99年度会议系列)的65-72页上发表的“Artdefo，accurate real time deformable objects”(Artdefo，精确实时的形变物体)中提出了一个基于边界元法的形变物体实时模拟系统Artdefo。该系统从边界元计算的角度对物体形变进行处理，有效的解决了模拟精确性的问题。但是，在系统的实时性能方面该模拟系统还存在着计算效率低，系统交互性不够的不足。

综上所述，设计一个同时兼顾精确性和实时性的形变模拟系统，对虚拟手术中的器官形变模拟至关重要。

发明内容

本发明的目的在于克服现有技术中的不足，提供一种虚拟手术系统中形变物体的实时模拟系统，使其能够在实时的条件下精确模拟出软组织器官的受力情况和形变情况。

本发明是通过以下技术方案实现的，本发明包括：边界元模块、形状匹配模块、有限状态机，其中：边界元模块和形状匹配模块共同模拟软组织物体的形变，有限状态机分析物体的受力状态并且切换形变的计算模式，其中：

所述边界元模块首先为形变物体建立物理模型，划分物体表面网格单元，对边界积分方程进行离散化并求解出每个网格单元的形变向量；边界元模块是形变模拟系统中负责系统精确性的模块，它由有限状态机在需要精确形变模拟的情况下进行调用，并接受形变物体的当前的空间信息作为输入信息，通过模块中的偏微分方程处理和形变向量计算，对形变物体进行精确模拟。

所述形状匹配模块根据划分好的网格单元，对每一个网格单元的形变位置和初始位置建立一一映射关系，根据网格单元的位置和物体材料属性建立运动学方程，计算物体的形变回复形状；形状匹配处理模块是形变模拟系统中负责系统实时性的模块，它由系统状态机在不需要精确形变模拟的情况下进行调用，并接受形变物体的当前位置状态作为输入信息，通过模块中的点集映射和运动模拟处理，对形变物体进行实时处理。

所述有限状态机在软组织物体形变过程中，对物体的受力状态进行分析，并通过基本状态之间的切换来控制当前使用的边界元模块和形状匹配模块。有限状态机通过调用边界元模块，把系统的当前状态输入该模块中进行形变的精确计算。在系统资源不足或者需要交互性很高的情况下，状态机切换到形状匹配模块，以该模块的对内存和CPU的低消耗来实现实时性的特点。

所述的有限状态机控制系统模块之间的通信和协作。该状态机把一个软组织物体的受力形变过程划分为三个基本状态：平衡状态、形变状态和回复状态。在平衡状态中，软组织物体的形状处于受力稳定的状态，形状不发生改变。在形变状态中，作用在物体表面某边界单元的外力被检测到，并根据高斯分布函数分布在表由外部输入指定的局部范围的单元内。有限状态机此时调用边界元模块开始计算物体的表面形变和受力大小。在回复状态中，物体不再受外力作用，但物体的表面形状不处于平衡状态，形变物体在内力的作用下回复成平衡状态的形状。有限状态机通过分析形变物体表面单元的受力情况，在三个基本状态之间切换，并调用该状态所对应的边界元模块和形状匹配模块计算表面形变。

所述有限状态机在检测到一个新外力加载在物体表面时把当前状态转到形变状态，并调用边界元模块精确计算形变和反馈力。在外力从物体表面移除时，有限状态机的当前状态转到回复状态，并调用形状匹配模块计算形变的恢复过程。当形变物体回到平衡位置时，有限状态机的当前状态转到平衡状态，此时不调用计算模块，释放之前计算所用的内存资源。

所述边界元模块和形状匹配模块的系统复杂度不同，在系统上层调度过程中所分配的时间也不同。边界元模块是建立在实时求解偏微分方程基础上的计算模块，模拟形变的准确度高，反馈力分析精确，所耗费的系统CPU资源和内存资源都很高。形状匹配模块是一个建立在动力学弹簧基础上的处理模块，可以快速匹配物体的形变形状和平衡形状，模拟形变中的形状恢复过程，系统处理速度很快，不需要使用额外的内存空间。

所述的边界元模块首先为形变物体建立物理模型。物体形变的数学模型是基于弹性力学中线弹性物体的纳维方程来描述了物体的形变位移和所受外力之间的关系。对作用在物体整个体区域内的纳维方程，把Navier方程进行边界化，从而把弹性物体整个体空间的三维问题转化为局限于物体表面的二维问题。

所述的边界元模块在为形变物体建立物理模型之后，接收从系统外部读取的医学数据模型文件，把物体的表面划分为一系列互不重叠的三角片。这种表面划分方式和虚拟手术所需要用的图形渲染接口诸如VTK，OpenGL或者DirectX的模型划分方式是一致的，这样也避免了额外的数据结构和存储开销。在划分了物体表面网格之后，就把一个连续的物体表面划分为了一系列个离散的边界单元。然后对其中的每一个单元应用上面提到的边界积分方程，根据不同的边界条件，把边界单元的位移或受力作为未知量，从而建立起了一个可由计算机进行求解的线性方程组。

所述边界元模块使用形变基向量的方法来提高系统处理方程组的速度。形变基向量方法的关键点是在实时计算之前对系数矩阵进行预计算，从而把计算复杂性很高的解线性方程组操作转化为矩阵向量的乘法操作。在刷新形变向量基的时候，只需把矩阵向量中的非零元素或值发生变化的元素提取出来进行矩阵向量乘法操作，在维数很大的情况下可以大大的提高计算效率。在求得了形变基向量之后，根据形变物体的线性特点，对这些向量进行简单的线性缩放操作，就可以快速求得不同外力作用于物体同一位置的情况下物体不同的形变情况。

所述形状匹配模块，当外力从软组织物体上移除时，软组织需要回复到从前的形状，形状匹配模块对这一过程的模拟使用了无网格形状匹配技术。从空间结构来说，这一技术不需要对形变向量通过方程组的方式求解，也不需要额外知道整个网格的拓扑信息，只需要对每一个边界单元建立一个从当前空间位置到初始空间的一一映射。在这一过程中，把形变物体的表面单元看作是一个没有相互作用的粒子系统。

所述形状匹配模块，为了模拟回复状态下的物体形变过程，在对所有边界单元做了映射匹配之后，对每一个映射建立一个含阻尼的弹簧系统，每一个形变后的元素都被映射中的弹簧拉向它的初始位置。因为阻尼力的存在，每一个单元都将最终在空间上匹配到形变物体的初始位置，而回复过程的剧烈程度和快慢程度则可以通过调节刚度系数和阻尼力来实现。

在本发明中，有限状态机被引入了形变模拟系统，和以往单纯依赖边界元模拟技术或者表面形状匹配技术的模拟系统相比，状态机很好的结合了两种技术的优点(边界元技术的精确性和形状匹配技术的实时性)，通过不同模块间的调度实现了系统实时性和精确性之间的平衡。状态机通过调用边界元模块，把系统的当前状态输入该模块中进行形变的精确计算。在系统资源不足或者需要交互性很高的情况下，状态机可以很快切换到形状匹配模块，以该模块的对内存和CPU的低消耗来实现实时性的特点。通过状态机的这种自动调度，系统可以根据当前的具体环境在精确和实时方面进行不同的侧重，该系统应用在虚拟手术中，可以显著的提高虚拟手术系统中器官形变的真实性，同时提升手术系统的交互性能。

附图说明

图1本发明系统结构框图。

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

如图1所示，在本发明的形变模拟系统中，系统的基本结构包括两个不同的形变处理模块(边界元模块和形状匹配模块)和一个控制模块间调度并且管理信息输入输出的状态机。

所述边界元模块建立物体弹性形变偏微分方程，划分物体表面网格单元，对边界积分方程进行离散化并求解出每个网格单元的形变向量；边界元模块是形变模拟系统中负责系统精确性的模块，它由有限状态机在需要精确形变模拟的情况下进行调用，并接受形变物体的当前的空间信息作为输入信息，通过模块中的偏微分方程处理和形变向量计算，对形变物体进行精确模拟。

所述形状匹配模块根据划分好的网格单元，对每一个网格单元的形变位置和初始位置建立一一映射关系，根据网格单元的位置和物体材料属性建立运动学方程，计算物体的形变回复形状；形状匹配处理模块是形变模拟系统中负责系统实时性的模块，它由系统状态机在不需要精确形变模拟的情况下进行调用，并接受形变物体的当前位置状态作为输入信息，通过模块中的点集映射和运动模拟处理，对形变物体进行实时处理。

所述有限状态机在软组织物体形变过程中，对物体的受力状态进行分析，并通过基本状态之间的切换来控制当前使用的边界元模块和形状匹配模块。有限状态机通过调用边界元模块，把系统的当前状态输入该模块中进行形变的精确计算。在系统资源不足或者需要交互性很高的情况下，状态机切换到形状匹配模块，以该模块的对内存和CPU的低消耗来实现实时性的特点。

所述的边界元模块，完成以下的处理：

a)建立物体弹性形变偏微分方程

边界元模块首先为形变物体建立物理模型。物体形变的数学模型是基于弹性力学中线弹性物体的Navier方程：

Nu+b＝0

其中N是二次偏微分算子，u表示物体区域内点的位移，b表示物体所受的体力(比如重力)。这个方程描述了物体的形变位移和所受外力之间的关系。对作用在物体整个体区域内的Navier方程，使用偏微分方程求解中常用的Green-Gauss理论和Kelvin基本解把Navier方程进行边界化(见《边界元理论及应用》第二章，北京理工大学出版社，2002年9月)，从而把弹性物体整个体空间的三维问题转化为局限于物体表面的二维问题。

$C_{\mathrm{lk}}\left(P^{\prime }\right)u_k\left(P^{\prime }\right)=\int \Gamma [u_{\mathrm{lk}}^{*}(P^{\prime },Q^{\prime })-u_k\left(Q^{\prime }\right)p_{\mathrm{lk}}^{*}(P^{\prime },Q^{\prime })]\mathrm{d\Gamma }\left(Q^{\prime }\right)+\int \Omega u_{\mathrm{lk}}^{*}(P^{\prime },Q)f_k\left(Q\right)\mathrm{d\Omega }\left(Q\right)$

在这里u*项和p*项是Navier方程的基本解。

b)划分物体表面网格单元

边界元模块把物体的表面划分为一系列互不重叠的三角片。这种表面划分方式和虚拟手术所需要用的图形渲染接口诸如VTK，OpenGL或者DirectX的模型划分方式是一致的，这样也避免了额外的数据结构和存储开销。

在划分了物体表面网格之后，就把一个连续的物体表面划分为了n个离散的边界单元。然后对其中的每一个单元应用上面提到的边界积分方程，根据不同的边界条件，把边界单元的位移u(x)或受力p(x)作为未知量，从而建立起了一个n个方程，n个未知数的线性方程组。

HU＝GP

其中H和G是3n×3n的系数矩阵，U和P是包含3n个元素的向量。

对这个方程组进行规范化的移项处理，把u(x)和p(x)中的未知量移到方程等号的左边，已知量移到右边，同时系数矩阵H和G中的相应列也需要进行交换。得到新的线性方程组

AX＝CX

其中A和C分别是交互后的H和G矩阵，X在交换后包含的是所有的未知量，X包含的是所有的已知量。

c)实时求解网格单元的形变向量

为了能够实时的求解所有边界单元在每一个时间步内的形变向量，在本发明中我们使用了形变基向量的方法来提高方程组的求解速度。形变基向量方法的关键点是在实时计算之前对系数矩阵进行预计算，从而把计算复杂性很高的解线性方程组操作转化为矩阵向量的乘法操作。预计算首先计算出系数矩阵A的逆矩阵A^-1。然后把方程组变换为

X＝A^-1CX

这样，在刷新形变向量基的时候，只需把X中的非零元素或值发生变化的元素提取出来进行矩阵向量乘法操作，在X的维数很大的情况下可以大大的提高计算效率。

在求得了形变基向量之后，根据形变物体的线性特点，对这些向量进行简单的线性缩放操作，就可以快速求得不同外力作用于物体同一位置的情况下物体不同的形变情况。

2.形状匹配模块，做如下处理：

a)点集的映射

当外力从软组织物体上移除时，软组织需要回复到从前的形状，形状匹配模块对这一过程的模拟使用了无网格形状匹配技术。从空间结构来说，这一技术不需要对形变向量通过方程组的方式求解，也不需要额外知道整个网格的拓扑信息，只需要对每一个边界单元建立一个从当前空间位置到初始空间的一一映射。在这一过程中，把形变物体的表面单元看作是一个没有相互作用的粒子系统。

b)回复状态下的运动模拟

为了模拟回复状态下的物体形变过程，在对所有边界单元做了映射匹配之后，我们对每一个映射建立一个含阻尼的弹簧系统，每一个形变后的元素都被映射中的弹簧拉向它的初始位置。弹簧的回复力可以通过如下公式计算得出

f＝-k[p(x)-p(O)]-f_damp

在公式中，p(O)为物体的初始位置，p(x)为物体的当前位置，k为所设弹簧的劲度系数，f_damp为弹簧的阻尼力。这个方程描述了物体形变回复的时的受力情况。因为阻尼力f_damp的存在，每一个单元都将最终在空间上匹配到形变物体的初始位置，而回复过程的剧烈程度和快慢程度则可以通过调节刚度系数k和阻尼力f_damp来实现。

本实施例在CPU为Pentuim M 1.5GHz，内存为1.0GB的计算机中实现，编程语言为C++，具体参数情况如表1所示：

表1实施例中各种参数情况

模型编号 形状属性 四面体个数 顶点数A Tumor 960 482B Kidney 2584 1294

1、首先从vtk文件中导入模型A、B的网格数据，在虚拟手术系统的渲染模块中通过vtk渲染通道渲染在虚拟场景中。为了模拟模型的形变，为两个模型分别建立前文所述的有限状态机，并把状态都设为平衡状态。

2、分别对边界元模块和形状匹配模块进行预计算处理。在边界元模块的预计算中，对每个物体通过上文所述的边界元方法建立物体弹性形变偏微分方程，划分物体表面网格单元，对边界积分方程进行离散化并求解出系数矩阵的逆矩阵。在形状匹配模块的预计算中，记录模型所有网格单元的初始位置。

3、对模型的某一个部分区域内的边界单元施加外力，计算状态机改变到形变状态，并调用边界元计算模块。在边界元计算模块中，通过外力作用的大小和边界单元的位置确定形变条件向量X，再由边界元模块中预先计算出的系数矩阵逆矩阵A^-1和条件矩阵C，根据X＝A^-1CX求得模型每一个边界单元的形变及受力向量。并根据这些形变向量刷新模型显示在场景中的形状，根据受力刷新场景力反馈设备中的受力。

4、当外力从模型表面移开时，有限状态机从形变状态转换到回复状态，并调用形状匹配模块。在形状匹配模块中，根据预计算过程中记录的初始位置如上文所述建立网格当前位置和初始位置的一一映射。并通过模块中的回复弹簧计算回复形变位置。

5、当物体完全恢复到初始位置是，有限状态机把模块切换回平衡状态。当有新的外力被施加在模型上时，重新从步骤3开始形变的模拟过程。

在该环境下的精确形变模拟能够达到实时的效果，与传统方法相比，性能也有了很大的提高，分别对应于上述情况A和B，具体结果见表2。

表2实施例结果及与传统方法的比较情况

情况本发明系统处理时间(ms)传统系统处理时间(ms)性能提高率(％)A315038.0B12517327.7