% function [ output_args ] = demo3DDiffeoSynth( input_args )
%This demo illustrates different ways to sample from points in a
%diffeomorphic model.

scriptname = mfilename;
scriptpath = which(scriptname);
tokens = regexp(scriptpath, filesep);
cd(scriptpath(1:tokens(end)));

load('helamodel.mat')

p = model.cellShapeModel.positions;

rminds = any(isnan(p),2);

pos = p(~rminds,1:2);

c = 1;
gmmfitfile = ['gmmfit.mat'];  
if ~exist(gmmfitfile, 'file')
    objgmm = {};
    while c < 5
        tryagain = 0;
        while tryagain < 2
            try
                objgmm{c} = gmdistribution.fit(pos,c, 'replicates', 200);
                break
            catch
                tryagain = tryagain + 1;    
            end
        end
        c = c+1;
    end
    save(gmmfitfile, 'objgmm')
else
    load(gmmfitfile)
end


[~, aicind] = min(cellfun(@(x) x.AIC, objgmm));
[~, bicind] = min(cellfun(@(x) x.BIC, objgmm));

gmm_aic = objgmm{aicind};
gmm_bic = objgmm{bicind};



tri = delaunay(pos);

linx = linspace(min(pos(:,1)), max(pos(:,1)), 100);
liny = linspace(min(pos(:,2)), max(pos(:,2)), 100);

[x,y] = meshgrid(linx,liny);

for i = 1:length(x(:))
    gmm_pdf(i) = gmm_aic.pdf([x(i), y(i)]);
end

figure('color', 'w')
axis('tight')

scatter(x(:), y(:), 300, gmm_pdf, '.')
hold on

scatter(pos(:,1), pos(:,2), 300, '.', 'k')
triplot(tri, pos(:,1), pos(:,2), 'k')
saveas(gcf, 'gmmdemo.png', 'png')


positions = nan(size(model.cellShapeModel.positions,1),2);
positions(~rminds,:) = pos;


model.cellShapeModel.positions = pos;
model.nuclearShapeModel.positions = pos;
model.cellShapeModel.tessellation = tri;
model.nuclearShapeModel.tessellation = tri;

numimgs = 5;
for i = 1:numimgs
    %sample from the gmm ANDed to the convex hull of the triangulation
    simplex_idx = nan;
    while isnan(simplex_idx)
        point = gmm_aic.random;
        simplex_idx = tsearchn(model.cellShapeModel.positions,model.cellShapeModel.tessellation,point);
    end
    points(i,:) = point;
    
    param.position = point;
    
    
    savedir = [scriptpath '_img' num2str(i)];
    if ~exist(savedir, 'dir')
        mkdir(savedir)
    end
    cd(savedir);
    
    param.targetDirectory = [scriptname filesep 'img' num2str(i)];
    
    imgs{i} = model2img({model}, param);
end





