function [gridinterpX,gridinterpY,gridinterpZ]= growth_fromgrid_ag3D(leaflineg,gridX1,gridY1,gridZ1,gridXW,gridYW,gridZW,opt)
%opt=1: display
if opt==1
ha2=figure;
h=line(leaflineg(:,1),leaflineg(:,2));
hold on
end
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)]';
        zcoord1=[gridZ1(r,c) gridZ1(r,c+1) gridZ1(r+1,c+1) gridZ1(r+1,c)]';
        %area1=polyarea(xcoord1,ycoord1);
        %decompose quadrilateral (which may be non planar) into 2 triangles
        %and add areas of the 2 triangles
        area1=quadriarea3D(xcoord1,ycoord1,zcoord1);
        
        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)]';
        zcoord=[gridZW(r,c) gridZW(r,c+1) gridZW(r+1,c+1) gridZW(r+1,c)]';
        %area=polyarea(xcoord,ycoord);
        area=quadriarea3D(xcoord,ycoord,zcoord);
        
        RG=100*(area-area1)/area1;
        if opt==1
        h=patch(xcoord1,ycoord1,RG);
        end
        [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

if opt==1
%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');

end