Merge branch 'devel' into multi_image_fit_devel

CMakeLists.txt

 cmake_minimum_required(VERSION 3.1.3)
 project(eos)
 set(eos_VERSION_MAJOR 0)
-set(eos_VERSION_MINOR 11)
-set(eos_VERSION_PATCH 0)
+set(eos_VERSION_MINOR 12)
+set(eos_VERSION_PATCH 1)
 set(eos_VERSION ${eos_VERSION_MAJOR}.${eos_VERSION_MINOR}.${eos_VERSION_PATCH})
 
 set_property(GLOBAL PROPERTY USE_FOLDERS ON)

MANIFEST.in

 include README.md LICENSE
 global-include CMakeLists.txt *.cmake
-recursive-include 3rdparty/cereal-1.1.1/include *
+recursive-include 3rdparty/cereal/include *
+recursive-include 3rdparty/eigen/Eigen *
 recursive-include 3rdparty/eigen3-nnls/src *.h
 recursive-include 3rdparty/glm/glm *
 recursive-include 3rdparty/nanoflann/include *
 recursive-include 3rdparty/pybind11 *
-recursive-include cmake *
 recursive-include include *
 recursive-include python *

README.md

@@ -32,9 +32,10 @@ At the moment, it mainly provides the following functionality:
 To use the library in your own project, just add the following directories to your include path:
 
 * `eos/include`
-* `eos/3rdparty/cereal-1.1.1/include`
+* `eos/3rdparty/cereal/include`
 * `eos/3rdparty/glm`
 * `eos/3rdparty/nanoflann/include`
+* `eos/3rdparty/eigen/Eigen`
 * `eos/3rdparty/eigen3-nnls/src`
 
 **Make sure to clone with `--recursive` to download the required submodules!**
@@ -84,9 +85,10 @@ The full model is available at [http://www.cvssp.org/facemodel](http://www.cvssp
 
 ## Python bindings
 
-eos includes python bindings for some of its functionality (and more can be added!). Set `-DEOS_GENERATE_PYTHON_BINDINGS=on` when running `cmake` to build them (and optionally set `PYTHON_EXECUTABLE` to point to your python interpreter if it's not found automatically).
+eos includes python bindings for some of its functionality (and more can be added!). An experimental package is on PyPI: Try `pip install eos-py`. You will still need the data files from this repository.
+In case of issues, build the bindings manually: Clone the repository and set `-DEOS_GENERATE_PYTHON_BINDINGS=on` when running `cmake` (and optionally set `PYTHON_EXECUTABLE` to point to your python interpreter if it's not found automatically).
 
-After building the bindings, they can be used like any python module:
+After having obtained the bindings, they can be used like any python module:
 
 ```
 import eos

examples/fit-model-ceres.cpp

@@ -380,7 +380,7 @@ int main(int argc, char *argv[])
 	// Colour model fitting (this needs a Morphable Model with colour (albedo) model, see note above main()):
 	if (!morphable_model.has_color_model())
 	{
-		cout << "The MorphableModel used does not contain a colour (albedo) model. ImageCost requires a model that contains a colour PCA model. You may want to use the full Surrey Face Model or remove this section.";
+		cout << "Error: The MorphableModel used does not contain a colour (albedo) model. ImageCost requires a model that contains a colour PCA model. You may want to use the full Surrey Face Model or remove this section.";
 		return EXIT_FAILURE;
 	}
 	std::vector<double> colour_coefficients;

include/eos/fitting/blendshape_fitting.hpp

@@ -24,10 +24,10 @@
 
 #include "eos/morphablemodel/Blendshape.hpp"
 
-#include "Eigen/Core" // for nnls.h
+#include "Eigen/Core"
+#include "Eigen/QR"
 #include "nnls.h"
 
-#include "Eigen/Core"
 #include "opencv2/core/core.hpp"
 
 #include <vector>
@@ -61,62 +61,56 @@ inline std::vector<float> fit_blendshapes_to_landmarks_linear(const std::vector<
 	assert(landmarks.size() == vertex_ids.size());
 
 	using cv::Mat;
+	using Eigen::VectorXf;
 	using Eigen::MatrixXf;
 
 	const int 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:
-	Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> blendshapes_as_basis = morphablemodel::to_matrix(blendshapes);
-	// Above converts to a RowMajor matrix on return - for now, since the core algorithm still uses cv::Mat (and OpenCV stores data in row-major memory order).
-	Mat blendshapes_basis_as_mat = Mat(blendshapes_as_basis.rows(), blendshapes_as_basis.cols(), CV_32FC1, blendshapes_as_basis.data());
+	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}$:
-	Mat V_hat_h = Mat::zeros(4 * num_landmarks, num_blendshapes, CV_32FC1);
+	MatrixXf V_hat_h = MatrixXf::Zero(4 * num_landmarks, num_blendshapes);
 	int row_index = 0;
 	for (int i = 0; i < num_landmarks; ++i) {
-		Mat basis_rows = blendshapes_basis_as_mat.rowRange(vertex_ids[i] * 3, (vertex_ids[i] * 3) + 3);
-		basis_rows.copyTo(V_hat_h.rowRange(row_index, row_index + 3));
+		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:
-	Mat P = Mat::zeros(3 * num_landmarks, 4 * num_landmarks, CV_32FC1);
+	MatrixXf P = MatrixXf::Zero(3 * num_landmarks, 4 * num_landmarks);
 	for (int i = 0; i < num_landmarks; ++i) {
-		Mat submatrix_to_replace = P.colRange(4 * i, (4 * i) + 4).rowRange(3 * i, (3 * i) + 3);
-		affine_camera_matrix.copyTo(submatrix_to_replace);
+		using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+		P.block(3 * i, 4 * i, 3, 4) = Eigen::Map<RowMajorMatrixXf>(affine_camera_matrix.ptr<float>(), affine_camera_matrix.rows, affine_camera_matrix.cols);
 	}
 
 	// The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
-	Mat y = Mat::ones(3 * num_landmarks, 1, CV_32FC1);
+	VectorXf y = VectorXf::Ones(3 * num_landmarks);
 	for (int i = 0; i < num_landmarks; ++i) {
-		y.at<float>(3 * i, 0) = landmarks[i][0];
-		y.at<float>((3 * i) + 1, 0) = landmarks[i][1];
-		//y.at<float>((3 * i) + 2, 0) = 1; // already 1, stays (homogeneous coordinate)
+		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
-	Mat v_bar = Mat::ones(4 * num_landmarks, 1, CV_32FC1);
+	VectorXf v_bar = VectorXf::Ones(4 * num_landmarks);
 	for (int i = 0; i < num_landmarks; ++i) {
-		//cv::Vec4f model_mean = morphable_model.get_shape_model().get_mean_at_point(vertex_ids[i]);
-		cv::Vec4f model_mean(face_instance(vertex_ids[i]*3), face_instance(vertex_ids[i]*3 + 1), face_instance(vertex_ids[i]*3 + 2), 1.0f);
-		v_bar.at<float>(4 * i, 0) = model_mean[0];
-		v_bar.at<float>((4 * i) + 1, 0) = model_mean[1];
-		v_bar.at<float>((4 * i) + 2, 0) = model_mean[2];
-		//v_bar.at<float>((4 * i) + 3, 0) = 1; // already 1, stays (homogeneous coordinate)
-		// note: now that a Vec4f is returned, we could use copyTo?
+		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:
-	Mat A = P * V_hat_h; // camera matrix times the basis
-	Mat b = P * v_bar - y; // camera matrix times the mean, minus the landmarks.
+	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
 	
-	Mat AtAReg = A.t() * A + lambda * Mat::eye(num_blendshapes, num_blendshapes, CV_32FC1);
-	// Solve using OpenCV:
-	Mat c_s;
-	bool non_singular = cv::solve(AtAReg, -A.t() * b, c_s, cv::DECOMP_SVD); // DECOMP_SVD calculates the pseudo-inverse if the matrix is not invertible.
-	// Because we're using SVD, non_singular will always be true. If we were to use e.g. Cholesky, we could return an expected<T>.
+	const MatrixXf AtAReg = A.transpose() * A + lambda * Eigen::MatrixXf::Identity(num_blendshapes, num_blendshapes);
+	const MatrixXf rhs = -A.transpose() * b;
 
-	return std::vector<float>(c_s);
+	const VectorXf coefficients = AtAReg.colPivHouseholderQr().solve(rhs);
+
+	return std::vector<float>(coefficients.data(), coefficients.data() + coefficients.size());
 };
 
 /**
@@ -139,65 +133,53 @@ inline std::vector<float> fit_blendshapes_to_landmarks_nnls(const std::vector<eo
 {
 	assert(landmarks.size() == vertex_ids.size());
 
-	using cv::Mat;
+	using Eigen::VectorXf;
 	using Eigen::MatrixXf;
 
 	const int 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:
-	Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> blendshapes_as_basis = morphablemodel::to_matrix(blendshapes);
-	// Above converts to a RowMajor matrix on return - for now, since the core algorithm still uses cv::Mat (and OpenCV stores data in row-major memory order).
-	Mat blendshapes_basis_as_mat = Mat(blendshapes_as_basis.rows(), blendshapes_as_basis.cols(), CV_32FC1, blendshapes_as_basis.data());
+	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}$:
-	Mat V_hat_h = Mat::zeros(4 * num_landmarks, num_blendshapes, CV_32FC1);
+	MatrixXf V_hat_h = MatrixXf::Zero(4 * num_landmarks, num_blendshapes);
 	int row_index = 0;
 	for (int i = 0; i < num_landmarks; ++i) {
-		Mat basis_rows = blendshapes_basis_as_mat.rowRange(vertex_ids[i] * 3, (vertex_ids[i] * 3) + 3);
-		basis_rows.copyTo(V_hat_h.rowRange(row_index, row_index + 3));
+		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:
-	Mat P = Mat::zeros(3 * num_landmarks, 4 * num_landmarks, CV_32FC1);
+	MatrixXf P = MatrixXf::Zero(3 * num_landmarks, 4 * num_landmarks);
 	for (int i = 0; i < num_landmarks; ++i) {
-		Mat submatrix_to_replace = P.colRange(4 * i, (4 * i) + 4).rowRange(3 * i, (3 * i) + 3);
-		affine_camera_matrix.copyTo(submatrix_to_replace);
+		using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+		P.block(3 * i, 4 * i, 3, 4) = Eigen::Map<RowMajorMatrixXf>(affine_camera_matrix.ptr<float>(), affine_camera_matrix.rows, affine_camera_matrix.cols);
 	}
-
 	// The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
-	Mat y = Mat::ones(3 * num_landmarks, 1, CV_32FC1);
+	VectorXf y = VectorXf::Ones(3 * num_landmarks);
 	for (int i = 0; i < num_landmarks; ++i) {
-		y.at<float>(3 * i, 0) = landmarks[i][0];
-		y.at<float>((3 * i) + 1, 0) = landmarks[i][1];
-		//y.at<float>((3 * i) + 2, 0) = 1; // already 1, stays (homogeneous coordinate)
+		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
-	Mat v_bar = Mat::ones(4 * num_landmarks, 1, CV_32FC1);
+	VectorXf v_bar = VectorXf::Ones(4 * num_landmarks);
 	for (int i = 0; i < num_landmarks; ++i) {
-		cv::Vec4f model_mean(face_instance(vertex_ids[i]*3), face_instance(vertex_ids[i]*3 + 1), face_instance(vertex_ids[i]*3 + 2), 1.0f);
-		v_bar.at<float>(4 * i, 0) = model_mean[0];
-		v_bar.at<float>((4 * i) + 1, 0) = model_mean[1];
-		v_bar.at<float>((4 * i) + 2, 0) = model_mean[2];
-		//v_bar.at<float>((4 * i) + 3, 0) = 1; // already 1, stays (homogeneous coordinate)
-		// note: now that a Vec4f is returned, we could use copyTo?
+		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:
-	Mat A = P * V_hat_h; // camera matrix times the basis
-	Mat b = P * v_bar - y; // camera matrix times the mean, minus the landmarks.
-
+	// 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:
-	using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
-	Eigen::Map<RowMajorMatrixXf> A_Eigen(A.ptr<float>(), A.rows, A.cols);
-	Eigen::Map<RowMajorMatrixXf> b_Eigen(b.ptr<float>(), b.rows, b.cols);
-
-	Eigen::VectorXf x;
-	bool non_singular = Eigen::NNLS<Eigen::MatrixXf>::solve(A_Eigen, -b_Eigen, x);
-	Mat c_s(x.rows(), x.cols(), CV_32FC1, x.data()); // create an OpenCV Mat header for the Eigen data
+	VectorXf coefficients;
+	bool non_singular = Eigen::NNLS<MatrixXf>::solve(A, -b, coefficients);
 
-	return std::vector<float>(c_s);
+	return std::vector<float>(coefficients.data(), coefficients.data() + coefficients.size());
 };
 
 	} /* namespace fitting */

include/eos/fitting/fitting.hpp

@@ -538,7 +538,7 @@ inline std::pair<core::Mesh, fitting::RenderingParameters> fit_shape_and_pose(co
  * @param[in] lambda Regularisation parameter of the PCA shape fitting.
  * @return The fitted model shape instance and the final pose.
  */
-inline std::pair<core::Mesh, fitting::RenderingParameters> fit_shape_and_pose(const morphablemodel::MorphableModel& morphable_model, const std::vector<morphablemodel::Blendshape>& blendshapes, const core::LandmarkCollection<cv::Vec2f>& landmarks, const core::LandmarkMapper& landmark_mapper, int image_width, int image_height, const morphablemodel::EdgeTopology& edge_topology, const fitting::ContourLandmarks& contour_landmarks, const fitting::ModelContour& model_contour, int num_iterations = 5, boost::optional<int> num_shape_coefficients_to_fit = boost::none, float lambda = 30.0f)
+inline std::pair<core::Mesh, fitting::RenderingParameters> fit_shape_and_pose(const morphablemodel::MorphableModel& morphable_model, const std::vector<morphablemodel::Blendshape>& blendshapes, const core::LandmarkCollection<cv::Vec2f>& landmarks, const core::LandmarkMapper& landmark_mapper, int image_width, int image_height, const morphablemodel::EdgeTopology& edge_topology, const fitting::ContourLandmarks& contour_landmarks, const fitting::ModelContour& model_contour, int num_iterations = 5, boost::optional<int> num_shape_coefficients_to_fit = boost::none, float lambda = 50.0f)
 {
 	std::vector<float> pca_coeffs;
 	std::vector<float> blendshape_coeffs;

include/eos/fitting/linear_shape_fitting.hpp

@@ -24,7 +24,7 @@
 
 #include "eos/morphablemodel/MorphableModel.hpp"
 
-//#include "Eigen/LU"
+#include "Eigen/QR"
 
 #include "opencv2/core/core.hpp"
 
@@ -57,10 +57,14 @@ namespace eos {
  * @param[in] model_standard_deviation The standard deviation of the 3D vertex points in the 3D model, projected to 2D (so the value is in pixels).
  * @return The estimated shape-coefficients (alphas).
  */
-inline std::vector<float> fit_shape_to_landmarks_linear_multi(morphablemodel::MorphableModel morphable_model, std::vector<cv::Mat> affine_camera_matrix, std::vector<std::vector<cv::Vec2f>>& landmarks, std::vector<std::vector<int>>& vertex_ids, std::vector<Eigen::VectorXf> base_face=std::vector<Eigen::VectorXf>(), float lambda=3.0f, boost::optional<int> num_coefficients_to_fit=boost::optional<int>(), boost::optional<float> detector_standard_deviation=boost::optional<float>(), boost::optional<float> model_standard_deviation=boost::optional<float>())
+inline std::vector<float> fit_shape_to_landmarks_linear_multi(const morphablemodel::MorphableModel& morphable_model, std::vector<cv::Mat> affine_camera_matrix, std::vector<std::vector<cv::Vec2f>>& landmarks, std::vector<std::vector<int>>& vertex_ids, std::vector<Eigen::VectorXf> base_face=std::vector<Eigen::VectorXf>(), float lambda=3.0f, boost::optional<int> num_coefficients_to_fit=boost::optional<int>(), boost::optional<float> detector_standard_deviation=boost::optional<float>(), boost::optional<float> model_standard_deviation=boost::optional<float>())
 {
-	using cv::Mat;
+
 	assert(affine_camera_matrix.size() == landmarks.size() && landmarks.size() == vertex_ids.size()); // same number of instances (i.e. images/frames) for each of them
+	assert(landmarks.size() == vertex_ids.size());
+
+	using Eigen::VectorXf;
+	using Eigen::MatrixXf;
 
 	int num_coeffs_to_fit = num_coefficients_to_fit.get_value_or(morphable_model.get_shape_model().get_num_principal_components());
 	int num_images = affine_camera_matrix.size();
@@ -72,20 +76,27 @@ inline std::vector<float> fit_shape_to_landmarks_linear_multi(morphablemodel::Mo
 	for (auto&& l : landmarks) {
 		total_num_landmarks_dimension += l.size();
 	}
+
 	// $\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}$:
-	Mat V_hat_h = Mat::zeros(4 * total_num_landmarks_dimension, num_coeffs_to_fit, CV_32FC1);
+    MatrixXf V_hat_h = MatrixXf::Zero(4 * total_num_landmarks_dimension, num_coeffs_to_fit);
 	int V_hat_h_row_index = 0;
 	// 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:
-	Mat P = Mat::zeros(3 * total_num_landmarks_dimension, 4 * total_num_landmarks_dimension, CV_32FC1);
+    MatrixXf P = MatrixXf::Zero(3 * total_num_landmarks_dimension, 4 * total_num_landmarks_dimension);
 	int P_index = 0;
-    Mat Omega = Mat::zeros(3 * total_num_landmarks_dimension, 3 * total_num_landmarks_dimension, CV_32FC1);
-    int Omega_index = 0; // this runs the same as P_index
+	// The variances: Add the 2D and 3D standard deviations.
+	// If the user doesn't provide them, we choose the following:
+	// 2D (detector) standard deviation: In pixel, we follow [1] and choose sqrt(3) as the default value.
+	// 3D (model) variance: 0.0f. It only makes sense to set it to something when we have a different variance for different vertices.
+	// The 3D variance has to be projected to 2D (for details, see paper [1]) so the units do match up.
+	float sigma_squared_2D = std::pow(detector_standard_deviation.get_value_or(std::sqrt(3.0f)), 2) + std::pow(model_standard_deviation.get_value_or(0.0f), 2);
+	// We use a VectorXf, and later use .asDiagonal():
+	VectorXf Omega = VectorXf::Constant(3 * total_num_landmarks_dimension, 1.0f / sigma_squared_2D);
 	// The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
-	Mat y = Mat::ones(3 * total_num_landmarks_dimension, 1, CV_32FC1);
+	VectorXf y = VectorXf::Ones(3 * total_num_landmarks_dimension);
 	int y_index = 0; // also runs the same as P_index. Should rename to "running_index"?
 	// The mean, with an added homogeneous coordinate (x_1, y_1, z_1, 1, x_2, ...)^t
-	Mat v_bar = Mat::ones(4 * total_num_landmarks_dimension, 1, CV_32FC1);
+	VectorXf v_bar = VectorXf::Ones(4 * total_num_landmarks_dimension);
 	int v_bar_index = 0; // also runs the same as P_index. But be careful, if I change it to be only 1 variable, only increment it once! :-)
 						 // Well I think that would make it a bit messy since we need to increment inside the for (landmarks...) loop. Try to refactor some other way.
 
@@ -105,82 +116,50 @@ inline std::vector<float> fit_shape_to_landmarks_linear_multi(morphablemodel::Mo
         // And we insert a row of zeros after every third row, resulting in matrix $\hat{V}_h \in R^{4N\times m-1}$:
         //Mat V_hat_h = Mat::zeros(4 * num_landmarks, num_coeffs_to_fit, CV_32FC1);
         for (int i = 0; i < num_landmarks; ++i) {
-            Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> basis_rows_ = morphable_model.get_shape_model().get_rescaled_pca_basis_at_point(vertex_ids[k][i]); // In the paper, the orthonormal basis might be used? I'm not sure, check it. It's even a mess in the paper. PH 26.5.2014: I think the rescaled basis is fine/better.
-            // Above converts to a RowMajor matrix on return - for now, since the core algorithm still uses cv::Mat (and OpenCV stores data in row-major memory order).
-            Mat basis_rows = Mat(basis_rows_.rows(), basis_rows_.cols(), CV_32FC1, basis_rows_.data());
-            //basisRows.copyTo(V_hat_h.rowRange(rowIndex, rowIndex + 3));
-            basis_rows.colRange(0, num_coeffs_to_fit).copyTo(V_hat_h.rowRange(V_hat_h_row_index, V_hat_h_row_index + 3));
+            MatrixXf basis_rows_ = morphable_model.get_shape_model().get_rescaled_pca_basis_at_point(vertex_ids[k][i]); // In the paper, the orthonormal basis might be used? I'm not sure, check it. It's even a mess in the paper. PH 26.5.2014: I think the rescaled basis is fine/better.
+            V_hat_h.block(V_hat_h_row_index, 0, 3, V_hat_h.cols()) = basis_rows_.block(0, 0, basis_rows_.rows(), num_coeffs_to_fit);
             V_hat_h_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:
         //Mat P = Mat::zeros(3 * num_landmarks, 4 * num_landmarks, CV_32FC1);
         for (int i = 0; i < num_landmarks; ++i) {
-            Mat submatrix_to_replace = P.colRange(4 * P_index, (4 * P_index) + 4).rowRange(3 * P_index, (3 * P_index) + 3);
-            affine_camera_matrix[k].copyTo(submatrix_to_replace);
+            using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+            P.block(3 * P_index, 4 * P_index, 3, 4) = Eigen::Map<RowMajorMatrixXf>(affine_camera_matrix[k].ptr<float>(), affine_camera_matrix[k].rows, affine_camera_matrix[k].cols);
             ++P_index;
         }
-        // The variances: Add the 2D and 3D standard deviations.
-        // If the user doesn't provide them, we choose the following:
-        // 2D (detector) standard deviation: In pixel, we follow [1] and choose sqrt(3) as the default value.
-        // 3D (model) variance: 0.0f. It only makes sense to set it to something when we have a different variance for different vertices.
-        // The 3D variance has to be projected to 2D (for details, see paper [1]) so the units do match up.
-        float sigma_squared_2D = std::pow(detector_standard_deviation.get_value_or(std::sqrt(3.0f)), 2) + std::pow(model_standard_deviation.get_value_or(0.0f), 2);
-        //Mat Sigma = Mat::zeros(3 * num_landmarks, 3 * num_landmarks, CV_32FC1);
-        for (int i = 0; i < 3 * num_landmarks; ++i) {
-            // Sigma(i, i) = sqrt(sigma_squared_2D), but then Omega is Sigma.t() * Sigma (squares the diagonal) - so we just assign 1/sigma_squared_2D to Omega here:
-            Omega.at<float>(Omega_index, Omega_index) = 1.0f / sigma_squared_2D; // the higher the sigma_squared_2D, the smaller the diagonal entries of Sigma will be
-            ++Omega_index;
-        }
 
         // The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
         //Mat y = Mat::ones(3 * num_landmarks, 1, CV_32FC1);
         for (int i = 0; i < num_landmarks; ++i) {
-            y.at<float>(3 * y_index, 0) = landmarks[k][i][0];
-            y.at<float>((3 * y_index) + 1, 0) = landmarks[k][i][1];
-            //y.at<float>((3 * i) + 2, 0) = 1; // already 1, stays (homogeneous coordinate)
+            y(3 * y_index) = landmarks[k][i][0];
+            y((3 * y_index) + 1) = landmarks[k][i][1];
+            //y((3 * i) + 2) = 1; // already 1, stays (homogeneous coordinate)
             ++y_index;
         }
         // The mean, with an added homogeneous coordinate (x_1, y_1, z_1, 1, x_2, ...)^t
         //Mat v_bar = Mat::ones(4 * num_landmarks, 1, CV_32FC1);
         for (int i = 0; i < num_landmarks; ++i) {            
-            //cv::Vec4f model_mean = morphable_model.get_shape_model().get_mean_at_point(vertex_ids[i]);
-            cv::Vec4f model_mean(base_face[k](vertex_ids[k][i] * 3), base_face[k](vertex_ids[k][i] * 3 + 1), base_face[k](vertex_ids[k][i] * 3 + 2), 1.0f);
-            v_bar.at<float>(4 * v_bar_index, 0) = model_mean[0];
-            v_bar.at<float>((4 * v_bar_index) + 1, 0) = model_mean[1];
-            v_bar.at<float>((4 * v_bar_index) + 2, 0) = model_mean[2];
-            //v_bar.at<float>((4 * i) + 3, 0) = 1; // already 1, stays (homogeneous coordinate)
+            v_bar(4 * v_bar_index) = base_face[k](vertex_ids[k][i] * 3);
+            v_bar((4 * v_bar_index) + 1) = base_face[k](vertex_ids[k][i] * 3 + 1);
+            v_bar((4 * v_bar_index) + 2) = base_face[k](vertex_ids[k][i] * 3 + 2);
+            //v_bar.at<float>((4 * i) + 3) = 1; // already 1, stays (homogeneous coordinate)
             ++v_bar_index;
-            // note: now that a Vec4f is returned, we could use copyTo?
         }
     }
 
-	// Bring into standard regularised quadratic form with diagonal distance matrix Omega
-	Mat A = P * V_hat_h; // camera matrix times the basis
-	Mat b = P * v_bar - y; // camera matrix times the mean, minus the landmarks.
-	//Mat c_s; // The x, we solve for this! (the variance-normalised shape parameter vector, $c_s = [a_1/sigma_{s,1} , ..., a_m-1/sigma_{s,m-1}]^t$
-	//int numShapePc = morphableModel.getShapeModel().getNumberOfPrincipalComponents();
-	const int num_shape_pc = num_coeffs_to_fit;
-	Mat AtOmegaA = A.t() * Omega * A;
-	Mat AtOmegaAReg = AtOmegaA + lambda * Mat::eye(num_shape_pc, num_shape_pc, CV_32FC1);
-
-	// Invert (and perform some sanity checks) using Eigen:
-/*	using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
-	Eigen::Map<RowMajorMatrixXf> AtOmegaAReg_Eigen(AtOmegaAReg.ptr<float>(), AtOmegaAReg.rows, AtOmegaAReg.cols);
-	Eigen::FullPivLU<RowMajorMatrixXf> luOfAtOmegaAReg(AtOmegaAReg_Eigen); // Calculate the full-pivoting LU decomposition of the regularized AtA. Note: We could also try FullPivHouseholderQR if our system is non-minimal (i.e. there are more constraints than unknowns).
-	auto rankOfAtOmegaAReg = luOfAtOmegaAReg.rank();
-	bool isAtOmegaARegInvertible = luOfAtOmegaAReg.isInvertible();
-	float threshold = std::abs(luOfAtOmegaAReg.maxPivot()) * luOfAtOmegaAReg.threshold(); // originaly "2 * ..." but I commented it out
-	RowMajorMatrixXf AtARegInv_EigenFullLU = luOfAtOmegaAReg.inverse(); // Note: We should use ::solve() instead
-	Mat AtOmegaARegInvFullLU(AtARegInv_EigenFullLU.rows(), AtARegInv_EigenFullLU.cols(), CV_32FC1, AtARegInv_EigenFullLU.data()); // create an OpenCV Mat header for the Eigen data
-*/
-	// Solve using OpenCV:
-	Mat c_s; // Note/Todo: We get coefficients ~ N(0, sigma) I think. They are not multiplied with the eigenvalues.
-	bool non_singular = cv::solve(AtOmegaAReg, -A.t() * Omega.t() * b, c_s, cv::DECOMP_SVD); // DECOMP_SVD calculates the pseudo-inverse if the matrix is not invertible.
-	// Because we're using SVD, non_singular will always be true. If we were to use e.g. Cholesky, we could return an expected<T>.
-
-	return std::vector<float>(c_s);
+	// Bring into standard regularised quadratic form with diagonal distance matrix Omega:
+	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 AtOmegaAReg = A.transpose() * Omega.asDiagonal() * A + lambda * Eigen::MatrixXf::Identity(num_coeffs_to_fit, num_coeffs_to_fit);
+	const MatrixXf rhs = -A.transpose() * Omega.asDiagonal() * b; // It's -A^t*Omega^t*b, but we don't need to transpose Omega, since it's a diagonal matrix, and Omega^t = Omega.
+	// c_s: The 'x' that we solve for. (The variance-normalised shape parameter vector, $c_s = [a_1/sigma_{s,1} , ..., a_m-1/sigma_{s,m-1}]^t$.)
+	// We get coefficients ~ N(0, 1), because we're fitting with the rescaled basis. The coefficients are not multiplied with their eigenvalues.
+	const VectorXf c_s = AtOmegaAReg.colPivHouseholderQr().solve(rhs);
+
+	return std::vector<float>(c_s.data(), c_s.data() + c_s.size());
 };
 
+
 /**
  * Fits the shape of a Morphable Model to given 2D landmarks (i.e. estimates the maximum likelihood solution of the shape coefficients) as proposed in [1].
  * It's a linear, closed-form solution fitting of the shape, with regularisation (prior towards the mean).
@@ -209,7 +188,6 @@ inline std::vector<float> fit_shape_to_landmarks_linear(const morphablemodel::Mo
     return fit_shape_to_landmarks_linear_multi(morphable_model, { affine_camera_matrix }, all_landmarks, all_vertex_ids, { base_face }, lambda, num_coefficients_to_fit, detector_standard_deviation, model_standard_deviation );
 }
 
-
 	} /* namespace fitting */
 } /* namespace eos */
 

include/eos/morphablemodel/Blendshape.hpp

@@ -102,6 +102,20 @@ inline Eigen::MatrixXf to_matrix(const std::vector<Blendshape>& blendshapes)
 	return blendshapes_as_basis;
 };
 
+/**
+ * @brief Maps an std::vector of coefficients with Eigen::Map, so it can be multiplied
+ * with a blendshapes matrix.
+ *
+ * Note that  the resulting Eigen::Map only lives as long as the data given lives and is in scope.
+ *
+ * @param[in] coefficients Vector of blendshape coefficients.
+ * @return An Eigen::Map pointing to the given coefficients data.
+ */
+inline Eigen::Map<const Eigen::VectorXf> to_vector(const std::vector<float>& coefficients)
+{
+	return Eigen::Map<const Eigen::VectorXf>(coefficients.data(), coefficients.size());
+};
+
 	} /* namespace morphablemodel */
 } /* namespace eos */
 

matlab/+eos/+fitting/fit_shape_and_pose.m

@@ -33,7 +33,7 @@ end
 % We'll use default values to the following arguments, if they're not
 % provided:
 if (~exist('edge_topology', 'var')), edge_topology = '../share/sfm_3448_edge_topology.json'; end
-if (~exist('contour_landmarks', 'var')), contour_landmarks = '../share/ibug2did.txt'; end
+if (~exist('contour_landmarks', 'var')), contour_landmarks = '../share/ibug_to_sfm.txt'; end
 if (~exist('model_contour', 'var')), model_contour = '../share/model_contours.json'; end
 if (~exist('num_iterations', 'var')), num_iterations = 5; end
 if (~exist('num_shape_coefficients_to_fit', 'var')), num_shape_coefficients_to_fit = -1; end

matlab/demo.m

@@ -3,7 +3,7 @@
 %% Set up some required paths to files:
 model_file = '../share/sfm_shape_3448.bin';
 blendshapes_file = '../share/expression_blendshapes_3448.bin';
-landmark_mappings = '../share/ibug2did.txt';
+landmark_mappings = '../share/ibug_to_sfm.txt';
 
 %% Load an image and its landmarks in ibug format:
 image = imread('../bin/data/image_0010.png');

python/demo.py

@@ -10,9 +10,9 @@ def main():
 
     model = eos.morphablemodel.load_model("../share/sfm_shape_3448.bin")
     blendshapes = eos.morphablemodel.load_blendshapes("../share/expression_blendshapes_3448.bin")
-    landmark_mapper = eos.core.LandmarkMapper('../share/ibug2did.txt')
+    landmark_mapper = eos.core.LandmarkMapper('../share/ibug_to_sfm.txt')
     edge_topology = eos.morphablemodel.load_edge_topology('../share/sfm_3448_edge_topology.json')
-    contour_landmarks = eos.fitting.ContourLandmarks.load('../share/ibug2did.txt')
+    contour_landmarks = eos.fitting.ContourLandmarks.load('../share/ibug_to_sfm.txt')
     model_contour = eos.fitting.ModelContour.load('../share/model_contours.json')
 
     (mesh, pose, shape_coeffs, blendshape_coeffs) = eos.fitting.fit_shape_and_pose(model, blendshapes,

python/generate-python-bindings.cpp

@@ -105,6 +105,8 @@ PYBIND11_PLUGIN(eos) {
 		.def("get_shape_model", [](const morphablemodel::MorphableModel& m) { return m.get_shape_model(); }, "Returns the PCA shape model of this Morphable Model.") // Not sure if that'll really be const in Python? I think Python does a copy each time this gets called?
 		.def("get_color_model", [](const morphablemodel::MorphableModel& m) { return m.get_color_model(); }, "Returns the PCA colour (albedo) model of this Morphable Model.") // (continued from above:) We may want to use py::overload, but in any case, we need to tell pybind11 if it should use the const or non-const overload.
 		.def("get_mean", &morphablemodel::MorphableModel::get_mean, "Returns the mean of the shape- and colour model as a Mesh.")
+		.def("draw_sample", (core::Mesh(morphablemodel::MorphableModel::*)(std::vector<float>, std::vector<float>) const)&morphablemodel::MorphableModel::draw_sample, "Returns a sample from the model with the given shape- and colour PCA coefficients.", py::arg("shape_coefficients"), py::arg("color_coefficients"))
+		.def("has_color_model", &morphablemodel::MorphableModel::has_color_model, "Returns true if this Morphable Model contains a colour model, and false if it is a shape-only model.")
 		;
 
 	morphablemodel_module.def("load_model", &morphablemodel::load_model, "Load a Morphable Model from a cereal::BinaryInputArchive (.bin) from the harddisk.", py::arg("filename"));

setup.py

@@ -88,7 +88,7 @@ class CMakeBuild(build_ext):
 
 setup(
     name='eos-py',
-    version='0.11.0',
+    version='0.12.1',
     author='Patrik Huber',
     author_email='patrikhuber@gmail.com',
     description='Python bindings for eos - A lightweight 3D Morphable Face Model fitting library in modern C++11/14',

share/readme.txt

@@ -10,7 +10,7 @@ Files in this directory:
 - sfm_shape_3448.bin:
 	The public shape-only Surrey 3D Morphable Face Model.
 	To obtain a full 3DMM and higher resolution levels, follow the instructions
-	at <todo: add link to the page of the Uni>.
+	at cvssp.org/facemodel.
 	Details about the different models can be found in:
 	"A Multiresolution 3D Morphable Face Model and Fitting Framework",
 	P. Huber, G. Hu, R. Tena, P. Mortazavian, W. Koppen, W. Christmas, M. Rätsch, J. Kittler,
@@ -22,7 +22,9 @@ Files in this directory:
 
 - sfm_3448_edge_topology.json:
 	Contains a precomputed list of the model's edges, and the two faces and vertices that are
-	adjacent to each edge. Used in the edge-fitting.
+	adjacent to each edge. Uses 1-based indexing ("0" has a special meaning of "no adjacent
+	vertex/edge") - this may change to 0-based in the future to be consistent with the rest of
+	the library. The file is used in the edge-fitting.
 
 - model_contours.json:
 	Definition of the model's contour vertices of the right and left side of the face.

utils/CMakeLists.txt

@@ -28,13 +28,11 @@ target_link_libraries(scm-to-cereal eos ${OpenCV_LIBS} ${Boost_LIBRARIES})
 add_executable(bfm-binary-to-cereal bfm-binary-to-cereal.cpp)
 target_link_libraries(bfm-binary-to-cereal eos ${OpenCV_LIBS} ${Boost_LIBRARIES})
 
-
 # Reads an edgestruct CSV file created from Matlab, and converts it to json:
 add_executable(edgestruct-csv-to-json edgestruct-csv-to-json.cpp)
 target_link_libraries(edgestruct-csv-to-json eos ${Boost_LIBRARIES})
 
-
-# install target:
+# Install targets:
 install(TARGETS scm-to-cereal DESTINATION bin)
 install(TARGETS bfm-binary-to-cereal DESTINATION bin)
 install(TARGETS edgestruct-csv-to-json DESTINATION bin)