function [stage_grid] = calcRegionalGrowthTensors(stage_grid, sectors_pre, sectors_post, min_number_sectors, smooth)

% for each region that is not empty, merge surrounding regions if they are
% empty.

centroids_prev = [sectors_pre.Centroid]';
centroids_prev = reshape(centroids_prev, 2, length(centroids_prev)/2);
centroids_post = [sectors_post.Centroid]';
centroids_post = reshape(centroids_post, 2, length(centroids_post)/2);

[density.prev] = calcDeficientRegions(stage_grid.tri_info.T, stage_grid.tri_info.pts, stage_grid.region_tri_ind, centroids_prev, min_number_sectors);
[density.post] = calcDeficientRegions(stage_grid.tri_info.T, stage_grid.tri_info.pts, stage_grid.region_tri_ind, centroids_post, min_number_sectors);

stage_grid.density = density;

G_t = stage_info2average_region_info(stage_grid, sectors_pre);
G_t_delta = stage_info2average_region_info(stage_grid, sectors_post);


num_regions = size(G_t,3);
growth_tensors = zeros(3,3,num_regions);

for j = 1:num_regions
    [Emaj1, Emin1, theta1] = gtlib_growthTensor2Params(G_t(:,:,j));
    [Emaj2, Emin2, theta2] = gtlib_growthTensor2Params(G_t_delta(:,:,j));
    Smaj = (Emaj2/Emaj1);
    Smin = (Emin2/Emin1);
    [Smaj, Smin] = checkShrinking(Smaj, Smin);
    if Smin>Smaj
        theta2 = theta2 + pi/2;
    end
    growth_tensors(:,:,j) = gtlib_growthParams2Tensor(Smaj, Smin, theta2);
end
if smooth
    growth_tensors = smooth_tensors(growth_tensors, stage_grid);
end
stage_grid.ungrowth_tensors = growth_tensors;
stage_grid.average_sectors.prev = G_t;
stage_grid.average_sectors.post = G_t_delta;
return
%%
function [Smaj, Smin] = checkShrinking(Smaj, Smin)
if Smaj<1
    Smaj = 1 + realmin;
end
if Smin<1
    Smin = 1 + realmin;
end
return
%%
function smooth_growth_tensors = smooth_tensors(growth_tensors, stage_grid)
num_regions = size(growth_tensors,3);
T = stage_grid.tri_info.T;
regions_ind = stage_grid.region_tri_ind;
smooth_growth_tensors = growth_tensors;
for i=1:num_regions
    neigh_tri = findTouchingTriangles(T,  stage_grid.region_tri_ind{i});
    neighbours = unique(findRegionContainingTri(regions_ind, neigh_tri));

        neighbours(neighbours == i) = [];
    if ~isempty(neighbours)
        neighbours = cat(1, neighbours, i);
        smooth_growth_tensors(:,:,i) = mean(growth_tensors(:,:,neighbours),3);
    end
end
return