CN106934189A - 外科手术软组织变形的仿真方法及装置 - Google Patents

Abstract

本发明适用于生物力学仿真技术领域，提供了外科手术软组织变形的仿真方法及装置，包括：在外科手术仿真过程中，对软组织中的每个粒子i，计算其形变梯度对所述形变梯度F_i进行对角化，得到奇异值矩阵并计算对角化后所述粒子i的第一应力计算所述粒子i的每个邻域粒子j对所述粒子i的作用力根据计算结果对所述软组织进行形变处理。本发明实施例基于点的方法来重现外科手术交互过程中人体器官的生物力学反应，在软组织的形变中，其每个粒子的位置变化仅仅取决于材料性质，能够达到很好的模型精确度。

Description

外科手术软组织变形的仿真方法及装置

技术领域

本发明属于生物力学仿真技术领域，尤其涉及外科手术软组织变形的仿真方法及装置。

背景技术

近年来，尽管各种医疗手段都取得了巨大的进步，但是对于情况较为严重的病患，手术依然是治愈率最高的一种。虚拟现实技术的手术仿真可以有效地缩短实习医生进行手术训练的周期，其主要应用包括术前医学影像处理、解剖结构可视化、生物力学模拟等。

在虚拟现实技术的手术仿真领域，交互式生物力学仿真是一项关键技术，近年来发展起来的无网格方法阐述了不使用网格在传统意义上的优势，这对于几何形状复杂的研究对象来说处理方式更为灵活。此外，由于拉格朗日框架对质点分布的广泛兼容性，无网格法可以较容易地控制有限应力及大尺度形变，自五十年代有限元法出现之后，无网格方法已成为最具扩展性的仿真方法。然而，无网格方法仍然存在着高精度几何模型带来的高计算复杂度的问题同时，在外科手术过程中软组织变形的情况下，将不得不使用网格定义以控制可能出现的大尺度变形。

发明内容

有鉴于此，本发明实施例提供了交互式生物力学的仿真方法及装置，以解决无网格法无法处理外科手术仿真中软组织大尺度变形的问题。

第一方面，提供了一种外科手术软组织变形的仿真方法，包括：

在外科手术仿真过程中，对软组织中的每个粒子i，计算其形变梯度其中，为粒子i的基本矩阵A_i，所述j为所述粒子i的邻域粒子，所述w_ij为权值，所述u_ij为初始时刻粒子i和邻域粒子j的间距，所述x_ij为当前时刻粒子i和邻域粒子j的间距；

对所述形变梯度F_i进行奇异值分解得到奇异值矩阵并计算对角化后所述粒子i的第一应力所述I为单位矩阵，所述λ和μ均为拉梅系数；

计算所述粒子i的每个邻域粒子j对所述粒子i的作用力根据计算结果对所述软组织进行形变处理，所述v_i为所述粒子i的体积。

第二方面，提供了一种外科手术软组织变形的仿真装置，包括：

第一计算单元，用于在外科手术仿真过程中，对软组织中的每个粒子i，计算其形变梯度其中，为粒子i的基本矩阵A_i，所述j为所述粒子i的邻域粒子，所述w_ij为权值，所述u_ij为初始时刻粒子i和邻域粒子j的间距，所述x_ij为当前时刻粒子i和邻域粒子j的间距；

第二计算单元，用于对所述形变梯度F_i进行奇异值分解得到奇异值矩阵并计算对角化后所述粒子i的第一应力 所述I为单位矩阵，所述λ和μ均为拉梅系数；

第三计算单元，用于计计算所述粒子i的每个邻域粒子j对所述粒子i的作用力根据计算结果对所述软组织进行形变处理，所述v_i为所述粒子i的体积。

本发明实施例基于点的方法来重现外科手术交互过程中人体器官的生物力学反应，在软组织的形变中，其每个粒子的位置变化仅仅取决于材料性质，能够达到很好的模型精确度。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例提供的外科手术软组织变形的仿真方法的实现流程图；

图2是本发明另一实施例提供的外科手术软组织变形的仿真方法的实现流程图；

图3是本发明实施例提供的实验中采用的带有纹理的表面渲染模型图；

图4是本发明实施例提供的实验中采用的连续介质粒子化表达示意图；

图5是有限元法和本发明实施例方案的变形效果图；

图6是基于点的方法与本发明实施例方案在体积上的比较示意图；

图7是基于点的方法与本发明实施例方案在运行时间上的比较示意图；

图8是基于点的方法与本发明实施例方案在视觉效果上的对比示意图；

图9和图10分别是肝和胃在大变形下的仿真效果示意图；

图11是本发明实施例提供的外科手术软组织变形的仿真装置的结构框图。

具体实施方式

以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透切理解本发明实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本发明。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本发明的描述。

图1示出了本发明实施例提供的外科手术软组织变形的仿真方法的实现流程，在本方案中，采用了基于点的方法来重现外科手术交互过程中人体器官的生物力学反应，其实现流程详述如下：

在S101，在外科手术仿真过程中，对软组织中的每个粒子i，计算其形变梯度其中，为粒子i的基本矩阵A_i，所述j为所述粒子i的邻域粒子，所述w_ij为权值，所述u_ij为初始时刻粒子i和邻域粒子j的间距，所述x_ij为当前时刻粒子i和邻域粒子j的间距。

在S102中，对所述形变梯度F_i进行奇异值分解得到对角矩阵并计算所述粒子i的应力其中，所述I为单位矩阵，所述λ和μ均为拉梅系数，由材料属性决定。

在S103中，计算所述粒子i的每个邻域粒子j对所述粒子i的作用力根据计算结果对所述软组织进行形变处理，其中，所述v_i为所述粒子i的体积。

通过图1所示方案，就可以基于点的方法来重现外科手术交互过程中人体器官的生物力学反应，在软组织的形变中，其每个粒子的位置变化仅仅取决于材料性质，能够达到很好的模型精确度。

此外，作为本发明的一个实施例，为了使现有模型高度可控和无条件稳定，可以采用基于位置的方法去快速处理交互过程，如图2所示，所述方法还包括：

在S201中，计算所述粒子i的位置修正值其中，所述M为质量矩阵，所述C_i(x)为所述粒子i周围约束方程。

在S202中，根据所述位置修正值对所述粒子i的位置进行修正。

在本发明实施例中，对手术器械表面网格进行简化采样，利用这些点代替手术器械模型来实现快速碰撞检测，在碰撞过程中，基于现有技术中解决保证物体不相互渗透问题的研究工作，在所述粒子i和所述粒子j的碰撞过程中，通过C(x_i，x_j)＝|x_ij|-r≥0来进行碰撞限制，其中，所述r为在所述粒子i和所述粒子j的碰撞方向上，为所述粒子i和所述粒子j预设的碰撞限制距离，来防止物体间互相渗透现象的发生。此外，对于加速物体间的碰撞检测，本方案中采用空间哈希网格来加速粒子邻域信息的获取速度。

作为本发明的一个实施例，为了防止相互渗透现象的发生，在粒子i和粒子j的碰撞过程中，还需要计算因摩擦力而产生的正切方向的位移：

其中，和分别为碰撞的粒子i和粒子j当前可能到达的位置，x_i和x_j分别为碰撞开始时刻粒子i和粒子j所在的位置，是粒子i和粒子j接触点的法线，所述

由此，对每个粒子i来说，因摩擦而产生的位移如下：

其中，u_k和u_s分别是动摩擦力系数和静摩擦力系数，所述w_i为粒子i权值，所述w_j为粒子j权值，所述d为任意常量，当粒子相对速度小于牵引力阈值，上式可以直接忽略正切方向位移。

对于与粒子i碰撞的粒子j来说，其因摩擦而产生的位移如下：

由于大多数的生物软组织是不可压缩的，因此，体积守恒在可变形物体的动态仿真中有重要意义。对于四面体网格模型，在本发明实施例中，在每一个四面体网格上，对网格顶点(p₁，p₂，p₃，p₄)定义一个体积限制：

其中，p_i，j＝p_i-p_j，所述V₀为四面体的初始体积。

接下来，通过一系列实验来验证本发明实施例提供的生物力学仿真方法的准确性、鲁棒性以及不可压缩性。首先，对实验环境的说明如下：

下述实验均在Windows 7操作平台上完成，由C++实现并且通过OpenGL渲染。实验运行在Intel Xeon CPU，2.40GHz，8GB memory，NVIDIA GTX 650Ti的硬件平台上。采用的模型为带有纹理的表面渲染模型以及粒子化表达的连续介质力学计算模型，其中，图3为实验中采用的带有纹理的表面渲染模型，图4为实验中采用的连续介质粒子化表达。

1、对模型仿真的准确性验证：

通过对比SOFA(Simulation Open Framework Architecture)中有限元仿真方法同等条件下肝脏变形标记点位移结果，来验证本方案所提模型在动力学计算方面的准确性。图5左图为有限元法的变形效果，右图为本方案的变形效果，可以看出，本方案在模型仿真的精确度上逼迫有限元方法。

2、对模型的不可压缩性的验证：

体积守恒可以用来保证手术对象的不可压缩性。图6是仿真过程中，基于点的方法与本方案在体积方面的比较，图7是这两种方法在运行时间上的对比，通过对比可以看出，本方案提出的模型在没有增加运行时间的情况下，具有较好的体积守恒特点，从而更加逼近真实值。图8为两种方法视觉效果对比渲染图。

3、对模型鲁棒性的验证：

可变形固体的仿真方法，如质点弹簧法，FEM法等在大变形或大时间步长的条件下经常会出现数值计算不稳定的情况。在本方案的模型中，全部点的位置是由基于点的方法来计算的，这仅仅取决于材料性质。此外，在本方案中引入基于位置的方法，来实现交互以及体积守恒，上述方案的采用使得本方案模型在大时间下以及大形变下也能够达到稳定性要求。图9和图10是肝和胃分别在大变形下仿真的情况。虽然超大变形很少在实际情况中出现，但这是一种检验稳定性的重要手段。

应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本发明实施例的实施过程构成任何限定。

对应于上文实施例所述的外科手术软组织变形的仿真方法，图11示出了本发明实施例提供的外科手术软组织变形的仿真装置的结构框图。为了便于说明，仅示出了与本实施例相关的部分。

参照图11，该装置包括：

第一计算单元1101，在外科手术仿真过程中，对软组织中的每个粒子i，计算其形变梯度其中，为粒子i的基本矩阵A_i，所述j为所述粒子i的邻域粒子，所述w_ij为权值，所述u_ij为初始时刻粒子i和邻域粒子j的间距，所述x_ij为当前时刻粒子i和邻域粒子j的间距；

第二计算单元1102，对所述形变梯度F_i进行奇异值分解得到奇异值矩阵并计算对角化后所述粒子i的第一应力 所述I为单位矩阵，所述λ和μ均为拉梅系数；

第三计算单元1103，计算所述粒子i的每个邻域粒子j对所述粒子i的作用力根据计算结果对所述软组织进行形变处理，所述v_i为所述粒子i的体积。

可选地，所述装置还包括：

第四计算单元，计算所述粒子i的位置修正值其中，

修正单元，根据所述位置修正值对所述粒子i的位置进行修正。

可选地，所述装置还包括：

碰撞限制单元，在所述粒子i和所述粒子j的碰撞过程中，通过|x_ij|-r≥0进行碰撞限制，其中，所述r为在所述粒子i和所述粒子j的碰撞方向上，为所述粒子i和所述粒子j预设的碰撞限制距离。

可选地，所述装置还包括：

第五计算单元，在所述粒子i和所述粒子j的碰撞过程中，计算所述粒子i因摩擦而产生的位移：

并计算所述粒子j因摩擦而产生的位移：

其中，u_k和u_s分别是动摩擦力系数和静摩擦力系数， xi*和xj*分别为所述粒子i和所述粒子j当前可能到达的位置，xi和xj分别为碰撞开始时刻所述粒子i和所述粒子j所在的位置，n＝xij*/xij*是所述粒子i和所述粒子j接触点的法线，所述

可选地，所述装置还包括：

体积限制单元，对于四面体网格模型，在每一个四面体网格上，对网格顶点(p₁，p₂，p₃，p₄)定义体积限制其中，p_i，j＝p_i-p_j，所述V₀为四面体的初始体积。

所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。实施例中的各功能单元、模块可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。另外，各功能单元、模块的具体名称也只是为了便于相互区分，并不用于限制本申请的保护范围。上述系统中单元、模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

在本发明所提供的实施例中，应该理解到，所揭露的装置和方法，可以通过其它的方式实现。例如，以上所描述的系统实施例仅仅是示意性的，例如，所述模块或单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通讯连接可以是通过一些接口，装置或单元的间接耦合或通讯连接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。

所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实施例的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)或处理器(processor)执行本发明实施例各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例各实施例技术方案的精神和范围。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种外科手术软组织变形的仿真方法，其特征在于，包括：

在外科手术仿真过程中，对软组织中的每个粒子i，计算其形变梯度其中，为粒子i的基本矩阵A_i，所述j为所述粒子i的邻域粒子，所述w_ij为权值，所述u_ij为初始时刻粒子i和邻域粒子j的间距，所述x_ij为当前时刻粒子i和邻域粒子j的间距；

对所述形变梯度F_i进行奇异值分解得到奇异值矩阵并计算对角化后所述粒子i的第一应力所述I为单位矩阵，所述λ和μ均为拉梅系数；

计算所述粒子i的每个邻域粒子j对所述粒子i的作用力根据计算结果对所述软组织进行形变处理，所述v_i为所述粒子i的体积。

2.如权利要求1所述的方法，其特征在于，所述方法还包括：

计算所述粒子i的位置修正值其中， ${\mathrm{\λ}}_i=\frac{C_{i\left(x\right)}}{\▿C_{i\left(x\right)}{M}^{-1}\▿C_i{\left(x\right)}^T};$

根据所述位置修正值对所述粒子i的位置进行修正。

3.如权利要求1所述的方法，其特征在于，所述方法还包括：

在所述粒子i和所述粒子j的碰撞过程中，通过|x_ij|-r≥0进行碰撞限制，其中，所述r为在所述粒子i和所述粒子j的碰撞方向上，为所述粒子i和所述粒子j预设的碰撞限制距离。

4.如权利要求1所述的方法，其特征在于，所述方法还包括：

在所述粒子i和所述粒子j的碰撞过程中，计算所述粒子i因摩擦而产生的位移：

${\mathrm{\&Delta;x}}_i=\frac{w_i}{w_i+w_j}\left\{\begin{array}{c}\mathrm{\&Delta;}{x}_{\&perp;},\left|\mathrm{\&Delta;}{x}_{\&perp;}\right|<u_{S}d\\ \mathrm{\&Delta;}{x}_{\&perp;}\&CenterDot;\min (\frac{u_{k}d}{\left|{\mathrm{\&Delta;x}}_{\&perp;}\right|},1),\end{array}\right.$ 其他，

并计算所述粒子j因摩擦而产生的位移：

${\mathrm{\&Delta;x}}_j=\frac{w_i}{w_i+wj}{\mathrm{\&Delta;x}}_i,$

其中，u_k和u_s分别是动摩擦力系数和静摩擦力系数， xi*和xj*分别为所述粒子i和所述粒子j当前可能到达的位置，xi和xj分别为碰撞开始时刻所述粒子i和所述粒子j所在的位置，n＝xiij*/*是所述粒子i和所述粒子j接触点的法线，所述

5.如权利要求1所述的方法，其特征在于，所述方法还包括：

对于四面体网格模型，在每一个四面体网格上，对网格顶点(p₁,p₂,p₃,p₄)定义体积限制 $C(p_1,p_2,p_3,p_4)=\frac{1}{6}(p_{1,2}\&times;p_{3,1})\&CenterDot;p_{4,1}-V_0,$ 其中，p_i,j＝p_i-p_j，所述V₀为四面体的初始体积。

6.一种外科手术软组织变形的仿真装置，其特征在于，包括：

第一计算单元，用于在外科手术仿真过程中，对软组织中的每个粒子i，计算其形变梯度其中，为粒子i的基本矩阵A_i，所述j为所述粒子i的邻域粒子，所述w_ij为权值，所述u_ij为初始时刻粒子i和邻域粒子j的间距，所述x_ij为当前时刻粒子i和邻域粒子j的间距；

第二计算单元，用于对所述形变梯度F_i进行奇异值分解得到奇异值矩阵并计算对角化后所述粒子i的第一应力 所述I为单位矩阵，所述λ和μ均为拉梅系数；

第三计算单元，用于计计算所述粒子i的每个邻域粒子j对所述粒子i的作用力根据计算结果对所述软组织进行形变处理，所述v_i为所述粒子i的体积。

7.如权利要求6所述的装置，其特征在于，所述装置还包括：

第四计算单元，用于计算所述粒子i的位置修正值其中， ${\mathrm{\&lambda;}}_i=\frac{C_{i\left(x\right)}}{\&dtri;C_{i\left(x\right)}{M}^{-1}\&dtri;C_i{\left(x\right)}^T};$

修正单元，用于根据所述位置修正值对所述粒子i的位置进行修正。

8.如权利要求6所述的装置，其特征在于，所述装置还包括：

碰撞限制单元，用于在所述粒子i和所述粒子j的碰撞过程中，通过|x_ij|-r≥0进行碰撞限制，其中，所述r为在所述粒子i和所述粒子j的碰撞方向上，为所述粒子i和所述粒子j预设的碰撞限制距离。

9.如权利要求6所述的装置，其特征在于，所述装置还包括：

第五计算单元，用于在所述粒子i和所述粒子j的碰撞过程中，计算所述粒子i因摩擦而产生的位移：

${\mathrm{\&Delta;x}}_i=\frac{w_i}{w_i+w_j}\left\{\begin{array}{c}\mathrm{\&Delta;}{x}_{\&perp;},\left|\mathrm{\&Delta;}{x}_{\&perp;}\right|<u_{S}d\\ \mathrm{\&Delta;}{x}_{\&perp;}\&CenterDot;\min (\frac{u_{k}d}{\left|{\mathrm{\&Delta;x}}_{\&perp;}\right|},1),\end{array}\right.$ 其他，

并计算所述粒子h因摩擦而产生的位移：

${\mathrm{\&Delta;x}}_j=\frac{w_i}{w_i+wj}{\mathrm{\&Delta;x}}_i,$

其中，u_k和u_s分别是动摩擦力系数和静摩擦力系数， xi*和xj*分别为所述粒子i和所述粒子j当前可能到达的位置，xi和xj分别为碰撞开始时刻所述粒子i和所述粒子j所在的位置，n＝xij*/xij*是所述粒子i和所述粒子j接触点的法线，所述

10.如权利要求6所述的装置，其特征在于，所述装置还包括：

体积限制单元，用于对于四面体网格模型，在每一个四面体网格上，对网格顶点(p₁,p₂,p₃,p₄)定义体积限制 $C(p_1,p_2,p_3,p_4)=\frac{1}{6}(p_{1,2}\&times;p_{3,1})\&CenterDot;p_{4,1}-$ $V_0,$ 其中，p_i,j＝p_i-p_j，所述V₀为四面体的初始体积。