1 引言

图像的边缘是图像的基本特征之一,清晰准确的边缘提取结果有利于提高特征提取和模式识别等的准确性^[1-3]。但在实际情况中,噪声会影响图像边缘检测结果,使得边界不明显,所以国内外很多学者都对这个难题进行了研究^[4]。

目前,包含Canny算法在内的不少经典图像边缘检测方法都是采用整数阶的微分算法^[5],它们对弱边缘与噪声的分辨能力较差,在较大噪声的图像中容易检测出伪边缘^[6]。为了解决这个问题,Li等^[7]提出改进模板的分数阶微分来实现对医学图像边缘的检测,刁智华等^[8]提出了一种基于改进模糊边缘的边缘提取算法,用于小麦病斑图像的边缘提取。这些改进方法都是基于单一尺度下图像的空间域,较难兼顾图像的边缘、纹理和平滑区域。近些年,随着以小波变换^[9]为代表的多尺度几何理论的快速发展,contourlet变换^[10]、非下采样contourlet变换(NSCT)^[11-12]等新型多尺度算法被提出,它们为图像边缘检测提供新的思路。其中contourlet 变换和NSCT是小波变换的改进,它们的方向分解更充分、选择更灵活,二者弥补了小波变换方向分解不够的缺陷,在边缘检测过程中可以提取更多的方向信息和细节,从而弥补了小波变换对某些方向边缘的检测能力差的缺陷。例如,肖易寒等^[11]提出在NSCT域使用双阈值的模极大值方法获取图像的高频分量,有效地减少了图像的伪边缘,这说明图像高频分量的提取是边缘准确检测的关键之一;苏金凤等^[13]结合了鲁棒主成分和NSCT,实现了多波段图像的压缩融合,该方法较传统算法可以获得更多的清晰度。同时廖一鹏等^[14]通过结合分数阶微分的方法,可以有效地兼顾图像中高频段纹理信息的增强与低频边缘信息的保留,而多尺度图像的融合也进一步优化了边缘信息的保留与提取能力^[15]。付华栋等^[16]将NSCT算法与改进的Canny算法相结合,能够较好地发挥NSCT算法的优势,从而能够识别并提取较为精确的图像边缘,但该方法没有很好地弥补Canny算法的缺陷,导致该方法处理噪声图像的效果一般。

综上分析,可以将分数阶微分运用于NSCT域高频子带的边缘检测中,结合高频子带的相关性,通过识别更多的边缘信息,可获得更精确的边缘提取结果。为此,本文将轮廓提取算法、NSCT和分数阶微分结合起来,提出了一种在NSCT域结合分数阶微分的噪声图像边缘检测方法。在使用NSCT对噪声图像进行分解后,先在低频图像中使用改进的轮廓提取算法,该算法引入动态梯度自适应分数阶微分并基于自适应高、低阈值连接边缘,从而完成低频子带的边缘轮廓提取;其次,通过利用NSCT域高频子带具有方向性和尺度间具有相关性的特点,结合具有方向的分数阶微分模板,可实现强、弱边缘信息的提取与噪声的抑制,从而获得高频边缘提取结果;最后,将NSCT域尺度内和尺度间的检测结果进行融合,从而得到最终的边缘图像检测结果。通过多张图像的边缘提取结果验证本文算法的有效性,并将其与其他几种边缘检测方法进行比较,以验证本文方法的有效性。

2 相关理论

2.1 分数阶微分及掩模算子

Grümwald-Letnikov(G-L)的分数阶微分定义是从连续函数整数阶导数的经典定义出发, 将微积分的阶数由整数向非整数推衍而来的^[17],即

aGDtv=limh→01hv∑m=0t-vh(-1)mΓ(v+1)m!Γ(v-m+1)×f(t-mh),(1)

式中: aGDtv为分数阶微分算子,其中v为分数阶微分的阶数,t为自变量,h为积分步长;f(·)为源函数;Γ(n)为Gamma函数, Γ(n)=∫0∞e-ttn-1dt=(n-1)!。

假设现在一个信号中t的范围为:t∈[a,b],在t的范围里可以采用间隔进行划分,于是可以得出:n= b-ahh=1=[b-a],由此可以推导出该信号分数阶微分的差分表达式为

dvf(t)dtv≈f(t)+(-v)f(t-1)+(-v)(-v+1)2f(t-2)+…+Γ(-v+1)n!Γ(-v+n+1)f(t-n)。(2)

由(2)式可以直接获得差分方程的第n项的系数为

an=Γ(-v+1)n!Γ(-v+n+1)。(3)

对二维数字图像进行分数阶微分增强就是利用分数阶微分的掩模算子和图像像素点进行卷积处理,本文方法使用了前3项系数,采用5×5的分数阶掩模算子,选择8个方向对图像进行计算,得到的Tiansi算子如图1所示,它可以兼顾结果的误差与方向性。

图 1. Tiansi算子

Fig. 1. Tiansi operator

下载图片 查看所有图片

2.2 NSCT

原理上,具备了平移不变性的多尺度变换可以对图像中的噪声进行更好的抑制^[18],但contourlet变换过程中需要进行下采样,导致了不同尺度上的分解结果不同,使coutourlet变换在一定程度上失去了平移不变性的优势,因此da Cunha等^[19]提出了NSCT,该变换保留了contourlet变换具有的多尺度、多方向的频域拆分能力等优势的同时,强化了处理过程中的平移不变性,因此该变换能够提高边缘信息的分离与识别能力。

NSCT流程如图2所示,首先通过非下采样金字塔滤波器(NSP)对源图像进行多尺度分解,将图像的低频和高频信息分解到多张图像中;再通过非下采样方向滤波器组(NSDFB)将高频图像进行多方向分解,得到多张具有方向特征的方向子带图像,从而获得了高频图像的方向信息^[19]。

图 2. NSCT分解图

Fig. 2. Decomposing process of NSCT

下载图片 查看所有图片

3 NSCT域多尺度边缘检测

NSCT域多尺度边缘检测算法就是在不同尺度的图像上,根据其特点,选用不同的针对性算法提取边缘信息,最后进行融合,得到最终的边缘图像。

3.1 低频图像轮廓提取

为了更好地提取低频子带中的边缘轮廓信息,使用适合低频子带边缘检测的改进轮廓提取算法进行提取。处理流程为:

1) 对图像进行NSCT后,噪声大多会集中在高频子带,低频图像的噪声较小,为了避免高斯平滑导致的边缘模糊问题,直接使用动态梯度的自适应分数阶微分对低频图像进行增强处理,对于一幅m×n大小的低频图像f(i,j),计算像素点(i,j)的8邻域梯度的平均值M(i,j):

M(i,j)=18∑m=-11 ∑n=-11f(i,j)-f(i+m,j+n)。(4)

对于每个像素点(i,j),根据其梯度平均值M(i,j)使用不同阶数的分数阶微分模板进行增强处理:

v(i,j)=0.6,20<M(i,j)0.4,2<M(i,j)≤200.3,M(i,j)≤2,(5)

式中:v(i,j)是像素对应的分数阶微分的阶数。

2) 记增强后的低频图像为S(i,j),通过水平与垂直方向上的一阶有限差分近似式来计算两个方向的偏导数P_x(i,j)与 P_y(i,j),并进一步通过计算获得梯度幅值I(i,j):

Px(i,j)≈∑k=01S(i+1,j+k)-S(i,j+k),(6)Py(i,j)≈∑k=01S(i+k,j+1)-S(i+1,j),(7)I(i,j)=Px(i,j)2+Py(i,j)2。(8)

3) 对梯度幅值结果进行非极大值抑制,以细化梯度幅值屋脊。

4)为了提高阈值处理的鲁棒性,使用自适应高低阈值进行边缘提取。将所有的像素根据梯度幅值的大小进行升序排列后,得到具有相同的梯度幅值的像素的数量为r_i(i=0,1…),将像素点根据梯度幅值分为非边缘点、疑似边缘点和边缘点三类,分别记为A₁、A₂和A₃,其中A₁和A₂之间的梯度幅值阈值为T_L,而A₂和A₃之间的梯度幅值阈值为T_H。

图像各梯度幅值的分布概率p_i为

pi=rim×n,i=(0,1,2…)。(9)

A₁、A₂、A₃的平均相对熵ε₁、ε₂和ε₃为

εk=-∑ri∈AkpiPk×lnpiPk(k=1,2,3),(10)

其中

Pk=∑ri∈Akpi(k=1,2,3)。(11)

根据熵值最大的原理确定最佳的阈值T_L和T_H,确定T_L和T_H后,如果满足条件:

ε1+ε2+ε3TL;TH=maxε1+ε2+ε3,(12)

此时,所取的T_L和T_H即为所求的最佳高、低阈值,从而可获得三类阈值点。之后通过连通性判断,确定A₂类像素中哪些属于边缘点,就可以获得轮廓提取结果。

3.2 多尺度高频子带边缘检测

高频方向子带均包含具有方向特征的边缘信息和不具有方向特征的噪声,NSCT对噪声有较大的抑制,因此可以通过阈值处理方式,初步区分有用信息与噪声。假设第n个尺度上第m个子带存在一个阈值 tnm^[14]:

tnm=121m×n∑i=1m∑j=1nWnm(i,j)-W̅nm2,(13)

式中: Wnm(i,j)表示第n个尺度上第m个子带中(i,j)位置的系数; W̅nm为子带系数的均值。考虑到同一尺度的不同方向子带中,噪声的强度大致相同,而有效信息的强度会随着图像的方向性特征的不同而发生改变,如果某个方向子带的能量较低,则说明该子带上噪声较多,这时应当提高阈值、提高抗噪能力,因此本文根据能量特征改进阈值:

Tnm=∑i=1LEniL×Enmtnm,(14)

式中: Enm为第n个尺度上第m个子带的能量;L为第n个尺度的方向总数。

在同一子带中,不同区域的特征也有所不同,因此单一阈值的去噪方法会因为某些弱边缘强度不足而将其错误去除,从而导致间断点产生。经过实验观察发现,弱边缘大部分都集中在[0.8 Tnm,1.2 Tnm],因此本文方法结合子带具有方向相关性的特点,对各个子带进行自适应增强,以提升弱边缘和噪声的区分度。

定义第n个尺度上第m个子带上的尺度相关系数:

rnm(i,j)=∏k=1KWkm(i,j)EkmE∏k=1KWkm(i,j),(15)

式中: ∏k=1KWkmi,j为同一位置不同尺度的系数的乘积,而 Ekm/E∏k=1KWkm(i,j)是归一化系数,该系数越大,说明该位置各尺度间相关性越好。根据以上结果,计算各个像素的增强系数 vnm(i,j),该系数为分数阶微分的阶次,可表示为

vnm(i,j)=0,Mnm(i,j)<0.8Tnm0,0.8Tnm≤Mnm(i,j)<1.2Tnmrnm(i,j)<Mnm(i,j)0.7,0.8Tnm≤Mnm(i,j)<1.2Tnmrnm(i,j)≥Mnm(i,j)0.7,Mnm(i,j)≥1.2Tnm。(16)

根据NSCT的时频特性,高频信息具有较强的方向性,例如分解方向个数为4的分解结果中包含了水平、竖直、左对角、右对角4个方向的尺度分解图,该方向性对特定的子带是确定的,因此通过对Tiansi算子进行改进,可使其能够根据方向性对图像进行选择性增强,即

wh=v2-v120v2-v120v2-v120-v6-v6-v60006000-v6-v6-v60v2-v120v2-v120v2-v12,(17)

式中:v为增强系数 vnm(i,j),这个矩阵去除了水平方向的系数,使算子可以获取更多的垂直、左对角、右对角方向的信息,从而可增强水平方向的边缘信息强度,抑制其他方向的噪声。对于其他方向的尺度图像,先进行一次旋转操作,使其中的方向信息变为水平方向的信息,之后通过该矩阵进行增强处理,最后执行一次反向旋转操作进行还原。

3.3 本文方法流程

本文提出的结合分数阶微分的噪声图像NSCT域边缘检测方法的具体步骤如图3所示,具体步骤如下:

1) 对输入图像执行NSCT,得到低频图像和高频各子带分量信息;

2) 使用改进轮廓提取算法对低频图像进行处理,通过分数阶微分算子对图像进行增强操作后,对图像的梯度幅值进行非极大值抑制,之后基于熵值最大的原理,从中提取出边缘点与疑似边缘点,最后结合连通性得到低频轮廓提取结果;

3) 对于高频子带,先根据子带系数与同位置各尺度间的相关性,利用不同阶次的分数阶微分模板对系数进行增强,并通过自适应阈值处理来进行噪声抑制,再对各个尺度的高频信息进行融合,最后进行二值化处理,得到高频边缘提取结果;

4) 对高频和低频结果进行相加融合,生成最终的边缘图像,并输出结果。

图 3. 方法流程框图

Fig. 3. Flow chart of the proposed method

下载图片 查看所有图片

4 实验结果与分析

所有实验都是采用Windows 10系统的MATLAB R2018a环境。本文方法中NSCT采用“maxflat”塔型分解和“dmaxflat7”方向滤波器进行3个尺度的分解。其中第1个尺度为低频子带,对其不进行方向分解;第2、3个尺度则为高频子带,将其对应的方向数目分别设定为4、8。因篇幅限制,文中只展示出几幅纹理边缘信息较丰富的图像处理的效果。

4.1 原图像边缘检测效果

理论上来说,分数阶微分的引入可以帮助提取更多高频纹理边缘信息,故本文选用纹理边缘信息较丰富的256 pixel×256 pixel的动漫人物图像和标准测试灰度图像Lena进行边缘检测实验,原图如图4所示。

图 4. 实验图像原图。 (a) Lena; (b)动漫卡通图像

Fig. 4. Original images for experiment. (a) Lena; (b) cartoon image

下载图片 查看所有图片

利用文中提出的方法对两幅图像进行处理,得出低频分量、高频尺度1、高频尺度2、高频分量和最终融合的边缘结果,结果如图5和图6所示。可以发现本文方法所获取的图像的低频轮廓信息比较完整,其中动漫人物图像中的图案轮廓信息明显,基于Lena图像能大体检测出帽子的轮廓。而高频信息包括边缘和纹理信息,很多细节信息都可以体现出来。如动漫人物图像中图案纹理更加清晰完整,Lena图像中头发信息更多并且有层次感。最终把两幅图像的低频轮廓信息和高频边缘纹理信息加以融合,从而可以观察到边缘完整、纹理清晰的动漫人物图像和Lena图像。

图 5. Lena图像的边缘检测结果。(a)低频; (b)高频尺度1; (c)高频尺度2; (d)高频; (e)融合结果

Fig. 5. Edge detection results of Lena image. (a) Low frequency; (b) high frequency scale 1; (c) high frequency scale 2; (d) high frequency; (e) fusion result

下载图片 查看所有图片

图 6. 动漫图像的边缘检测结果。(a)低频; (b)高频尺度1; (c)高频尺度2; (d)高频; (e)融合结果

Fig. 6. Edge detection results of cartoon image. (a) Low frequency; (b) high frequency scale 1; (c) high frequency scale 2; (d) high frequency; (e) fusion result

下载图片 查看所有图片

4.2 噪声图像边缘检测效果

本实验选择纹理边缘信息较丰富的标准测试灰度图像Lena和两幅医学图像作为边缘检测图像进行测试,先利用Matlab中的imnoise函数对图像进行处理,添加不同比例的高斯噪声,再通过本文方法进行轮廓提取,其结果如图7、图8和图9所示,可以发现,在不同的高斯噪声密度下,本文方法都能较好地提取出三幅含噪图像的边缘信息,并且噪声密度的大小对提取结果的影响不大。

图 7. Lena图像的边缘提取结果。 (a)原图; (b)加噪10%的图像; (c)加噪10%的提取结果; (d)加噪50%的图像; (e)加噪50%的提取结果

Fig. 7. Edge extraction results of Lena image. (a) Original image; (b) image with 10% noise; (c) extraction result with 10% noise; (d) image with 50% noise; (e) extraction result with 50% noise

下载图片 查看所有图片

图 8. 医学图像1的边缘提取结果。 (a)原图; (b)加噪10%的图像; (c)加噪10%的提取结果; (d)加噪50%的图像; (e)加噪50%的提取结果

Fig. 8. Edge extraction results of medical image 1. (a) Original image; (b) image with 10% noise; (c) extraction result with 10% noise; (d) image with 50% noise; (e) extraction result with 50% noise

下载图片 查看所有图片

图 9. 医学图像2的边缘提取结果。 (a)原图; (b)加噪10%的图像; (c)加噪10%的提取结果; (d)加噪50%的图像; (e)加噪50%的提取结果

Fig. 9. Edge extraction results of medical image 2. (a) Original image; (b) image with 10% noise; (c) extraction result with 10% noise; (d) image with 50% noise; (e) extraction result with 50% noise

下载图片 查看所有图片

4.3 原始图像多种边缘提取方法比较

针对无噪声图像,采用Canny 算子、文献[ 7]中的改进分数阶微分方法、文献[ 15]中的基于NSST的多尺度边界检测方法、文献[ 20]中的NSCT结合KFCM聚类方法对上述两幅无噪图像进行处理,这些方法的结果与本文方法的检测结果如图10与图11所示。

图 10. 各方法对Lena图像的提取结果。(a) Canny; (b)文献[ 7]中方法; (c)文献[ 15]中方法; (d)文献[ 20]中方法; (e)本文方法

Fig. 10. Extracted results of methods for Lena image. (a) Canny; (b) method in Ref. [7]; (c) method in Ref. [15]; (d) method in Ref. [20]; (e) proposed method

下载图片 查看所有图片

图 11. 各方法对动漫图像的提取结果。(a) Canny; (b)文献[ 7]中方法; (c)文献[ 15]中方法; (d)文献[ 20]中方法; (e)本文方法

Fig. 11. Extracted results of methods for cartoon image. (a) Canny; (b) method in Ref. [7]; (c) method in Ref. [15]; (d) method in Ref. [20]; (e) proposed method

下载图片 查看所有图片

从图10和图11中可以观察到,基于Canny算子的方法得到的结果中边缘信息比较丰富,提取了过多的细节信息,其中包含了一些伪边缘信息,这导致有些边缘轮廓不是很清晰;从文献[ 7]中方法得到的结果中可以看出,较多的边缘信息没有被检测出来,边缘曲线会出现断裂的情况。文献[ 15]中方法容易获得较多的边缘信息,但连贯性不是很好;采用文献[ 20]中方法检测出的图像边缘效果较文献[ 15]中方法得到的图像细节更丰富,但有些边缘地方断裂、轮廓不完整。而本文方法与其他边缘检测方法相比,除了能够获得较好的边缘提取结果外,图10(e)所标区域中轮廓更完整连贯,细节边缘的定位更准确,可以获得更多的边缘轮廓和纹理信息,这些结果证明本文方法对于无噪图像具有更好的边缘提取效果。

为了可以客观地评价各个方法的边缘检测结果,本文选用文献[ 21]中的连续边缘像素比(R)来分析以下几种检测方法的性能,R可定义为

R=α/β,(18)

式中:α表示提取结果中连续边缘像素点的数量;β表示提取结果中所有边缘像素的数量;R可以用来表征边缘提取结果的连续性,其值越大,则连续性越好。

表1和表2中给出了5种方法的α、β和R值,从表中可以发现,各个方法的处理结果中,本文方法的α、β和R值都是最大的,这就表明本文方法对图像的强、弱边缘信息的识别性最好,同时所得结果具有较好的边缘连续性,因此本文方法可在一定程度上避免边缘断裂的问题。

4.4 噪声图像多种边缘提取方法比较

将本文方法与文献[ 7]中的改进分数阶微分方法、文献[ 15]中的基于NSST的多尺度边界检测方法、文献[ 16]中的基于NSCT与改进Canny相结合的边缘检测算法与文献[ 20]中的NSCT结合KFCM聚类方法对含噪图像的处理结果进行了比较。具体地,采用上面5种方法对上述高斯噪声密度为10%的三幅图像进行边缘检测,结果如图12~14所示。

图 12. 各方法对加噪Lena图像的处理效果图。 (a)文献[ 7]中方法; (b)文献[ 15]中方法; (c)文献[ 16]中方法; (d)文献[ 20]中方法; (e)本文方法

Fig. 12. Processing results of different methods for noisy Lena image. (a) Method in Ref. [7]; (b) method in Ref. [15]; (c) method in Ref. [16]; (d) method in Ref. [20]; (e) proposed method

下载图片 查看所有图片

图 13. 各方法对加噪医学图像1的处理效果图。 (a)文献[ 7]中方法; (b)文献[ 15]中方法; (c)文献[ 16]中方法; (d)文献[ 20]中方法; (e)本文方法

Fig. 13. Processing results of different methods for noisy medical image 1. (a) Method in Ref. [7]; (b) method in Ref. [15]; (c) method in Ref. [16]; (d) method in Ref. [20]; (e) proposed method

下载图片 查看所有图片

图 14. 各方法对加噪医学图像2的处理效果图。(a)文献[ 7]中方法; (b)文献[ 15]中方法; (c)文献[ 16]中方法; (d)文献[ 20]中方法; (e)本文方法

Fig. 14. Processing results of different methods for noisy medical image 2. (a) Method in Ref. [7]; (b) method in Ref. [15]; (c) method in Ref. [16]; (d) method in Ref. [20]; (e) proposed method

下载图片 查看所有图片

从图12~14中可知,文献[ 7]中方法虽然有较好的抗噪能力,但检测到的边缘信息比较少,边缘断裂情况明显。文献[ 15]中方法有较好的去噪能力,其连续性较文献[ 7]中的方法好,但仍然无法较好地检测出图像的边缘信息,且依然存在断裂的情况。文献[ 16]中方法虽然能够较好地提取有效的边缘信息,但检测出来的伪边缘也较多,这说明该方法容易受到噪声影响。文献[ 20]中方法相对于文献[ 15]中方法来说,检测到的图像边缘信息更清楚,但是还是存在一些噪声信息。本文方法较前面几种方法获取的边缘信息更丰富、更清晰,同时边缘更为连续,边缘位置更为准确,去噪效果更好,因此与其他方法相比,本文方法具有很大的优势。

为了客观地评价针对含噪图像的实验中得到的边缘检测结果,先使用各个方法对具有不同强度的噪声进行处理,并使用R值来分析,以对各个方法进行比较,具体结果如图15所示。

图 15. 不同噪声下各方法的结果对比。(a) Lena图像;(b)医学图像1;(c)医学图像2

Fig. 15. Comparison of results of different methods for different noise. (a) Lena image; (b) medical image 1; (c) medical image 2

下载图片 查看所有图片

图15中给出了三幅图像不同高斯噪声密度下5种方法的R值变化图,从图中可以看到:本文方法的R值都是最大的,这证明了本文方法获取到的边缘连续性最好。对于Lena图像、医学图像1、医学图像2,从无噪到加入5%噪声,本文方法的R值分别会下降5.3%、6.0%、4.4%,文献[ 7]中方法的R值分别会下降24.6%、30.3%、29.1%,文献[ 15]中方法的R值分别会下降10.7%、10.1%、7.4%,文献[ 16]中方法的R值分别下降10.1%、9.7%、7.6%,文献[ 20]中方法的R值分别下降6.8%、9.3%、5.7%。本文方法相对其他4种方法的R值的下降幅度最小,而且随着噪声密度的不断加大,本文方法检测的边缘R值下降趋势相对缓慢,并且R值一直是最大的。这说明了本文方法相对其他几种方法来说抗噪性能会更好一些。

5 结论

提出了一种结合分数阶微分的噪声图像NSCT域边缘检测方法。首先对原始图像进行NSCT分解,得到低频与高频分量,利用改进轮廓提取算法提取低频的边缘轮廓,然后利用NSCT域的多尺度积和边缘点的相关性对高频系数进行自适应阈值去噪,使用不同阶次的方向分数阶微分模板对不同类型的系数进行自适应增强,并将其融合成高频边缘图像,再通过阈值处理去除噪声并进行二值化,最后将处理完成的高、低频图像进行相加融合,以获得完整的边缘图像。对不同类型的图像进行实验,将本文方法与其他现有方法进行对比,结果表明:采用本文方法对图像进行边缘检测得到的无噪图像的平均连续边缘像素比为0.931,含10%高斯噪声图像的平均连续边缘像素比为0.861。本文方法检测得到的边缘位置更准确,图像轮廓信息更完整,既能针对无噪灰度和彩色图像、也可以针对含噪的各类图像进行检测,且随着噪声水平的上升,本文方法保持较高的平均连续边缘像素比,抗噪性能好,边缘完整、连续。