edge.m

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

      bx = abs(filter2(op',a)); by = abs(filter2(op,a));
      b = kx*bx.*bx + ky*by.*by;
      if isempty(thresh), % Determine cutoff based on RMS estimate of noise
         cutoff = 4*sum(sum(b(rr,cc)))/prod(size(b(rr,cc))); thresh = sqrt(cutoff);
      else                   % Use relative tolerance specified by the user
         cutoff = (thresh).^2;
      end
      rows = 1:blk(1);
      for i=0:mblocks,
         if i==mblocks, rows = (1:nrem(1)); end
         for j=0:nblocks,
            if j==0, cols = 1:blk(2); elseif j==nblocks, cols=(1:nrem(2)); end
            if ~isempty(rows) & ~isempty(cols)
               r = rr(i*mb+rows); c = cc(j*nb+cols);
               e(r,c) = (b(r,c)>cutoff) & ...
               ( ( (bx(r,c) >= (kx*by(r,c)-eps*100)) & ...
               (b(r,c-1) <= b(r,c)) & (b(r,c) > b(r,c+1)) ) | ...
               ( (by(r,c) >= (ky*bx(r,c)-eps*100 )) & ...
               (b(r-1,c) <= b(r,c)) & (b(r,c) > b(r+1,c))));
            end
         end
      end
      
   elseif strcmp(method,'prewitt')
      op = [-1 -1 -1;0 0 0;1 1 1]/6; % Prewitt approximation to derivative
      bx = abs(filter2(op',a)); by = abs(filter2(op,a));
      b = kx*bx.*bx + ky*by.*by;
      if isempty(thresh), % Determine cutoff based on RMS estimate of noise
         cutoff = 4*sum(sum(b(rr,cc)))/prod(size(b(rr,cc))); thresh = sqrt(cutoff);
      else                   % Use relative tolerance specified by the user
         cutoff = (thresh).^2;
      end
      rows = 1:blk(1);
      for i=0:mblocks,
         if i==mblocks, rows = (1:nrem(1)); end
         for j=0:nblocks,
            if j==0, cols = 1:blk(2); elseif j==nblocks, cols=(1:nrem(2)); end
            if ~isempty(rows) & ~isempty(cols)
               r = rr(i*mb+rows); c = cc(j*nb+cols);
               e(r,c) = (b(r,c)>cutoff) & ...
               ( ( (bx(r,c) >= (kx*by(r,c)-eps*100) ) & ...
               (b(r,c-1) <= b(r,c)) & (b(r,c) > b(r,c+1)) ) | ...
               ((by(r,c) >= (ky*bx(r,c)-eps*100) )  & ...
               (b(r-1,c) <= b(r,c)) & (b(r,c) > b(r+1,c)) ) );
            end
         end
      end
      
   elseif strcmp(method, 'roberts')
      op = [1 0;0 -1]/sqrt(2); % Roberts approximation to diagonal derivative
      bx = abs(filter2(op,a)); by = abs(filter2(rot90(op),a));
      b = kx*bx.*bx + ky*by.*by;
      if isempty(thresh), % Determine cutoff based on RMS estimate of noise
         cutoff = 6*sum(sum(b(rr,cc)))/prod(size(b(rr,cc))); thresh = sqrt(cutoff);
      else                   % Use relative tolerance specified by the user
         cutoff = (thresh).^2;
      end
      rows = 1:blk(1);
      for i=0:mblocks,
         if i==mblocks, rows = (1:nrem(1)); end
         for j=0:nblocks,
            if j==0, cols = 1:blk(2); elseif j==nblocks, cols=(1:nrem(2)); end
            if ~isempty(rows) & ~isempty(cols)
               r = rr(i*mb+rows); c = cc(j*nb+cols);
               e(r,c) = (b(r,c)>cutoff) & ...
               ( ( (bx(r,c) >= (kx*by(r,c)-eps*100)) & ...
               (b(r-1,c-1) <= b(r,c)) & (b(r,c) > b(r+1,c+1)) ) | ...
               ( (by(r,c) >= (ky*bx(r,c)-eps*100)) & ...
               (b(r-1,c+1) <= b(r,c)) & (b(r,c) > b(r+1,c-1)) ) );
            end
         end
      end
   else
      error([method,' is not a valid method.']);
   end
end

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.  

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,Method,Thresh,Sigma,H,kx,ky] = parse_inputs(varargin)
% OUTPUTS:
%   I      Image Data
%   Method Edge detection method
%   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
Method='sobel';
Thresh=[];
Direction='both';
Sigma=2;
H=[];
K=[1 1];

methods = {'canny','prewitt','sobel','marr-hildreth','log','roberts','zerocross'};
directions = {'both','horizontal','vertical'};

% Now parse the nargin-1 remaining input arguments

% First get the strings - we do this because the intepretation of the 
% rest of the arguments will depend on the method.
nonstr = [];   % ordered indices of non-string arguments
for i = 2:nargin
   if ischar(varargin{i})
      str = lower(varargin{i});
      j = strmatch(str,methods);
      k = strmatch(str,directions);
      if ~isempty(j)
         Method = methods{j(1)};
         if strcmp(Method,'marr-hildreth')  
            warning('''Marr-Hildreth'' has been grandfathered, use ''LoG'' instead.');
         end
      elseif ~isempty(k)
         Direction = directions{k(1)};
      else
         error(['Invalid input string: ''' varargin{i} '''.']);
      end
   else
      nonstr = [nonstr i];
   end
end

% Now get the rest of the arguments 

switch Method
   
case {'prewitt','sobel','roberts'}
   threshSpecified = 0;  % Threshold is not yet specified
   for i = nonstr
      if prod(size(varargin{i}))<=1 & ~threshSpecified % Scalar or empty
         Thresh = varargin{i};
         threshSpecified = 1;
      elseif prod(size(varargin{i}))==2  % The dreaded K vector
         K=varargin{i};     
      else
         error('Invalid input arguments');
      end
   end
   
case 'canny'
   Sigma = 1.0;          % Default Std dev of gaussian for canny
   threshSpecified = 0;  % Threshold is not yet specified
   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
      
case 'log'
   threshSpecified = 0;  % Threshold is not yet specified
   for i = nonstr
      if prod(size(varargin{i}))<=1  % Scalar or empty
         if ~threshSpecified
            Thresh = varargin{i};
            threshSpecified = 1;
         else
            Sigma = varargin{i};
         end
      else
         error('Invalid input arguments');
      end
   end
   
case 'zerocross'
   threshSpecified = 0;  % Threshold is not yet specified
   for i = nonstr
      if prod(size(varargin{i}))<=1 & ~threshSpecified % Scalar or empty
         Thresh = varargin{i};
         threshSpecified = 1;
      elseif prod(size(varargin{i})) > 1 % The filter for zerocross
         H = varargin{i};
      else
         error('Invalid input arguments');
      end
   end

case 'marr-hildreth'
   for i = nonstr
      if prod(size(varargin{i}))<=1  % Scalar or empty
         Thresh = varargin{i};
      elseif prod(size(varargin{i}))==2  % The dreaded K vector 
         warning('The [kx ky] direction factor has no effect for ''Marr-Hildreth''.');
      elseif prod(size(varargin{i})) > 2 % The filter for zerocross
         H = varargin{i};
      else
         error('Invalid input arguments');
      end
   end
   
otherwise
   error('Invalid input arguments');
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