% This performs autocorrelation and mean grain-size analyses on a series of images 
% All images should be saved as .TIF and stored in a single folder.
% 
% This folder must also contain a file (mm_pix.txt), which contains the
% linear scaling factors for each image in alphanumeric order.
% To build the mm_pix.txt, a Matlab function (measure_mm2pix.m) has been
% written by Josh Logan and is available for 
% 
% Calibration data is provide in the files:
%     JAW_calibdata_150.txt & ypix_JAW.txt  
% These data extend through a max pixel offset of only 150 pixels. Note
% that these data were generated from and for the Elwha River delta, and
% although they may be applied to other sites, there may be systematic bias
% in this application that should be investigated and corrected for.  
% 
% Uses Calibrations.m and AutoCorrS_Vert.m.
% Developed by David Rubin (USGS), Erin Todd (UCSC), Josh Logan (USGS), and
% Jon Warrick (USGS)
%
% contact:  jwarrick@usgs.gov
%
% Version 2.0
% March 30, 2009 (Date of Last Modification)



%Preliminary clean up
close; fclose('all');
clear

%Determine platform, adjust path dividers (jlogan).
if ismac 
    slashdir = '/';
elseif ispc
    slashdir = '\';
end

%Load calibration data
load JAW_calibdata_150.txt;
CalibData=JAW_calibdata_150;
MaxOffset = 150;

load ypix_JAW.txt;
ymm=ypix_JAW;
str_out='pix';

% run_type=input('What type of run?  O for vary pixel step, 1 for vary autocorr limits [pix step at 3]:');

run_type=1;

%++++++++++++++++++++++++++++
% These vairables will allow users to modify the kind of analysis
% conducted:
% 'PixelStep' influences the analysis speed and (1 = slowest, but most
%      complete analysis.
% 'limits' assigns a value for R below which the analysis will stop.
%      Warrick et al. (in press) ESPL found an optimized value at 0.25.
% 'bia_correct' will modify the final grain-size results by a systematic
%      bias (Warrick et al found 2% for Elwha (i.e., = 1.02).  Note that this
%      value may change significantly with application to other sites/settings 
%      as shown by the Kachemak data in Warrick et al.  

PixelStep = 1

limits=0.25

bias_correct = 1.02
%+++++++++++++++++++++++++++


run_mode = 0;

d=uigetdir('','Select folder location of photo batch');
dir_index=strfind(d, slashdir);  fldr_index=length(dir_index);
exp_name=strcat(d(((dir_index(fldr_index))+1):length(d)),'_out.csv');
if d
else
    return
end
filecheck = exist(exp_name);
if filecheck==2;
    display('An output file with this name already exists');
    return
end

%Find the scaling data and TIF images
path(path,d);
cd(d);
load mm_pix.txt;    

dir_data=dir(fullfile(d,'*.tif'));
num_files=length(dir_data);

%DOUBLE CHECK THAT THERE IS A MM_PIX FOR EACH PHOTO (jlogan)
if length(mm_pix) ~= length(dir_data)
    disp('Length of mm_pix.txt file does not match number of images')
    return
end


%%%%%%%%%%%%%%%%
% This is a for-loop that computes an autocorrelation R for each pixel step
% of each photo using the AutoCorrS_Vert.m function.  Each R value is fit to 
% a grain-size (by pixels) using linear regression with the ypix_JAW data. 
% The analysis is run using both horizontal and vertical shifting of
% pixels.

for f=1:length(dir_data);
    FileName=dir_data(f).name;
    ImageData=imread(FileName);
    
    %%%%%%%%%%%%%
    %       PlotGen_AC_Call;
    %        INSERT THE PlotGen_AC_Call COMMANDS:
    [M,N] = size(FileName);
    FileName = FileName(1,1:N-4);
    %         data=double(ImageData);
    data=ImageData;
    if run_type == 0;
        AutoCorrPixStep; % Calculate mean grain size.3
    elseif run_type == 1;
        %     AutoCorrS;
        cd ..
        AutoCorrS_Vert;
        cd(d)
    else
        break
        return
    end;
    
    %%%%%%%%%%%%
    
    FilesIndex(f,1)=cellstr(FileName);
    GS_pix(f,1)=GrainSizePix;
    GS_pix_V(f,1)=GrainSizePix_V;
    GS_pix_av(f,1)=GrainSizePix_av;
    
    sd_dev(f,1)=std_pix;
    sd_dev_V(f,1)=std_pix_V;
    sd_dev_av(f,1)=std_pix_av;
    
end

% Convert all results in pixels to mm and correct for bias.
GS_mm    =  GS_pix   .* mm_pix .* bias_correct;  %horizontal
GS_mm_V  =  GS_pix_V .* mm_pix .* bias_correct;  %vertical
GS_mm_av =  GS_pix_av.* mm_pix .* bias_correct;  %average of H and V

std_mm    = sd_dev   .* mm_pix;  
std_mm_V  = sd_dev_V .* mm_pix;
std_mm_av = sd_dev_av.* mm_pix;
SE_av     = sd_dev_av./ GS_pix_av;  % the SE is related to grain sorting (See Warrick et al.)


FilesIndex=char(FilesIndex);

%HeaderStr='Photo, GS pix(h),StDev(h),GS_mm(h),GS_Phi(h),GS pix(v),StDev(v),GS_mm(v),GS_Phi(v),GS pix(Ave),StDev(Ave),GS_mm(av)';
HeaderStr='Photo, GS pix(h),StDev(h),GS_mm(h),GS pix(v),StDev(v),GS_mm(v),GS pix(Ave),StDev(Ave),SE(Ave),GS_mm(av)';
exp_name

fid=fopen(exp_name,'a');

fprintf(fid,HeaderStr);

fprintf(fid,'\r');
for o=1:length(dir_data);
%    fprintf(fid,'%s,%6.2f,%6.2f,%5.1f,%5.3f,%6.2f,%6.2f,%5.1f,%5.3f,%6.2f,%6.2f,%5.1f,%5.3f',FilesIndex(o,:),GS_pix(o,:),std_mm(o,:),GS_mm(o,:),GS_phi(o,:),GS_pix_V(o,:),std_mm_V(o,:),GS_mm_V(o,:),GS_phi_V(o,:),GS_pix_av(o,:),std_mm_av(o,:),GS_mm_av(o,:),GS_phi_av(o,:));%,%5.1f,%5.3f,%5.2f,%5.2f,Meas_Baxis_mm(o,:),Meas_Baxis_phi(o,:),prcnt_AC_meas_mm(o,:),prcnt_AC_meas_phi(o,:)
    fprintf(fid,'%s,%6.2f,%6.2f,%5.1f,%6.2f,%6.2f,%5.1f,%6.2f,%6.2f,%6.2f,%5.1f',FilesIndex(o,:),GS_pix(o,:),std_mm(o,:),GS_mm(o,:),GS_pix_V(o,:),std_mm_V(o,:),GS_mm_V(o,:),GS_pix_av(o,:),std_mm_av(o,:),SE_av(o,:),GS_mm_av(o,:));
    fprintf(fid,'\r');
end
fclose(fid);