function FEMGrowthTool(varargin)
%function FEMGrowthTool(varargin)
%
% Examples:
%
% FEMGrowthTool;
%
% FEMGrowthTool('growth_filename', 'Data/a.mat');
%
%
% Dr. A. I. Hanna (2007)
ad.data.meshInfo = [];
ad.data.X = [];
ad.data.growth_dir = pwd;
ad.data.showgrowthtensors = 1;
ad.data.growthparams = [];
ad.data.plothandles.meshhandle = [];
error(nargchk(0,inf,nargin));

if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error('%s %s %s', 'Arguments to', mfilename, 'must come param/value in pairs.')
end
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'growth_filename'
            ad = doLoadGrowthData(ad, varargin{i+1});
        case 'growth_directory'
            ad.data.growth_dir = varargin{i+1};
        otherwise
            error('%s\n%s %s', 'Unknown parameter name passed to show_region_areal_growth.', 'Name was', varargin{i});
    end
end

% Open, move, and get the handles from the figure file.
fig = openfig(mfilename, 'reuse');

set(fig, 'Visible', 'on');
movegui(fig, 'west');
handles = guihandles(fig);
%set(fig, 'Color', get(handles.controls_panel, 'BackgroundColor'));


% Set all the callback functions
set(handles.load_growth, 'callback', {@doLoadGrowth});
set(handles.startstopbtn, 'callback', {@doStartStop});
set(handles.exportgrowthmenu, 'callback', {@doExportData});
set(handles.fixnodesbtn, 'callback', {@doFixNodes});
set(handles.inittempslider, 'callback', {@doUpdateTempSlider});
set(handles.youngsslider, 'callback', {@doUpdateYoungsSlider});
set(handles.poissonslider, 'callback', {@doUpdatePoissonSlider});
set(handles.freeallnodes, 'callback', {@doFreeAllNodes});
set(handles.quitmenu, 'callback', {@doQuit});



% Initialize the application data structure
ad.figMain = fig;

ad.handles = handles;
ad.handles.growthfig = [];
ad.handles.growthaxes = [];
ad.data.plothandles.circles = [];
ad.data.plothandles.imh  = [];
setappdata(0,'FEMGrowthToolData',ad);
doUpdateTempSlider(ad.handles.inittempslider);
doUpdateYoungsSlider(ad.handles.youngsslider);
doUpdatePoissonSlider(ad.handles.poissonslider);
try
    uiwait(fig);
catch
    if ishandle(fig)
        delete(fig)
    end
end
return;
%%%
%
%%%
function doQuit(quitmenu, evd) %#ok<INUSD,INUSD>
ad = getappdata(0, 'FEMGrowthToolData');
delete(ad.figMain);
setappdata(0, 'FEMGrowthToolData', ad);
%%%
%
%%%
function doUpdateTempSlider(tempslider, evd) %#ok<INUSD,INUSD>
ad = getappdata(0, 'FEMGrowthToolData');
set(ad.handles.temptxt, 'String', num2str(get(tempslider, 'value')));
setappdata(0, 'FEMGrowthToolData', ad);
%%%
%
%%%
function doUpdateYoungsSlider(youngsslider, evd) %#ok<INUSD,INUSD>
ad = getappdata(0, 'FEMGrowthToolData');
set(ad.handles.youngstxt, 'String', num2str(get(youngsslider, 'value')));
setappdata(0, 'FEMGrowthToolData', ad);
%%%
%
%%%
function doUpdatePoissonSlider(poissonslider, evd) %#ok<INUSD,INUSD>
ad = getappdata(0, 'FEMGrowthToolData');
set(ad.handles.poissontxt, 'String', num2str(get(poissonslider, 'value')));
setappdata(0, 'FEMGrowthToolData', ad);

%%%
%
%%%
function doLoadGrowth(loadbtn, evd) %#ok<INUSD,INUSD>
ad = getappdata(0, 'FEMGrowthToolData');
[filename, pathname] = uigetfile('*.mat', 'Load Growth Data', [ad.data.growth_dir, filesep]);
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel');
    return;
else
    disp(['User selected ', fullfile(pathname, filename)])
end
ad = doLoadGrowthData(ad, [pathname, filesep, filename]);
setappdata(0, 'FEMGrowthToolData', ad);
%%%
%
%%%
function ad = doLoadGrowthData(ad, filename)
growthparams = load(filename);
if ~isfield(growthparams, 'GROWTHSTRUCT');
    fprintf('Invalid growth file\n');
    return;
end
growthparams = growthparams.GROWTHSTRUCT;
growthparams = calcShapeAreas(growthparams);
ad.data.growthparams = growthparams;
gparams = growthparams{1};
pts = gparams.vertices;
meshInfo = femlib_createBasicMesh(pts, gparams.tri);
X = meshInfo.X;
meshInfo = updateMeshGUIParams(ad, meshInfo, X);
ad.data.meshInfo = meshInfo;
ad.data.X = X;
if ishandle(ad.handles.growthfig)
    delete(ad.handles.growthfig);
end
ad = newFEMGrowthFig(ad);
return
%%%
%
%%%
function ph = plotCircles(circles, axish)
hold(axish, 'on');
ph = cell(length(circles),1);
for i=1:length(circles)
    z = circles{i};
    ph{i} = fill3(z(:,1), z(:,2), 10*ones(size(z,1),1),'r', 'Parent', axish);
end
hold(axish, 'off');
%%%
%
%%%
function meshInfo = updateMeshGUIParams(ad, meshInfo, X)
N = size(X ,1);
inittemp = str2double(get(ad.handles.temptxt, 'String'));
meshInfo.vertex_temp = -inittemp*ones(N,1);
youngsMod = str2double(get(ad.handles.youngstxt, 'String'));
poissonRat = str2double(get(ad.handles.poissontxt, 'String'));
meshInfo.D = femlib_getMaterialStiffnessMatrix(youngsMod, poissonRat);
return;
%%%%%
%
%
%%%%%
function circles = getCircles(X, tri)
theta = linspace(0,2*pi,20);
r = 2;
z = [r*cos(theta)', r*sin(theta)'];
circles = cell(length(tri),1);
for i=1:length(tri)
    ind = tri(i).ind;
    pts = X(ind,:);
    mu = mean(pts(:,1:2));
    circles{i} = z + ones(size(z,1),1)*mu;
end
return;
%%%%%
%
%
%%%%%
function doStartStop(startstopbtn, evd) %#ok<INUSD>
ad = getappdata(0,'FEMGrowthToolData');
if isempty(ad.data.growthparams)
    set(ad.handles.startstopbtn, 'Value', 0);
    return;
end
if get(startstopbtn, 'Value') == 1
    ad = newFEMGrowthFig(ad);
    ad.data.circles = getCircles(ad.data.X, ad.data.meshInfo.element);
    set(startstopbtn, 'String', 'Stop');
    ad = simulate(ad);
else
    set(startstopbtn, 'String', 'Start');
end
setappdata(0,'FEMGrowthToolData',ad);
%%%
%
%%%
% function setAxisLimits(shape, axish)
% minX = min(shape(:,1))-20;
% maxX = max(shape(:,1))+20;
% minY = min(shape(:,2))-20;
% maxY = max(shape(:,2))+20;
% axis(axish, [minX maxX minY maxY -20 20]);
% view(axish, 2);
% return;

%%%
%
%%%
function ad = simulate(ad)
ad = newFEMGrowthFig(ad);
growthparams = ad.data.growthparams;
meshInfo = ad.data.meshInfo;
X = ad.data.X;
meshInfo = updateMeshGUIParams(ad, meshInfo, X);
axish = ad.handles.growthaxes;
numgrowthlevels = length(growthparams);
s = 1;
indx = 1;

makemovie = get(ad.handles.makemoviechk, 'Value');
while (get(ad.handles.startstopbtn, 'value') == 1) && (s<numgrowthlevels)
    target_area_reached = 0;
    gparams = growthparams{s};
    fprintf('Reading next growth data\n');
    if ~isempty(gparams.tridata)
        for i=1:length(meshInfo.element)
            epsilon_1 = gparams.tridata(i,1);
            epsilon_2 = gparams.tridata(i,2);
            theta = gparams.tridata(i,3);
            tensor = gtlib_growthParams2Tensor(epsilon_1, epsilon_2, theta);
            meshInfo.element(i).preStrain = [tensor(1,1); tensor(2,2); 0; 0; 0; 2*tensor(1,2)];
        end
        fprintf('%d number of elements in mesh\n', length(meshInfo.element));
        maxIter = 1000;
        currentArea = growthparams{s}.initArea;
        nextArea = growthparams{s+1}.initArea;
        iter = 1;
        while (~target_area_reached) && (iter<maxIter) && (get(ad.handles.startstopbtn, 'value') == 1)
            Xold = X;
            meshInfo.X = X;
            meshInfo = femlib_wedgeIteration(meshInfo);
            X = meshInfo.X;
            pts = X(1:meshInfo.dorsalDivide,1:2);
            K = convhull(pts(:,1),pts(:,2),{'Qt','Pp'}) ;
            currentArea = polyarea(pts(K,1), pts(K,2));
            target_area_reached = currentArea>=nextArea;
            fprintf(' - current area %2.2f, Next Area %2.2f\n', currentArea, nextArea);
            if target_area_reached==0
                ad.data.circles = updateCircles(Xold, X, meshInfo.element, ad.data.circles, gparams.tridata);
                updateCirclesPlot(ad.data.plothandles.circles, ad.data.circles);
                femlib_updateMeshPlot(X, meshInfo.element, meshInfo.dorsalDivide, ad.data.plothandles.meshhandle, axish);
                axis(axish, 'image', 'vis3d');
                drawnow;
                if makemovie
                    if sum(sum((X(:,1:2) - Xold(:,1:2)).^2,2))>.1
                        createMovieFiles(ad.handles.growthfig , indx);
                        fprintf('-added to movie\n');
                    end
                end
            end
            iter = iter + 1;
            indx = indx + 1;
        end
        X = Xold;
        s = s + 1;
    else
        set(ad.handles.startstopbtn, 'value', 0);
        doStartStop(ad.handles.startstopbtn);
    end
end
axis(axish, 'image', 'vis3d');
if get(ad.handles.startstopbtn, 'value')
    set(ad.handles.startstopbtn, 'value', 0)
    doStartStop(ad.handles.startstopbtn);
end
return;
%%%
%
%%%
function createMovieFiles(figh, indx)
%%%
%
%%% This is where we capture images for gifs and movies %%%
% I = getframe(figh);
% I = frame2im(I);
% moviefilesdir = [pwd, filesep, 'MovieImages'];
% if ~exist(moviefilesdir, 'dir')
%     mkdir(moviefilesdir);
% end
% imwrite(I, [moviefilesdir, filesep, sprintf('%2.3d.tiff', indx)], 'TIFF', 'Compression', 'none');

if indx==1
    I = getframe(figh);
    I = frame2im(I);
    [Xd, map] = rgb2ind(I, 128);
    imwrite(Xd, map, 'test.gif', 'GIF', 'WriteMode', 'overwrite', 'DelayTime', .5, 'LoopCount', Inf);
else
    I = getframe(figh);
    I = frame2im(I);
    [Xd, map] = rgb2ind(I, 128);
    imwrite(Xd, map, 'test.gif', 'GIF', 'WriteMode', 'append', 'DelayTime', .5);
end
%
%%%

%%%
%
%%%
function updateCirclesPlot(circh, circles)
for i=1:length(circh)
    ph = circh{i};
    z = circles{i};
    set(ph, 'XData', z(:,1), 'YData', z(:,2));
end
%%%
%
%%%
function circles = updateCircles(Xold, Xnew, tri, circles, tri_tensors)
for i=1:length(tri)
    ind = tri(i).ind;
    xold = Xold(ind,1:2);
    xnew = Xnew(ind,1:2);
    g1 = interpTriangle(xold');
    g2 = interpTriangle(xnew');
    T = inv(g1'*g1)*g1'*g2;
    z = circles{i};
    z = z - ones(size(z,1),1)*mean(z);
    znew = z*T;
    znew = znew + ones(size(znew,1),1)*mean(xnew);
    circles{i} = znew;
end
return;
%%%
%
%%%
function t = interpTriangle(t)
t = t - mean(t,2)*ones(1,3);
t = interp_tri_pts(t');
return;
%%%
%
%%%
function GP = calcShapeAreas(growthparams)
GP = cell(length(growthparams),1);
for i=1:length(growthparams)
    gparams = growthparams(i);
    pts = gparams.vertices;
    K = convhull(pts(:,1),pts(:,2),{'Qt','Pp'}) ;
    gparams.initArea = polyarea(pts(K,1), pts(K,2));
    %gparams.initArea = femlib_trisurfacearea(gparams.tri, pts);
    GP{i} = gparams;
end
return;
%%%%%
%
%%%%%
function ad = newFEMGrowthFig(ad)
if (isempty(ad.handles.growthfig) || isempty(ad.handles.growthaxes)) || (~ishandle(ad.handles.growthfig) || ~ishandle(ad.handles.growthaxes))
    ad.handles.growthfig = figure;
    cameratoolbar(ad.handles.growthfig);
    set(ad.handles.growthfig, 'Name', 'FEM Mesh Window');
    ad.handles.growthaxes = gca;
    movegui(ad.handles.growthfig, 'center');
    set(ad.handles.growthfig, 'DoubleBuffer', 'on');
    set(ad.handles.growthaxes, 'DrawMode', 'fast');
    final_tri = ad.data.growthparams{end}.tri;
    final_pts = ad.data.growthparams{end}.vertices;

    tr = mean(final_pts) - mean(ad.data.X(:,1:2));

    ad.data.X(:,1:2) = ad.data.X(:,1:2) + ones(size(ad.data.X(:,1:2),1),1)*tr;

    th = trisurf(final_tri, final_pts(:,1), final_pts(:,2), -200*ones(size(final_pts,1),1), 'Parent', ad.handles.growthaxes, 'FaceAlpha', 0.5);
    hold(ad.handles.growthaxes, 'on');
    % Now plot the mesh
    ad.data.plothandles.meshhandle = femlib_plotMesh('pts', ad.data.X, 'mesh', ad.data.meshInfo, 'axish', ad.handles.growthaxes, 'edgecolor', [0 .5 0]);
    axis(ad.handles.growthaxes, 'image', 'ij');
    %setAxisLimits(final_pts, ad.handles.growthaxes);
    set(ad.handles.growthfig, 'Color', [1 1 1]);
    set(ad.handles.growthaxes, 'visible', 'off');
    set(ad.handles.growthaxes, 'Color', [1 1 1], 'XColor', [1 1 1],  'YColor',  [1 1 1],'ZColor', [1 1 1], 'box', 'on');
    grid(ad.handles.growthaxes, 'on');
    view(ad.handles.growthaxes, 2);
    if ad.data.showgrowthtensors
        ad.data.circles = getCircles(ad.data.X, ad.data.meshInfo.element);
        ad.data.plothandles.circles = plotCircles(ad.data.circles, ad.handles.growthaxes);
    end
end
return
%%%
%
%%%
function doFreeAllNodes(freenodes, evd)
ad = getappdata(0,'FEMGrowthToolData');
ad.data.meshInfo.fixedEquations = [];
setappdata(0,'FEMGrowthToolData', ad);

%%%
%
%%%
function doFixNodes(fixnodes, evd)
ad = getappdata(0,'FEMGrowthToolData');
axes(ad.handles.growthaxes);
[x, y] = ginput(1);
mesh = ad.data.meshInfo;
X = mesh.X;
d = matdistance(X(:,1:2), [x;y]');
[val, ind] = min(d);
mesh.fixedEquations = [mesh.fixedEquations, ind*3-2, ind*3-1, ind*3];
ad.data.meshInfo = mesh;
setappdata(0,'FEMGrowthToolData', ad);
%%%
%
%%%
function doExportData(exportbtn, evd) %#ok<INUSD,INUSD>
ad = getappdata(0,'FEMGrowthToolData');
[filename, pathname] = uigetfile('*.mat', 'Open growth data file...');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
    return;
end
objdirectory = uigetdir(pwd, 'Select directory to save obj files in.');
if objdirectory==0
    disp('User pressed cancel')
    return;
end
growthstruct2obj([pathname, filesep, filename], objdirectory);
setappdata(0,'FEMGrowthToolData',ad);
return;
%%%
%
%%%
function growthstruct2obj(growthparamsfilename, objdirectory, t)
% function growthstruct2obj(growthparamsfilename)
%
%
% Dr. A. I. Hanna
growthparams = load(growthparamsfilename);
growthparams = growthparams.GROWTHSTRUCT;
s = length(growthparams);
filename = [objdirectory, filesep, 'stage_'];
if nargin<3
    t = 20;
end
for s = 1:s
    gparams = growthparams(s);
    npts = gparams.vertices;
    ntri = gparams.tri;
    tri_info = gparams.tridata;
    for i=1:size(tri_info, 1)
        e1 = gparams.tridata(i,1);
        e2 = gparams.tridata(i,2);
        theta = gparams.tridata(i,3);
        tri_info(i,:) = [e1, e2, theta];
    end
    filenamenum = sprintf('%.4d', s);
    save_data([filename, filenamenum, '.obj'], npts, ntri, tri_info);
end
%%%%%
%
%
%%%%
function save_data(filename, d, T, tri_info)
fprintf('Saving as %s.\n', filename);
if ~isempty(filename)
    fid = fopen(filename,'w');
    for i=1:size(d,1)
        fprintf(fid, 'v %5.5f %5.5f %5.5f\n', d(i,1), d(i,2), 0);
    end
    for i=1:size(T,1)
        fprintf(fid, 'f %d %d %d\n', T(i,1), T(i,2), T(i,3));
    end
    for i=1:size(tri_info, 1)
        fprintf(fid, 'x-growth %d %5.5f %5.5f %5.5f\n', i, tri_info(i,1), tri_info(i,2), tri_info(i,3));
    end
    fclose(fid);
end
return;

