%TextureEQ    Texture Equalisation Algorithm
%  TextureEQ(X) performs texture equalisation on image X
%  Returns Output
%
%  Comprising Sobel Edge detection, Hough Line Detection and Texture equalisation
%  Copyright Henry Bruce, University of Nottingham, 2003

function[output]=Hough(img)

figure;imagesc(img);colorbar;title('Input Image');	%Display Input image



%Local Variables----------------------------------------
dimens=size(img);

height=dimens(1);		%Image height and width
width=dimens(2);

matrixGx=[-1 0 1;-2 0 2;-1 0 1];		%Sobel edge detector kernels
matrixGy=[1 2 1;0 0 0;-1 -2 -1];

hfactorAC=50;			%Theta multiplication factor
heightAC=round(2*pi*hfactorAC);		%Height and width of Accumulator Array
widthAC=round(sqrt(width*width+height*height));

%Initialise Image Matrices

gradientmatrix=zeros(height,width);		%Edge gradient matrix
accumulator=zeros(heightAC,widthAC);		%First accumulator matrix
accumulator2=zeros(heightAC,widthAC);		%Second accumulator matrix
output=zeros(height,width);		%Output matrix
summatrix=zeros(height,width);		%Sum from CIA technique matrix
sumareas=zeros(height,width);			%Sum of all individual areas matrix

nolocs=10;		%No of peaks to locate, ie no of lines
maxlocs=zeros(nolocs,2);

locsize=25;		%Size of block to zero after locating peak in accumulator
locoff=floor(locsize/2);

%Initialise Cell Variables 
areas=cell(1,1);			%Separate detected areas cell
divmatrix=cell(1,(2*nolocs));			%Areas either side of lines cell
linematrix=cell(1,nolocs);			%Lines cell




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

for x=2:(width-1)		%Loop through image pixels calculating gradient magnitude and place in gradientmatrix
   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;title('Edge Gradient Image');		%Display Edge Detection array

edgeThreshold=input('Enter Gradient Threshold: ');		%Get user input for threshold




%Accumulator Calculation----------------------------------------
disp('Calculating Accumulator Array');

for x=2:(width-1)		%Loop though gradient pixels 
   for y=2:(height-1)
      if gradientmatrix(y,x)>edgeThreshold
         
         curvevalues=zeros(1,heightAC);	%Zero vector used to determine second accumulator
         
         for row=1:heightAC		%For all possible theta values calculate r
            theta=pi*2*(row/heightAC);
            r=(x*cos(theta))+(y*sin(theta));
            
            p=uint16(theta*hfactorAC);		%Convert to integers
            q=uint16(r);
            
            if p>0		%If r,theta exist in accumulator, accumulate by gradient at that point
               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);		%Calculate r and theta for max point along curve from 1st accumulator
         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);		%Increment second accumulator
      end
   end
end

figure;imagesc(accumulator);colorbar;title('First Accumulator Array');		%Display 1st and 2nd accumulators 
figure;imagesc(accumulator2);colorbar;title('Second Accumulator Array');




%Retrieve top points from Accumulator----------------------------------------
disp('Retrieving Lines');

for n=1:nolocs		%Get nolocs top peaks from 2nd accumulator
   [a,b]=max(accumulator2);
   [c,d]=max(a);   
   j=b(d);
   k=d;
   
   maxlocs(n,1)=j;		%Save r and theta in maxlocs array
   maxlocs(n,2)=k;
   
   accumulator2(j,k)=0;		%Zero area around peak to prevent inclusion of surrounding noise
   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




%Draw lines and areas----------------------------------------
disp('Calculating Areas');

for n=1:nolocs
   divmatrix{1,2*n-1}=zeros(height,width);	%Zero relevant line and areas matrices
   divmatrix{1,2*n}=zeros(height,width);
   linematrix{1,n}=zeros(height,width);

   x=maxlocs(n,2);	%Get peak location
   y=maxlocs(n,1);
   
   for b=1:width		%Draw line by working across and then down - to cover all possibilities
      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   

   theta=y/hfactorAC;		%Determine angle of line for creating area
   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	%'Fill' area to one side of line
      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}+linematrix{1,n});   	%Calculate area other side of line
   
end

for n=1:nolocs		%Perform cumulative intersection addition technique
   summatrix=summatrix+(n*divmatrix{1,2*n});
end

figure;imagesc(summatrix);colorbar;title('CIA Matrix');		%Display CIA matrix

s=1;
areas{1,1}=zeros(height,width);		%Get maximum in CIA matrix
summax=max(max(summatrix));




%Split CIA matrix into individual matrices----------------------------------------
disp('Spliting Areas');

for n=1:(summax+1)		%Work up to summax traversing the CIA array diagonally creating new matrices in areas cell
   
   if n>1		%Keep track of where the last area was found
   	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

for n=1:s-1			%Display areas found
   sumareas=sumareas+n*(areas{1,n});
end
figure;imagesc(sumareas);colorbar;title('Separate areas identified through CIA technique');




%Perform Equalisation----------------------------------------
disp('Equalising Image');

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

for n=1:(s-1)
   means(1,n)=mean(nonzeros(areas{1,n}.*img));		%Get first mean
   if n>1
      diff=means(1,1)/means(1,n);		%Compare subsequent means to first and equalise accordingly
      outimg=diff*(areas{1,n}.*img);
   else
      outimg=(areas{1,n}.*img);		%Set area in output image
   end
   output=output+outimg;
end

figure;imagesc(output);colorbar;title('Output Image');		%Display Output