How to implement a gradient based Hough transform

I am trying to use the Hough transform for edge detection, and would like to use gradient images as the basis.

What I have done so far, given the image I of size [M,N] and its partial derivatives gx, gy, is to calculate the gradient angle in each pixel as thetas = atan(gy(x,y) ./ gx. Similarly I calculate the gradient magnitude as magnitudes = sqrt(gx.^2+gy.^2).

To build the Hough transform, I use the following MATLAB code:

max_rho = ceil(sqrt(M^2 + N^2));
hough = zeros(2*max_rho, 101);
for x=1:M
for y=1:N
theta = thetas(x,y);
rho = x*cos(theta) + y*sin(theta);

rho_idx = round(rho)+max_rho;
theta_idx = floor((theta + pi/2) / pi * 100) + 1;
hough(rho_idx, theta_idx) = hough(rho_idx, theta_idx) + magnitudes(x,y);
end
end

The resulting Hough transform looks plausible (see http://i.stack.imgur.com/hC9mP.png), but when I try to use its maxima as edge parameters in the original image, the results look more or less random. Did I do something wrong in constructing the Hough transform?

UPDATE: I had a stupid mistake in my code: rho was calculated as x*cos(theta)+y*cos(theta) instead of x*cos(theta)+y*sin(theta). That is, I was using two cosines instead of a cosine and a sine. I have edited the code above and the new resulting image is below. This did not give much better edges though.

@endolith: To plot an edge, given a maximal value in the hough-matrix at rho_idx, theta_idx, I translate the indices to rho,theta values:

theta = (theta_idx -1) / 100 * pi - pi / 2;
rho = rho_idx - max_rho;

Finally I plot the edge as y= (rho - x*cos(theta)) / sin(theta). • "when I try to use its maxima as edge parameters in the original image" How are you doing that? – endolith Oct 27 '11 at 19:41
• @Jonas Due Vesterheden Just wondering is this a time VS doppler frequency image?... – Spacey Nov 6 '11 at 22:25
• @Mohammad: Not that I know of. The original image is of some circuit board. What do you mean by "VS"? – Jonas Due Vesterheden Nov 9 '11 at 14:02
• @JonasDueVesterheden 'VS' just means 'versus'. (Time versus doppler frequency') :-) – Spacey Nov 10 '11 at 16:29
• You should smooth your hough map before applying Non max suppression to it. – user791 Jan 12 '12 at 11:10