function Export2ObjFormat()
[data, ungrown_grids] = getUngrowthData;
% pathname = uigetdir([pwd, filesep], 'Pick a Directory');
% if pathname == 0
%     return;
% end
[filename, pathname] = uiputfile('*.mat', 'Save as...');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel');
    return;
end



max_interp_iter = 0;
%grown_S = ungrown_grid(1).S;
%grown_E = ungrown_grid(1).Edgs;
%order = ungrown_grid(1).order;
%shape = double(ungrown_grid(1).shape);
%R = buildregions(ungrown_grid(1));
%shape = getshape(R);
tri_method=questdlg('How to trianulate?', 'Triangulation Question', 'Auto','Manual', 'Load', 'Auto');

if strcmp(tri_method, 'Manual')
    figure(1); clf; hold on;
    imagesc(R);
    axis image ij;
end
tri_ratio = 20;
[T, d] = getTriangulation(grown_S, grown_E, order, max_interp_iter, shape, tri_method, tri_ratio);
exportData(data, ungrown_grid, T, [pathname, filesep, filename], max_interp_iter);
%%%%%%
%
%
%%%%%%
function [T, d] = getTriangulation(S, E, order, max_interp_iter, shape, tri_method, tri_ratio)
T = [];
[d, P, C] = interp_points(S, E, order, max_interp_iter);

switch tri_method
    case 'Manual'
        T = triangulate_manual(d(:,1), d(:,2));
        [filename, pathname] = uiputfile('*.mat', 'Save as...');
        if isequal(filename,0) || isequal(pathname,0)
            disp('User pressed cancel')
            return;
        else
            save([pathname, filesep, filename], 'T');
            disp(['User selected ', fullfile(pathname, filename)])
        end
    case 'Load'
        [filename, pathname] = uigetfile('*.mat', 'Load triangulation file');
        if isequal(filename,0) || isequal(pathname,0)
            disp('User pressed cancel')
            return;
        else
            T = load([pathname, filesep, filename]);
            if ~isfield(T, 'T')
                return;
            end
            T = T.T;
            disp(['User selected ', fullfile(pathname, filename)]);
        end
    case 'Auto'
        T = delaunay(d(:,1), d(:,2));
        T = removebadtriangles(T, d, shape, tri_ratio);
    case 'Centroids'
        display_grid(S, E, order, gca, 'r', 1);

        d = [P; C];
        T = [];

        for i=1:length(order)
            o = order{i};
            rpts = P(o, :);
            [orpts, oind] = order_xy_pts(rpts);

            cind = [];
            for j =1:size(C,1)
                c = C(j,:);
                if inpolygon(c(1),c(2),orpts(:,1),orpts(:,2))
                    cind = j;
                end
            end
            oind = [oind; oind(1)];
            for j=2:length(oind)
                T = [T; o(oind(j-1)) o(oind(j)) size(P,1)+cind];
            end
        end
        fprintf('Using centroids\n');
end
%[T, d] = interp_tri(T, d, 1);
%T = removeOutsideFaces(T, d, shape);

return;
%%%%%%
%
%
%%%%%%
function [d, P, C] = interp_points(S, E, order, max_interp_iter)
d = [];
x = E(:, end-1); y = E(:, end);
xc = S(:, end-1); yc = S(:, end);

% Interpolate the points
for i=1:length(order)
    o = order{i};
    dx = x(o); dy = y(o);
    for ii=1:max_interp_iter
        [dx, dy] = interp_vertex_points(dx, dy);
    end
    x = [x; dx]; y = [y; dy];
end
P = [x, y];
P = unique(P, 'rows');
C = [xc, yc];
C = unique(C, 'rows');

x = [x; xc]; y = [y; yc];
[x, i, j] = unique(x, 'rows');
x = x(j);
[y, i, j] = unique(y, 'rows');
y = y(j);
d = [x, y];
%%%%%%
%
%
%%%%%%
function [l, G] = getUngrowthData
[filename, pathname] = uigetfile('*.mat', 'Pick some ungrowth data', [pwd, filesep, 'Staged\Data\UngrowthData\']);
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
    return;
else
    disp(['User selected ', fullfile(pathname, filename)]);
    l = load(fullfile(pathname, filename));
    G = l.ungrown_grids;
end
%%%%%%
%
%
%%%%%%
function exportData(l, G, T, filename, max_interp_iter)
% Remember here that the growth data and the growth parameters are the
% wrong way up. Because of the way we did the ungrow, we had to first
% ungrow with the small sector data and then finally work our way up to the
% large sector data. So the growth data at index, 1, is actually the last
% bit of growth to apply. Remember that we also need to include the final
% shape of the organ we are growing. This will give us a target area to
% reach in functions that use this data.

waitbar_h = waitbar(0, 'Exporting...');
T_orig = T;
GROWTHSTRUCT = struct([]);
g = G(end);
d = interp_points(g.S, g.Edgs, g.order, max_interp_iter);

    
% GROWTHSTRUCT(1).vertices = d;
% GROWTHSTRUCT(1).tri = T;
% GROWTHSTRUCT(1).tridata = [];

%for stage_indx=1:length(G)-1
for stage_indx=length(G):-1:2
    waitbar(stage_indx/length(G), waitbar_h);
    g = G(stage_indx);
    S = g.S;
    E = g.Edgs;
    R = buildregions(g);
    shape = getshape(R);
    if nargin<7
        d = interp_points(S, E, g.order, max_interp_iter);
    end

    growth_rate = l.growth_parameters.scaling_factors(stage_indx-1).G;
    anisotropy = l.growth_parameters.scaling_factors(stage_indx-1).A;
    theta = l.growth_parameters.scaling_factors(stage_indx-1).theta;
    remove_ind = [];
    
    tri_info = zeros(size(T,1),3);
    for t = 1:size(T,1)
        tri_center = sum(d(T(t,:), :))/3;
        region_number = R(round(tri_center(2)), round(tri_center(1)));
        if region_number == 0
            tri_info(t, :) = [0, 0 ,0];
            disp('Setting to unity');
        else
            GR = growth_rate(region_number);
            A = anisotropy(region_number);
            doublex = sqrt(2*A);
            doubley = doublex/A;
            rate_x = log2(doublex)/GR;
            rate_y = log2(doubley)/GR;
            tri_info(t, :) = [rate_x, rate_y, theta(region_number)];
            %tri_info(t, :) = [growth_rate(region_number), anisotropy(region_number), theta(region_number)];
        end
    end
    figure(1); clf; hold on;
    imagesc(R);
    triplot(T, d(:,1), d(:,2), 'w', 'LineWidth', 1);
    axis image ij;
    drawnow
    GROWTHSTRUCT(end+1).vertices = d;
    GROWTHSTRUCT(end).tri = T;
    GROWTHSTRUCT(end).tridata = tri_info;
end
g = G(1);
S = g.S;
E = g.Edgs;
d = interp_points(S, E, g.order, max_interp_iter);
GROWTHSTRUCT(end+1).vertices = d;
GROWTHSTRUCT(end).tri = T;
GROWTHSTRUCT(end).tridata = [];
save(filename, 'GROWTHSTRUCT');
close(waitbar_h);
fprintf('Done\n');
%%%%%
%
%
%%%%%
function R = buildregions(g)
E = g.Edgs(:,end-1:end);
R = zeros(ceil(max(E(:,2))), ceil(max(E(:,1))));
for i=1:length(g.order)
    p=g.Edgs(g.order{i}, end-1:end);
    r = roipoly(R, p(:,1), p(:,2));
    R = R + r*i;
end
return;
%%%%%
%
%
%%%%%
function shape = getshape(R)
B = R>0;
shape = fliplr(binim2SortedEdge(B));
return;


