%HOUGH TRANSFORM****************************************

%Hough Line Detection
%Arguments input image

function[output]=Hough(img)
figure;imagesc(img);colorbar;

%Local Variables---------------------
dimens=size(img);
height=dimens(1);
width=dimens(2);
matrixGx=[-1 0 1;-2 0 2;-1 0 1];
matrixGy=[1 2 1;0 0 0;-1 -2 -1];

hfactorAC=50;

heightAC=round(2*pi*hfactorAC);
widthAC=round(sqrt(width*width+height*height));

gradientmatrix=zeros(height,width);
accumulator=zeros(heightAC,widthAC);
accumulator2=zeros(heightAC,widthAC);
output=zeros(height,width);
summatrix=zeros(height,width);

nolocs=3;
maxlocs=zeros(nolocs,2);
locsize=25;
locoff=floor(locsize/2);

inters=cell(1,1);
newinters=cell(1,1);
finalinters=cell(1,1);

%-----------------------------------

for x=2:(width-1)
for y=2:(height-1)
gx=(sum(sum(matrixGx.*img((y-1):(y+1),(x-1):(x+1)))))/8;
gy=(sum(sum(matrixGy.*img((y-1):(y+1),(x-1):(x+1)))))/8;

gradientmatrix(y,x)=sqrt(gx*gx+gy*gy);
end
end

figure;imagesc(gradientmatrix); colorbar;

edgeThreshold=input('Enter Gradient Threshold: ');
%edgeThreshold=500;

disp('Calculating Accumulator Array');

for x=2:(width-1)
for y=2:(height-1)
t=gradientmatrix(y,x);
if gradientmatrix(y,x)>edgeThreshold

curvevalues=zeros(1,heightAC);

for row=1:heightAC
theta=pi*2*(row/heightAC);
r=(x*cos(theta))+(y*sin(theta));
p=uint16(theta*hfactorAC);
q=uint16(r);
if p>0
if p<=heightAC
if q>0
if q<=widthAC
accumulator(p,q)=accumulator(p,q)+gradientmatrix(y,x);
curvevalues(row)=accumulator(p,q);
end
end
end
end
end
[curvemax,peak]=max(curvevalues);
theta=pi*2*(peak/heightAC);
r=(x*cos(theta))+(y*sin(theta));
p=uint16(theta*hfactorAC);
q=uint16(r);
accumulator2(p,q)=accumulator2(p,q)+gradientmatrix(y,x);

end
end
end

%figure;imagesc(accumulator);colorbar;figure;imagesc(accumulator2);colorbar;

disp('Retrieving Lines');

for n=1:nolocs
[a,b]=max(accumulator2);
[c,d]=max(a);
j=b(d);
k=d;
maxlocs(n,1)=j;
maxlocs(n,2)=k;

accumulator2(j,k)=0;

if (k-locoff>0) & (k+locoff<=widthAC) & (j-locoff>0) & (j+locoff<=heightAC)
accumulator2(j-locoff:j+locoff,k-locoff:k+locoff)=0;
else
for x=-locoff:locoff
for y=-locoff:locoff
try
accumulator2(j-locoff+y:j+locoff+y,k-locoff+x:k+locoff+x)=0;
catch
end
end
end

end
end

disp('Calculating Areas');

figure;imagesc(accumulator2); colorbar;
divmatrix=cell(1,(2*nolocs));
linematrix=cell(1,nolocs);

for n=1:nolocs
divmatrix{1,2*n-1}=zeros(height,width);
divmatrix{1,2*n}=zeros(height,width);
linematrix{1,n}=zeros(height,width);

x=maxlocs(n,2);
y=maxlocs(n,1);

for b=1:width
innards=y/hfactorAC;
a=uint16((x-(b*cos(innards)))/sin(innards));
if a>0
if a<=height
linematrix{1,n}(a,b)=1;
end
end
end

for a=1:height
innards=y/hfactorAC;
b=uint16((x-(a*sin(innards)))/cos(innards));
if b>0
if b<=width
linematrix{1,n}(a,b)=1;
end
end
end
%figure;imagesc(linematrix{1,n});colorbar;

theta=y/hfactorAC;
if((theta>(1.5*pi)) & (theta<=(2*pi))) | ((theta>0) & (theta<=(pi/2)))
cstart=1;
cend=height;
dir=1;
else
if (theta>(pi/2)) & (theta<=pi)
cstart=height;
cend=1;
dir=-1;
end
end

for i=cstart:dir:cend
if linematrix{1,n}(i,1)==1
break
end
for j=1:width
if linematrix{1,n}(i,j)==0
divmatrix{1,2*n-1}(i,j)=1;
else
break
end
end
end

divmatrix{1,2*n}=~divmatrix{1,2*n-1};
%divmatrix{1,2*n}=~(divmatrix{1,2*n-1}+linematrix{1,n});
%figure;imagesc(divmatrix{1,2*n-1});colorbar; figure;imagesc(divmatrix{1,2*n});colorbar;

end

disp('Calculating Intersections');

s=1;
for r=1:(2*nolocs)

for q=1:(2*nolocs)

ident=0;

if q~=r
temp=divmatrix{1,r}.*divmatrix{1,q};
%figure;imagesc(temp);colorbar;
if nnz(temp)~=0
%if isempty(nonzeros(temp))==0
%intersection
if isempty(inters{1,1})==0 & s>1
for p=1:(s-1)
if inters{1,p}==temp
ident=1;
break
end
end
end

if ident==0
inters{1,s}=temp;
%figure;imagesc(inters{1,s});colorbar;
s=s+1;
end

end
end

end
end

disp('Removing subintersections');

%newinters=inters;
a=1;
for r=1:(s-1)
for q=1:(s-1)
if q~=r

subinters=inters{1,q}.*inters{1,r};
if nnz(subinters)~=0 %& nnz(newinters{1,r})~=0
%figure;imagesc(inters{1,q});colorbar;figure;imagesc(inters{1,r});colorbar;
%figure;imagesc(subinters);colorbar;
temp=double(uint16(inters{1,r}-subinters));
if nnz(temp)~=0
newinters{1,a}=temp;
%figure;imagesc(newinters{1,p});colorbar;
a=a+1;
end

end

end
end
end

disp('Removing duplicates');

c=1;
for n=1:(a-1)

ident=0;

if n>1
for p=1:(n-1)
if newinters{1,p}==newinters{1,n}
%figure;imagesc(newinters{1,p});colorbar;
%figure;imagesc(newinters{1,n});colorbar;

ident=1;
break
end
end
end

if ident==0 %& nnz(newinters{1,n})~=0
finalinters{1,c}=newinters{1,n};
%figure;imagesc(finalinters{1,c});colorbar;

c=c+1;
end
end

suminters=zeros(height,width);
for n=1:c-1
suminters=suminters+(n*finalinters{1,n});

end
figure;imagesc(suminters);

disp('Equalising Image');

means=zeros(1,(c-1));

for n=1:(c-1)
means(1,n)=mean(nonzeros(finalinters{1,n}.*img));
if n>1
diff=means(1,1)/means(1,n);
outimg=diff*(finalinters{1,n}.*img);
else
outimg=(finalinters{1,n}.*img);
end
output=output+outimg;
%figure;imagesc(outimg);caxis([0 5e4]);colorbar;

end

figure;imagesc(output);colorbar;