canny检测出的多个边缘合成一个_Canny边缘检测算法

PS：计算机视觉课程的一项大作业，作为门外汉，在网上查找资料学习和复现时踩了一些坑，在完成作业后分享一下，希望后面的小白同学少走弯路。

1.问题描述

在计算机视觉/机器视觉领域，图像分割的应用十分普遍，它是指将数字图像细分为多个图像子区域的过程。图像分割的目的是简化或改变图像的表示形式，使得图像更容易理解和分析。它的基本形式通常为定位图像中的物体和边界，结果是图像上子区域或轮廓线的集合。图像分割方法主要包括阈值处理(二值化)、聚类法、边缘检测和区域生长等。对于图像分割问题的求解没有统一范式，通常要与领域知识充分结合，才能更为有效地解决。

边缘检测是基于灰度突变来分割图像的常用方法，其实质是提取图像中不连续部分的特征。目前常见边缘检测算子有差分算子、Roberts算子、Sobel算子、Prewitt算子、Log算子以及Canny算子等。其中，Canny算子是由计算机科学家John F. Canny于1986年提出的一种边缘检测算子，是目前理论上相对最完善的一种边缘检测算法。在MATLAB、OpenCV等常用图像处理工具中已有内置的Canny算子API，本次大作业将依据Canny算子的算法原理在MATLAB环境中进行复现，并与MATLAB内置Canny算子的效果进行对比。

2.算法原理

基于Canny算子的边缘检测主要有5个步骤，依次是高斯滤波、像素梯度计算、非极大值抑制、滞后阈值处理和孤立弱边缘抑制。

(1) 高斯滤波

高斯滤波使用的高斯核是具有

和

两个维度的高斯函数，且两个维度上标准差一般取相同，形式为：

高斯滤波，即使用即使用某一尺寸的二维高斯核与图像进行卷积。由于数字图像的数据形式为离散矩阵，高斯核是对连续高斯函数的离散近似，通过对高斯曲面进行离散采样和归一化得出。例如，尺寸为

，标准差为

的高斯核为：

在确定高斯核后，将其与图像进行离散卷积即可。

(2) 使用Sobel算子计算像素梯度

Sobel算子是两个

的矩阵，分别为

和

。前者用于计算图像

方向像素梯度矩阵

，后者用于计算图像

方向像素梯度矩阵

。具体形式为：

其中，

为灰度图像矩阵，且此处的

表示互相关运算(卷积运算可视为将卷积核旋转180°后的互相关运算)。需要说明的是，图像矩阵坐标系原点在左上角，且

正方向为从左到右，

正方向为从上到下。则由

可计算得到梯度强度矩阵

。

(3) 非极大值像素梯度抑制

非极大值像素梯度抑制的目的在于消除边缘检测带来的杂散响应，起到将边缘“瘦身”的作用。其基本方法是将当前像素梯度强度与沿正负梯度方向上的相邻像素的梯度强度进行比较，若其最大(即为极值)，则保留该像素为边缘点，若不是最大，则对其进行抑制，不将其作为边缘点。为了更精确计算，通常在跨越梯度方向的两个相邻像素之间使用线性插值来得到要参与比较的像素梯度。图1. 像素梯度方向线性插值示意图

如图1所示，可将像素的邻接情况划分为4个区域，其中每个区域包含上下两部分。若中心像素点

方向梯度强度为

，

方向梯度强度为

，梯度强度为

，则根据

和

的正负和大小可判断出其梯度方向所属区域，进而根据其像素进而根据其像素梯度方向以及相邻点的像素梯度线性插值得到正负梯度方向的两个参与比较的梯度强度

和

。公式如下：

其他三个区域的计算方法类似。需要注意的是，若

，说明像素点无像素梯度，其为非边缘。

(4) 阈值滞后处理

定义一个高阈值和一个低阈值。梯度强度低于低阈值的像素点被抑制，不作为边缘点；高于高阈值的像素点被定义为强边缘，保留为边缘点；处于高低阈值之间的定义为弱边缘，留待进一步处理。

(5) 孤立弱边缘抑制

通常而言，由真实边缘引起的弱边缘像素点将连接到强边缘像素点，而噪声响应则未连接。通过查看弱边缘像素及其8个邻域像素，可根据其与强边缘的连接情况来进行判断。一般，可定义只要其中邻域像素其中一个为强边缘像素点，则该弱边缘就可以保留为强边缘，即真实边缘点。

3.程序流程图图2. 程序流程图

4.实验结果

在MATLAB R2016a环境下对实现了基于Canny算子的边缘检测，并与MATLAB内置的Canny边缘检测API进行了四组对比实验。其中“Lena”和“房屋”为计算机视觉领域的常用测试图片，“机翼”为出差时用手机自行拍摄的照片，“小刘鸭”为微信表情包图片。图3.实验结果-Lena

、图4.实验结果-房屋图5.实验结果-机翼图6.实验结果-小刘鸭

为了达到较好的检测效果，需要反复调整自行编写的边缘检测程序中的滤波器参数和滞后阈值处理时的高阈值和低阈值。MATLAB的Canny边缘检测API则采用其默认参数，猜测其中高、低阈值可能由自适应算法自行计算。图3和图4的结果表明，相较API，自行编写的Canny检测程序提取到了更多的细节信息，但也保留了较多的局部噪点；API提取的轮廓特征更为清晰，但局部失真情况相对较重。图5结果表明，相较于API，自行编写的程序提取到了更为清晰的机翼轮廓，并保留了较少的云层的纹理。图6效果相差不大。从实验过程可知Canny 算子对于滤波参数和高、低阈值的选取还是较为敏感，使得实际应用过程中需要反复调试。囿于时间和能力，没有对参数的自适应选择展开进一步学习和尝试。

5.源代码

clc;

clear;

%--------------------------------------------------%图片导入与预处理

% f_original = imread('lena.png'); %导入lena图片

% f_original = imread('house.png'); %导入房屋图片

% f_original = imread('wing.tif'); %导入机翼图片

f_original = imread('xiaoliuya.jpg'); %导入小刘鸭

f_grey = rgb2gray(f_original); %转换为灰度图像

%------------------------------------------------- %图片导入与预处理

%--------------------------------------------------%自行编写的CANNY算子边缘检测

%1-高斯滤波

gw = fspecial('gaussian',[5,5],0.5); %高斯滤波设置核，5*5，标准差为0.5

f_filter = imfilter(f_grey,gw,'replicate'); %高斯滤波

f = f_filter;

%2-利用Sobel算子计算像素梯度

Sobel_X = [-1,0,1;-2,0,2;-1,0,1]; %X方向Sobel算子(互相关算子，非卷积算子)

Sobel_Y = [-1,-2,-1;0,0,0;1,2,1]; %Y方向Sobel算子(互相关算子)

[rowNum,columnNum] = size(f);

f_extend = zeros(rowNum+2,columnNum+2); %图像扩充，边界补充为0

for i = 2:rowNum+1

for j = 2:columnNum+1

f_extend(i,j) = f(i-1,j-1);

end

end

Gx = zeros(rowNum,columnNum);

Gy = zeros(rowNum,columnNum);

for i = 2:rowNum+1 %计算x向和y向梯度

for j = 2:columnNum+1

window = [f_extend(i-1,j-1),f_extend(i-1,j),f_extend(i-1,j+1);...

f_extend(i,j-1),f_extend(i,j),f_extend(i,j+1);...

f_extend(i+1,j-1),f_extend(i+1,j),f_extend(i+1,j+1)];

Gx(i-1,j-1) = sum(sum(Sobel_X .* window)); %计算x向梯度

Gy(i-1,j-1) = sum(sum(Sobel_Y .* window)); %计算y向梯度

end

end

Sxy = sqrt(Gx.*Gx + Gy.*Gy); %梯度强度矩阵计算

%3-非极大值抑制

indexD = zeros(rowNum,columnNum);

for i = 1:rowNum %判断梯度方向所属区间，Gx=Gy=0，则令其为5，肯定不是边界点

for j = 1:columnNum

ix = Gx(i,j);

iy = Gy(i,j);

if (iy<=0 && ix>-iy) || (iy>=0 && ix

indexD(i,j) = 1;

elseif (ix>0 && ix<=-iy) || (ix<0 && ix>=-iy) %梯度方向属于区间2

indexD(i,j) = 2;

elseif (ix<=0 && ix>iy) || (ix>=0 && ix

indexD(i,j) = 3;

elseif (iy<0 && ix<=iy) || (iy>0 && ix>=iy) %梯度方向属于区间4

indexD(i,j) = 4;

else %Gx和Gy均为0，无梯度，肯定非边缘

indexD(i,j) = 5;

end

end

end

Gup = zeros(rowNum,columnNum);

Gdown = zeros(rowNum,columnNum);

for i = 2:rowNum-1 %计算非边界处的插值梯度强度

for j = 2:columnNum-1

ix = Gx(i,j);

iy = Gy(i,j);

if indexD(i,j) == 1 %计算区间1内插值梯度，Gup为上方区间的梯度，Gdown为下方区间的梯度

t = abs(iy./ix);

Gup(i,j) = Sxy(i,j+1).*(1-t) + Sxy(i-1,j+1).*t;

Gdown(i,j) = Sxy(i,j-1).*(1-t) + Sxy(i+1,j-1).*t;

elseif indexD(i,j) == 2 %计算区间2内插值梯度

t = abs(ix./iy);

Gup(i,j) = Sxy(i-1,j).*(1-t) + Sxy(i-1,j+1).*t;

Gdown(i,j) = Sxy(i+1,j).*(1-t) + Sxy(i+1,j-1).*t;

elseif indexD(i,j) == 3 %计算区间3内插值梯度

t = abs(ix./iy);

Gup(i,j) = Sxy(i-1,j).*(1-t) + Sxy(i-1,j-1).*t;

Gdown(i,j) = Sxy(i+1,j).*(1-t) + Sxy(i+1,j+1).*t;

elseif indexD(i,j) == 4 %计算区间4内插值梯度

t = abs(iy./ix);

Gup(i,j) = Sxy(i,j-1).*(1-t) + Sxy(i-1,j-1).*t;

Gdown(i,j) = Sxy(i,j+1).*(1-t) + Sxy(i+1,j+1).*t;

end

end

end

Sxy_NMX = zeros(rowNum,columnNum); %判断是否为梯度方向极大值

for i = 1:rowNum

for j = 1:columnNum

if Sxy(i,j) >= Gup(i,j) && Sxy(i,j) >= Gdown(i,j) %若为梯度方向极大值，则保留；

Sxy_NMX(i,j) = Sxy(i,j); %否则，进行抑制(置0)

end

end

end

%4-滞后阈值法+5-抑制孤立的弱边缘

f_final = zeros(rowNum,columnNum);

%Tl为高阈值，Th为低阈值，connectNum为联通参数，一般为1

% Tl = 60; %lena

% Th = 120;

% connectNum = 1;

% Tl = 40; %房屋

% Th = 120;

% % connectNum = 1;

% Tl = 22; %机翼

% Th = 40;

% connectNum = 1;

Tl = 15; %小刘鸭

Th = 35;

connectNum = 1;

for i = 2:rowNum-1

for j = 2:columnNum-1

if Sxy_NMX(i,j) >= Th %高于高阈值的像素为强边缘

f_final(i,j) = 1;

elseif Sxy_NMX(i,j) <= Tl %低于低阈值的像素为非边缘

f_final(i,j) = 0;

else %位于高低阈值之间的像素为弱边缘，进行孤立性检测

count = 0;

if Sxy_NMX(i-1,j-1)~=0 %左上方像素

count = count+1;

end

if Sxy_NMX(i-1,j)~=0 %上方像素

count = count+1;

end

if Sxy_NMX(i-1,j+1)~=0 %右上方像素

count = count+1;

end

if Sxy_NMX(i,j-1)~=0 %左方像素

count = count+1;

end

if Sxy_NMX(i,j+1)~=0 %右方像素

count = count+1;

end

if Sxy_NMX(i+1,j-1)~=0 %左下方像素

count = count+1;

end

if Sxy_NMX(i+1,j)~=0 %下方像素

count = count+1;

end

if Sxy_NMX(i+1,j+1)~=0 %右下方像素

count = count+1;

end

if count >= connectNum %弱边缘非孤立，则为边缘

f_final(i,j) = 1;

end

end

end

end

%--------------------------------------------------%自行编写的CANNY算子边缘检测

%--------------------------------------------------%MATLAB的CANNY算子边缘检测

[fCanny_dafault,tc] = edge(f_grey,'canny'); %使用默认参数

%--------------------------------------------------%MATLAB的CANNY算子边缘检测

%--------------------------------------------------%测试

subplot(1,3,1);

imshow(f_original);

title('原始图像','fontsize',20);

subplot(1,3,2);

imshow(f_final);

title('自行编写的Canny边缘检测','fontsize',20);

subplot(1,3,3);

imshow(fCanny_dafault);

title('MATLAB内置的Canny边缘检测','fontsize',20);

%--------------------------------------------------%测试

6.References

[3] 维基百科