function [eout] = edgedetector(varargin)
%EDGEDETECTOR Find edges in intensity image.
%
%   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 = EDGEDETECTOR(I,'sobel') specifies the Sobel method.
%
%   BW = EDGE(I,'sobel',DIRECTION) specifies directionality for the
%   Sobel method. DIRECTION is a string specifying whether to look for
%   'horizontal' or 'vertical' edges, or 'both' (the default)
%
%   Prewitt Method
%   --------------
%   BW = EDGEDETECTOR(I,'prewitt') specifies the Prewitt method.
%
%   BW = EDGE(I,'prewitt',DIRECTION) specifies directionality for
%   the Prewitt method. DIRECTION is a string specifying whether to look
%   for 'horizontal' or 'vertical' edges, or 'both' (the default).
%
%   Roberts Method
%   --------------
%   BW = EDGEDETECTOR(I,'roberts') specifies the Roberts method.
%
%   BW = EDGE(I,'roberts',OPTION) specifies the specific formuala for
%   the Roberts method. OPTION is a string specifying whether to use
%   'absolute' or 'square' formula.'absolute' (the default).
%
%   Laplacian
%   ----------------------------
%   BW = EDGEDETECTOR(I,'laplacian') specifies the Laplacian of Gaussian method.
%
%   Kirsch Method
%   -----------------
%   BW = EDGEDETECTOR(I,'kirsch') specifies the kirsch method,
%
%
%   Robinson Method
%   -----------------
%   BW = EDGEDETECTOR(I,'robison') specifies the robinson method,
%
%   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 sobel
%   method:
%
%       I = imread('rice.tif');
%       BW = edgedetector(I,'sobel');
%       figure; colormap(gray(256)); image(BW);
%

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

% Transform to a double precision intensity image
if isa(a, 'uint8') | isa(a, 'uint16') 
   a = double(a);
end
if strcmp(method,'roberts')
   op = [1 0;0 -1];
   bx = abs(filter2(op,a)); by = abs(filter2(rot90(op),a));
   if strcmp(roption,'absolute')
       b = bx + by;
   else
       b = sqrt(kx*bx.*bx + ky*by.*by);
   end
   strtitle = 'Roberts Edge Map';
      
elseif any(strcmp(method, {'sobel','prewitt'}))
   if strcmp(method,'sobel')
      op = [-1 -2 -1;0 0 0;1 2 1];
      strtitle = 'Sobel Edge Map';
   else 
      op = [-1 -1 -1;0 0 0;1 1 1];
      strtitle = 'Prewitt Edge Map';
   end
   bx = abs(filter2(op',a)); by = abs(filter2(op,a));
   b = sqrt(kx*bx.*bx + ky*by.*by);
   
elseif any(strcmp(method, {'kirsch','robinson'}))
   m = size(a,1);
   n = size(a,2);
   if strcmp(method,'kirsch')
       op0 = [-3 -3 5;-3 0 5;-3 -3 5]/8;
       op1 = [-3 5 5; -3 0 5 ; -3 -3 -3]/8;
       strtitle = 'Kirsch Edge Map';
   else
       op0 = [-1 0 1;-2 0 2;-1 0 1]/2;
       op1 = [0 1 2; -1 0 1; -2 -1 0]/2;
       strtitle = 'Robinson Edge Map';
   end
   op2 = rot90(op0);
   op4 = rot90(op2);
   op6 = rot90(op4);
   op3 = rot90(op1);
   op5 = rot90(op3);
   op7 = rot90(op5);
   b0 = abs(filter2(op0,a));
   b1 = abs(filter2(op1,a));
   b2 = abs(filter2(op2,a));
   b3 = abs(filter2(op3,a));
   b4 = abs(filter2(op4,a));
   b5 = abs(filter2(op5,a));
   b6 = abs(filter2(op6,a));
   b7 = abs(filter2(op7,a));
   for i=1:m
     for j=1:n
       temp = 0;
       if b0(i,j) > temp
         temp = b0(i,j);
       end
       if b1(i,j) > temp
         temp = b1(i,j);
       end
       if b2(i,j) > temp
         temp = b2(i,j);
       end
       if b3(i,j) > temp
         temp = b3(i,j);
       end
       if b4(i,j) > temp
         temp = b4(i,j);
       end
       if b5(i,j) > temp
         temp = b5(i,j);
       end
       if b6(i,j) > temp
         temp = b6(i,j);
       end
       if b7(i,j) > temp
         temp = b7(i,j);
       end
       b(i,j) = temp;
     end
   end
  
elseif strcmp(method,'laplacian')
   op = [-1 -1 -1; -1 8 -1; -1 -1 -1];
   b=filter2(op,a);
   strtitle = 'Laplacian Edge Map';
else
   error([method,' is not a valid method.']);
end

if nargout==0,
   % Plot the original image
    subplot(1,2,1);
    imagesc(a);
    colormap(gray);
    axis square;
    axis off;
    grid off;
    title('Original Image','FontSize',12,'Color','b')

    % Plot the reversed image
    subplot(1,2,2);
    imagesc(b);
    colormap(gray);
    axis square;
    axis off;
    grid off;
    title(strtitle,'FontSize',12,'Color','b')
%    figure; colormap(gray(256)); image(a); title('Original Image');
%    figure; colormap(gray(256)); image(b); title(strtitle);
else
   eout = b;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Local Function : parse_inputs
%
function [I,Method,kx,ky,Roption] = parse_inputs(varargin)
% OUTPUTS:
%   I      Image Data
%   Method Edge detection method
%   kx,ky  From Directionality vector
%   Roption  For Roberts method(absolute or square formula)

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

I = varargin{1};

% Defaults
Method='roberts';
Direction='both';
Roption='absolute';

methods = {'roberts','prewitt','sobel','kirsch','robinson','laplacian'};
directions = {'both','horizontal','vertical'};
option = {'absolute','square'};

% Now parse the nargin-1 remaining input arguments

for i = 2:nargin
   if ischar(varargin{i})
      str = lower(varargin{i});
      j = strmatch(str,methods);
      k = strmatch(str,directions);
      l = strmatch(str,option);
      if ~isempty(j)
         Method = methods{j(1)};
      elseif ~isempty(k)
         Direction = directions{k(1)};
      elseif ~isempty(l)
         Roption = option{l(1)};
      else
         error(['Invalid input string: ''' varargin{i} '''.']);
      end
   end
end

% Now get the direction arguments

switch Direction
case 'both',
   kx = 1; ky = 1; 
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