function [shape] = warp_organ_shape(Enew, Snew, Eold, Sold, order, shape, REGS)

%function [warped_edge] = warp_organ_shape(Enew, Snew, Eold, Sold, order, edge_petal)
% l = load('Staged\Data\UngrowthData\3_1_5_5.mat');
% g = load('Staged\Data\Grids\grid_data_5_by_5.mat');
% g = g.region_info;
% 
% Enew = l.ungrown_grid(2).E;
% Snew = l.ungrown_grid(2).S;
% Eold = l.ungrown_grid(1).E;
% Sold = l.ungrown_grid(1).S;
for i=1:size(shape,1)
    ii = 2;
    r = [];
    while isempty(r) & (ii<=length(REGS))
      region = REGS{ii};
      if any(region.edge==i)
          r = ii;
      end
      ii = ii + 1;
    end
    if iscell(order) & (length(order)<=r)
        order = order{r};
    else
        order = [];
    end
    if isempty(r)
        cent = Sold(:, end-1:end);
        p = shape(i,:);
        [val, r] = min(sqrt(sum((cent - ones(size(cent,1),1)*p).^2,2)));
        r = r+1;
    end
        
    [xnew, ynew] = get_region_vertices(Enew, Snew, r-1, order);
    [xnew, ynew] = interp_vertex_points(xnew, ynew);
    [xold, yold] = get_region_vertices(Eold, Sold, r-1, order);
    [xold, yold] = interp_vertex_points(xold, yold);
    T = cp2tform([xnew, ynew], [xold, yold],'polynomial', 2);
    u = double(shape(i,1));
    v = double(shape(i,2));
    shape(i,:) = [ones(length(u),1)  u  v  u.*v  u.^2  v.^2]*T.tdata;
end
return;
figure(1); clf; hold on;
plot(edge_petal(:,1), edge_petal(:,2), '.')
warped_edge = edge_petal;
indexes = edge_petal(:,3);
%[edge_petal, I, J] = unique(edge_petal(:,1:2), 'rows');
%edge_petal = [edge_petal, indexes(I)];
warped_ind = [];
for r=1:max(edge_petal(:,3))
    if iscell(order) & (length(order)<=r)
        order = order{r};
    else
        order = [];
    end
    [xnew, ynew] = get_region_vertices(Enew, Snew, r, order);
    [xnew, ynew] = interp_vertex_points(xnew, ynew);
    [xold, yold] = get_region_vertices(Eold, Sold, r, order);
    [xold, yold] = interp_vertex_points(xold, yold);
    
    
    ind = find(edge_petal(:,3)==r);

    
    T = cp2tform([xnew, ynew], [xold, yold],'polynomial', 2);
    u = double(edge_petal(ind,1));
    v = double(edge_petal(ind,2));
    warped_edge(ind,1:2) = [ones(length(u),1)  u  v  u.*v  u.^2  v.^2]*T.tdata;
    warped_edge(ind,3) = edge_petal(ind,3);
    %T = cp2tform([xold, yold],[xnew, ynew],'linear conformal');


    %[warped_edge(ind,1), warped_edge(ind,2)] = tformfwd(T, edge_petal(ind,1), edge_petal(ind,2));
    plot(edge_petal(ind,1), edge_petal(ind,2), 'ro');
    plot(warped_edge(ind,1), warped_edge(ind,2), 'go');
end
%[warped_edge, toind] = image_sort(warped_edge(:,1), warped_edge(:,2));
%warped_edge = [warped_edge, indexes(toind)];