一种用于模拟软组织形变和路径切割的方法及装置

1.本发明涉及一种软组织力触觉模拟方法，尤其涉及一种用于模拟软组织形变和路径切割的方法及装置。


背景技术：

2.近年来远程医疗会诊的需求日益凸显，远程医疗会诊中的一个关键技术难点就是患者虚拟病变器官的精确模拟，从而促使会诊专家提出正确的、科学的、合适的治疗方案，最终实现精确指导。在现有的虚拟器官模拟方法中，有限元模型将软组织器官离散为如四面体、六面体等有限个单元体，单元体之间通过网格节点相连，单元体上任一点的位移是用各个单元节点的位移作为变量的函数来表示，通过求解节点位移，得到软组织形变量，最终模拟软组织器官的形变行为。它依据弹性力学，将复杂的求解域离散化求解，因此具有仿真精度高的特点。但由于该模型在形变过程中网格单元的拓扑结构会不断重组，导致计算量较大、计算效率较低，极大地影响仿真实时性，因而在远程会诊中实施困难。


技术实现要素：

3.发明目的：针对有限元模型存在的问题，本发明提出一种既能降低形变计算量又能保证仿真实时性，并当形变达到极限时，形变模型从该处开始断裂，由此形成一种用于模拟软组织形变和路径切割的方法及装置。
4.技术方案：本发明提供一种用于模拟软组织形变和路径切割的方法，具体包括以下步骤：
5.(1)使用有限元模型对待模拟软组织进行形变模拟，并在形变仿真过程中应用本征正交分解与伽辽金投影相结合的模型降阶法实现模型降阶；
6.(2)根据手术机械与软组织的碰撞检测形成的交点位置形成切割路径；
7.(3)使用贝塞尔曲线来绘制软组织被切割后产生的表面切口。
8.进一步地，步骤(1)所述的使用有限元模型对待模拟软组织进行形变模拟的过程如下：
[0009][0010]
式中，m表示质量矩阵，u表示位置向量，v表示速度向量，g表示软组织所受外力，t表示迭代时间，f表示软组织所受内力，h
t
λ表示约束力，由手术器械作用于软组织表面时产生，h表示时间间隔矩阵，λ表示约束参数，
′
表示转置；
[0011]
将软组织离散化为一系列四面体网格单元，使用隐式欧拉法对有限元模型的动态表达式进行数值计算；考虑将连续的迭代时间离散为区间在时间间隔[t
n
，t
n+1
]上的具体数值计算公式为：
[0012][0013]
式中，令h表示时间间隔，即h＝t
n+1
‑
t
n
，表示求导符号，d表示微分符号，表示f在t
n
时刻的值，表示g在t
n+1
时刻的值。
[0014]
进一步地，步骤(1)所述的模型降阶法实现过程如下：
[0015]
利用本征正交分解对有限元全阶模型进行降阶，从而求得一组正交基函数φ＝(φ1，φ2，...，φ
n
)，使之能够最大程度地近似原始样本，获得正交基函数的具体计算公式如下：
[0016][0017]
式中，j表示最小二乘意义下全阶模型的解向量样本与基函数之间的误差函数，为参数空间λ的一个离散子集，表示中的某个分量，t0和分别表示迭代时间t的两个离散时刻，表示存储在快照矩阵的数据，φ
i
表示正交基函数的第i个分量；
[0018]
通过得到的一组本征正交分解基函数，并将这一组基函数张成的空间作为降阶模型解向量所在的函数空间，然后使用伽辽金投影将全阶模型投影至降阶空间；使用伽辽金投影求解降阶模型解向量的具体计算公式如下：
[0019][0020]
式中，α表示正交基函数φ的系数向量，
[0021]
进一步地，所述步骤(2)包括以下步骤：
[0022]
(21)将手术器械简化抽象为一条线段，手术器械与软组织的碰撞检测即为线段与三角单元的交点检测；假设线段的两端点分别为i和l，它们的位置向量分别为u
i
(x
i
，y
i
，z
i
)和u
l
(x
l
，y
l
，z
l
)，则手术器械的空间直线方程具体表示如下：
[0023][0024]
式中，x、y和z分别表示空间直线方程的三个未知数，k表示空间直线方程的参数；则线段上的任一节点坐标为：
[0025]
[0026]
(22)确定与手术器械相交的三角单元的平面方程：假设碰撞的三角单元的三个端点分别为o、p、q，它们的位置坐标分别为u
o
(x
o
，y
o
，z
o
)、u
p
(x
p
，y
p
，z
p
)和u
q
(x
q
，y
q
，z
q
)，且三角单元的法向量为n(n
x
，n
y
，n
z
)，故三角单元的平面方程为：
[0027]
ux+vy+wz+t＝0
[0028]
式中，u、v、w和t表示平面方程的四个待求参数，x、y和z表示平面方程的未知数；
[0029]
根据计算平面方程的点法式方法，通过点o和法向量n的坐标可得平面方程：
[0030]
n
x
(x
‑
x
o
)+n
y
(y
‑
y
o
)+n
z
(z
‑
z
o
)＝0
[0031]
式中，u＝n
x
，v＝n
y
，w＝n
z
，t＝
‑
n
x
x
o
‑
n
y
x
o
‑
n
z
z
o
；
[0032]
(23)联立直线方程和平面方程，从而获得任一离散时刻手术器械与三角单元平面的交点s的位置u
s
(x
s
，y
s
，z
s
)，其具体计算公式如下：
[0033][0034]
进一步地，步骤(3)所述的贝塞尔曲线是通过三个控制点并采用2次bemstein基函数构成；所述2次bernstein基函数r
i，2
(t)的具体表达式为：
[0035][0036]
进一步地，步骤(3)所述的表面切口方程为：
[0037][0038]
式中，e(t)表示绘制切割路径的控制点在迭代时间t时的位置向量，r
i，2
(t)表示2次bemstein基函数。
[0039]
基于相同的发明构思，本发明还提供一种用于模拟软组织形变和路径切割的装置，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述计算机程序被加载至处理器时实现上述的用于模拟软组织形变和路径切割的方法。
[0040]
有益效果：与现有技术相比，本发明的有益效果：(1)使用基于模型降阶法的有限元模型构建软组织形变模型，利用本征正交分解获得一组正交基函数，然后通过伽辽金投影将有限元全阶模型的解向量投影至这一组正交基函数所张成的降阶空间，以此获得软组织节点的形变位移向量，从而减少了形变计算量，提高了计算效率并保证了仿真实时性；(2)采用直线模型来模拟切割工具，将手术器械简化为一条直线。把手术器械与软组织的碰撞转化为该直线与虚拟软组织表面三角单元的相交检测，从而实现模拟切割过程，既保证了仿真实时性，同时降低了切割模拟的复杂程度；(3)采用贝塞尔曲线来绘制软组织被切割
后产生的表面切口，通过规定曲线端点和增加曲线中间的控制点来控制曲线首末以及确定曲线曲度，大大降低切割程序的空间复杂度，并且能够促使切口效果更加平滑、自然。
附图说明
[0041]
图1是本发明的流程图；
[0042]
图2是本发明所述的降阶模型示意图；
[0043]
图3是本发明所述的手术器械与三角单元的交点示意图；
[0044]
图4是本发明所述的贝塞尔曲线切口绘制示意图；
[0045]
图5是基于不同方法切割胃部软组织的切口效果图，其中，(a)为基于本发明的模拟效果图；(b)为基于b样条曲线的模拟效果图；(c)为基于非均匀有理b样条曲线的模拟效果图。
具体实施方式
[0046]
下面结合附图对本发明作进一步详细说明。
[0047]
本发明提供一种用于模拟软组织形变和路径切割的方法，如图1所示，具体包括以下步骤：
[0048]
步骤1：使用有限元模型对待模拟软组织进行形变模拟，并在形变仿真过程中应用本征正交分解与伽辽金投影相结合的模型降阶方法实现模型降阶。
[0049]
根据医学影像数据，使用基于模型降阶法的有限元模型对人体器官构建形变模型，将一组有限元全阶模型的求解结果作为样本，并在此基础上通过本征正交分解和伽辽金投影相结合的模型降阶方法，将全阶模型降阶为降阶模型，然后把这一组解向量投影至降阶模型所在的降阶空间，来计算网格单元节点的形变位移向量。使用opengl将ct扫描采集的胃部医学图像数据进行三维几何重现，将胃部软组织离散化为一系列四面体网格单元，然后使用降阶模型进行形变建模，如图2所示。
[0050]
通过将有限元全阶模型的解向量投影至降阶空间中来计算节点的形变位移向量，具体包括以下步骤：
[0051]
(a)使用有限元模型来模拟软组织形变，其动态表达式为：
[0052][0053]
式中，m表示质量矩阵，u表示位置向量，v表示速度向量，g表示软组织所受外力，t表示迭代时间，f表示软组织所受内力，h
t
λ表示约束力，由手术器械作用于软组织表面时产生，h表示时间间隔矩阵，λ表示约束参数，
′
表示转置。
[0054]
(b)考虑将连续的迭代时间t离散为区间在时间间隔[t
n
，t
n+1
]上对有限元模型动态表达式进行积分：
[0055][0056][0057]
式中，h表示时间间隔，即h＝t
n+1
‑
t
n
，表示求导符号，d表示微分符号，表示f在t
n
时刻的值，表示g在t
n+1
时刻的值。
[0058]
(c)将软组织离散化为一系列四面体网格单元，使用隐式欧拉法对有限元模型的动态表达式进行数值计算，可以得到下式：
[0059][0060][0061]
式中，表示所构成的矩阵。
[0062]
(d)由于软组织所受内力f(u，t)是关于位置和速度的非线性函数，因此使用newton
‑
raphson迭代对其进行迭代计算：
[0063][0064]
式中，中，表示软组织所受内力f在t
n
时刻的值。
[0065]
(f)最后，得到一组有限元模型全阶模型的求解结果，具体数值计算公式为：
[0066][0067]
式中，令
[0068]
(g)将通过步骤(f)得到的一组解向量组成快照矩阵并对u(t，λ)进行本征正交分解：
[0069][0070]
式中，φ＝(φ1，φ2，...，φ
n
)表示表示基于本征正交分解的正交基函数，α
i
表示基函数φ
i
的对应系数。
[0071]
(h)对求得的正交基函数进行截断，以在最大程度降低全阶模型维度的情况下近似原始样本。定义误差函数j，来最终选取表征全阶模型解向量的正交基函数φ＝(φ1，φ2，...，φ
n
)：
[0072][0073]
式中，j表示最小二乘意义下全阶模型的解向量样本与基函数之间的误差函数，为参数空间λ的一个离散子集，λ
*
表示中的某个分量，t0和分别表示迭代时间t的两个离散时刻，表示存储在快照矩阵
的数据，φ
i
表示正交基函数的第i个分量。
[0074]
(i)通过步骤(h)得到一组本征正交分解基函数，并将这一组基函数张成的空间作为降阶模型解向量所在的函数空间，然后使用伽辽金投影将全阶模型投影至降阶空间。使用伽辽金投影求解降阶模型解向量的具体计算公式如下：
[0075][0076]
式中，
[0077]
步骤2：根据手术机械与软组织的碰撞检测形成的交点位置形成切割路径。
[0078]
切割路径可以理解为手术器械在离散时间内在软组织表面划过的一段连续的折线。因此，首先将手术器械简化抽象为一条线段，手术器械与软组织的碰撞检测即为线段与三角单元的交点检测，如图3所示。首先，假设线段的两端点分别为i和l，它们的位置向量分别为u
i
(x
i
，y
i
，z
i
)和u
l
(x
l
，y
l
，z
l
)，则手术器械的空间直线方程具体表示如下：
[0079][0080]
式中，x、y和z分别表示空间直线方程的三个未知数，k表示空间直线方程的参数；则线段上的任一节点坐标为：
[0081][0082]
其次，确定与手术器械相交的三角单元的平面方程。假设碰撞的三角单元的三个端点分别为o、p、q，它们的位置坐标分别为u
o
(x
o
，y
o
，z
o
)、u
p
(x
p
，y
p
，z
p
)和u
q
(x
q
，y
q
，z
q
)，且三角单元的法向量为n(n
x
，n
y
，n
z
)，故三角单元的平面方程为：
[0083]
ux+vy+wz+t＝0
[0084]
式中，u、v、w和t表示平面方程的四个待求参数，x、y和z表示平面方程的未知数。
[0085]
然后，根据计算平面方程的点法式方法，通过点o和法向量n的坐标可得平面方程：
[0086]
n
x
(x
‑
x
o
)+n
y
(y
‑
y
o
)+n
z
(z
‑
z
o
)＝0
[0087]
式中，u＝n
x
，v＝n
y
，w＝n
z
，t＝
‑
n
x
x
o
‑
n
y
x
o
‑
n
z
z
o
。
[0088]
最后，联立直线方程和平面方程，从而获得任一离散时刻手术器械与三角单元平面的交点s的位置u
s
(x
s
，y
s
，z
s
)，其具体计算公式如下：
[0089][0090]
步骤3：使用贝塞尔曲线来绘制软组织被切割后产生的表面切口。
[0091]
首先，确定切割路径的开始位置a和终点位置b，并获取切割路径经过的所有表面节点，记为e1＝{e
10
，e
11
，e
12
，e
13
，e
14
，e
15
}，切割路径周边的节点记为点集e2＝{e
20
，e
21
，e
22
，
e
23
，e
24
，e
25
，e
26
，e
27
，e
28
，e
29
，e
30
，e
31
}，分布于切割路径两侧。选取节点c和d为规定表面切口曲度的控制点，如图4所示。然后根据a、b、c三个顶点和a、b、d三个顶点分别绘制两条二次贝塞尔曲线来表示由切割操作产生的切口。
[0092]
其中，假设n+1个控制点的位置为e
i
＝(x
i
，y
i
，z
i
)，i＝0，1，
…
，n，由这些控制点混合产生位置向量e(t)，用以描述e0与e
n
之间逼近贝塞尔多项式函数的路径：
[0093][0094]
式中，n表示多项式次数，通常由切口曲线的n+1个控制点确定，r
i，n
(t)为多项表达式，叫做n次bernstein基函数，定义为：
[0095][0096]
本文选取二次贝塞尔曲线实现切口的再现，其中二次贝塞尔曲线通过3个控制点产生，把n＝2代入bernstein基函数中，由此得到2次bernstein基函数集合：
[0097][0098]
然后，基于二次贝塞尔曲线得到绘制表面切口的方程为：
[0099][0100]
最后，分别采用本发明、b样条曲线以及非均匀有理b样条曲线对虚拟胃部软组织进行切割模拟对比。图5给出了基于不同方法切割胃部软组织的切口效果图。从图5中可以看出，基于本发明的模拟切割的逼真度较高，切口更加真实、光滑和自然。
[0101]
基于相同的发明构思，本发明还提供一种用于模拟软组织形变和路径切割的装置，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述计算机程序被加载至处理器时实现上述的用于模拟软组织形变和路径切割的方法。
[0102]
本发明可以减少形变计算量的同时保证计算效率和模拟胃部软组织形变的计算精度，以及提高仿真实时性，让操作者在人机交互的过程中可以感受到软组织逼真的切割过程。