Appendix B

Matlab code used to do the Hough Transform on an image containing a straight line.

% Hough Transform MatLab code.
% Version 2.0: Looks for feature points in an image which
% contains a line and replots the line.

% Written by Vicky Safouris and Steven Manos.
% Date last updated: Friday 26/5/2000.
% filename: /home/smanos/hought/htline.m

% Define gradient (orig_m) and y-intercept (orig_b) of a line.

orig_m=0.75
orig_b=20.0

% Set up array (x,y) which consists of all zeros and add a 
% straight line to it according to the input m and b.

im=zeros(100,100);

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

colormap(gray)

% Plot image of array (x,y)

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

% Set up of the accumulator array A(r,theta). r is the longest 
% line possible in the 100 by 100 pixel image, which is given 
% by floor(sqrt(100^2+100^2)) = 141. 360 gives the range of 
% theta we want to examine, being 1 to 360.
 
A=zeros(141,360);

% Here all the values of x and y are searched through to find
% any feature points. For each (x,y) feature point all the
% corresponding (r,theta) points are found. r negative is not
% 'voted for' in the accumulator cell.

for x=1:100;
 for y=1:100;
  if im(y,x)>=200;
   for theta=1:360;
    r=round(x*cos(theta*pi/180) + y*sin(theta*pi/180));
       if r<1
          r=1;
          A(r,theta)=A(r,theta);
       else
          A(r,theta)=A(r,theta)+1;
       end;
     end;
  end;
 end;
end;
A;

% Plot the accumulator array.

subplot(312), image(A*18)
axis equal
axis([1 360 1 141])
axis xy
xlabel('theta')
ylabel('r')
title('Accumulator array')

% The maximum value of the A(r,theta) array is found and the
% corresponding (r,theta) values of this maximum value are
% found. These are defined as maxtheta and maxr.

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

maxtheta
maxr

% The slope m and y-intercept b of the found line is
% calculated.

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

% Defining the array for the output image.

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 the output image.

subplot(313), image(im2)
axis equal
axis([1 100 1 100])
axis xy
xlabel('x')
ylabel('y')
title('Line found')