function [A, B, R, Votes] = hough (I, min_r, max_r, num_circles)

% Hough transform to detect circles in gray image I.
% Finds top candidates given by num_circles.
% A, B are vectors listing centers.
% R lists the respective radii.
% The vectors will have length num_circles.
% min_r and max_r are min and max radius to expect.
% Circle is given in the form
%    (y-a)^2 + (x-b)^2 = r^2
% For simplicity, we will assume a,b,r are integers.


I = edge(I,'sobel');
figure, imshow(I)

%I = double(I);
[m,n] = size(I);

num_r = max_r - min_r + 1;   % Number of radii to consider.
max_a = m;    % Max possible value of a.  Assume a starts at 1.
max_b = n;    % Max possible value of b.  Assume b starts at 1.

Votes = zeros(max_a, max_b, num_r);


for a = 1:max_a
  for b = 1:max_b
    for r = 1:num_r
      % We know both x and y will be within radius of (a,b).
      % This allows us to cut down our search space for x & y.
      % But we need to we stay within the matrix bounds [1:m,1:n].
      radius = r + min_r - 1;
      min_x = max(1,b-radius);
      max_x = min(n,b+radius);
      for x = min_x:max_x
         min_y = max(1,a-radius);
         max_y = min(m,a+radius);
         for y = min_y:max_y
	 
	    % Determine if a circle given by a,b,r passes through the pixel:
	    %   [ x-0.5 : x+0.5 , y-0.5 : y+0.5 ]
	    % We check each side of the pixel.
	    % First check the pixel is bright enough.
	    if I(y,x) > 0
	    
	      xL = x - 0.5;  % Left side of pixel square.
	      xR = x + 0.5;  % Right side.
	      yT = y + 0.5;  % Top.
	      yB = y - 0.5;  % Bottom.
	      
	      % We generate the corresponding x & y values if the circle
	      % passed through any of the 4 sides of the square.  Note this
	      % gives us 8 values since there is a +/- square root.
	      
	      y1 = a + sqrt(radius^2 - (xL - b)^2);
	      y2 = a - sqrt(radius^2 - (xL - b)^2);
	      y3 = a + sqrt(radius^2 - (xR - b)^2);
	      y4 = a - sqrt(radius^2 - (xR - b)^2);
	      
	      x1 = b + sqrt(radius^2 - (yT - a)^2);
	      x2 = b - sqrt(radius^2 - (yT - a)^2);
	      x3 = b + sqrt(radius^2 - (yB - a)^2);
	      x4 = b - sqrt(radius^2 - (yB - a)^2);	      
	      
	      % Check if any of these 8 values falls within the square.
	      if y1 < yT & y1 > yB
	        Votes(a,b,r) = Votes(a,b,r) + 1;
	      elseif y2 < yT & y2 > yB
	        Votes(a,b,r) = Votes(a,b,r) + 1;	
	      elseif y3 < yT & y3 > yB
	        Votes(a,b,r) = Votes(a,b,r) + 1;
	      elseif y4 < yT & y4 > yB
	        Votes(a,b,r) = Votes(a,b,r) + 1;
	      elseif x1 < xR & x1 > xL
	        Votes(a,b,r) = Votes(a,b,r) + 1;	
	      elseif x2 < xR & x2 > xL
	        Votes(a,b,r) = Votes(a,b,r) + 1;
	      elseif x3 < xR & x3 > xL
	        Votes(a,b,r) = Votes(a,b,r) + 1;
	      elseif x4 < xR & x4 > xL
	        Votes(a,b,r) = Votes(a,b,r) + 1;
	      end;															      
	    end;
	    
    %        if I(y,x) >=128 & (y-a)^2 + (x-b)^2 == radius^2
    %            Votes(a,b,r) = Votes(a,b,r) + 1;
    %       end;
         end;  % end y for loop
      end;  % x for loop
    end;   % r for loop
  end;   % b for loop
end;   % a for loop


% Locate maximum in Votes matrix.
A = zeros(1,num_circles);
B = zeros(1,num_circles);
R = zeros(1,num_circles);
Most_votes = zeros(1,num_circles);  % Contains counts of cells with most votes.

for a = 1:max_a
  for b = 1:max_b
    for r = 1:num_r
          for i = 1:num_circles
            if Votes(a,b,r) > Most_votes(i)
                  % We've found a new circle.
                  % So bump all circles with less votes down the list.
                  for j = num_circles:-1:i+1
                    Most_votes(j) = Most_votes(j-1);
                    A(j) = A(j-1);
                    B(j) = B(j-1);
                    R(j) = R(j-1);
                  end;
		  
                  % Write our new (a,b,r) candidate in.
                  Most_votes(i) = Votes(a,b,r);
                  A(i) = a;
                  B(i) = b;
                  R(i) = r + min_r - 1;
                  break;   % Breaks out of i for loop
                end;
          end;
    end;
  end;
end;
