function [olw,wpt,xyW,Tr_rec,ldiv]=maketrsfo_by_section_ag3(mat1W,mat2W,nd,olX,olY,xl,yl,Tr_rec,t,polorder,numpoints);

if nargin<11
    numpoints=round(size(mat1W,1)*50/100);
end
if nargin<10
polorder=3;
end
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 
    if nargin<10
        polorder=3;
    end
elseif minbeadsn>=6 
    polorder=2;
    disp('insufficient number of beads to estimate a 3rd order polynomial, or user choice: using 2d order polynomial');
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 or lwm transformation');
    end
end

if nd==1 || polorder==4
    
    if polorder<=3
    Tr=cp2tform(mat2W(:,1:2), mat1W(:,1:2),'polynomial',polorder); %derive a 3rd order polynomial
    elseif polorder==4 % apply lwm trsfo
        Tr=cp2tform(mat2W(:,1:2), mat1W(:,1:2),'lwm',numpoints); 
    end
    Tr_rec{1,t}=Tr;%To record trsfos - AG 2013
    %% 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);
    elseif polorder==2
        [xyW] = trans2nd(Tr.tdata,[xl yl]); %Apply transformation
        [olw] = trans2nd(Tr.tdata,[olX olY]);
        [wpt] = trans2nd(Tr.tdata,mat1W(:,1:2));
    elseif polorder==4
        [xyW] = lwm_forward(Tr,[xl yl]); %Apply transformation
        [olw] = lwm_forward(Tr,[olX olY]);
        [wpt] = lwm_forward(Tr,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

    Tr_rec{d,t}=Tr;%To record trsfos - AG 2013
    
    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
        Tr_rec{d,t}=Tr;%To record trsfos - AG 2013
        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
    Tr_rec{d,t}=Tr;%To record trsfos - AG 2013
    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