CN106485000A - 基于虚拟弹簧和球面调和函数力反馈模型建模方法 - Google Patents

Abstract

本发明公开了一种基于虚拟弹簧和球面调和函数力反馈模型建模方法，包括以下步骤：一，基于球面调和函数构建重构函数对目标软组织进行重构；二，利用Nyquist采样定律简化重构函数系数的求解，计算获得重构函数系数；三，在受力点处构建弹簧‑质点模型，基于各质点的法向形变量由重构函数曲面半径的变化来确定，计算形变区域面积；四，引入应力‑应变曲线求解虚拟弹簧弹性系数，通过虚拟弹簧弹性系数求解反馈力。本发明可以在较好表现模型轮廓和细节的前提下实现快速而逼真的力反馈计算；球面调和函数建模可保证模型的稳定性；弹簧‑质点模型使力反馈的求解方便快捷。

Description

基于虚拟弹簧和球面调和函数力反馈模型建模方法

技术领域

本发明属于软组织力触觉建模技术领域，具体涉及一种基于虚拟弹簧和球面调和函数力反馈模型建模方法。

背景技术

研究结果表明，近80％的医疗临床差错由人为因素造成，因此手术训练至关重要。传统的医疗培训成本高、风险高，通过虚拟手术训练可以有效地解决这一问题。

力反馈模型是虚拟触觉交互过程的关键。常用的物理模型中，弹簧-质点模型具有建模简单、计算灵活、速度快的优点，但该模型精度有限，稳定性差。有限元模型精度高，但计算相对复杂。因此，具有良好的交互功能和形变效果的柔体模型对于虚拟手术系统来说至关重要。

发明内容

本发明的目的在于克服现有技术中的不足，提供了一种基于虚拟弹簧和球面调和函数力反馈模型建模方法，这种算法可以在计算简单的同时确保高精度的变形模拟，在交互过程中，操作者可以感受到逼真的触觉信息。

为解决上述技术问题，本发明提供了一种基于虚拟弹簧和球面调和函数力反馈模型建模方法，其特征是，包括以下步骤：

步骤一，基于球面调和函数构建重构函数对目标软组织进行重构；

球面参数化后，用l阶m次的实值球面调和函数Y_l ^m(θ,φ)作为基函数表达重构函数f(θ,φ)为：

其中是重构函数系数，l为从0开始的正整数，|m|≤l；求解重构函数系数公式如下：

其中，是的共轭；

步骤二，利用Nyquist采样定律简化重构函数系数的求解，计算得到重构函数系数；

基于Nyquist采样定律，将球面调和函数系数的计算简化为函数采样加权和的形式，公式如下，

其中，

其中为m阶l次连带勒让德多项式，B为带宽；

通过以上公式计算获得重构函数系数，进而确定了重构函数；

步骤三，在受力点处构建弹簧-质点模型，基于各质点的法向形变量由重构函数曲面半径的变化来确定，计算形变区域面积；

其具体计算过程为，由弹簧-质点模型已知其中ε_j为第j层弹簧切向表面弹簧形变率，Δr_j为第j层表面弹簧形变量，r为表面弹簧长度，s为表面弹簧同心圆层数，是边界参数，切向表面弹簧形变率ε_j的计算公式为，

法向虚拟弹簧形变量L_j与切向表面弹簧形变率ε_j的关系可表示为，

将ε_j代入此公式可得下式，

其中，L_j是第j层虚拟弹簧形变量；

令j＝s，可简化为：

其中，代表区域边界的形变量，也是边界参数。

由上式可推导出边界参数s表达如下，

由上式可知，要计算边界参数s，需计算出接触点形变率ε₁；根据下式联立计算接触点形变率ε₁，

其中ΔZ_j是质点N_j处法向软组织形变量，可由N_j质点处对应的重构函数曲面半径的变化获得，已知边界参数s和表面弹簧长度r得出受影响软组织面积A，

A＝π(sr)²；

步骤四，引入应力-应变曲线求解虚拟弹簧弹性系数，通过虚拟弹簧弹性系数求解反馈力；其具体计算过程为：

以接触点N₁作为考察点，软组织弹性系数E可表示为，

其中，λ为应力，μ为应变，F_ext为施加的外力，A为步骤三中计算出的受影响软组织面积，Z₁为软组织形变前接触点的曲面半径，ΔZ₁为软组织形变后接触点改变的曲面半径；则虚拟弹簧弹性系数K：

由应力—应变的曲线图可以得出软组织弹性系数E，从而计算出虚拟弹簧弹性系数K；

通过已知的虚拟弹簧弹性系数K对反馈力F进行求解，根据已知技术可知第s层质点的力反馈可忽略不计，则反馈力的计算公式为：

式中，K为虚拟弹簧弹性系数，ΔZ_j为N_j质点处法向软组织形变量，可由重构函数对应点曲面半径的变化获得。

与现有技术相比，本发明所达到的有益效果是：本发明可以在较好表现模型轮廓和细节的前提下实现快速而逼真的力反馈计算；球面调和函数建模可保证模型的稳定性；弹簧-质点模型使力反馈的求解方便快捷。本发明算法可以在计算简单的同时确保高精度的变形模拟，在交互过程中，操作者可以感受到逼真的触觉信息。

附图说明

图1是本发明的流程示意图。

图2是虚拟弹簧的表面弹簧结构图。

图3是形变平面某一方向示例图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

如图1所示，本发明的一种基于虚拟弹簧和球面调和函数力反馈模型建模方法，包括以下步骤：

步骤一，基于球面调和函数构建重构函数对目标软组织进行重构。

基于球面调和函数的三维物体表达属于现有技术，就是将目标软组织表面用球面调和函数作为基函数表达出，使用球面调和函数进行建模可以更好地表现软组织的表面细节。球面参数化后软组织上的任意一个顶点都可以用f(θ,φ)表示，用l阶m次的实值球面调和函数(其中l为从0开始的正整数，表示频率波段指数；m为每个频带下的频率分量，且|m|≤l)展开球面参数化后，得到重构函数为：

球面调和函数的数学定义如下：

其中，为m阶l次连带勒让多项式，i为虚数单位；

其中是球面调和函数系数，即重构函数系数。求解球面调和函数系数公式如下：

其中，是的共轭。

步骤二，利用Nyquist采样定律简化重构函数系数的求解过程，计算得到重构函数系数。

由于球面调和函数对应每个固定频率的展开系数的模的平方和是旋转不变量，如果存在某个正数B>0，使得所有的l≥B时球面调和函数系数为0，则称B为带宽。由现有技术中Nyquist采样定律，将球面调和函数系数的计算简化为函数采样加权和的形式，公式如下，

其中，

其中为m阶l次连带勒让德多项式，f(θ_x,φ_y)在θ_x和φ_y均已知情况下，可由目标软组织的医学影像获得，i为虚数单位。

这种算法简化了繁琐的积分过程，加快了模型重建的速度。通过此简化过程计算获得重构函数系数，进而确定了重构函数，因此目标软组织受到外力产生的法向形变量则由球面调和函数曲面半径f(θ,φ)的变化获得。

步骤三，在受力点处建立虚拟弹簧模型，基于弹簧-质点模型计算形变区域面积。

基于弹簧-质点模型进行建模属于现有技术，参见图2，在受力接触面的接触点N₁处动态产生一个有虚拟弹簧的正六边形，以三角形为拓扑结构设计表面弹簧，以接触点为中心，用辐射状的表面同心圆分层。虚拟弹簧形变量由目标软组织模型球面调和函数曲面半径的改变给出，动态产生受力区域可以节省时间，同时减少计算量，方便力的求解。

在虚拟弹簧模型中，已知形变区域面积由s得出，s是一个边界参数，它表示在虚拟物体上以受力点为中心划分的表面弹簧同心圆层数，即外力影响的形变区域范围，共有s层质点(受力点N₁为第1层)、s层法向虚拟弹簧和s层切向表面弹簧。这样，第1层质点有1个质点，第2层质点有6个质点，第3层质点有12个质点，第s层质点有6(s-1)个质点。通过计算边界参数s可以算出外力影响的软组织表面区域面积：即共有s层表面弹簧，单个表面弹簧的长度为r，根据圆的面积公式可算出受影响表面区域面积。面积是之后应力计算中所需要的重要参数，同时能更直观地展现形变的影响范围。

形变平面内某一方向的模型如图3所示，Δr_j为第j层表面弹簧形变量，r为表面弹簧长度，令为第j层弹簧切向表面弹簧形变率，j＝1,2,3,...,s，j为自然数。根据现有技术已知ε_j表达如下，

法向虚拟弹簧形变量L_j与切向表面弹簧形变率ε_j的关系可表示为，

将ε_j代入L_j公式可得下式，

其中，L_j是第j层虚拟弹簧形变量；

令j＝s，可简化为，

其中，代表区域边界的形变量，也是边界参数。

由上式可推导出边界参数s表达如下，

由上式可知，要计算边界参数s，需计算出接触点形变率ε₁。由于碰撞点处法向软组织形变量，可由球面调和函数在碰撞点处曲面半径f(θ,φ)的变化确定，从而根据下式联立计算接触点形变率ε₁，

其中ΔZ_j是质点N_j处虚拟弹簧的形变量，也即质点N_j处法向软组织形变量，可由N_j质点处对应的球面调和函数曲面半径的变化获得，根据边界参数s和表面弹簧长度r得出受影响软组织面积A，

A＝π(sr)²

计算形变区域可以让力反馈模型更加逼真。

步骤四，引入应力-应变曲线求解虚拟弹簧弹性系数，通过虚拟弹簧弹性系数求解反馈力。

应力是外力和受力面积的比，应变是曲面半径形变量和原始曲面半径的比，由应力-应变曲线可得出软组织弹性系数E，从而计算出虚拟弹簧弹性系数K。

以接触点N₁作为考察点，软组织弹性系数E可表示为，

其中，λ—应力，μ—应变，F_ext—施加的力，A—上面计算出的受影响软组织面积，Z₁—软组织形变前接触点的曲面半径，ΔZ₁—软组织形变后接触点改变的曲面半径。计算虚拟弹簧弹性系数K，

由应力—应变的曲线图可以得出软组织弹性系数E，从而计算出虚拟弹簧弹性系数K。通过计算得到弹性系数K而不是人为规定，减少了模型中的人为干涉，使力反馈更自然。

通过已知的虚拟弹簧弹性系数K对反馈力F进行求解，根据已知技术可得第s层质点的力反馈可忽略不计，则反馈力的计算公式为：

式中，K—虚拟弹簧弹性系数，ΔZ_j—N_j质点处法向软组织形变量，可由球面调和函数对应点曲面半径的变化获得。

本发明算法可以在计算简单的同时确保高精度的变形模拟，在交互过程中，操作者可以感受到逼真的触觉信息。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。

Claims (2)

1.一种基于虚拟弹簧和球面调和函数力反馈模型建模方法，其特征是，包括以下步骤：

步骤一，基于球面调和函数构建重构函数对目标软组织进行重构；

步骤二，利用Nyquist采样定律简化重构函数系数的求解，计算得到重构函数系数；

步骤三，在受力点处构建弹簧-质点模型，基于各质点的法向形变量由重构函数曲面半径的变化来确定，计算形变区域面积；

其具体计算过程为，由弹簧-质点模型已知j＝1,2,3,...,s，其中ε_j为第j层弹簧切向表面弹簧形变率，Δr_j为第j层表面弹簧形变量，r为表面弹簧长度，s为表面弹簧同心圆层数，是边界参数，切向表面弹簧形变率ε_j的计算公式为，

${\mathrm{\ϵ}}_j=\frac{{\mathrm{\ϵ}}_1}{j+(j-1){\mathrm{\ϵ}}_1}$

法向虚拟弹簧形变量L_j与切向表面弹簧形变率ε_j的关系可表示为，

$L_j=r\sqrt{(2+{\mathrm{\ϵ}}_j){\mathrm{\ϵ}}_j}$

将ε_j代入此公式可得下式，

$\frac{L_j}{r}=\frac{1}{j+(j-1){\mathrm{\ϵ}}_1}\sqrt{\[2j+(2j-1){\mathrm{\ϵ}}_1\]{\mathrm{\ϵ}}_1}$

其中，L_j是第j层虚拟弹簧形变量；

令j＝s，可简化为：

$\frac{L_s}{r}\≈\sqrt{\frac{2{\mathrm{\ϵ}}_1}{(1+{\mathrm{\ϵ}}_1)s}}$

其中，代表区域边界的形变量，也是边界参数。

由上式可推导出边界参数s表达如下，

$s=\frac{2{\mathrm{\ϵ}}_1}{1+{\mathrm{\ϵ}}_1}{\left(\frac{r}{L_s}\right)}^2$

由上式可知，要计算边界参数s，需计算出接触点形变率ε₁；根据下式联立计算接触点形变率ε₁，

$\left\{\begin{array}{c}{\mathrm{\ΔZ}}_1=L_1+L_2+...+L_s\\ {\mathrm{\ΔZ}}_2=L_2+...+L_s\\ .\\ .\\ .\\ {\mathrm{\ΔZ}}_s=L_s\\ L_j=\frac{r}{j+(j-1){\mathrm{\ϵ}}_1}\sqrt{\[2j+(2j-1){\mathrm{\ϵ}}_1\]{\mathrm{\ϵ}}_1}\end{array}\right.$

其中ΔZ_j是质点N_j处法向软组织形变量，可由N_j质点处对应的重构函数曲面半径的变化获得，已知边界参数s和表面弹簧长度r得出受影响软组织面积A，

A＝π(sr)²；

步骤四，引入应力-应变曲线求解虚拟弹簧弹性系数，通过虚拟弹簧弹性系数求解反馈力；其具体计算过程为：

以接触点N₁作为考察点，软组织弹性系数E可表示为，

$E=\frac{\mathrm{\λ}}{\mathrm{\μ}}=\frac{\frac{F_{ext}}{A}}{\frac{{\mathrm{\ΔZ}}_1}{Z_1}}$

其中，λ为应力，μ为应变，F_ext为施加的外力，A为步骤三中计算出的受影响软组织面积，Z₁为软组织形变前接触点的曲面半径，ΔZ₁为软组织形变后接触点改变的曲面半径；则虚拟弹簧弹性系数K：

$K=\frac{F_{ext}}{{\mathrm{\ΔZ}}_1}=\frac{{\mathrm{\ΔZ}}_{1}AE}{Z_1}\frac{1}{{\mathrm{\ΔZ}}_1}=\frac{AE}{Z_1}$

由应力—应变的曲线图可以得出软组织弹性系数E，从而计算出虚拟弹簧弹性系数K；

通过已知的虚拟弹簧弹性系数K对反馈力F进行求解，根据已知技术可知第s层质点的力反馈可忽略不计，则反馈力的计算公式为：

$F=K\·{\mathrm{\ΔZ}}_1+K\·\underset{j=2}{\overset{s-1}{\Σ}}6\·\left(j-1\right)\·{\mathrm{\ΔZ}}_j$

式中，K为虚拟弹簧弹性系数，ΔZ_j为N_j质点处法向软组织形变量，可由重构函数对应点曲面半径的变化获得。

2.根据权利要求1所述的一种基于虚拟弹簧和球面调和函数力反馈模型建模方法，其特征是，在步骤二中，基于Nyquist采样定律，将重构函数系数的计算简化为函数采样加权和的形式，公式如下，

$C_l^m=\frac{\sqrt{2\mathrm{\π}}}{2B}\underset{x=0}{\overset{2B-1}{\Σ}}\underset{y=0}{\overset{2B-1}{\Σ}}{a}_x^{\left(B\right)}f({\mathrm{\θ}}_x,{\mathrm{\φ}}_y)e^{-{\mathrm{im\φ}}_y}{P}_l^m\left({\mathrm{cos\θ}}_x\right)$

其中，

$\begin{array}{cc}{\mathrm{\θ}}_x=\frac{\mathrm{\π}(2x+1)}{4B}& x=0,1,\mathrm{...},2B-1\end{array}$

$\begin{array}{cc}{\mathrm{\φ}}_y=\frac{y\mathrm{\π}}{B}& y=0,1,\mathrm{...2}B-1\end{array}$

其中f(θ,φ)为重构函数，P_l ^m为m阶l次连带勒让德多项式，B为带宽。