function Sample_cross_sections
%Program to compute cross sections for individual samples (samples are warped to the mean shape of their relevant time point for alignment purposes)
%Next program to run: Mean_cross_sections

close all;
all=[];
Path1=uigetdir(cd,'Select Reference data folder');%AG July 2013: for saving mean leaf data
Pathname=uigetdir(cd,'Select main project folder in  which all shape results and figures are stored');%AG July 2013: for saving mean leaf data
resdir='Cross Sections of all samples (warped)';
if ~exist([Pathname '\' resdir],'dir')
    mkdir(Pathname,resdir);
end
Pathname=[Pathname '\' resdir];



resfold=[Path1 '\SourceData'];

cd(resfold);
listpots=dir;
for  pot=3:size(listpots,1); %for each pot
    cd([resfold '\' listpots(pot).name]);
    listplts=dir;
    for plt=3:size(listplts,1); %for each plant
        cd([resfold '\' listpots(pot).name '\' listplts(plt).name]);
        listdas=dir;
        LWplant=[];
        for day=3:size(listdas,1) %for each successive pair of DAS
            cd([resfold '\' listpots(pot).name '\' listplts(plt).name '\' listdas(day).name]);
            disp(['Processing ' listpots(pot).name '\' listplts(plt).name '\' listdas(day).name])
            load([pwd '\3Daligned_data.mat']); XYZ=XYZg;
            load('warpedleaf3D.mat'); %load xyzgW 
            XYZg(:,1)=xyzgW(:,1); XYZg(:,2)=xyzgW(:,2);
            
            %Make set of points along the length of the leaf
            yi=(min(XYZg(:,2)):(max(XYZg(:,2))-min(XYZg(:,2)))/100:max(XYZg(:,2)))';
            xi=zeros(size(yi,1),1);
            F = TriScatteredInterp(XYZg(:,1),XYZg(:,2),XYZg(:,3)); %Interpolate the leaf surface based on the mesh points
            zi=F(xi,yi); %Use this to get the z value of each point down the midline
            xi(isnan(zi))=[]; %The min and max may not necessarily be at the midline, so there ends up being NaNs where the line goes past the surface
            yi(isnan(zi))=[];
            zi(isnan(zi))=[];
            yl=yi; %To save
            xl=xi;
            zl=zi;
            
            %Make set of points across the width of the leaf
            xi=(min(XYZg(:,1)):(max(XYZg(:,1))-min(XYZg(:,1)))/100:max(XYZg(:,1)))';
            yi=zeros(size(xi,1),1);
            F = TriScatteredInterp(XYZg(:,1),XYZg(:,2),XYZg(:,3)); %Interpolate the leaf surface based on the mesh points
            zi=F(xi,yi); %Use this to get the z value of each point down the midline
            xi(isnan(zi))=[]; %The min and max may not necessarily be at the midline, so there ends up being NaNs where the line goes past the surface
            yi(isnan(zi))=[];
            zi(isnan(zi))=[];
            yw=yi; %To save
            xw=xi;
            zw=zi;

             save([Pathname '\cross_sections_' listpots(pot).name  listplts(plt).name listdas(day).name '.mat'],'xl','yl','zl','xw','yw','zw'); 
        end
    end
end
disp(['Cross sections for all samples have been saved to ' Pathname]);

