function [cnt, E] = houghCircles(im, points);

I = im;
%get size of image
[m n] = size(I);

% initialize counter
a = 0;
% 3D matrix (coordinates + radius) 
% to store the result of transform
acc = zeros(m, n, 9);

% search for different radii values
for r = 6 : 2 : 22,
    % update counter
    a = a + 1;
    % circle kenrel  
    ckernel = zeros(2 * r + 1, 2 * r + 1);
  
    for i = 1 : 4 : 360,     
        
        % get circle dimensions
        xo = round(1 + r + r * cos(i * (pi / 180)));
        yo = round(1 + r + r * sin(i * (pi / 180)));
        % update circle kernel
        ckernel(xo, yo) = 1;
        
    end
  
    % find only the whites of I
    [i, j] = find(I);
    np = length(i);
  
    % main hough transform loop - update acc
    for ii = 1 : np,    
        if (i(ii) - r) > 0 & (j(ii) - r) > 0 & (i(ii) + r) < m & (j(ii) + r) < n
            acc(i(ii) - r : i(ii) + r , j(ii) - r : j(ii) + r, a) = ...
            acc(i(ii) - r : i(ii) + r , j(ii) - r : j(ii) + r , a ) + ckernel;
        end
    end
 
end

% get the max values of the 3D matrix and their indices
[vrad irad] = max(acc, [], 3);
acc = vrad;
% initial threshold in an attempt to 
% locate the transoform's local maxima
acc = acc > max(max(acc)) / 2.7;
% prepare our regions of interest 
% for the next step by smoothing
% our blobs -minimize possible interference
x = [-5 : 5];
y = normpdf(x, 0, 2.6);
kernel = y' * y;
acc = conv2(acc, kernel, 'same');

% array to store the center of the circles
cnt = zeros(m,n);
% array to store the final image
E = zeros(m, n);

for t = 1 : 10
    
    % check to see whether all the balls where
    % detected and terminate this iterat - we
    % do that by detecting all the maxima and
    % using a tolerance level of .01 for the max
    if max(max(acc)) < 0.026 | max(max(acc')) < 0.026
        return
    end
    
    % get the indeces of the maximum pt in acc
    [val, yind] = max(max(acc));
    [val, xind] = max(max(acc'));
    
    % In case there are different max values
    % for the rows and cols do it again
    if acc(xind, yind) == 0
        [p, q] = size(acc);
        val = max(max(acc));
        for i = 1 : p
            for j = 1 : q
                if acc(i, j) == val
                    xind = i;
                    yind = j;
                    break;
                end
            end
        end
    end
    
    % erase all the neighboring pixels such that
    % when the loop runs again we will locate
    % the next highest maximum - repeat until
    % all the circles are detected
    acc(xind - 15 : xind + 15, yind - 15 : yind + 15) = 0;
     
    % update cnt matrix
    cnt(xind, yind) = 1;
    % match the counter a to the
    % corresponding radius value
    rr = 2 * irad(xind, yind) + 6;
     
    % loop which computes the circles to
    % be displayed on the final image
    for i = 0 : 2 : 360,
      
        xf  = round(xind + rr * cos(i * (pi / 180)));
        yf = round(yind + rr * sin(i * (pi / 180)));
        E(xf - 1 : xf + 1, yf - 1 : yf + 1) = 1;
        
    end
end
 
return