function figure_compare_selected_pcs
% function figure_compare_selected_pcs
%
% Example:
%
% figure_compare_selected_pcs;
%
% Dr. A. I. Hanna (2007)

[filename, pathname] = uigetfile('*.xls', 'Pick an Excel File Generated By a T-Test');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
    return;
end
[numeric, txt] = xlsread([pathname, filesep, filename], 'sheet2');
if isempty(numeric)
    str = sprintf('There was no information regarding T-Test found \nin this Excel file.\n\nPlease try again.');
    uiwait(msgbox(str,'Warning','modal'));
    return;
end

[filename, pathname] = uigetfile('*.mat', 'Pick a Groups File');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
    return;
end
groups = load([pathname, filesep, filename], '-mat');
groups = groups.groups;


[filename, pathname] = uigetfile('*.aam_dat', 'Pick a Model File');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
    return;
end
model = load([pathname, filesep, filename], '-mat');
model = model.aam;
if exist([pathname, filesep, 'mod_pdm.mat'], 'file')
    pdm = load([pathname, filesep, 'mod_pdm.mat'], '-mat');
    pdm = pdm.pdm;
else
    fprintf('No pdm found\n');
    return;
end

figuredirectory = uigetdir(pwd, 'Pick a Directory To Save The Figures');
if figuredirectory==0
    return;
end

template = get(model, 'PointModelTemplate');
templatename = get(template, 'name');
[templatepath, templatename] = fileparts(templatename);
edges = get(template, 'loops');
if iscell(edges) && (length(edges)>0)
edges = edges{1};
end
% scaling, rotation, translation
opts_params = get(model, 'procrustesParams');
opts = [opts_params.scaling opts_params.rotation opts_params.translation];
figh = figure;
M = 1;
N = 3;
ptcol = uisetcolor('Pick the Landmark Colour');
fillcol = uisetcolor('Pick the Fill Colour');

prompt={'Minimum x value: ','Maximum x value: ', 'Minimum y value: ','Maximum y value: '};
name='Limits';
numlines=1;
defaultanswer={'0','400', '0', '400'};
answer=inputdlg(prompt,name,numlines,defaultanswer);
xlim = [str2double(answer{1}), str2double(answer{2})];
ylim = [str2double(answer{3}), str2double(answer{4})];

   
%a1 = input('Minimum x value: ');
%a2 = input('Minimum x value: ');
%xlim = [a1 a2];
%a1 = input('Minimum y value: ');
%a2 = input('Maximum y value: ');
%ylim = [a1 a2];

for j=1:length(txt)
    g = txt{j};
    sig_pcs = numeric(:,j);
    sig_pcs = abs(sig_pcs(sig_pcs~=0));
    if ~isempty(sig_pcs)
        for i=1:length(groups)
            if strcmpi(g, groups(i).name)
                clf(figh);
                % This is the mean of the raw data
                [X] = load_group_pointmodels(templatename, groups(i).elements);
                [pointmodels] = pmalign(X, opts, pdm.Xm);
                Xm = mean(pointmodels,2);
                %xlim = [min(Xm(1:2:end))-30, max(Xm(1:2:end))+30];
                %ylim = [min(Xm(2:2:end))-30, max(Xm(2:2:end))+30];
                axish = subplot(M,N,2);
                plotshape(Xm, edges, axish, 1, ptcol, fillcol);
                axis(axish, 'image', 'ij');
                axis(axish, [xlim(1) xlim(2) ylim(1), ylim(2)]);
                title(sprintf('Mean Shape - Raw Points'));
                
                % This is the shape of the mean of the groups using the
                % first N PC's
                b = pdm.pca.P'*(Xm - pdm.Xm);
                Xm_model = pdm.Xm + pdm.P*b(1:size(pdm.P,2));
                mse_model = mean((Xm - Xm_model).^2,1);
                axish = subplot(M,N,1);
                plotshape(Xm_model, edges, axish, 0, ptcol);
             
                axis(axish, 'image', 'ij');
                axis(axish, [xlim(1) xlim(2) ylim(1), ylim(2)]);
                title(sprintf('Mean Shape - First %d PCs\nMSE = %3.3f', size(pdm.P,2), mse_model));
                
                % This is the shape of the mean of the group using the
                % significant PCs from the ttest or ranked ttest
                Xm_model_sigpcs = pdm.Xm + pdm.pca.P(:,sig_pcs) *b(sig_pcs);
                mse_sigpcs = mean((Xm - Xm_model_sigpcs).^2,1);
                axish = subplot(M,N,3);
                plotshape(Xm_model_sigpcs, edges, axish, 0, ptcol);
                axis(axish, 'image', 'ij');
                axis(axish, [xlim(1) xlim(2) ylim(1), ylim(2)]);
                title(sprintf('Mean Shape - Sig. PCs\n MSE = %3.3f', mse_sigpcs));
                saveas(figh, [figuredirectory, filesep, groups(i).name, '.fig'], 'fig');
                print(figh, '-dtiffnocompression', [figuredirectory, filesep, groups(i).name, '.tiff']);
            end
        end
    end
end
return;

function plotshape(pts, edges, axish, fillflag, ptcol, fillcol)
hold(axish, 'on');
if fillflag
	fill(pts(1:2:end), pts(2:2:end), fillcol, 'Parent', axish);
end
plot(axish, pts(1:2:end), pts(2:2:end), '.', 'MarkerFaceColor', ptcol, 'MarkerEdgeColor', ptcol,'Color', ptcol);
pts = reshape(pts, 2, length(pts)/2);
for k=1:size(edges,1)
    plot(axish, [pts(1,edges(k,1)), pts(1,edges(k,2))] , [pts(2,edges(k,1)), pts(2,edges(k,2))], '-', 'Color', ptcol); 
end
