How to implement a gradient based Hough transform

Question

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 https://i.sstatic.net/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).

Answer 1

I have no idea whether this an issue but atan() only gives you angles from -90 to +90 degree because of the quadrant ambiguity. To get the full gradient angle (from -180 to 180) you need to use atan2().

Answer 2

I am a bit confused by your question. Hough transform is used to detect lines, not edges.

If all you want is an edge map, you should simply threshold the gradient magnitude, or use something more fancy like the Canny edge detector.

If you want to detect straight lines, you would be better off starting with an edge map, and then using the hough function if the Image Processing toolbox, if you have access to it. The problem with doing a Hough transform on the gradient is that edge pixels forming a straight line might have opposite gradient orientations. For example, consider a checkerboard pattern. A edge between two rows of squares flips orientation depending on whether you have a black square above and a white square below, or the other way around.

As to your implementation, I think the problem is that the bins in your Hough matrix are too small. Essentially the bin size in the rho dimension is 1, and the bin size in the theta dimension is less than 2 degrees. That means that the gradient orientations must alight very precisely to form a line, which rarely happens in practice. If you caculate rho_idx and theta_idx so that the bins are larger, that will make your line detector more tolerant to errors, and you might get better lines.