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

%Hough Line Detection
%Arguments input image,edge threshold,Hough threshold

function[maxlocs]=Hough(img,edgeThreshold)

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);
linematrix=zeros(height,width);
areamatrix1=zeros(height,width);

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;

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

if gradient>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)+gradient;
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)+gradient;

end
end
end

%figure;imagesc(gradientmatrix); colorbar;

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

nolocs=1;
maxlocs=zeros(nolocs,2);
locsize=11;
locoff=floor(locsize/2);

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

%figure;imagesc(accumulator2); colorbar;

for n=1:nolocs
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(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(a,b)=1;
end
end
end

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(i,1)==1
break
end

for j=1:width
if linematrix(i,j)==0
areamatrix1(i,j)=1;
else
break
end
end
end
areamatrix2=~(areamatrix1+linematrix);

mean1=mean(nonzeros(areamatrix1.*img))
mean2=mean(nonzeros(areamatrix2.*img))

end

figure;imagesc(linematrix); colorbar;