edge.m

有关matlab的电子书籍有一定的帮助希望有用

function [eout,thresh] = edge(varargin)
%EDGE Find edges in intensity image.
%   EDGE 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.
%
%   EDGE supports six different edge-finding methods:
%
%      The Sobel method finds edges using the Sobel approximation to the
%      derivative. It returns edges at those points where the gradient of
%      I is maximum.
%
%      The Prewitt method finds edges using the Prewitt approximation to
%      the derivative. It returns edges at those points where the gradient
%      of I is maximum.
%
%      The Roberts method finds edges using the Roberts approximation to
%      the derivative. It returns edges at those points where the gradient
%      of I is maximum.
%
%      The Laplacian of Gaussian method finds edges by looking for zero
%      crossings after filtering I with a Laplacian of Gaussian filter.
%
%      The zero-cross method finds edges by looking for zero crossings
%      after filtering I with a filter you specify.
%
%      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.
%
%   The parameters you can supply differ depending on the method you
%   specify. If you do not specify a method, EDGE uses the Sobel method.
%
%   Sobel Method
%   ------------
%   BW = EDGE(I,'sobel') specifies the Sobel method.
%
%   BW = EDGE(I,'sobel',THRESH) specifies the sensitivity threshold for 
%   the Sobel method. EDGE ignores all edges that are not stronger than 
%   THRESH.  If you do not specify THRESH, or if THRESH is empty ([]), 
%   EDGE chooses the value automatically.
%
%   BW = EDGE(I,'sobel',THRESH,DIRECTION) specifies directionality for the
%   Sobel method. DIRECTION is a string specifying whether to look for
%   'horizontal' or 'vertical' edges, or 'both' (the default).
%
%   [BW,thresh] = EDGE(I,'sobel',...) returns the threshold value.
%
%   Prewitt Method
%   --------------
%   BW = EDGE(I,'prewitt') specifies the Prewitt method.
%
%   BW = EDGE(I,'prewitt',THRESH) specifies the sensitivity threshold for
%   the Prewitt method. EDGE ignores all edges that are not stronger than
%   THRESH. If you do not specify THRESH, or if THRESH is empty ([]),
%   EDGE chooses the value automatically.
%
%   BW = EDGE(I,'prewitt',THRESH,DIRECTION) specifies directionality for
%   the Prewitt method. DIRECTION is a string specifying whether to look
%   for 'horizontal' or 'vertical' edges, or 'both' (the default).
%
%   [BW,thresh] = EDGE(I,'prewitt',...) returns the threshold value.
%
%   Roberts Method
%   --------------
%   BW = EDGE(I,'roberts') specifies the Roberts method.
%
%   BW = EDGE(I,'roberts',THRESH) specifies the sensitivity threshold for
%   the Roberts method. EDGE ignores all edges that are not stronger than
%   THRESH. If you do not specify THRESH, or if THRESH is empty ([]),
%   EDGE chooses the value automatically.
%
%   [BW,thresh] = EDGE(I,'roberts',...) returns the threshold value.
%
%   Laplacian of Gaussian Method
%   ----------------------------
%   BW = EDGE(I,'log') specifies the Laplacian of Gaussian method.
%
%   BW = EDGE(I,'log',THRESH) specifies the sensitivity threshold for the
%   Laplacian of Gaussian method. EDGE ignores all edges that are not
%   stronger than THRESH. If you do not specify THRESH, or if THRESH is 
%   empty ([]), EDGE chooses the value automatically.
%
%   BW = EDGE(I,'log',THRESH,SIGMA) specifies the Laplacian of Gaussian
%   method, using SIGMA as the standard deviation of the LoG filter. The
%   default SIGMA is 2; the size of the filter is N-by-N, where
%   N=CEIL(SIGMA*3)*2+1. 
%
%   [BW,thresh] = EDGE(I,'log',...) returns the threshold value.
%
%   Zero-cross Method
%   -----------------
%   BW = EDGE(I,'zerocross',THRESH,H) specifies the zero-cross method,
%   using the specified filter H. If THRESH is empty ([]), EDGE chooses 
%   the sensitivity threshold automatically.
%
%   [BW,THRESH] = EDGE(I,'zerocross',...) returns the threshold value.
%
%   Canny Method
%   ----------------------------
%   BW = EDGE(I,'canny') specifies the Canny method.
%
%   BW = EDGE(I,'canny',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 = EDGE(I,'canny',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] = EDGE(I,'canny',...) returns the threshold values as a
%   two-element vector.
%
%   Class Support
%   -------------
%   I can be of class uint8 or double. BW is of class uint8.
%
%   Remarks
%   -------
%   For the 'log' and 'zerocross' methods, if you specify a
%   threshold of 0, the output image has closed contours, because
%   it includes all of the zero crossings in the input image.
%
%   Example
%   -------
%   Find the edges of the rice.tif image using the Prewitt and Canny
%   methods:
%
%       I = imread('rice.tif');
%       BW1 = edge(I,'prewitt');
%       BW2 = edge(I,'canny');
%       imshow(BW1)
%       figure, imshow(BW2)
%
%   See also FSPECIAL.

%   Grandfathered syntax
%   --------------------
%   BW = EDGE(... ,K) allows the specification of a directionality
%   factor, K.  This only works for the 'sobel', 'prewitt', and
%   'roberts' methods.   K must be a 1-by-2 vector, K = [kx ky].
%   For Sobel and Prewitt, K=[1 0] looks for vertical edges,
%   K=[0 1] looks for horizontal edges, and K=[1 1], the default,
%   looks for non-directional edges.   For the Roberts edge detector,
%   K=[1 0] looks for 135 degree diagonal edges, K=[0 1] looks
%   for 45 degree diagonal edges, and K=[1 1], the default, looks
%   for non-directional edges.
%
%   Clay M. Thompson 10-8-92
%   Revised by Chris Griffin, 1996,1997
%   Copyright 1993-1998 The MathWorks, Inc.  All Rights Reserved.
%   $Revision: 5.12 $  $Date: 1997/11/24 15:34:36 $

[a,method,thresh,sigma,H,kx,ky] = parse_inputs(varargin{:});

% Transform to a double precision intensity image
if isa(a, 'uint8')  
   a = im2double(a);
end

m = size(a,1);
n = size(a,2);
rr = 2:m-1; cc=2:n-1;

% The output edge map:
e = repmat(logical(uint8(0)), m, n);

if strcmp(method,'canny')
   % 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
   
   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: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
   
   % 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
   
elseif any(strcmp(method, {'log','marr-hildreth','zerocross'}))
   % We don't use image blocks here
   if isempty(H),
      fsize = ceil(sigma*3) * 2 + 1;  % choose an odd fsize > 6*sigma;
      op = fspecial('log',fsize,sigma); 
   else 
      op = H; 
   end
   
   op = op - sum(op(:))/prod(size(op)); % make the op to sum to zero
   b = filter2(op,a);
   
   if isempty(thresh)
      thresh = .75*mean2(abs(b(rr,cc)));
   end
  
   % Look for the zero crossings:  +-, -+ and their transposes 
   % We arbitrarily choose the edge to be the negative point
   [rx,cx] = find( b(rr,cc) < 0 & b(rr,cc+1) > 0 ...
      & abs( b(rr,cc)-b(rr,cc+1) ) > thresh );   % [- +]
   e((rx+1) + cx*m) = 1;
   [rx,cx] = find( b(rr,cc-1) > 0 & b(rr,cc) < 0 ...
      & abs( b(rr,cc-1)-b(rr,cc) ) > thresh );   % [+ -]
   e((rx+1) + cx*m) = 1;
   [rx,cx] = find( b(rr,cc) < 0 & b(rr+1,cc) > 0 ...
      & abs( b(rr,cc)-b(rr+1,cc) ) > thresh);   % [- +]'
   e((rx+1) + cx*m) = 1;
   [rx,cx] = find( b(rr-1,cc) > 0 & b(rr,cc) < 0 ...
      & abs( b(rr-1,cc)-b(rr,cc) ) > thresh);   % [+ -]'
   e((rx+1) + cx*m) = 1;
   
   % Most likely this covers all of the cases.   Just check to see if there
   % are any points where the LoG was precisely zero:
   [rz,cz] = find( b(rr,cc)==0 );
   if ~isempty(rz)
      % Look for the zero crossings: +0-, -0+ and their transposes
      % The edge lies on the Zero point
      zero = (rz+1) + cz*m;   % Linear index for zero points
      zz = find(b(zero-1) < 0 & b(zero+1) > 0 ...
         & abs( b(zero-1)-b(zero+1) ) > 2*thresh);     % [- 0 +]'
      e(zero(zz)) = 1;
      zz = find(b(zero-1) > 0 & b(zero+1) < 0 ...
         & abs( b(zero-1)-b(zero+1) ) > 2*thresh);     % [+ 0 -]'
      e(zero(zz)) = 1;
      zz = find(b(zero-m) < 0 & b(zero+m) > 0 ...
         & abs( b(zero-m)-b(zero+m) ) > 2*thresh);     % [- 0 +]
      e(zero(zz)) = 1;
      zz = find(b(zero-m) > 0 & b(zero+m) < 0 ...
         & abs( b(zero-m)-b(zero+m) ) > 2*thresh);     % [+ 0 -]
      e(zero(zz)) = 1;
   end

else  % one of the easy methods (roberts,sobel,prewitt)
   
   % Determine edges in blocks for easy methods 
   nr = length(rr); nc = length(cc);
   
   blk = bestblk([nr nc]);
   nblks = floor([nr nc]./blk); nrem = [nr nc] - nblks.*blk;
   mblocks = nblks(1); nblocks = nblks(2);
   mb = blk(1); nb = blk(2);
   
   if strcmp(method,'sobel')
      op = [-1 -2 -1;0 0 0;1 2 1]/8; % Sobel approximation to derivative