blendshape_fitting.hpp

/*
 * eos - A 3D Morphable Model fitting library written in modern C++11/14.
 *
 * File: include/eos/fitting/blendshape_fitting.hpp
 *
 * Copyright 2015 Patrik Huber
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
 
#ifndef BLENDSHAPEFITTING_HPP_
#define BLENDSHAPEFITTING_HPP_
 
#include "eos/morphablemodel/Blendshape.hpp"
 
#include "Eigen/Core"
#include "Eigen/QR"
#include "Eigen/Sparse"
#include "nnls.h"
 
#include <vector>
#include <cassert>
 
namespace eos {
namespace fitting {
 
inline std::vector<float> fit_blendshapes_to_landmarks_linear(
    const std::vector<morphablemodel::Blendshape>& blendshapes, const Eigen::VectorXf& face_instance,
    Eigen::Matrix<float, 3, 4> affine_camera_matrix, const std::vector<Eigen::Vector2f>& landmarks,
    const std::vector<int>& vertex_ids, float lambda = 500.0f)
{
    assert(landmarks.size() == vertex_ids.size());
 
    using Eigen::MatrixXf;
    using Eigen::VectorXf;
 
    const auto num_blendshapes = blendshapes.size();
    const int num_landmarks = static_cast<int>(landmarks.size());
 
    // Copy all blendshapes into a "basis" matrix with each blendshape being a column:
    const MatrixXf blendshapes_as_basis = morphablemodel::to_matrix(blendshapes);
 
    // $\hat{V} \in R^{3N\times m-1}$, subselect the rows of the eigenvector matrix $V$ associated with the $N$ feature points
    // And we insert a row of zeros after every third row, resulting in matrix $\hat{V}_h \in R^{4N\times m-1}$:
    MatrixXf V_hat_h = MatrixXf::Zero(4 * num_landmarks, num_blendshapes);
    int row_index = 0;
    for (int i = 0; i < num_landmarks; ++i)
    {
        V_hat_h.block(row_index, 0, 3, V_hat_h.cols()) =
            blendshapes_as_basis.block(vertex_ids[i] * 3, 0, 3, blendshapes_as_basis.cols());
        row_index += 4; // replace 3 rows and skip the 4th one, it has all zeros
    }
    // Form a block diagonal matrix $P \in R^{3N\times 4N}$ in which the camera matrix C (P_Affine, affine_camera_matrix) is placed on the diagonal:
    MatrixXf P = MatrixXf::Zero(3 * num_landmarks, 4 * num_landmarks);
    for (int i = 0; i < num_landmarks; ++i)
    {
        P.block<3, 4>(3 * i, 4 * i) = affine_camera_matrix;
    }
 
    // The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
    VectorXf y = VectorXf::Ones(3 * num_landmarks);
    for (int i = 0; i < num_landmarks; ++i)
    {
        y(3 * i) = landmarks[i][0];
        y((3 * i) + 1) = landmarks[i][1];
        // y((3 * i) + 2) = 1; // already 1, stays (homogeneous coordinate)
    }
    // The mean, with an added homogeneous coordinate (x_1, y_1, z_1, 1, x_2, ...)^t
    VectorXf v_bar = VectorXf::Ones(4 * num_landmarks);
    for (int i = 0; i < num_landmarks; ++i)
    {
        v_bar(4 * i) = face_instance(vertex_ids[i] * 3);
        v_bar((4 * i) + 1) = face_instance(vertex_ids[i] * 3 + 1);
        v_bar((4 * i) + 2) = face_instance(vertex_ids[i] * 3 + 2);
        // v_bar((4 * i) + 3) = 1; // already 1, stays (homogeneous coordinate)
    }
 
    // Bring into standard regularised quadratic form:
    const MatrixXf A = P * V_hat_h;   // camera matrix times the basis
    const MatrixXf b = P * v_bar - y; // camera matrix times the mean, minus the landmarks
 
    const MatrixXf AtAReg =
        A.transpose() * A + lambda * Eigen::MatrixXf::Identity(num_blendshapes, num_blendshapes);
    const MatrixXf rhs = -A.transpose() * b;
 
    const VectorXf coefficients = AtAReg.colPivHouseholderQr().solve(rhs);
 
    return std::vector<float>(coefficients.data(), coefficients.data() + coefficients.size());
};
 
inline std::vector<float> fit_blendshapes_to_landmarks_nnls(
    const std::vector<morphablemodel::Blendshape>& blendshapes, const Eigen::VectorXf& face_instance,
    Eigen::Matrix<float, 3, 4> affine_camera_matrix, const std::vector<Eigen::Vector2f>& landmarks,
    const std::vector<int>& vertex_ids)
{
    assert(landmarks.size() == vertex_ids.size());
 
    using Eigen::MatrixXf;
    using Eigen::VectorXf;
 
    const auto num_blendshapes = blendshapes.size();
    const int num_landmarks = static_cast<int>(landmarks.size());
 
    // Copy all blendshapes into a "basis" matrix with each blendshape being a column:
    const MatrixXf blendshapes_as_basis = morphablemodel::to_matrix(blendshapes);
 
    // $\hat{V} \in R^{3N\times m-1}$, subselect the rows of the eigenvector matrix $V$ associated with the $N$ feature points
    // And we insert a row of zeros after every third row, resulting in matrix $\hat{V}_h \in R^{4N\times m-1}$:
    MatrixXf V_hat_h = MatrixXf::Zero(4 * num_landmarks, num_blendshapes);
    int row_index = 0;
    for (int i = 0; i < num_landmarks; ++i)
    {
        V_hat_h.block(row_index, 0, 3, V_hat_h.cols()) =
            blendshapes_as_basis.block(vertex_ids[i] * 3, 0, 3, blendshapes_as_basis.cols());
        row_index += 4; // replace 3 rows and skip the 4th one, it has all zeros
    }
 
    // Form a block diagonal matrix $P \in R^{3N\times 4N}$ in which the camera matrix C (P_Affine, affine_camera_matrix) is placed on the diagonal:
    Eigen::SparseMatrix<float> P(3 * num_landmarks, 4 * num_landmarks);
    std::vector<Eigen::Triplet<float>> P_coefficients; // list of non-zeros coefficients
    for (int i = 0; i < num_landmarks; ++i)
    { // Note: could make this the inner-most loop.
        for (int x = 0; x < affine_camera_matrix.rows(); ++x)
        {
            for (int y = 0; y < affine_camera_matrix.cols(); ++y)
            {
                P_coefficients.push_back(
                    Eigen::Triplet<float>(3 * i + x, 4 * i + y, affine_camera_matrix(x, y)));
            }
        }
    }
    P.setFromTriplets(P_coefficients.begin(), P_coefficients.end());
 
    // The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
    VectorXf y = VectorXf::Ones(3 * num_landmarks);
    for (int i = 0; i < num_landmarks; ++i)
    {
        y(3 * i) = landmarks[i][0];
        y((3 * i) + 1) = landmarks[i][1];
        // y_((3 * i) + 2) = 1; // already 1, stays (homogeneous coordinate)
    }
 
    // The mean, with an added homogeneous coordinate (x_1, y_1, z_1, 1, x_2, ...)^t
    VectorXf v_bar = VectorXf::Ones(4 * num_landmarks);
    for (int i = 0; i < num_landmarks; ++i)
    {
        v_bar(4 * i) = face_instance(vertex_ids[i] * 3);
        v_bar((4 * i) + 1) = face_instance(vertex_ids[i] * 3 + 1);
        v_bar((4 * i) + 2) = face_instance(vertex_ids[i] * 3 + 2);
        // v_bar((4 * i) + 3) = 1; // already 1, stays (homogeneous coordinate)
    }
 
    // Bring into standard least squares form:
    const MatrixXf A = P * V_hat_h;   // camera matrix times the basis
    const MatrixXf b = P * v_bar - y; // camera matrix times the mean, minus the landmarks
 
    // Solve using NNLS:
    VectorXf coefficients;
    Eigen::NNLS<MatrixXf> nnls(A, 100); // In the multi-image fitting, the solver sometimes got stuck. This limits it to 100 iterations. It normaly converges within <10 iterations. This should be fine for the single-image fitting as well.
    bool non_singular = nnls.solve(-b);
    coefficients.noalias() = nnls.x();
 
    return std::vector<float>(coefficients.data(), coefficients.data() + coefficients.size());
};
 
} /* namespace fitting */
} /* namespace eos */
 
#endif /* BLENDSHAPEFITTING_HPP_ */