function [ungrown_grids] = satlib_ungrow(varargin)
% function [ungrown_grids] = satlib_ungrow(varargin)
%
%
% Dr. A. I. Hanna (2007)

[curr_grid, extern_iter, spr_const, plot_data, lambda, damp, dt, tol] = parseargs(varargin);

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);

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}.deltaSt = 0;
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;
    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
            [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);
                [Smaj, Smin, theta] = gtlib_growthTensor2Params(GT);
                [junk, junk, dtheta] = gtlib_growthTensor2Params(region_rot(:,:,r));
                GT = region_rot(:,:,r)'*GT;

                [Smaj, Smin, theta] = gtlib_growthTensor2Params(GT);
                
                stage_grids{j}.ungrowth_tensors(:,:,r) =  gtlib_growthParams2Tensor(Smaj, Smin, theta);
                
                %stage_grids{j}.ungrowth_tensors(:,:,r) = region_rot(:,:,r)'*GT;
            end
        end
    end
 

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

    %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
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;
return
%%
%
function [curr_grid, maxiter, K, plot_data, lambda, damp, dt, tol] = parseargs(varargin)
varargin = varargin{1};
error(nargchk(0,inf,nargin));
if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error(['Arguments to ', mfilename, ' must come param/value in pairs.'])
end
curr_grid = [];
maxiter = 200;
K = 1;
plot_data = 0;
lambda = 0;
damp = 1;
dt = 0.01;
tol = 1e-5;
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'curr_grid'
            curr_grid = varargin{i+1};
        case 'maxiter'
            maxiter = varargin{i+1};
        case 'k'
            K = varargin{i+1};
        case 'plot_data'
            plot_data = varargin{i+1};
        case 'lambda'
            lambda = abs(varargin{i+1});
        case 'damp'
            damp = varargin{i+1};
        case 'dt'
            dt = varargin{i+1};
        case 'tol'
            tol = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to ', mfilename, '.  Name was ' varargin{i}])
    end
end
if isempty(curr_grid)
    error('Parameter curr_grid must be supplied');
end
return
