Commit bbc82d0d authored by Patrik Huber's avatar Patrik Huber

Moved shape and camera fitting to headers, renamed to snake_case

parent 77f271f1
......@@ -48,8 +48,6 @@ set(SOURCES
src/eos/morphablemodel/PcaModel.cpp
src/eos/morphablemodel/MorphableModel.cpp
src/eos/morphablemodel/io/cvssp.cpp
src/eos/fitting/AffineCameraEstimation.cpp
src/eos/fitting/LinearShapeFitting.cpp
src/eos/render/Mesh.cpp
src/eos/render/utils.cpp
)
......
......@@ -174,7 +174,7 @@ int main(int argc, char *argv[])
}
// Estimate the camera from the 2D - 3D point correspondences
Mat affineCam = fitting::estimateAffineCamera(imagePoints, modelPoints);
Mat affineCam = fitting::estimate_affine_camera(imagePoints, modelPoints);
// Draw the mean-face landmarks projected using the estimated camera:
for (auto&& vertex : modelPoints) {
......@@ -184,7 +184,7 @@ int main(int argc, char *argv[])
// Estimate the shape coefficients by fitting the shape to the landmarks:
float lambda = 5.0f; ///< the regularisation parameter
vector<float> fittedCoeffs = fitting::fitShapeToLandmarksLinear(morphableModel, affineCam, imagePoints, vertexIndices, lambda);
vector<float> fittedCoeffs = fitting::fit_shape_to_landmarks_linear(morphableModel, affineCam, imagePoints, vertexIndices, lambda);
// Obtain the full mesh and draw it using the estimated camera:
render::Mesh mesh = morphableModel.drawSample(fittedCoeffs, vector<float>());
......
/*
* Eos - A 3D Morphable Model fitting library written in modern C++11/14.
*
* File: include/eos/fitting/AffineCameraEstimation.hpp
*
* Copyright 2014, 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 AFFINECAMERAESTIMATION_HPP_
#define AFFINECAMERAESTIMATION_HPP_
#include "eos/render/utils.hpp"
#include "opencv2/core/core.hpp"
#include <vector>
namespace eos {
namespace fitting {
/**
* The Gold Standard Algorithm for estimating an affine
* camera matrix from world to image correspondences.
* See Algorithm 7.2 in Multiple View Geometry, Hartley &
* Zisserman, 2nd Edition, 2003.
*
* Requires >= 4 corresponding points.
*
* @param[in] imagePoints A list of 2D image points.
* @param[in] modelPoints Corresponding points of a 3D model.
* @return A 3x4 affine camera matrix (the third row is [0, 0, 0, 1]).
*/
cv::Mat estimateAffineCamera(std::vector<cv::Vec2f> imagePoints, std::vector<cv::Vec4f> modelPoints);
/**
* Projects a point from world coordinates to screen coordinates.
* First, an estimated affine camera matrix is used to transform
* the point to clip space. Second, the point is transformed to
* screen coordinates using the window transform. The window transform
* also flips the y-axis (the image origin is top-left, while in
* clip space top is +1 and bottom is -1).
*
* Note: Assumes the affine camera matrix only projects from world
* to clip space, because a subsequent window transform is applied.
*
* @param[in] vertex A vertex in 3D space. vertex[3] = 1.0f.
* @param[in] affine_camera_matrix A 3x4 affine camera matrix.
* @param[in] screen_width Width of the screen or window used for projection.
* @param[in] screen_height Height of the screen or window used for projection.
* @return A vector with x and y coordinates transformed to screen coordinates.
*/
inline cv::Vec2f project_affine(cv::Vec4f vertex, cv::Mat affine_camera_matrix, int screen_width, int screen_height)
{
// Transform to clip space:
cv::Mat clip_coords = affine_camera_matrix * cv::Mat(vertex);
// Take the x and y coordinates in clip space and apply the window transform:
cv::Vec2f screen_coords = render::clip_to_screen_space(cv::Vec2f(clip_coords.rowRange(0, 2)), screen_width, screen_height);
return screen_coords;
};
} /* namespace fitting */
} /* namespace eos */
#endif /* AFFINECAMERAESTIMATION_HPP_ */
/*
* Eos - A 3D Morphable Model fitting library written in modern C++11/14.
*
* File: include/eos/fitting/LinearShapeFitting.hpp
*
* Copyright 2014, 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 LINEARSHAPEFITTING_HPP_
#define LINEARSHAPEFITTING_HPP_
#include "eos/morphablemodel/MorphableModel.hpp"
#include "opencv2/core/core.hpp"
#include "boost/optional.hpp"
#include <vector>
namespace eos {
namespace fitting {
/**
* 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).
*
* [1] O. Aldrian & W. Smith, Inverse Rendering of Faces with a 3D Morphable Model, PAMI 2013.
*
* Note: Using less than the maximum number of coefficients to fit is not thoroughly tested yet and may contain an error.
* Note: Returns coefficients following standard normal distribution (i.e. all have similar magnitude). Why? Because we fit using the normalised basis?
* Note: The standard deviations given should be a vector, i.e. different for each landmark. This is not implemented yet.
*
* @param[in] morphableModel The Morphable Model whose shape (coefficients) are estimated.
* @param[in] affineCameraMatrix A 3x4 affine camera matrix from world to clip-space (should probably be of type CV_32FC1 as all our calculations are done with float).
* @param[in] landmarks 2D landmarks from an image, given in clip-coordinates.
* @param[in] vertexIds The vertex ids in the model that correspond to the 2D points.
* @param[in] lambda The regularisation parameter (weight of the prior towards the mean).
* @param[in] numCoefficientsToFit How many shape-coefficients to fit (all others will stay 0). Not tested thoroughly.
* @param[in] detectorStandardDeviation The standard deviation of the 2D landmarks given (e.g. of the detector used).
* @param[in] modelStandardDeviation The standard deviation of the 3D vertex points in the 3D model.
* @return The estimated shape-coefficients (alphas).
*/
std::vector<float> fitShapeToLandmarksLinear(morphablemodel::MorphableModel morphableModel, cv::Mat affineCameraMatrix, std::vector<cv::Vec2f> landmarks, std::vector<int> vertexIds, float lambda=20.0f, boost::optional<int> numCoefficientsToFit=boost::optional<int>(), boost::optional<float> detectorStandardDeviation=boost::optional<float>(), boost::optional<float> modelStandardDeviation=boost::optional<float>());
} /* namespace fitting */
} /* namespace eos */
#endif /* LINEARSHAPEFITTING_HPP_ */
/*
* Eos - A 3D Morphable Model fitting library written in modern C++11/14.
*
* File: src/eos/fitting/AffineCameraEstimation.cpp
* File: include/eos/fitting/AffineCameraEstimation.hpp
*
* Copyright 2014, 2015 Patrik Huber
*
......@@ -17,49 +17,55 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include "eos/fitting/AffineCameraEstimation.hpp"
#pragma once
#ifndef AFFINECAMERAESTIMATION_HPP_
#define AFFINECAMERAESTIMATION_HPP_
#include "eos/render/utils.hpp"
#include "opencv2/core/core.hpp"
#include "opencv2/core/core_c.h" // for CV_REDUCE_AVG
#include <iostream>
#include <exception>
#include <vector>
#include <cassert>
using cv::Mat;
using cv::Vec2f;
using cv::Vec3f;
using cv::Vec4f;
using std::vector;
namespace eos {
namespace fitting {
Mat estimateAffineCamera(vector<Vec2f> imagePoints, vector<Vec4f> modelPoints)
/**
* The Gold Standard Algorithm for estimating an affine
* camera matrix from world to image correspondences.
* See Algorithm 7.2 in Multiple View Geometry, Hartley &
* Zisserman, 2nd Edition, 2003.
*
* Requires >= 4 corresponding points.
*
* @param[in] image_points A list of 2D image points.
* @param[in] model_points Corresponding points of a 3D model.
* @return A 3x4 affine camera matrix (the third row is [0, 0, 0, 1]).
*/
cv::Mat estimate_affine_camera(std::vector<cv::Vec2f> image_points, std::vector<cv::Vec4f> model_points)
{
assert(imagePoints.size() == modelPoints.size());
using cv::Mat;
assert(image_points.size() == model_points.size());
assert(image_points.size() >= 4); // Number of correspondence points given needs to be equal to or larger than 4
const auto numCorrespondences = imagePoints.size();
if (numCorrespondences < 4) {
std::string errorMsg("AffineCameraEstimation: Number of points given needs to be equal to or larger than 4.");
std::cout << errorMsg << std::endl;
throw std::runtime_error(errorMsg);
}
const auto numCorrespondences = image_points.size();
Mat matImagePoints; // will be numCorrespondences x 2, CV_32FC1
Mat matModelPoints; // will be numCorrespondences x 3, CV_32FC1
for (int i = 0; i < imagePoints.size(); ++i) {
for (int i = 0; i < image_points.size(); ++i) {
Mat imgPoint(1, 2, CV_32FC1);
imgPoint.at<float>(0, 0) = imagePoints[i][0];
imgPoint.at<float>(0, 1) = imagePoints[i][1];
imgPoint.at<float>(0, 0) = image_points[i][0];
imgPoint.at<float>(0, 1) = image_points[i][1];
matImagePoints.push_back(imgPoint);
Mat mdlPoint(1, 3, CV_32FC1);
mdlPoint.at<float>(0, 0) = modelPoints[i][0];
mdlPoint.at<float>(0, 1) = modelPoints[i][1];
mdlPoint.at<float>(0, 2) = modelPoints[i][2];
mdlPoint.at<float>(0, 0) = model_points[i][0];
mdlPoint.at<float>(0, 1) = model_points[i][1];
mdlPoint.at<float>(0, 2) = model_points[i][2];
matModelPoints.push_back(mdlPoint);
}
......@@ -76,20 +82,20 @@ Mat estimateAffineCamera(vector<Vec2f> imagePoints, vector<Vec4f> modelPoints)
}
averageNorm /= matImagePoints.rows;
// 2) multiply every vectors coordinate by sqrt(2)/avgnorm
float scaleFactor = std::sqrt(2)/averageNorm;
float scaleFactor = std::sqrt(2) / averageNorm;
matImagePoints *= scaleFactor; // add unit homogeneous component here
// The points in matImagePoints now have a RMS distance from the origin of sqrt(2).
// The normalisation matrix so that the 2D points are mean-free and their norm is as described above.
// The points in matImagePoints now have a RMS distance from the origin of sqrt(2).
// The normalisation matrix so that the 2D points are mean-free and their norm is as described above.
Mat T = Mat::zeros(3, 3, CV_32FC1);
T.at<float>(0, 0) = scaleFactor; // s_x
T.at<float>(1, 1) = scaleFactor; // s_y
T.at<float>(0, 2) = -imagePointsMean.at<float>(0, 0) * scaleFactor; // t_x
T.at<float>(1, 2) = -imagePointsMean.at<float>(0, 1) * scaleFactor; // t_y
T.at<float>(2, 2) = 1;
// center the model points to the origin:
Mat tmpOrigMdlPoints = matModelPoints.clone(); // Temp for testing: Save the original coordinates.
// translate the centroid of the model points to the origin:
// translate the centroid of the model points to the origin:
Mat modelPointsMean; // use non-homogeneous coords for the next few steps? (less submatrices etc overhead)
cv::reduce(matModelPoints, modelPointsMean, 0, CV_REDUCE_AVG);
modelPointsMean = cv::repeat(modelPointsMean, matModelPoints.rows, 1);
......@@ -104,8 +110,8 @@ Mat estimateAffineCamera(vector<Vec2f> imagePoints, vector<Vec4f> modelPoints)
// 2) multiply every vectors coordinate by sqrt(3)/avgnorm
scaleFactor = std::sqrt(3) / averageNorm;
matModelPoints *= scaleFactor; // add unit homogeneous component here
// The points in matModelPoints now have a RMS distance from the origin of sqrt(3).
// The normalisation matrix so that the 3D points are mean-free and their norm is as described above.
// The points in matModelPoints now have a RMS distance from the origin of sqrt(3).
// The normalisation matrix so that the 3D points are mean-free and their norm is as described above.
Mat U = Mat::zeros(4, 4, CV_32FC1);
U.at<float>(0, 0) = scaleFactor; // s_x
U.at<float>(1, 1) = scaleFactor; // s_y
......@@ -123,16 +129,16 @@ Mat estimateAffineCamera(vector<Vec2f> imagePoints, vector<Vec4f> modelPoints)
//Mat p_8(); // p_8 is 8 x 1. We are solving for it.
Mat b(numCorrespondences * 2, 1, CV_32FC1);
for (int i = 0; i < numCorrespondences; ++i) {
A_8.at<float>(2*i, 0) = matModelPoints.at<float>(i, 0); // could maybe made faster by assigning the whole row/col-range if possible?
A_8.at<float>(2*i, 1) = matModelPoints.at<float>(i, 1);
A_8.at<float>(2*i, 2) = matModelPoints.at<float>(i, 2);
A_8.at<float>(2*i, 3) = 1;
A_8.at<float>((2*i)+1, 4) = matModelPoints.at<float>(i, 0);
A_8.at<float>((2*i)+1, 5) = matModelPoints.at<float>(i, 1);
A_8.at<float>((2*i)+1, 6) = matModelPoints.at<float>(i, 2);
A_8.at<float>((2*i)+1, 7) = 1;
b.at<float>(2*i, 0) = matImagePoints.at<float>(i, 0);
b.at<float>((2*i)+1, 0) = matImagePoints.at<float>(i, 1);
A_8.at<float>(2 * i, 0) = matModelPoints.at<float>(i, 0); // could maybe made faster by assigning the whole row/col-range if possible?
A_8.at<float>(2 * i, 1) = matModelPoints.at<float>(i, 1);
A_8.at<float>(2 * i, 2) = matModelPoints.at<float>(i, 2);
A_8.at<float>(2 * i, 3) = 1;
A_8.at<float>((2 * i) + 1, 4) = matModelPoints.at<float>(i, 0);
A_8.at<float>((2 * i) + 1, 5) = matModelPoints.at<float>(i, 1);
A_8.at<float>((2 * i) + 1, 6) = matModelPoints.at<float>(i, 2);
A_8.at<float>((2 * i) + 1, 7) = 1;
b.at<float>(2 * i, 0) = matImagePoints.at<float>(i, 0);
b.at<float>((2 * i) + 1, 0) = matImagePoints.at<float>(i, 1);
}
Mat p_8 = A_8.inv(cv::DECOMP_SVD) * b;
Mat C_tilde = Mat::zeros(3, 4, CV_32FC1);
......@@ -142,7 +148,35 @@ Mat estimateAffineCamera(vector<Vec2f> imagePoints, vector<Vec4f> modelPoints)
Mat P_Affine = T.inv() * C_tilde * U;
return P_Affine;
}
};
/**
* Projects a point from world coordinates to screen coordinates.
* First, an estimated affine camera matrix is used to transform
* the point to clip space. Second, the point is transformed to
* screen coordinates using the window transform. The window transform
* also flips the y-axis (the image origin is top-left, while in
* clip space top is +1 and bottom is -1).
*
* Note: Assumes the affine camera matrix only projects from world
* to clip space, because a subsequent window transform is applied.
*
* @param[in] vertex A vertex in 3D space. vertex[3] = 1.0f.
* @param[in] affine_camera_matrix A 3x4 affine camera matrix.
* @param[in] screen_width Width of the screen or window used for projection.
* @param[in] screen_height Height of the screen or window used for projection.
* @return A vector with x and y coordinates transformed to screen coordinates.
*/
inline cv::Vec2f project_affine(cv::Vec4f vertex, cv::Mat affine_camera_matrix, int screen_width, int screen_height)
{
// Transform to clip space:
cv::Mat clip_coords = affine_camera_matrix * cv::Mat(vertex);
// Take the x and y coordinates in clip space and apply the window transform:
cv::Vec2f screen_coords = render::clip_to_screen_space(cv::Vec2f(clip_coords.rowRange(0, 2)), screen_width, screen_height);
return screen_coords;
};
} /* namespace fitting */
} /* namespace eos */
#endif /* AFFINECAMERAESTIMATION_HPP_ */
/*
* Eos - A 3D Morphable Model fitting library written in modern C++11/14.
*
* File: src/eos/fitting/LinearShapeFitting.cpp
* File: include/eos/fitting/LinearShapeFitting.hpp
*
* Copyright 2014, 2015 Patrik Huber
*
......@@ -17,54 +17,81 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include "eos/fitting/LinearShapeFitting.hpp"
#pragma once
#include "eos/render/utils.hpp"
#ifndef LINEARSHAPEFITTING_HPP_
#define LINEARSHAPEFITTING_HPP_
#include <cassert>
#include "eos/morphablemodel/MorphableModel.hpp"
#include "Eigen/LU"
using eos::morphablemodel::MorphableModel;
using cv::Mat;
using std::vector;
#include "opencv2/core/core.hpp"
#include "boost/optional.hpp"
#include <vector>
#include <cassert>
namespace eos {
namespace fitting {
vector<float> fitShapeToLandmarksLinear(MorphableModel morphableModel, Mat affineCameraMatrix, vector<cv::Vec2f> landmarks, std::vector<int> vertexIds, float lambda/*=20.0f*/, boost::optional<int> numCoefficientsToFit/*=boost::optional<int>()*/, boost::optional<float> detectorStandardDeviation/*=boost::optional<float>()*/, boost::optional<float> modelStandardDeviation/*=boost::optional<float>()*/)
/**
* 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).
*
* [1] O. Aldrian & W. Smith, Inverse Rendering of Faces with a 3D Morphable Model, PAMI 2013.
*
* Note: Using less than the maximum number of coefficients to fit is not thoroughly tested yet and may contain an error.
* Note: Returns coefficients following standard normal distribution (i.e. all have similar magnitude). Why? Because we fit using the normalised basis?
* Note: The standard deviations given should be a vector, i.e. different for each landmark. This is not implemented yet.
*
* @param[in] morphable_model The Morphable Model whose shape (coefficients) are estimated.
* @param[in] affine_camera_matrix A 3x4 affine camera matrix from world to clip-space (should probably be of type CV_32FC1 as all our calculations are done with float).
* @param[in] landmarks 2D landmarks from an image, given in clip-coordinates.
* @param[in] vertex_ids The vertex ids in the model that correspond to the 2D points.
* @param[in] lambda The regularisation parameter (weight of the prior towards the mean).
* @param[in] num_coefficients_to_fit How many shape-coefficients to fit (all others will stay 0). Not tested thoroughly.
* @param[in] detector_standard_deviation The standard deviation of the 2D landmarks given (e.g. of the detector used).
* @param[in] model_standard_deviation The standard deviation of the 3D vertex points in the 3D model.
* @return The estimated shape-coefficients (alphas).
*/
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, float lambda=20.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>())
{
assert(landmarks.size() == vertexIds.size());
int numCoeffsToFit = numCoefficientsToFit.get_value_or(morphableModel.getShapeModel().getNumberOfPrincipalComponents());
using cv::Mat;
assert(landmarks.size() == vertex_ids.size());
int num_coeffs_to_fit = num_coefficients_to_fit.get_value_or(morphable_model.getShapeModel().getNumberOfPrincipalComponents());
// $\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 * landmarks.size(), numCoeffsToFit, CV_32FC1);
int rowIndex = 0;
Mat V_hat_h = Mat::zeros(4 * landmarks.size(), num_coeffs_to_fit, CV_32FC1);
int row_index = 0;
for (int i = 0; i < landmarks.size(); ++i) {
Mat basisRows = morphableModel.getShapeModel().getNormalisedPcaBasis(vertexIds[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.
//basisRows.copyTo(V_hat_h.rowRange(rowIndex, rowIndex + 3));
basisRows.colRange(0, numCoeffsToFit).copyTo(V_hat_h.rowRange(rowIndex, rowIndex + 3));
rowIndex += 4; // replace 3 rows and skip the 4th one, it has all zeros
Mat basis_rows = morphable_model.getShapeModel().getNormalisedPcaBasis(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.
//basisRows.copyTo(V_hat_h.rowRange(rowIndex, rowIndex + 3));
basis_rows.colRange(0, num_coeffs_to_fit).copyTo(V_hat_h.rowRange(row_index, row_index + 3));
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, affineCam) is placed on the diagonal:
Mat P = Mat::zeros(3 * landmarks.size(), 4 * landmarks.size(), CV_32FC1);
for (int i = 0; i < landmarks.size(); ++i) {
Mat submatrixToReplace = P.colRange(4 * i, (4 * i) + 4).rowRange(3 * i, (3 * i) + 3);
affineCameraMatrix.copyTo(submatrixToReplace);
Mat submatrix_to_replace = P.colRange(4 * i, (4 * i) + 4).rowRange(3 * i, (3 * i) + 3);
affine_camera_matrix.copyTo(submatrix_to_replace);
}
// The variances: Add the 2D and 3D standard deviations.
// If the user doesn't provide them, we choose the following:
// 2D (detector) variance: Assuming the detector has a standard deviation of 3 pixels, and the face size (IED) is around 80px. That's 3.75% of the IED. Assume that an image is on average 512x512px so 80/512 = 0.16 is the size the IED occupies inside an image.
// Now we're in clip-coords ([-1, 1]) and take 0.16 of the range [-1, 1], 0.16/2 = 0.08, and then the standard deviation of the detector is 3.75% of 0.08, i.e. 0.0375*0.08 = 0.003.
// 3D (model) variance: 0.0f. It only makes sense to set it to something when we have a different variance for different vertices.
float sigma_2D_3D = detectorStandardDeviation.get_value_or(0.003f) + modelStandardDeviation.get_value_or(0.0f);
float sigma_2D_3D = detector_standard_deviation.get_value_or(0.003f) + model_standard_deviation.get_value_or(0.0f);
// Note: Isn't it a bit strange to add these as they have different units/normalisations? Check the paper.
Mat Sigma = Mat::zeros(3 * landmarks.size(), 3 * landmarks.size(), CV_32FC1);
for (int i = 0; i < 3 * landmarks.size(); ++i) {
Sigma.at<float>(i, i) = 1.0f / sigma_2D_3D; // the higher the sigma_2D_3D, the smaller the diagonal entries of Sigma will be
}
Mat Omega = Sigma.t() * Sigma; // just squares the diagonal
// The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
// The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
Mat y = Mat::ones(3 * landmarks.size(), 1, CV_32FC1);
for (int i = 0; i < landmarks.size(); ++i) {
y.at<float>(3 * i, 0) = landmarks[i][0];
......@@ -74,10 +101,10 @@ vector<float> fitShapeToLandmarksLinear(MorphableModel morphableModel, Mat affin
// The mean, with an added homogeneous coordinate (x_1, y_1, z_1, 1, x_2, ...)^t
Mat v_bar = Mat::ones(4 * landmarks.size(), 1, CV_32FC1);
for (int i = 0; i < landmarks.size(); ++i) {
cv::Vec4f modelMean = morphableModel.getShapeModel().getMeanAtPoint(vertexIds[i]);
v_bar.at<float>(4 * i, 0) = modelMean[0];
v_bar.at<float>((4 * i) + 1, 0) = modelMean[1];
v_bar.at<float>((4 * i) + 2, 0) = modelMean[2];
cv::Vec4f model_mean = morphable_model.getShapeModel().getMeanAtPoint(vertex_ids[i]);
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?
}
......@@ -85,27 +112,30 @@ vector<float> fitShapeToLandmarksLinear(MorphableModel morphableModel, Mat affin
// 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();
int numShapePc = numCoeffsToFit;
//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();
int num_shape_pc = num_coeffs_to_fit;
Mat AtOmegaA = A.t() * Omega * A;
Mat AtOmegaAReg = AtOmegaA + lambda * Mat::eye(numShapePc, numShapePc, CV_32FC1);
Mat AtOmegaAReg = AtOmegaA + lambda * Mat::eye(num_shape_pc, num_shape_pc, CV_32FC1);
// Invert using OpenCV:
Mat AtOmegaARegInv = AtOmegaAReg.inv(cv::DECOMP_SVD); // DECOMP_SVD calculates the pseudo-inverse if the matrix is not invertible.
// Invert (and perform some sanity checks) using Eigen:
/*
Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> AtOmegaAReg_Eigen(AtOmegaAReg.ptr<float>(), AtOmegaAReg.rows, AtOmegaAReg.cols);
Eigen::FullPivLU<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> 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).
// 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).
float rankOfAtOmegaAReg = luOfAtOmegaAReg.rank();
bool isAtOmegaARegInvertible = luOfAtOmegaAReg.isInvertible();
Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> AtARegInv_EigenFullLU = luOfAtOmegaAReg.inverse();
float threshold = /*2 * */ std::abs(luOfAtOmegaAReg.maxPivot()) * luOfAtOmegaAReg.threshold();
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
*/
Mat AtOmegatb = A.t() * Omega.t() * b;
Mat c_s = -AtOmegaARegInv * AtOmegatb; // Note/Todo: We get coefficients ~ N(0, sigma) I think. They are not multiplied with the eigenvalues.
return vector<float>(c_s);
}
return std::vector<float>(c_s);
};
} /* namespace fitting */
} /* namespace eos */
#endif /* LINEARSHAPEFITTING_HPP_ */
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment