function [Xm, P, b, pcaDat] = pca(AAM, indices, filenames,FractionPCs, in_context_flag, lda_shape_flag)
if nargin<4
    FractionPCs=0.95;
end
%v = 0.95;
v = FractionPCs;
% 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);
modelDirec=get(AAM,'modelDirec');
for e =1:length(elements)
    name = filenames{e};
    name = name(1:length(name)-4);
    %pts = load(['PointModels', filesep, templatename, filesep,  name, '_aligned']);
    pts = load(fullfile(modelDirec,[name,'_aligned']));
    pts = pts.pts;
    pts = reshape(pts, 2, length(pts)/2);

    if in_context_flag
        ind = setdiff(1:size(pts,2), indices);
        pts(:, ind) = meanshape(:, ind);
    else
        pts = pts(:, indices);
    end

    X(:, indx) = pts(:);
    indx = indx+1;
end
[Xm, P, b, pcaDat] = principle_component_analysis(X, v);

if lda_shape_flag==1
    % This is to sort out some LDA stuff...we might need to project back using
    % the PCA matrix.
    X = [];
    % We must reduce the dimensionality of the original data.
    groups_file = [modelDirec, filesep, 'groups.mat'];
    if exist(groups_file)
        groups = load(groups_file);
        groups = groups.groups;
        for i=1:length(groups)
            elements = groups(i).elements;
            indx = 1;
            for ii=1:length(elements)
                name=  elements{ii};
                name =name(1:end-7);
                pts = load(fullfile(modelDirec,[name,'_aligned']));
                pts = pts.pts;
                pts = reshape(pts, 2, length(pts)/2);
                pts = pts(:, indices);
                X(:,indx) = pts(:)';
                indx = indx + 1;
            end
            G{i} = (P'*X)';
        end
        [Xm, P, b, pcaDat] = lda(G, 3);
    end
end
return;
