canny边缘检测学习笔记

Canny边缘检测的步骤：

1、高斯平滑（gaussian smoothing）

2、图像求导（derivation）

3、非极大值抑制（non-local maxima suppression）

4、滞后门限（hysteresis threshold）

* 非极大值抑制 和 滞后门限 是canny边缘检测的特色步骤

具体介绍：

http://blog.csdn.net/likezhaobin/article/details/6892176

http://www.cnblogs.com/blue-lg/archive/2011/12/25/2300899.html

算法实现（matlab源代码）：

matlab中边缘检测的函数是edge(gray_image, 'canny');我把里面关于canny算法的部分提出来写成一个单独的函数。

 

0、首先是预处理

变量a是一个m*n的矩阵，是输入的灰度图。变量e是用来保存边缘的矩阵，是一个二值图。

function e = canny_edge_matlab(a)
% Transform to a double precision intensity image if necessary
if ~isa(a,'double') && ~isa(a,'single') 
  a = im2single(a);
end
[m,n] = size(a);
% The output edge map:
e = false(m,n);

 

1、高斯平滑

GaussianDieOff是高斯衰变值，小于该值的高斯函数值都视为0.

PercentOfPixelsNotEdges代表非边缘像素占总像素个数的百分比，默认为70%。ThresholdRatio是一个系数，lowThreshhold=ThresholdRatio*highThreshold。

sigma是高斯函数里的标准差，默认为1。高斯函数的模板大小最大为60*60，实际大小由GaussianDieOff和sigma共同决定

gau保存了1D的高斯模板，用来做高斯平滑。分别沿x和y轴做平滑

dgau2D保存了2D高斯函数对x的一阶偏导，为了做图像微分。

% Magic numbers
GaussianDieOff = .0001;  
PercentOfPixelsNotEdges = .7; % Used for selecting thresholds
ThresholdRatio = .4;          % Low thresh is this fraction of the high.
% Design the filters - a gaussian and its derivative
sigma=1; % sigma是一个入参，默认是1，这句是我加的，为了代码可以正常运行
pw = 1:30; % possible widths
ssq = sigma^2;

% 为了确定高斯核大小，当函数值小于衰减系数时，就把他视为0
% 窗口大小为power(2*width)
width = find(exp(-(pw.*pw)/(2*ssq))>GaussianDieOff,1,'last');
if isempty(width)
width = 1;  % the user entered a really small sigma
end
t = (-width:width);
gau = exp(-(t.*t)/(2*ssq))/(2*pi*ssq);     % the gaussian 1D filter

% Find the directional derivative of 2D Gaussian (along X-axis)
% Since the result is symmetric along X, we can get the derivative along
% Y-axis simply by transposing the result for X direction.
[x,y]=meshgrid(-width:width,-width:width);
dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq))/(pi*ssq);

% Convolve the filters with the image in each direction
% The canny edge detector first requires convolution with
% 2D gaussian, and then with the derivitave of a gaussian.
% Since gaussian filter is separable, for smoothing, we can use 
% two 1D convolutions in order to achieve the effect of convolving
% with 2D Gaussian.  We convolve along rows and then columns.

% smooth the image out
% 高斯模糊
aSmooth=imfilter(a,gau,'conv','replicate');   % run the filter across rows
aSmooth=imfilter(aSmooth,gau','conv','replicate'); % and then across columns

 

2、图像微分

dgau2D是高斯函数对x的偏导，其转置dgau2D‘是高斯函数对y的偏导。

ax是图像在x方向上与高斯微分卷积的结果（x方向上的图像微分），ay是y方向的图像微分。用ax和ay可以求出图像各个像素的梯度值和梯度方向。mag是梯度值。最后要对mag进行归一化。

%apply directional derivatives
%一阶高斯
ax = imfilter(aSmooth, dgau2D, 'conv','replicate');
ay = imfilter(aSmooth, dgau2D', 'conv','replicate');
mag = sqrt((ax.*ax) + (ay.*ay));

 

3、非极大值抑制

• 如果用户指定了阈值，就用用户指定的值，否则，用默认值计算阈值。阈值的计算方法是：

1）将图像按灰度值转化为64阶的直方图。

2）计算高阈值。高阈值就是图像中边缘像素灰度值的下限，这里默认有70%的像素不是边缘，现在要求出这70%的像素的灰度值水平。例如：

图片共45个像素，做成7阶直方图，70%像素不是边缘，即31个像素不是边缘，所以高阈值要淘汰掉这些像素。由图可知，前5阶像素个数和为31，所以将灰度值的水平设置为6，则70%的像素就被淘汰了。由于阈值取值范围是（0，1），所以该图片的阈值就是6/7。

3）计算低阈值。计算公式是lowThresh = ThresholdRatio*highThresh;

% Select the thresholds
thresh = ''; % 我自己加的，为了程序可以进入if
if isempty(thresh) 
  % 这里的直方图阶数可任意，理论上大于1小于255都可以，64是个比较好的值
counts=imhist(mag, 64);
highThresh = find(cumsum(counts) > PercentOfPixelsNotEdges*m*n,...
                  1,'first') / 64;
lowThresh = ThresholdRatio*highThresh;
thresh = [lowThresh highThresh];
% 用户指定的阈值，我把他注释掉了
%   elseif length(thresh)==1
%     highThresh = thresh;
%     if thresh>=1
%       eid = sprintf('Images:%s:thresholdMustBeLessThanOne', mfilename);
%       msg = 'The threshold must be less than 1.'; 
%       error(eid,'%s',msg);
%     end
%     lowThresh = ThresholdRatio*thresh;
%     thresh = [lowThresh highThresh];
%   elseif length(thresh)==2
%     lowThresh = thresh(1);
%     highThresh = thresh(2);
%     if (lowThresh >= highThresh) || (highThresh >= 1)
%       eid = sprintf('Images:%s:thresholdOutOfRange', mfilename);
%       msg = 'Thresh must be [low high], where low < high < 1.'; 
%       error(eid,'%s',msg);
%     end
end

 

• 非极大值抑制主要使用了cannyFindLocalMaxima()这个函数来查找梯度方向上，与该像素相邻的区域内梯度的最大值。其实就是线性插值。这里要用到梯度的方向。

根据对称性，只需要分别求出4个方向区间的极大值。idxWeak保存了弱边缘的像素编号，idxStrong保存了强边缘的像素编号。代码如下：

% The next step is to do the non-maximum suppression.  
% We will accrue indices which specify ON pixels in strong edgemap
% The array e will become the weak edge map.
idxStrong = [];  
for dir = 1:4
idxLocalMax = cannyFindLocalMaxima(dir,ax,ay,mag);
idxWeak = idxLocalMax(mag(idxLocalMax) > lowThresh);
e(idxWeak)=1; % 先把弱边缘视为边缘，后面将用迟滞算法筛选与强边缘相连的弱边缘
idxStrong = [idxStrong; idxWeak(mag(idxWeak) > highThresh)]; %#ok<AGROW>
end

 

• 某个方向上的极大值查找主要使用了Local Function : cannyFindLocalMaxima。如果该像素是其梯度方向上的局部极大值，则返回该像素编号。注释写在代码中了：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Local Function : cannyFindLocalMaxima
%
function idxLocalMax = cannyFindLocalMaxima(direction,ix,iy,mag)
%
% This sub-function helps with the non-maximum suppression in the Canny
% edge detector.  The input parameters are:
% 
%   direction - the index of which direction the gradient is pointing, 
%               read from the diagram below. direction is 1, 2, 3, or 4.
%   ix        - input image filtered by derivative of gaussian along x 
%   iy        - input image filtered by derivative of gaussian along y
%   mag       - the gradient magnitude image
%
%    there are 4 cases:
%
%                         The X marks the pixel in question, and each
%         3     2         of the quadrants for the gradient vector
%       O----0----0       fall into two cases, divided by the 45 
%     4 |         | 1     degree line.  In one case the gradient
%       |         |       vector is more horizontal, and in the other
%       O    X    O       it is more vertical.  There are eight 
%       |         |       divisions, but for the non-maximum suppression  
%    (1)|         |(4)    we are only worried about 4 of them since we 
%       O----O----O       use symmetric points about the center pixel.
%        (2)   (3)        


[m,n] = size(mag);

% Find the indices of all points whose gradient (specified by the 
% vector (ix,iy)) is going in the direction we're looking at.  
% 找到各个区间的梯度，其中ix.iy是沿坐标方向的两个梯度，按照上面的划分
% 第一个区间的梯度应满足(iy<=0 & ix<iy) | (iy>0 & ix>iy)，其他区间以此类推
% ix和iy都是m*n的矩阵，find会把他们都当作一个向量来处理（按列存储！）
% find将在他们中找出满足条件的梯度对应的下标，返回的结果是一个列向量，重复的下标只保留一个
switch direction
 case 1
  idx = find((iy<=0 & ix>-iy)  | (iy>=0 & ix<-iy));
 case 2
  idx = find((ix>0 & -iy>=ix)  | (ix<0 & -iy<=ix));
 case 3
  idx = find((ix<=0 & ix>iy) | (ix>=0 & ix<iy));
 case 4
  idx = find((iy<0 & ix<=iy) | (iy>0 & ix>=iy));
end

% Exclude the exterior pixels
% 去除矩阵最外面的一圈
%  1  5  9 13
%  2  6 10 14
%  3  7 11 15
%  4  8 12 16
% v==1和v==0 可以去除第一行和最后一行
% idx<=m可以去除第一列，idx>(n-1)*m可以去除最后一列
if ~isempty(idx)
  v = mod(idx,m);
  extIdx = (v==1 | v==0 | idx<=m | (idx>(n-1)*m));
  idx(extIdx) = [];
end

ixv = ix(idx);  
iyv = iy(idx);   
gradmag = mag(idx);

% Do the linear interpolations for the interior pixels
% 线性插值。对于下面这个矩阵：
%  1  5  9 13
%  2  6 10 14
%  3  7 11 15
%  4  8 12 16
% 比较idx与其梯度方向上相邻两个像素梯度的大小，由于idx梯度方向上不一定存在像素，
% 所以要进行线性插值（注意matlab中，矩阵按列存储）
% 对于case1，假设idx=10，则mag(idx+m)得到的是idx右临像素14的梯度值，
% mag(idx+m-1)得到的是像素13的梯度值。
% d确定了梯度的方向，根据这个方向可以求出插值点与两个参考点的距离，
% 距离越小，这个参考点的权重越大
switch direction
 case 1
  d = abs(iyv./ixv);
  gradmag1 = mag(idx+m).*(1-d) + mag(idx+m-1).*d; 
  gradmag2 = mag(idx-m).*(1-d) + mag(idx-m+1).*d; 
 case 2
  d = abs(ixv./iyv);
  gradmag1 = mag(idx-1).*(1-d) + mag(idx+m-1).*d; 
  gradmag2 = mag(idx+1).*(1-d) + mag(idx-m+1).*d; 
 case 3
  d = abs(ixv./iyv);
  gradmag1 = mag(idx-1).*(1-d) + mag(idx-m-1).*d; 
  gradmag2 = mag(idx+1).*(1-d) + mag(idx+m+1).*d; 
 case 4
  d = abs(iyv./ixv);
  gradmag1 = mag(idx-m).*(1-d) + mag(idx-m-1).*d; 
  gradmag2 = mag(idx+m).*(1-d) + mag(idx+m+1).*d; 
end
idxLocalMax = idx(gradmag>=gradmag1 & gradmag>=gradmag2); 

 

4、滞后门限

在强边缘的八个领域中，如果存在弱边缘，则将弱边缘加入到强边缘中，直到没有新的像素被加入到强边缘。

matlab使用的是bwselect()函数实现这个功能的，该函数可以从某些特定的像素开始，在二值图中查找4连通或者8联通区域。

程序最后，使用了形态学函数bwmorph(e, 'thin', 1)。我们知道，一个高质量的边缘检测算法应该满足两个条件：good localization (edges marked should be as close as possible to the edge in the real image) 以及 good detection(the algorithm should mark as many real edges in the images as possible)。这个函数可以把一条“粗”线，取中间值变成一条“细线”，得到 good localization。非极大值抑制可以得到 good localization，滞后门限可以得到 good detection。

%  1  5  9 13
%  2  6 10 14
%  3  7 11 15
%  4  8 12 16
if ~isempty(idxStrong) % result is all zeros if idxStrong is empty
rstrong = rem(idxStrong-1, m)+1; % 余数，代表强边缘的行号
cstrong = floor((idxStrong-1)/m)+1; % 商，代表强边缘的列号
% 从e中选出某些像素开始的8连通像素
e = bwselect(e, cstrong, rstrong, 8);
%      'thin'         With N = Inf, remove pixels so that an object
%                       without holes shrinks to a minimally
%                       connected stroke, and an object with holes
%                       shrinks to a ring halfway between the hold
%                       and outer boundary
e = bwmorph(e, 'thin', 1);  % Thin double (or triple) pixel wide contours
end