function figure_show_original_and_modelled(modeldir)

% this function generates a figure where each of the original shapes are
% shown and next to that is the modelled shape using a number of PC's


% if nargin<1
%     modeldir = '';
% end
% if isempty(modeldir)
%     modeldir = uigetdir(pwd, 'Pick your model directory (it must contain a mod_pmd.mat)');
% end

[group_filename, group_pathname] = uigetfile('*.mat', 'Pick a group file');
if isequal(group_filename,0) || isequal(group_pathname,0)
    disp('User pressed cancel');
    return;
end
groups = load([group_pathname, filesep, group_filename]);
groups = groups.groups;

pdmfilename = [modeldir, filesep, 'mod_pdm.mat'];
sfamfilename = [modeldir, filesep, 'mod_sfam.mat'];
aammodel = [modeldir, filesep, 'model.aam_dat'];
if ~exist(aammodel, 'file')
    fprintf('There was no model in this directory\n');
    return;
end
aammodel = load(aammodel, '-mat');
aammodel = aammodel.aam;

if ~exist(pdmfilename, 'file')
    fprintf('There was no model in this directory\n');
    return;
end
pdm = load(pdmfilename);
pdm = pdm.pdm;

b = pdm.pca.b;
for i=1:length(b); str{i} = num2str(i); end;
[pcs,v] = listdlg('PromptString','Select PC''s to display:', 'SelectionMode','multiple', 'ListString',str)
if isempty(pcs)
    return;
end

figure; set(gcf, 'Color', [1 1 1]);
M = 3; N = 4;
indx =1;
template = get(aammodel, 'PointModelTemplate');
primaries = get(template, 'primaries');
edge = get(template, 'loops');
edge = edge{1};

dim = 2;
fillcol = 'r';
fillflag = 1;
edgecol = 'k';
P = pdm.pca.P(:, pcs);
all_axish = [];
for i=1:length(groups)
     if (indx>M*N)
        figure; set(gcf, 'Color', [1 1 1]);
        M = 3; N = 4;
        indx =1;
     end
    all_axish = cat(2, all_axish, subplot(M, N, indx)); indx = indx + 1;
     if isfield(groups(i).shapedata, 'pts_aligned')
    group_mean =  mean([groups(i).shapedata.pts_aligned],2);
    plotshape(group_mean(:), primaries, edge, [0 0], gca, fillcol, 0, edgecol);
     end
    p = get(get(gca, 'xlabel'), 'position');
    text(p(1), p(2), p(3), groups(i).name, 'FontAngle', 'italic', 'FontName', 'times', 'HorizontalAlignment', 'center', 'Interpreter', 'none')
    all_axish = cat(2, all_axish, subplot(M, N, indx)); indx = indx + 1;
    model_group_mean = pdm.Xm + P*(P'*(group_mean - pdm.Xm));
    plotshape(model_group_mean(:), primaries, edge, [0 0], gca, groups(i).color, fillflag, edgecol);
end
normalise_axis_limits(all_axish);
set(all_axish, 'Visible', 'off');
%%%%%
%
%
%%%%%
function plotshape(pts, pri, edge, offset, axish, fillcol, fillflag, edgecol)
[s, o] = isAugmentedProject;
pts = reshape(pts, 2, length(pts)/2);
pts = splitPoints(pts', s, o)';
pts = pts(:);
facecol = [1 1 1];
priedgecol = [0 0 0];
prifacecol = [0 0 0];
if nargin<8
    edgecol = [0 0 0];
end
pts = reshape(pts, 2, length(pts)/2);
pts(1,:) = pts(1,:) + offset(1);
pts(2,:) = pts(2,:) + offset(2);

hold(axish, 'on');
if fillflag
   % fillcol = [0.7529    0.7529    0.7529];
    fill(pts(1,:),pts(2,:),fillcol)
end 
for i=1:length(edge)
    xp = [pts(1,edge(i,1)) pts(1,edge(i,2))];
    yp = [pts(2,edge(i,1)) pts(2,edge(i,2))];
    plot(axish, xp, yp, '-', 'Color', edgecol, 'LineWidth', 2);
end
axis(axish, 'image', 'ij');
%%%
%
%%%
function normalise_axis_limits(all_axis_h)
newxlim = [realmax, 0];
newylim = [realmax, 0];
for i=1:length(all_axis_h)
   ah = all_axis_h(i); 
   xlim = get(ah, 'Xlim');
   ylim = get(ah, 'Ylim');
   if xlim(1)<newxlim(1)
       newxlim(1) = xlim(1);
   end
   if xlim(2)>newxlim(2)
       newxlim(2) = xlim(2);
   end
   if ylim(1)<newylim(1)
       newylim(1) = ylim(1);
   end
   if ylim(2)>newylim(2)
       newylim(2) = ylim(2);
   end
end
set(all_axis_h, 'Xlim', newxlim, 'Ylim', newylim);
return
