function [tall, fall, tauvalues, Dvalues]=EK_FRAPbatch

% batch processes FRAP data using OME-tiff files
%
% 1-Prompts file listing 
% 2-loads files in the list
% 3-Detects area that is to be bleached, calculates radius
% 4-Calculates intensity vs. time for the bleached area
% 5-calls EK_SoumpassisFit(t,f) to fit FRAP data to the Soumpassis equation
% and returns the fit parameter 'tau' (requires Curve Fitting Toolbox)
% 6-Displays distribution of tau parameters and statistics
% 7-Averages individual data files
% 7-Plots all data together and makes a fit to avg. data
% 
% Erdem Karatekin
% Last updated 15 December, 2015


% CHECK THESE VALUES, CHANGE AS NEEDED

% size of a pixel, in um 
pixSize = 0.265; % Olympus IX81, 60x/1.45 TIRF obj

tauStart= 30; % initial guess for fit, in sec. Look at a few curves to guess. 


        
[fn, fp]=uigetfile('*.txt','file listing data to be batch processed');
fullf=fullfile(fp,fn);
fid=fopen(fullf,'r'); %open read only
C = textscan(fid, '%s %f32', 'HeaderLines', 0);
%Cm = [C{:,1} C{:,2}]; %convert cell array to matrix
fclose(fid);

files = C{1,1}; % file names in a cell array
% relBleachInt = C{1,2}; %column vector of values

Nfile = size(files,1);
msg_nfiles=sprintf('Total # of Files: %d', Nfile) ;
disp(msg_nfiles)


[fndata, fpdata]=uigetfile('*.tif','select any file in the directory where original data are');
%fullfdata=fullfile(fpdata, fndata);
% filepath used later below

% -----------------------------------------
% --- Copy file names into variable nene (check the character # if it is not =20 change it ):

%ff=fopen([fp fn],'r');
%nene=textscan(ff,'%s');   %read file names into a 1x1 cell array. The only element in the array is a list of filenames. 
%fclose(ff);



fileindex=1;
NN=Nfile+1;

% initialize vectors into which data will be concatenated
tauvalues=zeros(length(Nfile),1);
Dvalues=zeros(length(Nfile),1);
legendtext=cell(1, length(Nfile)); % text cell array to use in fig legends
tall=cell(1,length(Nfile));
fall=cell(1,length(Nfile));

% Set up figure to receive datasets and fits
f_ = clf; % clear current fig, get handles
figure(f_);
%set(f_,'Units','Pixels','Position',[815 281 814 578]);
set(f_,'Units','Normalized','Position',[0.2 0.2 0.6 0.6]);
% legh_ = []; legt_ = {};   % handles and text for legend
% xlim_ = [Inf -Inf];       % limits of x axis
ax_ = axes;
set(ax_,'Units','normalized','OuterPosition',[0 0 1 1]);
set(ax_,'Box','on', 'FontSize', 14);
axes(ax_); 


xlabel('time, s', 'FontSize', 14)
ylabel('normalized fluorescence', 'FontSize', 14)
hold all; % toggles linestyles and markers as new plots are added


while NN > fileindex;
  %fname=char(nene{1}(fileindex));
  fname=char(files(fileindex));
     %name1 = sscanf(name1_, '%s') %read non-white-space characters, i.e. get rid of white space
    % name1=nene{1}(ind_file);
    msg_working=sprintf('Working on file: %s',fname);
    disp(msg_working);
    legendtext(fileindex)={fname};


     %___CALL THE PROGRAM/FILE WITH NAME name1...

     % open image file + metadata
     fullfdata=fullfile(fpdata,fname);
     
     reader = bfGetReader(fullfdata);
     omeMeta=reader.getMetadataStore();
     
     stackSizeX = omeMeta.getPixelsSizeX(0).getValue(); % image width, pixels
     stackSizeY = omeMeta.getPixelsSizeY(0).getValue(); % image height, pixels
     stackSizeT = omeMeta.getPixelsSizeT(0).getValue(); % time stack # of frames
     
     ch=1; % set channel
     z=1; % set z-plane
     
     tStack = zeros(stackSizeX,stackSizeY,stackSizeT); % not sure if X, Y shoudl be inverted
     frTimes = zeros(stackSizeT,1);
     
     for t=1:stackSizeT
         % get plane
         iPlane = reader.getIndex(z-1, ch-1, t-1)+1;
         tStack(:,:,t) = bfGetPlane(reader, iPlane); % this loads the entire time-stack
         frTimes(t) = omeMeta.getPlaneDeltaT(0,t-1).value;% args are image index, plane index
     end
   
    
     %% GET NOISE LEVEL
          
     % position has the form [xmin ymin width height].
     rectpos= [1, 1, floor(stackSizeX./4), floor(stackSizeY./4)];
     
     % create mask
     BW=false(stackSizeX,stackSizeY);
     
     for m=1:rectpos(3)
         for n=1:rectpos(4)
             BW(m,n)=true;
         end
     end
     
     noiseStats = regionprops(BW, tStack(:,:,1), 'MeanIntensity'); % get mean int. in the ROI box
     Fnoise = noiseStats.MeanIntensity;

     %% AUTO-DETECTED BLEACH AREA AND PIXEL VALUES
     %[Abl, intValues] = EK_getFRAPvals(tStack, relBleachInt, Fnoise);
     [Abl, intValues] = EK_getFRAPvals(tStack);
        
        t=frTimes;  % in sec
        y=intValues;
 
%         t=t-t(y==min(y)); % shift time axis such that t=0 is when y is min
%         %t=t.*dt; %convert t to time in sec (dt is period between successive acquisitions)
%         %t=double(t); %convert to double for use in curve fitting toolbox
%         %y=double(y); %convert y to double precision
% 
%         f=(y-min(y))./(y(end)-min(y)); % frac recovery, normalized to y_inf
%         %f2=(y-min(y))./(y(1)-min(y)); % frac recovery, normalized to 
%                  %pre-bleach value (takes first reading as pre-bleach value)
        
                 
%ignoring the frame acquired before bleach ended:                  
        t = t - t( find( y==min(y)) + 1 ); % shift time axis such that t=0 is when y is min, 
                % ignoring the frame acquired before the end of bleach period                 
        f = ( y - y(t==0) )./( y(end) - y(t==0) ); % frac recovery, normalized to y_inf         
                
                
         tall{:,fileindex}=t;
         fall{:,fileindex}=f;


        if isempty(Fnoise)~=1
            K=(y(1)-min(y))./(y(1)-Fnoise);
            disp({'K=', K})
        else %do nothing
        end
        
        plot(t,f,'LineWidth',2);
        
        % Make Soumpassis type fit
        [cf_, gof]=EK_SoumpassisFit(t,f,tauStart);
        
        disp(cf_);
        disp(gof);
        
        coefn=coeffnames(cf_); % get coeff names from cf object
        coefv=coeffvalues(cf_); % get coeff value from cf object
        tau=coefv(3);
        
        msg=sprintf('best fit coefficients are %s = %g, %s = %g, %s = %g', coefn{1}, coefv(1), coefn{2}, coefv(2), coefn{3}, coefv(3));
        disp(msg);
        disp('                                                    ')
        disp('****************************************************')
        
        tauvalues(fileindex)=tau;  %fill vector of tau values       
        
        %calculate xy diff coeff Dxy from tau, given
        %R_bleach=17 um on Nikon Ti Eclipse w/ 60x,1x
        % cf. On Olympus IX71 R_bleach was 9.25 um!!
        % Dxy=w^2/(4 tau), with w (radius of bleached spot) in um and tau s, get um^2/s
        Ablu = Abl .* (pixSize).^2; % convert bleach area from pixels to um^2
        D=Ablu /(4.* pi .*tau);
        Dvalues(fileindex)=D;  %fill vector of diff. coeffs.
                
        fileindex=fileindex+1; % increment file index
end

figure(f_);
legend(legendtext,'FontSize', 10);

% averaged data

avgtau=mean(tauvalues);
stdtau=std(tauvalues);
stdetau=stdtau/sqrt(Nfile);

disp('                                                           ')
disp('-----------------------------------------------------------')
disp('                  recovery times from fits                 ')
disp('-----------------------------------------------------------')
msg_res=sprintf('<tau>=%g s, std(tau)=%g s, stderr(tau)=%g s',avgtau,stdtau,stdetau);
disp(msg_res);

avgD=mean(Dvalues);
stdD=std(Dvalues);
stdeD=stdD/sqrt(Nfile);


disp('                                                           ')
disp('-----------------------------------------------------------')
disp('                     diffusion coefficients                ')
disp('-----------------------------------------------------------')
msg_res=sprintf('<D>=%g um^2/s, std(D)=%g um^2/s, stderr(D)=%g um^2/s',avgD, stdD, stdeD);
msg_res2=sprintf('using bleach area A=%g um^2', Ablu);
disp(msg_res);
disp(msg_res2);


%plot(D, f_D, 'ko','MarkerSize',6, 'Linewidth', 2);


% % call cftool functions to make double expt fit and display results:
% ok_ = isfinite(D) & isfinite(f_D);
% st_ = [0.5 -10 0.5 -0.1 ];  % initial guesses
% ft_ = fittype('exp2');
% 
% % Fit this model
% [cf_, gof, output] = fit(D(ok_),f_D(ok_),ft_,'Startpoint',st_);
% 
% plot(cf_)
% 
% disp(cf_)
% disp(gof)  ; % display goodness of fit parameters
% disp(output); % display fit parameters
% 
% 
% % Done plotting data and fits.  Now finish up loose ends.
% hold off;
% % leginfo_ = {'Orientation', 'vertical', 'Location', 'NorthEast'}; 
% % h_ = legend(ax_,legh_,legt_,leginfo_{:});  % create legend
% % set(h_,'Interpreter','none');
% xlabel(ax_,'time since docking, s', 'FontSize',20);               % put x label
% ylabel(ax_,'surviving fraction', 'FontSize',20);               % put y label
% % title(ax_, plottitle);



disp('-----------------------------------------------------------')
disp('                           PROGRAM IS OVER                 ')
disp('-----------------------------------------------------------')