function Display_deformations_and_growth_daytonext3D

%Display program to show grid deformation from one day to the next in 3D and/or to 
%show growth maps from 3D data in 2D or 3D – Used to make growth maps and to make 3D figures. 
%User is prompted to answer a few questions regarding display options.
%
%Next program to run: Display_deformations_forgridmovie3D
%
days=7:19;

close all
maingrid=5;
mainthick=1;%thickness of main grid lines
secondthick=0.1;
secondcol=[0.7 0.7 0.7];
scalegrid=0; % if scalegrid is set to 1, the pre-deformation shape is set to a constant size, so axes do not represent true size. A vertical scale bar is added to the figure;

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


opt1=input('press 0 to show grids only, 2 to show grid growth maps: ');
if opt1==2
    gridl=input('Show grid lines? y/n','s') ;
end
if opt1==0
    gridl='y';
end

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']); % we load this because predef and postdef mean outlines are stored there
    addname2='pol3';
else
    load([Pathname '\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([Pathname '\Grid deformations in 2D\Meangrid_day_by_day_deform07-19.mat']);
load([Pathname '\3D surfaces for 3D deformations\Mean3Dsurfaces_for_3D_deformations.mat']);

if opt1==2
    opt2=input('Press 2 to display results from 3D grid deformation in 2D, 3 to display in 3D');
else %if displaying grids, we display in 3D (2D display already done in previous steps of the analysis)
    opt2=3;
end

if opt2==2
    addname='_dispin2D';
else
    addname='_dispin3D';
end
if opt1==2 && gridl=='y'
    resdir='Surface deformations from day to day in 3D_withgrowth_withgridlines';
elseif opt1==2 && gridl=='n'
    resdir='Surface deformations from day to day in 3D_withgrowth_nogridlines';
elseif opt1==0
    resdir='Grid deformations from day to day in 3D';
end

resdir=[resdir addname addname2];

if ~exist([Pathname '\' resdir],'dir')
    mkdir(Pathname,resdir);
end
Pathname=[Pathname '\' resdir];

cd(Pathname)

fignum=0;

for t=days(1):days(end-1)
    figure; hold on
    set(gcf, 'color', 'white');
    axis equal;
    
    if scalegrid==1
        rangex=max(max(avxpredef{t}))-min(min(avxpredef{t}));
        scfac=1000/rangex;
    else scfac=1;
    end
    if opt2==3 % display mean 3D outine based on mode
        plot3([Mode_outlineZ{t,1}(:,1);Mode_outlineZ{t,1}(1,1)*scfac],[Mode_outlineZ{t,1}(:,2);Mode_outlineZ{t,1}(1,2)]*scfac,[Mode_outlineZ{t,1}(:,3);Mode_outlineZ{t,1}(1,3)]*scfac,'color',[0 0 0],'LineWidth',1)
    else
        
        
        plot([avxpredef{t};avxpredef{t}(1)]*scfac ,[avypredef{t};avypredef{t}(1)]*scfac,'color',[0 0 0],'LineWidth',mainthick)
        hold on
    end
    
    title(['DAS' num2str(t,'%02.0f') 'before deformation']);
    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');
    
    
    %NOTE: here can add grid colour-coded for growth
    %grid pre-deformation = gX{t}, gY{t}
    %grid post deformation = xm{t+1}, ym{t+1}
    
    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','none');
    gZ{t}=F(gX{t},gY{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);
    
    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','none');
    F2=scatteredInterpolant(pts(:,1),pts(:,2),pts(:,3),'linear','linear');%to extrapolate warped outline later
    
    zm{t+1}=F(xm{t+1},ym{t+1});
    [gridinterpX,gridinterpY,gridinterpZ]= growth_fromgrid_ag3D([avxpredef{t} avypredef{t}],gX{t},gY{t},gZ{t},xm{t+1},ym{t+1},zm{t+1},0);
    
    if opt2==2
        gZ{t}=gY{t}*0;
    end
    
    if opt1==0
        
        
        %%%%%%%5
        
        for k=1:size(gX{t},1) %drawing lines instead of mesh b/c mesh won't save to eps files
            if rem(k-1,maingrid)==0
                plot3(gX{t}(k,:)*scfac,gY{t}(k,:)*scfac,gZ{t}(k,:)*scfac,'k','LineWidth',mainthick);
            else
                hh=plot3(gX{t}(k,:)*scfac,gY{t}(k,:)*scfac,gZ{t}(k,:)*scfac,'LineWidth',secondthick);
                set(hh,'Color',secondcol);
            end
        end
        for j=1:size(gX{t},2)
            
            if rem(j-1,maingrid)==0
                plot3(gX{t}(:,j)*scfac,gY{t}(:,j)*scfac,gZ{t}(:,j)*scfac,'k','LineWidth',mainthick);
            else
                hh=plot3(gX{t}(:,j)*scfac,gY{t}(:,j)*scfac,gZ{t}(:,j)*scfac,'k','LineWidth',secondthick);
                set(hh,'Color',secondcol);
            end
        end
        %%%%%%%%
        
    elseif opt1==2 && gridl=='y'
        % h=pcolor(gridinterpX,gridinterpY,gridinterpZ);
        patchsurf3D(gX{t}*scfac,gY{t}*scfac,gZ{t}*scfac,gridinterpZ);
    elseif opt1==2 && gridl=='n'
        %else
        %  h=surfc(gX{t},gY{t},gY{t}*0,gridinterpZ,'LineStyle','none');
        
        %h=pcolor(gridinterpX,gridinterpY,gridinterpZ);
        %set(h,'LineStyle','none');
        patchsurf3D(gX{t}*scfac,gY{t}*scfac,gZ{t}*scfac,gridinterpZ,0);
    end
    %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 off;
    %AG: removed 'if' Aug 29
    %     if loaddata==3 && opt2==3
    
    %     end
    %if opt1~=0
    if scalegrid==1
        nscale=1000;
        title('click where you want the scale bar to be')
        [x1,y1]=ginput(1);
        hold on,
        h=line([x1;x1],[y1;y1+nscale*scfac]);
        set(h,'Color',[0 0 0]);
        
    end
    saveas(gcf,['Figure' num2str(fignum)],'fig');
    close gcf %if opt1==0, we will display the 2 grids side by side
    %end
    
    hold on
    if  opt2==3
        plot3([Mode_outlineZ{t+1,2}(:,1);Mode_outlineZ{t+1,2}(1,1)]*scfac,[Mode_outlineZ{t+1,2}(:,2);Mode_outlineZ{t+1,2}(1,2)]*scfac,[Mode_outlineZ{t+1,2}(:,3);Mode_outlineZ{t+1,2}(1,3)]*scfac,'color',[0 0 0],'LineWidth',1)
        olzw{t+1}=F2(olxw{t+1},olyw{t+1});
        %plot3([olxw{t+1};olxw{t+1}(1)]*scfac,[olyw{t+1};olyw{t+1}(1)]*scfac,[olzw{t+1};olzw{t+1}(1)]*scfac,'color', [0.7 0.7 1],'LineWidth',mainthick);
        plot3([olxw{t+1};olxw{t+1}(1)]*scfac,[olyw{t+1};olyw{t+1}(1)]*scfac,[olzw{t+1};olzw{t+1}(1)]*scfac,'color', [0 191/255 1],'LineWidth',mainthick);
    else
        cst=0;
        plot([avxpostdef{t+1};avxpostdef{t+1}(1)]*scfac+cst ,[avypostdef{t+1};avypostdef{t+1}(1)]*scfac,'color',[0 0 0],'LineWidth',mainthick);
        plot([olxw{t+1};olxw{t+1}(1)]*scfac+cst,[olyw{t+1};olyw{t+1}(1)]*scfac,'color', [0.7 0.7 0.7],'LineWidth',mainthick);
        
    end
    
    if opt2==2
        zm{t+1}=ym{t+1}*0;
    end
   
    if opt1==0
        for k=1:size(xm{t+1},1) %drawing lines instead of mesh b/c mesh won't save to eps files
            if rem(k-1,maingrid)==0
                h=plot3(xm{t+1}(k,:)*scfac,ym{t+1}(k,:)*scfac,zm{t+1}(k,:)*scfac,'k','LineWidth',mainthick);
            else
                hh=plot3(xm{t+1}(k,:)*scfac,ym{t+1}(k,:)*scfac,zm{t+1}(k,:)*scfac,'LineWidth',secondthick);
                set(hh,'Color',secondcol);
            end
            
        end
        for j=1:size(xm{t+1},2)
            if rem(j-1,maingrid)==0
                h=plot3(xm{t+1}(:,j)*scfac,ym{t+1}(:,j)*scfac,zm{t+1}(:,j)*scfac,'k','LineWidth',mainthick);
            else
                hh=plot3(xm{t+1}(:,j)*scfac,ym{t+1}(:,j)*scfac,zm{t+1}(:,j)*scfac,'LineWidth',secondthick);
                set(hh,'Color',secondcol);
            end
        end
    else
        patchsurf3D(xm{t+1}*scfac,ym{t+1}*scfac,zm{t+1}*scfac,gridinterpZ,2);%to only show deformed grid
    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 off;
    
    saveas(gcf,['Figure' num2str(fignum)],'fig');
    %        end
    close gcf;
    
    %keyboard
end

disp(['Figures of day to day deformations in 3D are saved in ' Pathname]);
