function FVC2 = blend_textures_mesh(FVC, textures, values, lambda)
%BLEND_TEXTURES_MESH    Blend Textures on a Mesh.
%   FVC2 = blend_textures_mesh(FVC, textures, values, lambda) outputs a 
%   textured mesh by blending some input textures on a triangular base mesh
%   via a screened Poisson equation with least-criterion gradient selection 
%   (or null gradient for missing data).
%
%   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.
%       * textures: array of size #vertices x 3 x #textures representing the textures to blend on the mesh.
%       * values: array of size #faces x #textures representing the criterion values of the respective faces.
%       * lambda: non-negative scalar representing the regularization with respect to the base mesh texture.
%
%   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.
%
%   Copyright (c) 2013, Arnaud Dessein (University of York)

% Get the sizes
Nfaces      = size(FVC.faces, 1);
Nvertices  	= size(FVC.vertices, 1);

% Compute the triangle areas
T = 0.5 * sqrt(sum(cross(                                                   ...
    FVC.vertices(FVC.faces(:, 1), :) - FVC.vertices(FVC.faces(:, 2), :),    ...
    FVC.vertices(FVC.faces(:, 1), :) - FVC.vertices(FVC.faces(:, 3), :)     ...
    ).^2, 2));

% Initialize the sparse linear systems of equations per color channel
I	= zeros(9 * Nfaces, 1);
J	= zeros(9 * Nfaces, 1);
A	= zeros(9 * Nfaces, 1);
B1  = zeros(9 * Nfaces, 1);
B2  = zeros(9 * Nfaces, 1);
B3  = zeros(9 * Nfaces, 1);
l 	= 1;

% For each triangle
for k = 1: Nfaces
    
    % Get the vertex indices
    i1  = FVC.faces(k, 1);
    i2  = FVC.faces(k, 2);
    i3  = FVC.faces(k, 3);
    
    % Get the vertex positions
    a   = FVC.vertices(i1, :);
    b   = FVC.vertices(i2, :);
    c   = FVC.vertices(i3, :);
    
    % Compute the gradients
    ab  = a - b;
    bc  = b - c;    
    ca  = c - a;
    v1	= ab - sum(ab .* bc) / sum(bc .* bc) * bc;
    v1 	= v1 / sum(v1 .* v1);      
    v2 	= bc - sum(bc .* ca) / sum(ca .* ca) * ca;
    v2 	= v2 / sum(v2 .* v2);
    v3  = -(v1 + v2);
      
    % Get the less foreshortened observed textures
    tmp_indices     = find(all(all(~isnan(textures([i1 i2 i3], :, :)), 2), 1));
    tmpNtextures    = length(tmp_indices);
    if tmpNtextures == 0
        TMP = zeros(3);
    elseif tmpNtextures == 1
        TMP = textures([i1 i2 i3], :, tmp_indices);
    else
        tmp = values(k, tmp_indices);
        idx = find(tmp == min(tmp));
        if length(idx) ~= 1
            idx = idx(randi(length(idx)));
        end
        TMP = textures([i1 i2 i3], :, tmp_indices(idx));
    end
   
    % Update the linear systems of equations per color channel (can speed up by simplifying some redundant computations)
    
    % 1-1
    tmp     = sum(v1 .* v1) * T(k);
    I(l) 	= i1;
    J(l) 	= i1;
    A(l)	= tmp;
    B1(l)   = tmp * TMP(1, 1);
    B2(l)   = tmp * TMP(1, 2);
    B3(l)   = tmp * TMP(1, 3);
    l       = l + 1;
    
    % 2-2
    tmp     = sum(v2 .* v2) * T(k);
    I(l) 	= i2;
    J(l) 	= i2;
    A(l)	= tmp;
    B1(l)   = tmp * TMP(2, 1);
    B2(l)   = tmp * TMP(2, 2);
    B3(l)   = tmp * TMP(2, 3);
    l       = l + 1;
    
    % 3-3
    tmp     = sum(v3 .* v3) * T(k);
    I(l) 	= i3;
    J(l) 	= i3;
    A(l)	= tmp;
    B1(l)   = tmp * TMP(3, 1);
    B2(l)   = tmp * TMP(3, 2);
    B3(l)   = tmp * TMP(3, 3);
    l       = l + 1;
    
    % 1-2
    tmp     = sum(v1 .* v2) * T(k);
    I(l) 	= i1;
    J(l) 	= i2;
    A(l)	= tmp;
    B1(l)   = tmp * TMP(2, 1);
    B2(l)   = tmp * TMP(2, 2);
    B3(l)   = tmp * TMP(2, 3);
    l       = l + 1;
    I(l) 	= i2;
    J(l) 	= i1;
    A(l)	= tmp;
    B1(l)   = tmp * TMP(1, 1);
    B2(l)   = tmp * TMP(1, 2);
    B3(l)   = tmp * TMP(1, 3);
    l       = l + 1;
    
    % 1-3
    tmp     = sum(v1 .* v3) * T(k);
    I(l) 	= i1;
    J(l) 	= i3;
    A(l)	= tmp;
    B1(l)   = tmp * TMP(3, 1);
    B2(l)   = tmp * TMP(3, 2);
    B3(l)   = tmp * TMP(3, 3);
    l       = l + 1;
    I(l) 	= i3;
    J(l) 	= i1;
    A(l)	= tmp;
    B1(l)   = tmp * TMP(1, 1);
    B2(l)   = tmp * TMP(1, 2);
    B3(l)   = tmp * TMP(1, 3);
    l       = l + 1;
    
    % 2-3
    tmp     = sum(v2 .* v3) * T(k);
    I(l) 	= i2;
    J(l) 	= i3;
    A(l)	= tmp;
    B1(l)   = tmp * TMP(3, 1);
    B2(l)   = tmp * TMP(3, 2);
    B3(l)   = tmp * TMP(3, 3);
    l       = l + 1;
    I(l) 	= i3;
    J(l) 	= i2;
    A(l)	= tmp;
    B1(l)   = tmp * TMP(2, 1);
    B2(l)   = tmp * TMP(2, 2);
    B3(l)   = tmp * TMP(2, 3);
    l       = l + 1;
    
end

% Format the linear systems of equations per color channel
A           = sparse(I, J, A,  Nvertices, Nvertices);
B(:, 1)     = sum(sparse(I, J, B1,  Nvertices, Nvertices), 2);
B(:, 2)     = sum(sparse(I, J, B2,  Nvertices, Nvertices), 2);
B(:, 3)     = sum(sparse(I, J, B3,  Nvertices, Nvertices), 2);

clear B1 B2 B3 I J T TMP a b c ab bc ca i1 i2 i3 v1 v2 v3 idx k l tmp tmp_indices tmpNtextures

% Add the regularization terms
I       = find(all(~isnan(FVC.facevertexcdata), 2));
Atmp    = sparse(1: length(I), I, sqrt(lambda), length(I), Nvertices);
Btmp    = sqrt(lambda) * FVC.facevertexcdata(I, :);
A0      = cat(1, A, Atmp);
B0      = cat(1, B, Btmp);

% Solve the linear systems of equations per color channel
state = warning('query', 'MATLAB:rankDeficientMatrix');
if strcmp(state.state, 'on')
    warning('off', 'MATLAB:rankDeficientMatrix')
    X0 = A0 \ B0;
    warning('on', 'MATLAB:rankDeficientMatrix')
else
    X0 = A0 \ B0;
end

clear A B A0 B0 Atmp Btmp I Nfaces Nvertices state

% Format the output
FVC2                    = FVC;
FVC2.facevertexcdata    = max(0, min(1, X0));

clear X0

end