function MeanLeaf_extract
%Program to compute a true mean and a characteristic leaf shape for each time point from all samples available at that time point.
%Uses samples warped to the mean leaf shape computed from all samples available at each time point.
%From the mean and characteristic leaf  fit shapes we can then calculate cross sections along proximodistal or mediolateral 
%axis and compare those to cross sections of all samples to see how mean and best fit leaf shapes fit the dataset. 
%We can also compute spatial patterns of local curvature. All those are computed in the next programs.
%Next program to run : Sample_cross_sections.


days=7:19;
close all;
Path1=uigetdir(cd,'Select Reference data folder');
load([Path1 '\meanoutlines_DAS7-19.mat']);
Pathdata=[Path1 '\SourceData'];

Pathname=uigetdir(cd,'Select (or make and select) main project folder in  which all shape results and figures will be stored, e.g. name it Shape Results');%AG July 2013: for saving mean leaf data

surfdir='3D surface and figures';
if ~exist([Pathname '\' surfdir],'dir')
    mkdir(Pathname,surfdir);
end
Pathname=[Pathname '\' surfdir];
fname=fullfile(Pathname,'Mean3Dsurface.mat');

%Pathdata='S:\Camille Girard-Bock\Shape paper_Summer2013\RESULTSgood';
%Savedata=['S:\Camille Girard-Bock\Shape paper_Summer2013\figures\Camille\Mode_Leaf'];
leaf_XIg={};
leaf_YIg={};
leaf_ZIg={};
Mode_leaf_ZIg={};
%numberDAS=zeros(13,1);
numberDAS={};

for i=min(days):max(days)
    kk=0;
    %step=10;
    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}=XIg;
    leaf_YIg{i}=YIg;
    leaf_ZIg{i}=zeros(size(XIg,1),size(XIg,2));
    lineW3dat={};
    Mode_leaf_ZIg{i}=zeros(size(XIg,1),size(XIg,2));
    numberDAS{i}=zeros(size(XIg,1),size(XIg,2));
    [IN ON] = inpolygon(leaf_XIg{i},leaf_YIg{i},avx{i},avy{i});
    KP=IN+ON;
    KPnot=~KP;
    leaf_ZIg{i}(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 exist([Pathdata '\' pot(x).name '\' plt(y).name '\DAS' num2str(i,'%02d') ])
                    cd([Pathdata '\' pot(x).name '\' plt(y).name '\DAS' num2str(i,'%02d') ])
                    load('warpedleaf3D.mat'); %load xyzgW %note AG: loads lineW too, has 99 pts
                    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)];
                        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}(ii,jj)=numberDAS{i}(ii,jj)+1;
                                leaf_ZIg{i}(ii,jj,kk)=ZIg(ii,jj);
                            end
                        end
                        lineW3dat{i}(:,:,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},3)
                temps_ZIg(kk)=leaf_ZIg{i}(ii,jj,kk);
            end
            if sum(isnan(temps_ZIg))<=numel(temps_ZIg)/2
                B=[1:numberDAS{i}(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}(ii,jj)=1;
                else
                    Mode_leaf_ZIg{i}(ii,jj)=0; %all values are at 0
                end
            else
                Mode_leaf_ZIg{i}(ii,jj)=NaN; %if more than the half of the values at this point are NaN
            end
            %             [IN ON] = inpolygon(leaf_XIg{i}(ii,jj),leaf_YIg{i}(ii,jj),avx{i},avy{i});
            %             if IN+ON==0
            %                 Mode_leaf_ZIg{i}(ii,jj)=NaN;
            %             end
        end
    end
    
    Mode_outlineZ{i}=mean(lineW3dat{i}(:,:,index),3);
    mean_xy_allsamples=mean(lineW3dat{i}(:,1:2,:),3);
    Mode_outlineZ{i}=cat(2,mean_xy_allsamples,Mode_outlineZ{i}(:,3));
    Mean_outlineZ{i}=mean(lineW3dat{i}(:,:,:),3);
    
     Truemean_leaf_ZIg{i}=zeros(size(Mode_leaf_ZIg{i}));
    RMSE_truemeanleaf{i}=zeros(size(Mode_leaf_ZIg{i}));
    RMSE_modeleaf{i}=zeros(size(Mode_leaf_ZIg{i}));
    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}(ii,jj)) || Mode_leaf_ZIg{i}(ii,jj)==0;
                Truemean_leaf_ZIg{i}(ii,jj)=mean(leaf_ZIg{i}(ii,jj,:));
                RMSE_truemeanleaf{i}(ii,jj)=sum(sqrt((leaf_ZIg{i}(ii,jj,:) - Truemean_leaf_ZIg{i}(ii,jj)).^2))/size(leaf_ZIg{i},3);
                RMSE_modeleaf{i}(ii,jj)=sum(sqrt((leaf_ZIg{i}(ii,jj,:) - Mode_leaf_ZIg{i}(ii,jj)).^2))/size(leaf_ZIg{i},3);
                
            else
                Mode_leaf_ZIg{i}(ii,jj)=mean(leaf_ZIg{i}(ii,jj,index));%AG: See Camille: index is overwritten all the time in the previous loop?
                
                Truemean_leaf_ZIg{i}(ii,jj)=mean(leaf_ZIg{i}(ii,jj,:));
                RMSE_truemeanleaf{i}(ii,jj)=sum(sqrt((leaf_ZIg{i}(ii,jj,:) - Truemean_leaf_ZIg{i}(ii,jj)).^2))/size(leaf_ZIg{i},3);
                RMSE_modeleaf{i}(ii,jj)=sum(sqrt((leaf_ZIg{i}(ii,jj,:) - Mode_leaf_ZIg{i}(ii,jj)).^2))/size(leaf_ZIg{i},3);
                
                
            end
            [IN ON] = inpolygon(leaf_XIg{i}(ii,jj),leaf_YIg{i}(ii,jj),avx{i},avy{i});
            if IN+ON==0
                Mode_leaf_ZIg{i}(ii,jj)=NaN;
                Truemean_leaf_ZIg{i}(ii,jj)=NaN;
                RMSE_truemeanleaf{i}(ii,jj)=NaN;
                RMSE_modeleaf{i}(ii,jj)=NaN;
            end
        end
    end
    
    tmp=reshape(RMSE_truemeanleaf{i},prod(size(RMSE_truemeanleaf{i})),1);
    tmp2=tmp((~isnan(tmp)));
    RMSE_truemeanleaftot(i)=sum(tmp2);
    
    tmp=reshape(RMSE_modeleaf{i},prod(size(RMSE_modeleaf{i})),1);
    tmp2=tmp((~isnan(tmp)));
    RMSE_modeleaftot(i)=sum(tmp2);

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

%AG July 2013: saving mean leaf data
%FileName=input(['File name: MeanleavesDAS7_19'], 's');
save(fname,'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(['Characteristic (and true mean) 3D surface data and figures have been saved to ' Pathname]);

%cd(['S:\Camille Girard-Bock\Shape paper_Summer2013\programs'])






