Hough transform in MATLAB without using hough function

Question

I found an implementation of the Hough transform in MATLAB at Rosetta Code, but I'm having trouble understanding it. Also I would like to modify it to show the original image and the reconstructed lines (de-Houghing).

Any help in understanding it and de-Houghing is appreciated. Thanks

1. Why is the image flipped?

   theImage = flipud(theImage);

2. I can't wrap my head around the norm function. What is its purpose, and can it be avoided?

EDIT: norm is just a synonym for euclidean distance: sqrt(width^2 + height^2)

rhoLimit = norm([width height]);

1. Can someone provide an explanation of how/why rho, theta, and houghSpace is calculated?

   rho = (-rhoLimit:1:rhoLimit);          
   theta = (0:thetaSampleFrequency:pi);
   
   numThetas = numel(theta);
   houghSpace = zeros(numel(rho),numThetas);
   

2. How would I de-Hough the Hough space to recreate the lines?

Calling the function using a 10x10 image of a diagonal line created using the identity (eye) function

theImage = eye(10)
thetaSampleFrequency = 0.1
[rho,theta,houghSpace] = houghTransform(theImage,thetaSampleFrequency)

The actual function

function [rho,theta,houghSpace] = houghTransform(theImage,thetaSampleFrequency)

    %Define the hough space
    theImage = flipud(theImage);
    [width,height] = size(theImage);

    rhoLimit = norm([width height]);
    rho = (-rhoLimit:1:rhoLimit);          
    theta = (0:thetaSampleFrequency:pi);

    numThetas = numel(theta);
    houghSpace = zeros(numel(rho),numThetas);

    %Find the "edge" pixels
    [xIndicies,yIndicies] = find(theImage);

    %Preallocate space for the accumulator array
    numEdgePixels = numel(xIndicies);
    accumulator = zeros(numEdgePixels,numThetas);

    %Preallocate cosine and sine calculations to increase speed. In
    %addition to precallculating sine and cosine we are also multiplying
    %them by the proper pixel weights such that the rows will be indexed by 
    %the pixel number and the columns will be indexed by the thetas.
    %Example: cosine(3,:) is 2*cosine(0 to pi)
    %         cosine(:,1) is (0 to width of image)*cosine(0)
    cosine = (0:width-1)'*cos(theta); %Matrix Outerproduct  
    sine = (0:height-1)'*sin(theta); %Matrix Outerproduct

    accumulator((1:numEdgePixels),:) = cosine(xIndicies,:) + sine(yIndicies,:);

    %Scan over the thetas and bin the rhos 
    for i = (1:numThetas)
        houghSpace(:,i) = hist(accumulator(:,i),rho);
    end

    pcolor(theta,rho,houghSpace);
    shading flat;
    title('Hough Transform');
    xlabel('Theta (radians)');
    ylabel('Rho (pixels)');
    colormap('gray');

end

Answer 1

The Hough Transform is a "voting" approach where each image point casts a vote on the existence of a certain line (not a line segment) in an image. The voting is carried out in the parameter space for a line: the polar coordinate representation of normal vectors.

We discretize the parameter space and allow each image point to suggest parameters which would be compatible with a line through the point. Each of your questions can be addressed in terms of how the parameter space is treated in code. Wikipedia has a good article with worked examples that might clarify things (if you are having any conceptual troubles).

For your specific questions:

1. The image is flipped so the origin is the bottom right corner. As far as I can tell this step is not technically necessary. It does change the outcome somewhat due to discretization issues. The other implementations on Rosetta Code do not flip the image.
2. rhoLimit holds the maximum radius of an image point in polar coordinates (recall the norm of a vector is its magnitude).
3. rho and theta are discretizations of the polar coordinate plane according to a sampling rate. houghSpace creates a matrix with an element for each possible combination of the discrete rho/theta values.
4. The Hough Transform does not specify the lengths of putative lines; the peaks in the voting space just specify the polar coordinates of the normal vector of the line. You can "de-Hough" by selecting the peaks and drawing the corresponding lines, or perhaps by drawing every possible line and using the number of votes as a grayscale weight. It is not possible to re-create the original image from the Hough Transform, just the lines identified by the transform (and your thresholding scheme on the votes).

Following the example from the question produces the following graph. The placement of grid lines and the datatips cursor can be a bit misleading (though the variable values in the 'tip are correct). Since this is an image of the parameter space and not the image space the sampling rate we chose is determining the number of bins in each variable. At this sampling rate, the image points are compatible with more than one possible line; in other words our lines have subpixel resolution, in the sense that they cannot be drawn without overlap in a 10x10 image.

Once we have chosen a peak, such as that corresponding to the line with normal (rho,theta) = (6.858,0.9), we can draw that line in an image however we choose. Automated peak picking, that is thresholding to find the highly up-voted lines, is its own problem - you could ask a another question about the topic in DSP or about a particular algorithm here.

For example methods see the code and documentation of MATLAB's houghpeaks and houghlines functions.