function Mean_cross_sections
%Program to compute cross sections for the mean and the mode leaves.
%next program to run: Show_cross_sections

days=7:19;

%load mean outline
Path1=uigetdir(cd,'Select Reference data folder');%AG July 2013: for saving mean leaf data
load([Path1 '\meanoutlines_DAS7-19.mat']);

% [fn2, pn2]=uigetfile(pwd,['Load mean outlines file']);
% load([pn2 '\' fn2]);

%load mean and mode leaf data
%load('C:\Documents and Settings\arolland\My Documents\Dropbox\Lauren_shapepaper\Sep 25 - new mean leaf data\MeanleavesDAS7_19cSep24v5testborder.mat');
%SEP 29: GOOD ONE: load('C:\Users\Administrator\Desktop\TESTMEANLEAF_SEP27_noextrap_borderclassok_V9\MeanleavesDAS7_19cSep27v9noextra.mat')
Pathname=uigetdir(cd,'Select main project folder in  which all shape results and figures are stored');%AG July 2013: for saving mean leaf data
load([Pathname '\3D surface and figures\Mean3Dsurface.mat']);
%[fn2, pn2]=uigetfile(pwd,['Load mode leaf surface file']);
%load([pn2 '\' fn2]);
%this loads 'leaf_XIg','leaf_YIg','Mode_leaf_ZIg' and 'Truemean_leaf_ZIg'

%PathName=uigetdir(cd,'Select folder in  which to save mean and mode leave PD and Transverse sections data');%AG July 2013: for saving mean leaf data
resdir='Cross sections of characteristic leaf shape and of true mean 3D surface';
if ~exist([Pathname '\' resdir],'dir')
    mkdir(Pathname,resdir);
end
Pathname=[Pathname '\' resdir];

for k=1:2
    
    for di=min(days):max(days)
        %reshape into columns
        Xl=reshape(leaf_XIg{di},prod(size(leaf_YIg{di})),1);
        Yl=reshape(leaf_YIg{di},prod(size(leaf_YIg{di})),1);
        
        
        if k==1 %mean first
            varname='truemean';
            Zlmean=reshape(Truemean_leaf_ZIg{di},prod(size(Truemean_leaf_ZIg{di})),1);
            %Zlmean=reshape(Mode_leaf_ZIg{di},prod(size(Mode_leaf_ZIg{di})),1);
    
            XYZg_mean=[Xl Yl Zlmean];
            XYZg=XYZg_mean;
            
        else %then mode
            varname='characteristic';
            Zlmode=reshape(Mode_leaf_ZIg{di},prod(size(Mode_leaf_ZIg{di})),1);
            XYZg_mode=[Xl Yl Zlmode];
            XYZg=XYZg_mode;
         
        end
        
        %remove all NaN values (values from grid outside the leaf outline)
        XYZg(isnan(XYZg(:,3)),:)=[];
        
        %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;
        
        xl_all{di}=xl;
        yl_all{di}=yl;
        zl_all{di}=zl;
        xw_all{di}=xw;
        yw_all{di}=yw;
        zw_all{di}=zw;

        
    end
    save([Pathname '\cross_sections_' varname '.mat'],'xl_all','yl_all','zl_all','xw_all','yw_all','zw_all'); %save within each folder
    
end

disp(['Cross sections for characteristic leaf shape and true mean leaf have been saved to ' Pathname]);
