2.4 图像边缘及其检测
边缘检测是底层计算机视觉中最重要的问题之一,但到目前为止,如何有效地从真实场景中进行边缘检测仍然非常具有挑战性。底层好的边缘检测结果对后续高层的计算机视觉任务的成功具有非常重要的作用。在大量的原始图像数据中,如何有效地获取计算机视觉任务所需要的特征(如边缘、角点等)并去除无关的背景和各种干扰,成为底层计算机视觉图像处理技术的核心问题。
例如,在立体视觉技术中,一种常见的操作就是需要获取图像的深度信息(Depth Information),如图2.17所示。该操作一般采用双目视觉的基本原理来进行,如图2.18所示,采用来自不同位置的图像来获取相关的深度信息。而这正好与人类双眼获取图像后,可以感知一定的深度信息相对应。在视觉技术中,这种信息的获取都需要做左右图像之间的匹配对应,即如图2.17的原始图像所示,求出两幅图像中的匹配点,然后通过视差来确定该点的深度信息。而做这种匹配对应(配准),图像的一些特征信息起着决定性作用。此外,在视觉跟踪中,运动目标的检测和关联实际上也是采用目标的表观特征并结合时序信息来进行处理的。来自运动的结构恢复技术、来自形状的结构恢复技术和遮挡的对象识别与推理等,都利用了对象的边缘、角点等具有区分能力的特征来进行处理。因此,底层图像处理结果对视觉任务至关重要。
图2.17 BookArrival序列的纹理图(左)和其对应的深度图(右)
图2.18 基于双目视觉的深度测量原理
ITU与MPEG成立的联合视频专家组在3DAV编码标准的研究与制定中,将立体视频数据定义为纹理加深度的表示格式。其中,纹理图是采集的一个视点图像;深度图是对纹理图中每个区域相对深度的描述。在多视点视频编码标准中,多视点视频由少数几个纹理视频联合深度视频表示,深度视频图像中的像素值表示与其对应纹理图像中像素点位置的目标与相机的相对距离,像素值按照深度的范围(Zmin,Zmax)映射,在多视点视频编码标准中,深度取值按照8bit采样,即将深度范围分为256个均匀间隔,距离摄像机最远的Zmax取值为0,最近的Zmin取值为255。这里的图像由摄像机(视点为12)采集的BookArrival序列的纹理图和其对应的深度图来表示。
图2.18中P为物体上的一个点,其在左右两个相机成像平面上的像点分别为Pl、Pr、Cl、Cr,它们分别为左右两个相机的光学中心位置。注意,图中左右两个光学相机像平面在同一个平面上,这里希望求出深度d与相机位置及焦距之间的关系,f为焦距,la-lb为点P在左右两个相机的像平面上的视差。由三角形的相似关系△PPrB≈△PCrD可知,得?,并且可得(d−f)(a+lb)=ad,即;由△PPlB≈△PClD可知,得,从而有,计算出,从而可得。因此,一旦对齐相机后,只需求出对应点在两幅图像中的位置,然后求出其视差,就可以估计出物点的深度信息。
2.4.1 边缘类型
从概念上来看,边缘一般指像素值发生突变的区域,但由于数字图像成像过程中无法形成理想的突变状态,因此一般在图像中很少出现突变的边缘,而是有一个渐变的过程。一般在图像处理中,有4种边缘类型,当然这4种类型也可以上下翻转,分别对应阶梯状、斜坡状、脉冲状和屋顶状4种边缘类型,如图2.19所示。
图2.19 常见的4种边缘类型
从这4种边缘类型可以看出一些明显特性,若采用曲线的形式对这些边缘进行拟合,然后对相应的曲线进行一阶导数和二阶导数的计算,则可以计算出不同类型的边缘应该满足的条件。在实际的边缘检测中,可以采用这些条件来进行边缘的判定。在实际处理数字图像时,采用差分近似导数进行处理。
2.4.2 边缘检测的三个阶段
边缘检测可以分为滤波处理、差分处理和最终检测三个阶段。
在滤波处理阶段,实际上是对图像进行滤波处理。例如,通过滤波来去除图像的噪声,这些噪声可能来自成像过程中的采样、量化、模糊、散焦及对象表面结构的不规则性。通过逆滤波来去除各种模糊,通过双边滤波来强化边缘等。
在差分处理阶段,通过利用边缘的特性采取一阶和二阶差分来获取边缘区域。一般情况下,通过计算图像像素值的梯度峰值来确定边缘。显然,梯度值有大有小,必须对其进行阈值处理,从而抑制伪边缘,而阈值的选取又依赖于图像的特性。如何有效地获取恰当的阈值,仍然是一个非常困难的问题。目前常见的阈值方法包括根据前景背景分布来获取最优阈值法,OTSU大津阈值法,以及基于熵的阈值法等。
Haralick边缘检测方法的差分处理看起来与此稍有不同,该方法首先通过分段连续多项式来拟合邻域像素值,然后通过对多项式求偏导数来进行差分处理。但实际上,该方法也包括滤波处理,只不过这里的滤波通过拟合方式来达到。
而常见的Canny边缘检测在滤波阶段使用高斯滤波器,在差分阶段沿梯度方向计算其一阶方向导数,然后在检测阶段,通过检测上一步导数输出的峰值来定位边缘点。
Marr和Hildreth也使用高斯滤波器,通过Laplacian来进行差分阶段的操作。实际上,Laplacian操作相当于沿横轴和纵轴方向二阶偏导数的和。
下面我们分别介绍滤波和相关的差分处理,以及几种常见的边缘检测算法。
2.4.3 滤波操作及双边滤波器
在信号处理中,滤波与预测是紧密相关的概念,实质上就是对输入的数据进行处理,然后产生输出。若输入的数据为以前的数据,而产生的数据为未来的数据,则称为预测;若产生的数据只是对当前数据的校正,则称为滤波。在数字图像处理中,假设原始空域图像用f(x,y)表示,其频域(傅里叶变换)用F(u,v)表示,滤波器用h(x,y)表示,滤波器的频域用H(u,v)表示,图像f(x,y)的滤波输出用g(x,y)表示,频域用G(u,v)表示,则滤波过程可以写为
两边同时进行傅里叶变换
显然,滤波H(u,v)的性质决定了滤波的结果。而若要确定滤波器的效果,则只需要对滤波器h(x,y)进行频谱分析即可。通常所说的低通、高通和带通滤波器,实际上就是通过对H(u,v)进行分析得出明确的结论。
注意,这里一般滤波器中系数的值只与处理像素的距离有关。后面还会谈到,滤波器中的系数不仅与当前处理像素的位置有关,而且与处理位置的像素值大小有关,双边滤波器在滤波的同时还能保持边缘不被平滑。
这里以一维滤波器为例,多维滤波器分析与此类似。假设h(n)={1,1,1},即只有3个值为1,其余值均为0。输入信号和输出信号的关系为g(x)=1/3[f(x−1)+f(x)+f(x+1)]。根据傅里叶变换的性质则有
G(u)=1/3[e-j2πu/N+1+ej2πu/N]F(u)=H(u)F(u)
其频率响应函数为
当u从0到N-1变化时,对应从低频到高频变化,其频谱图如图2.20所示,从中可以看出,低频部分得以保留,而高频部分被削弱,因此该滤波器可以认为是低通滤波器,实际上就是将连续3个输入值的平均值作为输出。
图2.20 滤波器h(n)={1,1,1}对应的频谱图
从图2.20中可以明显看出,低频时其幅度大,高频时其幅度小,产生低通滤波的效果。
最简单的滤波器为均值滤波器,此时滤波器系数均相同,且系数为。此外,另一种常用的低通滤波器为高斯滤波器,其响应函数为。若设定σ=2,则滤波器大小为13×13,并将滤波器系数进行调整,使得最大值为255,则可以得到如图2.21所示的滤波器。
图2.21 一个13×13的高斯低通滤波器
其频谱响应可用Matlab中的freqz2函数画出如图2.22所示的频谱响应。注意,一般,生成的高斯滤波器会对其滤波器系数进行归一化,即保证所有滤波器系数之和为1。思考:如何设计图2.21中的滤波器?
图2.22 图2.21中的高斯低通滤波器的频谱响应
在实际处理中,若形成滤波器模板,则只需要对原始图像用滤波器模板做滑动卷积处理,即可生成滤波后的图像。
从上述的介绍中可知,在滤波器形成滤波器模板后,直接与原始图像做卷积即可达到对图像进行滤波的目的。并且可以知道,滤波器具有局部响应特性,因此通过采用滤波器来获取图像的局部特征,如Gabor特征。此外,基于滤波器的卷积操作在CNN(卷积神经网络,Convolutional Neural Networks)中有着广泛的应用。其中卷积的操作和形式也存在许多变化,包括变形的卷积横板、1×1的卷积横板等。
在图像处理中,假设输入图像用I(i,j)(0≤i<width,0≤j<height)表示,输出图像用表示,卷积核用表示,一般,核的大小为方阵,此时width=height。此时卷积操作可以写为
注意,卷积操作具有可交换性,该式可等价写为
显然,随着(i,j)的移动,可以生成输入图像局部过滤的结果图像。随着卷积核的不同,对应不同的图像。很多常见的图像处理操作都可以表示为这种卷积核的操作。
而常见的用于边缘检测的模板算子包括以下4种。
(1)Prewitt算子:水平算子为,竖直算子为;
(2)Sobel算子:水平算子为,竖直算子为;
(3)Roberts算子:分为45°和135°两种方向,为2×2模板,45°为,135°为;
(4)Laplacian算子:。
前面提到双边滤波器,这里简单进行介绍。一般滤波器只考虑当前处理像素的位置与邻域像素位置距离之间的关系,一般考虑图像的局部平滑特性,距离近的像素所贡献的权值大,也就是滤波器系数值大,反之则小。但前述的滤波操作有一个缺陷,即在去除噪声的同时也会弱化边界。这对边界检测相当不利,因此后来有人提出了双边滤波器。分别综合利用几何空间距离和像素差值来共同决定滤波器系数。
与前述滤波类似,双边滤波器中输出像素与输入像素关系为
可以明显看出,式(2.10)与前述公式的差别就在于滤波器系数的确定形式上。一般,权值h(i,j,x,y)=d(i,j,x,y)·r(i,j,x,y)。其中,称为定义域核,而称为值域核。注意,一般要求h(i,j,x,y)的所有滤波器系数之和为1。可以看出h(i,j,x,y)中的d考虑了像素空间距离的差别,而r则考虑了像素邻域内像素值之间的差别。注意,双边滤波器在滤波的同时考虑了边缘的保持能力。
2.4.4 差分操作
众所周知,连续函数的导数在离散情况下采用差分来近似。这从连续函数的导数定义:可知,若令Δx=1,则成为差分形式。实际上在图像处理中,该形式是相邻像素值的差值,近似为一阶导数。而根据前面介绍的边缘类型可以知道,利用这种相邻像素值的差可以检测边缘。
将其推广到二维的情况,这时图像f(x,y)在点(x,y)处的梯度定义为向量(fx,fy)。其梯度幅度M和方向θ分别定义为
而在方向θ上的方向导数定义为
为了完成对差分的计算,有多种方法来进行运算,如2.4.3节谈及的4种类型的边缘检测算法,实际上就是利用差分进行计算的。
此外,上述梯度幅度值的计算在实际处理过程中也有很多的简化算法。显然,上述求水平、竖直方向梯度的平方和再开方是最准确的,但其涉及平方和开方的操作,因此计算复杂度较高。在实际计算中,若对其要求精度不高,则水平和竖直方向的梯度可以直接采用|Δx|和|Δy|进行简化,进一步可将M简化为M=max{fx,fy}或者M=|Δx|+|Δy|也可达到类似效果。
2.4.5 边缘检测操作
令为图像f(x,y)在像素点(x,y)处的梯度幅度,则M称为图像f的梯度图像。若对梯度图像M进行局部峰值阈值化操作,则可以确定边缘像素的位置为
其中,E(x,y)为边缘图像,Th为阈值,该阈值既可以全局决定,又可以根据局部特性来确定,在应用时可以自适应选择。
2.4.6 非极大值抑制操作
2.4.5节直接利用梯度幅度进行阈值化操作来检测边缘,但并未用到梯度的方向信息。梯度方向表示函数值增加的方向,因此若函数值在某个方向上没有任何变化,则其梯度值为0,如图2.23所示。图中若存在一条竖直边缘线,则边缘线上的任意一点的梯度值都只有水平方向的分量,且方向水平向右,正好与垂直的边缘线垂直。因此,若图像中的点的梯度值变化不显著,则该点很可能不是边缘点。在进行边缘检测时,我们需要抑制这样的非极大值点。该过程可表示为
注意,局部邻域可以有很多的选取方法,如常见的四邻域、八邻域等。选取方法可以根据实际情况进行选择。
图2.23 均匀灰度图像的梯度值
图2.23的每列像素值都一样,每行像素值分别从0~255均匀变化。显然,图像中每个像素点都可以根据前面的梯度公式计算梯度向量,但其竖直方向上像素值没有变化,因此其梯度的竖直方向分量为0,只剩下水平梯度,其方向为水平向右,指向灰度值增加的方向。显然,梯度的方向与边缘方向垂直,而梯度方向表示灰度值变化最剧烈的方向。
若当前点(x,y)处的梯度方向两边的点的梯度幅度表示为M′(x,y)和M′(x,y),则可以只对这两个点进行比较来实现非极大值抑制,可以用图 2.24 表示,只需要比较当前点与边缘两点的幅度大小就可以确定当前点是否是候选边缘点。当然,我们也可以在2.4.5节之前做这个操作,从而获取更好的边缘检测结果。
图2.24 局部非极大值抑制
2.4.7 几种典型的边缘检测算子
在真实边缘检测算子中,为了避免噪声的干扰,一般先对图像进行平滑等预处理,然后再采用上述的检测过程进行检测。应用最为广泛的边缘检测算子之一就是John Canny在1986年提出的Canny算子,它与Marr(LoG)边缘检测方法类似,也属于先平滑后求导数的方法。其次采用连续的函数来逼近图像的局部区域,然后利用连续函数的偏导数来获取其不连续点(即边缘点),典型的方式为Haralick算子。第三类就是LoG算子,实际上相当于对图像进行不同尺度下的平滑操作,然后在对其求差,从而凸显边缘操作。
1.Canny边缘检测算子
在Canny边缘检测中,最显著的特征在于其采用两个阈值来进行判断,高阈值可以按照2.4.5节中的方式来检测那些更具有可能性的真实边缘点,然后根据低阈值来检测高阈值后剩余的点的梯度值,判断其是否可能将高阈值点连接成闭合边缘。这样可以减少噪声对边缘检测的影响。
2.Haralick边缘检测算子
另外一类常见的边缘检测算子是基于图像曲面模型,如Haralick用双三次的多项式来拟合像素的邻域关系。在图像的梯度方向上,根据多项式的系数来计算二阶和三阶导数。在给定像素位置后,可以通过最小二乘法来拟合当前像素邻域内的像素值。当该点处的二阶导数为0且三阶导数为负时,判定该点为边缘点。
例如,若令x,y表示图像中的行列坐标,其像素值f(x,y)采用坐标的双三次多项式来进行表示,则有
此时在点(x,y)处的梯度fx与fy可以求出,若将点(0,0)代入,则可得在点(0,0)处的梯度方向θ的正弦与余弦值分别为
在给定方向向量(sinθ,cosθ)后,可求出其一阶和二阶方向导数,分别为
若令x=ρsinθ,y=ρcosθ代入图像的双三次插值公式,则原公式变换为
可知C0=k00,C1=k10sinθ+k01cosθ,C2=k20sin2θ+k11sinθcosθ+k02cos2θ,C3=k30sin3θ+k21sin2θcosθ+k12sinθcos2θ+k03cos3θ。此时若对幅度ρ求导数,则有
由其判断边缘点的两个条件:和,可得C3<0,2C2+6C3ρ=0,从而有,注意这里ρ0为阈值。
总结Haralick边缘检测算子如下:
(1)根据局部区域的像素值,采用最小二乘法来拟合双三次多项式的系数kij,0≤i+j≤3;
(2)计算θ、sinθ、cosθ;
(3)计算C2、C3;
(4)判断是否为边缘点:若C3<0,则是边缘点,否则不是。
其中第一步拟合多项式的系数中以当前待处理点为点(0,0),周围像素点坐标为
若采用一次多项式来拟合像素位置与像素值的关系,则有
这里只需要确定三个系数,将当前像素邻域坐标代入可得一个超定的线性方程组。若将上述9个像素点的坐标分别用(x1,y1),(x2,y2),…,(x9,y9)表示,其像素点的值分别用(f1,f2,…,f9)表示,则该方程组可表示为
将其表示为f=Ak,采用最小二乘法,可求得k
令B=(ATA)−1AT=(bij),则其为3×9的矩阵,k的每个分量都可以通过卷积的形式来进行计算,如图2.25所示,其中k00、k10、k01正好表示为卷积核与像素值的卷积操作。
图2.25 Haralick边缘检测算子中通过卷积操作来求拟合系数
应用这样的方式,对5×5的局部邻域进行类似的操作,可以求得双三次拟合中的所有系数k(0ij≤i+j≤3),这里实际上可以形成10个卷积核,分别对应这10个系数。实际上,使用这种检测边缘的方式在于其采用连续的曲面模型对图像的局部区域进行建模,从而形成连续的逼近,然后再采用连续函数的偏导数来检测灰度级别的不连续性,从而获取图像的边缘。
3.LoG边缘检测算子
LoG(Laplacian of Gaussian)算子也称Marr算子,在计算机视觉相关领域有着广泛的应用。拉普拉斯(Laplacian)算子是最简单的各向同性微分算子,具有一定的旋转不变性。在一维情况下,类似于二阶导数操作,在二维时其定义为
可以通过前面边缘检测的基本原理来理解,当出现边缘时,灰度值会发生突变,因此对其一阶导数出现较大幅值的概率较大,但需要阈值来进行取舍。而再对一阶梯度再次求导数,则这时往往会出现0交叉的情况,然后就可以用0交叉来判断是否为边缘点,但该过程容易受噪声的影响。为了保证检测免受噪声影响,首先进行高斯低通滤波去除噪声。在不同的尺度上获得不同的表示,然后再应用Laplacian算子进行处理,这就是LoG算子。
假设输入图像为f(x,y),尺度为σ的高斯核定义为
LoG算子可表示为
注意:对卷积的导数等于对其中之一的导数求导后再做卷积,并且高斯函数的导数仍然是高斯函数,因此可以定义LoG算子的核函数为
且有
由此可知,两个不同尺度的高斯函数的差值可以用来近似LoG算子,进而用来检测边缘。
4.Difference of Gaussian算子
从2.4.7节介绍可知,采用两个相邻尺度的高斯低通滤波对图像做处理后,然后做差,可以用来近似LoG算子。可以形式化表示为
其中,DoGk=(Gσ1−Gσ2)称为DoG算子的核函数,然后通过前述的卷积操作可以获得图像的边缘图。这相当于在需要进行连续尺度操作时极大地简化了计算过程,该算子在图像处理中有着广泛的应用,如角点检测。