pca.asv

matlab aamtool box

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'); end
elements = 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;
end
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;