function FVC = stitch_scans_mesh(FV, scans, patches, options)
%STITCH_SCANS_MESH  Stitch Scans on a Mesh.
%   FVC = stitch_scans_mesh(FV, scans, patches, options) outputs a textured
%   mesh by stitching some input scans on a base shape mesh via texture 
%   averaging as well as optional least-angle selection of vertices or 
%   patches as pre-processing and texture blending as post-processing.
%
%   Input arguments:
%       * FV: structure representing a 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.
%       * scans: cell of length #scans where entries are structures containing scans to stich on the base shape mesh with the following fields:
%           - facevertexcdata: array of size #vertices x 3 where rows contain the 3 color coordinates of the respective vertices.
%           - facenormals: array of size #faces x 3 where rows contain the 3 spatial coordinates of the respective face normals;
%           - vertexnormals: array of size #vertices x 3 where rows contain the 3 spatial coordinates of the respective vertex normals.
%       * patches: cell of length #patches where entries are vectors containing the indices of vertices that are included in the respective patches (can be left empty for vertex- instead of patch-based stitching).
%       * options: structure of options with the following fields:
%           - selec: boolean to perform least-angle selection of vertices or patches (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:
%       * 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.
%
%   Copyright (c) 2013, Arnaud Dessein (University of York)

% Patches
if nargin < 3
    patches = {};
end
Npatches = length(patches);

% Options
if nargin < 4
    options = struct;
end

% Selection
selec = false;
if isfield(options, 'selec')
    selec = options.selec;
end

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

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

% Initialization
Nvertices   = size(FV.vertices, 1);
Nfaces      = size(FV.faces, 1);
Nscans      = length(scans);
textures    = NaN(Nvertices, 3, Nscans);

% Get the textures
for i = 1: Nscans
    textures(:, :, i) = scans{i}.facevertexcdata;
end

% Least-angle selection
if selec
    
    % Find the observed vertices
    test = ~squeeze(any(isnan(textures), 2));
    
    % Compute the angles
    angles = 0.5 * pi * ones(Nvertices, Nscans);
    for i = 1: Nscans
        angles(test(:, i), i) = acos(scans{i}.vertexnormals(test(:, i), 3));
    end
    angles = pi - angles;
    
    % Initialize the fixed textures
    textures2 = NaN(Nvertices, 3, Nscans);
    
    % Vertex-based
    if isempty(patches)
        
        % Find the scan with least-angle per vertex
        [~, jmin] = min(angles, [], 2);
        
        % Update the fixed textures
        for j = 1: Nscans
            test                    = find(jmin == j);
            textures2(test, :, j)	= textures(test, :, j);
        end
        
        % Patch-based
    else
        
        % For each patch
        for i = 1: Npatches
            
            % Sort the scans by increasing angle per vertex
            [~, jmin] = sort(mean(angles(patches{i}, :), 1), 'ascend');
            
            % Initialize the vertices to fill
            test = patches{i};
            
            % For each scan
            for j = 1: Nscans
                
                % Fill in the vertices with observed data from the current scan
                texture                     = textures(test, :, jmin(j));
                textures2(test, :, jmin(j))	= texture;
                
                % Update the vertices to fill in
                test = test(any(isnan(texture), 2));
                
                % Check whether filling is finished
                if isempty(test)
                    break
                end
                
            end
            
        end
        
    end

% No selection
else
    
    % Copy the textures
    textures2 = textures;
    
end

clear textures texture test angles jmin i j

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

% Merge the textures
if blend
    
    % Compute the face angles if needed
    if selec && isempty(patches)
        angles = [];
    else
        angles = zeros(Nfaces, Nscans);
        for i = 1: Nscans
            angles(:, i) = acos(scans{i}.facenormals(:, 3));
        end
        angles = pi - angles;
    end
    
    % Perform texture blending
    FVC = blend_textures_mesh(FVC, textures2, angles, lambda);
    
end

clear textures2 angles test norm blend selec lambda i Nfaces Nvertices Nscans Npatches

end