function [ungrown_grids] = ungrow_aih(curr_grid, extern_iter, spr_const, plot_data, lambda, damp, dt, tol)
% function [R, ungrowth_params] = ungrow_aih(R, ungrowth_params, extern_iter, spr_const, recalc_flag, animate, axis_handle, error_axis_handle)
%
%
% Dr. A. I. Hanna (2007)

animate =  get(plot_data.animate_edge_chk, 'value');
recalc_flag = get(plot_data.recalc_chk, 'Value');
axis_handle = plot_data.imageAxis;
error_axis_handle = plot_data.error_axis;
tensor_axis = plot_data.tensor_axis;
cla(error_axis_handle);
stage_grids = curr_grid.stage_grids;
num_growth_stages = length(stage_grids);
%num_regions = size(ungrowth_params.ungrowth_tensors{1},3);


base_grid = curr_grid.base_grid;

reinit_grid = 1;
original_stage_grids = stage_grids;
updated_tensor_colour = [0 0.2510 0.5020];
original_tensor_colour = [1 0 0];

ungrown_grids = cell(num_growth_stages+1,1);

ungrown_grids{1} = base_grid;
ungrown_grids{1}.ungrowth_tensors = zeros(3,3, length(base_grid.region_tri_ind));

for i=1:num_growth_stages

    cla(tensor_axis);
    hold(tensor_axis, 'on');
    drawStageGrowthTensors(stage_grids{i}, updated_tensor_colour, tensor_axis);
    drawStageGrowthTensors(original_stage_grids{i}, original_tensor_colour, tensor_axis);
    set(tensor_axis, 'visible', 'off');

    stage_ungrowth_tensors = stage_grids{i}.ungrowth_tensors;

    xy_before = base_grid.tri_info.pts;

    stage_springs = stage_grids{i}.tri_info.spring_edge;

    new_spring_lengths = warp_springs(base_grid, stage_ungrowth_tensors, stage_springs);

    base_grid = satlib_shrinkMesh(stage_ungrowth_tensors, base_grid);
    

    % relax the mesh according to the calculated ungrowth data and user
    % specified spring parameters
    [base_grid, iter, plotE] = satlib_relaxSprings('lambda', lambda, 'damp', damp, 'dt', dt, 'tol', tol, 'S', base_grid, ...
        'stage_springs', stage_springs, 'new_spring_lengths', new_spring_lengths, ...
        'maxiter', extern_iter, 'K', spr_const, 'animate', animate, 'image_axis', axis_handle, 'error_axis', error_axis_handle);

    xy_after = base_grid.tri_info.pts;
    % xy_after = R.region_info.spring.vert;
    stage_region_tri_ind = stage_grids{i}.region_tri_ind;
    base_grid = warp_and_update(base_grid, xy_before, xy_after);

    if recalc_flag
        for j=i:num_growth_stages
            %grid_id = getGridId(organ_grid, ungrowth_params.stages{j});
            [region_rot] = calc_strain_cross_matrix(xy_before, xy_after, stage_grids{j}.region_tri_ind, base_grid.tri_info.T);
            fprintf(' - rotating tensors for stages %0.2d - %0.2d\n', i, num_growth_stages);
            num_regions = size(region_rot,3);
            for r = 1:num_regions
                GT = stage_grids{j}.ungrowth_tensors(:,:,r);
                stage_grids{j}.ungrowth_tensors(:,:,r) = region_rot(:,:,r)'*GT;
            end
        end
    end
    %     if reinit_grid
    %         R.region_info = init_grid(R.region_info);
    %     end
    %[R.region_info] = makeExternalSprings( R.region_info);
    %[R.region_info] = makeRegionSprings( R.region_info);

    ungrown_grids{i+1} = base_grid;
    ungrown_grids{i+1}.region_tri_ind = stage_grids{i}.region_tri_ind;
    ungrown_grids{i+1}.ungrowth_tensors = stage_grids{i}.ungrowth_tensors;

    %grids{i+1} = organ_grid;
    %grids{i+1}.R = [];
    %plot_springs(organ_grid, axis_handle);
    %plot_centroids(R.region_info, axis_handle);
    plot_shape(base_grid, axis_handle);
    plot_grid(base_grid, axis_handle);
    axis(axis_handle, 'ij');
    drawnow;
end
return;
%%%
%
%%%
function R = warp_and_update(R, xy_before, xy_after)
% warp the shape of the organ
%R.region_info.region_vertices = tpslib_ptswarp(R.region_info.region_vertices, xy_before, xy_after);
T = cp2tform(xy_after, xy_before,'polynomial', 2);
u = R.shape(:,1);
v = R.shape(:,2);
R.shape = [ones(length(u),1)  u  v  u.*v  u.^2  v.^2]*T.tdata;
% for i=1:length(R.centroids)
%     c = R.centroids{i};
%     u = c(:,1);
%     v = c(:,2);
%     c = [ones(length(u),1)  u  v  u.*v  u.^2  v.^2]*T.tdata;
%    R.centroids{i} = c;
% end
return
%%%
%
%%%
function curr_grid = init_grid(curr_grid)
curr_grid = segmentOrgan(curr_grid);

[curr_grid] = makeInternalSprings(curr_grid);
[curr_grid] = makeExternalSprings(curr_grid);
[curr_grid] = makeRegionSprings(curr_grid);

%[curr_grid] = makeInternalSprings(curr_grid);
%[curr_grid] = makeExternalSprings(curr_grid);
%[curr_grid] = makeShapeSprings(curr_grid);
return
