function [growth_params_bar] = calculate_growth_params(av_info, N, smooth)
% av_info, is a cell structure, where av_info{1} corresponds to smaller
% sectors, and so corresponds to later development, this should be used
% first when ungrowing. I.e. we use the growth factors that were calculated
% at the end of development to start the ungrow. In this way we will use
% the data from the larger sectors when then organ is nearly fully ungrown.
growth_params = cell(1, length(av_info)-1);
% remember that columns 4 and 5 correspond to L and W, and 6 is theta
if isempty(av_info)
    return;
end
L = zeros(length(av_info{1}(:,4)),1);
W = zeros(length(av_info{1}(:,4)),1);

% Recall that L(k-delta) is pulsed earlier and so will (hopefully have
% larger sectors) this means that L and W should > 0.
for i=2:length(av_info)
 
    L_k_delta = av_info{i}(:,4);
    L_k = av_info{i-1}(:,4);
    W_k_delta = av_info{i}(:,5);
    W_k = av_info{i-1}(:,5);
    
    for ii=1:length(L_k)
        if L_k(ii)>0
            L(ii) = L_k_delta(ii)./L_k(ii);
        end
    end
    for ii=1:length(W_k)
        if W_k(ii)>0
            W(ii) = W_k_delta(ii)./W_k(ii);
        end
    end
    growth_params{i-1} = [L, W, av_info{i-1}(:,6)];
end
% Now we need to smooth the length and width parameters with their
% respective neighbours.
% g = load([grid_path, filesep, grid_name]);
% g = g.region_info;
% N = g.neighbours;
%G1 = gtlib_growthParams2Tensor(av_info{2}(1, 4), av_info{2}(1, 5), av_info{2}(1, 6));
%G2 = gtlib_growthParams2Tensor(av_info{1}(1, 4), av_info{1}(1, 5), av_info{1}(1, 6));
%[a, b, c] = gtlib_growthTensor2Params(G1 - G2);
%growth_params{1}(1,:)
clear g;
growth_params = correct_principle_direction_of_growth(growth_params, N);
if smooth
    growth_params_bar = smooth_growth_params_tensor(growth_params, N);
else
    growth_params_bar = growth_params;
end
%growth_params_bar = correct_principle_direction_of_growth(growth_params_bar, N);
return;
%%%%%
%
%
%%%%%
function [growth_params_bar] = correct_principle_direction_of_growth(growth_params_bar, N)
for i=1:size(N,1)
    ind = find(N(i,:));
    for ii=1:length(growth_params_bar)
        if (growth_params_bar{ii}(i,1) > growth_params_bar{ii}(i,2))
            theta = growth_params_bar{ii}(i,3);
        else
            scratch = growth_params_bar{ii}(i,1);
            growth_params_bar{ii}(i,1) = growth_params_bar{ii}(i,2);
            growth_params_bar{ii}(i,2) = scratch;
            theta = growth_params_bar{ii}(i,3)+pi/2;
        end
        theta = mod(theta,pi);
        growth_params_bar{ii}(i,3) = theta;
    end
end
%%%%%
%
%
%%%%%
function [growth_params_bar] = smooth_growth_params_tensor(growth_params, N)
growth_params_bar = growth_params;
for i=1:size(N,1)
    ind = find(N(i,:));
    for ii=1:length(growth_params)
        L = growth_params{ii}(:,1);
        W = growth_params{ii}(:,2);
        theta = growth_params{ii}(:,3);
        M = [zeros(length(ind),3), L(ind), W(ind), theta(ind)];
        av = average_ellipse_av_growth_tensor(M);
  
        growth_params_bar{ii}(i,1) = av(4);
        growth_params_bar{ii}(i,2) = av(5);
       growth_params_bar{ii}(i,3) = av(6);
    end
end
return;
