% 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] = HTEllipse2(x_dir, y_dir, edges, radius_range)

    [rows, columns] = size(edges);
    
    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) = [];
%     alpha_range = pi/8;
        
    indexes_edges = find(edges)';
    s = size(edges);
    
    for alpha = alpha_range
        rot_mat = [cos(alpha) -sin(alpha); sin(alpha) cos(alpha)];
        rot_mat_inv = [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_edges
                    [y, x] = ind2sub(s, i);
                    v = [-x_dir(i); -y_dir(i)];
                    u = rot_mat_inv * v;
                    if u(1) == 0
                        if u(2) > 0
                            theta = pi / 2.;
                        else
                            theta = 3. * pi / 2.;
                        end
                    elseif u(2) == 0
                        if u(1) > 0
                            theta = 0;
                        else
                            theta = pi;
                        end
                    elseif u(2) < 0
                        t = u(2) / u(1);
                        theta = 2*pi + 2 * atan(...
                            (-r1 - r2*t*sqrt((r1^2+r2^2*t^2)/(r2^2*t^2)))...
                            /(r2*t));
                    else
                        t = u(2) / u(1);
                        theta = 2 * atan(...
                            (-r1 + r2*t*sqrt((r1^2+r2^2*t^2)/(r2^2*t^2)))...
                            /(r2*t));
                    end
                    C = round([x; y] - rot_mat * [r1*cos(theta); r2*sin(theta)]);
                    if C(1) < columns && C(1) > 0 && C(2) < rows && C(2) > 0
                        acc(C(2),C(1)) = acc(C(2),C(1)) + 1;
                    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

