function FVC2 = render_texture_mesh(FVC, oglp)
%RENDER_TEXTURE_MESH    Render Texture on a Mesh.
%   FVC2 = render_texture_mesh(FVC, oglp) outputs a rendered texture on a
%   mesh according to the specified lighting parameters.
%
%   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.
%       * oglp: structure representing the lighting parameters in OpenGL format.
%
%   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 rendering parameters
T               = oglp.T;                   % homogeneous extrinsic transformation matrix
i_amb_light     = oglp.ambient_light;   	% ambient light intensity
i_dir_light     = oglp.dir_light.intens;    % directed light intensity
d_dir_light     = oglp.dir_light.dir;       % directed light direction
shininess       = oglp.phong_exp;           % shininess
specularity     = oglp.specular / 255;      % specularity
do_cast_shadows = oglp.do_cast_shadows;     % flag for casting shadows
shadow_buf_size = oglp.sbufsize;            % buffer size for casting shadows

% Get the vertices
V = FVC.vertices;

% Build the rendered mesh
FVC2.faces      = FVC.faces;
FVC2.vertices   = FVC.vertices;

% Compute the transformed vertices
M       = T(1: 3, :);   % the output does not need to be in homogeneous coordinates
V(:, 4)	= 1;            % use homogeneous coordinates for input
V2   	= V * M.';      % the vertices are transposed

clear V M

% Get the triangle vertices
v1 = FVC.faces(:, 1);
v2 = FVC.faces(:, 2);
v3 = FVC.faces(:, 3);

% Compute the edge vectors
e1s = V2(v2, :) - V2(v1, :);
e2s = V2(v3, :) - V2(v1, :);
e3s = V2(v2, :) - V2(v3, :);

% Normalize the edge vectors
e1s_norm = e1s ./ repmat(sqrt(sum(e1s.^2, 2)), 1, 3);
e2s_norm = e2s ./ repmat(sqrt(sum(e2s.^2, 2)), 1, 3);
e3s_norm = e3s ./ repmat(sqrt(sum(e3s.^2, 2)), 1, 3);

% Compute the angles
angles(:, 1) = acos(sum(e1s_norm .* e2s_norm, 2));
angles(:, 2) = acos(sum(e3s_norm .* e1s_norm, 2));
angles(:, 3) = pi - (angles(:, 1) + angles(:, 2));

% Compute the triangle weighted normals
triangle_normals    = cross(e1s, e3s, 2);
w1_triangle_normals = triangle_normals .* repmat(angles(:, 1), 1, 3);
w2_triangle_normals = triangle_normals .* repmat(angles(:, 2), 1, 3);
w3_triangle_normals = triangle_normals .* repmat(angles(:, 3), 1, 3);

clear e1s e2s e3s e1s_norm e2s_norm e3s_norm angles triangle_normals

% Initialize the vertex normals
normals = zeros(size(FVC.vertices, 1), 3);

% Determine the visibility from the directed light source
if do_cast_shadows && isposintscalar(shadow_buf_size)
  
    % Compute the rotation from the viewer to the directed light source
    u = cross(d_dir_light, [0; 0; 1]) ; % unnormalized axis
    u = u / norm(u, 2);                 % normalized axis
    t = acos([0 0 1] * d_dir_light);    % angle
    R = [   u(1) * u(1) * (1 - cos(t)) +        cos(t), u(2) * u(1) * (1 - cos(t)) - u(3) * sin(t), u(3) * u(1) * (1 - cos(t)) + u(2) * sin(t)
            u(1) * u(2) * (1 - cos(t)) + u(3) * sin(t), u(2) * u(2) * (1 - cos(t)) +        cos(t), u(3) * u(2) * (1 - cos(t)) - u(1) * sin(t)
            u(1) * u(3) * (1 - cos(t)) - u(2) * sin(t), u(2) * u(3) * (1 - cos(t)) + u(1) * sin(t), u(3) * u(3) * (1 - cos(t)) +        cos(t)
        ];
    
    % Compute the transformed vertices
    V3 = V2 * R.'; % the vertices are transposed
    
    % Store the vertex depths for z-buffering
    Z = V3(:, 3);
    
    % Compute the projected vertices in the image plane
    width       = shadow_buf_size;
    height      = shadow_buf_size;
    UV(:, 1)    = (V3(:, 1) - min(V3(:, 1))) / (max(V3(:, 1)) - min(V3(:, 1))) * width;
    UV(:, 2)    = (V3(:, 2) - min(V3(:, 2))) / (max(V3(:, 2)) - min(V3(:, 2))) * height;
    
    clear V2 V3 R u t
    
    % Transform to the pixel plane
    UV = UV + 0.5;
    
    % Construct the pixel grid
    [rows, cols] = meshgrid(0: width + 1, 0: height + 1); % pad to avoid boundary problems when interpolating
    
    % Compute bounding boxes for the projected triangles
    x       = [UV(v1, 1) UV(v2, 1) UV(v3, 1)];
    y       = [UV(v1, 2) UV(v2, 2) UV(v3, 2)];
    minx    = max(0,            ceil (min(x, [], 2)));
    maxx    = min(width + 1,    floor(max(x, [], 2)));
    miny    = max(0,            ceil (min(y, [], 2)));
    maxy    = min(height + 1,   floor(max(y, [], 2)));
    
    clear x y
    
    % Initialize the z-buffer
    zbuffer = -inf(height + 2, width + 2); % pad to avoid boundary problems when interpolating
    
    % For each triangle (can speed up by comparing the triangle depths to the z-buffer and priorly sorting the triangles by increasing depth)
    for i = 1: size(FVC.faces, 1)
        
        % Update the vertex normals
        normals(v1(i), :) = normals(v1(i), :) + w1_triangle_normals(i, :);
        normals(v2(i), :) = normals(v2(i), :) + w2_triangle_normals(i, :);
        normals(v3(i), :) = normals(v3(i), :) + w3_triangle_normals(i, :);
        
        % If some pixels lie in the bounding box
        if minx(i) <= maxx(i) && miny(i) <= maxy(i)
            
            % Get the pixels lying in the bounding box
            px = rows(miny(i) + 1: maxy(i) + 1, minx(i) + 1: maxx(i) + 1);
            py = cols(miny(i) + 1: maxy(i) + 1, minx(i) + 1: maxx(i) + 1);
            px = px(:);
            py = py(:);
            
            % Compute the edge vectors
            e0 = UV(v1(i), :);
            e1 = UV(v2(i), :) - e0;
            e2 = UV(v3(i), :) - e0;
            
            % Compute the barycentric coordinates (can speed up by first computing and testing a solely)
            det     = e1(1) * e2(2) - e1(2) * e2(1);
            tmpx    = px - e0(1);
            tmpy    = py - e0(2);
            a       = (tmpx * e2(2) - tmpy * e2(1)) / det;
            b       = (tmpy * e1(1) - tmpx * e1(2)) / det;
            
            % Test whether the pixels lie in the triangle
            test = a >= 0 & b >= 0 & a + b <= 1;
            
            % If some pixels lie in the triangle
            if any(test)
                
                % Get the pixels lying in the triangle
                px = px(test);
                py = py(test);
                
                % Interpolate the triangle depth for each pixel
                w2 = a(test);
                w3 = b(test);
                w1 = 1 - w2 - w3;
                pz = Z(v1(i)) * w1 + Z(v2(i)) * w2 + Z(v3(i)) * w3;
                
                % Update the z-buffer
                for j = 1: length(pz)
                    zbuffer(py(j) + 1, px(j) + 1) = max(zbuffer(py(j) + 1, px(j) + 1), pz(j));
                end
                
            end
            
        end
        
    end
    
    clear det e0 e1 e2 i j minx maxx miny maxy px py pz test a b w1 w2 w3 tmpx tmpy
    
    % Normalize the vertex normals
    normals = normals ./ repmat(sqrt(sum(normals.^2, 2)), 1, 3);
    
    clear v1 v2 v3 w1_triangle_normals w2_triangle_normals w3_triangle_normals
    
    % Compute the dot products between the vertex normals and the directed light source direction
    NdotL = d_dir_light(1) * normals(:, 1) + d_dir_light(2) * normals(:, 2) + d_dir_light(3) * normals(:, 3);
    
    % Test attached shadows
    test = NdotL >= 0;
    
    % Interpolate the z-buffer at the projected vertices
    Z2          = inf(size(FVC.vertices, 1), 1);
    Z2(test)    = interp2(rows, cols, zbuffer, UV(test, 1), UV(test, 2), '*linear');
    
    % Fix the projected outline
    tmp_test        = Z2 == -inf;
    Z2(tmp_test)	= interp2(rows, cols, zbuffer, UV(tmp_test, 1), UV(tmp_test, 2), '*nearest');
    
    clear cols rows
    
    % Fix the projected outline
    tmp_ind         = find(Z2 == -inf);
    [i j]           = find(isfinite(zbuffer));
    i               = i - 1;
    j               = j - 1;
    if exist('knnsearch', 'file')
        tmp_idx = knnsearch([j i], UV(tmp_ind, :),          ...
                            'K',            1,              ...
                            'NSMethod',     'kdtree',       ...
                            'IncludeTies',  true,           ...
                          	'Distance',     'euclidean',    ...
                            'BucketSize',   50              ...
                            );
        for k = 1: length(tmp_ind)
            Z2(tmp_ind(k)) = max(zbuffer(i(tmp_idx{k}) + 1 + j(tmp_idx{k}) * (height + 2)));
        end
    else
        for k = 1: length(tmp_ind)
            d2              = (UV(tmp_ind(k), 1) - j).^2 + (UV(tmp_ind(k), 2) - i).^2;
            tmp_idx         = find(d2 == min(d2));
            Z2(tmp_ind(k))  = max(zbuffer(i(tmp_idx) + 1 + j(tmp_idx) * (height + 2)));
        end
    end
    
    % Test cast shadows
    test(test) = abs(Z2(test) - Z(test)) <= 0.02 * (max(Z) - min(Z));
    
    clear UV Z2 Z zbuffer tmp_idx tmp_ind tmp_test i j k d2 width height
    
else
    
    % Update the vertex normals
    for i = 1: size(FVC.faces, 1)
        normals(v1(i), :) = normals(v1(i), :) + w1_triangle_normals(i, :);
        normals(v2(i), :) = normals(v2(i), :) + w2_triangle_normals(i, :);
        normals(v3(i), :) = normals(v3(i), :) + w3_triangle_normals(i, :);
    end
    
    % Normalize the vertex normals
    normals = normals ./ repmat(sqrt(sum(normals.^2, 2)), 1, 3);
    
    clear v1 v2 v3 w1_triangle_normals w2_triangle_normals w3_triangle_normals i V2
    
    % Compute the dot products between the vertex normals and the directed light source direction
    NdotL = d_dir_light(1) * normals(:, 1) + d_dir_light(2) * normals(:, 2) + d_dir_light(3) * normals(:, 3);
    
    % Test attached shadows
    test = NdotL >= 0;
    
end

% Compute the halfway vector between the viewer and the directed light source directions (both assumed at infinity)
H = [0; 0; 1] + d_dir_light;
H = H / norm(H, 2);

% Compute the dot products between the vertex normals and the halfway vector
NdotH = H(1) * normals(test, 1) + H(2) * normals(test, 2) + H(3) * normals(test, 3);

% Fix the dot products
NdotH = max(0, NdotH);

% Get the texture
texture1 = FVC.facevertexcdata(:, 1);
texture2 = FVC.facevertexcdata(:, 2);
texture3 = FVC.facevertexcdata(:, 3);

% Compute the ambient component
ambient1 = i_amb_light(1) * texture1;
ambient2 = i_amb_light(2) * texture2;
ambient3 = i_amb_light(3) * texture3;

% Compute the diffuse component
diffuse1 = i_dir_light(1) * texture1(test) .* NdotL(test);
diffuse2 = i_dir_light(2) * texture2(test) .* NdotL(test);
diffuse3 = i_dir_light(3) * texture3(test) .* NdotL(test);

% Compute the specular component
specular    = specularity * NdotH.^shininess;
specular1   = i_dir_light(1) * specular;
specular2   = i_dir_light(2) * specular;
specular3   = i_dir_light(3) * specular;

clear specular normals NdotH NdotL H

% Render the texture
texture1        = ambient1;
texture2        = ambient2;
texture3        = ambient3;
texture1(test)  = texture1(test) + diffuse1 + specular1;
texture2(test)  = texture2(test) + diffuse2 + specular2;
texture3(test)  = texture3(test) + diffuse3 + specular3;

clear test ambient1 diffuse1 specular1 ambient2 diffuse2 specular2 ambient3 diffuse3 specular3

% Build the texture
FVC2.facevertexcdata = [texture1, texture2, texture3];

clear texture1 texture2 texture3

% Clamp the texture
FVC2.facevertexcdata = max(0, min(1, FVC2.facevertexcdata));

clear T i_amb_light i_dir_light d_dir_light shininess specularity do_cast_shadows shadow_buf_size

end