function [Xm, P, b, pcaDat] = principal_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 (2006).
if isempty(X)
    Xm = [];
    P = [];
    b = [];
    pcaDat = [];
    return;
end
X = X';
[m, n] = size(X);
% Find mean shape and the covariance matrix
Xm = mean(X,1);
d = X - repmat(Xm, m, 1);
S = cov(d);
% Eigen-decomposition of covariance matrix
[V D] = eig(S);
% Sort in ascending order
D = diag(D);
[evals idx] = sort(D, 'descend'); 
evecs = V(:,idx);
% Find surviving eigenvectors
vr = cumsum(evals) / sum(evals);
if v==1
    t = length(evecs)
else
    t = min(find(vr >= v));
end
% Keep only modes 1:t
Xm = Xm(:);
P = evecs(:,1:t);
b = evals(1:t);
pcaDat.P = evecs;
pcaDat.b = evals;
pcaDat.v = vr;