Merge branch 'multi_image_fit_devel' into multi-frame-richard

examples/CMakeLists.txt

@@ -30,6 +30,11 @@ add_executable(fit-model fit-model.cpp)
 target_link_libraries(fit-model eos ${OpenCV_LIBS} ${Boost_LIBRARIES})
 target_link_libraries(fit-model "$<$<CXX_COMPILER_ID:GNU>:-pthread>$<$<CXX_COMPILER_ID:Clang>:-pthreads>")
 
+# Model fitting example that fits orthographic camera, shape, blendshapes, and contours to multiple images:
+add_executable(fit-model-multi fit-model-multi.cpp)
+target_link_libraries(fit-model-multi eos ${OpenCV_LIBS} ${Boost_LIBRARIES})
+target_link_libraries(fit-model-multi "$<$<CXX_COMPILER_ID:GNU>:-pthread>$<$<CXX_COMPILER_ID:Clang>:-pthreads>")
+
 # Generate random samples from the model:
 add_executable(generate-obj generate-obj.cpp)
 target_link_libraries(generate-obj eos ${OpenCV_LIBS} ${Boost_LIBRARIES})
@@ -42,6 +47,7 @@ target_link_libraries(fit-model-multi-frame "$<$<CXX_COMPILER_ID:GNU>:-pthread>$
 # Install these targets:
 install(TARGETS fit-model-simple DESTINATION bin)
 install(TARGETS fit-model DESTINATION bin)
+install(TARGETS fit-model-multi DESTINATION bin)
 install(TARGETS generate-obj DESTINATION bin)
 install(TARGETS fit-model-multi-frame DESTINATION bin)
 install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/data DESTINATION bin)

examples/fit-model-multi.cpp

+/*
+ * eos - A 3D Morphable Model fitting library written in modern C++11/14.
+ *
+ * File: examples/fit-model.cpp
+ *
+ * Copyright 2016 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.
+ */
+#include "eos/core/Landmark.hpp"
+#include "eos/core/LandmarkMapper.hpp"
+#include "eos/morphablemodel/MorphableModel.hpp"
+#include "eos/morphablemodel/Blendshape.hpp"
+#include "eos/fitting/fitting.hpp"
+#include "eos/render/utils.hpp"
+#include "eos/render/texture_extraction.hpp"
+#include "eos/render/render.hpp"
+
+#include "opencv2/core/core.hpp"
+#include "opencv2/highgui/highgui.hpp"
+
+#include "boost/program_options.hpp"
+#include "boost/filesystem.hpp"
+
+#include <vector>
+#include <iostream>
+#include <fstream>
+
+using namespace eos;
+namespace po = boost::program_options;
+namespace fs = boost::filesystem;
+using eos::core::Landmark;
+using eos::core::LandmarkCollection;
+using cv::Mat;
+using cv::Vec2f;
+using cv::Vec3f;
+using cv::Vec4f;
+using std::cout;
+using std::endl;
+using std::vector;
+using std::string;
+
+/**
+ * Reads an ibug .pts landmark file and returns an ordered vector with
+ * the 68 2D landmark coordinates.
+ *
+ * @param[in] filename Path to a .pts file.
+ * @return An ordered vector with the 68 ibug landmarks.
+ */
+LandmarkCollection<cv::Vec2f> read_pts_landmarks(std::string filename)
+{
+	using std::getline;
+	using cv::Vec2f;
+	using std::string;
+	LandmarkCollection<Vec2f> landmarks;
+	landmarks.reserve(68);
+
+	std::ifstream file(filename);
+	if (!file.is_open()) {
+		throw std::runtime_error(string("Could not open landmark file: " + filename));
+	}
+
+	string line;
+	// Skip the first 3 lines, they're header lines:
+	getline(file, line); // 'version: 1'
+	getline(file, line); // 'n_points : 68'
+	getline(file, line); // '{'
+
+	int ibugId = 1;
+	while (getline(file, line))
+	{
+		if (line == "}") { // end of the file
+			break;
+		}
+		std::stringstream lineStream(line);
+
+		Landmark<Vec2f> landmark;
+		landmark.name = std::to_string(ibugId);
+		if (!(lineStream >> landmark.coordinates[0] >> landmark.coordinates[1])) {
+			throw std::runtime_error(string("Landmark format error while parsing the line: " + line));
+		}
+		// From the iBug website:
+		// "Please note that the re-annotated data for this challenge are saved in the Matlab convention of 1 being
+		// the first index, i.e. the coordinates of the top left pixel in an image are x=1, y=1."
+		// ==> So we shift every point by 1:
+		landmark.coordinates[0] -= 1.0f;
+		landmark.coordinates[1] -= 1.0f;
+		landmarks.emplace_back(landmark);
+		++ibugId;
+	}
+	return landmarks;
+};
+
+/**
+ * Draws the given mesh as wireframe into the image.
+ *
+ * It does backface culling, i.e. draws only vertices in CCW order.
+ *
+ * @param[in] image An image to draw into.
+ * @param[in] mesh The mesh to draw.
+ * @param[in] modelview Model-view matrix to draw the mesh.
+ * @param[in] projection Projection matrix to draw the mesh.
+ * @param[in] viewport Viewport to draw the mesh.
+ * @param[in] colour Colour of the mesh to be drawn.
+ */
+void draw_wireframe(cv::Mat image, const core::Mesh& mesh, glm::mat4x4 modelview, glm::mat4x4 projection, glm::vec4 viewport, cv::Scalar colour = cv::Scalar(0, 255, 0, 255))
+{
+	for (const auto& triangle : mesh.tvi)
+	{
+		const auto p1 = glm::project({ mesh.vertices[triangle[0]][0], mesh.vertices[triangle[0]][1], mesh.vertices[triangle[0]][2] }, modelview, projection, viewport);
+		const auto p2 = glm::project({ mesh.vertices[triangle[1]][0], mesh.vertices[triangle[1]][1], mesh.vertices[triangle[1]][2] }, modelview, projection, viewport);
+		const auto p3 = glm::project({ mesh.vertices[triangle[2]][0], mesh.vertices[triangle[2]][1], mesh.vertices[triangle[2]][2] }, modelview, projection, viewport);
+		if (render::detail::are_vertices_ccw_in_screen_space(glm::vec2(p1), glm::vec2(p2), glm::vec2(p3)))
+		{
+			cv::line(image, cv::Point(p1.x, p1.y), cv::Point(p2.x, p2.y), colour);
+			cv::line(image, cv::Point(p2.x, p2.y), cv::Point(p3.x, p3.y), colour);
+			cv::line(image, cv::Point(p3.x, p3.y), cv::Point(p1.x, p1.y), colour);
+		}
+	}
+};
+
+/**
+ * @brief Merges isomaps from a live video with a weighted averaging, based
+ * on the view angle of each vertex to the camera.
+ *
+ * An optional merge_threshold can be specified upon construction. Pixels with
+ * a view-angle above that threshold will be completely discarded. All pixels
+ * below the threshold are merged with a weighting based on its vertex view-angle.
+ * Assumes the isomaps to be 512x512.
+ */
+class WeightedIsomapAveraging
+{
+public:
+    /**
+     * @brief Constructs a new object that will hold the current averaged isomap and
+     * be able to add frames from a live video and merge them on-the-fly.
+     *
+     * The threshold means: Each triangle with a view angle smaller than the given angle will be used to merge.
+     * The default threshold (90°) means all triangles, as long as they're a little bit visible, are merged.
+     *
+     * @param[in] merge_threshold View-angle merge threshold, in degrees, from 0 to 90.
+     */
+    WeightedIsomapAveraging(float merge_threshold = 90.0f)
+    {
+        assert(merge_threshold >= 0.f && merge_threshold <= 90.f);
+
+        visibility_counter = cv::Mat::zeros(512, 512, CV_32SC1);
+        merged_isomap = cv::Mat::zeros(512, 512, CV_32FC4);
+
+        // map 0° to 255, 90° to 0:
+        float alpha_thresh = (-255.f / 90.f) * merge_threshold + 255.f;
+        if (alpha_thresh < 0.f) // could maybe happen due to float inaccuracies / rounding?
+            alpha_thresh = 0.0f;
+        threshold = static_cast<unsigned char>(alpha_thresh);
+    };
+
+    /**
+     * @brief Merges the given new isomap with all previously processed isomaps.
+     *
+     * @param[in] isomap The new isomap to add.
+     * @return The merged isomap of all images processed so far, as 8UC4.
+     */
+    cv::Mat add_and_merge(const cv::Mat& isomap)
+    {
+        // Merge isomaps, add the current to the already merged, pixel by pixel:
+        for (int r = 0; r < isomap.rows; ++r)
+        {
+            for (int c = 0; c < isomap.cols; ++c)
+            {
+                if (isomap.at<cv::Vec4b>(r, c)[3] <= threshold)
+                {
+                    continue; // ignore this pixel, not visible in the extracted isomap of this current frame
+                }
+                // we're sure to have a visible pixel, merge it:
+                // merged_pixel = (old_average * visible_count + new_pixel) / (visible_count + 1)
+                merged_isomap.at<cv::Vec4f>(r, c)[0] = (merged_isomap.at<cv::Vec4f>(r, c)[0] * visibility_counter.at<int>(r, c) + isomap.at<cv::Vec4b>(r, c)[0]) / (visibility_counter.at<int>(r, c) + 1);
+                merged_isomap.at<cv::Vec4f>(r, c)[1] = (merged_isomap.at<cv::Vec4f>(r, c)[1] * visibility_counter.at<int>(r, c) + isomap.at<cv::Vec4b>(r, c)[1]) / (visibility_counter.at<int>(r, c) + 1);
+                merged_isomap.at<cv::Vec4f>(r, c)[2] = (merged_isomap.at<cv::Vec4f>(r, c)[2] * visibility_counter.at<int>(r, c) + isomap.at<cv::Vec4b>(r, c)[2]) / (visibility_counter.at<int>(r, c) + 1);
+                merged_isomap.at<cv::Vec4f>(r, c)[3] = 255; // as soon as we've seen the pixel visible once, we set it to visible.
+                ++visibility_counter.at<int>(r, c);
+            }
+        }
+        cv::Mat merged_isomap_uchar;
+        merged_isomap.convertTo(merged_isomap_uchar, CV_8UC4);
+        return merged_isomap_uchar;
+    };
+
+private:
+    cv::Mat visibility_counter;
+    cv::Mat merged_isomap;
+    unsigned char threshold;
+};
+
+/**
+ * This app demonstrates estimation of the camera and fitting of the shape
+ * model of a 3D Morphable Model from an ibug LFPW image with its landmarks.
+ * In addition to fit-model-simple, this example uses blendshapes, contour-
+ * fitting, and can iterate the fitting.
+ *
+ * 68 ibug landmarks are loaded from the .pts file and converted
+ * to vertex indices using the LandmarkMapper.
+ */
+int main(int argc, char *argv[])
+{
+    fs::path modelfile, isomapfile, mappingsfile, contourfile, edgetopologyfile, blendshapesfile, outputfilebase;
+    vector<fs::path> imagefiles, landmarksfiles;
+	try {
+		po::options_description desc("Allowed options");
+		desc.add_options()
+			("help,h",
+				"display the help message")
+			("model,m", po::value<fs::path>(&modelfile)->required()->default_value("../share/sfm_shape_3448.bin"),
+				"a Morphable Model stored as cereal BinaryArchive")
+            //("image,i", po::value<vector<fs::path>>(&imagefiles)->required()->default_value("data/image_0010.png"),
+            ("image,i", po::value<vector<fs::path>>(&imagefiles)->multitoken(),
+				"an input image")
+            ("landmarks,l", po::value<vector<fs::path>>(&landmarksfiles)->multitoken(),
+				"2D landmarks for the image, in ibug .pts format")
+            ("mapping,p", po::value<fs::path>(&mappingsfile)->required()->default_value("../share/ibug_to_sfm.txt"),
+				"landmark identifier to model vertex number mapping")
+			("model-contour,c", po::value<fs::path>(&contourfile)->required()->default_value("../share/model_contours.json"),
+				"file with model contour indices")
+			("edge-topology,e", po::value<fs::path>(&edgetopologyfile)->required()->default_value("../share/sfm_3448_edge_topology.json"),
+				"file with model's precomputed edge topology")
+			("blendshapes,b", po::value<fs::path>(&blendshapesfile)->required()->default_value("../share/expression_blendshapes_3448.bin"),
+				"file with blendshapes")
+            ("output,o", po::value<fs::path>(&outputfilebase)->required()->default_value("out"),
+				"basename for the output rendering and obj files")
+			;
+		po::variables_map vm;
+		po::store(po::command_line_parser(argc, argv).options(desc).run(), vm);
+		if (vm.count("help")) {
+			cout << "Usage: fit-model [options]" << endl;
+			cout << desc;
+			return EXIT_SUCCESS;
+		}
+		po::notify(vm);
+	}
+	catch (const po::error& e) {
+		cout << "Error while parsing command-line arguments: " << e.what() << endl;
+		cout << "Use --help to display a list of options." << endl;
+		return EXIT_SUCCESS;
+	}
+
+    if (landmarksfiles.size() != imagefiles.size()) {
+        cout << "Number of landmarksfiles not equal to number of images given: "<<landmarksfiles.size() <<"!=" <<imagefiles.size()<< endl;
+        return EXIT_SUCCESS;
+    }
+
+    if (landmarksfiles.empty()) {
+        cout << "Please give at least 1 image and landmarkfile" << endl;
+        return EXIT_SUCCESS;
+    }
+	// Load the image, landmarks, LandmarkMapper and the Morphable Model:
+    vector<Mat> images;
+    for (auto& imagefile : imagefiles){
+        images.push_back(cv::imread(imagefile.string()));
+    }
+    vector<LandmarkCollection<cv::Vec2f>> landmarkss;
+	try {
+        for (auto& landmarksfile : landmarksfiles){
+            landmarkss.push_back(read_pts_landmarks(landmarksfile.string()));
+        }
+	}
+	catch (const std::runtime_error& e) {
+		cout << "Error reading the landmarks: " << e.what() << endl;
+		return EXIT_FAILURE;
+	}
+	morphablemodel::MorphableModel morphable_model;
+	try {
+		morphable_model = morphablemodel::load_model(modelfile.string());
+	}
+	catch (const std::runtime_error& e) {
+		cout << "Error loading the Morphable Model: " << e.what() << endl;
+		return EXIT_FAILURE;
+	}
+	// The landmark mapper is used to map ibug landmark identifiers to vertex ids:
+	core::LandmarkMapper landmark_mapper = mappingsfile.empty() ? core::LandmarkMapper() : core::LandmarkMapper(mappingsfile);
+
+	// The expression blendshapes:
+	vector<morphablemodel::Blendshape> blendshapes = morphablemodel::load_blendshapes(blendshapesfile.string());
+
+	// These two are used to fit the front-facing contour to the ibug contour landmarks:
+	fitting::ModelContour model_contour = contourfile.empty() ? fitting::ModelContour() : fitting::ModelContour::load(contourfile.string());
+	fitting::ContourLandmarks ibug_contour = fitting::ContourLandmarks::load(mappingsfile.string());
+
+	// The edge topology is used to speed up computation of the occluding face contour fitting:
+	morphablemodel::EdgeTopology edge_topology = morphablemodel::load_edge_topology(edgetopologyfile.string());
+
+	// Fit the model, get back a mesh and the pose:
+    vector<core::Mesh> meshs;
+    vector<fitting::RenderingParameters> rendering_paramss;
+    vector<int> image_widths;
+    vector<int> image_heights;
+    for (auto& image : images) {
+        image_widths.push_back(image.cols);
+        image_heights.push_back(image.rows);
+    }
+
+    std::vector<float> pca_shape_coefficients;
+    std::vector<std::vector<float>> blendshape_coefficients;
+    std::vector<std::vector<cv::Vec2f>> fitted_image_points;
+
+    std::tie(meshs, rendering_paramss) = fitting::fit_shape_and_pose_multi(morphable_model, blendshapes, landmarkss, landmark_mapper, image_widths, image_heights, edge_topology, ibug_contour, model_contour, 50, boost::none, 30.0f, boost::none, pca_shape_coefficients, blendshape_coefficients, fitted_image_points);
+    //std::tie(meshs, rendering_paramss) = fitting::fit_shape_and_pose_multi(morphable_model, blendshapes, landmarks, landmark_mapper, image_width, image_height, edge_topology, contour_landmarks, model_contour, num_iterations, num_shape_coefficients_to_fit, lambda, boost::none, pca_shape_coefficients, blendshape_coefficients, fitted_image_points);
+    //fit_shape_and_pose_multi(const morphablemodel::MorphableModel& morphable_model, const std::vector<morphablemodel::Blendshape>& blendshapes, const std::vector<core::LandmarkCollection<cv::Vec2f>>& landmarks, const core::LandmarkMapper& landmark_mapper, std::vector<int> image_width, std::vector<int> image_height, const morphablemodel::EdgeTopology& edge_topology, const fitting::ContourLandmarks& contour_landmarks, const fitting::ModelContour& model_contour, int num_iterations, boost::optional<int> num_shape_coefficients_to_fit, float lambda, boost::optional<fitting::RenderingParameters> initial_rendering_params, std::vector<float>& pca_shape_coefficients, std::vector<std::vector<float>>& blendshape_coefficients, std::vector<std::vector<cv::Vec2f>>& fitted_image_points)
+
+    std::cout<<"final pca shape coefficients: ";
+    for (auto i: pca_shape_coefficients)
+      std::cout << i << ' ';
+    std::cout << std::endl;
+
+    WeightedIsomapAveraging isomap_averaging(60.f); // merge all triangles that are facing <60° towards the camera
+    Mat merged_isomap;
+
+    for (unsigned i =0; i <images.size(); ++i) {
+
+        // The 3D head pose can be recovered as follows:
+        float yaw_angle = glm::degrees(glm::yaw(rendering_paramss[i].get_rotation()));
+        // and similarly for pitch and roll.
+
+        // Extract the texture from the image using given mesh and camera parameters:
+        Mat affine_from_ortho = fitting::get_3x4_affine_camera_matrix(rendering_paramss[i], images[i].cols, images[i].rows);
+        Mat isomap = render::extract_texture(meshs[i], affine_from_ortho, images[i]);
+
+        // Draw the loaded landmarks:
+        Mat outimg = images[i].clone();
+        for (auto&& lm : landmarkss[i]) {
+            cv::rectangle(outimg, cv::Point2f(lm.coordinates[0] - 2.0f, lm.coordinates[1] - 2.0f), cv::Point2f(lm.coordinates[0] + 2.0f, lm.coordinates[1] + 2.0f), { 255, 0, 0 });
+        }
+
+        // Draw the fitted mesh as wireframe, and save the image:
+        draw_wireframe(outimg, meshs[i], rendering_paramss[i].get_modelview(), rendering_paramss[i].get_projection(), fitting::get_opencv_viewport(images[i].cols, images[i].rows));
+        fs::path outputfile = outputfilebase;
+        outputfile += fs::path(imagefiles[i].stem());
+        outputfile += fs::path(".png");
+        cv::imwrite(outputfile.string(), outimg);
+
+        //save frontal rendering with texture:
+        Mat frontal_rendering;
+        glm::mat4 modelview_frontal = glm::mat4( 1.0 );
+        core::Mesh neutral_expression = morphablemodel::sample_to_mesh(morphable_model.get_shape_model().draw_sample(pca_shape_coefficients), morphable_model.get_color_model().get_mean(), morphable_model.get_shape_model().get_triangle_list(), morphable_model.get_color_model().get_triangle_list(), morphable_model.get_texture_coordinates());
+        std::tie(frontal_rendering, std::ignore) = render::render(neutral_expression, modelview_frontal, glm::ortho(-130.0f, 130.0f, -130.0f, 130.0f), 256, 256, render::create_mipmapped_texture(isomap), true, false, false);
+        outputfile.replace_extension(".frontal.png");
+        cv::imwrite(outputfile.string(), frontal_rendering);
+        outputfile.replace_extension("");
+
+        // And save the isomap:
+        outputfile.replace_extension(".isomap.png");
+        cv::imwrite(outputfile.string(), isomap);
+
+        // merge the isomaps:
+        merged_isomap = isomap_averaging.add_and_merge(isomap);
+    }
+
+    // save the merged isomap:
+    fs::path outputfile = outputfilebase;
+    outputfile +=fs::path("merged.isomap.png");
+    cv::imwrite(outputfile.string(), merged_isomap);
+    outputfile.replace_extension("");
+
+    // save the frontal rendering with merged isomap:
+    Mat frontal_rendering;
+    glm::mat4 modelview_frontal = glm::mat4( 1.0 );
+    core::Mesh neutral_expression = morphablemodel::sample_to_mesh(morphable_model.get_shape_model().draw_sample(pca_shape_coefficients), morphable_model.get_color_model().get_mean(), morphable_model.get_shape_model().get_triangle_list(), morphable_model.get_color_model().get_triangle_list(), morphable_model.get_texture_coordinates());
+    std::tie(frontal_rendering, std::ignore) = render::render(neutral_expression, modelview_frontal, glm::ortho(-130.0f, 130.0f, -130.0f, 130.0f), 512, 512, render::create_mipmapped_texture(merged_isomap), true, false, false);
+    outputfile.replace_extension(".frontal.png");
+    cv::imwrite(outputfile.string(), frontal_rendering);
+    outputfile.replace_extension("");
+
+    // Save the mesh as textured obj:
+    outputfile.replace_extension(".obj");
+    core::write_textured_obj(morphable_model.draw_sample(pca_shape_coefficients,std::vector<float>()), outputfile.string());
+
+    cout << "Finished fitting and wrote result mesh and isomap to files with basename " << outputfilebase << "." << endl;
+
+	return EXIT_SUCCESS;
+}

include/eos/fitting/fitting.hpp

@@ -190,7 +190,7 @@ inline auto concat(const std::vector<T>& vec_a, const std::vector<T>& vec_b)
 
 /**
  * @brief Fit the pose (camera), shape model, and expression blendshapes to landmarks,
- * in an iterative way.
+ * in an iterative way. Can fit to more than one set of landmarks, thus multiple images.
  *
  * Convenience function that fits pose (camera), the shape model, and expression blendshapes
  * to landmarks, in an iterative (alternating) way. It fits both sides of the face contour as well.
@@ -229,21 +229,27 @@ inline auto concat(const std::vector<T>& vec_a, const std::vector<T>& vec_b)
  * @param[out] fitted_image_points Debug parameter: Returns all the 2D points that have been used for the 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, boost::optional<int> num_shape_coefficients_to_fit, float lambda, boost::optional<fitting::RenderingParameters> initial_rendering_params, std::vector<float>& pca_shape_coefficients, std::vector<float>& blendshape_coefficients, std::vector<cv::Vec2f>& fitted_image_points)
+inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParameters>> fit_shape_and_pose_multi(const morphablemodel::MorphableModel& morphable_model, const std::vector<morphablemodel::Blendshape>& blendshapes, const std::vector<core::LandmarkCollection<cv::Vec2f>>& landmarks, const core::LandmarkMapper& landmark_mapper, std::vector<int> image_width, std::vector<int> image_height, const morphablemodel::EdgeTopology& edge_topology, const fitting::ContourLandmarks& contour_landmarks, const fitting::ModelContour& model_contour, int num_iterations, boost::optional<int> num_shape_coefficients_to_fit, float lambda, boost::optional<fitting::RenderingParameters> initial_rendering_params, std::vector<float>& pca_shape_coefficients, std::vector<std::vector<float>>& blendshape_coefficients, std::vector<std::vector<cv::Vec2f>>& fitted_image_points)
 {
     assert(blendshapes.size() > 0);
-	assert(landmarks.size() >= 4);
-	assert(image_width > 0 && image_height > 0);
+    assert(landmarks.size() > 0 && landmarks.size() == image_width.size() && image_width.size() == image_height.size());
     assert(num_iterations > 0); // Can we allow 0, for only the initial pose-fit?
     assert(pca_shape_coefficients.size() <= morphable_model.get_shape_model().get_num_principal_components());
+    int num_images = static_cast<int>(landmarks.size());
+    for (int j = 0; j < num_images; ++j) {
+        assert(landmarks[j].size() >= 4);
+        assert(image_width[j] > 0 && image_height[j] > 0);
+    }
     // More asserts I forgot?
 
 	using std::vector;
 	using cv::Vec2f;
 	using cv::Vec4f;
+	using cv::Mat;
 	using Eigen::VectorXf;
 	using Eigen::MatrixXf;
 
+
 	if (!num_shape_coefficients_to_fit)
 	{
 		num_shape_coefficients_to_fit = morphable_model.get_shape_model().get_num_principal_components();
@@ -258,111 +264,252 @@ inline std::pair<core::Mesh, fitting::RenderingParameters> fit_shape_and_pose(co
 
     if (blendshape_coefficients.empty())
     {
-		blendshape_coefficients.resize(blendshapes.size());
+        for (int j = 0; j < num_images; ++j) {
+            std::vector<float> current_blendshape_coefficients;
+            current_blendshape_coefficients.resize(blendshapes.size());
+            blendshape_coefficients.push_back(current_blendshape_coefficients);
+        }
     }
 
     MatrixXf blendshapes_as_basis = morphablemodel::to_matrix(blendshapes);
 
     // Current mesh - either from the given coefficients, or the mean:
     VectorXf current_pca_shape = morphable_model.get_shape_model().draw_sample(pca_shape_coefficients);
-	VectorXf current_combined_shape = current_pca_shape + blendshapes_as_basis * Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients.data(), blendshape_coefficients.size());
-	auto current_mesh = morphablemodel::sample_to_mesh(current_combined_shape, morphable_model.get_color_model().get_mean(), morphable_model.get_shape_model().get_triangle_list(), morphable_model.get_color_model().get_triangle_list(), morphable_model.get_texture_coordinates());
+    vector<VectorXf> current_combined_shapes;
+    vector<eos::core::Mesh> current_meshs;
+    for (int j = 0; j < num_images; ++j) {
+        VectorXf current_combined_shape = current_pca_shape + blendshapes_as_basis * Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients[j].data(), blendshape_coefficients[j].size());
+        current_combined_shapes.push_back(current_combined_shape);
+
+        eos::core::Mesh current_mesh = morphablemodel::sample_to_mesh(current_combined_shape, morphable_model.get_color_model().get_mean(), morphable_model.get_shape_model().get_triangle_list(), morphable_model.get_color_model().get_triangle_list(), morphable_model.get_texture_coordinates());
+        current_meshs.push_back(current_mesh);
+    }
 
     // The 2D and 3D point correspondences used for the fitting:
-	vector<Vec4f> model_points; // the points in the 3D shape model
-	vector<int> vertex_indices; // their vertex indices
-	vector<Vec2f> image_points; // the corresponding 2D landmark points
+    vector<vector<Vec4f>> model_points; // the points in the 3D shape model of all frames
+    vector<vector<int>> vertex_indices; // their vertex indices of all frames
+    vector<vector<Vec2f>> image_points; // the corresponding 2D landmark points of all frames
+
+    for (int j = 0; j < num_images; ++j) {
+        vector<Vec4f> current_model_points;
+        vector<int> current_vertex_indices;
+        vector<Vec2f> current_image_points;
 
         // Sub-select all the landmarks which we have a mapping for (i.e. that are defined in the 3DMM),
         // and get the corresponding model points (mean if given no initial coeffs, from the computed shape otherwise):
-	for (int i = 0; i < landmarks.size(); ++i) {
-		auto converted_name = landmark_mapper.convert(landmarks[i].name);
+        for (int i = 0; i < landmarks[j].size(); ++i) {
+            auto converted_name = landmark_mapper.convert(landmarks[j][i].name);
             if (!converted_name) { // no mapping defined for the current landmark
                 continue;
             }
             int vertex_idx = std::stoi(converted_name.get());
-		Vec4f vertex(current_mesh.vertices[vertex_idx].x, current_mesh.vertices[vertex_idx].y, current_mesh.vertices[vertex_idx].z, current_mesh.vertices[vertex_idx].w);
-		model_points.emplace_back(vertex);
-		vertex_indices.emplace_back(vertex_idx);
-		image_points.emplace_back(landmarks[i].coordinates);
+            Vec4f vertex(current_meshs[j].vertices[vertex_idx].x, current_meshs[j].vertices[vertex_idx].y, current_meshs[j].vertices[vertex_idx].z, current_meshs[j].vertices[vertex_idx].w);
+            current_model_points.emplace_back(vertex);
+            current_vertex_indices.emplace_back(vertex_idx);
+            current_image_points.emplace_back(landmarks[j][i].coordinates);
+        }
+
+        model_points.push_back(current_model_points);
+        vertex_indices.push_back(current_vertex_indices);
+        image_points.push_back(current_image_points);
     }
 
     // Need to do an initial pose fit to do the contour fitting inside the loop.
     // We'll do an expression fit too, since face shapes vary quite a lot, depending on expressions.
-	fitting::ScaledOrthoProjectionParameters current_pose;
-	current_pose = fitting::estimate_orthographic_projection_linear(image_points, model_points, true, image_height);
-	fitting::RenderingParameters rendering_params(current_pose, image_width, image_height);
+    vector<fitting::RenderingParameters> rendering_params;
+    for (int j = 0; j < num_images; ++j) {
+        fitting::ScaledOrthoProjectionParameters current_pose = fitting::estimate_orthographic_projection_linear(image_points[j], model_points[j], true, image_height[j]);
+        fitting::RenderingParameters current_rendering_params(current_pose, image_width[j], image_height[j]);
+        rendering_params.push_back(current_rendering_params);
 
-	cv::Mat affine_from_ortho = fitting::get_3x4_affine_camera_matrix(rendering_params, image_width, image_height);
-	blendshape_coefficients = fitting::fit_blendshapes_to_landmarks_nnls(blendshapes, current_pca_shape, affine_from_ortho, image_points, vertex_indices);
+        cv::Mat affine_from_ortho = fitting::get_3x4_affine_camera_matrix(current_rendering_params, image_width[j], image_height[j]);
+        blendshape_coefficients[j] = fitting::fit_blendshapes_to_landmarks_nnls(blendshapes, current_pca_shape, affine_from_ortho, image_points[j], vertex_indices[j]);
 
         // Mesh with same PCA coeffs as before, but new expression fit (this is relevant if no initial blendshape coeffs have been given):
-	current_combined_shape = current_pca_shape + morphablemodel::to_matrix(blendshapes) * Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients.data(), blendshape_coefficients.size());
-	current_mesh = morphablemodel::sample_to_mesh(current_combined_shape, morphable_model.get_color_model().get_mean(), morphable_model.get_shape_model().get_triangle_list(), morphable_model.get_color_model().get_triangle_list(), morphable_model.get_texture_coordinates());
+        current_combined_shapes[j] = current_pca_shape + morphablemodel::to_matrix(blendshapes) * Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients[j].data(),blendshape_coefficients[j].size());
+        current_meshs[j] = morphablemodel::sample_to_mesh(current_combined_shapes[j], morphable_model.get_color_model().get_mean(), morphable_model.get_shape_model().get_triangle_list(), morphable_model.get_color_model().get_triangle_list(), morphable_model.get_texture_coordinates());
+    }
 
     // The static (fixed) landmark correspondences which will stay the same throughout
     // the fitting (the inner face landmarks):
-	auto fixed_image_points = image_points;
-	auto fixed_vertex_indices = vertex_indices;
+    vector<vector<int>> fixed_vertex_indices (vertex_indices);
+    vector<vector<Vec2f>> fixed_image_points (image_points);
 
     for (int i = 0; i < num_iterations; ++i)
     {
+        std::vector<cv::Mat> affine_from_orthos;
+        std::vector<VectorXf> mean_plus_blendshapes;
+
         image_points = fixed_image_points;
         vertex_indices = fixed_vertex_indices;
+
+        for (int j = 0; j < num_images; ++j) {
+
             // Given the current pose, find 2D-3D contour correspondences of the front-facing face contour:
             vector<Vec2f> image_points_contour;
             vector<int> vertex_indices_contour;
-		auto yaw_angle = glm::degrees(glm::eulerAngles(rendering_params.get_rotation())[1]);
+            auto yaw_angle = glm::degrees(glm::eulerAngles(rendering_params[j].get_rotation())[1]);
             // For each 2D contour landmark, get the corresponding 3D vertex point and vertex id:
-		std::tie(image_points_contour, std::ignore, vertex_indices_contour) = fitting::get_contour_correspondences(landmarks, contour_landmarks, model_contour, yaw_angle, current_mesh, rendering_params.get_modelview(), rendering_params.get_projection(), fitting::get_opencv_viewport(image_width, image_height));
+            std::tie(image_points_contour, std::ignore, vertex_indices_contour) = fitting::get_contour_correspondences(landmarks[j], contour_landmarks, model_contour, yaw_angle, current_meshs[j], rendering_params[j].get_modelview(), rendering_params[j].get_projection(), fitting::get_opencv_viewport(image_width[j], image_height[j]));
             // Add the contour correspondences to the set of landmarks that we use for the fitting:
-		vertex_indices = fitting::concat(vertex_indices, vertex_indices_contour);
-		image_points = fitting::concat(image_points, image_points_contour);
+            vertex_indices[j] = fitting::concat(vertex_indices[j], vertex_indices_contour);
+            image_points[j] = fitting::concat(image_points[j], image_points_contour);
 
             // Fit the occluding (away-facing) contour using the detected contour LMs:
             vector<Eigen::Vector2f> occluding_contour_landmarks;
             if (yaw_angle >= 0.0f) // positive yaw = subject looking to the left
             { // the left contour is the occluding one we want to use ("away-facing")
-			auto contour_landmarks_ = core::filter(landmarks, contour_landmarks.left_contour); // Can do this outside of the loop
+                auto contour_landmarks_ = core::filter(landmarks[j], contour_landmarks.left_contour); // Can do this outside of the loop
                 std::for_each(begin(contour_landmarks_), end(contour_landmarks_), [&occluding_contour_landmarks](auto&& lm) { occluding_contour_landmarks.push_back({ lm.coordinates[0], lm.coordinates[1] }); });
             }
             else {
-			auto contour_landmarks_ = core::filter(landmarks, contour_landmarks.right_contour);
+                auto contour_landmarks_ = core::filter(landmarks[j], contour_landmarks.right_contour);
                 std::for_each(begin(contour_landmarks_), end(contour_landmarks_), [&occluding_contour_landmarks](auto&& lm) { occluding_contour_landmarks.push_back({ lm.coordinates[0], lm.coordinates[1] }); });
             }
-		auto edge_correspondences = fitting::find_occluding_edge_correspondences(current_mesh, edge_topology, rendering_params, occluding_contour_landmarks, 180.0f);
-		image_points = fitting::concat(image_points, edge_correspondences.first);
-		vertex_indices = fitting::concat(vertex_indices, edge_correspondences.second);
+            auto edge_correspondences = fitting::find_occluding_edge_correspondences(current_meshs[j], edge_topology, rendering_params[j], occluding_contour_landmarks, 180.0f);
+            image_points[j] = fitting::concat(image_points[j], edge_correspondences.first);
+            vertex_indices[j] = fitting::concat(vertex_indices[j], edge_correspondences.second);
 
             // Get the model points of the current mesh, for all correspondences that we've got:
-		model_points.clear();
-		for (const auto& v : vertex_indices)
+            model_points[j].clear();
+            for (const auto& v : vertex_indices[j])
             {
-			model_points.push_back({ current_mesh.vertices[v][0], current_mesh.vertices[v][1], current_mesh.vertices[v][2], current_mesh.vertices[v][3] });
+                model_points[j].push_back({ current_meshs[j].vertices[v][0], current_meshs[j].vertices[v][1], current_meshs[j].vertices[v][2], current_meshs[j].vertices[v][3] });
             }
 
             // Re-estimate the pose, using all correspondences:
-		current_pose = fitting::estimate_orthographic_projection_linear(image_points, model_points, true, image_height);
-		rendering_params = fitting::RenderingParameters(current_pose, image_width, image_height);
+            fitting::ScaledOrthoProjectionParameters current_pose = fitting::estimate_orthographic_projection_linear(image_points[j], model_points[j], true, image_height[j]);
+            rendering_params[j] = fitting::RenderingParameters(current_pose, image_width[j], image_height[j]);
 
-		cv::Mat affine_from_ortho = fitting::get_3x4_affine_camera_matrix(rendering_params, image_width, image_height);
+            cv::Mat affine_from_ortho = fitting::get_3x4_affine_camera_matrix(rendering_params[j], image_width[j], image_height[j]);
+            affine_from_orthos.push_back(affine_from_ortho);
 
             // Estimate the PCA shape coefficients with the current blendshape coefficients:
-		VectorXf mean_plus_blendshapes = morphable_model.get_shape_model().get_mean() + blendshapes_as_basis * Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients.data(), blendshape_coefficients.size());
-		pca_shape_coefficients = fitting::fit_shape_to_landmarks_linear(morphable_model, affine_from_ortho, image_points, vertex_indices, mean_plus_blendshapes, lambda, num_shape_coefficients_to_fit);
+            VectorXf current_mean_plus_blendshapes = morphable_model.get_shape_model().get_mean() + blendshapes_as_basis * Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients[j].data(),blendshape_coefficients[j].size());
+            mean_plus_blendshapes.push_back(current_mean_plus_blendshapes);
+        }
+        pca_shape_coefficients = fitting::fit_shape_to_landmarks_linear_multi(morphable_model, affine_from_orthos, image_points, vertex_indices, mean_plus_blendshapes, lambda, num_shape_coefficients_to_fit);
 
         // Estimate the blendshape coefficients with the current PCA model estimate:
         current_pca_shape = morphable_model.get_shape_model().draw_sample(pca_shape_coefficients);
-		blendshape_coefficients = fitting::fit_blendshapes_to_landmarks_nnls(blendshapes, current_pca_shape, affine_from_ortho, image_points, vertex_indices);
 
-		current_combined_shape = current_pca_shape + blendshapes_as_basis * Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients.data(), blendshape_coefficients.size());
-		current_mesh = morphablemodel::sample_to_mesh(current_combined_shape, morphable_model.get_color_model().get_mean(), morphable_model.get_shape_model().get_triangle_list(), morphable_model.get_color_model().get_triangle_list(), morphable_model.get_texture_coordinates());
+        for (int j = 0; j < num_images; ++j) {
+            blendshape_coefficients[j] = fitting::fit_blendshapes_to_landmarks_nnls(blendshapes, current_pca_shape, affine_from_orthos[j], image_points[j], vertex_indices[j]);
+
+            current_combined_shapes[j] = current_pca_shape + blendshapes_as_basis * Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients[j].data(),blendshape_coefficients[j].size());
+            current_meshs[j] = morphablemodel::sample_to_mesh(current_combined_shapes[j], morphable_model.get_color_model().get_mean(), morphable_model.get_shape_model().get_triangle_list(), morphable_model.get_color_model().get_triangle_list(), morphable_model.get_texture_coordinates());
+        }
     }
 
     fitted_image_points = image_points;
-	return { current_mesh, rendering_params }; // I think we could also work with a Mat face_instance in this function instead of a Mesh, but it would convolute the code more (i.e. more complicated to access vertices).
+    return { current_meshs, rendering_params }; // I think we could also work with a Mat face_instance in this function instead of a Mesh, but it would convolute the code more (i.e. more complicated to access vertices).
 };
 
+/**
+ * @brief Fit the pose (camera), shape model, and expression blendshapes to landmarks,
+ * in an iterative way.  Can fit to more than one set of landmarks, thus multiple images.
+ *
+ * Convenience function that fits pose (camera), the shape model, and expression blendshapes
+ * to landmarks, in an iterative (alternating) way. It fits both sides of the face contour as well.
+ *
+ * If you want to access the values of shape or blendshape coefficients, or want to set starting
+ * values for them, use the following overload to this function:
+ * std::pair<render::Mesh, fitting::RenderingParameters> fit_shape_and_pose(const morphablemodel::MorphableModel&, const std::vector<morphablemodel::Blendshape>&, const core::LandmarkCollection<cv::Vec2f>&, const core::LandmarkMapper&, int, int, const morphablemodel::EdgeTopology&, const fitting::ContourLandmarks&, const fitting::ModelContour&, int, boost::optional<int>, float, boost::optional<fitting::RenderingParameters>, std::vector<float>&, std::vector<float>&, std::vector<cv::Vec2f>&)
+ *
+ * Todo: Add a convergence criterion.
+ *
+ * \p num_iterations: Results are good for even a single iteration. For single-image fitting and
+ * for full convergence of all parameters, it can take up to 300 iterations. In tracking,
+ * particularly if initialising with the previous frame, it works well with as low as 1 to 5
+ * iterations.
+ * \p edge_topology is used for the occluding-edge face contour fitting.
+ * \p contour_landmarks and \p model_contour are used to fit the front-facing contour.
+ *
+ * @param[in] morphable_model The 3D Morphable Model used for the shape fitting.
+ * @param[in] blendshapes A vector of blendshapes that are being fit to the landmarks in addition to the PCA model.
+ * @param[in] landmarks 2D landmarks from an image to fit the model to.
+ * @param[in] landmark_mapper Mapping info from the 2D landmark points to 3D vertex indices.
+ * @param[in] image_width Width of the input image (needed for the camera model).
+ * @param[in] image_height Height of the input image (needed for the camera model).
+ * @param[in] edge_topology Precomputed edge topology of the 3D model, needed for fast edge-lookup.
+ * @param[in] contour_landmarks 2D image contour ids of left or right side (for example for ibug landmarks).
+ * @param[in] model_contour The model contour indices that should be considered to find the closest corresponding 3D vertex.
+ * @param[in] num_iterations Number of iterations that the different fitting parts will be alternated for.
+ * @param[in] num_shape_coefficients_to_fit How many shape-coefficients to fit (all others will stay 0). Should be bigger than zero, or boost::none to fit all coefficients.
+ * @param[in] lambda Regularisation parameter of the PCA shape fitting.
+ * @return The fitted model shape instance and the final pose.
+ */
+inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParameters>> fit_shape_and_pose_multi(const morphablemodel::MorphableModel& morphable_model, const std::vector<morphablemodel::Blendshape>& blendshapes, const std::vector<core::LandmarkCollection<cv::Vec2f>>& landmarks, const core::LandmarkMapper& landmark_mapper, std::vector<int> image_width, std::vector<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)
+{
+    std::vector<float> pca_shape_coefficients;
+    std::vector<std::vector<float>> blendshape_coefficients;
+    std::vector<std::vector<cv::Vec2f>> fitted_image_points;
+
+    return fit_shape_and_pose_multi(morphable_model, blendshapes, landmarks, landmark_mapper, image_width, image_height, edge_topology, contour_landmarks, model_contour, num_iterations, num_shape_coefficients_to_fit, lambda, boost::none, pca_shape_coefficients, blendshape_coefficients, fitted_image_points);
+};
+
+
+/**
+ * @brief Fit the pose (camera), shape model, and expression blendshapes to landmarks,
+ * in an iterative way.
+ *
+ * Convenience function that fits pose (camera), the shape model, and expression blendshapes
+ * to landmarks, in an iterative (alternating) way. It fits both sides of the face contour as well.
+ *
+ * If \p pca_shape_coefficients and/or \p blendshape_coefficients are given, they are used as
+ * starting values in the fitting. When the function returns, they contain the coefficients from
+ * the last iteration.
+ *
+ * Use render::Mesh fit_shape_and_pose(const morphablemodel::MorphableModel&, const std::vector<morphablemodel::Blendshape>&, const core::LandmarkCollection<cv::Vec2f>&, const core::LandmarkMapper&, int, int, const morphablemodel::EdgeTopology&, const fitting::ContourLandmarks&, const fitting::ModelContour&, int, boost::optional<int>, float).
+ * for a simpler overload with reasonable defaults and no optional output.
+ *
+ * \p num_iterations: Results are good for even a single iteration. For single-image fitting and
+ * for full convergence of all parameters, it can take up to 300 iterations. In tracking,
+ * particularly if initialising with the previous frame, it works well with as low as 1 to 5
+ * iterations.
+ * \p edge_topology is used for the occluding-edge face contour fitting.
+ * \p contour_landmarks and \p model_contour are used to fit the front-facing contour.
+ *
+ * Todo: Add a convergence criterion.
+ *
+ * @param[in] morphable_model The 3D Morphable Model used for the shape fitting.
+ * @param[in] blendshapes A vector of blendshapes that are being fit to the landmarks in addition to the PCA model.
+ * @param[in] landmarks 2D landmarks from an image to fit the model to.
+ * @param[in] landmark_mapper Mapping info from the 2D landmark points to 3D vertex indices.
+ * @param[in] image_width Width of the input image (needed for the camera model).
+ * @param[in] image_height Height of the input image (needed for the camera model).
+ * @param[in] edge_topology Precomputed edge topology of the 3D model, needed for fast edge-lookup.
+ * @param[in] contour_landmarks 2D image contour ids of left or right side (for example for ibug landmarks).
+ * @param[in] model_contour The model contour indices that should be considered to find the closest corresponding 3D vertex.
+ * @param[in] num_iterations Number of iterations that the different fitting parts will be alternated for.
+ * @param[in] num_shape_coefficients_to_fit How many shape-coefficients to fit (all others will stay 0). Should be bigger than zero, or boost::none to fit all coefficients.
+ * @param[in] lambda Regularisation parameter of the PCA shape fitting.
+ * @param[in] initial_rendering_params Currently ignored (not used).
+ * @param[in,out] pca_shape_coefficients If given, will be used as initial PCA shape coefficients to start the fitting. Will contain the final estimated coefficients.
+ * @param[in,out] blendshape_coefficients If given, will be used as initial expression blendshape coefficients to start the fitting. Will contain the final estimated coefficients.
+ * @param[out] fitted_image_points Debug parameter: Returns all the 2D points that have been used for the 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, boost::optional<int> num_shape_coefficients_to_fit, float lambda, boost::optional<fitting::RenderingParameters> initial_rendering_params, std::vector<float>& pca_shape_coefficients, std::vector<float>& blendshape_coefficients, std::vector<cv::Vec2f>& fitted_image_points)
+{
+    //we have to create a new vector here if the blendshape_coefficients are empty, as otherwise the new vector is not empty anymore and contains one element
+    std::vector<std::vector<float>> all_blendshape_coefficients;
+    if (!blendshape_coefficients.empty())
+    {
+        all_blendshape_coefficients = {blendshape_coefficients};
+    }
+    std::vector<std::vector<cv::Vec2f>> all_fitted_image_points;
+    if (!fitted_image_points.empty())
+    {
+        all_fitted_image_points = {fitted_image_points};
+    }
+    std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParameters>>  all_meshs_and_params = fit_shape_and_pose_multi( morphable_model, blendshapes, { landmarks }, landmark_mapper, { image_width }, { image_height }, edge_topology, contour_landmarks, model_contour, num_iterations, num_shape_coefficients_to_fit, lambda, initial_rendering_params, pca_shape_coefficients, all_blendshape_coefficients, all_fitted_image_points);
+    blendshape_coefficients = all_blendshape_coefficients[0];
+    fitted_image_points = all_fitted_image_points[0];
+
+    return {all_meshs_and_params.first[0], all_meshs_and_params.second[0] };
+}
+
 /**
  * @brief Fit the pose (camera), shape model, and expression blendshapes to landmarks,
  * in an iterative way.

include/eos/fitting/linear_shape_fitting.hpp

@@ -184,42 +184,40 @@ inline std::vector<float> fit_shape_to_landmarks_linear_multi(const morphablemod
  * @param[in] landmarks 2D landmarks from an image to fit the model to.
  * @param[in] vertex_ids The vertex ids in the model that correspond to the 2D points.
  * @param[in] base_face The base or reference face from where the fitting is started. Usually this would be the models mean face, which is what will be used if the parameter is not explicitly specified.
- * @param[in] lambda The regularisation parameter (weight of the prior towards the mean).
+ * @param[in] lambda The regularisation parameter (weight of the prior towards the mean). Gets normalized by the number of images given.
  * @param[in] num_coefficients_to_fit How many shape-coefficients to fit (all others will stay 0). Should be bigger than zero, or boost::none to fit all coefficients.
  * @param[in] detector_standard_deviation The standard deviation of the 2D landmarks given (e.g. of the detector used), in pixels.
  * @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(const morphablemodel::MorphableModel& morphable_model, cv::Mat affine_camera_matrix, const std::vector<cv::Vec2f>& landmarks, const std::vector<int>& vertex_ids, Eigen::VectorXf base_face=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>())
 {
+
+	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_landmarks = static_cast<int>(landmarks.size());
+	int num_images = affine_camera_matrix.size();
 
-	if (base_face.size() == 0)
-	{
-		base_face = morphable_model.get_shape_model().get_mean();
+    // the regularisation has to be adjusted when more than one image is given
+    lambda *= num_images;
+
+	int total_num_landmarks_dimension = 0;
+	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}$:
-	MatrixXf V_hat_h = MatrixXf::Zero(4 * num_landmarks, num_coeffs_to_fit);
-	int row_index = 0;
-	for (int i = 0; i < num_landmarks; ++i) {
-		MatrixXf basis_rows_ = morphable_model.get_shape_model().get_rescaled_pca_basis_at_point(vertex_ids[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(row_index, 0, 3, V_hat_h.cols()) = basis_rows_.block(0, 0, basis_rows_.rows(), num_coeffs_to_fit);
-		row_index += 4; // replace 3 rows and skip the 4th one, it has all zeros
-	}
+    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:
-	MatrixXf P = MatrixXf::Zero(3 * num_landmarks, 4 * num_landmarks);
-	for (int i = 0; i < num_landmarks; ++i) {
-		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);
-	}
+    Eigen::SparseMatrix<float> P (3 * total_num_landmarks_dimension, 4 * total_num_landmarks_dimension);
+    std::vector<Eigen::Triplet<float>> P_coefficients; // list of non-zeros coefficients
+	int P_index = 0;
 	// 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.
@@ -227,26 +225,69 @@ inline std::vector<float> fit_shape_to_landmarks_linear(const morphablemodel::Mo
 	// 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 * num_landmarks, 1.0f / sigma_squared_2D);
-	// Earlier, we set Sigma in a for-loop and then computed Omega, but it was really unnecessary:
-	// 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.
+	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$
+	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
+	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.
+
+    for (int k = 0; k < num_images; ++k)
+    {
+        // For each image we have, set up the equations and add it to the matrices:
+        assert(landmarks[k].size() == vertex_ids[k].size()); // has to be valid for each img
+
+        int num_landmarks = static_cast<int>(landmarks[k].size());
+
+        if (base_face[k].size()==0)
+        {
+            base_face[k] = morphable_model.get_shape_model().get_mean();
+        }
+
+        // $\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_coeffs_to_fit, CV_32FC1);
+        for (int i = 0; i < num_landmarks; ++i) {
+            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:
+        for (int i = 0; i < num_landmarks; ++i) {
+            for (int x = 0; x < affine_camera_matrix[k].rows; ++x){
+                for (int y = 0; y < affine_camera_matrix[k].cols; ++y){
+                    P_coefficients.push_back(Eigen::Triplet<float>(3 * P_index + x, 4 * P_index + y, affine_camera_matrix[k].at<float>(x,y)));
+                }
+            }
+            ++P_index;
+        }
+
 
         // The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
-	VectorXf y = VectorXf::Ones(3 * num_landmarks);
+        //Mat y = Mat::ones(3 * num_landmarks, 1, CV_32FC1);
         for (int i = 0; i < num_landmarks; ++i) {
-		y(3 * i) = landmarks[i][0];
-		y((3 * i) + 1) = landmarks[i][1];
+            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
-	VectorXf v_bar = VectorXf::Ones(4 * num_landmarks);
+        //Mat v_bar = Mat::ones(4 * num_landmarks, 1, CV_32FC1);
         for (int i = 0; i < num_landmarks; ++i) {            
-		v_bar(4 * i) = base_face(vertex_ids[i] * 3);
-		v_bar((4 * i) + 1) = base_face(vertex_ids[i] * 3 + 1);
-		v_bar((4 * i) + 2) = base_face(vertex_ids[i] * 3 + 2);
-		//v_bar((4 * i) + 3) = 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;
+        }
     }
 
+    // fill P with coefficients
+    P.setFromTriplets(P_coefficients.begin(), P_coefficients.end());
+
 	// 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
@@ -259,6 +300,35 @@ inline std::vector<float> fit_shape_to_landmarks_linear(const morphablemodel::Mo
     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).
+ *
+ * [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 model to screen-space (should probably be of type CV_32FC1 as all our calculations are done with float).
+ * @param[in] landmarks 2D landmarks from an image to fit the model to.
+ * @param[in] vertex_ids The vertex ids in the model that correspond to the 2D points.
+ * @param[in] base_face The base or reference face from where the fitting is started. Usually this would be the models mean face, which is what will be used if the parameter is not explicitly specified.
+ * @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). Should be bigger than zero, or boost::none to fit all coefficients.
+ * @param[in] detector_standard_deviation The standard deviation of the 2D landmarks given (e.g. of the detector used), in pixels.
+ * @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(const morphablemodel::MorphableModel& morphable_model, cv::Mat affine_camera_matrix, std::vector<cv::Vec2f> landmarks, std::vector<int> vertex_ids, Eigen::VectorXf base_face=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>())
+{
+    std::vector<std::vector<cv::Vec2f>> all_landmarks = {landmarks};
+    std::vector<std::vector<int>> all_vertex_ids = {vertex_ids};
+    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/fitting/multiframe_linear.hpp

+/*
+ * eos - A 3D Morphable Model fitting library written in modern C++11/14.
+ *
+ * File: include/eos/fitting/multiframe_linear.hpp
+ *
+ * Copyright 2016 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 MULTIFRAME_LINEAR_PATRIK_HPP_
+#define MULTIFRAME_LINEAR_PATRIK_HPP_
+
+#include "eos/morphablemodel/MorphableModel.hpp"
+#include "eos/morphablemodel/Blendshape.hpp"
+#include "eos/fitting/linear_shape_fitting.hpp"
+#include "eos/fitting/blendshape_fitting.hpp"
+#include "eos/fitting/contour_correspondence.hpp"
+#include "eos/fitting/nonlinear_camera_estimation.hpp"
+
+//#include "glm/gtc/matrix_transform.hpp"
+
+#include "opencv2/core/core.hpp"
+
+#include <vector>
+#include <iostream>
+
+namespace eos {
+	namespace fitting {
+		/**
+		 *
+		 * @param morphable_model
+		 * @param affine_camera_matrix
+		 * @param landmarks
+		 * @param vertex_ids
+		 * @param base_face
+		 * @param lambda
+		 * @param num_coefficients_to_fit
+		 * @param detector_standard_deviation
+		 * @param model_standard_deviation
+		 *
+		 * @return
+		 */
+		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 <cv::Mat> base_face = std::vector<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>()) {
+		}
+	}
+}

include/eos/fitting/multiframe_linear_patrik.hpp

+/*
+ * eos - A 3D Morphable Model fitting library written in modern C++11/14.
+ *
+ * File: include/eos/fitting/multiframe_linear.hpp
+ *
+ * Copyright 2016 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 MULTIFRAME_LINEAR_PATRIK_HPP_
+#define MULTIFRAME_LINEAR_PATRIK_HPP_
+
+#include "eos/morphablemodel/MorphableModel.hpp"
+#include "eos/morphablemodel/Blendshape.hpp"
+#include "eos/fitting/linear_shape_fitting.hpp"
+#include "eos/fitting/blendshape_fitting.hpp"
+#include "eos/fitting/contour_correspondence.hpp"
+#include "eos/fitting/nonlinear_camera_estimation.hpp"
+
+//#include "glm/gtc/matrix_transform.hpp"
+
+#include "opencv2/core/core.hpp"
+
+#include <vector>
+#include <iostream>
+
+namespace eos {
+	namespace fitting {
+
+// I think this function can be split in half - the setup, and the actual solving of the system? And then reuse parts of it in _multi?
+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<cv::Mat> base_face=std::vector<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>())
+{
+	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
+
+	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();
+
+	int total_num_landmarks_dimension = 0;
+	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);
+	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);
+	int P_index = 0;
+	Mat Sigma = Mat::zeros(3 * total_num_landmarks_dimension, 3 * total_num_landmarks_dimension, CV_32FC1);
+	int Sigma_index = 0; // this runs the same as P_index
+	Mat Omega;
+	// The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
+	Mat y = Mat::ones(3 * total_num_landmarks_dimension, 1, CV_32FC1);
+	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);
+	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.
+
+	for (int k = 0; k < num_images; ++k)
+	{
+		// For each image we have, set up the equations and add it to the matrices:
+		assert(landmarks[k].size() == vertex_ids[k].size()); // has to be valid for each img
+		
+		int num_landmarks = static_cast<int>(landmarks[k].size());
+
+		if (base_face[k].empty())
+		{
+			base_face[k] = morphable_model.get_shape_model().get_mean();
+		}
+
+		// $\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_coeffs_to_fit, CV_32FC1);
+		for (int i = 0; i < num_landmarks; ++i) {
+			Mat basis_rows = morphable_model.get_shape_model().get_normalised_pca_basis(vertex_ids[k][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(V_hat_h_row_index, V_hat_h_row_index + 3));
+			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);
+			++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.at<float>(Sigma_index, Sigma_index) = 1.0f / std::sqrt(sigma_squared_2D); // the higher the sigma_squared_2D, the smaller the diagonal entries of Sigma will be
+			++Sigma_index;
+		}
+		//Mat Omega = Sigma.t() * Sigma; // just squares the diagonal
+		// => moved outside the loop
+		
+		// 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_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].at<float>(vertex_ids[k][i] * 3), base_face[k].at<float>(vertex_ids[k][i] * 3 + 1), base_face[k].at<float>(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_index;
+			// note: now that a Vec4f is returned, we could use copyTo?
+		}
+	}
+	Omega = Sigma.t() * Sigma; // moved outside the loop. But can do even more efficiently anyway.
+
+	// 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>.
+
+/*	using RowMajorMatrixXf = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+	Eigen::Map<RowMajorMatrixXf> P_Eigen(P.ptr<float>(), P.rows, P.cols);
+	Eigen::Map<RowMajorMatrixXf> V_hat_h_Eigen(V_hat_h.ptr<float>(), V_hat_h.rows, V_hat_h.cols);
+	Eigen::Map<RowMajorMatrixXf> v_bar_Eigen(v_bar.ptr<float>(), v_bar.rows, v_bar.cols);
+	Eigen::Map<RowMajorMatrixXf> y_Eigen(y.ptr<float>(), y.rows, y.cols);
+	Eigen::MatrixXf A_Eigen = P_Eigen * V_hat_h_Eigen;
+	Eigen::MatrixXf res = (A_Eigen.transpose() * A_Eigen + lambda * Eigen::MatrixXf::Identity(num_coeffs_to_fit, num_coeffs_to_fit)).householderQr().solve(-A_Eigen.transpose() * (P_Eigen * v_bar_Eigen - y_Eigen));
+	Mat c_s(res.rows(), res.cols(), CV_32FC1, res.data()); // Note: Could also use LDLT.
+*/
+	return std::vector<float>(c_s);
+};
+
+std::vector<cv::Mat> fit_shape_multi(std::vector<cv::Mat> affine_camera_matrix, eos::morphablemodel::MorphableModel morphable_model, std::vector<eos::morphablemodel::Blendshape> blendshapes, std::vector<std::vector<cv::Vec2f>> image_points, std::vector<std::vector<int>> vertex_indices, float lambda, boost::optional<int> num_coefficients_to_fit, std::vector<float>& pca_shape_coefficients, std::vector<std::vector<float>>& blendshape_coefficients)
+{
+	using cv::Mat;
+	using std::vector;
+	// asserts... copy from fit_shape_to_landmarks_linear_multi?
+	int num_images = affine_camera_matrix.size();
+	
+	Mat blendshapes_as_basis(blendshapes[0].deformation.rows, blendshapes.size(), CV_32FC1); // assert blendshapes.size() > 0 and all of them have same number of rows, and 1 col
+	for (int i = 0; i < blendshapes.size(); ++i)
+	{
+		blendshapes[i].deformation.copyTo(blendshapes_as_basis.col(i));
+	}
+
+	vector<vector<float>> last_blendshape_coeffs, current_blendshape_coeffs; // this might be vector<vector<float>>? One for each frame?
+	vector<float> last_pca_coeffs, current_pca_coeffs; // this will stay ONE single vector
+	current_blendshape_coeffs.resize(num_images);
+	for (auto&& cbc : current_blendshape_coeffs) {
+		cbc.resize(blendshapes.size()); // starting values t_0, all zeros
+	}
+	// no starting values for current_pca_coeffs required, since we start with the shape fitting, and cv::norm of an empty vector is 0.
+	vector<Mat> combined_shape;
+	combined_shape.resize(num_images);
+
+	double bs_error = std::numeric_limits<double>::max();
+	do // run at least once:
+	{
+		last_blendshape_coeffs = current_blendshape_coeffs;
+		last_pca_coeffs = current_pca_coeffs;
+		// Estimate the PCA shape coefficients with the current blendshape coefficients (0 in the first iteration):
+		vector<Mat> means_plus_blendshapes;
+		for (int img = 0; img < num_images; ++img)
+		{
+			Mat mean_plus_blendshapes = morphable_model.get_shape_model().get_mean() + blendshapes_as_basis * Mat(last_blendshape_coeffs[img]);
+			means_plus_blendshapes.push_back(mean_plus_blendshapes);
+		}
+		current_pca_coeffs = fit_shape_to_landmarks_linear_multi(morphable_model, affine_camera_matrix, image_points, vertex_indices, means_plus_blendshapes, lambda, num_coefficients_to_fit);
+
+		// Estimate the blendshape coefficients with the current PCA model estimate:
+		Mat pca_model_shape = morphable_model.get_shape_model().draw_sample(current_pca_coeffs);
+		// Multi: actually, the expression fitting we keep like that: Use the current global shape and estimate expression BS for each frame separately. There's no use in doing that jointly for all frames.
+		for (int img = 0; img < num_images; ++img)
+		{
+			current_blendshape_coeffs[img] = eos::fitting::fit_blendshapes_to_landmarks_nnls(blendshapes, pca_model_shape, affine_camera_matrix[img], image_points[img], vertex_indices[img], 0.0f);
+		}
+
+		for (int img = 0; img < num_images; ++img) { // just for debug purposes here.
+			combined_shape[img] = pca_model_shape + blendshapes_as_basis * Mat(current_blendshape_coeffs[img]); // Note/Todo: There's actually no need to do this in the loop here? We can just do it once, at the end?
+		}
+
+		vector<double> bs_diffs; // the abs error
+		for (int img = 0; img < num_images; ++img) {
+			auto nc = cv::norm(current_blendshape_coeffs[img]);
+			auto nl = cv::norm(last_blendshape_coeffs[img]);
+			auto d = std::abs(nc - nl);
+			bs_diffs.push_back(d);
+		}
+		double minv = 0; double maxv = 0;
+		cv::minMaxLoc(bs_diffs, &minv, &maxv);
+		bs_error = maxv;
+	} while (std::abs(cv::norm(current_pca_coeffs) - cv::norm(last_pca_coeffs)) >= 0.01 || bs_error >= 0.01);
+	// have to change this loop: norm over all frames or something like that? (for pca_coeffs, it doesn't change - only one set!)
+
+	pca_shape_coefficients = current_pca_coeffs;
+	blendshape_coefficients = current_blendshape_coeffs;
+	return combined_shape;
+};
+
+	} /* namespace fitting */
+} /* namespace eos */
+
+#endif /* MULTIFRAME_LINEAR_HPP_ */

include/eos/morphablemodel/MorphableModel.hpp

@@ -26,9 +26,10 @@
 
 #include "eos/core/Mesh.hpp"
 
-#include "eos/morphablemodel/io/mat_cerealisation.hpp"
+#include "eos/morphablemodel/io/eigen_cerealisation.hpp"
 #include "cereal/cereal.hpp"
 #include "cereal/access.hpp"
+#include "cereal/types/array.hpp"
 #include "cereal/types/vector.hpp"
 #include "cereal/archives/binary.hpp"
 
@@ -41,6 +42,7 @@
 #include <vector>
 #include <array>
 #include <cstdint>
+#include <fstream>
 
 namespace eos {
 	namespace morphablemodel {

include/eos/morphablemodel/PcaModel.hpp

@@ -26,6 +26,7 @@
 #include "cereal/access.hpp"
 #include "cereal/types/array.hpp"
 #include "cereal/types/vector.hpp"
+#include "cereal/archives/binary.hpp"
 
 #include "Eigen/Core"
 
@@ -34,6 +35,7 @@
 #include <array>
 #include <random>
 #include <cassert>
+#include <fstream>
 
 namespace eos {
 	namespace morphablemodel {
@@ -57,7 +59,7 @@ public:
 	PcaModel() = default;
 
 	/**
-	 * Construct a PCA model from given mean, normalised PCA basis, eigenvalues
+	 * Construct a PCA model from given mean, orthonormal PCA basis, eigenvalues
 	 * and triangle list.
 	 *
 	 * See the documentation of the member variables for how the data should
@@ -292,6 +294,45 @@ private:
 	};
 };
 
+/**
+ * Helper method to load a PCA model from
+ * a cereal::BinaryInputArchive from the harddisk.
+ *
+ * Usually, morphablemodel::load_model(std::string) should be used to directly
+ * load a MorphableModel. This function can be useful when it is necessary to just
+ * load a PCA model.
+ *
+ * @param[in] filename Filename to a model.
+ * @return The loaded PCA model.
+ * @throw std::runtime_error When the file given in \c filename fails to be opened (most likely because the file doesn't exist).
+ */
+inline PcaModel load_pca_model(std::string filename)
+{
+	PcaModel model;
+
+	std::ifstream file(filename, std::ios::binary);
+	if (file.fail()) {
+		throw std::runtime_error("Error opening given file: " + filename);
+	}
+	cereal::BinaryInputArchive input_archive(file);
+	input_archive(model);
+
+	return model;
+};
+
+/**
+ * Helper method to save a PCA model to the
+ * harddrive as cereal::BinaryOutputArchive.
+ *
+ * @param[in] model The model to be saved.
+ * @param[in] filename Filename for the model.
+ */
+inline void save_pca_model(PcaModel model, std::string filename)
+{
+	std::ofstream file(filename, std::ios::binary);
+	cereal::BinaryOutputArchive output_archive(file);
+	output_archive(model);
+};
 
 /**
  * Takes an orthonormal PCA basis matrix (a matrix consisting

python/generate-python-bindings.cpp

@@ -81,8 +81,8 @@ PYBIND11_PLUGIN(eos) {
 	 * Bindings for the eos::morphablemodel namespace:
 	 *  - PcaModel
 	 *  - MorphableModel
-	 *  - load_model()
-	 *  - save_model()
+	 *  - load_model(), save_model()
+	 *  - load_pca_model(), save_pca_model()
 	 */
 	py::module morphablemodel_module = eos_module.def_submodule("morphablemodel", "Functionality to represent a Morphable Model, its PCA models, and functions to load models and blendshapes.");
 
@@ -107,10 +107,13 @@ PYBIND11_PLUGIN(eos) {
 		.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.")
+		.def("get_texture_coordinates", &morphablemodel::MorphableModel::get_texture_coordinates, "Returns the texture coordinates for all the vertices in the model.")
 		;
 
 	morphablemodel_module.def("load_model", &morphablemodel::load_model, "Load a Morphable Model from a cereal::BinaryInputArchive (.bin) from the harddisk.", py::arg("filename"));
 	morphablemodel_module.def("save_model", &morphablemodel::save_model, "Save a Morphable Model as cereal::BinaryOutputArchive.", py::arg("model"), py::arg("filename"));
+	morphablemodel_module.def("load_pca_model", &morphablemodel::load_pca_model, "Load a PCA model from a cereal::BinaryInputArchive (.bin) from the harddisk.", py::arg("filename"));
+	morphablemodel_module.def("save_pca_model", &morphablemodel::save_pca_model, "Save a PCA model as cereal::BinaryOutputArchive.", py::arg("model"), py::arg("filename"));
 
 	/**
 	 *  - Blendshape

utils/CMakeLists.txt

@@ -39,7 +39,7 @@ add_executable(accuracy-evaluation accuracy-evaluation.cpp)
 target_link_libraries(accuracy-evaluation eos ${OpenCV_LIBS} ${Boost_LIBRARIES})
 target_link_libraries(accuracy-evaluation "$<$<CXX_COMPILER_ID:GNU>:-pthread>$<$<CXX_COMPILER_ID:Clang>:-pthreads>")
 
-# 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)