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

%Hough Line Detection
%Arguments input image, median filter size

function[output]=Hough(inimg,fsize)
dimens=size(inimg);
height=dimens(1);
width=dimens(2);
medianoffset=floor(fsize/2);
finalimg=inimg(medianoffset+1:height-medianoffset,medianoffset+1:width-medianoffset);

figure;imagesc(inimg);colorbar;
img=medianfilter(inimg,fsize);
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=10;
maxlocs=zeros(nolocs,2);
locsize=25;
locoff=floor(locsize/2);

areas=cell(1,1);
finalinters=cell(1,1);

%-----------------------------------
disp('Calculating Gradient');

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: ');


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

for n=1:nolocs
   summatrix=summatrix+(n*divmatrix{1,2*n});
end

figure;imagesc(summatrix);colorbar;

s=1;
areas{1,1}=zeros(height,width);
summax=max(max(summatrix));

disp('Spliting Areas');

for n=1:(summax+1)
   if n>1
   	foundlast=found(1,2*width-1);
   	if foundlast==1
      	s=s+1;
      	areas{1,s}=zeros(height,width);
   	end
   end
   found=zeros(1,2*width-1);
   
   
   for d=(-1*(width-1)):1:(width-1)
      
      col=diag(summatrix,d);
      [len,a]=size(col);
      for p=1:len
         if col(p,1)==(n-1)
            %is area in diag
            if d<1
               x=p;
               y=p+(-1*d);
            else   
               x=p+d;
               y=p;
            end
            
            areas{1,s}(y,x)=1;
            found(1,(d+width))=1;
         end
      end
      
      if (d+width)>1
         if found(1,(d+width))==0 & found(1,(d+width-1))==1
            s=s+1;
            areas{1,s}=zeros(height,width);
         end
               
      end
      
   end
   
end
sumareas=zeros(height,width);
%for m=1:s-1   figure;imagesc(areas{1,m});colorbar;end
for n=1:s-1
   sumareas=sumareas+n*(areas{1,n});
end

figure;imagesc(sumareas);colorbar;

disp('Equalising Image');

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


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


figure;imagesc(output);colorbar;