% Hough Transform for circles.
% img : logical image, zeros are the circles edges
% radius_range is a vector : [min, max]
function [acc_votes, acc_radius1, acc_radius2, acc_alpha] = HTEllipse(img, radius_range)

    [rows, columns] = size(img);
    
    acc_votes = zeros(rows, columns);
    acc_radius1 = zeros(rows, columns);
    acc_radius2 = zeros(rows, columns);
    acc_alpha = zeros(rows, columns);
    
    alpha_range = linspace(0, pi / 2, 8);
    alpha_range(end) = [];
    
    phi_range = linspace(0, 2 * pi, 32);
    phi_range(end) = [];
    
    indexes_non_zero = find(img)';
    s = size(img);
    
    for alpha = alpha_range
        transform_mat = [cos(alpha) -sin(alpha); sin(alpha) cos(alpha)];        
        for r1 = radius_range(1):radius_range(2)
            for r2 = radius_range(1):radius_range(2)
                acc = zeros(rows, columns);
                for i = indexes_non_zero
                    [y, x] = ind2sub(s, i);
                    for phi = phi_range
                        p0 = [x; y] - transform_mat * [r1 * cos(phi); r2 * sin(phi)];
                        x0 = round(p0(1));
                        y0 = round(p0(2));

                        if x0 < columns && x0 > 0 && y0 < rows && y0 > 0
                            acc(y0, x0) = acc(y0, x0) + 1;
                        end
                    end
                end

                for x = 1:columns
                    for y = 1:rows
                        if acc(y, x) > acc_votes(y, x)
                            acc_votes(y, x) = acc(y, x);
                            acc_radius1(y, x) = r1;
                            acc_radius2(y, x) = r2;
                            acc_alpha(y, x) = alpha;
                        end
                    end
                end      
            end
        end
    end
    
    acc_votes = imgaussfilt(acc_votes, 1.0);
end

