图像处理 之 Canny 算子

function [eout,thresh] = canny(varargin)
%canny:Find edges in intensity image.
%canny takes an intensity image I as its input, and returns a binary image
% BW of the same size as I, with 1's where the function finds edges in I
% and 0's elsewhere.
%
%
%
% The Canny method finds edges by looking for local maxima of the
% gradient of I. The gradient is calculated using the derivative of a
% Gaussian filter. The method uses two thresholds, to detect strong
% and weak edges, and includes the weak edges in the output only if
% they are connected to strong edges. This method is therefore less
% likely than the others to be "fooled" by noise, and more likely to
% detect true weak edges.
%
%
% Canny Method
% ----------------------------
% BW = canny(I) specifies the Canny method.
%
% BW = canny(I,THRESH) specifies sensitivity thresholds for the
% Canny method. THRESH is a two-element vector in which the first element
% is the low threshold, and the second element is the high threshold. If
% you specify a scalar for THRESH, this value is used for the high
% threshold and 0.4*THRESH is used for the low threshold. If you do not
% specify THRESH, or if THRESH is empty ([]), EDGE chooses low and high
% values automatically.
%
% BW = canny(I,THRESH,SIGMA) specifies the Canny method, using
% SIGMA as the standard deviation of the Gaussian filter. The default
% SIGMA is 1; the size of the filter is chosen automatically, based
% on SIGMA.
%
% [BW,thresh] = canny(I,...) returns the threshold values as a
% two-element vector.
%
% Class Support
% -------------
% I can be of class uint8, uint16, or double. BW is of class uint8.
%
% Example
% -------
% Find the edges of the rice.tif image using the Prewitt and Canny
% methods:
%
% I = imread('rice.tif');
% BW2 = canny(I);
%
% Clay M. Thompson 10-8-92
% Revised by Chris Griffin, 1996,1997
% Copyright 1993-1998 The MathWorks, Inc. All Rights Reserved.
% $Revision: 5.14 $ $Date: 1998/12/17 17:20:14 $

// 取得输入参数 --> 实际上，为参数分配空间
[a,thresh,sigma,H,kx,ky] = parse_inputs(varargin{:});

% %变成双精度强度图像

if isa(a, 'uint8') | isa(a, 'uint16')
a = im2double(a);
end

m = size(a,1); 宽
n = size(a,2); 高

rr = 2:m-1; cc=2:n-1; //没用到这二个

% 建一个 m*n 的矩阵，并置0
e = repmat(logical(uint8(0)), m, n);

% 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

// 设计高斯滤波器--具体方式自己查下吧！ 参数放在了 gau,dgau上面
pw = 1:30; % possible widths
ssq = sigma*sigma;
width = max(find(exp(-(pw.*pw)/(2*sigma*sigma))>GaussianDieOff));
if isempty(width)
width = 1; % the user entered a really small sigma
end
t = (-width);
len = 2*width+1;
t3 = [t-.5; t; t+.5]; % We will average values at t-.5, t, t+.5
gau = sum(exp(-(t3.*t3)/(2*ssq))).'/(6*pi*ssq); % the gaussian 1-d filter
dgau = (-t.* exp(-(t.*t)/(2*ssq))/ ssq).'; % derivative of a gaussian

% Convolve the filters with the image in each direction
% The canny edge detector first requires convolutions with
% the gaussian, and then with the derivitave of a gauusian.
% I convolve the filters first and then make a call to conv2
% to do the convolution down each column.

ra = size(a,1);
ca = size(a,2);
ay = 255*a;
ax = 255*a';

h = conv(gau,dgau); %得取了传递函数
ax = conv2(ax, h, 'same').';
ay = conv2(ay, h, 'same');
mag = sqrt((ax.*ax) + (ay.*ay));
magmax = max(mag(:));
if magmax>0
mag = mag / magmax; % normalize
end
%mag 中存的是 通过高斯核卷积, 后的值

//////////////////////////////////////
// 下面是高，低阀值的计算: 结果放到 thresh 中去

% Select the thresholds 选择阈值
if isempty(thresh)
[counts,x]=imhist(mag, 64);
highThresh = min(find(cumsum(counts) > PercentOfPixelsNotEdges*m*n)) / 64;
lowThresh = ThresholdRatio*highThresh;
thresh = [lowThresh highThresh];
elseif length(thresh)==1
highThresh = thresh;
if thresh>=1
error('The threshold must be less than 1.');
end
lowThresh = ThresholdRatio*thresh;
thresh = [lowThresh highThresh];
elseif length(thresh)==2
lowThresh = thresh(1);
highThresh = thresh(2);
if (lowThresh >= highThresh) | (highThresh >= 1)
error('Thresh must be [low high], where low < high < 1.');
end
end

// 下面进行，计算弱边缘

% The next step is to do the non-maximum supression.
% 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)];
end

rstrong = rem(idxStrong-1, m)+1;
cstrong = floor((idxStrong-1)/m)+1;
e = bwselect(e, cstrong, rstrong, 8);
e = bwmorph(e, 'thin', 1); % Thin double (or triple) pixel wide contours

if nargout==0,
imshow(e);
else
eout = e;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Local Function : cannyFindLocalMaxima
% 得到梯度各个方向上的一阶导数的局部极大值点
function idxLocalMax = cannyFindLocalMaxima(direction,ix,iy,mag);
%
% This sub-function helps with the non-maximum supression 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 supression
% (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,o] = 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.

// 查出四种情况的点，放到 idx 中去
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
if ~isempty(idx)
v = mod(idx,m);
extIdx = find(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
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);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Local Function : parse_inputs
%
function [I,Thresh,Sigma,H,kx,ky] = parse_inputs(varargin)
% OUTPUTS:
% I Image Data
% Thresh Threshold value 阀值
% Sigma standard deviation of Gaussian //用来计算高斯滤波器的
% H Filter for Zero-crossing detection 滤波器对二阶导数过零检测
% kx,ky From Directionality vector // 方向矢量..估计是外国人想的周全，把你会不会计算单方向都考虑了！

error(nargchk(1,4,nargin));

I = varargin{1};

% Defaults
Thresh=[];
Direction='both';
Sigma=2;
H=[];
K=[1 1];

directions = {'both','horizontal','vertical'};

% Now parse the nargin-1 remaining input arguments
% Get the rest of the arguments

Sigma = 1.0; % Default Std dev of gaussian for canny
threshSpecified = 0; % Threshold is not yet specified
nonstr=[];
for i = nonstr
if prod(size(varargin{i}))==2 & ~threshSpecified
Thresh = varargin{i};
threshSpecified = 1;
elseif prod(size(varargin{i}))==1
if ~threshSpecified
Thresh = varargin{i};
threshSpecified = 1;
else
Sigma = varargin{i};
end
elseif isempty(varargin{i}) & ~threshSpecified
% Thresh = [];
threshSpecified = 1;
else
error('Invalid input arguments');
end
end

if Sigma<=0
error('Sigma must be positive');
end

switch Direction
case 'both',
kx = K(1); ky = K(2);
case 'horizontal',
kx = 0; ky = 1; % Directionality factor
case 'vertical',
kx = 1; ky = 0; % Directionality factor
otherwise
error('Unrecognized direction string');
end

if isrgb(I)
error('RGB images are not supported. Call RGB2GRAY first.');
end