% Shahram Beheshti     10/20/02



im = imread('Shahram3.jpg'); % Name your own picture 
greyim = rgb2gray(im);
I = edge(greyim,'canny',[0.03 0.19]);
I = double(I)/255;  
 
%or threshold its edges
figure(1)
E = edge(I);  
imshow(I);
title ('Gradient Image')
imwrite(E,'HoughlineE.jpg','jpg','Quality',60)

 
%Use matlab's radon function to compute the hough transform:
theta = (0:179)';
[R,xp] = radon(E,theta);
figure(2)
imagesc(theta,xp,R), colorbar;
xlabel ('theta (deg)'), ylabel ('rho (pixels from center)')
title('Line Space');  

plt=1;
for peaks = 80:10:100 % Change the threshold check the diferent result 
   %Find peaks in the linespace:
	i = find(R>peaks);
	%Sort the output and pick the top lines
	[foo,ind] = sort(-R(i));
	k = i(ind(1:size(i)));  
	%Convert linear index into coordinates of the peaks.
	[y,x] = ind2sub(size(R),k); 
 
	%Find the theta and rho values for the peak coordinates.
	t = -theta(x)*pi/180;
	r = xp(y);  
	%The lines parameters are computed as follows.  The line parameters have the coefficents of the equation Ax + By + C = 0,
	%and are invarient to scale.  However this particular scaling will produce the distance of a point to the line with the dot
	%product: [A; B; C]' * [x; y; 1].
	lines = [cos(t) sin(t) -r];  

	%Transform the line from the center of the image to the upper left.  (The minus ones are for matlabs 1 based coordinates.)
	cx = size(I,2)/2-1;
	cy = size(I,1)/2-1;
	lines(:,3) = lines(:,3) - lines(:,1)*cx - lines(:,2)*cy;  

	%Here are the top lines drawn on the gradient image
	figure(3+plt)
   text1='Gradient image with lines - Accumulator > ';
   text2=num2str(peaks);
   text=strcat(text1,text2);
	imshow(E); title(text);
   draw_lines(lines);
   plt=plt+1;
end