一种基于无网格径向基数据拟合的软组织受力形变模型建模方法

摘要

一种基于无网格径向基数据拟合的软组织受力形变模型建模方法。针对现有生物软组织物理形变模型结构复杂、实时性不佳等缺陷，考虑生物软组织各向异性特征，以力为数据流，由径向基无网格法离线获得以软组织受力点为中心，面上有限个等间距点受力位移数据，使用麦夸特法构建软组织表面节点力‑位移曲面形变模型。该模型结合几何形变模型及物理形变模型特点，在保留形变模型物理性特性真实性的同时，大大简化了计算过程，保证了交互的高实时性，实现虚拟手术中生物软组织形变模型与操作者准确实时交互。

说明书

技术领域

本发明涉及虚拟手术软组织受力形变模型建模，用于实现在虚拟手术中生物软组织模型与操作者的精确实时交互。

背景技术

当进行虚拟手术操作时，理想的生物软组织模型能够在实时交互时将其形变情况准确即时地反馈给操作者，这就要求生物软组织模型必须具备高精确性实时形变特性。

软组织的形变模型按照软组织形变过程是否基于物理原理，可以分为几何形变模型和物理形变模型。最早的几何形变模型为Sederberg提出的FFD(Free Form Deformation)模型，该模型完全从几何角度对物体的形变进行仿真。此外，Gibson提出了Chain-Mail模型，用于实现对膝关节组织形变实时仿真，，通过限定结点间的距离极限，来模拟位移传递过程。大多数医学仿真是针对某一具体器官而言。人体器官的几何形态建模是各种虚拟手术仿真的基础，对于可以交互操作的物理仿真更是如此。几何形态一般是用三角面表示的面模型以及四面体等体单元表示体模型，或由二者组合而成。因此，虚拟器官的形态建模就是设法建立上述面模型或体模型。有些模型可以采用纯数学的方法如样条曲面建模，有些可以借助统计数据建立，如采用拟合多项式建立表面模型，但几乎都是通过医学影像数据或组织切片彩色图像进行三维表面重构。早期，由于计算机运算能力的限制，研究者为了获得比较高的运行效率，不得不基于纯粹的几何形变技术对形变对象的形变过程进行建模。尽管集合形变技术能够提高了模型的运行速度，但是形变体在形变过程中所表现出来的物理性质被忽略，例如，没有考虑形变体的材料性质，仅仅从外观形式对形变过程进行描述，缺乏对力反馈的物理支持，从而影响到虚拟软组织器官模型的反馈力与形变关系之间关联的准确性。

物理形变模型主要包括质点弹簧建模MSM(Mass Spring Modeling)、光滑粒子流模型SPH(Smoothed Paticle Hydrodynamics)、有限元方法FEM(Finite Elements Method)、长元素方法LEM(Long Elements Method)和无网格法(Mesh Free)。基于物理的形变模型建模技术涉及到牛顿动力学、连续介质机械理论、差分几何、矢量运算、逼近理论和数值计算等领域的一个新的综合应用领域，这个领域仍然在不断发展和创新中。基于物理的形变模型通过应变能量函数定量地描述用户所需形变结果的特性，能使得形变结果贴近用户所描述的各种几何、物理属性。应变能量函数可以基于弹性力学方法定义，计算上具有一定程度的复杂性，但是逼近实际软组织的形变效果比较好。它也是计算机图形学中一个新兴的研究领域，研究对象包括刚体、弹塑性和流体等类型材料的运动、形变、流动和破裂等现象的仿真，研究方法包括时间积分策略、方程求解方法、模态分析、多分辨建模型方法和实时计算技术。非线性的物理模型具有非常高的计算复杂性，给计算机实时仿真带来了挑战。

综上所述，几何形变模型虽然运行效率高，但大多不考虑软组织的实际力学本构方程，也不考虑形变过程中物体质量、力或其他物理现象的作用，而只考虑几何形态的变化牺牲了软组织的物理性质。物理形变模型则基于软组织的力学本构方程，通过相应的计算模型得出组织受力时的形变。相对于几何形变模型，物理模型能够更加真实的反映软组织在虚拟手术过程中的形变过程，但具有计算量大，实时仿真速度低的缺点。

发明内容

本发明公开了一种基于曲面拟合的径向基无网格软组织数据的力反馈模型建模方法。针对现有生物软组织物理形变模型结构复杂、实时计算性能不佳等缺陷，考虑生物软组织各向异性特征，以力为数据流，通过径向基无网格法离线获得软组织受力点为中心，面上有限个等间距点受力位移数据，使用麦夸特法构建软组织表面节点力-位移曲面形变混合模型。该模型综合了几何形变模型及物理形变模型的优点，在保留形变模型物理真实性的同时，大大简化计算过程，保证了交互的实时性，实现虚拟手术中操作者与生物软组织形变模型交互时的实时性和准确性。

本发明是通过以下技术方案实现的。

生物软组织受空间力F作用后表面形变有如下曲面系统：

其中v_o为软组织未受力时表面各节点基准坐标，Δv_x、Δv_y、Δv_z分别为空间力F沿空间直角坐标系分量F_x、F_y、F_z作用下的软组织表面曲面形变分量，为随机误差项。

步骤(1)围绕受力点O建立软组织表面m个等间距d奇数节点矩阵，取软组织表面基准坐标v_o。

步骤(2)以垂直水平面Z轴方向施力F_z，由径向基无网格法获得力与软组织表面各节点形变数据，采样当前输入输出数据。获取L组软组织数据：F_zi→(Δx_i,Δy_i,Δz_i)_j，i∈{1,2,…,L}，j∈{1,2,…,m}。

步骤(3)麦夸特算法拟合力F_z与Δz_n函数关系。选取中间节点即受力点偏移量：ΔO＝(Δx_i,Δy_i,Δz_i)_n，n＝(1+m)/2，i∈{1,2,…,L}。获得L组软组织力与形变数据，由于施力方向垂直水平面受力点水平方向，Δx_ni、Δy_ni偏移为0，只存在Z轴方向偏移Δz_ni。因此基于麦夸特算法建立数据流Δz_n与力F_z的位移与力关系：Δz_n＝g_z(F_z)。

步骤(4)通过步骤(3)获得L组中间节点形变Δz_ni对应软组织表面曲面各节点偏移量(Δx_i，Δy_i，Δz_i)_j数据：Δz_ni→(Δx_i,Δy_i,Δz_i)_j，i∈{1,2,…,L}，j∈{1,2,…,m}，利用麦夸特算法拟合当受力节点偏移量为Δz_ni时对应曲面形变函数：Δv_zi＝γ_zi(x,y；p_i),i∈{1,2,…,L}获得受力节点偏移量为Δz_ni时对应曲面形变函数参数p_i：Δz_ni→p_i,i∈{1,2,…,L}。

步骤(5)由L组受力节点偏移量Δz_n对应曲面形变函数参数p数据，基于麦夸特算法拟合受力节点偏移量Δz_n与曲面形变参数p的对应关系：p＝ρ_z(Δz_n)。

步骤(6)以平行X轴方向向O点施力F_x，由径向基无网格法获得力与软组织表面各节点形变数据，采样当前输入输出数据。获取L组软组织数据。同理步骤(3)、(4)、(5)，可获得O节点施力F_x下：Δx_n＝g_x(F_x)；Δv_x＝γ_x(x,y；p)，p＝ρ_x(Δx_n)。

步骤(7)以平行Y轴方向O点施力F_y，由径向基无网格法获得力与软组织表面各节点形变数据，采样当前输入输出数据。获取L组软组织数据。同理步骤(3)、(4)、(5)，可获得O节点施力F_y下：Δy＝g_y(F_y)；Δv_y＝γ_y(x,y；p)，p＝ρ_y(Δy_n)。

步骤(8)F_x、F_y、F_z合成空间力F，可得当空间力F作用于O点时有软组织表面形变Δv_x、Δv_y、Δv_z，则有软组织形变函数：

步骤(9)力反馈对步骤(2)、步骤(6)和步骤(8)获得的受力节点O位移与施力F函数：Δx_n＝g_x(F_x)、Δy＝g_y(F_y)、Δz_n＝g_z(F_z)取反可得各实力分量与受力节点O位移函数：

从而可得力反馈

本发明的优点：以曲面拟合方法快速拟合由径向基无网格法获得的软组织形变输入输出数据，建立生物软组织混合形变模型。该方法结构简洁、计算方便，综合了物理形变模型及几何形变模型的优点，在保留径向基无网格法物理真实性的同时，又大大简化计算量，保证了交互的高实时性，实现虚拟手术中生物软组织形变模型与操作者的实时准确交互。

附图说明

图1为软组织表面节点受力位移图。

图2为曲面拟合图。

图3为分别对O点施加F_x、F_y、F_z时软组织形变示意图。其中，(a)为对O点施加F_x时软组织形变示意图，(b)为对O点施加F_y时软组织形变示意图，(c)为对O点施加F_z时软组织形变示意图。

具体实施方式

本发明将通过以下实例作进一步说明。

生物软组织受空间力F作用后表面形变有如下曲面系统：

其中，v_o为软组织未受力时表面各节点基准坐标，Δv_x、Δv_y、Δv_z分别为空间力F沿空间直角坐标系分量F_x、F_y、F_z作用下的软组织表面曲面形变分量，为随机误差项。

步骤1围绕受力点O建立软组织表面m个等间距为d的奇数节点矩阵。取软组织表面基准坐标；在现实手术中生物组织受力形变的组织直径范围在1～4cm左右，选取48mm×48mm大小软组织进行离散取点，取121个等间距节点，间距d＝4.8mm，选取121个点的中间点作为受力O点，此时，受力点O为第61节点，即n＝61。

步骤2以垂直水平面Z轴方向施力F_z，0<F_z≤5N，0～5N的施力范围符合真实手术中对软组织所施力的范围，由径向基无网格法获得力与软组织表面各节点形变数据，以0.05N为步长，采样当前输入输出数据。获取L＝100组软组织数据：

F_zi→(Δx_i,Δy_i,Δz_i)_j，i∈{1,2,…,L}，j∈{1,2,…,m}。(1)

步骤3麦夸特算法拟合力F_z与受力点偏移量Δz_n函数关系。选取中间节点即受力点偏移量：ΔO＝(Δx_i,Δy_i,Δz_i)_n，i∈{1,2,…,L}。获得L组软组织力与受力点偏移数据，由于施力方向垂直水平面受力点水平方向Δx_ni、Δy_ni偏移为0，只存在Z轴方向偏移Δz_ni。由此基于麦夸特算法建立数据流Δz_n与力F_z的位移与力关系：Δz_n＝g_z(F_z)。拟合过程具体描述如下：

(a)对于L组力与受力点偏移数据：F_zi→Δz_ni，i∈{1,2,…,L}，设函数形式为：

其中a为系统自变量最高阶数，α＝[α₁,α₂,…α_h](h＝a+1)为系统参数，为系统误差项。

(b)为避免过度拟合，利用交叉验证法确定(1)式系统自变量阶数a，而确定了自变量阶数的同时也确定了参数α的个数。对F_z做L次拟合，将第i次观测值剔除，i∈{1,2,…,L}，计算残差平方和，寻找到使其最小的阶数A

其中b_-i,a为删掉第i个观测之后的估计值。

(c)初步描述软组织力与受力点偏移函数系统为：

Δz_n＝α₁+α₂F_z¹+…+α_a+1F_z^a

(d)设置待估参数α初值设定阻尼因子d＝10^ud⁽⁰⁾(d⁽⁰⁾>0)；设定迭代次数u＝0即d＝d⁽⁰⁾，麦夸特法由于阻尼因子的引入，在一般的初值选择条件下都可收敛于所求结果。

(e)设定迭代次数u＝0，进行第一次迭代，利用L组软组织数据、参数初值α⁽⁰⁾及阻尼因子初值d⁽⁰⁾通过麦夸特算法参数估计解得α值，将求得参数α代入系统原函数式，计算残差平方和J⁽⁰⁾。

(f)将g_z(F_zi,α)在α⁽⁰⁾处按泰勒级数展开并略去二次及二次以上的项得：

(g)式(2)中除待估参数(α₁,α₂,…α_n)外，都是已知。对此使用最小二乘残差平方原理：

其中d≥0为阻尼因子。

(h)通过将残差平方J分别对α₁,α₂,…α_h一阶偏导设为0，以达到将其最小化的目的：

(i)令

可将式(6)写成：

(j)解得

(k)将步骤(j)计算获得参数α赋予α⁽⁰⁾，即α⁽⁰⁾＝α，令u＝1即d＝10d⁽⁰⁾进行第二次迭代，回到(e)继续执行解得新的参数α并代入系统原函数式，计算新的残差平方和J⁽¹⁾。比较J⁽¹⁾与J⁽⁰⁾，增加迭代次数，直至J⁽¹⁾<J⁽⁰⁾为止。具体过程描述如下：

比较：若J⁽¹⁾<J⁽⁰⁾，(k)结束；若J⁽¹⁾≥J⁽⁰⁾，取u＝2即d＝100d⁽⁰⁾，回到(e)继续执行重解α并计算新残差平方J⁽¹⁾，若新计算的J⁽¹⁾<J⁽⁰⁾(k)结束，反之则再取u＝3即d＝1000d⁽⁰⁾回到(e)继续执行…，如此反复迭代直至J⁽¹⁾<J⁽⁰⁾为止。在式(6)中，a_ij为定值，所以当d越大时Φ＝|α-α⁽⁰⁾|越小，残差平方和J收敛，而d相对初值d⁽⁰⁾增加的越大即u越大迭代的次数越多。选择的界限是看残差平方和是否下降。

(l)继续增加阻尼因子d的数量级即u＝u+1，反复执行(k)，直至新计算获得参数α与上一次计算参数即α⁽⁰⁾之差绝对值|Φ|≤δ允许最小误差为止，使得估计参数波动范围忽略不计。此时拟合力F_z与Δz_n函数完成，获得关系式：

Δz_n＝α₁+α₂F_z¹+…+α_a+1F_z^a(10)

步骤4由式(1)和(10)可获得L组中间节点形变Δz_ni对应软组织表面曲面各节点偏移量(Δx_i,Δy_i,Δz_i)_j数据，即：Δz_ni→(Δx_i,Δy_i,Δz_i)_j，i∈{1,2,…,L},j∈{1,2,…,m}。利用麦夸特算法拟合当受力节点偏移量为Δz_ni时对应曲面形变函数：Δv_zi＝γ_zi(x,y；p)，i∈{1,2,…,L}，获得受力节点偏移量为Δz_ni时对应曲面形变函数参数p_i：Δz_ni→p_i，i∈{1,2,…,L}。拟合曲面形变函数过程具体描述如下：

(aa)空间曲面形变包含有纵向形变分量ΔY，横向形变分量ΔX，垂直方向形变分量ΔZ，因此当受力节点偏移量为Δz_ni时曲面形变函数形式为：

(bb)对于m个奇数节点软组织表面基准坐标v_o：(x_o,y_o,z_o)_j，j∈{1,2,…,m}，有m组曲面形变纵向形变分量ΔY、横向形变分量ΔX、垂直方向形变分量ΔZ与软组织表面基准水平坐标(x_o,y_o)_j，j∈{1,2,…,m}对应数据：

设受力节点偏移量为Δz_ni时曲面各分量函数形式为：

其中h₁＝a₁+b₁+1、h₂＝a₂+b₂+1、h₃＝a₃+b₃+1。

(cc)参照步骤3(b)利用交叉验证法可求得曲面形变阶数a₁、a₂、a₃、b₁、b₂、b₃。

(dd)求得阶数后初步描述受力节点偏移量为Δz_ni时曲面形变函数为：

由基准坐标(x_o,y_o)_j，j∈{1,2,…,m}，0≤x_o≤4cm,0≤y_o≤4cm，令(x,y)＝(x_o,y_o)_j，x,y∈(0,4cm)，式(11)可写成：

(ee)参照步骤3函数参数麦夸特算法求解过程(d)～(j)，可解得受力节点偏移量为Δz_ni时曲面形变函数参数为：

其中：为按照式(4)在p_i1⁽⁰⁾处按泰勒级数展开并略去二次及二次以上的项，

其中：为按照式(4)在p_i2⁽⁰⁾处按泰勒级数展开并略去二次及二次以上的项，

其中：为按照式(4)在p_i3⁽⁰⁾处按泰勒级数展开并略去二次及二次以上的项，

(ff)参照执行步骤3的(k)、(l)过程获得有效曲面参数，即当中间节点形变Δz_ni时对应曲面形变函数参数。此时曲面函数拟合完成，当中间节点形变Δz_ni时曲面函数为：

因此对于L组中间节点形变Δz_ni对应软组织表面曲面各节点偏移量数据，即：Δz_ni→(Δx_i,Δy_i,Δz_i)_j，i∈{1,2,…,L},j∈{1,2,…,m}连续求解参数可以获得L组中间节点形变Δz_ni对应软组织表面曲面形变函数参数数据：

Δz_ni→[p_i1,p_i2,p_i3]，i∈{1,2,…,L}

其中

步骤5由L组受力节点偏移量Δz_n对应曲面形变函数参数p数据，基于麦夸特算法拟合受力节点偏移量Δz_n与曲面形变参数p的对应关系：p＝ρ_z(Δz_n)。具体拟合过程描述如下：

(aaa)对于L组受力节点偏移量Δz_n对应曲面形变函数参数p数据Δz_ni→[p_i1,p_i2,p_i3]，设函数形式为：

其中

参照步骤4可解得β得到如下式：

步骤6以平行X轴方向向O点施力F_x，由径向基无网格法获得力与软组织表面各节点形变数据，采样当前输入输出数据。获取L组软组织数据。同理步骤3、4、5，可获得O节点施力F_x下：Δx_n＝g_x(F_x)，p＝ρ_x(Δx_n)，Δv_x＝γ_x(x,y；p)。

步骤7以平行Y轴方向向O点施力F_y，由径向基无网格法获得力与软组织表面各节点形变数据，采样当前输入输出数据。获取L组软组织数据。同理步骤3、4、5，可获得O节点施力F_y下：Δy_n＝g_y(F_y)，p＝ρ_y(Δy_n)，Δv_y＝γ_y(x,y；p)。

步骤8F_z、F_y、F_z合成空间力F，可得当空间力F作用于O点时有软组织表面形变Δv_x、Δv_y、Δv_z，则有软组织形变函数：

步骤9力反馈对步骤2、步骤6和步骤8获得的受力节点O位移与施力F函数：Δx_n＝g_x(F_x)、Δy＝g_y(F_y)、Δz_n＝g_z(F_z)取反可得各实力分量与受力节点O位移函数：

从而可得力反馈