function [Xm, P, b, pcaDat] = pca(AAM, indices, filenames,FractionPCs)if nargin<4    FractionPCs=0.95;end%v = 0.95;%v = 0.99;% if(nargin < 1) error('No data defined for PCA'); endelements = get(AAM, 'elements');activeElements = get(AAM, 'activeElements');indx = 1;pmt = get(AAM, 'PointModelTemplate');templatename = get(pmt, 'name');templatename = templatename(1:length(templatename)-9);for e =1:length(elements)    name = filenames{e};    name = name(1:length(name)-4);    pts = load(['PointModels', filesep, templatename, filesep,  name, '_aligned']);    pts = pts.pts;    pts = reshape(pts, 2, length(pts)/2);    pts = pts(:, indices);    X(:, indx) = pts(:);    indx = indx+1;endX = X';[m, n] = size(X);% Find mean shape and the covariance matrixXm = 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 lengthif(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    endendevecs = V(:,idx);evals = flipud(evals);evecs = fliplr(evecs);% Find surviving eigenvectorsvr = cumsum(evals) / sum(evals);t = min(find(vr >= v));% Keep only modes 1:tXm = Xm(:);P = evecs(:,1:t);b = evals(1:t);pcaDat.P = evecs;pcaDat.b = evals;pcaDat.v = vr;