function MeanLeaf_extract_fordeformation
%Program used to compute characteristic and mean leaf shapes based on sets of leaves pre and post deformation, respectively. 
%This saves the true average and characteristic leaf shape for each deformation grid at time t within a t_t+1 deformation pair, 
%and the true average and characteristic leaf shape for each deformation grid at time t+1 within a t_t+1 deformation pair.
%Results are stored in cell array format such that modeleaf{day,1} = leaf at day t within a t,t+1 deformation pair, 
%modeleaf{day,2}=leaf at time t within a t-1,t deformation pair.
%(Note both true mean and characteristic leaf shapes are saved but only the characteristic leaf shape in used in subsequent programs).
%
%The data can then be used to show, for each deformation pair, the mean leaf grid in 3D before and after deformation. 
%We can then compute, for each square of the grid, its growth from time t to time t+1 using the deformation of the square in 3D space. 
%We can then compare the resulting growth patterns to the growth patterns computed from 3D particles positions in Remmler and Rolland-Lagan, 2012.
%
%(Note: user is prompted to indicate whether Compute_deformation_patterns or Compute_deformation_patterns_lwm was used earlier –
%note that this is only needed because those programs both save the average outlines from samples pre and post deformation, 
%and the matrices in which those are saved have different names depending on which program was used. However the outlines are the same 
%with both methods, so if one wants to test the 2 types of transformations, MeanLeaf_extract_fordeformation does not have to be run twice).
%
%Next program to run: Display_deformations_and_growth_daytonext3D

days=7:19;

Path1=uigetdir(cd,'Select Reference data folder');
Pathdata=[Path1 '\SourceData'];

Pathname=uigetdir(cd,'Select main project folder in  which all shape results and figures are stored');

trsfo=input('If you used Compute_deformation_patterns Press 1, if you used Compute_deformation_patterns_lwm press 2');
if trsfo==1
    load([Pathname '\Grid deformations in 2D\Meangrid_day_by_day_deform07-19.mat']); 
else
    load([Pathname '\Grid deformations in 2D\Meangrid_day_by_day_deform_lwm07-19.mat']); 
end

%load([Pathname '\Grid deformations in 2D\Meangrid_day_by_day_deform07-19.mat']);

%Pathname=uigetdir(cd,'Select folder in  which to save mean leaves data and figures');%AG July 2013: for saving mean leaf data
surfdir='3D surfaces for 3D deformations';
if ~exist([Pathname '\' surfdir],'dir')
    mkdir(Pathname,surfdir);
end
Pathname=[Pathname '\' surfdir];

leaf_XIg={};
leaf_YIg={};
leaf_ZIg={};
Mode_leaf_ZIg={};
numberDAS={};

for prepost=1:2

if prepost==1
    rangedays=[days(1) days(end-1)];%predefshapes
else
    rangedays=[days(2) days(end)];%postdef shapes
end
for i=rangedays(1):rangedays(2)%
    kk=0;
    %step=10;
    
    if prepost==1 % shape predeformation at time t
        avx{i}=avxpredef{i};
        avy{i}=avypredef{i};
    else %prepost==2 % shape postdeformation (shapes at time t resulting from growth from t-1)
        avx{i}=avxpostdef{i};
        avy{i}=avypostdef{i};
    end
    step=(abs(floor(min(avx{i})))+ceil(max(avx{i})))/40;
    ix=floor(min(avx{i})):step:ceil(max(avx{i}));
    iy=floor(min(avy{i})):step:ceil(max(avy{i}));
    [XIg,YIg] = meshgrid(ix,iy);
    leaf_XIg{i,prepost}=XIg;
    leaf_YIg{i,prepost}=YIg;
    leaf_ZIg{i,prepost}=zeros(size(XIg,1),size(XIg,2));
    lineW3dat={};
    Mode_leaf_ZIg{i,prepost}=zeros(size(XIg,1),size(XIg,2));
    numberDAS{i,prepost}=zeros(size(XIg,1),size(XIg,2));
    [IN ON] = inpolygon(leaf_XIg{i,prepost},leaf_YIg{i,prepost},avx{i},avy{i});
    KP=IN+ON;
    KPnot=~KP;
    leaf_ZIg{i,prepost}(KPnot)=NaN;
    cd (Pathdata);
    pot=dir;
    npots=size(pot,1);
    for x=3:npots
        cd([Pathdata '\' pot(x).name]);
        plt=dir;
        nplt=size(plt,1);
        %curpot=str2num(pot(i).name(1,size(pot(i).name,2)-1:size(pot(i).name,2)));
        for y=3:nplt
            cd([Pathdata '\' pot(x).name '\' plt(y).name]);
            das=dir;
            ndas=size(das,1);
            %for z=3:ndas
                
                if prepost==1
                    condmet=exist([Pathdata '\' pot(x).name '\' plt(y).name '\DAS' num2str(i,'%02d') '\G_tonext_automatches.mat'], 'file');%making means of leaf samples use to compute growth predeformation
                    if condmet
                        load([Pathdata '\' pot(x).name '\' plt(y).name '\DAS' num2str(i,'%02d') '\G_tonext_automatches.mat']);
                        if isempty(mat3D1)
                            condmet=0;
                        end
                        clear mat3D1 mat3D2
                    end
                else %prepost=2
                    condmet=exist([Pathdata '\' pot(x).name '\' plt(y).name '\DAS' num2str(i-1,'%02d') '\G_tonext_automatches.mat'], 'file');%if the current sample had growth data stored on previous day, it is to be considered for mean post deformation shape
                    if condmet
                        load([Pathdata '\' pot(x).name '\' plt(y).name '\DAS' num2str(i-1,'%02d') '\G_tonext_automatches.mat']);
                        if isempty(mat3D1)
                            condmet=0;
                        end
                        clear mat3D1 mat3D2
                    end
                end
                if condmet 
                    cd([Pathdata '\' pot(x).name '\' plt(y).name '\DAS' num2str(i,'%02d') ])
                    if prepost==1
                        load('warpedleaf3D_predef.mat'); %load leaf at time t that will be used for growth from t to t+1, leaf warped to avx and avypredef{t}
                        %load('warpedleaf3D.mat');%if it is the last day, can load normal wraping - not from samples that have growth data
                    else % prepost==2
                        load('warpedleaf3D_postdef.mat'); %load leaf at time t that was used for growth from t-1 to t, leaf warped to avx and avypostdef{t}
                        %load leaf warped to  xyzgW %note AG: loads lineW too, has 99 pts %data warped to mean shape for samples that have growth data
                    end
                    load('3Daligned_data.mat');XYZ=XYZg;%NOTEAG: loads leaflineg too
                     load('Newbladeline3D.mat');%AG - Load NewBladeLine - use to recover 3D outline: this only has 2D outline!
                    %AG July 2013: 
                    %lineW and NewBladeLine have corresponding (2D) points
                    %leaflineg matches 2D shap0e of NewBladeLine but has
                    %many more points and has z coordinates.
                    %For each points in NewBladeLine, find closest (x,y) in
                    %leaflineg and get the z. This z can be assigned to
                    %lineW.
                    if size(NewBladeLine,2)==2
                        NewBladeLine=cat(2,NewBladeLine,zeros(size(NewBladeLine,1),1));%add column for z
                        for ii=1:length(NewBladeLine)
                            dist=distance(leaflineg(:,1:2),NewBladeLine(ii,1:2));
                            [val,ind]=min(dist);
                            NewBladeLine(ii,3)=leaflineg(ind,3);
                        end
                    end
                    %recovered 3D NewBladeLine (save it? See Camille and decide), now get warped 3D line:
                    lineW3=cat(2,lineW,NewBladeLine(:,3));
                    
% AG: This was used to double check results of recovering 3D outline and
% warped 3D outline.
%                     figure;
%                     h=line(NewBladeLine(:,1),NewBladeLine(:,2),NewBladeLine(:,3));
%                     hold on
%                     h=line(leaflineg(:,1),leaflineg(:,2),leaflineg(:,3));set(h,'Color',[0 0 0]);
%                     h=line(lineW3(:,1),lineW3(:,2),lineW3(:,3));set(h,'Color',[1 0 0]);
%                     axis equal
                    
                    %%Camille - TO DO: make mean warped outline based on most used
                    %%leaves as for surface.

                    %DAS=(das(z).name);
                    %DAS(:,1:3)=[];
                    %DAS=str2num(DAS);
                    %if DAS==i
                        kk=kk+1;
                        %numberDAS(i)=numberDAS(i)+1;
                        
                        pts=[xyzgW(:,1),xyzgW(:,2),XYZ(:,3)];%removed for
                        %v10
                        %pts=[[xyzgW(:,1);lineW3(:,1)],[xyzgW(:,2);lineW3(:,2)],[XYZ(:,3);lineW3(:,3)]];%used
                        %for v10
                        pts=unique(pts,'rows');
                        indx=find(~isnan(pts(:,3)));
                        pts=pts(indx,:);
                        %F=scatteredInterpolant(pts(:,1),pts(:,2),pts(:,3),'linear','linear');
                        F=scatteredInterpolant(pts(:,1),pts(:,2),pts(:,3),'linear','none');
                        %F=TriScatteredInterp(pts(:,1),pts(:,2),pts(:,3)); % changed for V6
                        ZIg=F(XIg,YIg);
                        
                        [IN ON] = inpolygon(XIg,YIg,avx{i},avy{i});
                        INON=IN+ON;
                        ZIg(INON==0)=NaN;
%                         INOnly=IN-ON;  %CHANGED for v6
%                         ZIg(INOnly==0)=NaN;
                        
                        for ii=1:size(ZIg,1) % go through grid rows
                            for jj=1:size(ZIg,2); % go through grid columns
                                numberDAS{i,prepost}(ii,jj)=numberDAS{i,prepost}(ii,jj)+1;
                                leaf_ZIg{i,prepost}(ii,jj,kk)=ZIg(ii,jj);
                            end
                        end
                        lineW3dat{i,prepost}(:,:,kk)=lineW3;%storing outline data
                    %end
                end
            %end
        end
    end
    nbclass=5;
    temps_ZIg=[];
    most_freq=[];
    for ii=1:size(ZIg,1) % go through grid rows
        for jj=1:size(ZIg,2); % go through grid columns
            for kk=1:size(leaf_ZIg{i,prepost},3)
                temps_ZIg(kk)=leaf_ZIg{i,prepost}(ii,jj,kk);
            end
            if sum(isnan(temps_ZIg))<=numel(temps_ZIg)/2
                B=[1:numberDAS{i,prepost}(ii,jj)]'; %indices
                %     if max(abs(temps_ZIg))~=sum(abs(temps_ZIg))
                %%AG - sep27: removed - should not need this: 
                %%cd('S:\Camille Girard-Bock\Shape paper_Summer2013\programs\Lauren')
                meanval=nanmean(temps_ZIg,2); %if the mean is 0 we register 0 for that point without calculating the mode
                border=(max(temps_ZIg)-min(temps_ZIg))/100;
                if border~=0
                    %border=abs(meanval/100); %calculate what we're going to add to both min and max in order to creates categories that contain all values
                    
                    maxM=max(temps_ZIg)+border;
                    minM=min(temps_ZIg)-border;
                    interv=(maxM-minM)/nbclass; %create interval between edges in order to get 'nbclass' classes
                    edges=minM:interv:maxM; %creates edges
                    for ct=1:nbclass-1 % ct is the index corresponding to the category we're in
                        indx=find(temps_ZIg>=edges(ct) & temps_ZIg<edges(ct+1)); %index of numbers in that category
                        list_ind{ct}=indx; % we register all the indez in list_ind
                        freq(ct)=length(indx); % freq is the matrix containing how many numbers are in different categories
                        classval(ct)=mean(edges(ct:ct+1)); %class val is the mean value for each category (not used in this program)
                    end
                    [valmax,indmode]=max(freq); %the numver of values int he most frequent category (valmax) and the index of the category
                    most_freq=cat(1,most_freq,B(list_ind{indmode})); %list of all the index of those in most frequent category
                    n=histc(most_freq,1:max(most_freq)); %calculates how many times each index number from 1:max(most_freq) appear in most_freq
                    nfreq=cat(2,n,[1:max(most_freq)]'); %in 'nfreq' at the i row 1st column is the number of time i apears in most_freq, the second column show the value of i
                    index=find(n>=(max(n)/2)); %index of most used plant in the calculation %AG: See Camille: at that position only?
                    Mode_leaf_ZIg{i,prepost}(ii,jj)=1;
                else
                    Mode_leaf_ZIg{i,prepost}(ii,jj)=0; %all values are at 0
                end
            else
                Mode_leaf_ZIg{i,prepost}(ii,jj)=NaN; %if more than the half of the values at this point are NaN
            end
            %             [IN ON] = inpolygon(leaf_XIg{i,prepost}(ii,jj),leaf_YIg{i,prepost}(ii,jj),avx{i},avy{i});
            %             if IN+ON==0
            %                 Mode_leaf_ZIg{i,prepost}(ii,jj)=NaN;
            %             end
        end
    end
    
    Mode_outlineZ{i,prepost}=mean(lineW3dat{i,prepost}(:,:,index),3);
    mean_xy_allsamples=mean(lineW3dat{i,prepost}(:,1:2,:),3);
    Mode_outlineZ{i,prepost}=cat(2,mean_xy_allsamples,Mode_outlineZ{i,prepost}(:,3));
    Mean_outlineZ{i,prepost}=mean(lineW3dat{i,prepost}(:,:,:),3);
    
     Truemean_leaf_ZIg{i,prepost}=zeros(size(Mode_leaf_ZIg{i,prepost}));
    RMSE_truemeanleaf{i,prepost}=zeros(size(Mode_leaf_ZIg{i,prepost}));
    RMSE_modeleaf{i,prepost}=zeros(size(Mode_leaf_ZIg{i,prepost}));
    for ii=1:size(ZIg,1) % go through grid rows
        for jj=1:size(ZIg,2); % go through grid columns
            if isnan(Mode_leaf_ZIg{i,prepost}(ii,jj)) || Mode_leaf_ZIg{i,prepost}(ii,jj)==0;
                Truemean_leaf_ZIg{i,prepost}(ii,jj)=mean(leaf_ZIg{i,prepost}(ii,jj,:));
                RMSE_truemeanleaf{i,prepost}(ii,jj)=sum(sqrt((leaf_ZIg{i,prepost}(ii,jj,:) - Truemean_leaf_ZIg{i,prepost}(ii,jj)).^2))/size(leaf_ZIg{i,prepost},3);
                RMSE_modeleaf{i,prepost}(ii,jj)=sum(sqrt((leaf_ZIg{i,prepost}(ii,jj,:) - Mode_leaf_ZIg{i,prepost}(ii,jj)).^2))/size(leaf_ZIg{i,prepost},3);
                
            else
                Mode_leaf_ZIg{i,prepost}(ii,jj)=mean(leaf_ZIg{i,prepost}(ii,jj,index));%AG: See Camille: index is overwritten all the time in the previous loop?
                
                Truemean_leaf_ZIg{i,prepost}(ii,jj)=mean(leaf_ZIg{i,prepost}(ii,jj,:));
                RMSE_truemeanleaf{i,prepost}(ii,jj)=sum(sqrt((leaf_ZIg{i,prepost}(ii,jj,:) - Truemean_leaf_ZIg{i,prepost}(ii,jj)).^2))/size(leaf_ZIg{i,prepost},3);
                RMSE_modeleaf{i,prepost}(ii,jj)=sum(sqrt((leaf_ZIg{i,prepost}(ii,jj,:) - Mode_leaf_ZIg{i,prepost}(ii,jj)).^2))/size(leaf_ZIg{i,prepost},3);
                
                
            end
            [IN ON] = inpolygon(leaf_XIg{i,prepost}(ii,jj),leaf_YIg{i,prepost}(ii,jj),avx{i},avy{i});
            if IN+ON==0
                Mode_leaf_ZIg{i,prepost}(ii,jj)=NaN;
                Truemean_leaf_ZIg{i,prepost}(ii,jj)=NaN;
                RMSE_truemeanleaf{i,prepost}(ii,jj)=NaN;
                RMSE_modeleaf{i,prepost}(ii,jj)=NaN;
            end
        end
    end
    
    tmp=reshape(RMSE_truemeanleaf{i,prepost},prod(size(RMSE_truemeanleaf{i,prepost})),1);
    tmp2=tmp((~isnan(tmp)));
    RMSE_truemeanleaftot(i)=sum(tmp2);
    
    tmp=reshape(RMSE_modeleaf{i,prepost},prod(size(RMSE_modeleaf{i,prepost})),1);
    tmp2=tmp((~isnan(tmp)));
    RMSE_modeleaftot(i)=sum(tmp2);

%     Truemean_leaf_ZIg{i,prepost}=zeros(size(Mode_leaf_ZIg{i,prepost}));
%     RMSE_truemeanleaf{i,prepost}=zeros(size(Mode_leaf_ZIg{i,prepost}));
%     RMSE_modeleaf{i,prepost}=zeros(size(Mode_leaf_ZIg{i,prepost}));
%     for ii=1:size(ZIg,1) % go through grid rows
%         for jj=1:size(ZIg,2); % go through grid columns
%             if isnan(Mode_leaf_ZIg{i,prepost}(ii,jj)) || Mode_leaf_ZIg{i,prepost}(ii,jj)==0;
%                 
%             else
%                 Mode_leaf_ZIg{i,prepost}(ii,jj)=mean(leaf_ZIg{i,prepost}(ii,jj,index));%AG: See Camille: index is overwritten all the time in the previous loop?
%                 
%                 Truemean_leaf_ZIg{i,prepost}(ii,jj)=mean(leaf_ZIg{i,prepost}(ii,jj,:));
%                 RMSE_truemeanleaf{i,prepost}(ii,jj)=sum(sqrt(leaf_ZIg{i,prepost}(ii,jj,:) - Truemean_leaf_ZIg{i,prepost}(ii,jj)).^2)/size(leaf_ZIg{i,prepost},3);
%                 RMSE_modeleaf{i,prepost}(ii,jj)=sum(sqrt(leaf_ZIg{i,prepost}(ii,jj,:) - Mode_leaf_ZIg{i,prepost}(ii,jj)).^2)/size(leaf_ZIg{i,prepost},3);
%                 [IN ON] = inpolygon(leaf_XIg{i,prepost}(ii,jj),leaf_YIg{i,prepost}(ii,jj),avx{i},avy{i});
%                 if IN+ON==0
%                     Mode_leaf_ZIg{i,prepost}(ii,jj)=NaN;
%                     Truemean_leaf_ZIg{i,prepost}(ii,jj)=NaN;
%                     RMSE_truemeanleaf{i,prepost}(ii,jj)=NaN;
%                     RMSE_modeleaf{i,prepost}(ii,jj)=NaN;
%                 end
%             end
%         end
%     end
    
    
    figure;
    surfc(leaf_XIg{i,prepost},leaf_YIg{i,prepost},Mode_leaf_ZIg{i,prepost},Mode_leaf_ZIg{i,prepost}); 
    %h=line(avx{i},avy{i}); %for 2D outline
    h=line([Mode_outlineZ{i,prepost}(:,1); Mode_outlineZ{i,prepost}(1,1)],[Mode_outlineZ{i,prepost}(:,2); Mode_outlineZ{i,prepost}(1,2)],[Mode_outlineZ{i,prepost}(:,3);Mode_outlineZ{i,prepost}(1,3)]);
    set(h,'Color',[1 0 0]); 
    
    shading interp;
    axis equal;
    view(0,90);
    if prepost==1
        namefig='predef';
    else %prepost=2
        namefig='postdef';
    end
    title(['Mean leaf for DAS' num2str(i) ' ' namefig]);
        cd(Pathname);
        saveas(gcf, ['CharacteristicLeaf_height_DAS' num2str(i) namefig 'adjusted_grid_40.fig']);
        saveas(gcf, ['CharacteristicLeaf_height_DAS' num2str(i) namefig 'adjusted_grid_40.png']);
end

%AG July 2013: saving mean leaf data
cd(Pathname)

end
%FileName=input(['File name: MeanleavesDAS7_19'], 's');
save([Pathname '\Mean3Dsurfaces_for_3D_deformations.mat'],'leaf_XIg','leaf_YIg','Mode_leaf_ZIg','Mode_outlineZ','Mean_outlineZ','RMSE_modeleaf','RMSE_truemeanleaf','Truemean_leaf_ZIg','RMSE_modeleaftot','RMSE_truemeanleaftot');
%FileName= 'K:\My Documents\Manuscripts in prep\Lauren_shape_paper\AG reanalysis - copy of some of Lauren''s data\Meanleaves_freqmode.mat';
%save(FileName,'leaf_XIg','leaf_YIg','Mode_leaf_ZIg');

disp(['3D surfaces to use for 3D deformations have been saved in ' Pathname '\Mean3Dsurfaces_for_3D_deformations.mat']);
%cd(['S:\Camille Girard-Bock\Shape paper_Summer2013\programs'])






