function [gridX1,gridY1,gridXW,gridYW,olw,olt2,olt1]=Getsingleaf_transfogrid(dirname,optfig)
%get transformation from single leaf basedon beads data. dirname is  DAS
%folder of the plant to track to next day
%program saving trsformations saved this:
% FileName=input(['File name: GridDAS' num2str(days(1),'%02.0f') '-'  num2str(days(end),'%02.0f') ], 's');
% save(['gridDASday_by_day' num2str(days(1),'%02.0f') '-'  num2str(days(end),'%02.0f') FileName],'xm','ym','olxw','olyw','NormError','gX','gY','Tr_rec',...
%     'ldiv_rec','Tr_n','wpt_rec','allbeads1');
%clear all
%close all

numregs=40;%number of regions along x or y for grid
nd=1; % number of sections for polynomial transformation - NOTE: will have to adapt so can have any trsfo type
%load trsfos data
%load('C:\Users\Anne-Gaelle\Dropbox\Lauren_shapepaper\July6WE\grids_deform_daybyday\gridDASday_by_day07-19_nofixes_withbeads.mat');


%dirname=uigetdir(pwd,'Select DAS folder of the plant you want to track to next day');
disp(dirname);
cur_day=str2num(dirname(1,end-1:end));
dirname2=[dirname(1,1:end-2) num2str(cur_day+1,'%02.0f')];

load([dirname '\3Daligned_data.mat']);
if optfig==1
    ha=figure;
    set(ha,'Name',['DAS' num2str(cur_day,'%02.0f') 'data']);
    h=line(leaflineg(:,1),leaflineg(:,2));
    hold on
    %beadsg: 1st 2 columns contain bead number and size, respectively.
    plot(beadsg(:,3),beadsg(:,4),'.');
end
%make and plot a grid
[gridX1, gridY1]=makegleafgrid(leaflineg(:,1),leaflineg(:,2),numregs,optfig);


%load('K:\My Documents\Manuscripts in prep\Lauren_shape_paper\AG reanalysis - copy of some of Lauren''s data\Copy of RESULTSgood\POT01\PLT1\DAS07\G_tonext_automatches.mat');
load([dirname '\G_tonext_automatches.mat']);

if optfig==1
    hold on,h=plot(mat3D1(:,1),mat3D1(:,2),'+');set(h,'Color',[1 0 0]);%shows beads used for matching- fits.
end

%b) show trsformed leaf based on mean trsfo
load([dirname '\Newbladeline3D.mat']);
%hb=figure;
olX=NewBladeLine(:,1);
olY=NewBladeLine(:,2);
t=7;%BC WORKING WITH DAS07
olt1=[olX,olY];

dimg=size(gridX1);
xl=reshape(gridX1,dimg(1)*dimg(2),1); %Make grid points into list
yl=reshape(gridY1,dimg(1)*dimg(2),1);

%get trsfo from beads for that leaf
%make and disp trsfo
[olw,wpt,xyW]=maketrsfo_by_section(mat3D1,mat3D2,nd,olX,olY,xl,yl);

gridXW=reshape(xyW(:,1),dimg);%*scx+(min(avx{t+1})-min(olw(:,1))); %Make back into matrix format for grid
gridYW=reshape(xyW(:,2),dimg);%*scy+(min(avy{t+1})-min(olw(:,2)));

if optfig==1
    hb=figure;
    set(hb,'Name',['DAS' num2str(cur_day,'%02.0f') 'trsformed to ' num2str(cur_day+1,'%02.0f') ' using polynomial trsfo']);
    h=line(olw(:,1),olw(:,2));
    hold on
    %beadsg: 1st 2 columns contain bead number and size, respectively.
    plot(wpt(:,1),wpt(:,2),'.');
    
    
    for k=1:size(gridXW,1) %drawing lines instead of mesh b/c mesh won't save to eps files
        plot(gridXW(k,:),gridYW(k,:),'k','LineWidth',.5);
    end
    for j=1:size(gridXW,2)
        plot(gridXW(:,j),gridYW(:,j),'k','LineWidth',.5);
    end
    
    
    %b2) show same as b but colour-coded from growth based on polynomial
    %drawing region by region
    
    
    ha2=figure;
    h=line(leaflineg(:,1),leaflineg(:,2));
    hold on
    for r=1:size(gridXW,1)-1
        for c=1:size(gridXW,2)-1
            xcoord1=[gridX1(r,c) gridX1(r,c+1) gridX1(r+1,c+1) gridX1(r+1,c)]';
            ycoord1=[gridY1(r,c) gridY1(r,c+1) gridY1(r+1,c+1) gridY1(r+1,c)]';
            area1=polyarea(xcoord1,ycoord1);
            gridinterpX(r,c)=mean(xcoord1);
            gridinterpY(r,c)=mean(ycoord1);
            
            xcoord=[gridXW(r,c) gridXW(r,c+1) gridXW(r+1,c+1) gridXW(r+1,c)]';
            ycoord=[gridYW(r,c) gridYW(r,c+1) gridYW(r+1,c+1) gridYW(r+1,c)]';
            area=polyarea(xcoord,ycoord);
            RG=100*(area-area1)/area1;
            h=patch(xcoord1,ycoord1,RG);
            [IN ON] = inpolygon(gridinterpX(r,c),gridinterpY(r,c),leaflineg(:,1),leaflineg(:,2));
            if IN+ON==0
                gridinterpZ(r,c)=NaN;
            else
                gridinterpZ(r,c)=RG;
            end;
            
            
        end
    end
    
    %showing growth data from transformation on leaf before deformation
    %first using grid
    ha3=figure;
    h=line(leaflineg(:,1),leaflineg(:,2));
    hold on
    h=surfc(gridX1,gridY1,gridY1*0,gridinterpZ,'LineStyle','none');
    
    %then using contourf and data at center of each grid square
    ha4=figure;
    h=line(leaflineg(:,1),leaflineg(:,2));
    hold on
    h=contourf(gridinterpX,gridinterpY,gridinterpZ,0:1:300,'LineStyle','none');
    
    
    %NOTE: we get data on grid on non warped leaf, we could technically get
    %trsfos from individual leaves and use that to get detailed growth info
    %across the leaf from all leaves. This could be used to make average growth
    %maps based on individual leaf transformations.
    
    caxis([0 300]);
    
    
end

olw=cat(1,olw,olw(1,:));

if optfig==1
    %d) overlay b and c
    hd=figure;
    set(hd,'Name',['DAS' num2str(cur_day,'%02.0f') 'trsformed to DAS' num2str(cur_day+1,'%02.0f') ' using polynomial trsfo overlayed with ' 'DAS' num2str(cur_day+1,'%02.0f') ' data']);
    
    hd2=line(olw(:,1),olw(:,2)); set(hd2,'Color',[1 0 0]);
    hold on
    %beadsg: 1st 2 columns contain bead number and size, respectively.
    hd2=plot(wpt(:,1),wpt(:,2),'.'); set(hd2,'Color',[1 0 0]);
    
end
%leaflineg1=leaflineg;

load([dirname2 '\Newbladeline3D.mat']);
olt2=NewBladeLine;
olt2=cat(1,olt2,olt2(1,:));

if optfig==1
    h=line(olt2(:,1),olt2(:,2)); set(h,'Color',[0 0 0]);
    hold on
    %beadsg: 1st 2 columns contain bead number and size, respectively.
    %h=plot(beadsg(:,3),beadsg(:,4),'.'); set(h,'Color',[0 0 0]);
    hold on,h=plot(mat3D2(:,1),mat3D2(:,2),'+');set(h,'Color',[0 0 0]);%shows beads matched from DAS7.


close(hd);
close(ha3);
close(ha4);
close(ha2);
close(hb);
close(ha);

end
function [olw,wpt,xyW]=trsfo_by_section(ldiv,nd,Tr_rec,olX,olY,mat1W,t,xl,yl);
    d=1;
    %Tr_n=nd;%number of trsfos to use - to record - AG 2013
    Tr=Tr_rec{d,t};
    %% Apply this to the grid points and deforming outline
    [olw(olY<ldiv(d,2),1:2)] = trans3rdpoly([olX(olY<ldiv(d,2)) olY(olY<ldiv(d,2))],Tr.tdata);
    [wpt(mat1W(:,2)<ldiv(d,2),1:2)] = trans3rdpoly(mat1W(mat1W(:,2)<ldiv(d,2),1:2),Tr.tdata);
    [xyW(yl<ldiv(d,2),1:2)] = trans3rdpoly([xl(yl<ldiv(d,2)) yl(yl<ldiv(d,2))],Tr.tdata); %Apply transformation

    
    for d=2:nd-1 % apply trsfo per section
        Tr=Tr_rec{d,t};
         [xyW(yl>=ldiv(d,1) & yl<ldiv(d,2),1:2)] = trans3rdpoly([xl(yl>=ldiv(d,1) & yl<ldiv(d,2)) yl(yl>=ldiv(d,1) & yl<ldiv(d,2))],Tr.tdata); %Apply transformation
        %[xyW(yl>=ldiv(d,1) & yl<ldiv(d,2),1:2)] = trans3rdpoly([xl(yl>=ldiv(d,1) & yl<ldiv(d,2)) yl(yl>=ldiv(d,1) & yl<ldiv(d,2))],Tr.tdata); %Apply transformation
        [olw(olY>=ldiv(d,1) & olY<ldiv(d,2),1:2)] = trans3rdpoly([olX(olY>=ldiv(d,1) & olY<ldiv(d,2)) olY(olY>=ldiv(d,1) & olY<ldiv(d,2))],Tr.tdata);
        [wpt(mat1W(:,2)>=ldiv(d,1) & mat1W(:,2)<ldiv(d,2),1:2)] = trans3rdpoly(mat1W(mat1W(:,2)>=ldiv(d,1) & mat1W(:,2)<ldiv(d,2),1:2),Tr.tdata);
    end
   
    d=nd;
    Tr=Tr_rec{d,t};
    %% Apply this to the grid points and deforming outline
    [xyW(yl>=ldiv(d,1),1:2)] = trans3rdpoly([xl(yl>=ldiv(d,1)) yl(yl>=ldiv(d,1))],Tr.tdata); %Apply transformation
    [olw(olY>=ldiv(d,1),1:2)] = trans3rdpoly([olX(olY>=ldiv(d,1)) olY(olY>=ldiv(d,1))],Tr.tdata);
    [wpt(mat1W(:,2)>=ldiv(d,1) ,1:2)] = trans3rdpoly(mat1W(mat1W(:,2)>=ldiv(d,1),1:2),Tr.tdata);
    



function [gridX, gridY]=makegleafgrid(listx,listy,num,optfig)
        d=(max(listy)-min(listy))-(max(listx)-min(listx));
        if d<0 % WAS if d>0 - AG changed June 2013
            %AG JUNE 2013: was /9, changed by /7 to match paper
            %sc=(max(avy{days(1)})-min(avy{days(1)}))/7; %find increment to divide leaf into 5 sections across the widest part
            sc=(max(listx)-min(listx))/num; %AG fix june 2013
            %[gridX gridY]=meshgrid(min(listx):sc:max(listx),min(avy{days(1)})+d/2+sc/2:sc:max(avy{days(1)})-d/2-sc/2); %make a grid
            [gridX gridY]=meshgrid(min(listx):sc:max(listx),min(listy)+d/2:sc:max(listy)-d/2); %AG fic June 2013
        else
            sc=(max(listy)-min(listy))/num; %find increment to divide leaf into 5 sections across the widest part
            %[gridX gridY]=meshgrid(min(listx)-d/2+sc/2:sc:max(listx)+d/2-sc/2,min(listy):sc:max(listy)); %make a grid
            [gridX gridY]=meshgrid(min(listx)-d/2:sc:max(listx)+d/2,min(listy):sc:max(listy)); %AG added to replace prev line june 2013
        end
        if optfig==1
            figure(gcf);
            for k=1:size(gridX,1) %drawing lines instead of mesh b/c mesh won't save to eps files
                plot(gridX(k,:),gridY(k,:),'k','LineWidth',.5);
            end
            for j=1:size(gridX,2)
                plot(gridX(:,j),gridY(:,j),'k','LineWidth',.5);
            end
        end;
function [olw,wpt,xyW]=maketrsfo_by_section(mat1W,mat2W,nd,olX,olY,xl,yl);

div=(max(yl)-min(yl))/nd;
ldiv=[min(yl) min(yl)+div];
nbeadspersection(1)=sum((mat1W(:,2)<ldiv(1,2)));
if nd>1
    for k=2:nd
        ldiv(k,1:2)=ldiv(k-1,:)+div;
        nbeadspersection(k)=sum((mat1W(:,2)>=ldiv(k,1))&(mat1W(:,2)<ldiv(k,2)));
    end
end
minbeadsn=min(nbeadspersection);

if minbeadsn>=10
    polorder=3;
elseif minbeadsn>=6
    polorder=2;
    disp('insufficient number of beads to estimate a 3rd order polynomial, using 2d order instead');
else
    if nd>1
        error('Insufficient number of beads to warp by section - try a single polynomial (nd=1)');
    else
        error('Insufficient number of beads to estimate a polynomial transformation');
    end
end

if nd==1
    
    Tr=cp2tform(mat2W(:,1:2), mat1W(:,1:2),'polynomial',polorder); %derive a 3rd order polynomial
    %% Apply this to the grid points and deforming outline
    
    if polorder==3
        [xyW] = trans3rdpoly([xl yl],Tr.tdata); %Apply transformation
        [olw] = trans3rdpoly([olX olY],Tr.tdata);
        [wpt] = trans3rdpoly(mat1W(:,1:2),Tr.tdata);
    else
        [xyW] = trans2nd(Tr.tdata,[xl yl]); %Apply transformation
        [olw] = trans2nd(Tr.tdata,[olX olY]);
        [wpt] = trans2nd(Tr.tdata,mat1W(:,1:2));
    end
else

    d=1;
    %Tr_n=nd;%number of trsfos to use - to record - AG 2013
    Tr=cp2tform(mat2W(mat1W(:,2)<ldiv(d,2),1:2), mat1W(mat1W(:,2)<ldiv(d,2),1:2),'polynomial',polorder); %derive a 3rd order polynomial

    if polorder==3
        %% Apply this to the grid points and deforming outline
        [olw(olY<ldiv(d,2),1:2)] = trans3rdpoly([olX(olY<ldiv(d,2)) olY(olY<ldiv(d,2))],Tr.tdata);
        [wpt(mat1W(:,2)<ldiv(d,2),1:2)] = trans3rdpoly(mat1W(mat1W(:,2)<ldiv(d,2),1:2),Tr.tdata);
        [xyW(yl<ldiv(d,2),1:2)] = trans3rdpoly([xl(yl<ldiv(d,2)) yl(yl<ldiv(d,2))],Tr.tdata); %Apply transformation
    else
        [olw(olY<ldiv(d,2),1:2)] = trans2nd(Tr.tdata,[olX(olY<ldiv(d,2)) olY(olY<ldiv(d,2))]);
        [wpt(mat1W(:,2)<ldiv(d,2),1:2)] = trans2nd(Tr.tdata,mat1W(mat1W(:,2)<ldiv(d,2),1:2));
        [xyW(yl<ldiv(d,2),1:2)] = trans2nd(Tr.tdata,[xl(yl<ldiv(d,2)) yl(yl<ldiv(d,2))]); %Apply transformation

    end

    for d=2:nd-1 % apply trsfo per section
        Tr=cp2tform(mat2W(mat1W(:,2)>=ldiv(d,1) & mat1W(:,2)<ldiv(d,2),1:2), mat1W(mat1W(:,2)>=ldiv(d,1) & mat1W(:,2)<ldiv(d,2),1:2),'polynomial',3); %derive a 3rd order polynomial
        
        if polorder==3
            [xyW(yl>=ldiv(d,1) & yl<ldiv(d,2),1:2)] = trans3rdpoly([xl(yl>=ldiv(d,1) & yl<ldiv(d,2)) yl(yl>=ldiv(d,1) & yl<ldiv(d,2))],Tr.tdata); %Apply transformation
            %[xyW(yl>=ldiv(d,1) & yl<ldiv(d,2),1:2)] = trans3rdpoly([xl(yl>=ldiv(d,1) & yl<ldiv(d,2)) yl(yl>=ldiv(d,1) & yl<ldiv(d,2))],Tr.tdata); %Apply transformation
            [olw(olY>=ldiv(d,1) & olY<ldiv(d,2),1:2)] = trans3rdpoly([olX(olY>=ldiv(d,1) & olY<ldiv(d,2)) olY(olY>=ldiv(d,1) & olY<ldiv(d,2))],Tr.tdata);
            [wpt(mat1W(:,2)>=ldiv(d,1) & mat1W(:,2)<ldiv(d,2),1:2)] = trans3rdpoly(mat1W(mat1W(:,2)>=ldiv(d,1) & mat1W(:,2)<ldiv(d,2),1:2),Tr.tdata);
        else
            [xyW(yl>=ldiv(d,1) & yl<ldiv(d,2),1:2)] = trans2nd(Tr.tdata,[xl(yl>=ldiv(d,1) & yl<ldiv(d,2)) yl(yl>=ldiv(d,1) & yl<ldiv(d,2))]); %Apply transformation
            %[xyW(yl>=ldiv(d,1) & yl<ldiv(d,2),1:2)] = trans3rdpoly([xl(yl>=ldiv(d,1) & yl<ldiv(d,2)) yl(yl>=ldiv(d,1) & yl<ldiv(d,2))],Tr.tdata); %Apply transformation
            [olw(olY>=ldiv(d,1) & olY<ldiv(d,2),1:2)] = trans2nd(Tr.tdata,[olX(olY>=ldiv(d,1) & olY<ldiv(d,2)) olY(olY>=ldiv(d,1) & olY<ldiv(d,2))]);
            [wpt(mat1W(:,2)>=ldiv(d,1) & mat1W(:,2)<ldiv(d,2),1:2)] = trans2nd(Tr.tdata,mat1W(mat1W(:,2)>=ldiv(d,1) & mat1W(:,2)<ldiv(d,2),1:2));
        end
    end

    d=nd;
    Tr=cp2tform(mat2W(mat1W(:,2)>=ldiv(d,1),1:2), mat1W(mat1W(:,2)>=ldiv(d,1),1:2),'polynomial',polorder); %derive a 3rd order polynomial
    
    if polorder==3
        %% Apply this to the grid points and deforming outline
        [xyW(yl>=ldiv(d,1),1:2)] = trans3rdpoly([xl(yl>=ldiv(d,1)) yl(yl>=ldiv(d,1))],Tr.tdata); %Apply transformation
        [olw(olY>=ldiv(d,1),1:2)] = trans3rdpoly([olX(olY>=ldiv(d,1)) olY(olY>=ldiv(d,1))],Tr.tdata);
        [wpt(mat1W(:,2)>=ldiv(d,1) ,1:2)] = trans3rdpoly(mat1W(mat1W(:,2)>=ldiv(d,1),1:2),Tr.tdata);
    else
        %% Apply this to the grid points and deforming outline
        [xyW(yl>=ldiv(d,1),1:2)] = trans2nd(Tr.tdata,[xl(yl>=ldiv(d,1)) yl(yl>=ldiv(d,1))]); %Apply transformation
        [olw(olY>=ldiv(d,1),1:2)] = trans2nd(Tr.tdata,[olX(olY>=ldiv(d,1)) olY(olY>=ldiv(d,1))]);
        [wpt(mat1W(:,2)>=ldiv(d,1) ,1:2)] = trans2nd(Tr.tdata,mat1W(mat1W(:,2)>=ldiv(d,1),1:2));
    end

end