function [FVC2, filenames] = synthesize_face_hallu(FVC, patches, db, options)
%SYNTHESIZE_FACE_HALLU  Synthesize a Face by Hallucination.
%   [FVC2, filenames] = synthesize_face_hallu(FVC, patches, db, options) 
%   outputs the textured mesh of a synthesized face by hallucination from a 
%   patch database, where overlapping patches are selected so as to 
%   minimize the Euclidean distance to the input face with optional texture 
%   blending as post-processing.
%
%   Input arguments:
%       * FVC: structure representing a textured mesh with the following fields:
%           - faces: array of size #faces x 3 where rows contain the 3 vertex indices of the respective faces;
%           - vertices: array of size #vertices x 3 where rows contain the 3 spatial coordinates of the respective vertices;
%           - facevertexcdata: array of size #vertices x 3 where rows contain the 3 color coordinates of the respective vertices.
%       * patches: cell of length #samples (or #symmetric sample pairs) where entries are vectors (or 2-column arrays) containing the vertex indices that are included in the respective patches;
%       * db: string representing the path of the database.
%       * options: structure of options with the following fields:
%           - verbose: boolean to display the current file being read (default: false);
%           - blend: boolean to blend the selected textures (default: false);
%           - lambda: non-negative scalar representing the blending regularization with respect to the averaged textures (default: 1e-6).
%
%   Output arguments:
%       * FVC2: structure representing a textured mesh with the following fields:
%           - faces: array of size #faces x 3 where rows contain the 3 vertex indices of the respective faces;
%           - vertices: array of size #vertices x 3 where rows contain the 3 spatial coordinates of the respective vertices;
%           - facevertexcdata: array of size #vertices x 3 where rows contain the 3 color coordinates of the respective vertices.
%       * filenames: cell of size #patches where entries contain the names of the respective files selected for the different patches.
%
%   Copyright (c) 2013, Arnaud Dessein (University of York)

% Options
if nargin < 4
    options = struct;
end

% Verbose
verbose = false;
if isfield(options, 'verbose')
    verbose = options.verbose;
end

% Blend
blend = false;
if isfield(options, 'blend')
    blend = options.blend;
end

% Regularization
lambda = 1e-6;
if isfield(options, 'lambda')
    lambda = options.lambda;
end

% Get the dimensions
Nvertices   = size(FVC.vertices, 1);
Nfaces      = size(FVC.faces, 1);
Npatches    = length(patches);

% Get the filenames
files   = struct2cell(dir([db '*.mat']));
files 	= files(1, :);
Nfiles  = length(files);

% Initialization
face_errors     = zeros(Nfaces, Nfiles);
patch_errors    = zeros(Npatches, Nfiles);

% Compute the errors
for f = 1: Nfiles
    
    if verbose
        disp(['File number ' num2str(f) '.']);
    end
    
    % Load the file
    load([db files{f}]);
    
    % Compute the vertex errors
    vertex_errors = sum((FVC.facevertexcdata - texture).^2, 2);
    
    % Update the face errors
    face_errors(:, f) = sum(vertex_errors(FVC.faces), 2);
    
    % Update the patch errors
    for p = 1: Npatches
        patch_errors(p, f) = sum(vertex_errors(patches{p}(:)));
    end
    
end

% Get the minimum patch errors
[~, fs] = min(patch_errors, [], 2);

clear f p patch_errors vertex_errors texture

% Initialization
fs2         = unique(fs);
Nfiles      = length(fs2);
textures    = NaN(Nvertices, 3, Nfiles);

% Filenames
if nargout > 1
   filenames = files(fs); 
end

% Get the corresponding textures
for f = 1: Nfiles
    
    % Load the file
    load([db files{fs2(f)}]);
    
    % Update the texture
    patch_ind = find(fs == fs2(f));
    for p = 1: length(patch_ind)
        textures(patches{patch_ind(p)}, :, f) = texture(patches{patch_ind(p)}, :);
    end
    
end

clear files f p texture patch_ind

% Average the textures
test                    = isnan(textures);
textures(test)          = 0;
norm                    = Nfiles - sum(test(:, 1, :), 3);
FVC2                    = FVC;
FVC2.facevertexcdata    = sum(textures, 3) ./ repmat(norm, 1, 3);
textures(test)          = NaN;

% Merge the textures
if blend
    
    % Compute the face errors
    face_errors = face_errors(:, fs2);
    
    % Perform texture blending
    FVC2 = blend_textures_mesh(FVC2, textures, face_errors, lambda);

end

clear face_errors blend fs fs2 lambda norm textures test verbose Nfaces Nfiles Npatches Nvertices

end