eos  0.7.1
linear_shape_fitting.hpp
1 /*
2  * Eos - A 3D Morphable Model fitting library written in modern C++11/14.
3  *
4  * File: include/eos/fitting/linear_shape_fitting.hpp
5  *
6  * Copyright 2014, 2015 Patrik Huber
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 #pragma once
21 
22 #ifndef LINEARSHAPEFITTING_HPP_
23 #define LINEARSHAPEFITTING_HPP_
24 
25 #include "eos/morphablemodel/MorphableModel.hpp"
26 
27 //#include "Eigen/LU"
28 
29 #include "opencv2/core/core.hpp"
30 
31 #include "boost/optional.hpp"
32 
33 #include <vector>
34 #include <cassert>
35 
36 namespace eos {
37  namespace fitting {
38 
60 inline std::vector<float> fit_shape_to_landmarks_linear(morphablemodel::MorphableModel morphable_model, cv::Mat affine_camera_matrix, std::vector<cv::Vec2f> landmarks, std::vector<int> vertex_ids, cv::Mat base_face=cv::Mat(), 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>())
61 {
62  using cv::Mat;
63  assert(landmarks.size() == vertex_ids.size());
64 
65  int num_coeffs_to_fit = num_coefficients_to_fit.get_value_or(morphable_model.get_shape_model().get_num_principal_components());
66  int num_landmarks = static_cast<int>(landmarks.size());
67 
68  if (base_face.empty())
69  {
70  base_face = morphable_model.get_shape_model().get_mean();
71  }
72 
73  // $\hat{V} \in R^{3N\times m-1}$, subselect the rows of the eigenvector matrix $V$ associated with the $N$ feature points
74  // And we insert a row of zeros after every third row, resulting in matrix $\hat{V}_h \in R^{4N\times m-1}$:
75  Mat V_hat_h = Mat::zeros(4 * num_landmarks, num_coeffs_to_fit, CV_32FC1);
76  int row_index = 0;
77  for (int i = 0; i < num_landmarks; ++i) {
78  Mat basis_rows = morphable_model.get_shape_model().get_normalised_pca_basis(vertex_ids[i]); // In the paper, the not-normalised 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 normalised basis is fine/better.
79  //basisRows.copyTo(V_hat_h.rowRange(rowIndex, rowIndex + 3));
80  basis_rows.colRange(0, num_coeffs_to_fit).copyTo(V_hat_h.rowRange(row_index, row_index + 3));
81  row_index += 4; // replace 3 rows and skip the 4th one, it has all zeros
82  }
83  // 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:
84  Mat P = Mat::zeros(3 * num_landmarks, 4 * num_landmarks, CV_32FC1);
85  for (int i = 0; i < num_landmarks; ++i) {
86  Mat submatrix_to_replace = P.colRange(4 * i, (4 * i) + 4).rowRange(3 * i, (3 * i) + 3);
87  affine_camera_matrix.copyTo(submatrix_to_replace);
88  }
89  // The variances: Add the 2D and 3D standard deviations.
90  // If the user doesn't provide them, we choose the following:
91  // 2D (detector) standard deviation: In pixel, we follow [1] and choose sqrt(3) as the default value.
92  // 3D (model) variance: 0.0f. It only makes sense to set it to something when we have a different variance for different vertices.
93  // The 3D variance has to be projected to 2D (for details, see paper [1]) so the units do match up.
94  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);
95  Mat Sigma = Mat::zeros(3 * num_landmarks, 3 * num_landmarks, CV_32FC1);
96  for (int i = 0; i < 3 * num_landmarks; ++i) {
97  Sigma.at<float>(i, i) = 1.0f / std::sqrt(sigma_squared_2D); // the higher the sigma_squared_2D, the smaller the diagonal entries of Sigma will be
98  }
99  Mat Omega = Sigma.t() * Sigma; // just squares the diagonal
100  // The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
101  Mat y = Mat::ones(3 * num_landmarks, 1, CV_32FC1);
102  for (int i = 0; i < num_landmarks; ++i) {
103  y.at<float>(3 * i, 0) = landmarks[i][0];
104  y.at<float>((3 * i) + 1, 0) = landmarks[i][1];
105  //y.at<float>((3 * i) + 2, 0) = 1; // already 1, stays (homogeneous coordinate)
106  }
107  // The mean, with an added homogeneous coordinate (x_1, y_1, z_1, 1, x_2, ...)^t
108  Mat v_bar = Mat::ones(4 * num_landmarks, 1, CV_32FC1);
109  for (int i = 0; i < num_landmarks; ++i) {
110  //cv::Vec4f model_mean = morphable_model.get_shape_model().get_mean_at_point(vertex_ids[i]);
111  cv::Vec4f model_mean(base_face.at<float>(vertex_ids[i] * 3), base_face.at<float>(vertex_ids[i] * 3 + 1), base_face.at<float>(vertex_ids[i] * 3 + 2), 1.0f);
112  v_bar.at<float>(4 * i, 0) = model_mean[0];
113  v_bar.at<float>((4 * i) + 1, 0) = model_mean[1];
114  v_bar.at<float>((4 * i) + 2, 0) = model_mean[2];
115  //v_bar.at<float>((4 * i) + 3, 0) = 1; // already 1, stays (homogeneous coordinate)
116  // note: now that a Vec4f is returned, we could use copyTo?
117  }
118 
119  // Bring into standard regularised quadratic form with diagonal distance matrix Omega
120  Mat A = P * V_hat_h; // camera matrix times the basis
121  Mat b = P * v_bar - y; // camera matrix times the mean, minus the landmarks.
122  //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$
123  //int numShapePc = morphableModel.getShapeModel().getNumberOfPrincipalComponents();
124  const int num_shape_pc = num_coeffs_to_fit;
125  Mat AtOmegaA = A.t() * Omega * A;
126  Mat AtOmegaAReg = AtOmegaA + lambda * Mat::eye(num_shape_pc, num_shape_pc, CV_32FC1);
127 
128  // Invert (and perform some sanity checks) using Eigen:
129 /* using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
130  Eigen::Map<RowMajorMatrixXf> AtOmegaAReg_Eigen(AtOmegaAReg.ptr<float>(), AtOmegaAReg.rows, AtOmegaAReg.cols);
131  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).
132  auto rankOfAtOmegaAReg = luOfAtOmegaAReg.rank();
133  bool isAtOmegaARegInvertible = luOfAtOmegaAReg.isInvertible();
134  float threshold = std::abs(luOfAtOmegaAReg.maxPivot()) * luOfAtOmegaAReg.threshold(); // originaly "2 * ..." but I commented it out
135  RowMajorMatrixXf AtARegInv_EigenFullLU = luOfAtOmegaAReg.inverse(); // Note: We should use ::solve() instead
136  Mat AtOmegaARegInvFullLU(AtARegInv_EigenFullLU.rows(), AtARegInv_EigenFullLU.cols(), CV_32FC1, AtARegInv_EigenFullLU.data()); // create an OpenCV Mat header for the Eigen data
137 */
138  // Solve using OpenCV:
139  Mat c_s; // Note/Todo: We get coefficients ~ N(0, sigma) I think. They are not multiplied with the eigenvalues.
140  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.
141  // 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>.
142 
143  return std::vector<float>(c_s);
144 };
145 
146  } /* namespace fitting */
147 } /* namespace eos */
148 
149 #endif /* LINEARSHAPEFITTING_HPP_ */
cv::Mat get_normalised_pca_basis() const
Definition: PcaModel.hpp:188
int get_num_principal_components() const
Definition: PcaModel.hpp:84
Namespace containing all of eos&#39;s 3D model fitting functionality.
std::vector< float > fit_shape_to_landmarks_linear(morphablemodel::MorphableModel morphable_model, cv::Mat affine_camera_matrix, std::vector< cv::Vec2f > landmarks, std::vector< int > vertex_ids, cv::Mat base_face=cv::Mat(), 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 >())
Definition: linear_shape_fitting.hpp:60
cv::Mat get_mean() const
Definition: PcaModel.hpp:119
PcaModel get_shape_model() const
Definition: MorphableModel.hpp:77
A class representing a 3D Morphable Model, consisting of a shape- and colour (albedo) PCA model...
Definition: MorphableModel.hpp:54