GNU Octave: Hough Transform

Question

I am attempting to use a hough transform, unfortunately it doesn't seem to be outputting r and theta that correspond with the lines Drawn. I've been trying to find the answer on this site and others but everything I've tried so far has failed.

I=zeros(80, 80);
for n=1:25;
  I(n+20, n+2)=1;
  I(n+1, n*2+17)=1;
  I(n+1, n*2+16)=1;
  end

hough = houghtf(I,"line", pi*[0:360]/180);
threshHough = hough>.9*max(hough(:));
[r, theta] = find(threshHough>0)
%theta = (theta-91)*pi/180
%r=r-size(hough,1)/2

imshow(I)

Answer 1

The houghtf function in Octave parameterizes a line as

r = x*cos(theta) + y*sin(theta)

The output is an NxM matrix where the first dimension represents r and the second one theta.

The documentation is not very clear on the r parameter. All it says is that N = 2*diag_length - 1, with diag_length the diagonal length of the image. But it doesn't even tell you where the origin is. After some experimentation, I found that the origin is at (N-1)/2. Thus, after finding a peak at indices (ii,jj),

r = ii - (size(hough,1)-1)/2;

The theta parameter is described much better in the documentation: the bin returned corresponds to the values given in the input. So theta is the index into the array pi*[0:360]/180 that you pass to houghft. You can write:

angles = pi*[0:360]/180;
theta = angles(jj);

Do note that the Hough transform deals with the orientation of the line, not the direction of the line. You can parameterize the same line using theta or theta+180. Thus, the hough image produced by the code above is redundant, with all orientations repeated twice (actually, the result for orientation 0 degrees is repeated twice, at 180 and 360 degrees). Instead, use

angles = pi*[0:179]/180;

Next, simply thresholding the Hough transform and getting all pixels is not the best approach: some peaks will be much stronger, leading to bow-tie--shaped detections, of which you'll want to use only the highest point. The 'shrink' method of the bwmorph function is a quick-and-dirty way of reducing the detections after the threshold to a single point, but it's not guaranteed to be on the highest point. Note also that for lines at 0 degrees, the bow-tie is cut in half at the edge of the image, with the other half appearing back at angle 179 degrees. This requires extra effort to fix. Below, I fixed it by extending angles a bit, and then removing duplicate points.

Finally, when parameterizing the line, it seems that Octave picks x as the fist index and y as the second index into the image, which is opposite of what imshow does.

Putting it all together, we get:

% Create test image with some lines
I = zeros(80, 80);
for n=1:25;
  I(n+20,n+2) = 1;
  I(n+1,n*2+17) = 1;
  I(n+1,n*2+16) = 1;
end
I(34:73,55) = 1;
I(60,15:64) = 1;

% Compute the Hough transform
angles = pi*[-10:189]/180;
hough = houghtf(I,"line",angles);

% Detect maxima in the Hough transform -- this is a bit of a hack
detect = hough>.5*max(hough(:));
detect = bwmorph(detect,'shrink',inf);
[ii, jj] = find(detect);
r = ii - (size(hough,1)-1)/2;
theta = angles(jj);

% Remove duplicate points by deleting any point with an angle 
% outside of the interval [0,180).
dup = theta<-1e-6 | theta>=pi-1e-6;
r(dup) = [];
theta(dup) = [];

% Compute line parameters (using Octave's implicit singleton expansion)
r = r(:)';
theta = theta(:)';
x = repmat([1;80],1,length(r)); % 2xN matrix, N==length(r)
y = (r - x.*cos(theta))./sin(theta); % solve line equation for y
% The above goes wrong when theta==0, fix that:
horizontal = theta < 1e-6;
x(:,horizontal) = r(horizontal);
y(:,horizontal) = [1;80];

% Plot
figure
imshow(I)
hold on
plot(y,x,'r-','linewidth',2) % note the reversed x and y!

The diagonal lines are off by one pixel because the way we detected them. We never looked for the location of the local maximum, we took all pixels above a threshold and picked the point in the middle of that.