hough_me.m

function [h, theta, rho] = hough_me(f,dtheta, drho)
% HOUGH_ME computes the Hough transform of the image. The parameters that
% takes as input are: 
% dtheta    space of the Hough transform along the theta axis.
% drho      spcae of the Hough transform along the rho axis.

f = double(f);
[M,N] = size(f);
theta = linspace(-90,0,ceil(90/dtheta) + 1);
theta = [theta -fliplr(theta(2:end-1))];
ntheta = length(theta);
% calculate the eucledian distance and round it
D = sqrt((M-1)^2 + (N-1)^2);
q = ceil(D/drho);
nrho = 2*q-1;
rho = linspace(-q*drho, q*drho, nrho);
[x,y,val]=find(f);
x = x-1;
y = y-1;

%initialize output
h = zeros(nrho, length(theta));

for k = 1:ceil(length(val)/1000)
    first = (k-1)*1000+1;
    last = min(first+999, length(x));
    % for each axis repeat and copies the array to create the xmatrix,
    % y_matrix, theta_matrix and rho_matrix.
    x_matrix = repmat(x(first:last), 1, ntheta);
    y_matrix = repmat(y(first:last), 1, ntheta);
    val_matrix = repmat(val(first:last),1,ntheta);
    theta_matrix = repmat(theta, size(x_matrix, 1), 1)*pi/180;
    rho_matrix = x_matrix.*cos(theta_matrix) + ... 
        y_matrix.*sin(theta_matrix);
    
slope = (nrho-1)/(rho(end) - rho(1));
rho_bin_index = round(slope*(rho_matrix - rho(1))+1);
theta_bin_index = repmat(1:ntheta, size(x_matrix,1),1);

h = h + full(sparse(rho_bin_index(:), theta_bin_index(:), ...
    val_matrix(:), nrho, ntheta));
end