function show_group_model_estimates(modeldirectory)
% Get the shape model
if nargin<1
    modeldirectory = uigetdir(pwd, 'Pick a Model Directory');
    if modeldirectory == 0
        fprintf('User pressed cancel');
        return;
    end
end
[filename, pathname] = uigetfile('*.mat', 'Pick an Group File');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
    return;
else
    groupsfilename = fullfile(pathname, filename);
end
pdmfilename = [modeldirectory, filesep, 'mod_pdm.mat'];
detailsfilename = [modeldirectory, filesep, 'details.mat'];
modelfilename = [modeldirectory, filesep, 'model.aam_dat'];
aam = load(modelfilename, '-mat');
aam = aam.aam;
template = get(aam, 'PointModelTemplate');
[templatepath, templatename, ext, vers] = fileparts(get(template, 'name'));
pmdir = ['PointModels', filesep, templatename];

pdm = load(pdmfilename);
pdm = pdm.pdm;
if exist(detailsfilename, 'file')
    details = load(detailsfilename);
    details = details.details;
    opts = [details.scaling, details.rotation, details.translation];
else
    opts = [1 1 1];
end
edges = get(template, 'loops');
edges = edges{1};
% Get the groups file
grps = load(groupsfilename);
grps = grps.groups;
% Set the PM directory
% For each group, take the first name in that group and estimate that shape
% by the model.
var_covered = pdm.FractionPCs;
Xm = pdm.Xm;
P = pdm.pca.P;
b = pdm.pca.b;
for i=1:length(b); str{i} = num2str(i); end;
[s,v] = listdlg('PromptString','Select PCs:', 'SelectionMode','multiple', 'ListString',str);
if v == 0
    fprintf('User pressed cancel\n');
    return;
end
total_var = sum(b);
P = P(:,s);
b = b(s);
fprintf('Variance covered: %f%%\n', 100*(sum(b)/total_var));


M = ceil(min(10, length(grps)/2));
N = 4;
indx = 1;
fig = figure;
fs = 8;
set(fig, 'Color', [0 0 0]);
waitbar_h = waitbar(0, 'Calculating group estimates...');
lengrps = length(grps);
for i=1:lengrps
    waitbar(i/lengrps, waitbar_h);
    if indx>M*N
        indx = 1;
        fig = figure;
        set(fig, 'Color', [0 0 0]);
    end
    figure(fig);
    pmfilename = [pmdir, filesep, grps(i).elements{1}];
    grpname = grps(i).name;
    pts = load(pmfilename);
    pts = pts.pts;
    pts= pmalign(pts,opts,Xm);
    model_b = P'*(pts - Xm);
    model_pts = Xm + P*model_b;
    axis_h = subplot(M, N, indx); indx = indx+1;
    plotshape(pts, edges, 'w', axis_h);
    title(axis_h, ['Original - ', grpname], 'Interpreter', 'none', 'Color', 'w', 'FontSize', fs);
    set(axis_h, 'Color', [0 0 0]);
    axis_h = subplot(M, N, indx); indx = indx+1;
    plotshape(model_pts, edges, grps(i).color, axis_h);
    title(axis_h, ['Model - ', grpname], 'Interpreter', 'none', 'Color', 'w', 'FontSize', fs);
    set(axis_h, 'Color', [0 0 0]);
end
close(waitbar_h);
%%%%%%%%%%%%
%
%
%%%%%%%%%%%%
function plotshape(pts, edges, col, axis_h)
pts = reshape(pts, 2,length(pts)/2);
fill(pts(1,:), pts(2,:), col);
hold(axis_h, 'on');
for i=1:size(edges,1)
    p1 = pts(:, edges(i,1));
    p2 = pts(:, edges(i,2));
    plot(axis_h, [p1(1) p2(1)],[p1(2) p2(2)], '-', 'Color', col);
end
axis(axis_h, 'image');
