function [sfam, options] = buildsfam(options, npix, v, critical_triangle_expansion, template_name, verb, boundary_index)
%BUILDSFAM builds a shape-free appearance model
%  [SFAM, OPTIONS] = BUILDSFAM(OPTIONS,NPIX,V, VERB) builds a
%  shape-free model of appearance given a labelled set of images.
%  OPTIONS is a set of options given by SOPTIONS and NPIX is the
%  (approximate) number of pixels to include in the model.
%
%  V defines the percentage variance the model should retain
%  (default 0.95 corresponding to 95%).
%
%  VERB indicates whether the warped images are displayed.
%
%  NOTE: At present, there is no image normalisation - this
%  still need to be implemented! This is a cutdown version of
%  sfaload.m from Matthews appearance model toolbox, adapted to use
%  the options from soptions, and handle movies as input. Could later
%  write a script to convert the options filenames to a filespec and
%  use the appearance code - must then ensure the directories are
%  correct etc!
%
%
% Barry-John Theobald, b.theobald@uea.ac.uk
% School of Information Systems, University of East Anglia,
% Norwich, NR4 7TJ, UK
%
% Version 1.2

% HISTORY:
%     01/01/2003 - 1.2 - Optionally return the options structure and set
%                        values such as filenames and paths etc.
%     31/12/2002 - 1.1 - Added field ftype, with default value 'l' to
%                        ensure compatability with AAM toolbox!
%     31/12/2002 - 1.0 - Initial revision



% clear mex;
% Check args
if(nargin<7) boundary_index = []; end
if(nargin < 5) verb = 1; end
if(nargin < 4) critical_triangle_expansion = 5; end
if(nargin < 3) v=0.85; end; %v = 0.95; end
if(nargin < 2) npix = -1; end
if(nargin < 1) options = ''; end


% Get, or check the options are valid!
[options, msg] = get_options('soptions', options);
if(isempty(options))
    error(msg);
    return;
end

% Get the PDM defined in the options
if isempty(options.pdm)
    [pdm, msg] = get_pdm(options.PDMFile);
else
    pdm=options.pdm;
end
if(isempty(pdm))
    error(msg);
    return;
end


% SFAM is only defined for 2D Landmarks!
if(options.OutputType ~= 0)
    error('Can only build an SFAM from 2D landmarks');
end

% Check if OPTS has an output directory defined - need to know
% where the input landmarks are lcoated.
%if(isempty(options.OutputDirectory))
if(isempty(options.PMdirectory))
    [fname, pname] = uigetfile('*_pm.mat', 'Locate landmarks');
    if(~any(fname))
        return;
    else
        options.PMdirectory  = pname;
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialise

% Get the mean points
mpts = reshape(pdm.Xm, 2, length(pdm.Xm)/2)';

% Get triangulation of the mean points
tri = delaunay(mpts(:,1), mpts(:,2),{'QJ','QJ','QJ'});
tri_filename=fullfile(options.imagesDir,'triangulation.dat');
save(tri_filename,'tri','-mat');
tri = tri';
tri = tri(:);


%check on the filenames
[filenames,imagelist]=CheckPMexist(options.PMdirectory,options.InputDirectory,options.InputFilenames, options.imagesDir);

% How big is the mean shape image?
msiz = round(max(mpts) - min(mpts));
if any(msiz<10)
    fprintf('Image is very small, might cause trouble, buildsfam line 103');
end
switch(options.InputType)
    case 0
        %img = imread(strcat(options.InputDirectory, options.InputFilenames(1,:)));
        img = imread(strcat(options.InputDirectory, imagelist(1,:)));
    case {1,2}
        options = soptions(options, 'InputFilenames', '');
        if(options.InputType == 1)
            mov = aviread(strcat(options.InputDirectory, options.InputFilenames(1,:)));
        else
            [mov, msg] = get_movie(strcat(options.InputDirectory, options.InputFilenames(1,:)));
            if(isempty(mov)) error(msg); end
        end
        img = mov(1).cdata;
end

if(length(size(img)) == 3)
    msiz = [msiz 3];
    npix = npix * 3;
end

% Scale the model to contain npix pixels of interest
if(npix > 0)
    % Ensure mean model is zero centred
    mpts = mpts - ((sum(mpts) / size(mpts,1))' * ones(1,size(mpts,1)))';

    % Offset to (0,0) min for warping
    ompts = mpts - (min(mpts)' * ones(1,size(mpts,1)))';

    % Find out how many pixels of interest are currently returned
    timg = triwarp(ones(msiz)*100, ompts, ompts);
    idx = find(timg > 0);

    % Apply scaling
    scale = sqrt(npix / length(idx));
    mpts = applytrm(reshape(mpts', length(pdm.Xm), 1), [scale; 0; 0; 0]);
    mpts = reshape(mpts, 2, length(pdm.Xm)/2)';

end

% Offset to (0,0) min for warping destination
ompts = mpts - (min(mpts)' * ones(1,size(mpts,1)))';

% How big are the images? Mean pts dimensions. Might be colour.
npix = prod(round(max(ompts)));
if(length(size(img)) == 3)
    npix = npix * 3;
end

% Find indicies of pixels of interest in mean shape.
msiz = round(max(ompts));
if(length(size(img)) == 3)
    msiz = [msiz 3];
end


timg = triwarp(ones(msiz)*100, ompts, ompts);
omsidx = find(timg > 10);


% Display?
% if(verb)
% 	fH = figure;
% 	iH = image(uint8(zeros(size(timg))));
% 	set(gcf,'DoubleBuffer','on');
% 	truesize
% 	omsidx = find(timg > 0);
% end


modelpath=options.PDMFile(1:findstr(options.PDMFile,'mod_pdm.mat')-2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gather the shape-free textures
if verb
    figno=figure;
end
fidx = 1;
switch options.InputType
    case 0 % Image sequences

        % Initialise
        %sfa = zeros(npix, size(options.InputFilenames,1));
        sfa = zeros(npix, size(imagelist,1));
        sfa_NaN=zeros(npix, size(imagelist,1));
        % Check each image in turn
        %for k = 1:size(options.InputFilenames)
        expected_number=size(imagelist,1);
        for k = 1:size(imagelist,1)

            % Is the current image labelled?
            %baseName = options.InputFilenames(k,:);
            baseName = imagelist(k,:);
            imgName = strcat(options.InputDirectory, baseName);
            dot = max(findstr(baseName, '.'));
            %ptsName = strcat(options.OutputDirectory, baseName(1:dot-1), '_pm.mat');
            ptsName = strcat(options.PMdirectory, baseName(1:dot-1), '_pm.mat');
            image_baseName=baseName(1:dot-1);

            if(exist(ptsName, 'file'))

                % Collect the image and landmarks
                img = imread(imgName);
                %eval(sprintf('load %s', ptsName));
                load(ptsName);
                template = load(['Templates', filesep, template_name, '.temp_dat'], '-mat');
                template = template.PMTemplate;
                loops = get(template, 'loops');
                edges=loops{1};
                %             imagesDir=options.imagesDir;
                %             modelname=fullfile(imagesDir,'mod_sfam');
                %             temp=load(modelname);
                allTRI=[];
                ptsx=pts(1:2:end);
                ptsy=pts(2:2:end);
                if ~isempty(boundary_index)
                    if iscell(boundary_index)
                        mask  = zeros(size(img,1), size(img,2));
                        for i=1:length(boundary_index)
                            mask_x = ptsx(boundary_index{i});
                            mask_y = ptsy(boundary_index{i});
                            b = roipoly(img, mask_x, mask_y);
                            mask = mask | b;
                        end
                        mask = uint8(mask);
                        for i=1:size(img,3)
                            img(:,:,i) = img(:,:,i).*mask;
                        end

                    end
                end
                selected_pts=[ptsx(options.indices),ptsy(options.indices)]';

                %[img2,sfimg_safe,img1_temp,allTRI_temp,upts,uompts]=warp_to_mean_shape(pts,ompts,img,options.indices, critical_triangle_expansion,edges,[msiz(2),msiz(1),msiz(3)],allTRI);
                [img2,sfimg_safe,img1_temp,allTRI_temp,upts,uompts]=warp_to_mean_shape(selected_pts(:),ompts,img,options.indices, critical_triangle_expansion,edges,[msiz(2),msiz(1),msiz(3)],allTRI);

                %sfimg=uint8(img2); % converts NaN into zero % IGNORE THIS FOR
                %THE MOMENT
                %subplot(2,1,1);imshow(uint8(img1))
                %subplot(2,1,2);
                if verb
                    imshow(uint8(img2));
                    axis image ij;
                    if ~isempty(allTRI_temp)
                        hold on
                        %triplot(allTRI_temp,ompts(:,1),ompts(:,2),'w');
                        triplot(allTRI_temp,upts(:,1),upts(:,2),'w');
                        hold off
                    end
                    %subplot(2,1,1);imshow(uint8(img2))
                    title('Image Being Warped')
                    drawnow
                end

                sfa(:, fidx) = sfimg_safe(:);
                sfa_NaN(:, fidx) = img2(:);
                fidx = fidx + 1;

            end

        end
        if fidx<expected_number
            % Remove entries that did not exist
            keyboard
            sfa = sfa(:, 1:fidx-1);
        end


    case {1,2}

        for k = 1:size(options.InputFilenames)

            % Get the current movie
            if(options.InputType == 1)
                mov = aviread(strcat(options.InputDirectory, options.InputFilenames(k,:)));
            else
                [mov, msg] = get_movie(strcat(options.InputDirectory, options.InputFilenames(k,:)));
                if(isempty(mov)) error(msg); end
            end

            baseName = deblank(options.InputFilenames(k,:));

            % Now deal with frames as images!
            for l = 1:length(mov)

                % Is the current image labelled?
                %ptsName = strcat(options.OutputDirectory, ...
                ptsName = strcat(options.PMdirectory, ...
                    baseName(1:end-4), sprintf('_%.3d', l), '_pm.mat');

                if(exist(ptsName, 'file'))

                    % Collect the image and landmarks
                    eval(sprintf('load %s', ptsName));
                    pts = reshape(pts, 2, length(pts)/2)';

                    sfimg = double(triwarp(mov(l).cdata, pts, ompts));

                    % Display?
                    if(verb)
                        set(iH, 'CData', uint8(sfimg))
                        drawnow
                    end

                    sfa(:, fidx) = sfimg(:);
                    fidx = fidx + 1;

                end

            end

        end

end
% if verb	close; end
% drawnow
clear idx
clear img
% Calculate SFAM
clear gltriwarp % release the warp code space
% first we need to find the mean appearance ignoring all NaN's

[Am, P, b, pcaDat] = pca(sfa_NaN, v,sfa);
%[Am, P, b, pcaDat] = pca(sfa, v);
if verb
    clf
    imshow(uint8(reshape(Am,[msiz(2),msiz(1),msiz(3)])))
    drawnow
end;
sfam.Am = Am;
sfam.P = P;
sfam.b = b;

% Extra PCA info - useful to edit number of modes later
sfam.pca = pcaDat;


% Copy across all info
sfam.siz = msiz;
sfam.boundary_index = boundary_index;
sfam.tri = tri;
sfam.ompts = ompts;
sfam.omsidx = omsidx;
sfam.gs = 0;
sfam.apoff = 0;
sfam.apscl = 1;
sfam.ftype = 'l';
sfam.nrm = 'none';
sfam.var = v;
sfam.pdm = pdm;
% Save or return the SFAM
if(nargout < 1)
    % Save the model accoring to the options
    if(isempty(options.SFAMFile))
        [fname, pname] = uiputfile('*_sfam.mat', 'Save SFAM As');
        if(~any(fname)) return; end

        % Check the file extension is _sfam.mat
        dot = max(findstr(fname, '.'));
        if(~dot)
            fname = strcat(fname, '_sfam.mat');
        else
            if(~strcmp(fname(dot-4:end), '_sfam.mat'))
                fname = strcat(fname(1:dot-1), '_sfam.mat');
            end
        end
    end

    options.SFAMFile = strcat(pname, fname);

    eval(sprintf('save %s sfam', opts.SFAMFile));
    return;
end
pause(1)
if verb
    close(figno)
end


function [Xm, P, b, pcaDat] = pca(X, v, V)
%function [Xm, P, b, pcaDat] = pca(X, v, Xsafe)
%PCA compute a princpal component analysis
%  [XM, P, B, pcaDat] = PCA(X, V) performs a principal components
%  analysis on the column-wise data set X. The output arguments are
%
%X   - one column per image, one row per pixel
%v   - verbose (default 0)
%
%   Xm     - the mean observation
%   P      - first t modes accounting for V percent of the total variation
%   b      - the corresponding eigenvalues
%   pcaDat - all eigenvectors, eigenvalues and the cumulative percentage
%            variance captured by the successive modes, also the covariance.
%
%
% Barry-John Theobald, b.theobald@uea.ac.uk
% School of Information Systems, University of East Anglia,
% Norwich, NR4 7TJ, UK
%
% Version 1.6

% HISTORY:
%     02/12/2002 - 1.6 - No longer return the inverse of the covariance matrix, the
%                        Mahalonabis distance is approximated in the fitting routines
%                        rather than computed using this. (Also removed ffr);
%     02/12/2002 - 1.5 - Have added the fast method back in, but also added an optional
%                        flag that allows the normal method to be used even if there are
%                        not enough examples.
%     01/12/2002 - 1.4 - Removed fast method - using inv(S) for fast method, the
%                        dimensions are incorrect when computing the Mahalonabis distance
%     01/12/2002 - 1.3 - Changed return variable S to inv(S), so it doesn't need
%                        inverted at each stage of the fit.
%     29/11/2002 - 1.2 - Returned covariance matrix in pcaDat, useful if Mahalonabis
%                        is used by the fitter.
%     27/11/2002 - 1.1 - Return PCA data as individual elements rather than structs.
%                        Different models use different variable names (e.g. Xm or Gm).
%     21/11/2002 - 1.0 - Initial revision



% Check input args
if(nargin < 2) v = 0.95; end
if(nargin < 1) error('No data defined for PCA'); end


% Transpose into row-wise data for covaraince
X = X';

[m, n] = size(X);

% Find mean shape and the covariance matrix
% the computationally safe approach is not ignore the effect of expanding
% triangles by large amounts
Xmsafe = mean(V',1);
% but this could be experimentally unsafe so try and replace all such
% triangles (identified by NaN) by the mean of remaining images
d=isnan(X); %NAN=isnan(X);
ind=find(d); %ind=find(NAN);
X(ind)=0; % remove all the NaN's from the sum
% and find the mean
Xm=sum(X)./sum(~d); % one row per image Xm=sum(X)./sum(~NAN); % one row per image
% however, if there is a misslabelling there may be triangles that expand
% too much that are common to all the data images, so allow the computation
% to continue by substituting the un-flagged images
ind=find(isnan(Xm(:)));
if ~isempty(ind)
    uiwait(msgbox(sprintf(...
        'So many triangles are expanding by more than the set threshold that it is necessary to use uncorrected data. There may be an inconsistent use of landmarks.')));
    Xm(ind)=Xmsafe(ind); % safe only in the computational sense, unsafe as an experiment
end
% now substitute the mean triangle values back into the original triangle
% so replacing the NaN by the means that will not influence the PCA
for i=1:size(X,1)
    ind=find(d(i,:));%ind=find(NAN(i,:));
    X(i,ind)=Xm(ind);
end
d=[];
Xmsafe=[];
% HOWEVER, there might still be some NaN's because all images might a
% coincident erroneous triangle
% In which case we just have to re-insert the original data
ind=find(isnan(X(:)));
if ~isempty(ind)
    uiwait(msgbox(sprintf(...
        'So many triangles are expanding by more than the set threshold that it is necessary to use uncorrected data. There may be an inconsistent use of landmarks.')));
    X(ind)=V(ind);
end
% free up some memory
%clear Xsafe Xmsafe
%clear NAN
% try % if memory allows
%     d = X - repmat(Xm, m, 1);
% catch

%mu = zeros(length(Xm),1);
%S = zeros(length(Xm));

for i=1:size(X,1)

    d(i,:) = X(i,:) - Xm;
    %[S, mu] = recursive_covariance(d(i,:)', S, mu, i);

end

% % end
% S2 = S;
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
            V(:,k) = V(:,k) / sqrt(V(:,k)' * V(:,k));
        end
    end
end
clear d

evecs = V(:,idx);

V=[];

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;



