function [a, b, c, rho, theta, N_lines] = get_line_hough(x, y, rho_min, rho_max, N_rho, theta_min, theta_max, N_theta, N_lines, graphic)

% [a, b, c, rho, theta, N_lines] = get_line_hough(x, y, rho_min, rho_max, N_rho, theta_min, theta_max, N_theta, N_lines, graphic)
%
% DESC:
% calculate the N most voted lines using the Hough transform. The lines are
% parameterized as: cos(theta)*x + sin(theta)*y = rho
%
% AUTHOR
% Marco Zuliani - zuliani@ece.ucsb.edu
%
% VERSION:
% 1.1
%
% INPUT:
% x, y                  = point coordinates
% rho_min, rho_max      = rho interval
% N_rho                 = number of bins for rho
% theta_min, theta_max  = theta interval
% N_theta               = number of bins fro theta
% N_lines               = number of desired lines
% graphic               = true to display the Hough map
% 
%
% OUTPUT:
% a, b, c               = line parametrization: a*x + b*y + c = 0
% rho, theta            = hough parametrization
%
%                          cos(theta)*x + sin(theta)*y = rho
%
% N_lines               = number of lines that were actually found

if nargin < 9
    N_lines = 1;
end;

if nargin < 10
    graphic = false;
end;

N = length(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Build the accumulator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% create the accumulator
A = uint16(zeros(N_theta, N_rho));

Delta_theta = (theta_max - theta_min) / (N_theta - 1);
Delta_rho = (rho_max - rho_min) / (N_rho - 1);

theta = theta_min;
for h = 1:N_theta

    for n = 1:N

        % get the index for rho
        rho = cos(theta)*x(n) + sin(theta)*y(n);
        k = floor( (rho - rho_min) / (rho_max - rho_min) * (N_rho-1) ) + 1;

        % update the accumulator
        if (k > 0) && (k <= N_rho)
            A(h, k) = A(h, k) + 1;
        end;
        %         else
        %             warning('values out of accumulator bounds')
        %         end;

    end;

    % associate angle to bin
    theta = theta + Delta_theta;

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the lines
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% pick the largest N maxima of the accumulator
maxima_locations = imregionalmax(A);
[H_max, K_max] = find(maxima_locations);
ind = sub2ind(size(A), H_max, K_max);
A_max = A(ind);
[A_max_sorted ind_max_sorted] = sort(A_max, 'descend');
N_lines = min(length(A_max_sorted), N_lines);

for h = 1:N_lines
    
    H = H_max(ind_max_sorted(h));
    K = K_max(ind_max_sorted(h));

    % and convert it to the line parametrization
    theta(h) = theta_min + Delta_theta*(H-1);
    rho(h) = rho_min + Delta_rho*(K-1);

    a(h) = cos(theta(h));
    b(h) = sin(theta(h));
    c(h) = -rho(h);

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display the results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (graphic)

    figure
    imagesc([rho_min, rho_max], [theta_min/pi*180, theta_max/pi*180], A);
    hold on
    colormap gray
    plot(rho, theta/pi*180, '+r')
    truesize
    xlabel('\rho')
    ylabel('\theta [degrees]')
    colorbar;

end;

return