function [ratio, D, projD] = aamlib_calcProjError(varargin)
%function [ratio, D, projD] = aamlib_calcProjError(varargin)
%
%
% [ratio, D, projD] = aamlib_calcProjError('verbose', 1);
%
% Dr. A. I. Hanna (2007)
if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error(['Arguments to ', mfilename,' must come param/value in pairs.'])
end
files = [];
model = [];
mse = [];
X = [];
pcind = [];
verbose = 1;
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'files'
            files = varargin{i+1};
        case 'model'
            model = varargin{i+1};
        case 'pointmodels'
            X = varargin{i+1};
        case 'pcind'
            pcind = varargin{i+1};
        case 'verbose'
            verbose = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to ', mfilename, '.  Name was ' varargin{i}])
    end
end
if isempty(model)
    [model] = aamio_loadModel;
    if isempty(model)
        return;
    end
    Xm = model.Xm;
    P = model.pca.P;
    if isempty(pcind)
        [P, pcind] = choosePCs(P);
        if isempty(P)
            fprintf('User cancelled operation\n');
            return;
        end
    end
end
if isempty(X)
    if isempty(files)
        [X, files] = aamio_loadPointModels;
        if isempty(X)
            return;
        end
    else
        [X, files] = aamio_loadPointModels('filenames', files);
    end
end

if size(Xm,1)~=size(X,1)
    error('Number of points must agree');
end
% Align the point models to the mean, ignoring scale, rotation and
% translation
opts.scaling = 1;
opts.rotation = 1;
opts.translation = 1;
[X] = pcalib_GPA('data',X, 'template', Xm, 'opts', opts);

% get the length of the vector in PCA space using all the PCs
fullP = model.pca.P;
N = size(X,2);
D = zeros(N,1);
for i=1:N
    B = fullP'*(X(:,i) - Xm);
    D(i) = norm(fullP'*(X(:,i) - Xm));
end

% get the length of the vector in PCA space using SELECTED PCs
projD = zeros(N,1);
Xstar = X;
for i=1:N
    b = P'*(X(:,i) - Xm);
    projD(i) = norm(b);
    Xstar(:,i) = Xm + P*b;
end


% calculate the ratio between the projected vector and the original vector
% as requested by Rico.

ratio = projD./D;
if verbose
    plotData(Xm, X, 'Original Data');
    plotData(Xm, Xstar, 'Data Projected onto Selected PCs');
end
excel = questdlg('Would you like to save the result to Excel?', ...
                         'Save?', 'Yes', 'No', 'Yes');
if isempty(excel)
    return;
end
if strcmpi(excel, 'yes')
    [filename, pathname] = uiputfile('*.xls', 'Save as...');
    if isequal(filename,0) || isequal(pathname,0)
        disp('User pressed cancel');
        return
    end
    sheet = createExcelFormat(ratio, D, projD, files);
    xlswrite(fullfile(pathname, filename), sheet);
end
return
%%%
%
%%%
function plotData(Xm, X, str)
figure;
hold on;
plot(Xm(1:2:end), Xm(2:2:end), 'r');
for i=1:size(X,2)
    plot(X(1:2:end, i), X(2:2:end, i), 'b:');
end
axis image xy;
title(str);
return

%%%
%
%%%
function [P, s] = choosePCs(P)
N = size(P,2);
str = cell(N,1);
for i=1:N
    str{i} = num2str(i);
end
[s,v] = listdlg('Name', 'PC Selection', 'PromptString','Select PCs to use in reconstruction:',...
    'SelectionMode','multiple',...
    'ListString',str, 'ListSize', [600, 300]);
P = P(:,s);
return
%%%
%
%%%
function sheet = createExcelFormat(ratio, D, projD, files)
sheet = {};
num_files = length(files);
sheet{1,1} = 'Filename';
sheet{1,2} = 'Ratio';
sheet{1,3} = 'Original Distance';
sheet{1,4} = 'Projected Distance';
for i=1:num_files
    [junk, filename] = fileparts(files{i});
    sheet{i+1,1} = filename;
    sheet{i+1,2} = ratio(i);
    sheet{i+1,3} = D(i);
    sheet{i+1,4} = projD(i);

end
return