1 引言

在基于激光三角测量原理的线结构光测量系统中,激光器发射出的一字线结构光投射至被测物体(如工件、印刷电路板等)表面,经物体表面高度信息调制后的激光光带含有物体的三维轮廓信息。通过分析光带图像可以获得被测物体的三维(3D)几何尺寸^[1-4],该方法可用于质量检验、逆向工程、产品建模等领域中。这种结构光视觉测量技术具有非接触、快速、灵活性强等优势,已发展成为最为有效的3D测量技术之一^[5]。

在工程实际中,有一类细长型、板条状产品,如雪糕棒、木地板等,其材质较软、外形结构特殊,容易发生弯曲变形。其中有一类弯曲在工程上称作死弯,即局部弯曲。局部弯曲的特点是在产品表面出现局部小面积的凹陷或凸起,较大区域仍为平面。对这类产品进行质量检测时,对局部弯曲程度有着严格的规定。含有局部弯曲信息的光带图像,绝大部分光带均保持直线特征,只是局部区域出现一定程度的弯曲变形。对于这类含有局部弯曲光带图像,一般可以通过传统直线检测技术进行分析处理,如光带中心提取^[6-7]结合最小二乘拟合法、Hough变换^[8]、Radon变换^[9]、边缘跟踪算法^[10]等。目前,这些处理算法一般需要对目标图像进行阈值分割、轮廓提取或其他预处理操作,处理过程繁琐,当图像信噪比偏低、光带强度不均匀、光带强度较低或光带宽度较宽时,它们的适应能力和检测精度将明显降低。

光带图像边缘沿法线方向的灰度具有最大的变化率,即灰度梯度最大,在频域中主要表现为高频成分,因此幅度谱能量分布总会沿光带的法线方向伸展。傅里叶变换具有平移不变性,光带目标位于图像中任意位置时,其幅度谱始终以频谱图像的中心(即零频位置)为原点。极坐标变换可将笛卡儿坐标下的旋转成分转化为平移成分。本文将这两者结合起来,对光带图像进行傅里叶变换,再以幅度谱中心为原点进行极坐标变换,由此得到光带目标的法线方向,进而对光带图像在法线方向上进行灰度投影,得到光带目标的局部弯曲量。

文献[ 11-12]提到的算法是对经极坐标变换后的空域图像进行傅里叶变换,严格意义上应称为极坐标-傅里叶变换,主要用到旋转不变特征,用于两幅图像间的配准或图像检索。本文采用傅里叶-极坐标变换算法检测光带的局部弯曲,该方法无需对图像进行滤波、阈值分割等预处理操作,利用图像频域的整体信息,对噪声、降质图像具有很好的适应性,检测步骤简捷,检测结果更为准确。

2 傅里叶-极坐标变换

2.1 二维傅里叶变换的平移与旋转特性

假设图像f₂(x,y)相对于图像f₁(x,y)而言,存在平移与旋转变换,表示为

f2(x,y)=f1(xcosϕ+ysinϕ-x0,-xsinϕ+ycosϕ-y0),(1)

式中x、y为空域坐标,x₀、y₀为平移量坐标,ϕ为旋转角度。

设F₁(u,v)与F₂(u,v)分别为f₁(x,y)与f₂(x,y)的傅里叶变换,二者之间的关系为

F2(u,v)=F1(ucosϕ+vsinϕ,-usinϕ+vcosϕ)·exp[-j2π(x0u+y0v)],(2)

式中u、v为频域坐标。

设M₁(u,v)与M₂(u,v)分别为F₁(u,v)与F₂(u,v)的幅度谱,则有

M2(u,v)=M1(ucosϕ+vsinϕ,-usinϕ+vcosϕ)。(3)

(3)式显示,幅度谱与平移(x₀,y₀)无关,且

图 1. 二维傅里叶变换的平移和旋转特性。(a) f1; (b) f2; (c) M1; (d) M2

Fig. 1. Translation and rotation characteristics of 2D Fourier transformation. (a) f1; (b) f2; (c) M1; (d) M2

下载图片 查看所有图片

与空间域图像具有相同的旋转角度ϕ。

二维傅里叶变换的平移与旋转特性可用图1所示的空间域图像(原始图像f₁,含有平移和旋转变换的图像f₂)及其对应的幅度谱来直观反映。可以看出,幅度谱M₁和M₂并不依赖于目标的平移量,因为空域中的平移仅影响傅里叶变换的相位值。其次,空域图像旋转某一角度ϕ(给定角度为-10.2°)时,其幅度谱相应地旋转相同的角度。

幅度谱的能量分布总是沿目标图像灰度变化的最大梯度方向,对于图1所示的两个直线目标,其幅度谱能量的分布与直线目标的法线方向相对应。空域目标位于图像中任何位置时,其幅度谱能量分布总是以频谱中心(零频)为原点。故可通过估计幅度谱中能量沿角度的分布情况来估计空域中含有类似直线边缘目标的法线方向或目标伸展方向。

图 2. 用于估计幅度谱中能量随角度变化关系时的扇形区域

Fig. 2. Fan-shaped area used to evaluate relationship between energy and angle in magnitude spectrum

下载图片 查看所有图片

计算幅度谱M中能量随角度的变化关系时,以幅度谱中心为原点O,如图2所示,对每一小扇形区域Δα内的能量值进行累加,得到区域内总能量与角度之间的关系曲线h(ϕ),曲线中峰值对应的角度即为目标的法线方向。通过沿轮辐线积分可得到角度为ϕ_k时扇形区内能量累积值为

h(ϕk)=∫ϕk-Δα/2ϕk+Δα/2∫0ρM(r,ϕ)drdϕ,(4)

式中r=u2+v2,ρ为最大半径。ϕ_k变化一周后可得到能量与角度间的关系曲线h(ϕ)。

2.2 极坐标变换

事实上,幅度谱为离散信号,(4)式中的h(ϕ_k)可通过对角度满足(ϕ_k-Δα/2)<ϕ<(ϕ_k+Δα/2)的矩形网格里的幅度值进行累加得到。将笛卡儿坐标系下的幅度谱转化至极坐标系中再计算h(ϕ),过程更为有效便捷。

笛卡儿坐标系中幅度谱的某一像素(u,v),表示为极坐标系中(r,ϕ),转化公式^[13-14]为

r=u2+v2ϕ=tan-1(v/u),(5)

式中r表示笛卡儿坐标系下某一像素(u,v)对应矢量(相对于原点)的幅值;ϕ为该矢量的幅角,取逆时针方向为正。

极坐标变换的示意图如图3所示,笛卡儿坐标系中的一点映射为极坐标系中对应点的关系式由(5)式确定。根据极坐标变换原理可知,笛卡儿坐标系中以原点为圆心的圆环线映射为极坐标系中平行于ϕ轴的横向线;笛卡儿坐标系中过原点的轮辐线映射为极坐标系中平行于r轴的纵向线,故极坐标系下沿r方向进行累加便可得到h(ϕ)。

图 3. 极坐标变换示意图。(a)笛卡儿坐标系; (b)极坐标系

Fig. 3. Schematic of polar transformation. (a) Cartesian coordinate; (b) polar coordinate

下载图片 查看所有图片

为了得到更准确的h(ϕ),对幅度谱进行极坐标变换时,仅分析图4所示的环形区域,即ερ'<r<ρ',其中ρ'为幅度谱图像的半径或图像尺寸(长或宽的最小值)的一半。一般情况下低频值很大,在实际处理时选择系数ε=0.1,因此半径范围为0.1ρ'<r<ρ'。由于二维傅里叶变换具有原点对称性,故只需针对幅度谱上半部分进行极坐标变换即可计算幅度谱的能量分布方向,故角度范围为0°≤ϕ<180°。

图 4. 极坐标变换时的环形区域

Fig. 4. Ring-shaped area used for polar transformation

下载图片 查看所有图片

图1所示的幅度谱M₁和M₂经极坐标变换后的结果如图5所示,可以看出,幅度谱中相对于原点的旋转成分被转化为极坐标系中沿ϕ方向的平移。

图 5. 图1所示幅度谱的极坐标变换。(a) M1; (b) M2

Fig. 5. Polar transformation of magnitude spectra shown in Fig. 1. (a) M1; (b) M2

下载图片 查看所有图片

2.3 光带法线方向的计算

通过傅里叶-极坐标变换,可将(4)式表示的幅度谱扇形区域的积分转化为极坐标系中沿r方向的累加求和,即

h(ϕk)=∑r=0.1ρ'ρ'M(r,ϕk)。(6)

当ϕ_k在0°~180°范围内变化时计算出h(ϕ)。

图 6. 能量与角度间的函数关系

Fig. 6. Energy as a function of angle

下载图片 查看所有图片

根据图5所示的极坐标变换结果计算出的h(ϕ)曲线如图6所示。两条曲线的峰值位置分别为ϕ₁=90.87°,ϕ₂=101.16°,该角度即为图1所示两幅图像中光带目标的法线方向。直线目标的伸展方向为θ=90°-ϕ,其中θ表示直线目标与x轴正方向之间的夹角,逆时针为正。因此两幅图像中光带目标的伸展角分别为θ₁=-0.87°、θ₂=-11.16°,二者之差为-10.29°,与给定角度相符。

3 基于方向投影的光带位置检测

一幅图像在给定角度φ上的方向投影表示为

R(x')=∫f(x'cosφ-y'sinφ,x'sinφ+y'cosφ)dy',(7)

式中φ为投影坐标轴x'与图像坐标x轴的夹角。投影坐标(x',y')与原图像坐标之间的变换关系为

x'y'=cosφsinφ-sinφcosφxy。(8)

φ=90°、0°分别表示水平、垂直投影。对于给定角度φ,(7)式表示沿这一方向上的线积分或灰度累加^[15],如图7所示,若沿矩形目标长边缘的法线方向进行投影,由此计算出的方向投影曲线R(x')中方波信号的边界将对应于矩形目标的长边缘轮廓。

图 7. 方向投影示意图

Fig. 7. Schematic of directional projection

下载图片 查看所有图片

根据计算出的法线角度(ϕ₁=90.87°,ϕ₂=101.16°),对图1中两幅图像f₁与f₂分别进行方向投影,结果如图8所示,曲线脉冲相对于峰值左右对称、接近高斯分布,脉冲宽度均为14 pixel,该数值对应于光带宽度,说明光带目标近似为直线,不存在局部弯曲起伏。图8中横坐标x'的原点对应于原图像(f₁与f₂)的中心位置,曲线峰值位置相当于光带中心线在图像坐标y轴上的截距。通过上述分析,可获得光带中心线的拟合直线参数(斜率与截距)。

图 8. 图1中图像在目标法线方向上的投影曲线。(a) f1; (b) f2

Fig. 8. Projection curves along normal direction of images shown in Fig.1. (a) f1; (b) f2

下载图片 查看所有图片

4 光带图像局部弯曲的工程测量实验

针对长度为114 mm的细长平面木制品建立了形状误差检测的激光三角测量系统,激光器发出线结构光投射在产品表面。若表面出现局部凸起或凹陷时,光带图像将发生变形。

实验中采集到的两幅大小为1250 pixel×280 pixel的光带图像如图9(a)、(b)所示。利用傅里叶-极坐标变换算法得到幅度谱能量与角度的关系曲线h(ϕ)如图9(c)、(d)所示,计算出光带法线方向角分别为88.18°和89.97°。实验图像在各自法线方向上的投影结果如图9(e)、(f)所示,曲线上除了明显的峰值外,还有一个距离峰值位置最远的拐点,如图9中小矩形框所示,二者位置之差即为局部弯曲量。由此得到两幅实验图像的局部弯曲量分别为30 pixel和24 pixel。

图 9. 光带图像局部弯曲的实验检测结果。 (a)(b) 实验采集的光带图像; (c)(d)经傅里叶-极坐标变换后得到的h(?); (e)(f) 光带图像在法线方向上的投影曲线

Fig. 9. Experimental detecting results of local bending of laser stripe images. (a) (b) laser stripe images recorded in experiment; (c) (d) h(?) obtained after Fourier-polar transformation; (e) (f) projection curves of laser stripe images along normal direction

下载图片 查看所有图片

5 方法的性能测试及讨论

为了进一步验证所提方法的性能,在图9(b)对应的被测木制品表面人为添加局部污染,使物体表面的反射率发生较大起伏变化,此时光带目标将出现宽度不均匀与强度不均匀特征。同时,通过数值方法在该光带图像中加入均值为0.1、方差为0.005的高斯随机噪声,实验图像如图10(a)所示,其中两个矩形区域即为污染区域的光带。经傅里叶-极坐标变换后的结果及h(ϕ)曲线如图10(b)所示,得到光带法线方向角为89.97°,与由图9(d)计算出的角度值相等。实验图像在该法线方向上投影后的结果如图10(c)所示,投影前通过将图像各像素的灰度值减去背景灰度后可以消除投影曲线中R(x')值的整体偏移。比较图10(c)与图9(f)发现,两条曲线相似度极高,且由图10(c)计算出的局部弯曲量与图9(f)的计算结果吻合。

图 10. 含有局部污染并加入高斯噪声的光带图像的实验结果。(a)光带图像; (b) 傅里叶-极坐标变换结果与h(?); (c) 光带图像在法线方向上的投影曲线

Fig. 10. Experimental results of laser stripe images with local contamination and Gaussian noises. (a) Laser stripe image; (b) Fourier-polar transformation result and h(?); (c) projection curve of laser stripe image along normal direction

下载图片 查看所有图片

图 11. 灰度重心法对图10(a)中心线的提取结果

Fig. 11. Extracting result of center-line in Fig. 10(a) by gray level gravity algorithm

下载图片 查看所有图片

传统方法一般需要提取光带中心线,采用直线拟合算法(如最小二乘法或者Hough变换法)得到直线参数,进而计算拟合直线与中心线上最远点之间的距离,由此得到弯曲量。光带中心线提取与直线拟合过程均会引入计算误差,尤其是光带目标偏低的信噪比会对中心线提取算法的稳健性提出更高要求。图11所示为利用文献[ 16]提出的灰度重心法对图10(a)进行中心线提取后的结果,在光带亮度偏暗区,中心坐标提取错误的概率显著上升,由此计算出的局部弯曲量也是不准确的,甚至是错误值。光带中心提取时还应考虑光带的方向,否则在光带曲率变化较大位置处必将产生较大的计算误差^[17],如图11放大显示区域所示。而阈值化的检测算法往往需要阈值分割、骨架化处理等预处理操作,过程繁琐且精度不可控。

上述分析表明,当光带目标的宽度及强度出现不均匀分布、图像中存在较大随机噪声时,所提方法仍能实现准确检测;而一些传统方法既不能直接获得被测结果,还容易受光带图像质量的影响。

6 结论

利用基于一字单线结构光的三角测量技术测量平板的局部起伏时,关键技术是对采集到的光带图像进行局部弯曲检测,重点对光带图像的处理方法进行了研究。采用傅里叶-极坐标变换算法估计光带目标的法线方向,对光带图像进行了灰度方向投影以实现光带目标检测,由投影曲线直接计算得到了光带局部弯曲量。检测算法只包括傅里叶变换、极坐标变换与灰度投影三个处理过程,无需对原始图像进行预处理操作,无需人工参与,算法简捷。通过三个实际光带图像的检测实验,证实了该方法应用于光带局部弯曲检测及直线类目标检测的有效性,该方法对随机噪声、光带强度不一致分布等情况具有较好的免疫力。