Appendix C

Matlab code used to implement the Hough Transform on an image containing two straight lines embedded in uniformly distributed noise.

% Hough Transform Matlab code.
% Version 3.0: Looks for lines in an image containing random
% noise.

% Written by Vicky Safouris and Steven Manos.
% Date last updated: Tuesday 30/5/2000.
% filename: /home/smanos/hought/htrand.m

colordef white

% Define gradient and y-intercept of first line in (x,y)-image.

orig_m=2.0
orig_b=10.0

% Create 100 by 100 image of 'noise'. This noise is made from
% uniformly distributed random numbers in the range [0,1].

im=rand(100,100)*65;

% Place a straight line in noisy image with slope m and y-intercept 
% b, where the line is 'added' to the noise.

for x=1:100;
  y=floor(orig_m*x + orig_b);
  if(y <= 100 & y >=1)
    im(y,x)=im(y,x)+110;
  end;
end;

% Place a second line in the noisy image.

orig_m1=-0.7
orig_b1=85.0

for x=1:100;
  y=floor(orig_m1*x + orig_b1);
  if(y <= 100 & y >=1)
    im(y,x)=im(y,x)+110;
  end;
end;

colormap(gray)

% Plot of input image.

subplot(311), image(im)
title('Input image')
axis equal
axis([1 100 1 100])
axis xy
xlabel('x')
ylabel('y')
title('Input image')

% set up of array to precalculate sin and cos values of theta
% to optimise code.

for thet=1:1:360;
        sinarray(thet)=sin(thet*pi/180.0);
        cosarray(thet)=cos(thet*pi/180.0);
end;

% Define accumulator array with maximum r of 141.

A=zeros(141,360);

% Loop through all of the (x,y) points, not just the feature
% points (as in versions 2 & 3).

for x=1:100;
 for y=1:100;
   for theta=1:360;
    r=round(x*cosarray(theta) + y*sinarray(theta));
         if r<1
              r=1;
              A(r,theta)=A(r,theta);
         else
              A(r,theta)=A(r,theta)+im(y,x);
     end;
  end;
 end;
end;

% Search for the maximum element in A(r,theta) and its (r,theta)
% coordinates.

m=max(max(A));
for r=1:141;
 for theta=1:360;
  if A(r,theta)==m
      maxr = r;
  maxtheta = theta;
  end;
 end;
end;

% calculate slope m and y-intercept b of the line found.

m=-1.0*((cos(maxtheta*pi/180.0))/(sin(maxtheta*pi/180.0)))
b=maxr/(sin(maxtheta*pi/180.0))

% Hough space (accumulator array) plot of lines obtained from 
% r and theta values.

subplot(312), image(A*0.008)
title('Hough transform space plot (r vs theta)')
axis equal
axis([1 360 1 141])
axis xy
ylabel('r')
xlabel('theta');
title('Accumulator array')

% Setting up and plotting the output image with the found line.

im2=zeros(100,100);
for x=1:100;
  y=round(m*x + b);
  if(y <= 100 & y >=1)
    im2(y,x)=255;
  end;
end;

% Plot of the output image.

subplot(313)
image(im2);
title('Plot of lines found in the input image')
axis equal
axis([1 100 1 100])
axis xy
ylabel('y');
xlabel('x');
title('Line found')

% 3D plot of the accumulator array.

figure;
colormap(bone);
subplot(212), surfl(A);
shading interp;
axis([1 360 1 141 1 10000])
ylabel('r');
xlabel('theta');
zlabel('value');
title('3D plot of the accumulator array')