CN103869984A - 一种基于纹理图片的力触觉再现方法 - Google Patents

Abstract

本发明公开了一种基于纹理图片的力触觉再现方法，该力触觉再现方法包括以下过程：步骤10）提取纹理特征，纹理特征包含纹理表面微观高度和表征纹理粗糙度的动摩擦系数；步骤20）基于步骤10）提取的纹理特征，进行力触觉建模，测算纹理力合力；步骤30）将步骤20）测算的纹理力合力，通过手控器输出反馈给操作者。该基于纹理图片的力触觉再现方法，无需特定的测量设备即可获取虚拟纹理的动摩擦系数，无需专用的纹理表达装置，即可将虚拟纹理的粗糙感传递给人。

Description

一种基于纹理图片的力触觉再现方法

技术领域

本发明属于图像处理技术领域，具体来说，涉及一种基于纹理图片的力触觉再现方法。

背景技术

力触觉再现是人机交互的重要组成部分，力触觉再现所能表达的信息包括人在与环境物体接触过程中所感受到的表面形状、纹理，粗糙度、温度等特征信息。实现力触觉再现一般包含三个环节：物体物理特征信息提取、特征信息的力触觉建模以及力触觉刺激的渲染表达。

在虚拟现实中，纹理的表面微观高度和粗糙度是力触觉建模中的重要特征。获得纹理表面的表面微观高度的方法主要分为两类：一类是接触式测量，另一类是非接触式测量。接触式测量一般是借助传感器测量探针在纹理表面扫描时的接触状态，这种方法会改变系统的连接状态和动态特性，导致测量结果和原始结果之间存在偏差。非接触式测量的方法主要有两种：（1）运用图像处理方法从图像中提取纹理表面的轮廓特征，基于图像法快速便捷，硬件成本较低，但纹理图片会受光照、颜色、分辨率等因素的影响，使得实际提取的几何特征和真实纹理之间存在系统误差；（2）运用测微距设备对纹理表面进行扫描，以得到纹理表面的轮廓特征，这种方法多运用光学特性，具有非破坏性、测量精度高等特点，但硬件成本极高，比较耗费时间。

纹理粗糙度提取的方法同样分为接触式测量和非接触式测量两类，接触式测量：Jochen Lang等提出运用无线触觉传感器WHaT测量并记录手持探针轻敲纹理表面时的加速度信号和力信号，并用记录的信号计算纹理表面的粗糙度和刚度等。非接触式测量主要采用图像处理的方法，较为经典的纹理粗糙度度量算法有：Rosenfeld提出的基于最佳尺寸纹理粗糙度算法，该方法不但克服了粗糙度对度量尺度的依赖，而且能有效地处理显微纹理图像；Karu和Tamura等学者基于Rosenfeld纹理粗糙度算法，提高了度量纹理粗糙度的算法的速度；Pentland P提出了用分形维数作为度量一幅图像纹理粗糙度的方法；Julesz B等人提出利用纹理的二阶或高阶统计特性提取纹理区域的粗糙度分量的方法；Haralick R提出了基于共生矩阵的纹理粗糙度算法；Kaizer E提出了图像的自相关函数法；Jain R等人提出了功率谱的定义来度量纹理粗糙度。

常用的纹理力模型主要分有随机信号模型，正弦模型。（1）随机信号模型，该模型方法首先计算纹理图像像素点的均值m和方差s，设x为随机数，在国内外研究中较为常用的随机数模型有：均匀分布随机数（0～1均匀分布随机数）、正态分布随机数、对数正态分布随机数。纹理力的幅值可以被描述为均值，方差和随机数的线性组合。这类模型方法计算快速简便，但不能反映纹理的局部特征，即对纹理的区分度不够好，因为随机模型仅仅对纹理力进行了统计意义上的描述，而两幅完全不同的纹理图像可能会拥有相近甚至相同的统计变量（不限于均值和方差），会造成不同类型纹理输出相近或者相同纹理力的情况。（2）正弦模型，国外有学者将纹理力公式表述为正弦函数，其中正弦函数的幅值为与纹理图像像素平均值m正相关的常数，正弦函数周期是与纹理图像像素方差s成反比的常数。正弦模型具有一定的数学意义，因为如果将纹理图像看作是随空间坐标变化的函数，那么这个函数可以通过傅里叶变换变成一系列三角函数的线性组合。因此采用正弦函数来描述纹理变化是较为合适的。纹理力是反映纹理表面变化的力信号，所以同样采用正弦函数来描述纹理力变化。从公式可以看出,纹理图像像素值均值越大即纹理表面平均高度越大，周期性信号的幅值就越大,纹理越粗糙。纹理图像像素值即纹理表面高度变化越剧烈，则周期信号的频率越高纹理越细密。纹理力的大小随位置的变化而变化。该方法在一定程度上能较好的反映出纹理的特性，但确定性模型能够表达的纹理种类有限，仅能使力感觉有明显的区分度,但缺乏真实感，不能模拟真实的纹理接触。

发明内容

技术问题：本发明所要解决的技术问题是：提供一种基于纹理图片的力触觉再现方法，无需特定的测量设备即可获取虚拟纹理的动摩擦系数，无需专用的纹理表达装置，即可将虚拟纹理的粗糙感传递给人。

技术方案：为解决上述技术问题，本发明采用以下的技术方案：

一种基于纹理图片的力触觉再现方法，该力触觉再现方法包括以下过程：

步骤10）提取纹理特征，纹理特征包含纹理表面微观高度和表征纹理粗糙度的动摩擦系数；

步骤20）基于步骤10）提取的纹理特征，进行力触觉建模，测算纹理力合力；

步骤30）将步骤20）测算的纹理力合力，通过手控器输出反馈给操作者。

进一步，所述的步骤10）中，纹理表面微观高度采用从明暗恢复形状的方法提取。

进一步，所述的步骤10）的过程为：首先根据从明暗恢复形状的方法，假设光源为无限远处的点光源，反射模型为朗伯体反射模型，成像几何关系为正交投影，建立如式（1）所示的纹理图像的辐射照度方程；

$E(x,y)=\frac{1+p{p}_s+{\mathrm{qq}}_s}{\sqrt{1+p^2+q^2}\sqrt{1+{p_s}^2{q_s}^2}}$     式(1)

其中，E(x,y)是(x,y)点的灰度值，p表示表示虚拟纹理中(x,y)点单位法向量第一个分量，

q表示虚拟纹理中(x,y)点单位法向量第二个分量，p_s表示单位光源方向向量的第一个分量，

q_s表示单位光源方向向量的第二个分量，σ为光源方向与Z轴的夹角，τ为光源方向与X轴的夹角；(x,y)表示纹理图像中的任意一点的坐标，该坐标系的坐标系原点为图像左上角，坐标系横轴正方向为水平向右，为X轴；坐标系纵轴正方向为垂直向下，为Y轴；坐标系中正方向为竖直向上，为Z轴，Z轴垂直于虚拟纹理平面；

然后对辐射照度方程进行泰勒展开，并只保留展开式中线性项，采用后向有限差分的方法把保留下来的线性项离散化，形成离散化后的线性方程；最后采用迭代的方法求解离散化后的线性方程，得到纹理图像中每个像素点对应高度值，为纹理表面微观高度值。

进一步，所述的步骤10）中，表征纹理粗糙度的动摩擦系数包括以下步骤：

步骤101）采用Tamura算法提取纹理材质粗糙系数；

步骤102）利用纹理表面微观高度计算纹理表面任一点的局部粗糙系数；

步骤103）测算图像纹理表面任一点的动摩擦系数。

进一步，所述的步骤101）的过程为：首先，计算图像中大小为2^k×2^k个像素的活动窗口中像素的平均强度值A_k(x,y)，其中，k为与图像大小相关的参数；然后，计算图像中的每个像素在横轴方向上互不重叠的窗口之间的平均强度差E_kh，以及每个像素在纵轴方向上互不重叠的窗口之间的平均强度差E_kv，选择每个像素的平均强度差最大值E_k对应的k值，

E_k=max(E_kh,E_kv)    式（2）

式（2）中，E_k对应的k值为k′，S_best(x,y)=2^k′，S_best(x,y)表示纹理图像中点(x,y)的最佳像素点个数，最后，依据式（3）测算纹理材质粗糙系数F_crs：

$F_{\mathrm{crs}}=\frac{1}{m\×n}\underset{i_1=1}{\overset{m}{\mathrm{\Σ}}}\underset{i_1=1}{\overset{n}{\mathrm{\Σ}}}S(i_1,j_1)$     式（3）

式（3）中，m表示纹理图像横轴方向上的像素个数，n表示纹理图像纵轴方向上的像素个数，i₁表示取值范围从1～m的正整数，j₁表示表示取值范围从1～n的正整数，S(i₁,j₁)表示(i₁,j₁)点的最佳像素点个数尺寸。

进一步，所述的步骤102）的过程为：依据式（4）测算纹理表面任一点的局部粗糙系数；

$r(x,y)=\underset{j_2=1}{\overset{1}{\mathrm{\Σ}}}\underset{i_2=1}{\overset{1}{\mathrm{\Σ}}}|h(x,y)-h(x+i_2,y+j_2)|$     式（4）

式（4）中，r(x,y)表示纹理图像中(x,y)点的局部粗糙系数；h(x,y)表示纹理图像中(x,y)点归一化高度，h(x+i₂,y+j₂)表示纹理图像中(x+i₂,y+j₂)点归一化高度；i₂表示变量，取值为-1、0或1；j₂表示变量，取值为-1、0或1。

进一步，所述的步骤103）的过程为，依据式（5）测算图像纹理表面任一点的动摩擦系数μ(x,y)：

μ(x,y)=ρ×r(x,y)    式（5）

其中，ρ表示纹理材质粗糙度，r(x,y)表示纹理图像中的(x,y)点的局部粗糙系数。

进一步，所述的步骤20）的过程为：

步骤201）采用沿坐标轴的包围盒碰撞检测算法，检测虚拟探针是否与虚拟纹理触碰，如果虚拟探针与虚拟纹理触碰，则进入步骤202）；如果虚拟探针不与虚拟纹理触碰，则停留在本步骤，继续检测；

步骤202）测算虚拟探针与虚拟纹理在触碰点(x,y)时的纹理力合力：假设图像纹理表面刚性，且纹理表面的力符合胡克定律，依据式（6）测算法向纹理力：

${\stackrel{\→}{F}}_n(x,y)=k_n\×h(x,y)\×\stackrel{\→}{n}$     式（6）

其中，k_n为虚拟纹理的刚度系数，h(x,y)为(x,y)点的归一化高度，

为垂直方向单位向量，

为(x,y)点的法向纹理力向量；

切向纹理力是摩擦力，依据式（7）测算切向纹理力：

$\frac{{\mathrm{dF}}_f}{\mathrm{dx}}=k_n{|1-\frac{F_f}{F_c}\mathrm{sgn}\left(v\right)|}^{\mathrm{\η}}\mathrm{sgn}(1-\frac{F_f}{F_c}\mathrm{sgn}\left(v\right))$     式（7）

其中，F_f表示虚拟纹理和虚拟探针之间的摩擦力，k_n表示虚拟纹理的刚度系数，F_c表示纹理和虚拟探针之间的库伦摩擦力，F_c=μF_n，μ表示纹理动摩擦系数，F_n表示虚拟探针与虚拟纹理之间的法向纹理力，η表示取值范围在1到2之间的常数，sgn(w)为符号函数，w为sgn(w)中的变量， $\mathrm{sgn}\left(w\right)=\left\{\begin{array}{cc}-1& w<0\\ 0& w=0\\ 1& w>0\end{array}\right.;$ 在式（7）中，w为v或者

v表示虚拟探针在虚拟纹理表面的滑动速度；

依据式（8）测算纹理力合力向量

${\stackrel{\&RightArrow;}{F}}_t={\stackrel{\&RightArrow;}{F}}_f+{\stackrel{\&RightArrow;}{F}}_n$     式（8）

其中，

表示虚拟探针与虚拟纹理之间的切向纹理力向量。

表示虚拟探针与虚拟纹理之间的的法向纹理力向量。

有益效果：与现有技术相比，本发明的技术方案具有以下有益效果：

本发明的基于纹理图片的力触觉再现方法，无需特定的测量设备即可获取虚拟纹理的动摩擦系数，无需专用的纹理表达装置，即可将虚拟纹理的粗糙感传递给人。在本发明提出的纹理力模型中，纹理力被分解为法向分量和切向分量。首先假设纹理表面为刚性，则法向纹理力遵循胡克定律，其幅值与表面微观高度呈正比，法向力的大小主要反映了纹理表面凹凸的程度；切向纹理力表现为摩擦力，摩擦力模型采用Dahl模型，Dahl模型是工程力学领域常用的摩擦力模型，它具有计算简单，能够同时描述物体预滑动和宏观滑动等优点。与随机模型和正弦模型相比，本发明提出的模型充分利用纹理图像信息，对纹理的局部特征具有较好的区分度。经实验证明本发明提出的纹理特征提取与纹理力建模表达的方法具有良好的效果。

附图说明

图1为本发明的系统框图。

图2为本发明的作用力示意图。

具体实施方式

下面结合附图，对本发明的技术方案进行具体说明。

如图1和图2所示，本发明的一种基于纹理图片的力触觉再现方法，该力触觉再现方法包括以下过程：

步骤10）提取纹理特征，纹理特征包含纹理表面微观高度和表征纹理粗糙度的动摩擦系数；

纹理表面微观高度采用从明暗恢复形状的方法提取。具体过程为：首先根据从明暗恢复形状的方法，假设光源为无限远处的点光源，反射模型为朗伯体反射模型，成像几何关系为正交投影，建立如式（1）所示的纹理图像的辐射照度方程；

$E(x,y)=\frac{1+p{p}_s+{\mathrm{qq}}_s}{\sqrt{1+p^2+q^2}\sqrt{1+{p_s}^2{q_s}^2}}$     式（1）

其中，E(x,y)是(x,y)点的灰度值，p表示虚拟纹理中(x,y)点单位法向量第一个分量，

q表示虚拟纹理中(x,y)点单位法向量第二个分量，p_s表示单位光源方向向量的第一个分量，

q_s表示单位光源方向向量的第二个分量，

σ为光源方向与Z轴的夹角，τ为光源方向与X轴的夹角；(x,y)表示纹理图像中的任意一点的坐标，该坐标系的坐标系原点为图像左上角，坐标系横轴正方向为水平向右，为X轴；坐标系纵轴正方向为垂直向下，为Y轴；坐标系中正方向为竖直向上，为Z轴，Z轴垂直于虚拟纹理平面；然后对辐射照度方程进行泰勒展开，并只保留展开式中线性项，采用后向有限差分的方法把保留下来的线性项离散化，形成离散化后的线性方程；最后采用迭代的方法求解离散化后的线性方程，得到纹理图像中每个像素点对应高度值，为纹理表面微观高度值。

表征纹理粗糙度的动摩擦系数包括以下步骤101）至步骤103）。

步骤101）采用Tamura算法提取纹理材质粗糙系数：首先，计算图像中大小为2^k×2^k个像素的活动窗口中像素的平均强度值A_k(x,y)，其中，k为与图像大小相关的参数；然后，计算图像中的每个像素在横轴方向上互不重叠的窗口之间的平均强度差E_kh，以及每个像素在纵轴方向上互不重叠的窗口之间的平均强度差E_kv，选择每个像素的平均强度差最大值E_k对应的k值，

E_k=max(E_kh,E_kv)    式（2）

式（2）中，E_k对应的k值为k′，S_best(x,y)=2^k′，S_best(x,y)表示纹理图像中点(x,y)的最佳像素点个数，最后，依据式（3）测算纹理材质粗糙系数F_crs：

$F_{\mathrm{crs}}=\frac{1}{m\&times;n}\underset{i_1=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{i_1=1}{\overset{n}{\mathrm{\&Sigma;}}}S(i_1,j_1)$     式（3）

式（3）中，m表示纹理图像横轴方向上的像素个数，n表示纹理图像纵轴方向上的像素个数，i₁表示取值范围从1～m的正整数，j₁表示表示取值范围从1～n的正整数，S(i₁,j₁)表示(i₁,j₁)点的最佳像素点个数尺寸。

步骤102）利用纹理表面微观高度计算纹理表面任一点的局部粗糙系数：依据式（4）测算纹理表面任一点的局部粗糙系数；

$r(x,y)=\underset{j_2=1}{\overset{1}{\mathrm{\&Sigma;}}}\underset{i_2=1}{\overset{1}{\mathrm{\&Sigma;}}}|h(x,y)-h(x+i_2,y+j_2)|$     式（4）

式（4）中，r(x,y)表示纹理图像中(x,y)点的局部粗糙系数；h(x,y)表示纹理图像中(x,y)点归一化高度，h(x+i₂,y+j₂)表示纹理图像中(x+i₂,y+j₂)点归一化高度；i₂表示变量，取值为-1、0或1；j₂表示变量，取值为-1、0或1。

步骤103）测算图像纹理表面任一点的动摩擦系数：依据式（5）测算图像纹理表面任一点的动摩擦系数μ(x,y)：

μ(x,y)=ρ×r(x,y)    式（5）

其中，ρ表示纹理材质粗糙系数，ρ=F_crs，r(x,y)表示纹理图像中的(x,y)点的局部粗糙系数。

步骤20）基于步骤10）提取的纹理特征，进行力触觉建模，测算纹理力合力。

步骤20）具体包括步骤201）和步骤202）。

步骤201）采用沿坐标轴的包围盒碰撞检测算法，检测虚拟探针是否与虚拟纹理触碰，如果虚拟探针与虚拟纹理触碰，则进入步骤202）；如果虚拟探针不与虚拟纹理触碰，则停留在本步骤，继续检测；

步骤202）测算虚拟探针与虚拟纹理在触碰点(x,y)时的纹理力合力：假设图像纹理表面刚性，且纹理表面的力符合胡克定律，依据式（6）测算法向纹理力：

${\stackrel{\&RightArrow;}{F}}_n(x,y)=k_n\&times;h(x,y)\&times;\stackrel{\&RightArrow;}{n}$     式（6）

其中，k_n为虚拟纹理的刚度系数，h(x,y)为(x,y)点的归一化高度，为垂直方向单位向量，

为(x,y)点的法向纹理力向量；

切向纹理力是摩擦力，依据式（7）测算切向纹理力：

$\frac{{\mathrm{dF}}_f}{\mathrm{dx}}=k_n{|1-\frac{F_f}{F_c}\mathrm{sgn}\left(v\right)|}^{\mathrm{\&eta;}}\mathrm{sgn}(1-\frac{F_f}{F_c}\mathrm{sgn}\left(v\right))$     式（7）

其中，F_f表示虚拟纹理和虚拟探针之间的摩擦力，k_n表示虚拟纹理的刚度系数，F_c表示纹理和虚拟探针之间的库伦摩擦力，F_c=μF_n，μ表示纹理动摩擦系数，F_n表示虚拟探针与虚拟纹理之间的法向纹理力，η表示取值范围在1到2之间的常数，sgn(w)为符号函数，w为sgn(w)中的变量， $\mathrm{sgn}\left(w\right)=\left\{\begin{array}{cc}-1& w<0\\ 0& w=0\\ 1& w>0\end{array}\right.;$ 在式（7）中，w为v或者

v表示虚拟探针在虚拟纹理表面的滑动速度；

依据式（8）测算纹理力合力向量

${\stackrel{\&RightArrow;}{F}}_t={\stackrel{\&RightArrow;}{F}}_f+{\stackrel{\&RightArrow;}{F}}_n$     式（8）

其中，

表示虚拟探针与虚拟纹理之间的切向纹理力向量。

表示虚拟探针与虚拟纹理之间的的法向纹理力向量。

步骤30）将步骤20）测算的纹理力合力，通过手控器输出反馈给操作者。

上述再现方法的步骤10）中，式（1）

$\begin{array}{c}E(x,y)=R(p,q)=\frac{1+{\mathrm{pp}}_s+{\mathrm{qq}}_s}{\sqrt{1+p^2+q^2}\sqrt{1+{p_s}^2+{q_s}^2}}\\ =\frac{\cos \mathrm{\&sigma;}+p\cos \mathrm{\&tau;}\sin \mathrm{\&sigma;}+q\sin \mathrm{\&tau;}\sin \mathrm{\&sigma;}}{\sqrt{1+p^2+q^2}}\end{array}$

其中E(x,y)是(x,y)点的灰度值， $p=\frac{\&PartialD;z}{\&PartialD;x},q=\frac{\&PartialD;z}{\&PartialD;x},p_s=\frac{\cos \mathrm{\&tau;}\sin \mathrm{\&sigma;}}{\cos \mathrm{\&sigma;}},$

σ为光源方向与Z轴的夹角，τ为光源方向与X轴的夹角。

求解式（1）的方法考虑到非线性项对重构结果的影响不大，可以将反射函数进行泰勒展开。Tsai和Shah首次利用后向有限差分方法将反射函数离散化：

$p=\frac{\&PartialD;Z}{\&PartialD;x}=Z(x,y)-Z(x-1,y)$

$q=\frac{\&PartialD;Z}{\&PartialD;y}=Z(x,y)-Z(x,y-1)$

Z(x,y)表示纹理图像高度矩阵在(x,y)点的值，Z(x,y-1)表示纹理图像高度矩阵在(x-1,y)点的值，Z(x-1,y)表示纹理图像高度矩阵在(x,y-1)点的值。

将p,q代入式（1）可得

f(E(x,y),Z(x,y))

=E(x,y)-R(Z(x,y)-Z(x-1,y),Z(x,y)-Z(x,y-1))=0

式（9）f(E(x,y),Z(x,y))表示纹理图像灰度值与理论高度值关系的函数。

对于确定的(x,y)点，其灰度E(x,y)是确定的，假设已知第n-1次迭代求出的高度矩阵在(x,y)点的值为Z^n-1(x,y)，对式（9）泰勒展开：

$Z^n(x,y)=Z^{n-1}(x,y)+\frac{-f(E(x,y),Z^{n-1}(x,y))}{\frac{\mathrm{df}(E(x,y),Z^{n-1}(x,y))}{\mathrm{dZ}(x,y)}}$     式(10)

Zⁿ(x,y)表示第n次迭代求出的高度矩阵在(x,y)点的值，f(E(x,y),Z^n-1(x,y))表示纹理图像灰度值与当前迭代出来的高度值关系的函数。由于纹理高度无法直接获得，以第n次迭代高度代替纹理高度。其中，

$\frac{\mathrm{df}(E(x,y),Z^{n-1}(x,y))}{\mathrm{dZ}(x,y)}=-1\&times;(\frac{(p_s+q_s)}{\sqrt{p^2+q^2+1}\sqrt{p_s^2{q}_s^2+1}}-\frac{(p+q)({\mathrm{pp}}_s+{\mathrm{qq}}_s+1)}{\sqrt{{(p^2+q^2+1)}^3}\sqrt{p_s^2+q_s^2+1}})$

由于局部表面的法向量可能与光照方向相同，即p_s=p，q_s=q，造成式（10）中的分母为0。因此，在式（2）的基础上构造出如下的迭代公式：

Zⁿ(x,y)=Z^n-1(x,y)+Kⁿ(-f(E(x,y),Z^n-1(x,y)))    式（11）

Kⁿ需要满足三个条件：（1）Kⁿ应该约等于的倒数；（2）当

逼近0时，Kⁿ等于0；（3）当Zⁿ(x,y)逼近真实Z(x,y)时，Kⁿ等于0。因此，定义Kⁿ：

$K^n=\frac{S_{x,y}^n{M}_{x,y}}{W_{x,y}+S_{x,y}^n{M}_{x,y}^2}$     式（12）

其中， $M_{x,y}=\frac{\mathrm{df}\left(Z^{n-1}(x,y)\right)}{\mathrm{dZ}(x,y)},S_{x,y}^n=E\left[{(Z^n(x,y)-Z(x,y))}^2\right],$ E是数学期望算子，W_x,y是一个非常小的非零数，W_x,y的值小于0.01×0.01。

根据式（10），经过数次迭代后，Zⁿ(x,y)就是对应像素(x,y)的高度值。

$h(x,y)=\frac{Z(x,y)}{Z_{\max }}$     式（13）

其中，h(x,y)为纹理表面归一化高度，Z_max为像素点最大高度。

在步骤101）纹理材质粗糙系数计算中：Tamura算法提取纹理材质粗糙系数可以分为以下几个步骤进行。首先，计算图像中大小为2^k×2^k个像素的活动窗口中像素的平均强度值，即有:

$A_k(x,y)=\underset{i=1z-2^{k-1}}{\overset{{x+2}^{k-1}-1}{\mathrm{\&Sigma;}}}\underset{i=y-2^{k-1}}{\overset{{y+2}^{k-1}-1}{\mathrm{\&Sigma;}}}g(i,j)/2^k$     式（14）

由于必须保证g(i+2^k-1-1,j+2^k-1-1),g(i+2^k-1,j+2^k-1)等像素点必须在给定图像的有效像素范围内，所以k不宜取过大，一般取k=1,2,×××,5。然后，对于每个像素分别计算它在水平和垂直方向上互不重叠的窗口之间的平均强度差

横轴方向：E_kh(x,y)=|A_k(x+2^k-1,y)-A_k(x-2^k-1,y)|    式（15）

纵轴方向：E_kv(x,y)=|A_k(x,y+2^k-1)-A_k(x,y-2^k-1)|    式（16）

其中对于每个像素，选择能使E值达到最大（无论方向）的k值，即E_k=max(E_kh,E_kv)，E_k对应的k值为k′，S_best(x,y)=2^k′，S_best(x,y)表示纹理图像中点(x,y)的最佳像素点个数。

最后纹理材质粗糙系数可以通过计算整幅图像中S_best的平均值来获得，表达式为

$F_{\mathrm{crs}}=\frac{1}{m\&times;n}\underset{i_1=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j_1=1}{\overset{n}{\mathrm{\&Sigma;}}}S(i_1,j_1).$

其中，m表示纹理图像横轴方向上的像素个数，n表示纹理图像纵轴方向上的像素个数，i₁表示取值范围从1～m的正整数，j₁表示表示取值范围从1～n的正整数，S(i₁,j₁)表示(i₁,j₁)点的最佳像素点个数尺寸。

动摩擦系数μ(x,y)用于表征纹理的粗糙程度。一般来说以丝绸等丝织品为代表的光滑纹理，表面微观高度变化微小，此类纹理根据本专利提出的方法获得的动摩擦系数较小，力触觉再现后会给人较为顺滑的感觉；反之如砂纸为代表的粗糙纹理，表面微观高度变化剧烈，此类纹理根据本专利提出的方法获得的动摩擦系数较大，力触觉再现后会给人较为粗糙的感觉。

上述基于纹理图片的力触觉再现方法，采用图像处理算法从二维纹理图像中提取纹理表面微观高度。采用Tamura算法计算纹理材质粗糙系数，根据纹理表面微观高度的梯度计算接触点的局部粗糙系数。将接触点的动摩擦系数定义为纹理材质粗糙系数和局部粗糙系数的乘积。然后按垂直和水平两个方向进行力触觉建模，其中法向力的大小符合胡克定律，切向力等效为摩擦力遵循Dahl摩擦力模型，最终将计算出的力由通用手控器传递到人手。该方法具有硬件成本低、再现方便快捷，真实感高的优点。所提出的方法可以广泛应用与虚拟手术仿真、网上商店和虚拟博物馆浏览，远程机器人操作等需要进行纹理触觉再现的应用领域。

下面例举一具体实施例。

（1）提取纹理表面微观高度：

以一副4×4的灰度图为例，该图左半边为黑色，右半边为白色，即该图像素矩阵第一列到第八列灰度值为0，第九列到第十六列灰度值为255。即 $E(x,y)=\left\{\begin{array}{cc}0& 1\&le;x\&le;2\\ 255& 3\&le;x\&le;4\end{array}\right.,$ 计算前把灰度归一化，所以归一化后 $E(x,y)=\left\{\begin{array}{cc}0& 1\&le;x\&le;2\\ 1& 3\&le;x\&le;4\end{array}\right..$ 由于SFS算法求解高度的过程是一个迭代过程，一般需要10次左右的迭代才能求解出高度，对16个像素点进行10次迭代计算计算量十分巨大。对于任意像素点的每一次迭代计算过程其实十分相似，所以本专利仅详细计算某个像素点的一轮迭代计算方法。现在以像素矩阵中第2行第2列的元素为例，介绍从灰度矩阵E(x,y)计算高度矩阵Z(x,y)方法的首轮迭代计算。首先对公式（1）～（7）出现的参数进行初始化。迭代开始前的高度矩阵初始为全0矩阵，即Z⁰(x,y)=0,x,y∈[1,4]

假设光源方向为表示为(0.01,0.01,-1)，所以p_s=0.01,q_s=0.01。p=Z⁰(2,2)-Z⁰(1,2)=0,q=Z⁰(2,2)-Z⁰(2,1)=0，W_2,2=0.001²。

$\begin{array}{c}\frac{\mathrm{df}(E\left(\mathrm{2,2}\right),Z^1\left(\mathrm{2,2}\right))}{\mathrm{dZ}\left(\mathrm{2,2}\right)}=-1\&times;(\frac{(p_s+q_s)}{\sqrt{p^2+q^2+1}\sqrt{p_s^2+q_s^2+1}}-\frac{(p+q)({\mathrm{pp}}_s+{\mathrm{qq}}_s+1)}{\sqrt{{(p^2+q^2+1)}^3}\sqrt{p_s^2+q_s^2+1}})\\ =-1\&times;(\frac{0.01+0.01}{\sqrt{0+0+1}\&times;\sqrt{{0.01}^2+{0.01}^2+1}}-\frac{(0+0)(0+0+1)}{\sqrt{1}\&times;\sqrt{{0.01}^2+{0.01}^2+1}})\\ =\frac{0.98}{\sqrt{1+2\&times;{0.01}^2}}\\ =0.5734\end{array}$

$\begin{array}{c}f(E\left(\mathrm{2,2}\right),Z^0\left(\mathrm{2,2}\right))=-(E\left(\mathrm{2,2}\right)-\frac{1+{\mathrm{pp}}_s+{\mathrm{qq}}_s}{\sqrt{1+p^2+q^2}\sqrt{1+{p_s}^2{q_s}^2}})\\ =(0-\frac{1+0+0}{\sqrt{1}\&times;\sqrt{1+2*{0.01}^2}})=0.9999\end{array}$

在计算完成前无法确定，所以在计算中我们

用迭代值不断逼近真值，迭代公式为：

$S_{x,y}^n=(1-\frac{\mathrm{df}(E(x,y),Z^{n-1}(x,y))}{\mathrm{dZ}(x,y)}{k}^{n-1})S_{x,y}^{n-1}$

假设所以

$\begin{array}{c}{K}^1=\frac{S_{\mathrm{2,2}}^1{M}_{\mathrm{2,2}}}{W_{\mathrm{2,2}}+S_{\mathrm{2,2}}^1{M}_{\mathrm{2,2}}\&times;M_{\mathrm{2,2}}}=\frac{1\&times;\frac{1}{\frac{\mathrm{df}(E\left(\mathrm{2,2}\right),Z^1\left(\mathrm{2,2}\right))}{\mathrm{dZ}\left(\mathrm{2,2}\right)}}}{0{.001}^2+\frac{1}{\frac{\mathrm{df}(E\left(\mathrm{2,2}\right),Z^1\left(\mathrm{2,2}\right))}{\mathrm{dZ}\left(\mathrm{2,2}\right)}}\&times;\frac{1}{\frac{\mathrm{df}(E\left(\mathrm{2,2}\right),Z^1\left(\mathrm{2,2}\right))}{\mathrm{dZ}\left(\mathrm{2,2}\right)}}}\\ =0.5728\end{array}$

Z¹(2,2)=Z⁰(x,y)+K¹(-f(E(2,2),Z⁰(2,2)))

=0+0.5728′(-0.9999)

=-0.5727

这样，对于像素矩阵第2行第2列像素点的求解高度的首轮迭代运算就完成了，Z¹(2,2)表示各像素点之间的相对高度，可以为负值，计算完成后会将所有像素点的高度归一化。对于其它像素点的各轮迭代计算只要仿照上述的计算过程就可完成计算。

最后迭代结果 $H=\left[\begin{array}{cccc}-2.9287& -2.9287& -2.9287& -2.9287\\ -2.9287& -2.9287& 1.2306& 0.003\\ -2.9287& -2.9287& 1.0173& 0.003\\ -2.9287& -2.9287& -2.9287& 2.9287\end{array}\right],$ 归一化后

$H=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 1& 0.7042\\ 0& 0& 0.9487& 0.7042\\ 0& 0& 0& 0\end{array}\right]$

（2）纹理粗糙度：

设k=1、2，

$A_1\left(\mathrm{2,2}\right)=\frac{\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{2}{\mathrm{\&Sigma;}}}g(i,j)}{2}=0,A_1\left(\mathrm{2,3}\right)=\frac{\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}\underset{j=2}{\overset{3}{\mathrm{\&Sigma;}}}g(i,j)}{2}=255,A_1\left(\mathrm{2,4}\right)=\frac{\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}\underset{j=3}{\overset{4}{\mathrm{\&Sigma;}}}g(i,j)}{2}=510$

$A_1(,3)=\frac{\underset{i=2}{\overset{3}{\mathrm{\&Sigma;}}}\underset{j=2}{\overset{1}{\mathrm{\&Sigma;}}}g(i,j)}{2}=0,A_1\left(\mathrm{3,3}\right)=\frac{\underset{i=2}{\overset{3}{\mathrm{\&Sigma;}}}\underset{j=2}{\overset{3}{\mathrm{\&Sigma;}}}g(i,j)}{2}=255,$

$A_1\left(\mathrm{3,4}\right)=\frac{\underset{i=2}{\overset{3}{\mathrm{\&Sigma;}}}\underset{j=3}{\overset{4}{\mathrm{\&Sigma;}}}g(i,j)}{2}=510,$

$A_1(4,2)=\frac{\underset{i=3}{\overset{4}{\mathrm{\&Sigma;}}}\underset{j=2}{\overset{1}{\mathrm{\&Sigma;}}}g(i,j)}{2}=0,A_1\left(\mathrm{4,3}\right)=\frac{\underset{i=3}{\overset{4}{\mathrm{\&Sigma;}}}\underset{j=2}{\overset{1}{\mathrm{\&Sigma;}}}g(i,j)}{2}=255,$

$A_1(4,4)=\frac{\underset{i=3}{\overset{4}{\mathrm{\&Sigma;}}}\underset{j=2}{\overset{1}{\mathrm{\&Sigma;}}}g(i,j)}{2}=510,A_2\left(\mathrm{3,3}\right)=\frac{\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{4}{\mathrm{\&Sigma;}}}g(i,j)}{2}=510,$

E_1h(3,3)=|A₁(2,3)-A₁(4,3)|=0,E_1v(3,3)=|A₁(3,2)-A₁(3,4)|=510，

所以E₁=510。对于4×4的图像，只有(3,3)点可以完成式（13）～（15）的计算，所以S_best(3,3)=2， $F_{\mathrm{crs}}=\frac{3}{4\&times;4}=0.1875.$

（3）计算虚拟力：计算(3,3)点的纹理力合力，

r(3,3)=|0.9487-1|+|0.9487-0|+|0.9487-0|+|0.9487-0.7042|

=2.1932

μ(3,3)=0.1875×2.1932

=0.4112

设纹理刚度k_n=1，F_n(3,3)=0.9487N，式（7）中取η=1，式（7）的计算公式可简化为库伦摩擦力公式，即F_f(x,y)=F_c(x,y)=μ(x,y)×F_n(x,y)，

F_f(3,3)=0.4112×0.9487=0.3901N,

$F_t\left(\mathrm{3,3}\right)=\sqrt{F_n{\left(\mathrm{3,3}\right)}^2+F_f{\left(\mathrm{3,3}\right)}^2}=\sqrt{{0.9487}^2+0{.4112}^2}=1.0340N,$ F_n(x,y)表示向量的大小，F_f(x,y)表示向量的大小，F_t(x,y)表示向量

的大小。

Claims (8)

1.一种基于纹理图片的力触觉再现方法，其特征在于，该力触觉再现方法包括以下过程：

步骤10）提取纹理特征，纹理特征包含纹理表面微观高度和表征纹理粗糙度的动摩擦系数；

步骤20）基于步骤10）提取的纹理特征，进行力触觉建模，测算纹理力合力；

步骤30）将步骤20）测算的纹理力合力，通过手控器输出反馈给操作者。

2.按照权利要求1所述的基于纹理图片的力触觉再现方法，其特征在于，所述的步骤10）中，纹理表面微观高度采用从明暗恢复形状的方法提取。

3.按照权利要求2所述的基于纹理图片的力触觉再现方法，其特征在于，所述的步骤10）的过程为：首先根据从明暗恢复形状的方法，假设光源为无限远处的点光源，反射模型为朗伯体反射模型，成像几何关系为正交投影，建立如式（1）所示的纹理图像的辐射照度方程；

$E(x,y)=\frac{1+p{p}_s+{\mathrm{qq}}_s}{\sqrt{1+p^2+q^2}\sqrt{1+{p_s}^2{q_s}^2}}$     式（1）

其中，E(x,y)是(x,y)点的灰度值，p表示表示虚拟纹理中(x,y)点单位法向量第一个分量，

q表示虚拟纹理中(x,y)点单位法向量第二个分量，

p_s表示单位光源方向向量的第一个分量，

q_s表示单位光源方向向量的第二个分量，σ为光源方向与Z轴的夹角，τ为光源方向与X轴的夹角；(x,y)表示纹理图像中的任意一点的坐标，该坐标系的坐标系原点为图像左上角，坐标系横轴正方向为水平向右，为X轴；坐标系纵轴正方向为垂直向下，为Y轴；坐标系中正方向为竖直向上，为Z轴，Z轴垂直于虚拟纹理平面；

然后对辐射照度方程进行泰勒展开，并只保留展开式中线性项，采用后向有限差分的方法把保留下来的线性项离散化，形成离散化后的线性方程；最后采用迭代的方法求解离散化后的线性方程，得到纹理图像中每个像素点对应高度值，为纹理表面微观高度值。

4.按照权利要求1所述的基于纹理图片的力触觉再现方法，其特征在于，所述的步骤10）中，表征纹理粗糙度的动摩擦系数包括以下步骤：

步骤101）采用Tamura算法提取纹理材质粗糙系数；

步骤102）利用纹理表面微观高度计算纹理表面任一点的局部粗糙系数；

步骤103）测算图像纹理表面任一点的动摩擦系数。

5.按照权利要求1所述的基于纹理图片的力触觉再现方法，其特征在于，所述的步骤101）的过程为：首先，计算图像中大小为2^k×2^k个像素的活动窗口中像素的平均强度值A_k(x,y)，其中，k为与图像大小相关的参数；然后，计算图像中的每个像素在横轴方向上互不重叠的窗口之间的平均强度差E_kh，以及每个像素在纵轴方向上互不重叠的窗口之间的平均强度差E_kv，选择每个像素的平均强度差最大值E_k对应的k值，

E_k=max(E_kh,E_kv)    式（2）

式（2）中，E_k对应的k值为k′，S_best(x,y)=2^k′，S_best(x,y)表示纹理图像中点(x,y)的最佳像素点个数，最后，依据式（3）测算纹理材质粗糙系数F_crs：

$F_{\mathrm{crs}}=\frac{1}{m\&times;n}\underset{i_1=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{i_1=1}{\overset{n}{\mathrm{\&Sigma;}}}S(i_1,j_1)$     式（3）

式（3）中，m表示纹理图像横轴方向上的像素个数，n表示纹理图像纵轴方向上的像素个数，i₁表示取值范围从1～m的正整数，j₁表示表示取值范围从1～n的正整数，S(i₁,j₁)表示(i₁,j₁)点的最佳像素点个数尺寸。

6.按照权利要求4所述的基于纹理图片的力触觉再现方法，其特征在于，所述的步骤102）的过程为：依据式（4）测算纹理表面任一点的局部粗糙系数；

$r(x,y)=\underset{j_2=1}{\overset{1}{\mathrm{\&Sigma;}}}\underset{i_2=1}{\overset{1}{\mathrm{\&Sigma;}}}|h(x,y)-h(x+i_2,y+j_2)|$     式（4）

式（4）中，r(x,y)表示纹理图像中(x,y)点的局部粗糙系数；h(x,y)表示纹理图像中(x,y)点归一化高度，h(x+i₂,y+j₂)表示纹理图像中(x+i₂,y+j₂)点归一化高度；i₂表示变量，取值为-1、0或1；j₂表示变量，取值为-1、0或1。

7.按照权利要求4所述的基于纹理图片的力触觉再现方法，其特征在于，所述的步骤103）的过程为，依据式（5）测算图像纹理表面任一点的动摩擦系数μ(x,y)：

μ(x,y)=ρ×r(x,y)    式（5）

其中，ρ表示纹理材质粗糙度，r(x,y)表示纹理图像中的(x,y)点的局部粗糙系数。

8.按照权利要求1所述的基于纹理图片的力触觉再现方法，其特征在于，所述的步骤20）的过程为：

步骤201）采用沿坐标轴的包围盒碰撞检测算法，检测虚拟探针是否与虚拟纹理触碰，如果虚拟探针与虚拟纹理触碰，则进入步骤202）；如果虚拟探针不与虚拟纹理触碰，则停留在本步骤，继续检测；

步骤202）测算虚拟探针与虚拟纹理在触碰点(x,y)时的纹理力合力：假设图像纹理表面刚性，且纹理表面的力符合胡克定律，依据式（6）测算法向纹理力：

${\stackrel{\&RightArrow;}{F}}_n(x,y)=k_n\&times;h(x,y)\&times;\stackrel{\&RightArrow;}{n}$     式（6）

其中，k_n为虚拟纹理的刚度系数，h(x,y)为(x,y)点的归一化高度，

为垂直方向单位向量，为(x,y)点的法向纹理力向量；

切向纹理力是摩擦力，依据式（7）测算切向纹理力：

$\frac{{\mathrm{dF}}_f}{\mathrm{dx}}=k_n{|1-\frac{F_f}{F_c}\mathrm{sgn}\left(v\right)|}^{\mathrm{\&eta;}}\mathrm{sgn}(1-\frac{F_f}{F_c}\mathrm{sgn}\left(v\right))$     式（7）

其中，F_f表示虚拟纹理和虚拟探针之间的摩擦力，k_n表示虚拟纹理的刚度系数，F_c表示纹理和虚拟探针之间的库伦摩擦力，F_c=μF_n，μ表示纹理动摩擦系数，F_n表示虚拟探针与虚拟纹理之间的法向纹理力，η表示取值范围在1到2之间的常数，sgn(w)为符号函数，w为sgn(w)中的变量， $\mathrm{sgn}\left(w\right)=\left\{\begin{array}{cc}-1& w<0\\ 0& w=0\\ 1& w>0\end{array}\right.;$ 在式（7）中，w为v或者

v表示虚拟探针在虚拟纹理表面的滑动速度；

依据式（8）测算纹理力合力向量

${\stackrel{\&RightArrow;}{F}}_t={\stackrel{\&RightArrow;}{F}}_f+{\stackrel{\&RightArrow;}{F}}_n$     式（8）

其中，

表示虚拟探针与虚拟纹理之间的切向纹理力向量。

表示虚拟探针与虚拟纹理之间的的法向纹理力向量。

Images