 function Show_cross_sections
%Program to display cross sections along the proximo-distal or medio-lateral axes 
%(user is prompted to choose which to display). Displays cross sections for all samples, 
%and superimpose cross sections for mean leaf shape and characteristic leaf shape. 
%Next program to run: Localcurvature

days=7:19;%days=input('Enter DAS numbers contained in this data set, e.g. [7:19] ');

close all
Path1=uigetdir(cd,'Select Reference data folder');%AG July 2013: for saving mean leaf data
result=[Path1 '\SourceData'];

% result=uigetdir(pwd,['Select the folder containing the RESULTS of ALL processed images']);
%
load([Path1 '\meanoutlines_DAS7-19.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
resdirall='Cross Sections of all samples (warped)';
resdirmean='Cross sections of characteristic leaf shape and of true mean 3D surface';

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

%
 opt=input('Press 1 for proximodistal or 2 for mediolateral cross sections: ');
 scale=input('Press 1 to scale all sections to the same 2D length for alignment purpose, 0 otherwise');
%alldisp=input('Press 1 to draw sections for all samples, 0 otherwise');
if opt==1,
    n1='PD';
else
    n1='ML';
end
if scale==1
    n2='scaled';
else
    n2=[];
end
resdirfigs=['Cross sections figures ' n1 n2];
if ~exist([Pathname '\' resdirfigs],'dir')
    mkdir(Pathname,resdirfigs);
end


for i=1:length(days)
    figure('Name',['DAS' sprintf('%02.0f',days(i))],'NumberTitle','off')
    title(['DAS' sprintf('%02.0f',days(i))])
    set(gcf,'WindowStyle','docked');
    rotate3d
    hold on
end

alldisp=1;
if alldisp==1
cd(result);
pot=dir;
npot=size(pot,1);
for i1=3:npot
    cd([result '\' pot(i1).name]);
    plt=dir;
    nplt=size(plt,1);
    for i2=3:nplt
        cd([result '\' pot(i1).name '\'  plt(i2).name]);
        das=dir;
        ndas=size(das,1);
        for i3=3:ndas
            cd([result '\' pot(i1).name '\' plt(i2).name '\' das(i3).name])
            if exist([Pathname '\' resdirall '\cross_sections_' pot(i1).name plt(i2).name das(i3).name '.mat']) 
                %load('allometry3D.mat');
                load([Pathname '\' resdirall '\cross_sections_' pot(i1).name plt(i2).name das(i3).name '.mat']);
                disp([pot(i1).name '\' plt(i2).name '\' das(i3).name])
                figure(str2num(das(i3).name(4:5))-days(1)+1)
                daynum=str2num(das(i3).name(4:5));
                col=[0.5 0.5 0.5];
                displaysection(avx,avy,daynum,yl,zl,xw,zw,col,scale,opt);
                
                
            end
        end
    end
end
end
%%

for k=1:2
    if k==1
        name='cross_sections_characteristic';
        col=[1 0 0];
    elseif k==2
        name='cross_sections_truemean';
        col=[0 0 1];
    end
    %[fn2, pn2]=uigetfile(pwd,['Load ' name ' leaf allometry file']);
    %load([pn2 '\' fn2]);
    load([Pathname '\' resdirmean '\' name]);
    for i=days(1):days(end)
        %figure(['DAS' num2str(i,'%02u')]);
        figure(i-days(1)+1)
        displaysection(avx,avy,i,yl_all{i},zl_all{i},xw_all{i},zw_all{i},col,scale,opt);
        if opt==1
            
            axis([-50 5500 -1000 1000]); %for PD display
            %
            
        else
            axis([-2750 2800 -1000 1000]) %for ML display
        end
        saveas(gcf,[Pathname '\' resdirfigs '\Figure' n1 n2 num2str(i)],'fig');
        
    end
end

disp(['Figures have been saved to ' Pathname '\' resdirfigs]);




function displaysection(avx,avy,daynum,yl,zl,xw,zw,col,scale,opt)
if opt==1
    if scale==1
    %                     yvals2=[min(avx{daynum}):(max(avx{daynum})-min(avx{daynum}))/(length(zl)-1):max(avx{daynum})]';
    yvals2=[0:(max(avy{daynum})-min(avy{daynum}))/(length(zl)-1):abs(min(avy{daynum}))+max(avy{daynum})]';
    
    maxyl=max(yl);
    minyl=min(yl);
    rangeyl=maxyl-minyl;
    rangeavy=max(avy{daynum})-min(avy{daynum});
    yvals3=yl*rangeavy/rangeyl;
    yvals4=(yl*rangeavy/rangeyl)-min(yl*rangeavy/rangeyl);
    else
        yvals4=yl-min(yl);%to show 0 at leaf base
    end
    plot(yvals4,zl,'DisplayName','xl,yl,zl','Color',col);
    hold on
    axis equal
    
elseif opt==2
    if scale==1
    %                     xvals=[min(avy{daynum}):(max(avy{daynum})-min(avy{daynum}))/(length(zw)-1):max(avy{daynum})]';
    xvals=[min(avx{daynum}):(max(avx{daynum})-min(avx{daynum}))/(length(zw)-1):max(avx{daynum})]';
    
    maxxw=max(xw);
    minxw=min(xw);
    rangexw=maxxw-minxw;
    rangeavx=max(avx{daynum})-min(avx{daynum});
    xvals3=xw*rangeavx/rangexw;
    xvals4=(xw*rangeavx/rangexw)-min(xw*rangeavx/rangexw)+min(avx{daynum});
    else
        xvals4=xw;
    end
    plot(xvals4,zw,'DisplayName','xw,yw,zw','Color',col);
    hold on
    axis equal
end


