原创文章，欢迎转载。转载请注明：转载自 祥的博客
原文链接：http://blog.csdn.net/humanking7/article/details/46606791
图象的边缘是指图象局部区域亮度变化显著的部分，该区域的灰度剖面一般可以看作是一个阶跃，既从一个灰度值在很小的缓冲区域内急剧变化到另一个灰度相差较大的灰度值。
(1) 必须满足两个条件：①能有效地抑制噪声；②必须尽量精确确定边缘的位置。
(2) 根据对信噪比与定位乘积进行测度，得到最优化逼近算子。这就是Canny边缘检测算子。
(3) 类似与Marr（LoG）边缘检测方法，也属于先平滑后求导数的方法。
步骤1：用高斯滤波器平滑处理原图像；
步骤2：用一阶偏导的有限差分进行计算梯度的幅值和方向；
步骤3：对梯度幅值进行非极大值抑制；
步骤4：用双阈值算法检测和连接边缘。
使用平滑滤波的原因从根本上来说是边缘检测算子中的导数计算。导数计算对噪声十分敏感，如果不提前使用滤波器加以改善，则在导数计算后，噪声将会被放大，使得检测出来的虚假边缘变多，不利于边缘的提取。
平滑滤波和边缘检测是一对矛盾的概念。一方面，平滑滤波能够有效的抑制噪声，而在此过程中会使得图像边缘模糊，增加了边缘定位的不确定性。另一方面，平滑滤波能够除去对边缘检测中导数运算敏感的噪声，有效的抑制了虚假边缘的产生。实际工程经验表明，高斯滤波器可以在抗噪声干扰和边缘检测精确定位之间提供一个较好的折中方案。
图像的边缘有方向和幅度两个属性，沿边缘方向像素变化平缓，垂直于边缘方向像素变化剧烈，边缘上的这种变化可以用微分算子检测出来，通常用一阶或二阶导数来检测边缘.
由上图可以看出，一阶导数可以用于检测图像中的一个点是否是边缘点（也就是判断一个点是否在斜坡上）。同样，二阶导数的符号可以用于判断一个边缘像素是否在亮的一边还是暗的一边。
用一阶偏导的有限差分来计算梯度的幅值和方向。
下图中，图a经过梯度计算后的得到梯度三维示意图b。图b中x和y代表图像像素位置，竖轴的数值反映了梯度幅值的大小。由公式(3-17)可知，边缘方向与梯度方向是相垂直的，如图b所示。
主函数代码文件main.m
clear all;
clear
clc;
%读进图像
[filename, pathname] = uigetfile({‘*.jpg‘; ‘*.bmp‘; ‘*.gif‘}, ‘选择图片‘);
%没有图像
if filename == 0
return;
end
imgsrc = imread([pathname, filename]);
[y, x, dim] = size(imgsrc);
%转换为灰度图
if dim>1
imgsrc = rgb2gray(imgsrc);
end
sigma = 1;
gausFilter = fspecial(‘gaussian‘, [3,3], sigma);
img= imfilter(imgsrc, gausFilter, ‘replicate‘);
zz = double(img);
%----------------------------------------------------------
%自己的边缘检测函数
[m theta sector canny1 canny2 bin] = canny1step(img, 22);
[msrc thetasrc sectorsrc c1src c2src binsrc] = canny1step(imgsrc, 22);
%Matlab自带的边缘检测
ed = edge(img, ‘canny‘, 0.5);
[xx, yy] = meshgrid(1:x, 1:y);
figure(1)
%mesh(yy, xx, zz);
surf(yy, xx, zz);
xlabel(‘y‘);
ylabel(‘x‘);
zlabel(‘Grayscale‘);
axis tight
figure(2)
subplot(4,2,1);
imshow(imgsrc);%原图
subplot(4,2,2);
imshow(img);%高斯滤波后
subplot(4,2,3);
imshow(uint8(m));%导数
subplot(4,2,4);
imshow(uint8(canny1));%非极大值抑制
subplot(4,2,5);
imshow(uint8(canny2));%双阈值
subplot(4,2,6);
imshow(ed);%Matlab自带边缘检测
subplot(4,2,8);
imshow(bin);%我自己的bin
figure(3)
edzz = 255*double(ed);
mesh(yy,xx,edzz);
xlabel(‘y‘);
ylabel(‘x‘);
zlabel(‘Grayscale‘);
axis tight
figure(4)
mesh(yy,xx,m);%画偏导数
xlabel(‘y‘);
ylabel(‘x‘);
zlabel(‘Derivative‘);
axis tight
figure(5)
mesh(yy,xx,theta);
xlabel(‘y‘);
ylabel(‘x‘);
zlabel(‘Theta‘);
axis tight
figure(6)
mesh(yy,xx,sector);
xlabel(‘y‘);
ylabel(‘x‘);
zlabel(‘Sector‘);
axis tight
figure(7)
mesh(yy,xx,canny2);
xlabel(‘y‘);
ylabel(‘x‘);
zlabel(‘Sector‘);
axis tight
Canny边缘检测函数文件canny1step.m
function [ m, theta, sector, canny1, canny2, bin] = canny1step( src, lowTh)
%canny函数第一步，求去x，y方向的偏导，模板如下：
% Gx
% 1 -1
% 1 -1
% Gy
% -1 -1
% 1 1
%------------------------------------
% 输入：
% src：图像，如果不是灰度图转成灰度图
% lowTh：低阈值
% 输出：
% m： 两个偏导的平方差，反映了边缘的强度
% theta：反映了边缘的方向
% sector：将方向分为3个区域，具体如下
% 2 1 0
% 3 X 3
% 0 1 2
% canny1：非极大值
% canny2：双阈值抑制
% bin ： 二值化
%---------------------------------------
[Ay Ax dim ] = size(src);
%转换为灰度图
if dim>1
src = rgb2gray(src);
end
src = double(src);
m = zeros(Ay, Ax);
theta = zeros(Ay, Ax);
sector = zeros(Ay, Ax);
canny1 = zeros(Ay, Ax);%非极大值抑制
canny2 = zeros(Ay, Ax);%双阈值检测和连接
bin = zeros(Ay, Ax);
for y = 1:(Ay-1)
for x = 1:(Ax-1)
gx = src(y, x) + src(y+1, x) - src(y, x+1) - src(y+1, x+1);
gy = -src(y, x) + src(y+1, x) - src(y, x+1) + src(y+1, x+1);
m(y,x) = (gx^2+gy^2)^0.5 ;
%--------------------------------
theta(y,x) = atand(gx/gy) ;
tem = theta(y,x);
%--------------------------------
if (tem<67.5)&&(tem>22.5)
sector(y,x) = 0;
elseif (tem<22.5)&&(tem>-22.5)
sector(y,x) = 3;
elseif (tem<-22.5)&&(tem>-67.5)
sector(y,x) = 2;
else
sector(y,x) = 1;
end
%--------------------------------
end
end
%-------------------------
%非极大值抑制
%------> x
% 2 1 0
% 3 X 3
%y 0 1 2
for y = 2:(Ay-1)
for x = 2:(Ax-1)
if 0 == sector(y,x) %右上 - 左下
if ( m(y,x)>m(y-1,x+1) )&&( m(y,x)>m(y+1,x-1) )
canny1(y,x) = m(y,x);
else
canny1(y,x) = 0;
end
elseif 1 == sector(y,x) %竖直方向
if ( m(y,x)>m(y-1,x) )&&( m(y,x)>m(y+1,x) )
canny1(y,x) = m(y,x);
else
canny1(y,x) = 0;
end
elseif 2 == sector(y,x) %左上 - 右下
if ( m(y,x)>m(y-1,x-1) )&&( m(y,x)>m(y+1,x+1) )
canny1(y,x) = m(y,x);
else
canny1(y,x) = 0;
end
elseif 3 == sector(y,x) %横方向
if ( m(y,x)>m(y,x+1) )&&( m(y,x)>m(y,x-1) )
canny1(y,x) = m(y,x);
else
canny1(y,x) = 0;
end
end
end%end for x
end%end for y
%---------------------------------
%双阈值检测
ratio = 2;
for y = 2:(Ay-1)
for x = 2:(Ax-1)
if canny1(y,x)<lowTh %低阈值处理
canny2(y,x) = 0;
bin(y,x) = 0;
continue;
elseif canny1(y,x)>ratio*lowTh %高阈值处理
canny2(y,x) = canny1(y,x);
bin(y,x) = 1;
continue;
else %介于之间的看其8领域有没有高于高阈值的，有则可以为边缘
tem =[canny1(y-1,x-1), canny1(y-1,x), canny1(y-1,x+1);
canny1(y,x-1), canny1(y,x), canny1(y,x+1);
canny1(y+1,x-1), canny1(y+1,x), canny1(y+1,x+1)];
temMax = max(tem);
if temMax(1) > ratio*lowTh
canny2(y,x) = temMax(1);
bin(y,x) = 1;
continue;
else
canny2(y,x) = 0;
bin(y,x) = 0;
continue;
end
end
end%end for x
end%end for y
end%end of function
图2所示的四张图片展示了整个边缘检测的过程。
图3展示了梯度幅值在非极大值抑制前后的三维效果，对应于图2中的图(a)和图(b)。
注： 本文所用的图像和代码源自于本人的毕设论文。
原文地址：http://blog.csdn.net/humanking7/article/details/46606791
