%% Hough transform example 1

clear all
close all

% Start by making an image consisting of 5 points.

im=zeros(16);
im(1,1)=1;
im(5,5)=1;
im(8,8)=1;
im(1,15)=1;
im(15,1)=1;
figure
imshow(im)

% The question is, are these point on no, one or multiple lines?
% That is, is there any line y=ax+b that passes through at least 3 of these
% points for some values of a and b?
%
% To investigate this, let's first take the (x,y) coordinates of every
% point in the given image and plug it into the equation b=-xa+y (note that
% (x,y) should be parameters in this equation), and then vary a over some
% predefined interval (we will use integers between -5 and 5). This would
% produce lines in terms of a and b, and we will visualise these lines.

figure

% Point 1, x=1, y=1, b=-a+1
a=-5:5
b=-a+1
plot(a,b,'r')
hold on

% Point 2, x=5, y=5, b=-5a+5
a=-5:5
b=-5*a+5
plot(a,b,'g')

% Point 3, x=8, y=8, b=-8a+8
a=-5:5
b=-8*a+8
plot(a,b,'b')

% Point 4, x=1, y=15; b=-a+15
a=-5:5
b=-a+15
plot(a,b,'k');

% Point 5 x=15, y=1; b=-15a + 1
a=-5:5
b=-15*a+1
plot(a,b,'y');

grid on
axis([-5 5 -30 30])
xlabel('a')
ylabel('b')

% So it appears that the points form two lines, one with parameters a=1 and
% b=0, the other with parameters a=-1 and b=16. This could be visually
% verified by looking at the image of the points. It could also be
% interpreted directly from the (x,y) coordinates of the 5 points.
%
% Let's now look at how the accumulator matrix looks like when using the
% Cartesian representation. We will quantize a as integers between -5 and 5
% and b as integers between -100 and 100.

acc=zeros(11,201);
figure

a=-5:5
b=-a+1
for x=1:length(a)
    acc(a(x)+6,b(x)+101)=acc(a(x)+6,b(x)+101)+1;
end
imshow(acc,[])
xlabel('b')
ylabel('a')

a=-5:5
b=-5*a+5
for x=1:length(a)
    acc(a(x)+6,b(x)+101)=acc(a(x)+6,b(x)+101)+1;
end
imshow(acc,[])
xlabel('b')
ylabel('a')

a=-5:5
b=-8*a+8
for x=1:length(a)
    acc(a(x)+6,b(x)+101)=acc(a(x)+6,b(x)+101)+1;
end
imshow(acc,[])
xlabel('b')
ylabel('a')

a=-5:5
b=-a+15
for x=1:length(a)
    acc(a(x)+6,b(x)+101)=acc(a(x)+6,b(x)+101)+1;
end
imshow(acc,[])
xlabel('b')
ylabel('a')

a=-5:5
b=-15*a+1
for x=1:length(a)
    acc(a(x)+6,b(x)+101)=acc(a(x)+6,b(x)+101)+1;
end
imshow(acc,[])
xlabel('b')
ylabel('a')

figure
as = repmat((-5:5)', 1, 201);
bs = repmat(-100:100, 11, 1);
mesh(as, bs, acc)

w_half=2;
acc(6+1+(-w_half:w_half), 101+ 0+(-w_half:w_half))
acc(6-1+(-w_half:w_half), 101+16+(-w_half:w_half))

% So we get maxima in the accumulator matrix at two locations, each
% corresponding to the two lines we previously identified (one with
% parameters a=1 and b=0, the other with parameters a=-1 and b=16).
%
% Let's finally look at how the accumulator matrix looks like when using
% polar representation. We will quantize theta as integers between -89 and
% 90 degrees, while rho as integers between -23 and 23 (this range for rho
% was chosen because sqrt(16^2+16^2) = sqrt(256+256) = sqrt(512) is between
% 22.5 and 23).

acc=zeros(47,180);
figure

theta=-89:90
rho=round(cosd(theta)+sind(theta))
for x=1:length(theta)
    acc(rho(x)+24,theta(x)+90)=acc(rho(x)+24,theta(x)+90)+1;
end
imshow(acc,[])
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')

theta=-89:90
rho=round(5*cosd(theta)+5*sind(theta))
for x=1:length(theta)
    acc(rho(x)+24,theta(x)+90)=acc(rho(x)+24,theta(x)+90)+1;
end
imshow(acc,[])
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')

theta=-89:90
rho=round(8*cosd(theta)+8*sind(theta))
for x=1:length(theta)
    acc(rho(x)+24,theta(x)+90)=acc(rho(x)+24,theta(x)+90)+1;
end
imshow(acc,[])
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')

theta=-89:90
rho=round(cosd(theta)+15*sind(theta))
for x=1:length(theta)
    acc(rho(x)+24,theta(x)+90)=acc(rho(x)+24,theta(x)+90)+1;
end
imshow(acc,[])
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')

theta=-89:90
rho=round(15*cosd(theta)+sind(theta))
for x=1:length(theta)
    acc(rho(x)+24,theta(x)+90)=acc(rho(x)+24,theta(x)+90)+1;
end
imshow(acc,[])
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')

figure
rhos = repmat((-23:23)', 1, 180);
thetas = repmat(-89:90, 47, 1);
mesh(rhos, thetas, acc)

w_half=5;
rho_line2=round(8*cosd(45)+8*sind(45))
acc(24+        0+(-w_half:w_half), 90-45+(-w_half:w_half))
acc(24+rho_line2+(-w_half:w_half), 90+45+(-w_half:w_half))

% So we get maxima in the accumulator matrix at two locations. The line
% with parameters rho=0 and theta=-45 degrees corresponds to the line with
% parameters a=1 and b=0, while the line with parameters
% rho=round(8*cosd(45)+8*sind(45))=round(8*sqrt(2)/2+8*sqrt(2)/2)=round(8*sqrt(2))=round(11.3)=11
% and theta=45 degrees corresponds to the line with parameters a=-1 and b=16.


%% Hough transform example 2

clear all
close all

% Start by making an image consisting of two lines with 7 connected points.

i=zeros(12);
for j=3:9
  i(j,j)=1;
  i(j+3,j)=1;
end

figure
colormap(gray(2))
imagesc(i);
grid on
set(gca,'xcolor','w')
set(gca,'ycolor','w')
axis image

% Use the MATLAB function HOUGH to find the parameters of these lines using
% polar representation.

[acc,thetas,rhos]=hough(i);
figure
imagesc(thetas,rhos,acc)
colormap(cool)
colorbar
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')
grid on

% We see that the accumulator matrix has two maxima, one corresponding to a
% line with parameters rho=0 and theta=-45 degrees, while the other
% corresponds to a line with parameters rho=-2 and theta=-45 degrees.
% Looking at the image of the points (or the (x,y) coordinates of the
% points), we see that this is correct except that the latter line would
% have rho=2 in our definition of the polar representation (because the
% point on the latter line closest to origo is (1.5,-1.5) and sqrt(1.5^2 +
% 1.5^2) is about 2.1). The reason for this difference in sign is that the
% HOUGH function uses a definition of the angle described by theta that
% differs from the definition used in the lecture on Hough transform;
% strictly speaking, it is not the definition of angle that differs, but
% rather what is defined to be the x-axis and y-axis (the HOUGH function
% uses that the x-axis is a horizontal line increasing from coordinate
% (1,1) to the right, while the y-axis is a vertical line increasing from
% coordinate (1,1) and down).
%
% Let's now remove 2 points from each line in the original image to
% investigate whether this hampers our ability to detect the lines.

i(5,5)=0;
i(7,7)=0;

i(4+3,4)=0;
i(7+3,7)=0;

figure
colormap(gray(2))
imagesc(i);
grid on
set(gca,'xcolor','w')
set(gca,'ycolor','w')
axis image

% Again use the MATLAB function HOUGH to find the parameters of these lines
% using polar representation.

[acc,thetas,rhos]=hough(i);
figure
imagesc(thetas,rhos,acc)
colormap(cool)
colorbar
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')
grid on

% We see that the accumulator matrix is very similar as in the case where
% the lines consisted of 7 connected points. In particular, it still has
% the same two maxima.


%% Hough transform example 3

clear all
close all

% Start by reading an image of a corridor.

i=imread('corridor.png');
ig=double(rgb2gray(i));

figure
imshow(ig,[])

% Threshold the gradient magnitude image to obtain an edge image.

h1=fspecial('sobel');
h2=h1';
igh=imfilter(ig,h1,'replicate');
igv=imfilter(ig,h2,'replicate');
igm=sqrt(igh.^2+igv.^2);
figure
imshow(igm,[])

igmT=igm>120;
figure
imshow(igmT)

% Use the MATLAB function HOUGH to find the parameters of lines in the
% image using polar representation.

[igsH,thetas,rhos]=hough(igmT);

figure
s=surf(igsH);
set(s,'edgecolor','none')
shading interp
lighting phong
camlight right
set(gcf,'color','k')
set(gca,'color','k')

figure
imagesc(thetas,rhos,igsH)
colormap(hot)
colorbar
xlabel('\theta (degrees)')
ylabel('\rho (pixels from center)')
grid on

% Let's first use the MATLAB function HOUGHPEAKS to find the peaks in the
% accumulator matrix and then the MATLAB function HOUGHLINES to detect line
% segments in the image using the detected peaks.

peaks=houghpeaks(igsH,8,'threshold',ceil(0.3*max(igsH(:))));
hold on
x=thetas(peaks(:,2));
y=rhos(peaks(:,1));
plot(x,y,'s','color','white');

lines=houghlines(igmT,thetas,rhos,peaks,'FillGap',20,'MinLength',40);
figure
imshow(ig,[])
axis on
hold on
for k = 1:length(lines)
   xy = [lines(k).point1; lines(k).point2];
   plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');
   plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
   plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
end

% Let's now instead try to manually detect the peaks and then draw these as
% full lines in the image. Note that we will include more peaks here (20
% instead of 8) because HOUGHPEAKS will by default suppress the
% neighborhood around each peak, but this is not done in the manual peak
% detection code below.

ind=find(igsH>ceil(0.3*max(igsH(:))));
[~,index]=sort(-igsH(ind));
k=ind(index(1:20));
[r,c]=ind2sub(size(igsH),k);
th=thetas(c);
rho=rhos(r);

figure
imshow(ig,[])
axis on
hold on
x = 1:size(ig,2);
y = (repmat(rho',1,size(x,2)) - (repmat(x,length(k),1)-1).*repmat(cosd(th)',1,size(x,2))) ./ repmat(sind(th)',1,size(x,2)) + 1;
plot(x,y)