CN103886642A - 一种钢板表面三维重建快速实现方法 - Google Patents

Abstract

本发明涉及一种钢板表面三维重建的快速实现方法，采用红、绿、蓝三个单色光源沿不同角度同时照明钢板表面，其入射光与钢板轧制速度方向的角度在±50度范围，通过垂直于钢板表面的彩色CCD摄像机拍摄光源照射的钢板表面区域，分离其彩色图像的R、G、B通道，得到近似于红、绿、蓝光源单独照明下获得的三幅图像。本方法的复杂度与图像像素数量呈线性关系，而全局优化算法则呈平方关系，因此算法的耗时大大减少，且本算法的误差要略小于全局优化算法。

Description

一种钢板表面三维重建快速实现方法

技术领域

本发明涉及一种钢板表面三维重建的快速实现方法，该方法可用于钢板表面缺陷三维在线检测。采用红、绿、蓝三个单色光和彩色面阵CCD摄像机得到钢板表面的梯度，根据表面梯度重建出表面三维信息。本发明设计了从表面梯度重建出表面三维信息的快速算法，该算法复杂度与图像像素数量成线性关系，而全局优化算法的复杂度与图像数量成平方关系，并且该算法的误差要略小于全局优化算法。

背景技术

表面缺陷是影响钢板质量的一个重要因素，随着社会的发展，用户对板带产品的表面质量提出更高的要求，表面在线检测技术在钢板生产中的应用变得越来越广泛。基于CCD摄像的图像检测技术具有高速、高精度、适应性好等优势，已经成为表面在线检测研究和应用的主流。但是目前在线应用的钢板表面在线检测系统，其检测原理大多基于二维灰度图像检测技术，即通过灰度图像提供的几何形状、明暗分布等二维信息来检测与识别钢板的表面缺陷。但是，由于钢板表面缺陷非常复杂，多数缺陷的几何形状不固定，如压痕、辊印等缺陷，通过几何形状等二维特征很难对缺陷进行正确分类。此外，一些氧化铁皮、油渍、水迹等“伪缺陷”与缺陷在形状、灰度上非常相似，容易识别成缺陷，造成大量的“误检”。

通过三维测量方法可获取表面深度信息，由于光度立体法具有分辨率高、受环境影响小、结构简单的优点，容易在钢板生产线上安装应用。采用光度立体法可以得到表面的三维形貌，通过三维形貌检测压痕、划伤、辊印、裂纹、夹杂等缺陷。由于氧化铁皮、油渍、水迹等“伪缺陷”在表面深度上几乎没有变化，因此通过表面三维形貌检测缺陷可避免“伪缺陷”的干扰，减少“误检”。

光度立体法的理论基础是朗伯漫反射模型，通过单台摄像机拍摄物体在不同方向光源单独照明下的一组图像，根据表面反射模型可计算出表面法向量和梯度分布。金属表面反射模型较复杂，同时含有镜面反射和漫反射成分，轧制加工的钢板表面具有均匀的纹理，其微观结构沿轧制方向排列。入射光接近轧制纹理方向时，反射光含有较多的漫反射成分；接近轧制纹理垂直方向时则含所有较多镜面反射成分。将入射光限制在轧制纹理方向±50度范围内时可以近似视为理想漫反射表面。

通过光度立体法可以得到表面法向量或表面梯度，然后通过二维积分或优化拟合算法根据梯度信息重建出表面三维信息，这一过程称为三维重建。三维重建算法是表面三维测量的关键，传统的三维重建算法有两种：全局优化算法和局部积分算法。全局优化算法的复杂度与图像像素数量成平方关系，随着图像像素数量的增加，算法耗时急剧增加，不适用于在线测量。局部积分算法的误差太大，不具有实用的价值。

为了将表面三维重建方法应用于钢板表面缺陷在线检测，本发明设计了一种基于Haar小波的快速三维重建算法，其算法复杂度与图像像素数量成线性关系，远小于全局优化算法的复杂度，而该算法的误差要略小于全局优化算法。

发明内容

由于钢板在生产线上处于高速运动状态，为了解决动态物体的测量问题，本发明设计了动态光度立体视觉检测方案。

本发明采用红、绿、蓝三个单色光源沿不同角度同时照明钢板表面同一区域，其入射光与钢板轧制速度方向的角度在±50度范围，通过垂直于钢板表面的彩色CCD摄像机拍摄光源照射的钢板表面区域，分离其彩色图像的R、G、B通道，得到近似于红、绿、蓝光源单独照明下获得的三幅图像I_R,I_G,I_B。I_R,I_G,I_B的高h=2ⁿ,宽w=2ⁿ，图像坐标系(x,y)中某点的相对深度为Z(x,y)，梯度矩阵P为深度矩阵Z沿x方向的差分，Q为深度矩阵Z沿y方向的差分，P、Q由I_R,I_G,I_B和光源方向矩阵L得到，P(x,y),Q(x,y)的初始值用P¹,Q¹表示，深度矩阵Z当作小波分解的初始矩阵LL⁰，则Z可以通过式(1)所示的小波分解算法，式(2)所示的递推公式和式(3)所示的小波重构算法得到：

$\left\{\begin{array}{c}{\mathrm{LH}}^k=D(Q^k*L_D)\\ {\mathrm{HL}}^k=D(P^k*{L_D}^T)\\ {\mathrm{HH}}^k=D(P^k*{H_D}^T)=D(Q^k*H_D)\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{c}{P}^k=D(P^{k-1}*L_D*{L_D}^T*L_D)\\ Q^k=D(Q^{k-1}*{L_D}^T*L_D*{L_D}^T)\end{array}\right.---\left(2\right)$

LL^k-1＝U(LL^k)*L_R ^T*L_R+U(LH^k)*H_R ^T*L_R+

U(HL^k)*L_R ^T*H_R+U(HH^k)*H_R ^T*H_R  (3)

式(1)、式(2)、式(3)中：

■“*”——二维卷积运算；

■D(M)——二维下采样，抽取矩阵奇数行、列，D(M)_(x,y)=M_(2x-1,2y-1)；

■U(M)——二维上采样，矩阵扩充为(2h+1)×(2w+1)，其中偶数行列：U(M)_(2x,2y)=M_(x,y)，其余项填充0；

■LL⁰——小波分解的初始矩阵，Z=LL⁰；

■LL^k、LH^k、HL^k、HH^k——第k层分解的低频及高频子带(1≤k≤n)，设LLⁿ中元素的值为0；

■P^k,Q^k——k层梯度矩阵(1≤k≤n)，其中P¹,Q¹为P,Q，其余通过递推求得；

■L_D=(1,1),H_D=(1,-1)——小波分解低通、高通滤波器；

■L_R=(0.5,0.5),H_R=(-0.5,0.5)——小波重构低通、高通滤波器。

本发明的另一个技术方案是上述的LL^k(1≤k<n)可以通过式(4)进行优化：

$\begin{array}{c}{\mathrm{LL}}_m^k(x,y)=[{\mathrm{LL}}_{m-1}^k(x-1,y)+P^k(x-1,y)]/4+\\ [{\mathrm{LL}}_{m-1}^k(x,y-1)+Q^k(x,y-1)]/4+\\ [{\mathrm{LL}}_{m-1}^k(x+1,y)-Q^k(x+1,y)]/4+\\ [{\mathrm{LL}}_{m-1}^k(x,y+1)-P^k(x,y+1)]/4\end{array}---\left(4\right)$

式(4)中，P^k,Q^k的定义见权利要求1，m为迭代次数。

本发明是这样实现的：

将深度矩阵Z当作小波分解的初始矩阵LL⁰，用LL^k、LH^k、HL^k、HH^k(0＜k≤n)表示Z经第k层二维Haar小波分解后的尺度分量、垂直分量、水平分量和对角分量。根据小波重构算法，只要知道最末层LLⁿ的值以及LH^k、HL^k、HH^k(0＜k≤n)的值，就可以依次重构得到LL^k(0≤k＜n)，LL⁰即为深度矩阵Z。因此本发明根据P、Q与二维Haar小波分解之间的关系，通过小波分解、递推算法和小波重构等步骤，得到通过P、Q求Z的快速算法。

进一步说，由于P、Q分别为深度矩阵Z沿x、y方向的差分，将P、Q代入二维Haar小波分解，可得LH¹、HL¹、HH¹。根据P、Q的递推公式得到P²,P³,…,Pⁿ和Q¹,Q³,…,Qⁿ，并依次求出LH^k、HL^k、HH^k(1＜k≤n)。设LLⁿ的值为0，由于LH^k、HL^k、HH^k(0＜k≤n)已知，因此根据二维Haar小波重构算法可依次重构得到LL^k(0≤k＜n)，LL⁰即为深度矩阵Z。

在重构LL^k(0≤k＜n)时会产生局部误差，可以用传统的全局优化算法对LL^k进行迭代优化，以减少误差。

本发明的有益效果是：基于Haar小波的快速三维重建算法的复杂度与图像像素数量呈线性关系，而全局优化算法则呈平方关系，因此耗时大大减少，且误差要略小于全局优化算法。

附图说明

图1为光源与摄像机的布置方式。图1中：1为待检测的物体，2为彩色CCD摄像机，3为红色光源，4为绿色光源，5为蓝色光源。

图2为基于Harr小波的三维重建算法流程图。

图3为误差修正后的算法流程。

具体实施方式

图1中，摄像机2放置在垂直于钢板1的表面方向上。红色光源3、绿色光源4、蓝色光源5围绕摄像机2布置在不同的位置，3台光源照射的区域与摄像机的采集区域重合。摄像机2为彩色CCD摄像机，其采集到的图像可分离出红（R）、绿（G）、蓝（B）三个通道图像，分别对应光源3、光源4、光源5的反射光图像。

根据光度立体学（Photometric Stereo）原理，对于发光强度为E的光源，其反射光的光强为：

I＝ρEl·n  (5)

式(5)中，n＝(n_x n_y n_z)^T，为表面上某点处的单位法向量；ρ为漫反射系数；l＝(l_x l_y l_z)，为光源单位方向向量；E为入射光强。

通过调试使光源3、光源4、光源5在摄像机2中的感光亮度相等，即ρE＝ρ_RE_R＝ρ_GE_G＝ρ_BE_B。则根据式(5)有：

$\left\{\begin{array}{c}{I}_R={\mathrm{\&rho;El}}_R\&CenterDot;n\\ I_G={\mathrm{\&rho;El}}_G\&CenterDot;n\\ I_B={\mathrm{\&rho;El}}_B\&CenterDot;n\end{array}\right.---\left(6\right)$

式(6)中，I_R,I_G,I_B分别为R、G、B通道的图像亮度，l_R,l_G,l_B分别为红、绿、蓝光源入射方向向量。

设光源方向矩阵为 $L=\left(\begin{array}{ccc}{l}_{\mathrm{Rx}}& l_{\mathrm{Ry}}& l_{\mathrm{Rz}}\\ l_{\mathrm{Gx}}& l_{\mathrm{Gy}}& l_{\mathrm{Gz}}\\ l_{\mathrm{Bx}}& l_{\mathrm{By}}& l_{\mathrm{Bz}}\end{array}\right),$ 反射光强为：I_RGB＝(I_R I_G I_B)^T，则表面法向量为：

n＝L^-1·I_RGB/|L^-1·I_RGB|  (7)

假设图像坐标系中某点的相对深度为Z(x,y)，定义其沿图像x方向和y方向的梯度分别为P(x,y),Q(x,y)：

$\left\{\begin{array}{c}P(x,y)={\mathrm{\&Delta;Z}}_x(x,y)\\ Q(x,y)={\mathrm{\&Delta;Z}}_y(x,y)\end{array}\right.---\left(8\right)$

整合图像范围内的单位法向量n为三个矩阵：N_x、N_y、N_z，则有：

(P,Q)＝(N_x·/N_z,N_y·/N_z)  (9)

根据式(7)和式(9)可以得到P、Q。传统从梯度矩阵P、Q重构三维表面的算法主要分为两种：局部积分算法和全局优化算法。局部积分算法整体误差较大，而全局优化算法运算时间较长，不适于在线检测。为了提高重建算法的速度，设计了一种基于Haar小波的快速三维重建算法。假定原始图像的高h=2ⁿ,宽w=2ⁿ,图像坐标范围：

1≤x≤w,1≤y≤h。定义以下的符号和运算：

■“*”——二维卷积运算；

■D(M)——二维下采样，抽取矩阵奇数行、列，D(M)_(x,y)=M_(2x-1,2y-1)；

■U(M)——二维上采样，矩阵扩充为(2h+1)×(2w+1)，其中偶数行列：U(M)_(2x,2y)=M_(x,y)，其余项填充0；

■LL⁰——小波分解的初始矩阵；

■LL^k、LH^k、HL^k、HH^k——第k层分解的低频及高频子带(1≤k≤n)；

■P^k,Q^k——k层梯度矩阵(1≤k≤n)，其中P¹,Q¹为P,Q，其余通过递推求得；

■L_D=(1,1),H_D=(1,-1)——小波分解低通、高通滤波器；

L_R=(0.5,0.5),H_R=(-0.5,0.5)——小波重构低通、高通滤波器。

传统二维小波分解算法如下：

$\left\{\begin{array}{c}{\mathrm{LL}}^k=D({\mathrm{LL}}^{k-1}*L_D*{L_D}^T)\\ {\mathrm{LH}}^k=D({\mathrm{LL}}^{k-1}*L_D*{H_D}^T)\\ {\mathrm{HL}}^k=D({\mathrm{LL}}^{k-1}*H_D*{L_D}^T)\\ {\mathrm{HH}}^k=D({\mathrm{LL}}^{k-1}*H_D*{H_D}^T)\end{array}\right.---\left(10\right)$

式(10)中，0<k≤n，n为最大分解层数。二维小波重构公式为：

LL^k-1＝U(LL^k)*L_R ^T*L_R+U(LH^k)*H_R ^T*L_R+

U(HL^k)*L_R ^T*H_R+U(HH^k)*H_R ^T*H_R  (11)

图2为基于Haar小波的快速三维重建算法的流程，将深度矩阵Z当作小波分解的初始矩阵LL⁰。该算法的核心是分解算法、递推算法和重构算法。

（1）分解算法

由式(8)可知，梯度矩阵P,Q分别为深度矩阵Z（即LL⁰）沿x,y方向的差分，式(8)同样适用于P^k,Q^k，因此有：

$\left\{\begin{array}{c}{P}^k={{\mathrm{\&Delta;LL}}^{k-1}}_x={\mathrm{LL}}^{k-1}*H_D\\ Q^k={{\mathrm{\&Delta;LL}}^{k-1}}_y={\mathrm{LL}}^{k-1}*{H_D}^T\end{array}\right.---\left(12\right)$

式(12)中，0<k≤n，将式(12)代入式(10)中，得到P^k,Q^k与高频子带LH^k、HL^k、HH^k的关系：

$\left\{\begin{array}{c}{\mathrm{LH}}^k=D(Q^k*L_D)\\ {\mathrm{HL}}^k=D(P^k*{L_D}^T)\\ {\mathrm{HH}}^k=D(P^k*{H_D}^T)=D(Q^k*H_D)\end{array}\right.---\left(13\right)$

（2）递推算法

由式(10)、式(12)可得：

P^k＝D(LL^k-2*L_D*L_D ^T)*H_D  (14)

已知H_D=(1,-1)，H_D*L_D=(1,0,-1)，以H_D与H_D*L_D为桥梁，可建立起下采样运算D(X)括号内外的运算关系：

D(X)*H_D＝D(X*H_D*L_D)  (15)

将式(15)应用于式(14)，并将式(12)代入，得到关于P的递推关系：

P^k＝D(P^k-1*L_D*L_D ^T*L_D)  (16)

同理，Q的递推关系为：

Q^k＝D(Q^k-1*L_D ^T*L_D*L_D ^T)  (17)

式(16)、式(17)中，1≤k<n。通过式(12)、式(13)可递推得到P²～Pⁿ和Q²～Qⁿ，根据式(13)可得LH²～LHⁿ、HL²～HLⁿ、HH²～HHⁿ。

（3）重构算法

设最终层LLⁿ=0，由式(11)递推可得LL^n-1,…,LL¹,LL⁰，其中LL⁰就是深度矩阵Z。

（4）修正算法

如图3所示，在式(1)的每一步重构过程后，加入一个迭代优化过程，以修正其局部误差。即对低频子带LL^k(0≤k<n,w=h=2ⁿ)使用以下的全局优化公式进行迭代：

$\begin{array}{c}{\mathrm{LL}}_m^k(x,y)=[{\mathrm{LL}}_{m-1}^k(x-1,y)+P^k(x-1,y)]/4+\\ [{\mathrm{LL}}_{m-1}^k(x,y-1)+Q^k(x,y-1)]/4+\end{array}$

$\begin{array}{c}[{\mathrm{LL}}_{m-1}^k(x+1,y)-Q^k(x+1,y)]/4+\\ [{\mathrm{LL}}_{m-1}^k(x,y+1)-P^k(x,y+1)]/4\end{array}---\left(18\right)$

式(18)中，m为迭代次数。

表1为全局优化算法和本算法的时间复杂度和误差统计。

表1算法时间复杂度与误差统计

(注：运行环境：主频2.6GHz，Windows7操作系统，Matlab计算软件；误差函数为：

$\mathrm{\&epsiv;}=\frac{1}{w*h}\mathrm{\&Sigma;}{({\mathrm{\&Delta;Z}}_x(x,y)-P(x,y))}^2+{({\mathrm{\&Delta;Z}}_y(x,y)-Q(x,y))}^2)$

由表1可见，本算法误差要略小于全局优化算法。在算法时间复杂度上，全局优化算法与像素总数呈平方关系：O(n²)…n=w=h；而本算法则呈线性关系：O(n)…n=w=h。

Claims (2)

1.一种钢板表面三维重建的快速实现方法，其特征在于：采用红、绿、蓝三个单色光源沿不同角度同时照明钢板表面同一区域，其入射光与钢板轧制速度方向的角度在±50度范围，通过垂直于钢板表面的彩色CCD摄像机拍摄光源照射的钢板表面区域，分离其彩色图像的R、G、B通道，得到近似于红、绿、蓝光源单独照明下获得的三幅图像I_R,I_G,I_B。I_R,I_G,I_B的高h=2ⁿ,宽w=2ⁿ，图像坐标系(x,y)中某点的相对深度为Z(x,y)，梯度矩阵P为深度矩阵Z沿x方向的差分，Q为深度矩阵Z沿y方向的差分，P、Q由I_R,I_G,I_B和光源方向矩阵L得到，P(x,y),Q(x,y)的初始值用P¹,Q¹表示，深度矩阵Z当作小波分解的初始矩阵LL⁰，则Z可以通过式(1)所示的小波分解算法，式(2)所示的递推公式和式(3)所示的小波重构算法得到：

$\left\{\begin{array}{c}{\mathrm{LH}}^k=D(Q^k*L_D)\\ {\mathrm{HL}}^k=D(P^k*{L_D}^T)\\ {\mathrm{HH}}^k=D(P^k*{H_D}^T)=D(Q^k*H_D)\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{c}{P}^k=D(P^{k-1}*L_D*{L_D}^T*L_D)\\ Q^k=D(Q^{k-1}*{L_D}^T*L_D*{L_D}^T)\end{array}\right.---\left(2\right)$

LL^k-1＝U(LL^k)*L_R ^T*L_R+U(LH^k)*H_R ^T*L_R+

U(HL^k)*L_R ^T*H_R+U(HH^k)*H_R ^T*H_R  (3)

式(1)、式(2)、式(3)中：

■“*”——二维卷积运算；

■D(M)——二维下采样，抽取矩阵奇数行、列，D(M)_(x,y)=M_(2x-1,2y-1)；

■U(M)——二维上采样，矩阵扩充为(2h+1)×(2w+1)，其中偶数行列：U(M)_(2x,2y)=M_(x,y)，其余项填充0；

■LL⁰——小波分解的初始矩阵，Z=LL⁰；

■LL^k、LH^k、HL^k、HH^k——第k层分解的低频及高频子带(1≤k≤n)，设LLⁿ中元素的值为0；

■P^k,Q^k——k层梯度矩阵(1≤k≤n)，其中P¹,Q¹为P,Q，其余通过递推求得；

■L_D=(1,1),H_D=(1,-1)——小波分解低通、高通滤波器；

■L_R=(0.5,0.5),H_R=(-0.5,0.5)——小波重构低通、高通滤波器。

2.一种钢板表面三维重建的快速实现方法，其特征在于：所述的LL^k可以通过式(4)进行优化：

$\begin{array}{c}{\mathrm{LL}}_m^k(x,y)=[{\mathrm{LL}}_{m-1}^k(x-1,y)+P^k(x-1,y)]/4+\\ [{\mathrm{LL}}_{m-1}^k(x,y-1)+Q^k(x,y-1)]/4+\\ [{\mathrm{LL}}_{m-1}^k(x+1,y)-Q^k(x+1,y)]/4+\\ [{\mathrm{LL}}_{m-1}^k(x,y+1)-P^k(x,y+1)]/4\end{array}---\left(4\right)$

式(4)中，P^k,Q^k的定义见权利要求1，m为迭代次数。

Images