% called by EyeballS_FINAL
% to compute the autocorrelation (R) curves for each image.
% by Dave Rubin (USGS), Erin Todd (UCSC), and Jon Warrick (USGS)
%
% contact:  jwarrick@usgs.gov
%
% Version 2.0
% March 30, 2009 (Date of Last Modification)


%Get image dimensions and data
ImageHeight = size(data,1);
ImageWidth = size(data,2);

ImageWidthToProcess = ImageWidth-MaxOffset;
ImageHeightToProcess = ImageHeight-MaxOffset;  

data1 = data(:,1:ImageWidthToProcess);
data_V1 = data(1:ImageHeightToProcess,:);   

MinAutoC = 1.0; %resets this variable for threshold test
clear autoc1;
clear autoc_V1;     

clear pixSizeFromImage;
clear pixSizeFromImage_V;   



for j = 1:length(limits);       
    MinAutoC = 1.0; %resets this variable for threshold test
    clear autoc1;
    clear pixSizeFromImage;
    PixelOffset = 0;
    i = 0;
    
% Horizontal Autocorrelation Analysis

    while MinAutoC > limits(j) & PixelOffset <= length(CalibData)-PixelStep ; %.3 normal Stop calculations when autocorrelation is too small.
        i = i+1;
        
        %resize the two samples of the image to compare with autocorr
        PixelOffset = 1 + (i-1) * PixelStep;     
        data1=1:ImageHeight*ImageWidthToProcess;
        data2=1:ImageHeight*ImageWidthToProcess;
        data1(1:ImageHeight*ImageWidthToProcess)= reshape (data(1:ImageHeight,1:ImageWidthToProcess),ImageHeight*ImageWidthToProcess,1);
        data2(1:ImageHeight*ImageWidthToProcess)= reshape (data(1:ImageHeight,PixelOffset+1:ImageWidthToProcess+PixelOffset),ImageHeight*ImageWidthToProcess,1);
       
        %Compute the autocorr (R)
        correl= corrcoef(data1,data2);
        autoc1(i)=correl(1,2);  
        
        %Use R value to estimate GS for this pixel step using interp.
        pixSizeFromImage(i) = interp1(CalibData(PixelOffset,:), ymm, autoc1(i),'linear');
        MinAutoC = min(MinAutoC,autoc1(i));
    end

        %correct for values that are underresolved
        %this is generally not necessary and can be deleted, unless the 'limits' chosen is ~0.
    for k=1:i
        if isnan(pixSizeFromImage(k))==1;
            pixSizeFromImage(k)=[];
        end
    end
    
    
% Vertical Autocorrelation Analysis

    MinAutoC = 1.0; %resets this variable for threshold test
    clear autoc_V1
    clear pixSizeFromImage_V
    PixelOffset = 0;
    i = 0;
    
    while MinAutoC > limits(j) & PixelOffset <= length(CalibData)-PixelStep ; %.3 normal Stop calculations when autocorrelation is too small.
        i = i+1;      
        
        %resize the two samples of the image to compare with autocorr
        PixelOffset = 1 + (i-1) * PixelStep; 
        data_V1=1:ImageWidth*ImageHeightToProcess;
        data_V2=1:ImageWidth*ImageHeightToProcess;
        data_V1(1:ImageWidth*ImageHeightToProcess)= reshape (data(1:ImageHeightToProcess,1:ImageWidth),ImageHeightToProcess*ImageWidth,1);
        data_V2(1:ImageWidth*ImageHeightToProcess)= reshape (data(PixelOffset+1:ImageHeightToProcess+PixelOffset,1:ImageWidth),ImageHeightToProcess*ImageWidth,1);
          
        %Compute the autocorr (R)
        correl_V= corrcoef(data_V1,data_V2);
        autoc1_V(i)=correl_V(1,2);

        %Use R value to estimate GS for this pixel step using interp.
        pixSizeFromImage_V(i) = interp1(CalibData(PixelOffset,:), ymm, autoc1_V(i),'linear');
        MinAutoC = min(MinAutoC,autoc1_V(i));
    end
    
        %correct for values that are underresolved
        %this is generally not necessary and can be deleted, unless the
        %'limits' chosen is ~0.
    for k=1:i
        if isnan(pixSizeFromImage_V(k))==1;
            pixSizeFromImage_V(k)=[];
        end
    end
    
    % Compute Summary statistics from the Analyses...
    GrainSizePix = mean(pixSizeFromImage); 
    std_pix=std(pixSizeFromImage);    
    
    GrainSizePix_V = mean(pixSizeFromImage_V); 
    std_pix_V=std(pixSizeFromImage_V);
    
    GrainSizePix_av(1)=GrainSizePix;
    GrainSizePix_av(2)=GrainSizePix_V;
    GrainSizePix_av = mean(GrainSizePix_av);
    
    std_pix_av(1)=std_pix;
    std_pix_av(2)=std_pix_V;
    std_pix_av = mean(std_pix_av);
       
    % Write calculated grain size in command window to let user know
    % about analysis progress
    
    ['For photo file' FileName ', calc GS(pix) is:  ' num2str(GrainSizePix_av(1,j))]
end