function [lineW,Tr]=warpleaves_ShapePaper
% Program used to warp (time t samples used to compute deformation patterns from t to t+1) to the average leaf shape 
% of those samples (Avxpredef{t},Avypredef{t}). It saves warped outlines and surfaces in warpedleaf3D_preddef.mat in each relevant sample’s 
% folder in the Source data or 'Results' folder of the Reference data
% folder [provided with the programs.
% This later makes it possible to compute the average (characteristic) leaf shape for each deformation grid at time t within a t_t+1 deformation pair.
% In the same way, this program also warps (time t samples used to compute deformation patterns from t-1 to t) to the average leaf shape 
% of those samples (Avxpostdef{t},Avypostdef{t}), and saves warped outline in warpedleaf3D_postdef.mat in each relevant sample’s folder in the source data folder of the 
% Reference dataset folder. 
% This later makes it possible to compute the average (characteristic) leaf shape for each deformation grid at time t+1 within a t_t+1 deformation pair.
%
% (NOTE: this warping is only done on outlines and leaf surfaces, not on
% beads as it is only used to align leaf surfaces in a subsequent program and obtain leaf deformations in 3D - at this point transformations have been computed already and beads 
% positions are not needed anymore - beads warping has been commented out of this warping program). 
%
%(Note: When running the program, the user is prompted to indicate whether Compute_deformation_patterns or Compute_deformation_patterns_lwm was used earlier).
%(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, warpleaves_ShapePaper does not have to be run twice).
%
% Next program to run: MeanLeaf_extract_fordeformation
%
optdisp=0; % this can be set to 1 to display a visual check of the warping.

daystart=7;
dayend=19;

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

PathName=uigetdir(cd,'Select main project folder in  which all shape results and figures are stored');
%pn0=uigetdir(pwd,'Folder where data for all pots are stored');
cd(pn0);
listd0=dir;
ndsize00=size(listd0,1);

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
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
end
%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

for prepost=1:2 % running the whole thing twice, once to make warpedleaf3D_predef.mat and once to make warpedleaf3D_postdef.mat for all relevant samples
if prepost==1
    rangedays=[daystart dayend-1];%predefshapes
else
    rangedays=[daystart+1 dayend];%postdef shapes
end
for day=rangedays(1):rangedays(2)%for v12
    % for ndays=1:numdays,
    %    day=days(ndays);
    dayproc=num2str(day,'%02.0f');
    if optdisp==1
        figure;
    end
    
    for npot=3:ndsize00,
        pn1=[pn0 '\' listd0(npot).name];
        cd(pn1);
        listd=dir;
        ndsize0=size(listd,1);
        
        for ndir=3:ndsize0,
            close all;
            if exist([pn1 '\' listd(ndir).name '\DAS' dayproc])
                cd([pn1 '\' listd(ndir).name '\DAS' dayproc]);
                %day= str2num(listd(ndir).name(4));
                disp(['Processing ' listd0(npot).name '\' listd(ndir).name '\DAS' dayproc]);
                listd1=dir;
                ndsize=size(listd1,1);
                figure;
                if prepost==1
                    condmet=exist('G_tonext_automatches.mat','file');
                    if condmet
                        load('G_tonext_automatches.mat');
                        if isempty(mat3D1)
                            condmet=0;
                        end
                        clear mat3D1 mat3D2
                    end
                else %prepost==2
                    condmet=exist([pn1 '\' listd(ndir).name '\DAS' num2str(day-1,'%02.0f') '\G_tonext_automatches.mat'],'file');%there was growth at the previous time point
                    if condmet
                        load([pn1 '\' listd(ndir).name '\DAS' num2str(day-1,'%02.0f') '\G_tonext_automatches.mat']);
                        if isempty(mat3D1)
                            condmet=0;
                        end
                        clear mat3D1 mat3D2
                    end
                end
                if condmet
                    
                    load([pwd '\Newbladeline3D.mat']);
                    load([pwd '\3Daligned_data.mat']);
                    xy=NewBladeLine;
                    
                    %Load relevant mean leaf outline
                    if prepost==1
                        avx{day}=avxpredef{day};
                        avy{day}=avypredef{day};
                    else %prepost==2
                        avx{day}=avxpostdef{day}(1:end-1);%because those mean outlines have been made to close
                        avy{day}=avypostdef{day}(1:end-1); % so last point is same as 1st point, but original leaf outline is not closed
                    end
                    
                    if optdisp==1
                        subplot(1,2,1);
                        %load([pwd '\growth3D.mat'],'tri');
                        
                        hold on, h=line([xy(:,1);xy(1,1)],[xy(:,2);xy(1,2)]); %this displays the leaf outline
                        axis equal
                    end;
                    stop=0;
                    iter=0;
                    g1=xy; % calculate the transformation based on leaf outline mapped to average outline
                    %g2=mat1; % before growth coordinates
                    g3=XYZg(:,1:2);
                    %g4=beadsg(:,3:4);
                    
                    while stop==0 && iter<10;%% index exceed matrix
                        Tr=get_2ndtrans(g1,[avx{day} avy{day}]);%% error
                        g1=trans2nd(Tr,[g1(:,1) g1(:,2)]);
                        %g2=trans2nd(Tr,[g2(:,1) g2(:,2)]);
                        g3=trans2nd(Tr,[g3(:,1) g3(:,2)]);
                        %g4=trans2nd(Tr,[g4(:,1) g4(:,2)]);
                        
                        iter=iter+1;
                        Tr_tot{iter}=Tr;
                        %Note: repeat until error minimized?
                        Averror=sum(sqrt((g1(:,1)-avx{day}).^2+(g1(:,2)-avy{day}).^2).^2);
                        if Averror<1000,
                            stop=1;
                        end;
                    end;
                    
                    lineW=g1;%outline
                    %mat1W=mat1;
                    %mat1W(:,1:2)=g2;
                    xyzgW=g3;
                    %beadsW=g4;
                    
                    if optdisp==1
                        subplot(1,2,2);
                        
                        hold on,h=line([g1(:,1);g1(1,1)],[g1(:,2);g1(1,2)]);set(h,'Color',[1 0 0]);
                        hold on,h=line([avx{day};avx{day}(1)],[avy{day}; avy{day}(1)]);set(h,'Color',[0 1 0]);
                        axis equal
                        
                    end
                    if prepost==1
                        save([pwd '\warpedleaf3D_predef.mat'],'lineW','xyzgW','Tr_tot');
                    else
                        save([pwd '\warpedleaf3D_postdef.mat'],'lineW','xyzgW','Tr_tot');
                    end
                    clear lineW xyzgW Tr_tot
                end
            end
        end;
    end;
end

end

