function [Xm, P, b, pcaDat] = principle_component_analysis(X, v)
% function [Xm, P, b, pcaDat] = principle_component_analysis(X)
% A function that performs PCA on a matrix of data, where each column is an
% experiment, and each row corresponds to points in that experiment.
%
% re-arranged Dr. A. I. Hanna (2005).
X = X';
[m, n] = size(X);
% Find mean shape and the covariance matrix
Xm = mean(X,1);
d = X - repmat(Xm, m, 1);

if(m >= n)
	% Compute normal covaraince	
	S = d' * d / (m - 1);
else
	% Covariance is not full-rank
	S = d * d' / (m - 1);
end
% Eigen-decomposition of covariance matrix
[V D] = eig(S);
% Sort in ascending order
[evals idx] = sort(diag(D)); 
% Note if fast method was used, eigengevectors
% aren't orthonormal - rescale to unit length
if(m < n)
    V = d' * V;
    for k = 1:size(V,2)
        if sqrt(V(:,k)' * V(:,k))==0
            uiwait(msgbox(sprintf('There is no variation in the shapes being analysed.'),'Oh Dear','modal'));
            %error('');
            Xm = [];
            P = [];
            b = []; 
            pcaDat = [];
            return;
        else
            V(:,k) = V(:,k) / sqrt(V(:,k)' * V(:,k));
        end
    end
end
evecs = V(:,idx);
evals = flipud(evals);
evecs = fliplr(evecs);
% Find surviving eigenvectors
vr = cumsum(evals) / sum(evals);
t = min(find(vr >= v));
% Keep only modes 1:t
Xm = Xm(:);
P = evecs(:,1:t);
b = evals(1:t);
pcaDat.P = evecs;
pcaDat.b = evals;
pcaDat.v = vr;