function [grids, growth_params] = ungrow_assumption_three(R, growth_params, extern_iter, spr_const, recalc_flag, animate, axis_handle, error_axis_handle)
% function [R, growth_params] = ungrow_assumption_three(R, growth_params, extern_iter, spr_const, recalc_flag, animate, axis_handle, error_axis_handle)
%
%
% Dr. A. I. Hanna (2007)
if nargin<2
    axis_handle = gca;
end
cla(error_axis_handle);
num_growth_stages = length(growth_params.growth_params);
grids = cell(num_growth_stages+1,1);
grids{1} = R.region_info;
grids{1}.R = [];
for i=1:num_growth_stages
    SF = growth_params.scaling_factors{i};
    xy_before = R.region_info.vert;
    R = warp_springs(R, SF);
    [R, iter, plotE] = relax_springs(R, extern_iter, spr_const, animate, axis_handle, error_axis_handle);
    R.region_info.centroids = R.region_info.vert(R.region_info.num_ctrl_pts+1:end,:);
    xy_after = R.region_info.vert;
    
    % warp the shape of the organ
    R.region_info.shape = tpslib_ptswarp(R.region_info.shape, xy_before, xy_after);    
    R.region_info.region_vertices = tpslib_ptswarp(R.region_info.region_vertices, xy_before, xy_after);
    
    if recalc_flag
        [region_rot_mat] = calc_strain_cross_matrix(xy_before, xy_after, R.region_info.region_ctrl_pt_ind);
        fprintf(' - rotating tensors for stages %0.2d - %0.2d\n', i, num_growth_stages);
        for j=i+1:num_growth_stages
            growth_params.scaling_factors{j}.theta = growth_params.scaling_factors{j}.theta - region_rot_mat;
        end
    end
    grids{i+1} = R.region_info;
    grids{i+1}.R = [];
    %plot_springs(R.region_info, axis_handle);
    plot_centroids(R.region_info, axis_handle);
    plot_shape(R.region_info, axis_handle);
    plot_grid(R.region_info, axis_handle);
    axis(axis_handle, 'ij');
    drawnow;
end
return;
