getsfam.m

matlab aamtool box

%
% function [sfam, options] = getsfam(image_folder, pt_model_folder, resolution, scale, rotate, trans, v, triangulation)
%
% Input Parameters
%     image_folder  - the location of your images (probably Cropped)
%     pt_model_folder - the folder containing your point models, _pt.mat
%     resolution - the number of pixels you want to include in the
%     appearance model
%     scale - include scale into the procrustes
%     rotation - include rotation into the procrustes
%     translation - include translation into the procrustes
%     v  - the amount of variance to be captured.
%
% Output Parameters:
%     sfam - the shape and appearance structure
%    options - the options file used to build the models.
%
% Example:
% 
% [sfam, options] = getsfam('Cropped', 'PointModels\PetalTempalte1', 200, 1, 1, 1, .95, []);
% 
% Dr. A. I. Hanna (2005)
function [sfam, options] = getsfam(image_folder, pt_model_folder, resolution, scale, rotate, trans, v, triangulation)
% image_folder = 'Cropped';
% pt_model_folder = 'PointModels\PetalTempalte1';
% res = 200;
% scale = 1;
% rotate = 1;
% trans = 1;
% v = 0.95;
if nargin<8
    triangulation = [];
end
imagenames = dir([image_folder, filesep, '*.JPG']);
imagenames = {imagenames.name};
pt_model_files = dir([pt_model_folder, filesep, '*_pm.mat']);
pt_model_files = {pt_model_files.name};
for i=1:length(pt_model_files)
   file = pt_model_files{i};
   l = load([pt_model_folder, filesep, file]);
   pointmodels(:,i) = l.pts;
end
[pointmodels,ProcrustesAlignment] = pmalign(pointmodels,[scale rotate trans]);

[pdm.Xm, pdm.P, pdm.b, pdm.pca] = pca(pointmodels, v);

options = BuildOptions(image_folder, imagenames, pdm, pt_model_folder, ProcrustesAlignment, triangulation)
[sfam, options] = buildsfam(options, resolution, v, 1);





%%%%%%%%%%%%%%%%%%%%
%
% Build Options
%
%%%%%%%%%%%%%%%%%%%%
function options = BuildOptions(imagedir, imagenames, pdm, pt_model_dir, ProcrustesAlignment, triangulation)
% We need to create the default list of options
p = pwd;
PMTDirec = '';
PMTName = '';
ModelName = '';
options.Type = 'soptions';
options.CWD = pwd;
options.InputDirectory = [imagedir, filesep];
options.InputFilenames = imagenames;
options.triangulation = triangulation;
options.InputType = 0;
options.InputStartFrame = 1;
options.InputStartMovie = 1;
options.PointModelFile = '';
options.MeshFile = '';
options.MeshMappingFile = '';
options.PDMFile = '';
options.GLDMFile = '';
options.SFAMFile = '';
options.SAMFile = '';
options.LabelMode = 0;
options.FitterType = 0;
options.OutputDirectory = '';
options.Output3DVertices = '';
options.Output3DVerticesD = '';
options.OutputShapePCA = '';
options.OutputTexturePCA = '';
options.OutputCombinedPCA = '';
options.OutputType = 0;
options.ShowLandmarkSets = 0;
options.ShowPicture = 1;
options.TryFit = 1;
options.IncludeRotate = 0;
options.GLDMType = 1;
options.FITPOINTTYPE = 1;
options.subject = '';
options.Tryfit = 0;
options.ProcrustesAlignment = ProcrustesAlignment;

options.img ='';
options.pdm = pdm;
options.sfam = '';
options.PMdirectory = [pt_model_dir, filesep];
options.dat = '';
options.indices = 1:length(pdm.Xm)/2;
options.imagesDir = imagedir;
return

function [Xm, P, b, pcaDat] = pca(X, v)
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;