结合空域噪声信息的小波脊提取算法

王 勇，邹 辉，饶勤菲，王李福

(重庆理工大学计算机科学与工程学院 重庆 巴南区 400054)

基于光学的物体三维形貌测量因其非接触、无损、测量速度快、自动化程度高、精度适中等特点而广泛应用于3D打印、机器视觉、虚拟现实等领域[1-2]。小波变换轮廓术[3-4]作为一种光学物体三维测量技术，因其抗噪性强、重建精度高、仅需一幅图像就能获取相位等优点而受到学者的广泛研究。

小波变换轮廓术的实质就是正确提取小波脊的位置，从而得到最佳伸缩尺度和对应相位信息。然而，在现实图像的采集过程，由于受到图像传感器、传输信道、解码处理等各种因素的影响，图像在形成和传输的过程会引入噪声，特别是椒盐噪声[5]。小波变换轮廓术在提取小波脊的过程中极易受此类噪声的影响，造成小波脊提取不准确，进而影响测量精度。因此，如何抑制噪声并准确提取小波脊，成为小波变换轮廓术的重点之一。

文献[6]利用连续小波变换平面上分布参数的相位信息提取脊线，但该方法仅能应用于没有噪声的信号，在噪声环境下，效果较差。文献[7]提出了一种在时间尺度空间利用小波变换平面上的相位信息迭代来提取小波脊的方法，由于小波变化的相位受噪声影响较大，该方法提取的脊线也不准确。文献[8]提出了模极大值算法，该方法利用小波系数模值的极值来搜寻脊点，并将极值的最大值做为脊点。但该方法抗噪性能较弱，因此，该研究者对该方法进行了改进，提出了Crazy-Climber算法。该算法优化了极值的搜索，强化了结果的全局最优化，降低了对噪声的敏感性，但是较为复杂。文献[9]提出了使用代价函数提取小波脊，该算法结合小波系数模极大值和小波脊的光滑特性，将小波脊的提取转化为最优解的问题，并使用动态规划思想求解小波脊。该算法具有良好的抗噪性和鲁棒性，但代价函数仍有缺陷。文献[10]指出基于代价函数的小波脊相位提取方法是目前较为有效、稳健的方法。

本文基于文献[9]提出的代价函数法，提出了一种新的代价函数提取小波脊的算法。该算法对原代价函数的小波变换系数幅值和尺度参数曲线梯度两项指标使用Feature Scaling方法进行调整，为避免两项权重过于拟合，造成小波脊受噪声影响波动的问题，本文结合空域噪声信息对调整后的尺度参数曲线梯度进行补偿，以获得更加合理的代价值，从而准确提取小波脊。实验证实了本文算法的有效性。

1 小波变换轮廓术基本原理

小波变换轮廓术的测量系统的光路原理和傅里叶变换轮廓术相同[2,4]。成像装置获得的变形条纹图一般可表示为：

式中，a(x,y)为背景分量；b(x,y)为条纹可见度；f0为空间载波频率；φ(x,y)为物体高度引起的调制相位；c(x,y)为随机噪声。其中，f0必须满足下式才能恢复相位[11]：

条纹图一维连续小波变换为[12]：

式中，a为尺度因子；b为平移因子；f(x)为取条纹图一行的信号；ψ(x)为母小波；“*”表示共轭。

Morlet小波在空域和频域中具有良好的局部性，文献[13]通过对几种小波实验对比，指出复Morlet相对于其他小波在相位恢复中更具优势，故本文选用复Morlet作为母小波，其表示如下：

式中，fb为带宽参数；fc为小波中心频率参数。

通过使用复Morlet小波对条纹图信号进行连续小波变换，可得到小波变换系数W(a,b)，其实部记为虚部为则幅值和相位分别为[4,13]：

由式(3)可以看出，小波变换系数W(a,b)为信号和小波序列函数内积的结果，它反映了信号和小波序列函数的相似程度。因此，当信号的局部频率和相应尺度的小波函数振荡频率相近或者相同时，较大，即相应的幅值A(a,b)较大。小波脊即为在小波平移方向上，不同位置的幅值最大值的连线。脊处对应的尺度即为最佳伸缩尺度，对应的相位即为包含物体高度信息的调制相位。

2 结合空域噪声信息的代价函数提取小波脊

条纹图受到噪声影响后，小波变换系数幅值会产生多个局部极大值点，若直接采用模极大值提取相位，易造成脊点的误选，导致小波脊定位不准确。根据条纹图瞬时频率的连续特性，文献[9]提出了基于代价函数的小波脊提取算法，该算法通过计算最小代价值来确定小波脊。

代价函数定义为：

式中，b为平移因子；φ(b)为位置b处的尺度因子；表示在尺度为φ(b)，位置为b处的小波变换系数幅值；C0和C1分别表示小波变换系数幅值和尺度参数曲线梯度的权值系数。

在实际应用中，处理的对象为数字图像，位置是离散化的，因此需要对式(7)进行离散化，虽然原理上可计算任意位置的代价值，但实际上这是不现实也是不必要的。式(7)离散化形式为：

式中，L表示小波变换系数幅值矩阵的宽度；b表示平移因子。

由式(8)可以看出，其选取的小波脊代价函数较小且较为光滑，避免了因噪声产生较大幅值造成的脊线不光滑。

如式(8)所示，在一般研究中，都默认设置C0和的值为1。小波变换系数幅值通常在0～1之间。小波变换尺度φ(b)无论选择步长为1还是步长为指数增长形式[10]，一般大于1。因此通常有：

由式(8)选择的小波脊过于偏向平滑，受到噪声影响后，易误选使脊线更加平滑的噪点作为脊点。因此，本文提出一种新的结合空域噪声信息和Feature Scaling方法的代价函数，其定义为：

式中，Maxs和Mins分别表示候选脊点对应小波变换系数幅值矩阵除第一列；的最大值和最小值；Maxb和Minb分别表示候选脊点对应的尺度因子矩阵梯度的最大值和最小值；K为λ的调整因子，一般为10n，n为整数，使得为噪声程度系数，具体判断标准可根据实际情况定义使用，本文中取噪声数与总像素比，其表达式为：

式中，Count为检测的噪声总数；Pixel_num为变形条纹的总像素数。

由式(10)可以看出，首先使用Feature Scaling方法对小波变换系数幅值和尺度参数曲线梯度两项特征因子进行特征缩放，权衡两项特征因子的影响比重，然后利用空域噪声信息对尺度参数曲线梯度因子进行补偿。在较低噪声情况下，小波脊受噪声的影响较小，对尺度参数曲线梯度因子补偿较小，使两项因子权重等阶，既使用模极大值提取小波脊，又利用小波脊线的光滑性过滤噪声。随着噪声的增加，噪声会产生较多属于极值的小波变换系数幅值，更易干扰小波脊的正确提取，此时增加对尺度参数曲线梯度因子的补偿，加大该因子的权重，即增强利用小波脊线光滑特性滤除噪声，从而正确定位小波脊。

本文结合条纹图灰度的连续特性，选取噪声检测中常用的基于梯度的判别法来判断噪点[15]，从而计算λ值。噪点分析示意图如图1所示，相关公式为：

式中，L(c)为噪点梯度和；c表示检测点x0的邻接区域；x0(c)表示检测点x0的灰度值；表示x0邻接点的灰度值。

如图1a所示，若L(c)大于某个阈值，可判定x0为噪点，对于边缘像素点，如图1b所示，L(c)远小于设定阈值。

图1 噪点分析示意图

获取λ步骤如下：

1)提取变形条纹图的一行。

2)读取第j个像素点，其中为条纹图一行的长度。设读取的像素点为A，其灰度值为pixel[i][j]。A后第2个像素点B的灰度值为计算A与B的灰度差值d=

3)分别计算A和B的噪点梯度和L1、L2。

小波脊提取步骤如下：

1)分别计算λ值和K值。

2)提取变形条纹图的一行对其进行连续小波变换，获得小波变换系数矩阵Coefficient(m,n)，由式(5)可求得幅值矩阵Modulus(m,n)。其中，m为尺度a的维数，n为条纹图一行的长度。

3)沿平移因子方向对Modulus(m,n)第i列做尺度-模值曲线，求取曲线局部极大值点和拐点，并将大于均值的点作为候选脊点。

4)重复步骤 3)直到i=n，获得候选脊点矩阵CandidateRidgePoints，进而获得对应的小波变换幅值矩阵Can_Magnitudes和尺度因子矩阵Can_Scales。

5)根据Can_Magnitudes和Can_Scales分别计算Mins、Maxs、Minb、Maxb。

6)设置和CandidateRidgePoints相同大小的代价值矩阵Costs，并将第一列代价值设置为0。

7)从第二列开始，分别选取CandidateRidgePoints第j列每个候选脊点i，根据式(10)代价函数以及Can_Magnitudes和Can_Scales，依次求得列j-1各个候选脊点m到选取脊点i的代价函数值Cost(m)，并求最小代价函数值将Cost(i,j)放入代价值矩阵Costs对应位置中。

8)重复步骤7)到最后一列处理完毕，选取最后一列最小值的点作为脊点，并反推出所有脊点，得到小波脊。

9)重复步骤2)～8)，得到整个变形条纹图的小波脊，并根据式(6)求解出小波脊相位。

3 计算机实验

为验证本文所提算法的有效性和鲁棒性，本文使用两种不同形状的物体来进行模拟实验验证。并对变形条纹图分别添加方差为0.02、0.04、0.06、0.08、0.1的斑点噪声来验证所提方法的抗噪性。本文使用常用的评价指标-均方根误差RMSE(root mean square error)来计算算法在不同噪声强度的重建误差，RMSE的表达式为：

3.1 peaks模拟

模拟的被测物体由Matlab函数库中的peaks函数生成。物体的高度设为模拟条纹图大小为512×512 pixels，其中标准正弦光栅周期为16 pixel，模拟物体和在方差为0.1的噪声下变形条纹图分别如图2、图3所示。

图2 模拟物体

图3 含噪变形条纹图

图4 原代价函数法

图5 本文算法

图6 各算法在不同噪声下模拟10次RMSE曲线图

表1 peaks实验各算法模拟10次RMSE均值

首先，本文选取中值滤波技术对含噪的变形光栅图进行滤波处理，设定其滤波窗口为3*3。为了便于对比分析，分别采用了模极大值法、原代价函数法和本文算法对模拟变形条纹图进行小波脊相位提取。并对提取的截断相位都采用相同方法进行相位展开。选取的母小波为复Morlet小波，带宽参数为0.5，小波中心频率参数为1，尺度因子为1～64，步长为1。对于每种强度的噪声，分别使用3种方法模拟10次。在不同噪声情况下，各算法获得的10次RMSE 均值如表1所示，RMSE曲线如图6所示，原代价函数法和本文算法在方差为0.1的噪声下重建形状分别如图4、图5所示。

可以看出，随着噪声的增加，3种方法RMSE值都在增加，本文算法RMSE较小，明显优于模极大值法和原代价函数法，具有更好的抗噪能力。本文算法在噪声低的情况下明显优于其他两种方法，随着噪声的增加，模极大值法波动增大，误差极不稳定，原代价函数法和文本算法则波动较小，较为稳定。噪声增加后，本文算法与原代价函数法差距逐渐减小，但本文算法仍优于其他两种方法。

3.2 圆锥体模拟

为进一步验证本文算法的有效性和鲁棒性，使用高度变化较大的圆锥体作为验证对象，模拟的圆锥体的表达式为：

模拟条纹图大小为250×250 pixels，标准正弦光栅周期为16 pixel，模拟物体和在方差为0.1的噪声下变形条纹图分别如图7、图8所示。

图8 含噪变形条纹图

图9 原代价函数法

图10 本文算法

图11 各算法在不同噪声下模拟10次RMSE曲线图

表2 圆锥体实验各算法模拟10次RMSE均值

其他步骤和参数设置与peaks模拟相同，3种方法在不同噪声情况下模拟10次，均方根误差RMSE均值如表2所示，RMSE曲线如图11所示，原代价函数法和本文算法在方差为0.1的噪声下重建形状分别如图9、图10所示。

可以看出，本文算法明显优于其他两种方法。在较低噪声情况下，模极大值法波动较大，原代价函数法和本文算法波动较小，较为稳定，随着噪声的增加，原代价函数出现波动情况，本文算法仍相对稳定，这是由于圆锥物体高度变化较大，原代价函数提取的小波脊趋于过于平滑，造成小波脊提取不准确。

4 结束语

本文提出了一种新的代价函数提取小波脊的算法。相比原代价函数，本文首先使用Feature Scaling方法对小波变换系数幅值和尺度参数曲线梯度的权值进行调整，然后结合空域噪声信息，对尺度参数曲线梯度进行补偿。避免了在低噪声情况下提取小波脊时，脊线过于趋向平滑，以及在高噪声情况下，脊线波动较大造成小波脊选取不准确的问题。实验表明，无论对于高度变化平滑还是陡峭的物体，本文算法具有良好的抗噪性、鲁棒性以及稳健性。

本文研究工作得到重庆理工大学研究生创新基金(YCX2016226)的资助，在此表示感谢。