function Deform_across_days
%To display grid deformations across time points without resetting the grid (the grid from the first time point 
%is deformed through all subsequent time points).
%Figures named Figdef0.fig, Figdef1.fig etc… show a display in 2D of the deformation across time points.
%The green outline shows how the original outline (of the 1st time point) deforms through subsequent time points. 
%It is normal that over time this outline shifts away from the outline of the leaf at each time point for 2 reasons: 
%1)small errors accumulate from time point to time point (i.e. in the day to day deformations the deformed outline is very 
%close to the actual outline of the next time point but it is not exactly the same); 3) the deformed outline originate from 
%the deformation of the sample outlines from the 1st time point. The samples at a later time point n may not be derived from the 
%same plants (unless the exact same set of leaves can be followed through all time points). Moreover, part of the leaf outline at 
%the 1st time point may be missing (e.g. at DAS 7 the very bottom of the blade may not have been extracted, this ‘missing part’ will 
%still be missing in the deformed outline across time points).
%Figures named Figure1.fig, Figure2.fig etc… will be displayed in 2D or 3D depending on the user’s choice. Figures with odd numbers 
%are colour-coded for growth and show relative growth to the next time point (relative growth is defined as in Remmler and Rolland-Lagan, 2012). Figures with even numbers show the grid post deformation and are not coloured but show the position of the deformed outline.
%Choices given to the user include: 
%- the choice of days to display deformations across (the more time points are included, the more errors will accumulate).
%- whether to display deformations in 2D or 3D (this will only apply to Figure1.fig, Figure2.fig… etc…).
%- if display is in 2D, whether to display all grid lines or a subset of grid lines on the odd-numbered figures (this is to make it easier to visualize growth data - all grid lines are displayed in the even-numbered figures as those are not coloured for growth).
%-if display is in 3D, whether to display grid lines or no grid lines on the odd-numbered figures.

close all
maingrid=5;
mainthick=1;%thickness of main grid lines
secondaries=0;%no display of small grid lines - if =1, will display them

%redef_regs=input('Press 1 to redefine regions to track, 0 otherwise');
redef_regs=0; 
trsfo=input('If you used Compute_deformation_patterns Press 1, if you used Compute_deformation_patterns_lwm press 2');


%polorder=input('Calculations options: Press 2 to use 2d order polynomials, 3 for 3rd order, 4 for lwm transformations');
%polorder=input('Calculations options for new computations: Press 3 for 3rd order polynomials, 4 for lwm transformations'); 
if trsfo==1
    polorder=3;
elseif trsfo==2
    polorder=4;
end
% note 2 is a valid option for polorder and will use 2d order polynomials, but is not suitable for our data

%loaddata=input('Press 2 to display 2D maps (if map data exist), 3 to display 3D maps (if maps and 3D mean leaf exist)');
loaddata=3;
opt1=2;


Path1=uigetdir(cd,'Select main project folder in  which all shape results and figures are stored');
if trsfo==1
    load([Path1 '\Grid deformations in 2D\Meangrid_day_by_day_deform07-19.mat']); % we load this because predef and postdef mean outlines are stored there
    addname2='pol3';
else
    load([Path1 '\Grid deformations in 2D\Meangrid_day_by_day_deform_lwm07-19.mat']); % we load this because predef and postdef mean outlines are stored there
    addname2='lwm';
end
%load([Path1 '\Grid deformations in 2D\Meangrid_day_by_day_deform07-19.mat']);

%[fn4 pn4]=uigetfile(pwd,['Load data to be displayed','*.mat']);%load avxpredef and avxpostdef
%load([pn4 '\' fn4]);
%PathName=uigetdir(pn4,'Select folder in  which to save mean grids');

%% Enter days to be processed and file locations
days=input('Enter DAS numbers to be processed e.g. [7:12],[11:19],or [7:19]');

resdir=['Figures of deformations across days' num2str(days(1)) '-' num2str(days(end))];
if ~exist([Path1 '\' resdir addname2],'dir')
    mkdir(Path1,[resdir addname2]);
end
PathName=[Path1 '\' resdir addname2];


cd(PathName)


if loaddata==3
    load([Path1 '\3D surfaces for 3D deformations\Mean3Dsurfaces_for_3D_deformations.mat']);
    opt2=input('Press 2 to display results from 3D grid deformation in 2D, 3 to display in 3D');
    if opt2==2
        gridl=input('Show grid lines of all regions? y/n (pressing n will only show a subset of grid lines)','s') ;
    elseif opt2==3
        gridl=input('Show grid lines ? y/n','s') ;
        
    end
    
end

fignum=0;


for t=days(1):days(end-1)
    figure; hold on
    set(gcf, 'color', 'white');
    axis equal;

    %if it is the first time point: load the grid and display it
    if t==days(1)
        
        cur_olx{t}=[avxpredef{t};avxpredef{t}(1)];%current outline
        cur_oly{t}=[avypredef{t};avypredef{t}(1)];
        
        if redef_regs==1
            numregs=input('How many regions along a side of the grid?');
            %make and plot a grid
            [gXcm{t}, gYcm{t}]=makeleafgrid_ag(cur_olx{t}, cur_oly{t},numregs,0);
        else
            %cm= cumulative deformations
            gXcm{t}=gX{t};
            gYcm{t}=gY{t};
            
        end
    
    end
    
    plot(cur_olx{t},cur_oly{t},'Color',[0 1 0],'LineWidth',1);
    hold on
    %plot(avx{t},avy{t},'Color',[0 0 1],'LineWidth',1);
            plot([avxpredef{t};avxpredef{t}(1)],[avypredef{t};avypredef{t}(1)],'Color',[0 0 1],'LineWidth',1);

    drawgrid(gXcm{t},gYcm{t},mainthick,maingrid,secondaries);
    
    axis equal
    axis([-3500 3500 -3500 3500]);
    saveas(gcf,['Figuredef' num2str(fignum)],'fig');
    fignum=fignum+1;
    %then: deform the initial grid based on transformations
    %calculated from day to day grid deformations
    %
    
    %now calculating tranformation based on mean grids
    %first reshape points data
    gridXr=reshape(gX{t},prod(size(gX{t})),1);
    gridYr=reshape(gY{t},prod(size(gY{t})),1);
    gridXWr=reshape(xm{t+1},prod(size(xm{t+1})),1);
    gridYWr=reshape(ym{t+1},prod(size(ym{t+1})),1);
    
    gridXrc=reshape(gXcm{t},prod(size(gXcm{t})),1);
    gridYrc=reshape(gYcm{t},prod(size(gYcm{t})),1);
    
    Tr_rec={};
    nd=1;
    npts=12;
   
    [olw,wpt,xyW,Tr_rec,ldiv]=maketrsfo_by_section_ag3([gridXr gridYr],[gridXWr gridYWr],nd,cur_olx{t},cur_oly{t},gridXrc,gridYrc,Tr_rec,t,polorder,npts);
   
    gXcm{t+1}=reshape(xyW(:,1),size(gXcm{t}));
    gYcm{t+1}=reshape(xyW(:,2),size(gXcm{t}));
    cur_olx{t+1}=olw(:,1);
    cur_oly{t+1}=olw(:,2);
    
    
    if t==days(end-1)
        figure
        plot(cur_olx{t+1},cur_oly{t+1},'Color',[0 1 0],'LineWidth',1);
        hold on
        plot([avxpostdef{t+1};avxpostdef{t+1}(1)],[avypostdef{t+1};avypostdef{t+1}(1)],'Color',[0 0 1],'LineWidth',1);
            drawgrid(gXcm{t+1},gYcm{t+1},mainthick,maingrid,secondaries);
        axis equal
        axis([-4000 4000 -4000 4000]);
    colorbar
    set(gcf, 'color', 'white');

        saveas(gcf,['Figuredef' num2str(fignum)],'fig');
        fignum=fignum+1;
    end
    
end
        
%%Now drawing figs with growth info
close all;
fignum=0;

for t=days(1):days(end-1)
    figure; hold on
    set(gcf, 'color', 'white');
    axis equal; 
    
    %     axis([-5000 5000 -5000 5000]) %axis([min(min(gYcm{days(end)})) max(max(gYcm{days(end)})) min(min(gYcm{days(end)})) max(max(gYcm{days(end)}))])
    if loaddata==3 && opt2==3 % display mean 3D outine based on mode
        plot3([Mode_outlineZ{t,1}(:,1);Mode_outlineZ{t,1}(1,1)],[Mode_outlineZ{t,1}(:,2);Mode_outlineZ{t,1}(1,2)],[Mode_outlineZ{t,1}(:,3);Mode_outlineZ{t,1}(1,3)],'color',[0 0 0],'LineWidth',1)
    else
        plot([avxpredef{t};avxpredef{t}(1)] ,[avypredef{t};avypredef{t}(1)],'color',[0 0 0],'LineWidth',1)
    end
    if opt1==0 || opt1==1
            drawgrid(gXcm{t},gYcm{t},mainthick,maingrid,secondaries);
    end

    title(['DAS' num2str(t,'%02.0f') 'before deformation']);
    axis([-3000 3000 -3000 3000]);
    if t<=9
        caxis([0 300]);
    elseif t<=11
        caxis([0 180]);
    elseif t<=13
        caxis([0 80]);
    elseif t<=15
        caxis([0 40]);
    else
        caxis([0 20]);
    end
    %set(gcf,'WindowStyle','Docked');
    colorbar
    set(gcf, 'color', 'white');

    
    %NOTE: here can add grid colour-coded for growth
    %grid pre-deformation = gXcm{t}, gYcm{t}
    %grid post deformation = xm{t+1}, ym{t+1}
    if opt1==2 || opt1==3
        if loaddata==3
            
            nrows=prod(size(leaf_XIg{t,1}));
            leafX=reshape(leaf_XIg{t,1},nrows,1);
            leafY=reshape(leaf_YIg{t,1},nrows,1);
            leafZ=reshape(Mode_leaf_ZIg{t,1},nrows,1);
            pts=[[leafX;Mode_outlineZ{t,1}(:,1)],[leafY; Mode_outlineZ{t,1}(:,2)],[leafZ;Mode_outlineZ{t,1}(:,3)]];
            pts=unique(pts,'rows');
            indx=find(~isnan(pts(:,3)));
            pts=pts(indx,:);
            F=scatteredInterpolant(pts(:,1),pts(:,2),pts(:,3),'linear','linear');
            gZcm{t}=F(gXcm{t},gYcm{t});

            
             %%
            
            nrows=prod(size(leaf_XIg{t+1,2}));
            leafX=reshape(leaf_XIg{t+1,2},nrows,1);
            leafY=reshape(leaf_YIg{t+1,2},nrows,1);
            leafZ=reshape(Mode_leaf_ZIg{t+1,2},nrows,1);
            %F=TriScatteredInterp(leafX,leafY,leafZ);
            %F=TriScatteredInterp([leafX;Mode_outlineZ{t+1}(:,1)],[leafY; Mode_outlineZ{t+1}(:,2)],[leafZ;Mode_outlineZ{t+1}(:,3)]);
           
            pts=[[leafX;Mode_outlineZ{t+1,2}(:,1)],[leafY; Mode_outlineZ{t+1,2}(:,2)],[leafZ;Mode_outlineZ{t+1,2}(:,3)]];
            pts=unique(pts,'rows');
            indx=find(~isnan(pts(:,3)));
            pts=pts(indx,:);
            F=scatteredInterpolant(pts(:,1),pts(:,2),pts(:,3),'linear','linear');
            
            gZcm{t+1}=F(xm{t+1},ym{t+1});
            [gridinterpX,gridinterpY,gridinterpZ]= growth_fromgrid_ag3D([avxpredef{t} avypredef{t}],gXcm{t},gYcm{t},gZcm{t},gXcm{t+1},gYcm{t+1},gZcm{t+1},0);
        else
            [gridinterpX,gridinterpY,gridinterpZ]= growth_fromgrid_ag([avxpredef{t} avypredef{t}],gXcm{t},gYcm{t},gXcm{t+1},gYcm{t+1},0);
        end
    end
    if opt1==2
        if gridl=='y'
            % h=surfc(gXcm{t},gYcm{t},gYcm{t}*0,gridinterpZ);
            if loaddata==3 && opt2==3
               % h=surfc(gXcm{t},gYcm{t},gZcm{t},gridinterpZ);
               patchsurf3D(gXcm{t},gYcm{t},gZcm{t},gridinterpZ); % note could also plot in 2D based on 3D grid data
            else
               % h=pcolor(gridinterpX,gridinterpY,gridinterpZ);
               patchsurf3D(gXcm{t},gYcm{t},0,gridinterpZ);
                   drawgrid(gXcm{t},gYcm{t},mainthick,maingrid,secondaries);

            end
        else
            %  h=surfc(gXcm{t},gYcm{t},gYcm{t}*0,gridinterpZ,'LineStyle','none');
            if loaddata==3 && opt2==3
                %h=surfc(gXcm{t},gYcm{t},gZcm{t},gridinterpZ,'LineStyle','none');
                patchsurf3D(gXcm{t},gYcm{t},gZcm{t},gridinterpZ,0);
            else
                %h=pcolor(gridinterpX,gridinterpY,gridinterpZ);
                %set(h,'LineStyle','none');
                patchsurf3D(gXcm{t},gYcm{t},0,gridinterpZ,0);
                drawgrid(gXcm{t},gYcm{t},mainthick,maingrid,secondaries);
            end
        end
    elseif opt1==3
        h=contourf(gridinterpX,gridinterpY,gridinterpZ,0:1:300,'LineStyle','none');
    end
    
    
    fignum=fignum+1;
%     print('-depsc2',['Grid DAS' num2str(days(1)) '-' num2str(days(end)) '   predefDAS' num2str(t,'%02.0f')]);
%     saveas(gcf,['Grid DAS' num2str(days(1)) '-' num2str(days(end)) '   predefDAS' num2str(t,'%02.0f')],'png');
    
      axis equal
        axis([-4000 4000 -4000 4000]);
    colorbar
    set(gcf, 'color', 'white');
%AG: removed 'if' Aug 29
%     if loaddata==3 && opt2==3
         saveas(gcf,['Figure' num2str(fignum)],'fig');
%     end
    close gcf
    
    if opt1==3 || opt1==3
        %set(gcf,'WindowStyle','Docked');
    else
        figure; hold on
        if loaddata==3 && opt2==3
            leafX=cat(1,leafX,[Mode_outlineZ{t+1,2}(:,1);Mode_outlineZ{t+1,2}(1,1)]);
            leafY=cat(1,leafY,[Mode_outlineZ{t+1,2}(:,2);Mode_outlineZ{t+1,2}(1,2)]);
            leafZ=cat(1,leafZ,[Mode_outlineZ{t+1,2}(:,3);Mode_outlineZ{t+1,2}(1,3)]);
            plot3([Mode_outlineZ{t+1,2}(:,1);Mode_outlineZ{t+1,2}(1,1)],[Mode_outlineZ{t+1,2}(:,2);Mode_outlineZ{t+1,2}(1,2)],[Mode_outlineZ{t+1,2}(:,3);Mode_outlineZ{t+1,2}(1,3)],'color',[0 0 0],'LineWidth',1)
            
            %to get 3D warped outline:
            %make 3D delaunay triangulation of grid points in 2D
            listpts=unique([leafX leafY leafZ],'rows');
            indx=find(~isnan(listpts(:,3)));
            listpts=listpts(indx,:);
            dtri=DelaunayTri(listpts(:,1:2));
            PI = nearestNeighbor(dtri,[cur_olx{t+1},cur_oly{t+1}]);%indx to closest point for each query point
            %this plots the triangle colosest to the ptn th point
            %hold on,h=plot(listpts(PI(ptn),1),listpts(PI(ptn),2),'*');set(h,'Color',[1 0 0])
            olzw3{t+1}=listpts(PI,3);
            
            xy=[[cur_olx{t+1};olxw{t+1}(1)],[cur_oly{t+1};olyw{t+1}(1)],[olzw3{t+1};olzw3{t+1}(1)]];            
            tt=1:length(xy);
            tts = 1: 0.1: length(xy);
            %xys = spline(t,xy,ts);
            
            xs = spline(tt,xy(:,1),tts);
            ys = spline(tt,xy(:,2),tts);
            zs = spline(tt,xy(:,3),tts);
            %a = 1;
            windowSize=50;
            zs2=cat(2,zs,zs,zs);
            zsnew=filter(ones(1,windowSize)/windowSize,1,zs2);
            zsnew=zsnew(length(zs)+1:2*length(zs));
%             b = [1/wndw 1/wndw 1/wndw 1/wndw];
%             zsnew=filter(b,a,zs);
            h=line(xs,ys,zsnew);set(h,'Color',[0 1 0]);%smoothed using averaging window
            
            
        else
            plot([avxpostdef{t+1};avxpostdef{t+1}(1)] ,[avypostdef{t+1};avypostdef{t+1}(1)],'color',[0 0 0],'LineWidth',1);
            plot([cur_olx{t+1};cur_olx{t+1}(1)],[cur_oly{t+1};cur_oly{t+1}(1)],'color', [0 1 0],'LineWidth',0.5);
        end
        
        
        if loaddata==3
            patchsurf3D(gXcm{t+1},gYcm{t+1},gZcm{t+1},gridinterpZ,2);%to only show deformed grid
        else
                drawgrid(gXcm{t+1},gYcm{t+1},mainthick,maingrid,secondaries);

        end
        
       
        title(['DAS' num2str(t+1,'%02.0f') 'deformed from DAS' num2str(t,'%02.0f')]);
        axis([-3000 3000 -3000 3000]);
    if t<=8
        caxis([0 300]);
    elseif t<=10
        caxis([0 180]);
    elseif t<=12
        caxis([0 80]);
    elseif t<=14
        caxis([0 40]);
    else
        caxis([0 20]);
    end
    set(gcf,'WindowStyle','Docked');
    colorbar
    set(gcf, 'color', 'white');
    axis equal; 
    
    
        fignum=fignum+1;
        %axis([-3000 3000 -3000 3000]);
%         print('-depsc2',['Grid DAS' num2str(days(1)) '-' num2str(days(end)) '   postdefDAS' num2str(t+1,'%02.0f')]);
%         saveas(gcf,['Grid DAS' num2str(days(1)) '-' num2str(days(end)) '   postdefDAS' num2str(t,'%02.0f')],'png');
        axis off;
% AG Aug 29 - removed 'if'
%        if loaddata==3 && opt2==3
        saveas(gcf,['Figure' num2str(fignum)],'fig');
%        end
        close gcf;
    end
    %keyboard
end

disp(['Figures of deformations across days have been saved in ' PathName]);

function drawgrid(gX,gY,mainthick,maingrid,second)

if nargin<5
    second=1;
end
for k=1:size(gX,1) %drawing lines instead of mesh b/c mesh won't save to eps files
    if rem(k-1,maingrid)==0
        plot(gX(k,:),gY(k,:),'k','LineWidth',mainthick);
    else
        if second==1
        plot(gX(k,:),gY(k,:),'k','LineWidth',.5);
        end
    end
end
for j=1:size(gX,2)
    
    if rem(j-1,maingrid)==0
        plot(gX(:,j),gY(:,j),'k','LineWidth',mainthick);
    else
        if second==1
        plot(gX(:,j),gY(:,j),'k','LineWidth',.5);
        end
    end
end


