Add various OpenMP pragmas, including changes to make it possible - speedup is gained

CMakeLists.txt

@@ -113,8 +113,17 @@ set(HEADERS
   ${CMAKE_CURRENT_SOURCE_DIR}/include/eos/video/Keyframe.hpp
 )
 
+find_package(OpenMP)
+if (OPENMP_FOUND)
+  set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+  set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+endif()
+
+message(STATUS ${OpenMP_CXX_FLAGS})
+message(STATUS ${OpenMP_INCLUDE_DIR})
 add_library(eos INTERFACE)
 target_compile_features(eos INTERFACE ${EOS_CXX_COMPILE_FEATURES})
+
 # Add header includes:
 target_include_directories(eos INTERFACE "include")
 target_include_directories(eos INTERFACE ${CEREAL_INCLUDE_DIR})
@@ -125,6 +134,7 @@ target_include_directories(eos INTERFACE ${glm_INCLUDE_DIR})
 target_include_directories(eos INTERFACE ${nanoflann_INCLUDE_DIR})
 target_include_directories(eos INTERFACE ${eigen3_nnls_INCLUDE_DIR})
 
+
 # Custom target for the library, to make the headers show up in IDEs:
 add_custom_target(eos-headers SOURCES ${HEADERS})
 source_group(core REGULAR_EXPRESSION include/eos/core/*)

include/eos/core/landmark_utils.hpp

@@ -8,7 +8,7 @@
 #include "opencv2/core/core.hpp"
 #include "opencv2/opencv.hpp"
 
-#include "eos/morphablemodel/morphablemodel.hpp"
+#include "eos/morphablemodel/MorphableModel.hpp"
 #include "eos/core/LandmarkMapper.hpp"
 #include "eos/core/Landmark.hpp"
 
@@ -256,15 +256,13 @@ namespace eos {
 		 * @param[in,out] image_points
 		 */
 		template <typename vec2f, typename vec4f>
-		inline void get_landmark_coordinates(core::LandmarkCollection<vec2f> landmarks,
-										const core::LandmarkMapper& landmark_mapper,
-										eos::core::Mesh& mesh,
-										vector<vector<vec4f>>& model_points,
-										vector<vector<int>>& vertex_indices,
-										vector<vector<vec2f>>& image_points) {
-			vector<vec4f> current_model_points;
-			vector<int> current_vertex_indices;
-			vector<vec2f> current_image_points;
+		std::tuple<vector<vec4f>, vector<int>, vector<vec2f>> get_landmark_coordinates(
+				core::LandmarkCollection<vec2f> landmarks,
+				const core::LandmarkMapper& landmark_mapper,
+				eos::core::Mesh& mesh) {
+			vector<vec4f> model_points;
+			vector<int> vertex_indices;
+			vector<vec2f> image_points;
 
 			for (auto &lm: landmarks) {
 				auto converted_name = landmark_mapper.convert(lm.name);
@@ -280,15 +278,62 @@ namespace eos {
 					mesh.vertices[vertex_idx].w
 				);
 
-				current_model_points.emplace_back(vertex);
-				current_vertex_indices.emplace_back(vertex_idx);
-				current_image_points.emplace_back(lm.coordinates);
+				model_points.emplace_back(vertex);
+				vertex_indices.emplace_back(vertex_idx);
+				image_points.emplace_back(lm.coordinates);
+			}
+
+			return std::make_tuple(model_points, vertex_indices, image_points);
+		}
+
+
+	/**
+ *
+ * Get the mesh coordinates, given a set of landmarks.
+ *
+ * @tparam vec2f
+ * @tparam vec4f
+ * @param[in] landmarks
+ * @param[in] landmark_mapper
+ * @param[in] mesh
+ * @param[in,out] model_points
+ * @param[in,out] vertex_indices
+ * @param[in,out] image_points
+ */
+	template <typename vec2f, typename vec4f>
+	inline void get_landmark_coordinates_inline(core::LandmarkCollection<vec2f> landmarks,
+										 const core::LandmarkMapper& landmark_mapper,
+										 eos::core::Mesh& mesh,
+										 vector<vector<vec4f>>& model_points,
+										 vector<vector<int>>& vertex_indices,
+										 vector<vector<vec2f>>& image_points) {
+		vector<vec4f> current_model_points;
+		vector<int> current_vertex_indices;
+		vector<vec2f> current_image_points;
+
+		for (auto &lm: landmarks) {
+			auto converted_name = landmark_mapper.convert(lm.name);
+			if (!converted_name) { // no mapping defined for the current landmark
+				continue;
 			}
 
-			model_points.push_back(current_model_points);
-			vertex_indices.push_back(current_vertex_indices);
-			image_points.push_back(current_image_points);
+			int vertex_idx = std::stoi(converted_name.get());
+			Vec4f vertex(
+				mesh.vertices[vertex_idx].x,
+				mesh.vertices[vertex_idx].y,
+				mesh.vertices[vertex_idx].z,
+				mesh.vertices[vertex_idx].w
+			);
+
+			current_model_points.emplace_back(vertex);
+			current_vertex_indices.emplace_back(vertex_idx);
+			current_image_points.emplace_back(lm.coordinates);
 		}
+
+		model_points.push_back(current_model_points);
+		vertex_indices.push_back(current_vertex_indices);
+		image_points.push_back(current_image_points);
+	}
 	}
 }
 

include/eos/fitting/closest_edge_fitting.hpp

@@ -46,12 +46,12 @@
 #include <cstddef>
 
 namespace eos {
-	namespace fitting {
+namespace fitting {
 
 /**
  * @brief Computes the intersection of the given ray with the given triangle.
  *
- * Uses the Möller-Trumbore algorithm algorithm "Fast Minimum Storage
+ * Uses the Möller-Trumbore algorithm algorithm "Fast Minimum Storage
  * Ray/
Triangle Intersection". Independent implementation, inspired by:
  * http://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/moller-trumbore-ray-triangle-intersection
  * The default eps (1e-6f) is from the paper.
@@ -391,6 +391,210 @@ inline std::pair<std::vector<cv::Vec2f>, std::vector<int>> find_occluding_edge_c
 	}
 	return { image_points, vertex_indices };
 };
+/**
+ * @brief Computes the vertices that lie on occluding boundaries, given a particular pose.
+ *
+ * This algorithm computes the edges that lie on occluding boundaries of the mesh.
+ * It performs a visibility text of each vertex, and returns a list of the (unique)
+ * vertices that make the boundary edges.
+ * An edge is defined as the line whose two adjacent faces normals flip the sign.
+ *
+ * @param[in] mesh The mesh to use.
+ * @param[in] edge_topology The edge topology of the given mesh.
+ * @param[in] R The rotation (pose) under which the occluding boundaries should be computed.
+ * @return A vector with unique vertex id's making up the edges.
+ */
+inline std::vector<int> occluding_boundary_vertices_parallel(const core::Mesh& mesh, const morphablemodel::EdgeTopology& edge_topology, glm::mat4x4 R)
+{
+	// Rotate the mesh:
+	std::vector<glm::vec4> rotated_vertices;
+	std::for_each(begin(mesh.vertices), end(mesh.vertices),
+				  [&rotated_vertices, &R](auto &&v) { rotated_vertices.push_back(R * v); });
+
+	// Compute the face normals of the rotated mesh:
+	std::vector<glm::vec3> facenormals;
+	for (auto &&f : mesh.tvi)
+	{ // for each face (triangle):
+		auto n = render::compute_face_normal(glm::vec3(rotated_vertices[f[0]]),
+											 glm::vec3(rotated_vertices[f[1]]),
+											 glm::vec3(rotated_vertices[f[2]]));
+		facenormals.push_back(n);
+	}
+
+	// Find occluding edges:
+	std::vector<int> occluding_edges_indices;
+	for (int edge_idx = 0; edge_idx < edge_topology.adjacent_faces.size(); ++edge_idx) // For each edge... Ef contains the indices of the two adjacent faces
+	{
+		const auto &edge = edge_topology.adjacent_faces[edge_idx];
+		if (edge[0] == 0) // ==> NOTE/Todo Need to change this if we use 0-based indexing!
+		{
+			// Edges with a zero index lie on the mesh boundary, i.e. they are only
+			// adjacent to one face.
+			continue;
+		}
+		// Compute the occluding edges as those where the two adjacent face normals
+		// differ in the sign of their z-component:
+		// Changing from 1-based indexing to 0-based!
+		if (glm::sign(facenormals[edge[0] - 1].z) != glm::sign(facenormals[edge[1] - 1].z))
+		{
+			// It's an occluding edge, store the index:
+			occluding_edges_indices.push_back(edge_idx);
+		}
+	}
+
+	// Select the vertices lying at the two ends of the occluding edges and remove duplicates:
+	// (This is what EdgeTopology::adjacent_vertices is needed for).
+	std::vector<int> occluding_vertices; // The model's contour vertices
+	for (auto &&edge_idx : occluding_edges_indices)
+	{
+		// Changing from 1-based indexing to 0-based!
+		occluding_vertices.push_back(edge_topology.adjacent_vertices[edge_idx][0] - 1);
+		occluding_vertices.push_back(edge_topology.adjacent_vertices[edge_idx][1] - 1);
+	}
+
+	// Remove duplicate vertex id's (std::unique works only on sorted sequences):
+	std::sort(begin(occluding_vertices), end(occluding_vertices));
+	occluding_vertices.erase(std::unique(begin(occluding_vertices), end(occluding_vertices)), end(occluding_vertices));
+
+	// Perform ray-casting to find out which vertices are not visible (i.e. self-occluded):
+	std::vector<bool> visibility;
+
+	int NUM_THREADS = 4;
+	auto t1 = std::chrono::high_resolution_clock::now();
+
+	 // we shoot the ray from the vertex towards the camera
+	glm::vec3 ray_direction(0.0f, 0.0f, 1.0f);
+
+#pragma omp parallel for schedule num_threads(NUM_THREADS)
+		for(int i = 0; i < occluding_vertices.size(); i++) {
+			const auto vertex_idx = occluding_vertices[i];
+			bool visible = true;
+
+			glm::vec3 ray_origin(rotated_vertices[vertex_idx]);
+
+			// For every tri of the rotated mesh:
+			for(int j = 0; j < mesh.tvi.size(); j++) {
+				auto tri = mesh.tvi[j];
+				auto &v0 = rotated_vertices[tri[0]];
+				auto &v1 = rotated_vertices[tri[1]];
+				auto &v2 = rotated_vertices[tri[2]];
+
+				auto intersect = ray_triangle_intersect(ray_origin, ray_direction, glm::vec3(v0), glm::vec3(v1), glm::vec3(v2), false);
+
+				// first is bool intersect, second is the distance t
+				if (intersect.first == true) {
+					// We've hit a triangle. Ray hit its own triangle. If it's behind the ray origin, ignore the intersection:
+					// Check if in front or behind?
+					if (intersect.second.get() <= 1e-4)
+					{
+						continue; // the intersection is behind the vertex, we don't care about it
+					}
+					// Otherwise, we've hit a genuine triangle, and the vertex is not visible:
+					visible = false;
+					break;
+				}
+			}
+
+			visibility.push_back(visible);
+		}
+
+	auto t2 = std::chrono::high_resolution_clock::now();
+	std::cout << "test derp 1: " << std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count() << "ms" << std::endl;
+
+	// Remove vertices from occluding boundary list that are not visible:
+	std::vector<int> final_vertex_ids;
+	for (int i = 0; i < occluding_vertices.size(); ++i)
+	{
+		if (visibility[i] == true)
+		{
+			final_vertex_ids.push_back(occluding_vertices[i]);
+		}
+	}
+	return final_vertex_ids;
+};
+
+/**
+ * @brief For a given list of 2D edge points, find corresponding 3D vertex IDs.
+ *
+ * This algorithm first computes the 3D mesh's occluding boundary vertices under
+ * the given pose. Then, for each 2D image edge point given, it searches for the
+ * closest 3D edge vertex (projected to 2D). Correspondences lying further away
+ * than \c distance_threshold (times a scale-factor) are discarded.
+ * It returns a list of the remaining image edge points and their corresponding
+ * 3D vertex ID.
+ *
+ * The given \c rendering_parameters camery_type must be CameraType::Orthographic.
+ *
+ * The units of \c distance_threshold are somewhat complicated. The function
+ * uses squared distances, and the \c distance_threshold is further multiplied
+ * with a face-size and image resolution dependent scale factor.
+ * It's reasonable to use correspondences that are 10 to 15 pixels away on a
+ * 1280x720 image with s=0.93. This would be a distance_threshold of around 200.
+ * 64 might be a conservative default.
+ *
+ * @param[in] mesh The 3D mesh.
+ * @param[in] edge_topology The mesh's edge topology (used for fast computation).
+ * @param[in] rendering_parameters Rendering (pose) parameters of the mesh.
+ * @param[in] image_edges A list of points that are edges.
+ * @param[in] distance_threshold All correspondences below this threshold.
+ * @return A pair consisting of the used image edge points and their associated 3D vertex index.
+ */
+inline std::pair<std::vector<cv::Vec2f>, std::vector<int>> find_occluding_edge_correspondences_parallel(const core::Mesh& mesh, const morphablemodel::EdgeTopology& edge_topology, const fitting::RenderingParameters& rendering_parameters, const std::vector<Eigen::Vector2f>& image_edges, float distance_threshold = 64.0f)
+{
+	assert(rendering_parameters.get_camera_type() == fitting::CameraType::Orthographic);
+	using std::vector;
+	using Eigen::Vector2f;
+
+	// Compute vertices that lye on occluding boundaries:
+	auto occluding_vertices = occluding_boundary_vertices_parallel(mesh, edge_topology, glm::mat4x4(rendering_parameters.get_rotation()));
+
+	// Project these occluding boundary vertices from 3D to 2D:
+	vector<Vector2f> model_edges_projected;
+	for (const auto& v : occluding_vertices)
+	{
+		auto p = glm::project({ mesh.vertices[v][0], mesh.vertices[v][1], mesh.vertices[v][2] }, rendering_parameters.get_modelview(), rendering_parameters.get_projection(), fitting::get_opencv_viewport(rendering_parameters.get_screen_width(), rendering_parameters.get_screen_height()));
+		model_edges_projected.push_back({ p.x, p.y });
+	}
+
+	// Find edge correspondences:
+	// Build a kd-tree and use nearest neighbour search:
+	using kd_tree_t = KDTreeVectorOfVectorsAdaptor<vector<Vector2f>, float, 2>;
+	kd_tree_t tree(2, image_edges); // dim, samples, max_leaf
+	tree.index->buildIndex();
+
+	vector<std::pair<std::size_t, double>> idx_d; // will contain [ index to the 2D edge in 'image_edges', distance (L2^2) ]
+	for (const auto& e : model_edges_projected)
+	{
+		std::size_t idx; // contains the indices into the original 'image_edges' vector
+		double dist_sq; // squared distances
+		nanoflann::KNNResultSet<double> resultSet(1);
+		resultSet.init(&idx, &dist_sq);
+		tree.index->findNeighbors(resultSet, e.data(), nanoflann::SearchParams(10));
+		idx_d.push_back({ idx, dist_sq });
+	}
+	// Filter edge matches:
+	// We filter below by discarding all correspondence that are a certain distance apart.
+	// We could also (or in addition to) discard the worst 5% of the distances or something like that.
+
+	// Filter and store the image (edge) points with their corresponding vertex id:
+	vector<int> vertex_indices;
+	vector<cv::Vec2f> image_points;
+	assert(occluding_vertices.size() == idx_d.size());
+	for (int i = 0; i < occluding_vertices.size(); ++i)
+	{
+		auto ortho_scale = rendering_parameters.get_screen_width() / rendering_parameters.get_frustum().r; // This might be a bit of a hack - we recover the "real" scaling from the SOP estimate
+		if (idx_d[i].second <= distance_threshold * ortho_scale) // I think multiplying by the scale is good here and gives us invariance w.r.t. the image resolution and face size.
+		{
+			auto edge_point = image_edges[idx_d[i].first];
+			// Store the found 2D edge point, and the associated vertex id:
+			vertex_indices.push_back(occluding_vertices[i]);
+			image_points.push_back(cv::Vec2f(edge_point[0], edge_point[1]));
+		}
+	}
+	return { image_points, vertex_indices };
+};
+
+
 
 	} /* namespace fitting */
 } /* namespace eos */

include/eos/fitting/fitting.hpp

@@ -43,6 +43,8 @@
 #include <vector>
 #include <algorithm>
 
+#include <omp.h>
+
 namespace eos {
 	namespace fitting {
 
@@ -903,7 +905,7 @@ inline eos::core::Mesh generate_new_mesh(
 }
 
 
-inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParameters>> fit_shape_and_pose_multi_2(
+inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParameters>> fit_shape_and_pose_multi_parallel(
 		const morphablemodel::MorphableModel& morphable_model,
 		const std::vector<morphablemodel::Blendshape>& blendshapes,
 		std::vector<eos::video::Keyframe>& key_frames,
@@ -953,6 +955,275 @@ inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParamete
 
 	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);
+	vector<VectorXf> current_combined_shapes(num_images);
+	vector<eos::core::Mesh> current_meshs(num_images);
+
+	// The 2D and 3D point correspondences used for the fitting:
+	vector<vector<Vec4f>> model_points(num_images); // the points in the 3D shape model of all frames.
+	vector<vector<int>> vertex_indices(num_images); // their vertex indices of all frames.
+	vector<vector<Vec2f>> image_points(num_images); // the corresponding 2D landmark points of all frames.
+
+	int NUM_THREADS = 4;
+
+	vector<fitting::RenderingParameters> rendering_params(num_images); // list of rendering params for all frames.
+
+#pragma omp parallel num_threads(NUM_THREADS)
+	{
+#pragma omp for
+	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());
+
+		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());
+
+		vector<Vec4f> curr_model_points;
+		vector<int> curr_vertex_indices;
+		vector<Vec2f> curr_image_points;
+
+		// Get the locations of the model locations of the meshes, vertex_indices and image points
+		// (equal to landmark coordinates), for every image / mesh.
+		std::tie(curr_model_points, curr_vertex_indices, curr_image_points) = eos::core::get_landmark_coordinates<Vec2f, Vec4f>(
+			key_frames[j].fitting_result.landmarks, landmark_mapper, current_mesh);
+
+		// Start constructing a list of rendering parameters needed for reconstruction.
+		// Get the current points from the last added image points and model points
+		auto current_pose = fitting::estimate_orthographic_projection_linear(
+			curr_image_points, curr_model_points, true, image_height);
+
+		fitting::RenderingParameters current_rendering_params(current_pose, image_width, image_height);
+		rendering_params[j] = current_rendering_params;
+
+		// update key frame rendering params
+		key_frames[j].fitting_result.rendering_parameters = current_rendering_params;
+
+		Mat affine_from_ortho = fitting::get_3x4_affine_camera_matrix(current_rendering_params, image_width, image_height);
+
+		blendshape_coefficients[j] = fitting::fit_blendshapes_to_landmarks_nnls(
+			blendshapes, current_pca_shape, affine_from_ortho, curr_image_points, curr_vertex_indices);
+
+		// 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[j].data(), blendshape_coefficients[j].size()
+		);
+
+		current_combined_shapes[j] = current_combined_shape;
+
+		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[j] = current_mesh;
+		model_points[j] = curr_model_points;
+		vertex_indices[j] = curr_vertex_indices;
+		image_points[j] = curr_image_points;
+	}
+}
+	// The static (fixed) landmark correspondences which will stay the same throughout
+	// the fitting (the inner face landmarks):
+	vector<vector<int>> fixed_vertex_indices (vertex_indices);
+	vector<vector<Vec2f>> fixed_image_points (image_points);
+
+	std::vector<cv::Mat> affine_from_orthos(num_images);
+	std::vector<VectorXf> mean_plus_blendshapes(num_images);
+
+	image_points = fixed_image_points;
+	vertex_indices = fixed_vertex_indices;
+
+#pragma omp parallel num_threads(NUM_THREADS)
+	{
+#pragma omp for
+		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 curr_keyframe = key_frames[j];
+			auto landmarks = curr_keyframe.fitting_result.landmarks;
+			float yaw_angle = curr_keyframe.yaw_angle;
+
+//			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_meshs[j],
+				rendering_params[j].get_modelview(),
+				rendering_params[j].get_projection(),
+				fitting::get_opencv_viewport(image_width, image_height)
+			);
+
+			// Add the contour correspondences to the set of landmarks that we use for the fitting:
+			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;
+
+			// positive yaw = subject looking to the left
+			if (yaw_angle >= 0.0f) {
+				// 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
+				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);
+				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_parallel(
+				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[j].clear();
+
+			for (const auto &v : vertex_indices[j])
+			{
+				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:
+			auto current_pose = fitting::estimate_orthographic_projection_linear(image_points[j], model_points[j], true, image_height);
+			rendering_params[j] = fitting::RenderingParameters(current_pose, image_width, image_height);
+
+			Mat affine_from_ortho = fitting::get_3x4_affine_camera_matrix(rendering_params[j], image_width, image_height);
+
+			affine_from_orthos[j] = affine_from_ortho;
+
+			// Estimate the PCA shape coefficients with the current blendshape coefficients:
+			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[j] = current_mean_plus_blendshapes;
+		}
+	}
+
+	pca_shape_coefficients = fitting::fit_shape_to_landmarks_linear_multi_parallel(
+		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);
+
+#pragma omp parallel num_threads(NUM_THREADS)
+	{
+#pragma omp for
+		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_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).
+};
+
+inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParameters>> fit_shape_and_pose_multi_2(
+	const morphablemodel::MorphableModel& morphable_model,
+	const std::vector<morphablemodel::Blendshape>& blendshapes,
+	std::vector<eos::video::Keyframe>& key_frames,
+	const core::LandmarkMapper& landmark_mapper,
+	int image_width,
+	int image_height,
+	int num_images,
+	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(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());
+
+	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();
+	}
+
+	if (pca_shape_coefficients.empty()) {
+		pca_shape_coefficients.resize(num_shape_coefficients_to_fit.get());
+	}
+
+	// TODO: This leaves the following case open: num_coeffs given is empty or defined, but the
+	// pca_shape_coefficients given is != num_coeffs or the model's max-coeffs. What to do then? Handle & document!
+	if (blendshape_coefficients.size() < num_images) {
+		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);
 	vector<VectorXf> current_combined_shapes;
@@ -969,7 +1240,7 @@ inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParamete
 	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());
+				Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients[j].data(), blendshape_coefficients[j].size());
 
 		eos::core::Mesh current_mesh = morphablemodel::sample_to_mesh(
 			current_combined_shape,
@@ -980,7 +1251,10 @@ inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParamete
 
 		// Get the locations of the model locations of the meshes, vertex_indices and image points
 		// (equal to landmark coordinates), for every image / mesh.
-		eos::core::get_landmark_coordinates<Vec2f, Vec4f>(
+
+//		std::tie(curr_model_points, curr_vertex_indices, curr_image_points) = eos::core::get_landmark_coordinates<Vec2f, Vec4f>(
+//			key_frames[j].fitting_result.landmarks, landmark_mapper, current_mesh);
+		eos::core::get_landmark_coordinates_inline<Vec2f, Vec4f>(
 			key_frames[j].fitting_result.landmarks,
 			landmark_mapper,
 			current_mesh,
@@ -1010,7 +1284,7 @@ inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParamete
 		current_combined_shape = current_pca_shape +
 			morphablemodel::to_matrix(blendshapes) *
 				Eigen::Map<const Eigen::VectorXf>(blendshape_coefficients[j].data(), blendshape_coefficients[j].size()
-		);
+				);
 
 		current_combined_shapes.push_back(current_combined_shape);
 
@@ -1065,7 +1339,7 @@ inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParamete
 
 		// positive yaw = subject looking to the left
 		if (yaw_angle >= 0.0f) {
-			 // the left contour is the occluding one we want to use ("away-facing")
+			// 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
 			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] });
@@ -1109,7 +1383,7 @@ inline std::pair<std::vector<core::Mesh>, std::vector<fitting::RenderingParamete
 		mean_plus_blendshapes.push_back(current_mean_plus_blendshapes);
 	}
 
-	pca_shape_coefficients = fitting::fit_shape_to_landmarks_linear_multi(
+	pca_shape_coefficients = fitting::fit_shape_to_landmarks_linear_multi_parallel(
 		morphable_model,
 		affine_from_orthos,
 		image_points,

include/eos/fitting/linear_shape_fitting.hpp

@@ -34,6 +34,10 @@
 #include <vector>
 #include <cassert>
 
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
 namespace eos {
 	namespace fitting {
 
@@ -110,9 +114,7 @@ inline std::vector<float> fit_shape_to_landmarks_linear_multi(
 	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 (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
 
@@ -127,6 +129,7 @@ inline std::vector<float> fit_shape_to_landmarks_linear_multi(
         // 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) {
+			// THIS ONE IS EXPENSE 12ms, but more the inner function generic_dense.. is expense
             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
@@ -206,6 +209,200 @@ inline std::vector<float> fit_shape_to_landmarks_linear(const morphablemodel::Mo
     return fit_shape_to_landmarks_linear_multi(morphable_model, { affine_camera_matrix }, all_landmarks, all_vertex_ids, { base_face }, lambda, num_coefficients_to_fit, detector_standard_deviation, model_standard_deviation );
 }
 
+/**
+ * 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). 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).
+ */
+#ifdef _OPENMP
+inline std::vector<float> fit_shape_to_landmarks_linear_multi_parallel(
+	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_images = affine_camera_matrix.size();
+
+	// the regularisation has to be adjusted when more than one image is given
+	lambda *= num_images;
+
+	int num_threads = 4;
+	int total_num_landmarks_dimension = 0;
+
+	std::vector<int> landmark_offset_index(landmarks.size());
+
+	for (int i = 0; i < landmarks.size(); i++)
+	{
+		landmark_offset_index[i] = total_num_landmarks_dimension;
+		total_num_landmarks_dimension += landmarks[i].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 * 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:
+	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
+
+	// The variances: Add the 2D and 3D standard deviations.
+	// If the user doesn't provide them, we choose the following:
+	// 2D (detector) standard deviation: In pixel, we follow [1] and choose sqrt(3) as the default value.
+	// 3D (model) variance: 0.0f. It only makes sense to set it to something when we have a different variance for different vertices.
+	// The 3D variance has to be projected to 2D (for details, see paper [1]) so the units do match up.
+	float sigma_squared_2D = std::pow(detector_standard_deviation.get_value_or(std::sqrt(3.0f)), 2)
+		+ std::pow(model_standard_deviation.get_value_or(0.0f), 2);
+	// We use a VectorXf, and later use .asDiagonal():
+	VectorXf Omega = VectorXf::Constant(3 * total_num_landmarks_dimension, 1.0f / sigma_squared_2D);
+	// The landmarks in matrix notation (in homogeneous coordinates), $3N\times 1$
+	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.
+
+	auto t1 = std::chrono::high_resolution_clock::now();
+	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());
+		int affine_camera_rows = affine_camera_matrix[k].rows;
+		int affine_camera_cols = affine_camera_matrix[k].cols;
+
+		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);
+
+#pragma omp parallel num_threads(num_threads)
+		{
+#pragma omp for
+			for (int i = 0; i < num_landmarks; ++i)
+			{
+				int index = (landmark_offset_index[k] * 4) + 4 * i;
+				int offset = landmark_offset_index[k] + i;
+				auto vertex_id = vertex_ids[k][i];
+
+				// THIS ONE IS EXPENSE 12ms, but more the inner function generic_dense.. is expense
+				MatrixXf basis_rows_ = morphable_model.get_shape_model()
+					.get_rescaled_pca_basis_at_point(vertex_id); // 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(index, 0, 3, V_hat_h.cols()) =
+					basis_rows_.block(0, 0, basis_rows_.rows(), num_coeffs_to_fit);
+
+				// Set y with landmark coordinates (will be used in linear lsq later)
+#pragma omp simd
+				for (int j = 0; j < 2; j++)
+				{
+					y(3 * offset + j) = landmarks[k][i][j];
+				}
+
+				// Set v_bar with face mesh (will be used in linear lsq later)
+#pragma omp simd
+				for (int j = 0; j < 3; j++)
+				{
+					v_bar(4 * offset + j) = base_face[k](vertex_id * 3 + j);
+				}
+			}
+
+			/* Note, cannot be executed in parallel yet due to push_back. Will need to find a way to assign instead of
+			 * using push_back, order matters.. : */
+#pragma omp single
+			{
+				for (int i = 0; i < num_landmarks; ++i)
+				{
+					int p_index = landmark_offset_index[k] + i;
+					for (int x = 0; x < affine_camera_rows; ++x)
+					{
+						for (int y = 0; y < affine_camera_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)));
+						}
+					}
+				}
+			}
+
+			// 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:
+		}
+	} // forloop
+
+	auto t2 = std::chrono::high_resolution_clock::now();
+
+	t1 = std::chrono::high_resolution_clock::now();
+	// fill P with coefficients
+	P.setFromTriplets(P_coefficients.begin(), P_coefficients.end());
+	MatrixXf one;
+	MatrixXf two;
+	MatrixXf A;
+	MatrixXf b;
+
+#pragma omp parallel sections num_threads(2)
+{
+#pragma omp section
+{
+	A = P * V_hat_h; // camera matrix times the basis
+	one = A.transpose() * Omega.asDiagonal();
+}
+
+#pragma omp section
+{
+	b = P * v_bar - y; // camera matrix times the mean, minus the landmarks
+	two = lambda * Eigen::MatrixXf::Identity(num_coeffs_to_fit, num_coeffs_to_fit);
+}
+}
+
+	const MatrixXf AtOmegaAReg = one * A + two;
+	const MatrixXf rhs = -one * b; // It's -A^t*Omega^t*b, but we don't need to transpose Omega, since it's a diagonal matrix, and Omega^t = Omega.
+	const VectorXf c_s = AtOmegaAReg.colPivHouseholderQr().solve(rhs);
+
+	t2 = std::chrono::high_resolution_clock::now();
+	std::cout << "matrix xf " << std::chrono::duration_cast<std::chrono::milliseconds>(t2-t1).count() << "ms" << std::endl;
+
+	return std::vector<float>(c_s.data(), c_s.data() + c_s.size());
+};
+#endif
+
+
 	} /* namespace fitting */
 } /* namespace eos */
 

include/eos/video/BufferedVideoIterator.hpp

@@ -93,40 +93,44 @@ public:
 	// TODO: build support for setting the amount of max_frames in the buffer.
 	BufferedVideoIterator(std::string videoFilePath, fitting::ReconstructionData reconstruction_data, boost::property_tree::ptree settings) {
 		std::ifstream file(videoFilePath);
-		std::cout << "video file path: " << videoFilePath << std::endl;
+		std::cout << "Opening video: " << videoFilePath << std::endl;
 
 		if (!file.is_open()) {
 			throw std::runtime_error("Error opening given file: " + videoFilePath);
 		}
 
 		cv::VideoCapture tmp_cap(videoFilePath); // open video file
-		if (!tmp_cap.isOpened()) { // check if we succeeded
+		if (!tmp_cap.isOpened()) {
 			throw std::runtime_error("Could not play video");
 		}
 
-		this->cap = tmp_cap;
+		cap = tmp_cap;
+
+		this->reconstruction_data = reconstruction_data;
 
 		// copy settings from gathered from a .ini file
-		this->min_frames = settings.get<int>("video.min_frames", 5);
-		this->drop_frames = settings.get<int>("video.drop_frames", 0);
-		this->skip_frames = settings.get<int>("video.skip_frames", 0);
-		this->frames_per_bin = settings.get<unsigned int>("video.frames_per_bin", 2);
-		this->num_shape_coefficients_to_fit = settings.get<unsigned int>("video.num_shape_coefficients_to_fit", 50);
+		min_frames = settings.get<int>("video.min_frames", 5);
+		drop_frames = settings.get<int>("video.drop_frames", 0);
+		skip_frames = settings.get<int>("video.skip_frames", 0);
+		frames_per_bin = settings.get<unsigned int>("video.frames_per_bin", 2);
 
-		this->reconstruction_data = reconstruction_data;
+		unsigned int num_shape_coeff = reconstruction_data.morphable_model.get_shape_model().get_num_principal_components();
+
+		this->num_shape_coefficients_to_fit = settings.get<unsigned int>(
+			"video.num_shape_coefficients_to_fit", num_shape_coeff);
 
 		// initialize bins
-		this->bins.resize(num_yaw_bins);
+		bins.resize(num_yaw_bins);
 
-		// reset all
-		this->n_frames = 0;
-		this->total_frames = 0;
+		// reset frame count
+		n_frames = 0;
+		total_frames = 0;
 
-		std::cout << "Buffered video iter: "
+		std::cout << "Settings: " << std::endl <<
 			"min_frames: " << min_frames << std::endl <<
-				  "drop_frames: " << drop_frames << std::endl <<
-				  "frames_per_bin: " << frames_per_bin << std::endl <<
-				  "num_shape_coefficients_to_fit: " <<  num_shape_coefficients_to_fit << std::endl;
+		    "drop_frames: " << drop_frames << std::endl <<
+		    "frames_per_bin: " << frames_per_bin << std::endl <<
+		    "num_shape_coefficients_to_fit: " << num_shape_coefficients_to_fit << std::endl;
 
 		std::cout << "total frames in video: " << cap.get(CV_CAP_PROP_FRAME_COUNT) << std::endl;
 	}
@@ -136,9 +140,11 @@ public:
 	}
 
 	/**
+	 * Generate a new keyframe containing information about pose and landmarks
+	 * These are needed to determine if we want the image in the first place.
 	 *
 	 * @param frame
-	 * @return
+	 * @return Keyframe
 	 */
 	Keyframe generate_new_keyframe(cv::Mat frame) {
 		int frame_height = frame.rows;
@@ -153,6 +159,7 @@ public:
 			return Keyframe();
 		}
 
+		// Get the necessary information for reconstruction.
 		auto landmarks = reconstruction_data.landmark_list[total_frames];
 		auto landmark_mapper = reconstruction_data.landmark_mapper;
 		auto blendshapes = reconstruction_data.blendshapes;
@@ -162,10 +169,6 @@ public:
 		vector<int> vertex_indices;
 		vector<cv::Vec2f> image_points;
 
-		if (num_shape_coefficients_to_fit == 0) {
-			num_shape_coefficients_to_fit = morphable_model.get_shape_model().get_num_principal_components();
-		}
-
 		// make a new one
 		std::vector<float> blend_shape_coefficients;
 
@@ -193,16 +196,6 @@ public:
 		fitting_result.rendering_parameters = rendering_params;
 		fitting_result.landmarks = landmarks;
 
-//		std::cout << face_roi << " ("<< frame_width << "," << frame_height << ")" << std::endl;
-//		for (auto &&lm : landmarks) {
-//			cv::rectangle(
-//					frame, cv::Point2f(lm.coordinates[0] - 2.0f, lm.coordinates[1] - 2.0f),
-//					cv::Point2f(lm.coordinates[0], lm.coordinates[1] + 2.0f), {255, 0, 0}
-//			);
-//		}
-//
-//		cv::imshow("frame", frame);
-//		cv::waitKey(0);
 		cv::Rect face_roi = core::get_face_roi(image_points, frame_width, frame_height);
 
 		float frame_laplacian_score = static_cast<float>(variance_of_laplacian(frame(face_roi)));
@@ -359,7 +352,7 @@ public:
 
 			// Setting that the buffer has changed:
 			frame_buffer_changed = true;
-			std::cout << "frame added(" << n_frames << "): " << keyframe.score << ", " << keyframe.yaw_angle << std::endl;
+			std::cout << "frame added(" << total_frames << "): " << keyframe.score << ", " << keyframe.yaw_angle << std::endl;
 		}
 
 		total_frames++;
@@ -384,8 +377,14 @@ public:
 		this->pca_shape_coefficients = pca_shape_coefficients;
 	}
 
-	// Converts a given yaw angle to an index in the internal bins vector.
-	// Assumes 9 bins and 20� intervals.
+	/**
+	 *
+	 * Converts a given yaw angle to an index in the internal bins vector.
+	 * Assumes 9 bins and 20� intervals.
+	 *
+	 * @param yaw_angle
+	 * @return
+	 */
 	static std::size_t angle_to_index(float yaw_angle)
 	{
 		if (yaw_angle <= -70.f)
@@ -436,15 +435,14 @@ public:
 
 
 	// todo: move to private if possible.
+	/**
+	 * Stop by releasing the VideoCapture.
+	 */
 	void __stop() {
 		cap.release();
 	};
 
 private:
-	int num_yaw_bins = 9;
-	unsigned int frames_per_bin;
-	bool frame_buffer_changed = false;
-
 	cv::VideoCapture cap;
 	std::deque<Keyframe> frame_buffer;
 	eos::fitting::ReconstructionData reconstruction_data;
@@ -455,6 +453,10 @@ private:
 	// latest pca_shape_coefficients
 	std::vector<float> pca_shape_coefficients;
 
+	std::size_t num_yaw_bins = 9;
+	bool frame_buffer_changed = false;
+	unsigned int frames_per_bin;
+
 	// TODO: make set-able
 	// total frames in processed, not persee in buffer (skipped frames)
 	int total_frames;
@@ -472,9 +474,6 @@ private:
 	int drop_frames = 0;
 
 	unsigned int num_shape_coefficients_to_fit = 0;
-
-	// laplacian threshold
-//	double laplacian_threshold = 10000000;
 };
 
 }

include/eos/video/Keyframe.hpp

@@ -227,6 +227,33 @@ public:
 		threshold = static_cast<unsigned char>(alpha_thresh);
 	};
 
+	/**
+	 *
+	 * @param frame
+	 * @param landmarks
+	 * @return
+	 */
+	bool has_eyes_open(cv::Mat frame, core::LandmarkCollection<Vec2f> landmarks) {
+		std::vector<cv::Vec2f> right_eye;
+		std::vector<cv::Vec2f> left_eye;
+
+		for(auto &&lm: landmarks) {
+			int landmark_index = std::stoi(lm.name);
+			if(landmark_index >= 37 && landmark_index <= 42) {
+				right_eye.push_back(lm.coordinates);
+			}
+
+			if(landmark_index >= 43 && landmark_index <= 48) {
+				left_eye.push_back(lm.coordinates);
+			}
+		}
+
+		cv::Rect right_box = cv::boundingRect(right_eye);
+		cv::Rect left_box = cv::boundingRect(left_eye);
+
+		return true;
+	}
+
 	/**
 	 * @brief Merges the given new isomap with all previously processed isomaps.
 	 *

utils/accuracy-evaluation.cpp

@@ -18,8 +18,8 @@
  */
 #include "eos/core/landmark_utils.hpp"
 #include "eos/render/render.hpp"
-#include "eos/morphablemodel/morphablemodel.hpp"
-#include "eos/morphablemodel/blendshape.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"
@@ -213,9 +213,6 @@ void render_output(
 		cv::imwrite("/tmp/eos/" +
 			std::to_string(frame_number) + "_" + std::to_string(yaw_angle) + ".png", frame);
 
-		int frame_width = frame.cols;
-		int frame_height = frame.rows;
-
 		// Extract the texture using the fitted mesh from this frame:
 		Mat affine_cam = fitting::get_3x4_affine_camera_matrix(rendering_paramss[i], frame.cols, frame.rows);
 		Mat isomap = render::extract_texture(
@@ -366,22 +363,28 @@ void evaluate_results(
 		Mat affine_from_ortho = fitting::get_3x4_affine_camera_matrix(
 			rendering_params, frame_width, frame_height);
 
-//		for (auto &&lm : landmarks) {
-//			cv::rectangle(
-//					outimg, cv::Point2f(lm.coordinates[0] - 2.0f, lm.coordinates[1] - 2.0f),
-//					cv::Point2f(lm.coordinates[0], lm.coordinates[1] + 2.0f), {255, 0, 0}
-//			);
-//		}
-//		// Draw the fitted mesh as wireframe, and save the image:
-//		draw_wireframe(
-//				outimg,
-//				meshs[i],
-//				rendering_params.get_modelview(),
-//				rendering_params.get_projection(),
-//				fitting::get_opencv_viewport(frame_width, frame_height));
-//
+		for (auto &&lm : landmarks) {
+			cv::rectangle(
+					outimg, cv::Point2f(lm.coordinates[0] - 2.0f, lm.coordinates[1] - 2.0f),
+					cv::Point2f(lm.coordinates[0], lm.coordinates[1] + 2.0f), {255, 0, 0}
+			);
+		}
+		// Draw the fitted mesh as wireframe, and save the image:
+		draw_wireframe(
+				outimg,
+				meshs[i],
+				rendering_params.get_modelview(),
+				rendering_params.get_projection(),
+				fitting::get_opencv_viewport(frame_width, frame_height));
+
+		bool eyes_open = isomap_averaging.has_eyes_open(frame, landmarks);
+
+		if (!eyes_open) {
+			continue;
+		}
+
 //		cv::imshow("Img", outimg);
-//		cv::waitKey(1);
+//		cv::waitKey(0);
 
 		// Draw the loaded landmarks:
 		Mat isomap = render::extract_texture(meshs[i], affine_from_ortho, frame);
@@ -660,7 +663,8 @@ int main(int argc, char *argv[]) {
 		if (vid_iterator.has_changed()) {
 			std::cout << "Going to reconstruct with " << key_frames.size() << " images."<< std::endl;
 			// Fit shape and pose:
-			std::tie(meshs, rendering_paramss) = fitting::fit_shape_and_pose_multi_2(
+			auto t1 = std::chrono::high_resolution_clock::now();
+			std::tie(meshs, rendering_paramss) = fitting::fit_shape_and_pose_multi_parallel(
 				morphable_model,
 				blendshapes,
 				key_frames,
@@ -680,6 +684,12 @@ int main(int argc, char *argv[]) {
 				fitted_image_points
 			);
 
+			auto t2 = std::chrono::high_resolution_clock::now();
+			std::cout << "Reconstruction took "
+					  << std::chrono::duration_cast<std::chrono::milliseconds>(t2-t1).count()
+					  << "ms, mean(" << std::chrono::duration_cast<std::chrono::milliseconds>(t2-t1).count() / key_frames.size()
+					  << "ms)" << std::endl;
+
 			vid_iterator.update_reconstruction_coeff(pca_shape_coefficients);
 
 			evaluate_results(