CN111341449B - 一种虚拟血管介入手术训练的模拟方法 - Google Patents

Abstract

本发明公开了一种虚拟血管介入手术训练的模拟方法，包括以下步骤：(1)使用基于点的方法对血管内部构建形变模型；(2)使用基于位置动力学法对血管内部节点施加距离约束、体积守恒约束和弹性势能守恒约束三种约束条件；(3)使用无网格移动最小二乘法构建从血管内部到血管表面的映射模型；本发明的方法不仅无需网格拓扑结构的初始划分和重构，甚至避免了采用网格法在产生大形变时的网格扭曲或畸形问题，而且考虑了血管的弹性特性和体积守恒特性，并提高了软组织形变模型的计算精度、稳定性和实时性。

Description

一种虚拟血管介入手术训练的模拟方法

技术领域

本发明涉及软组织力触觉模拟方法，具体涉及一种虚拟血管介入手术训练的模拟方法。

背景技术

心脑血管疾病是一种具有较高致残率和死亡率的疾病，主要通过血管介入手术进行诊断和治疗。传统血管介入手术训练是基于真实物体的解剖学，但该训练方法存在着些许不足与问题。随着计算机技术、仿真技术和虚拟现实技术的发展，虚拟血管介入手术训练逐渐成为一种新型的治疗心脑血管疾病的方法，为外科医生提供了良好的血管介入手术训练平台，他们可以根据具体需求进行反复练习，从而提高其手术水平，但搭建虚拟血管介入手术训练系统的一个重要技术难点就是对虚拟血管的高精度模拟。

在现有的虚拟软组织模拟方法中，有限元模型将问题求解域划分为若干个互不重叠的单元，通过单元形函数和节点插值函数建立刚度矩阵进而构建描述软组织材料力学特性的近似函数模型，由函数模型获得各单元内节点位移就可用来表现软组织的形变，因此其形变仿真精度最高，但在形变过程中网格拓扑结构会不断重组甚至扭曲，从而导致计算复杂度太高和计算量太大；质点弹簧模型将软组织离散化为一系列由弹簧相连的质点，质点同时受到弹簧的弹力作用和阻尼器的阻尼力约束，然后，根据牛顿第二定律对每个节点建立拉格朗日运动方程，最后通过求解此力学方程得到各节点的形变位移量，该模型结构简单易于实现，但模型稳定性差和仿真精度低。

综上，现有的有限元模型和质点弹簧模型等网格模型存在着当出现大形变时网格扭曲、拓扑结构重构时计算复杂、速度慢，且计算精度降低等问题。

发明内容

发明目的：本发明的目的在于提供一种能够避免采用网格法在产生大形变时的网格扭曲或畸形，提高软组织形变模型的计算精度、稳定性和实时性的虚拟血管介入手术训练的模拟方法。

技术方案：本发明的虚拟血管介入手术训练的模拟方法，包括以下步骤：

(1)使用基于点的方法对血管内部构建形变模型；

(2)使用基于位置动力学法对血管内部节点施加距离约束、体积守恒约束和弹性势能守恒约束三种约束条件；

(3)使用无网格移动最小二乘法构建从血管内部到血管表面的映射模型。

步骤(1)中，所述使用基于点的方法对血管内部构建形变模型，是指根据血管的医学影像数据，使用基于点的方法将血管内部离散化为一系列基于四面体的点云，并将血管的体积均匀分配给节点，通过获得节点的应力和应变来计算节点的形变位移向量。

所述通过获得节点的应力和应变来计算节点的形变位移向量，包括以下步骤：

(a)使用多项式核函数ω_ij来衡量中心节点i对其邻节点j的影响能力，具体计算公式为：

式中，h表示节点i的支撑半径，r表示节点i与其邻节点j之间的距离；

(b)计算节点i位移向量u_i的空间导数

具体计算公式为：

式中，

和/>

分别表示节点i的位移场u＝(m,n,p)^T中横坐标所对应的值m的空间导数、纵坐标所对应的值n的空间导数以及竖坐标所对应的值p的空间导数；

(c)通过求得的

可得到节点i处的应变ε_i和应力σ_i，具体计算公式为：

σ_i＝Cε_i

式中，J_i表示节点i的雅可比矩阵，I表示单位矩阵，C表示弹性矩阵，其值取决于弹性材料的杨氏模量和泊松比；

(d)根据连续介质力学理论，估计出节点i周围所存储的应变能U_i具体计算公式为：

式中，Ω表示节点i的支撑域，ν_i表示节点i的体积；

(e)通过对邻节点位移向量求导计算出每个邻节点的所受力f_j，从而得到中心节点i的所受内力f_i，其值为所有邻节点所受力f_j的负总和，具体计算公式为：

(f)通过对下式进行数值积分计算求解血管内部节点i的形变位移向量：

式中，M_i、u_i、v_i、a_i分别表示节点i的质量、位移、速度和加速度，t表示迭代时间，f_ext、f_i分别表示节点i所受的外力和内力。

步骤(b)中，所述节点i的位移场u＝(m,n,p)^T中横坐标所对应的值m的空间导数

的计算公式为：

式中，x_ij＝x_i-x_j，

表示矩量矩阵A_i的逆，m_i,m_j分别表示节点i和j的位移场u＝(m,n,p)^T横坐标所对应的值，ω_ij表示节点i和j间的权重，由多项式核函数计算，x_ij表示节点i的位置x_i与节点j的位置x_j之间的距离位移向量，||x_ij||表示节点i与节点j之间的距离。

步骤(b)中，节点i的位移场u＝(m,n,p)^T中纵坐标所对应的值n的空间导数

和节点i的位移场u＝(m,n,p)^T中竖坐标所对应的值p的空间导数/>

的计算方法与节点i的位移场u＝(m,n,p)^T中横坐标所对应的值m的空间导数/>

的计算方法相同。

步骤(2)中，所述距离约束C_distance(x₁,x₂)的函数公式为：

C_distance(x₁,x₂)＝|x₁-x₂|-d₀

式中，d₀表示节点x₁,x₂间的初始距离，根据距离约束的约束条件得到的节点修正因子Δx_i(i＝1,2)计算公式为：

式中，

表示节点x_i(i＝1,2)质量的倒数。

步骤(2)中，所述体积守恒约束C_volume(x₁,x₂,x₃,x₄)的函数公式为：

式中，V₀表示虚拟四面体单元(x₁,x₂,x₃,x₄)的初始体积，根据体积守恒约束的约束条件得到的节点修正因子Δx_i(i＝1,2,3,4)计算公式为：

式中，

表示节点x_i(i＝1,2,3,4)质量的倒数。

步骤(2)中，所述弹性势能守恒约束C_energy(x₁,x₂,x₃,x₄)的函数公式为：

式中，x₀表示四面体单元(x₁,x₂,x₃,x₄)的重心位置，k_i表示连接节点x_i,x₀的虚拟弹簧的弹性系数，d_i0表示节点x_i,x₀间的初始距离，根据弹性势能守恒约束的约束条件得到的节点修正因子Δx_i(i＝1,2,3,4)计算公式为：

式中，

表示节点x_i(i＝1,2,3,4)质量的倒数。

步骤(2)中，所述无网格移动最小二乘法中定义的映射函数为：

式中，u^h(X)为质点X的场函数u(X)的近似函数，Φ(X)为质点X支持域内的形函数，U_S为一个n维向量，用以描述支持域内所有节点的形变位移。

采用基函数p(X)和权函数W(X)来构造形函数，基函数p(X)、权函数W(X)的函数公式分别为：

p^T(X)＝[1,x,y,z]

式中，

为无量纲权函数影响半径，r_i为节点i的影响域半径。

有益效果：本发明与现有技术相比，其有益效果在于：(1)使用基于点的方法构建虚拟血管内部形变模型，通过计算应力和应变来获得节点的形变位移向量，在不损害计算效率的同时保证了计算精度和实时性；(2)使用基于位置动力学法对虚拟血管内部节点实施距离约束、体积守恒约束和弹性势能守恒约束三种约束条件，模拟了血管的弹性特性和体积守恒特性；(3)使用无网格移动最小二乘法构建从血管内部到血管表面的映射模型来渲染软组织形变，无需网格拓扑结构的初始划分和重构，甚至避免了采用网格法在产生大形变时的网格扭曲或畸形问题，提高了形变模型的稳定性。

附图说明

图1是本发明的流程图；

图2是本发明中虚拟血管模型结构图；

图3是本发明中距离约束投影示意图；

图4是本发明中弹性势能守恒约束投影示意图；

图5是本发明中从血管内部到血管表面的映射示意图。

具体实施方式

下面结合具体实施方式和说明书附图对本发明做进一步详细介绍。

如图1所示，本发明包括三个部分：使用基于点的方法对虚拟血管内部构建形变模型以此控制软组织的运动、使用基于位置动力学法对血管内部节点施加约束以及使用无网格移动最小二乘法构建从血管内部到血管表面的映射模型以此渲染软组织的形变；具体实施步骤如下：

(1)根据医学影像数据，使用基于点的方法对血管内部构建形变模型，将血管内部区域离散化为点云模型，并将血管的体积均匀分配给这些节点，通过获得节点的应力和应变计算其形变位移向量。使用OpenGL将CT扫描采集的血管医学图像数据进行三维几何重现，处理图像数据将血管表面离散化为一系列基于三角形网格的质点，将血管内部离散化一系列基于四面体的节点，如图2所示，然后，使用基于点的方法对血管内部建模；

其中，通过获得节点的应力和应变来计算节点的形变位移向量，包括以下步骤：

(a)使用多项式核函数ω_ij来衡量中心节点i对其邻节点j的影响能力，具体计算公式为：

式中，h表示节点i的支撑半径，r表示节点i与其邻节点j之间的距离；

(b)计算节点i位移向量u_i的空间导数

具体计算过程为：

x_ij＝x_i-x_j

式中，x_ij＝x_i-x_j，

表示节点i的位移场u＝(m,n,p)^T中横坐标所对应的值m的空间导数，/>

表示矩量矩阵A_i的逆，m_i,m_j分别表示节点i和j的位移场u＝(m,n,p)^T横坐标所对应的值，ω_ij表示节点i和j间的权重，由多项式核函数计算，x_ij表示节点i的位置x_i与节点j的位置x_j之间的距离位移向量，||x_ij||表示节点i与节点j之间的距离；

矩量矩阵A_i可由下式得到：

同样地，

也由上述方法计算求出，因此即可得到/>

式中，

分别表示节点i的位移场u＝(m,n,p)^T中纵坐标所对应的值n的空间导数以及竖坐标所对应的值p的空间导数；

(c)通过求得的

可得到节点i处的应变ε_i和应力σ_i，具体计算公式为：

σ_i＝Cε_i

式中，J_i表示节点i的雅可比矩阵，I表示单位矩阵，C表示弹性矩阵，其值取决于弹性材料的杨氏模量和泊松比；

(d)根据连续介质力学理论，估计出节点i周围所存储的应变能U_i，具体计算公式为：

式中，Ω表示节点i的支撑域，ν_i表示节点i的体积；

(e)应变能本质上是一个关于中心节点i和其所有邻节点j位移向量的函数，因此通过对邻节点位移向量求导就可计算出每个邻节点的所受力f_j，从而得到中心节点i的所受内力f_i，其值为所有邻节点所受力f_j的负总和：

(f)通过对下式进行数值积分计算求解血管内部节点i的形变位移向量：

式中，M_i、u_i、v_i、a_i分别表示节点i的质量、位移、速度和加速度，t表示迭代时间，f_ext、f_i分别表示节点i所受的外力和内力。

(2)使用基于位置动力学法对血管内部节点施加约束条件，对于由步骤(1)求出的血管内部节点的形变位置，将其投影到一个有效位置使其满足所定义的距离约束、体积守恒约束和弹性势能守恒约束三种约束条件，即找到某个修正因子Δx修正节点的形变位置以满足下式：

式中，C表示所定义的约束函数，节点的修正因子Δx_i可由下式表示：

其中，

表示节点x_i质量的倒数。

如图3所示，首先，定义任意两个节点之间的距离约束函数C_distance(x₁,x₂)：

C_distance(x₁,x₂)＝|x₁-x₂|-d₀

式中，d₀表示节点x₁与x₂间的初始距离；由距离约束得到节点位置的变化梯度

通过变化梯度求得各个节点的修正因子Δx_i(i＝1,2)：

然后，定义基于虚拟四面体单元的体积约束函数C_volume(x₁,x₂,x₃,x₄)以保证形变过后的软组织体积保持不变：

式中，V₀表示虚拟四面体单元(x₁,x₂,x₃,x₄)的初始体积；由体积约束得到节点位置的变化梯度

通过变化梯度求得各个节点的修正因子Δx_i(i＝1,2,3,4)：

如图4所示，最后，定义同样基于虚拟四面体单元的弹性势能守恒约束函数C_energy(x₁,x₂,x₃,x₄)，并假设四面体单元中任意两个节点间由虚拟弹簧连接，以此模拟血管的弹性特性：

式中，x₀表示四面体单元(x₁,x₂,x₃,x₄)的重心位置，k_i表示连接节点x_i，x₀的虚拟弹簧的弹性系数，d_i0表示节点x_i，x₀间的初始距离；由弹性势能守恒约束得到节点位置的变化梯度

通过变化梯度求得各个节点的修正因子Δx_i(i＝1,2,3,4)：

(3)使用无网格移动最小二乘法构建从血管内部到血管表面的映射模型来计算表面质点的形变位置，假设任意血管表面质点都可由其支持域内的一组内部节点来表示，如图5所示，假设某一表面质点为X，其支持域S内包含n个内部节点，定义从血管内部到血管表面的映射函数：

u(X)≈u^h(X)＝Φ(X)U_S

U_S＝[u₁,u₂,…,u_n]^T

式中，u(X)表示质点X处的场函数，u^h(X)表示u(X)的移动最小二乘近似函数，Φ(X)表示质点X支持域内的形函数，U_S为一个n维向量，用以描述支持域内所有节点的形变位移。

然后，利用无网格移动最小二乘法构造质点X形变位移的近似函数u^h(X)，故u^h(X)可近似表示为：

式中，p_j(X)为基函数，m为基函数的个数，a_j(X)为相应系数，其值为质点X所处空间坐标的函数，且

p^T(X)＝[p₁(X),p₂(X),…,p_m(X)]

a(X)＝[a₁(X),a₂(X),…,a_m(X)]^T

根据加权最小二乘法确定系数a(X)，使得u(X)的近似误差最小，定义Q：

式中，n为质点X支持域内的节点数，w(X-X_i)为权函数，X_i为支持域内的节点，u_i为节点X_i的形变位移；

上式可用矩阵形式表示为：

Q＝(Pa-U_s)^TW(X)(Pa-U_s)

式中：

为得到a(X)，对Q取极值，即得：

式中：A(X)＝P^TW(X)P，B(X)＝P^TW(X)。

因此，系数a(X)可表示为：

a(X)＝A^-1(X)B(X)U_s

故基于移动最小二乘法的近似函数u^h(X)为：

式中，形函数Φ(X)为：

Φ(X)＝[φ₁(X),φ₂(X),…,φ_n(X)]＝p^T(X)A^-1(X)B(X)

本发明采用如下基函数p(X)和权函数W(X)来构造形函数：

p^T(X)＝[1,x,y,z]

其中

为无量纲权函数影响半径，r_i为节点i的影响域半径。

本发明算法可以在不损害计算效率的同时保证了模拟血管形变的计算精度以及提高了形变模型的稳定性和实时性，让操作者在人机交互的过程中可以感受到软组织的弹性特性和体积守恒特性。

Claims (4)

1.一种虚拟血管介入手术训练的模拟方法，其特征在于，包括以下步骤：

(1)使用基于点的方法对血管内部构建形变模型；

(2)使用基于位置动力学法对血管内部节点施加距离约束、体积守恒约束和弹性势能守恒约束三种约束条件；

(3)使用无网格移动最小二乘法构建从血管内部到血管表面的映射模型；

步骤(1)中，所述使用基于点的方法对血管内部构建形变模型，是指根据血管的医学影像数据，使用基于点的方法将血管内部离散化为一系列基于四面体的点云，并将血管的体积均匀分配给节点，通过获得节点的应力和应变来计算节点的形变位移向量；

所述通过获得节点的应力和应变来计算节点的形变位移向量，包括以下步骤：

(a)使用多项式核函数ω_ij来衡量中心节点i对其邻节点j的影响能力，具体计算公式为：

式中，h表示节点i的支撑半径，r表示节点i与其邻节点j之间的距离；

(b)计算节点i位移向量u_i的空间导数▽u_i，具体计算公式为：

式中，

和/>

分别表示节点i的位移场u＝(m,n,p)^T中横坐标所对应的值m的空间导数、纵坐标所对应的值n的空间导数以及竖坐标所对应的值p的空间导数；

(c)通过求得的

可得到节点i处的应变ε_i和应力σ_i，具体计算公式为：

σ_i＝Cε_i

式中，J_i表示节点i的雅可比矩阵，I表示单位矩阵，C表示弹性矩阵，其值取决于弹性材料的杨氏模量和泊松比；

(d)根据连续介质力学理论，估计出节点i周围所存储的应变能U_i具体计算公式为：

式中，Ω表示节点i的支撑域，ν_i表示节点i的体积；

(e)通过对邻节点位移向量求导计算出每个邻节点的所受力f_j，从而得到中心节点i的所受内力f_i，其值为所有邻节点所受力f_j的负总和，具体计算公式为：

式中，

表示针对邻节点j的位移向量u_j求导；

(f)通过对下式进行数值积分计算求解血管内部节点i的形变位移向量：

式中，M_i、u_i、v_i、a_i分别表示节点i的质量、位移、速度和加速度，t表示迭代时间，f_ext、f_i分别表示节点i所受的外力和内力；

步骤(2)中，所述距离约束C_distance(x₁,x₂)的函数公式为：

C_distance(x₁,x₂)＝|x₁-x₂|-d₀

式中，d₀表示节点x₁,x₂间的初始距离，根据距离约束的约束条件得到的节点修正因子Δx_i计算公式为：

式中，

表示节点x_i质量的倒数，其中，i＝1,2；

步骤(2)中，所述体积守恒约束C_volume(x₁,x₂,x₃,x₄)的函数公式为：

式中，V₀表示虚拟四面体单元(x₁,x₂,x₃,x₄)的初始体积，根据体积守恒约束的约束条件得到的节点修正因子Δx_i计算公式为：

式中，

表示节点x_i质量的倒数，其中，i＝1,2,3,4；/>

表示针对节点x_i、x_j求导；

步骤(2)中，所述弹性势能守恒约束C_energy(x₁,x₂,x₃,x₄)的函数公式为：

式中，x₀表示四面体单元(x₁,x₂,x₃,x₄)的重心位置，k_i表示连接节点x_i，x₀的虚拟弹簧的弹性系数，d_i0表示节点x_i，x₀间的初始距离，根据弹性势能守恒约束的约束条件得到的节点修正因子Δx_i计算公式为：

式中，

表示节点x_i质量的倒数，其中，i＝1,2,3,4；

步骤(3)中，所述无网格移动最小二乘法中定义的映射函数为：

式中，u^h(X)为质点X的场函数u(X)的近似函数，Φ(X)为质点X支持域内的形函数，U_S为一个n维向量，用以描述支持域内所有节点的形变位移，

表示第i点的形函数，u_i为节点i的位移向量。

2.根据权利要求1所述的虚拟血管介入手术训练的模拟方法，其特征在于，步骤(b)中，所述节点i的位移场u＝(m,n,p)^T中横坐标所对应的值m的空间导数

的计算公式为：

式中，x_ij＝x_i-x_j，

表示矩量矩阵A_i的逆，m_i,m_j分别表示节点i和j的位移场u＝(m,n,p)^T横坐标所对应的值，ω_ij表示节点i和j间的权重，由多项式核函数计算，x_ij表示节点i的位置x_i与节点j的位置x_j之间的距离位移向量，||x_ij||表示节点i与节点j之间的距离。

3.根据权利要求2所述的虚拟血管介入手术训练的模拟方法，其特征在于，步骤(b)中，节点i的位移场u＝(m,n,p)^T中纵坐标所对应的值n的空间导数

和节点i的位移场u＝(m,n,p)^T中竖坐标所对应的值p的空间导数/>

的计算方法与节点i的位移场u＝(m,n,p)^T中横坐标所对应的值m的空间导数/>

的计算方法相同。

4.根据权利要求1所述的虚拟血管介入手术训练的模拟方法，其特征在于：

采用基函数p(X)和权函数W(X)来构造形函数，基函数p(X)、权函数W(X)的函数公式分别为：

p^T(X)＝[1,x,y,z]

式中，

为无量纲权函数影响半径，r_i为节点i的影响域半径。

Images