function Compute_deformation_patterns
% Program to compute mean deformation patterns, based on samples for which we have bead displacements from 1 day to the next. It also saves the mean outlines of samples used to compute the deformation: saves mean outlines predeformation and post deformation.
% Avxpredef{t} and Avypredef{t}: mean of samples at time t used to compute deformation from t to t+1
% Avxpostdef{t} and Avypredef{t}:mean of samples at time t used when computing deformation from t-1 to t.
%
% optfig=0: no figures displayed, optfig=1: figures are displayed (as the program runs) that show deformations of single samples as they are computed
%look at detailed code and comments therein to get more information on the
%figures that are displayed. Those are useful for checking what individual
%sample deformations are - program will take longer to run.
%when optfig is set to 1, single deformation figures are displayed but not
%saved (as there will be a very large number of figures), moreover figures
%are closed after each single transformation has been calculated, so needs to be at the compute while the program runs to visualize those.
%
%Next program: Display_Deformations_DaytoNext

if nargin<1
optfig=0; 
end

days=7:19;


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

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

%PathName=uigetdir(cd,'Select folder in  which to save mean grids'' data and figures');
resdir='Grid deformations in 2D';
if ~exist([PathName '\' resdir],'dir')
    mkdir(PathName,resdir);
end
PathName=[PathName '\' resdir];


cd(allresults);
pots=dir;

fignum=1;

%% Go through time points 
for t=days(1):days(end-1), %for each time point that user wants to process
    gridX1M=[];
    gridY1M=[];
    gridXWM=[];
    gridYWM=[];
    olWM=[];
    olt2m=[];
    olt1m=[];
    col=jet(20);
    testct=0;
    for npot=3:length(pots), %for each pot
        cd([allresults '\' pots(npot).name]);
        plants=dir;
        for nplt=3:length(plants), %for each plant
            pltdir=[allresults '\' pots(npot).name '\' plants(nplt).name];
            if exist([pltdir '\DAS' num2str(t,'%02.0f') '\G_tonext_automatches.mat'], 'file') %if it has data for the time point in question
                load([pltdir '\DAS' num2str(t,'%02.0f') '\G_tonext_automatches.mat']);
                if ~isempty(mat3D1)
                    clear mat3D1 mat3D2
                    disp(['Processing DAS' num2str(t,'%02.0f') ': ' pots(npot).name ' - ' plants(nplt).name]);
                    dirname=[pltdir '\DAS' num2str(t,'%02.0f')];
                    [gridX1,gridY1,gridXW,gridYW,olW,olt2,olt1]=Getsingleaf_transfogrid(dirname,optfig);%olw = warped outline
                    gridX1M=cat(3,gridX1M,gridX1);
                    gridY1M=cat(3,gridY1M,gridY1);
                    gridXWM=cat(3,gridXWM,gridXW);
                    gridYWM=cat(3,gridYWM,gridYW);
                    
                    olt2m=cat(3,olt2m,olt2);
                    olt1m=cat(3,olt1m,olt1);
                    olWM=cat(3,olWM,olW);
                    
                    % %                     %test grids - used for methods figure to explain grid
                    %                    % averaging
                    %                     figure;
                    %                     testct=testct+1;
                    %                     hold on,
                    %                     for k=1:10:size(gridX1,1) %drawing lines instead of mesh b/c mesh won't save to eps files
                    %                         plot(gridX1(k,:),gridY1(k,:),'LineWidth',.5,'Color',col(testct,:));
                    %                     end
                    %                     for j=1:10:size(gridX1,2)
                    %                         plot(gridX1(:,j),gridY1(:,j),'LineWidth',.5,'Color',col(testct,:));
                    %                     end
                    %                     olt1=cat(1,olt1,olt1(1,:));
                    %                     h=line(olt1(:,1),olt1(:,2)); set(h,'Color',col(testct,:));
                    %                      saveas(gcf,[PathName '\figpre' num2str(testct)],'fig');
                    %                     figure
                    %                     hold on,
                    %                     for k=1:10:size(gridX1,1) %drawing lines instead of mesh b/c mesh won't save to eps files
                    %                         plot(gridXW(k,:),gridYW(k,:),'LineWidth',.5,'Color',col(testct,:));
                    %                     end
                    %                     for j=1:10:size(gridXW,2)
                    %                         plot(gridXW(:,j),gridYW(:,j),'LineWidth',.5,'Color',col(testct,:));
                    %                     end
                    %                     h=line(olW(:,1),olW(:,2)); set(h,'Color',col(testct,:));
                    %                     %h=line(olt2(:,1),olt2(:,2)); set(h,'Color',col(testct,:));
                    %                     saveas(gcf,[PathName '\figpost' num2str(testct)],'fig');
                    %
                    %                     disp('')
                end
            end
        end
    end
    
    
    gridX1M=mean(gridX1M,3);
    gridY1M=mean(gridY1M,3);
    gridXWM=mean(gridXWM,3);
    gridYWM=mean(gridYWM,3);
    olt2m2=mean(olt2m,3);
    olWM2=mean(olWM,3);
    olt1m2=mean(olt1m,3);
    
    avxpredef{t}=olt1m2(:,1); %calculating mean shapes (from samples used) before deformation
    avypredef{t}=olt1m2(:,2);
    avxpostdef{t+1}=olt2m2(:,1); %recalculating mean shapes (from samples used) after deformation
    avypostdef{t+1}=olt2m2(:,2);
    
    if optfig==1
        fig1=figure;
        hold on
        for k=1:10:size(gridX1M,1) %drawing lines instead of mesh b/c mesh won't save to eps files
            plot(gridX1M(k,:),gridY1M(k,:),'k','LineWidth',.5);
        end
        for j=1:10:size(gridX1M,2)
            plot(gridX1M(:,j),gridY1M(:,j),'k','LineWidth',.5);
        end
        %line([avx{t};avx{t}(1)],[avy{t}; avy{t}(1)]);
        h=line(olt1m2(:,1),olt1m2(:,2));set(h,'Color',[1 0 0]);

        line([avxpredef{t};avxpredef{t}(1)],[avypredef{t}; avypredef{t}(1)]);
        
        fig2=figure;
        hold on
        for k=1:10:size(gridXWM,1) %drawing lines instead of mesh b/c mesh won't save to eps files
            plot(gridXWM(k,:),gridYWM(k,:),'k','LineWidth',.5);
        end
        for j=1:10:size(gridXWM,2)
            plot(gridXWM(:,j),gridYWM(:,j),'k','LineWidth',.5);
        end
        %line([avx{t+1};avx{t+1}(1)],[avy{t+1}; avy{t+1}(1)]);
        line([avxpostdef{t+1};avxpostdef{t+1}(1)],[avypostdef{t+1}; avypostdef{t+1}(1)]);
        hh=line(olWM2(:,1),olWM2(:,2));set(hh,'Color',[1 0 0]);
        
    end
    
    %store grid and mean of tranformed outliers in cell array
    gX{t}=gridX1M;
    gY{t}=gridY1M;
    xm{t+1}=gridXWM;
    ym{t+1}=gridYWM;
    olxw{t+1}=olWM2(:,1); % olxw and olyw: mean of single leaf warps of outlines
    olyw{t+1}=olWM2(:,2);
    
    %now calculating tranformation based on mean grids
    %first reshape points data
    gridXr=reshape(gridX1M,prod(size(gridX1M)),1);
    gridYr=reshape(gridY1M,prod(size(gridY1M)),1);
    gridXWr=reshape(gridXWM,prod(size(gridXWM)),1);
    gridYWr=reshape(gridYWM,prod(size(gridYWM)),1);
    
    %get transformation and apply it
    Tr_rec={};
    polorder=3;%using 3rd order polynomial - could use other method (2 = 2d order, 4 = lwm)
    nd=1;
    [olw,wpt,xyW,Tr_rec,ldiv]=maketrsfo_by_section_ag2([gridXr gridYr],[gridXWr gridYWr],nd,avxpredef{t},avypredef{t},gridXr,gridYr,Tr_rec,t,polorder);
    %Note: wpt = warped grid, xyW = warped grid also - rendundant bc program used
    %for other things also.
    
    %Save grid and outline warped from t to t+1 based on calculated transformation
    xmW{t+1}=reshape(xyW(:,1),size(gridX1M));
    ymW{t+1}=reshape(xyW(:,2),size(gridX1M));
    olxW2{t+1}=olw(:,1);
    olyW2{t+1}=olw(:,2);
    
    olw=cat(1,olw,olw(1,:));
    if optfig==1
    hold on, h=line(olw(:,1),olw(:,2));
    set(h,'Color',[0 1 0]);
    
    %     hold on
    %     for k=1:size(xmW{t+1},1) %drawing lines instead of mesh b/c mesh won't save to eps files
    %         plot(xmW{t+1}(k,:),ymW{t+1}(k,:),'r','LineWidth',.5);
    %     end
    %     for j=1:size(xmW{t+1},2)
    %         plot(xmW{t+1}(:,j),ymW{t+1}(:,j),'r','LineWidth',.5);
    %     end
    axis equal
    axis([min(min(xmW{t+1}))-40 max(max(xmW{t+1}))+40 min(min(ymW{t+1}))-40 max(max(ymW{t+1}))+40]);
    fignum=fignum+1;
    %saveas(gcf,[PathName '\Figure' num2str(fignum) 'DAS' num2str(t+1)],'fig');
    
    
    figure(fig1);
    axis equal
    axis([min(min(xmW{t+1}))-40 max(max(xmW{t+1}))+40 min(min(ymW{t+1}))-40 max(max(ymW{t+1}))+40]);
    %saveas(gcf,[PathName '\Figure' num2str(fignum-1) 'DAS' num2str(t)],'fig');
    
    end
    close all
    
end

%PathName=uigetdir(cd,'Select folder in  which to save mean grids');
cd(PathName)
save(['Meangrid_day_by_day_deform' num2str(days(1),'%02.0f') '-'  num2str(days(end),'%02.0f')],'xm','ym','gX','gY',...
    'avxpredef','avypredef', 'avxpostdef','avypostdef','olxw','olyw','xmW','ymW','olxW2','olyW2');

disp(['Deformation patterns in 2D have been saved in ' PathName]);

end

