2. 南京大学 地理与海洋科学学院 地理信息科学系，江苏 南京 210093
2. Department of Geographical Information Science,School of Geographic and Oceanographic Sciences, Nanjing University, Nanjing 210093, China
1 引 言
边缘检测是图像处理中重要的基础研究之一，它对于图像识别和理解具有重要意义。早期算法主要是利用边缘检测算子与图像的卷积运算来完成边缘检测工作，如Robert[1]、Prewitt[2]、Sobel[3]、Canny算子[4]和LoG算子[5]。此外学者们还提出基于小波变换[6, 7]、相位一致模型[8, 9, 10]、信息熵[11]、统计方法[12]和神经网络[13]进行边缘检测。它广泛应用各个领域，包括文字识别[14] 、人脸识别[15]、虹膜识别[16]、医学图像处理[17, 18]等。遥感图像处理领域中，地物边缘特征对于遥感图像信息提取至关重要。文献[19]利用小波变换实现SAR图像边缘检测；文献[20]利用Canny方法实现遥感图像的海岸线自动提取；文献[21]利用相位一致模型研究高分辨率遥感图像的多尺度边缘检测。
梯度、偏微分等算子和小波变换均是通过空域卷积运算进行边缘检测，卷积模板越大，则计算量也越大。虽然相位一致模型可以通过空域和频域两种方式完成，但都涉及很大计算量。这对于大数据量的图像尤其是遥感图像而言，边缘检测将耗费大量时间。同时传统方法是从水平和竖直方向分别进行运算，并通过平方和最终完成边缘检测，因此其他方向的边缘特征未能得到很好的响应。
本文从频域的角度出发，提出了基于二维Log Butterworth滤波器的全方向边缘检测的频域实现方法。该方法具有全方向边缘检测特点，同时可以更快地实现边缘检测。
2 改进Log Butterworth滤波器 2.1 改进一维Log Butterworth滤波器Butterworth滤波器最早由Butterworth在1930年提出，它具有通带内最大平坦的振幅特性，且在正频率内的振幅随频率升高而单调下降，常用于低通滤波，定义为
式中，Ωp为通带上限截止频率。通常，选择通带允许的最大衰减为3 dB，此时ε=1。
图 1给出了ε=1时，阶数N取不同值时的Butterworth滤波器形状。可见该滤波器的特点是N取值越大，通带内的频率响应曲线越平坦，则滤波器的衰减速度也越大。
|
| 图 1 N取不同值时的Butterworth滤波器 Fig. 1 The Butterworth low pass filter with different N |
Butterworth带通滤波器定义为
式中，Ω2为中心频率；Ω3为通带上限截止频率。
文献[22]从生物学的角度研究证明，人类识别纹理等图像特征的过程应涉及非线性的识别机制。因此，有学者将非线性的坐标系统引入到滤波器设计中，如Log坐标系。其中成功的改进有文献[23]提出的Log Gabor滤波器。它在Log坐标系下保持高斯形状，而在笛卡儿坐标系保持直流分量为0，解决了Gabor滤波器的带宽不能大于一倍频的局限性。由于自然图像的大部分能量分布在直流分量上，因此保持滤波器的直流分量为0，可以更好地突出纹理、边缘等中高频图像特征。Butterworth带通滤波器转换到Log坐标系后，由于它在高频部分带有一个明显的“尾巴”，所以，滤波结果可以更好地表现纹理和边缘等中高频特征。因此，Log Butterworth滤波器的频域响应为
不同参数的一维Log Butterworth滤波器如图 2所示。
|
| 图 2 不同参数设置的一维Log Butterworth滤波器的形状 Fig. 2 One-dimensional Log Butterworth filter with different parameters |
本文用模拟的一维离散信号(图 3(a))来验证Log Butterworth滤波器的边缘检测效果。信号长度为299；图 3(b)为幅度谱。根据式(2)和式(3)，给出具有相同参数的一维Butterworth带通滤波器(图 3(c))和Log Butterworth滤波器(图 3(d))，参数设置为N=2,Ω2=190,Ω3=210,ε=1。通过傅里叶变换、频域相乘和傅里叶反变换等完成滤波，结果分别如图 3(e)和图 3(f)所示。从最终的滤波结果可以看出：改进后的Log Butterworth滤波器具有更好的边缘检测效果。
|
| 图 3 相同参数下一维带通Butterworth和Log Butterworth滤波器的边缘检测结果 Fig. 3 Edge detection results based on 1-D Butterworth and Log Butterworth filter with the same parameters |
由于一维Log Butterworth滤波器在笛卡儿坐标系上是非线性的，且不以中心频率位置对称，无法使用旋转法来设计二维Log Butterworth滤波器。因此，本文利用极坐标的表达方式来描述二维Log Butterworth滤波器，公式如下
式中，HΩ和Hθ分别为滤波器的径向分量和角向分量；NΩ为滤波器径向分量的阶数；Nθ为滤波器角向分量的阶数。HΩ和Hθ两者分别控制径向和角向的滤波器形状。通过将径向和角向分量的相乘得到二维Log Butterworth滤波器。由于可以通过角向分量控制滤波器角向的分布，因此Log Butterworth滤波器呈扇环形分布。在ε、HΩ、Hθ、NΩ和Nθ确定的前提下，Ω2决定了滤波器的空间分布位置，Ω3和Ω2则共同决定了滤波器带宽。
由于频谱图以直流分量为中心对称，而本文提出的是全方向边缘检测方法，所以角向分量的带宽θ3-θ2设为π，同时角向分量的中心频率位置与提取效果无关。则二维Log Butterworth滤波器的表达式可改为
给定Lx和Ly为图像的长与宽，若Lx≠Ly，即图像行列数不一致，这也导致频谱图行列数的不一致。因此二维Log Butterworth滤波器的中心频率位置与直流分量间距离应随着角度变化而有所不同，从而满足不同行列数图像的边缘检测工作。中心频率应分布于在长短轴之比为图像长宽比的椭圆上，满足以下公式
可知中心频率最大值为
，当中心频率在x与y轴上的位置分别为lx和ly时，满足如下公式
其他方向上的中心频率分布于分别以lx和ly为短轴和长轴的椭圆上，如以下公式所示
由此公式(5)可化为
由图 4可以明显看出滤波器呈椭圆分布。
|
| 图 4 二维Log Butterworth滤波器(300×500)的形状 Fig. 4 Two dimensional Log Butterworth filter (300×500) |
针对不同大小图像，为了方便滤波器参数选取，有必要将中心频率进行归一化。同时，Ω3大于中心频率Ω2，且带宽Ω3-Ω2小于Ω2，因此1Ω3/Ω2＜2。为减少阻带分量造成的双边缘现象影响，Nθ应取较大值。令radius=，α=Ω2/radius，β=Ω3/Ω2，ε=1，NΩ=2，Nθ=18，则有
式中，0＜α＜1; 1＜β＜2。
3 算法实现及参数确定 3.1 算法实现设原图像为f(i,j)，边缘检测可由以下公式实现
边缘响应方向为
式中
式中，H(Ω,θ)y和H(Ω,θ)x为x和y方向上的边缘响应值。具体边缘检测步骤为：
(1) 对图像f(i,j)进行快速傅里叶变换，得频谱F(u,v)。
(2) 将Log Butterworth滤波器与F(u,v)进行相乘运算。
(3) 再对相乘结果进行快速反傅里叶变换，并取绝对值得到边缘检测结果g(i,j)。
(4) 利用非极大值抑制获得单像元边缘检测结果图。
3.2 滤波器参数确定滤波器参数α和β决定了最终边缘检测结果，因此有必要对参数进行优化选取。边缘定位准确性和抗噪声能力是两个主要考虑方面。因此本文利用F-measure[24]和PSNR(峰值信噪比)来描述不同参数下的定位准确性和去噪效果，并最终确定最佳参数范围。F-measure定义如下
式中，S1为标准参考图中的边缘集合；S2为待评价图的边缘集合；给定阈值td，当S1中边缘点p1与S2中所对应边缘点p2间的距离小于td，视参考图与待评价图中的两点匹配；Match(S1，S2)为S1中边缘点可以在S2中找匹配的点的总个数；Match(S2，S1)为S2中边缘点可以在S1中找到匹配点的总个数。本文将k设置为0.5。PSNR定义为
式中，图像灰阶为0~255；M和N为图像的行列数；I0(i,j)为原图边缘检测结果；I(i,j)为加噪后的边缘检测结果。
选取5幅不同类型的图像(图 5)作为试验数据来确定参数α和β的最佳参数范围。图 5(a)为加噪后的合成图像；图 5(b)为鸟的图像；图 5(c)为人类面部图像；图 5(d)为自然风景图像；图 5(e)为高分辨率遥感图像。图 5第2排给出了不同参数下边缘检测的F-measure曲面，F-measure越大定位效果越好。5幅图像检测结果的F-measure最大值如表 1所示。
|
| 图 5 试验数据及不同参数下边缘检测结果的F-measure曲面(第2排)和PSNR曲面(第3排) Fig. 5 Original image and F-measure (second row) and PSNR (third row) estimations of edge detection with different parameters |
| 图 5(a) | 图 5(b) | 图 5(c) | 图 5(d) | 图 5(e) | |
| (α,β) | (0.20,1.70) | (0.25,1.90) | (0.20,1.95) | (0.2,1.75) | (0.25,1.95) |
| F-measure | 0.905 1 | 0.865 3 | 0.705 9 | 0.649 8 | 0.889 8 |
从表 1可以看出，当α取值在[0.20,0.25]内，F-measure均具有较高值，且趋于稳定。同时当最大值时，β均大于1.70。
图 5第3排给出不同参数下边缘检测的PSNR曲面，PSNR越大则噪声影响越小。五幅图像检测结果的PSNR最大值如表 2所示。
| 图 5(a) | 图 5(b) | 图 5(c) | 图 5(d) | 图 5(e) | |
| (α,β) | (0.25,1.65) | (0.25,1.95) | (0.25,1.85) | (0.3,1.95) | (0.25,1.60) |
| PSNR | 26.303 3 | 31.497 0 | 23.610 7 | 25.717 2 | 17.225 9 |
从5个PSNR曲面可以看出，当α取值在[0.25,0.30]，PSNR均具有较高值，且趋于稳定。同时β均大于1.65。
F-measure和PSNR曲面的整体分布相近。F-measure和PSNR的最佳取值范围相似，而且最佳参数范围稳定，同时该取值范围内的F-measure和PSNR取值的差异较小。最终最优参数范围可选定为α为[0.20,0.30]，β为[1.70,0.95]。
4 结果与评价 4.1 耗时分析以下从计算耗时的角度分析本文方法的效率。时间复杂度仅代表运算次数的最高次，在特定的条件下，并不能准确描述算法效率。本文通过乘法和加法次数来分析算法的计算耗时。设图像大小为N×N，正反快速傅里叶变换的时间复杂度为O(N2log2N)，而实数乘法与加法次数分别为3N2log2N－5N2+8和11/2N2log2N－13/6N2+8/3。此外滤波还涉及相乘运算，乘法次数为N2。据式(11)，实数乘法与加法次数分别为6N2log2N－9N2+16和11N2log2N－13/3N2+16/3。
Canny方法(时间复杂度为O(N2))中，首先进行平滑滤波，再从x和y方向分别进行卷积运算，然后做平方和运算，并通过开方运算得到最终检测结果。平滑滤波中，通常是以9×9大小的模板来进行高斯平滑滤波，乘法次数为81N2，加法次数为80N2；卷积运算中，若采用最小的3×3模板，则在两个方向上的乘法次数均为9N2，加乘法次数均为8N2；平方运算中有2N2次乘法运算，N2次加法运算。此时，乘法次数总计为101N2，加法次数为97N2。同理可知模板为5×5时的乘法和加法次数分别为129N2次和133N2次(表 3)。
| Canny算法(3×3) | Canny算法(5×5) | 本文方法 | |
| 时间复杂度 | O(N2) | O(N2) | O(N2log2N) |
| 乘法次数 | 101N2 | 133N2 | 6N2log2N－9N2+16 |
| 加法次数 | 97 N2 | 129 N2 | 11N2log2N－13/3N2+16/3 |
定义RM和RA来描述两算法的乘法与加法次数比率
式中，Num_MULLog-Butterworh和Num_MULCanny分别为本文方法和Canny算法的乘法次数；Num_ADDLog-Butterworh和Num_ADDCanny分别为本文方法和Canny算法的加法次数。
表 4和表 5为卷积窗口为3×3和5×5时的RM3×3、RM5×5、RA3×3和RA5×5。图 6为该4个函数曲线。
| RM3×3 | RM5×5 | |
| 乘法次数比率 |  |  |
| AM3×3 | AM5×5 | |
| 加法次数比率 |  |  |
当Canny算子的模板大小为3×3，N小于330 281时，本文方法的乘法次数小于Canny方法(图 6(a))；当模板大小为5×5，N小于13 316 085时，本文方法的乘法次数小于Canny方法(图 6(b))；当模板大小为3×3，N小于593时，本文方法的加法次数小于Canny方法(图 6(c))；当模板大小为5×5，N小于4455时，本文方法的加法次数小于Canny方法(图 6(d))。由于乘法运算通常是由加法运算来实现，因此，乘法次数对起决定性作用。图像处理领域中，所用图像大小通常小于330 281像素×330 281像素和13 316 085像素×13 316 085像素。这也说明了本文方法能够更快地实现边缘检测。
|
| 图 6 Canny算子的分布曲线 Fig. 6 Diagram of Canny |
以1000像素×1000像素和3000像素×3000像素大小图像为试验数据，采用Matlab软件提供的快速傅里叶变换函数和Canny函数，来对比本文方法和Canny方法的运行效率。其中计算机处理器(CPU)为Inter(R)Core(TM)i5 3.20 GHz，内存(RAM)为8 GB，边缘检测耗时如表 6所示。可以看出本文方法效率速度明显优于传统的基于偏微分算子的Canny方法。
| 图像大小/像素 | Canny方法/s | 本文方法/s |
| 1000×1000 | 0.314 168 | 0.196 368 |
| 3000×3000 | 2.749 971 | 1.736 792 |
试验数据为BSDS图像。图 7给出了本文方法、Canny、Sobel、相位一致模型和数学形态学等方法的边缘检测结果。本文方法的滤波器参数为α=0.25，β=1.9，θ=π。Canny算子参数σ=1，卷积窗口为9×9；Sobel算子卷积窗口为3×3；相位一致模型中方向参数为6，尺度为4，σ=0.55；数学形态学方法中结构元素为3×3的方形窗口。
|
| 图 7 不同方法的边缘检测结果 Fig. 7 Edge detection with different algorithms |
从图 7(c)和7(d)中可以看出，本文方法的检测结果达到了较为理想的视觉效果。建筑物与背景间边缘，以及栅栏与建筑物间边缘在图 7(c)中得到了很好检测；图 7(d)中的高楼与人的轮廓边缘也得到了很好的检测，同时高楼内部的边缘纹理也有理想的检测效果。
对于Canny方法而言，由于方法本身加入了低通滤波，使得边缘更加平滑，但是对定位精度会有一定的影响。同时该算法中的非极大值抑制方法会产生过多的伪边缘，如图 7(e)右下角部分与图 7(f)中的树冠内部的纹理；Sobel方法在图 7(g)中获得了很好的检测效果，但是破碎点过多，同时图 7(h)中漏检测了部分高楼内部纹理边缘；相位一致模型的边缘检测结果(图 7(i)和7(j))同样具有破碎现象，并且边缘断裂严重；由于数学形态学方法不能获得边缘梯度方向，而无法利用非极大值抑制方法得到单像元边缘。所以采用阈值法，显然部分边界未能很好地检测，同时容易产生过多悬挂线等伪边缘(图 7(k)和7(l))。
利用文献[25]提出的边缘检测评价测度FOM(figure of merit)，分析4种算法的边缘检测效果
式中，N1和ND分别为真实和检测边缘数目；di为第i个边缘点与相应的检测边缘间的距离；ζ为尺度常量，设为1/9。
当FOM值越大时，图像边缘定位效果越好。根据公式(23)以及标准参考图，不同方法的FOM值如表 7所示。可明显看出，本文方法较于其他方法具有更好的定位效果。
| 本文方法 | Canny | Sobel | 相位一致模型 | 数学形态学 | |
| 图 8(a) | 0.678 1 | 0.585 3 | 0.630 5 | 0.519 5 | 0.562 2 |
| 图 8(b) | 0.667 4 | 0.607 3 | 0.353 1 | 0.593 0 | 0.634 1 |
此外，为验证该方法的应用的广泛性，选取WorldView (图 8(a))和QuickBird (图 8(b)) 高分辨率遥感图像进行边缘检测试验。图 8(a)中主要地物有水塘、建筑物、操场跑道等，而图 8(b)则包括了厂房等地物。图像大小分别为500×500和512×512。图 8(c)和8(d)给出了边缘幅度图，滤波器参数设置为α=0.25，β=1.90，θ=π。
|
| 图 8 高空间分辨率遥感图像的边缘检测结果 Fig. 8 Edge detection of high resolution remotely sensed imagery |
结果显示，图 8(a)中水塘、建筑物、操场跑道等地物边缘响应明显，尤其是操场边缘响应整体一致，充分显示出本文全方向边缘检测方法的优势；图 8(b)中厂房的边缘响应明显。同时检测时间稳定在0.033 s左右。因此，本文方法也可有效地用于高分辨率遥感图像的边缘检测工作。
表 8给出不同方法下两幅高遥感图像边缘检测结果的FOM值，可以明显看出本文具有更好的边缘检测效果。
| 本文方法 | Canny | Sobel | 相位一致模型 | 数学形态学 | |
| 图 8(c) | 0.615 6 | 0.578 0 | 0.260 5 | 0.487 1 | 0.419 1 |
| 图 8(d) | 0.601 7 | 0.574 0 | 0.264 7 | 0.538 5 | 0.565 9 |
为了更好地表现图像边缘特征的高频特性，本文以一维Butterworth滤波器为基础，将具有非线性识别机制的Log函数引入该滤波器。进而提出基于二维Log Butterworth滤波器的全方向图像边缘检测的频域实现方法。图像行列数不一致时，中心频率分布于椭圆之上，椭圆的长短轴之比为图像长宽比。为了方便滤波器参数选取，滤波器中心频率被归一化为α。同时最大上限截止频率与中心频率之比用取值范围为(1,2)的参数β代替。并借助F-measure和PSNR来实现滤波器最优参数范围的选定。
从本文方法和不同窗口大小的Canny算子的乘法和加法次数之比对来看，在加法次数上，本文方法并不占很大优势，但在乘法次数占有绝对优势，因此本文方法优于Canny算子等传统边缘检测方法。以不同大小图像作为试验数据进行边缘检测的结果也表明本文方法的运算速度更高。
对比分析本文方法与Canny、Sobel、相位一致模型和数学形态学方法的边缘检测效果的最终的FOM值显示出本文方法具有更好的定位效果。由于本文方法具有全方向边缘检测特点，因此检测结果中的边缘响应程度均不受边缘方向影响。此外，该方法不受图像行列数比例的影响，可以有效地应用于图像边缘检测。
| [1] | ROBERTS L G. Machine Perception of Three Dimensional Solids[R]. Massachusetts: Massachusetts Institute of Technology,1965. |
| [2] | SOBEL I. Neighbourhood Coding of Binary Images for Fast Contour Following and General Array Binary Processing[J]. Computer Graphics and Image Processing, 1978, 8: 127-135. |
| [3] | PREWITT J M S. Object Enhancement and Extraction[J]. Picture Processing and Psychopictorics, 1970,10(1):15-19. |
| [4] | CANNY J F. A Computational Approach to Edge Detection[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1986, 8 (6): 679-698. |
| [5] | MARR D, HILDRETH E C. Theory of Edge Detection[J]. The Royal Society, 1980, 207:187-217. |
| [6] | MALLAT S, ZHONG S. Characterization of Signals from Multiscale Edges[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1992, 14(7): 710-731. |
| [7] | JIANG W, LAM K M, SHEN T Z. Efficient Edge Detection Using Simplified Gabor Wavelets[J]. IEEE Transactions on Systems, 2009, 39(4): 1036-1047. |
| [8] | MORRONE M C, OWENS R A. Feature Detection from Local Energy[J]. Pattern Recognition Letters,1987, 6(5):303-313. |
| [9] | KOVESI P. Phase Congruency: A Low-level Image Invariant[J]. Psychological Research,2000, 64(2): 136-148. |
| [10] | WANG Ke, XIAO Pengfeng. The Algorithm of Image Features Detection from Phase Congruency Model Based on 2-D Hilbert Transform[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(6):605-611. (王珂, 肖鹏峰. 基于二维希尔伯特变换的相位一致模型图像特征检测方法[J]. 测绘学报, 2010, 39(6):605-611.) |
| [11] | MOHAMED A, TAREK A. New Edge Detection Technique Based on the Shannon Entropy in Gray Level Images[J]. International Journal on Computer Science and Engineering, 2011, 3(6): 2224-2232. |
| [12] | RAKESH R, CHAUDHURI P, MURTHY C A. Thresholding in Edge Detection: A Statistical Approach[J]. IEEE Transactions on Image Processing, 2004, 13:927-936. |
| [13] | KIM D S, LEE W H, KWEON L S. Automatic Edge Detection Using 3×3 Ideal Binary Pixel Patterns and Fuzzy-based Edge Thresholding[J].Pattern Recognition Letters, 2004, 25: 101-106. |
| [14] | FANG B, TANG Y Y. Improved Class Statistics Estimation for Sparse Data Problems in Offline Signature Verification[J]. IEEE Transactions on Systems, 2005, 35 (3): 276-286. |
| [15] | KARANDE K J, TALBAR S N. Independent Component Analysis of Edge Information for Face Recognition[J]. International Journal of Image Processing, 2009, 3(3): 120-130. |
| [16] | HUANG J, YOU X G, TANG Y Y, et al. A Novel Iris Segmentation Using Radial-suppression Edge Detection[J]. Signal Processing,2009, 89(4):2630-2643. |
| [17] | GUDMUNDSSON M, EI-KWAE E A, KABUKA M R. Edge Detection in Medical Images Using a Genetic Algorithm[J]. IEEE Transaction on Medical Imaging, 1998, 17(3):469-474. |
| [18] | BANDYOPADHYAY S K. Edge Detection in Brain Images[J]. International Journal of Computer Science and Information Technologies, 2011, 2(2):884-887. |
| [19] | TELLO M, LOPEZ M C, MALLORQUI J J, et al. Edge Enhancement Algorithm Based on the Wavelet Transform for Automatic Edge Detection in SAR Image[J]. IEEE Transaction on Geoscience and Remote Sensing, 2011, 49(1):222-235. |
| [20] | LIU H, JEZEK K C. Automated Extraction of Coastline from Satellite Imagery by Integrating Canny Edge Detection and Locally Adaptive Thresholding Methods[J]. International Journal of Remote Sensing, 2004, 25(5):937-958. |
| [21] | XIAO Pengfeng, FENG Xuezhi, ZHAO Shuhe, et al. Segmentation of High-resolution Remotely Sensed Imagery Based on Phase Congruency[J]. Acta Geodaetica et Cartographica Sinica, 2007, 36(2): 146-152. (肖鹏峰, 冯学智, 赵书河, 等. 基于相位一致的高分辨率遥感图像分割方法[J]. 测绘学报, 2007, 36(2): 146 -152.) |
| [22] | POLLEN D A, RONNER S E. Visual Cortical Neurons as Localized Spatial Frequency Filters[J]. IEEE Transactions on System, Man and Cybernetics, 1983, 13(5): 907-916. |
| [23] | FIELD D J. Relations between the Statistics of Natural Images and the Response Properties of Cortical Cells[J]. Journal of the Optical Society of America , 1987, 4(12): 2379-2394. |
| [24] | RIJSBERGEN C V. Information Retrieval[R]. Glasgow: University of Glasgow, 1979. |
| [25] | ABDOU I E, PRATT W K. Quantitative Design and Evaluation of Enhancement/ Thresholding Edge Detectors[J]. Proceedings of the IEEE ,1979, 67(5):753-763. |