function [h, theta, rho] = hough(varargin)
%HOUGH Hough transform.
%   HOUGH implements the Standard Hough Transform (SHT). HOUGH is designed
%   to detect lines. It uses the parametric representation of a line:
%
%                       rho = x*cos(theta) + y*sin(theta).
%   
%   The variable rho is the distance from the origin to the line along a 
%   vector perpendicular to the line.  Theta is the angle between
%   the x-axis and this vector. HOUGH generates parameter space matrix 
%   whose rows and columns correspond to rho and theta values respectively.
%   Peak values in this space represent potential lines in the input image.
% 
%   [H, THETA, RHO] = HOUGH(BW) computes the SHT of the binary image BW.
%   THETA (in degrees) and RHO are the arrays of rho and theta values over 
%   which the Hough transform matrix, H, was generated.
%   
%   [H, THETA, RHO] = HOUGH(BW,PARAM1,VAL1,PARAM2,VAL2) sets various
%   parameters.  Parameter names can be abbreviated, and case does not
%   matter. Each string parameter is followed by a value as indicated 
%   below:
%
%   'ThetaResolution' Real scalar between 0 and 90, exclusive.
%                     'ThetaResolution' specifies the spacing (in degrees)
%                     of the Hough transform bins along the theta axis.
%
%                     Default: 1.
%
%   'RhoResolution'   Real scalar between 0 and norm(size(BW)), exclusive.
%                     'RhoResolution' specifies the spacing of the Hough
%                     transform bins along the rho axis.
%
%                     Default: 1.
%
%   Notes
%   -----
%   The Hough transform matrix, H, is NRHO-by-NTHETA where 
%   NRHO = 2*ceil(norm(size(BW))/RhoResolution)-1, and 
%   NTHETA = 2*ceil(90/ThetaResolution). Theta angle values are in 
%   the range [-90, 90) degrees and rho values range from -DIAGONAL to 
%   DIAGONAL where DIAGONAL = RhoResolution*ceil(norm(size(BW))/
%   RhoResolution). Note that if 90/DTHETA is not an integer, 
%   the actual angle spacing will be 90/ceil(90/DTHETA).
%
%   Class Support
%   -------------
%   BW can be logical or numeric and it must be real, 2-D, and nonsparse.
%
%   Example
%   -------
%   Compute and display the Hough transform of the gantrycrane.png image
%
%      RGB = imread('gantrycrane.png');
%      I  = rgb2gray(RGB); % convert to intensity
%      BW = edge(I,'canny'); % extract edges
%      [H,T,R] = hough(BW,'RhoResolution',0.5,'ThetaResolution',0.5);
%
%      % display the original image
%      subplot(2,1,1);
%      imshow(RGB);
%      title('gantrycrane.png');
%
%      % display the hough matrix
%      subplot(2,1,2);
%      imshow(imadjust(mat2gray(H)),'XData',T,'YData',R,...
%             'InitialMagnification','fit');
%      title('Hough transform of gantrycrane.png');
%      xlabel('\theta'), ylabel('\rho');
%      axis on, axis normal, hold on;
%      colormap(hot);
%
%   See also HOUGHPEAKS and HOUGHLINES.

%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.8.1 $  $Date: 2004/08/10 01:39:43 $

%   References:
%   Rafael C. Gonzalez, Richard E. Woods, Steven L. Eddins, "Digital
%   Image Processing Using MATLAB", Prentice Hall, 2004

[bw, rho, theta] = parseInputs(varargin{:});

h = houghmex(bw,rho,theta*pi/180);

%-----------------------------------------------------------------------------

function [bw, rho, theta] = parseInputs(varargin)

iptchecknargin(1,5,nargin,mfilename);

bw     = varargin{1};
iptcheckinput(bw, {'numeric','logical'},...
              {'real', '2d', 'nonsparse', 'nonempty'}, ...
              mfilename, 'BW', 1);

if ~islogical(bw)
  bw = bw~=0;
end

% Set the defaults
theta_res = 1;
rho_res = 1;

% Process parameter-value pairs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
validStrings = {'ThetaResolution','RhoResolution'};

if nargin > 1
  done = false;
  
  idx = 2;
  while ~done
    input = varargin{idx};
    inputStr = iptcheckstrs(input, validStrings,mfilename,'PARAM',idx);
    
    idx = idx+1; %advance index to point to the VAL portion of the input 
    
    if idx > nargin
      eid = sprintf('Images:%s:valFor%sMissing', mfilename, inputStr);
      msg = sprintf('Parameter ''%s'' must be followed by a value.', inputStr);
      error(eid,'%s', msg);        
    end
    
    switch inputStr
      
     case 'ThetaResolution'
      theta_res = varargin{idx};
      iptcheckinput(theta_res, {'double'}, {'real', 'scalar', ...
                          'finite','positive'}, mfilename, inputStr, idx);

       if (theta_res >= 90)
         eid = sprintf('Images:%s:invalidThetaRes', mfilename);
         msg = sprintf(['Value of parameter,''%s'', must be between '...
                        '0 and 90.'], inputStr);
         error(eid,'%s', msg);
       end
      
     case 'RhoResolution'
      rho_res = varargin{idx};
      iptcheckinput(rho_res, {'double'}, {'real', 'scalar', ...
                          'finite','positive'}, mfilename, inputStr, idx);

       if (rho_res >= norm(size(bw)))
         eid = sprintf('Images:%s:invalidThetaRes', mfilename);
         msg = sprintf(['Value of parameter,''%s'', must be between '...
                        '0 and norm(size(BW)).'], inputStr);
         error(eid,'%s', msg);
       end
     
     otherwise
      eid = sprintf('Images:%s:internalError', mfilename);
      msg   = 'Unknown input string.'; %should never get here
      error(eid,'%s', msg);
    end
    
    if idx >= nargin
      done = true;
    end
    
    idx=idx+1;
  end
end

% Compute theta and rho
%%%%%%%%%%%%%%%%%%%%%%%
[M,N] = size(bw);
theta = linspace(-90, 0, ceil(90/theta_res) + 1);
theta = [theta -fliplr(theta(2:end - 1))];

D = sqrt((M - 1)^2 + (N - 1)^2);
q = ceil(D/rho_res);
nrho = 2*q - 1;
rho = linspace(-q*rho_res, q*rho_res, nrho);