Merge pull request #6 from patrikhuber/devel

v0.6.0: OpenGL-like rendering, texture rendering, texture extraction view-angle computation

CMakeLists.txt

@@ -25,6 +25,7 @@ if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
 			set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14")
 		endif()
 	endif()
+	# Note: gcc is fine without -pthreads.
 elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") # the quotes are needed here, maybe because "MSVC" seems to be a keyword
 	if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 19)
 		message(FATAL_ERROR "Visual Studio 2015 or newer is required.")
@@ -38,7 +39,10 @@ elseif(CMAKE_CXX_COMPILER_ID STREQUAL "Clang")
 	if(HAS_CXX14_FLAG)
 		set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14")
 	endif()
-else()
+	set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -pthreads")
+	# Eigen::LevenbergMarquardt probably needs -pthreads.
+	# Cleaner way would be to add it to fit-model's target_link_libraries, but that requires a CMake >2.8.10.
+else() # no GNU, no MSVC, no Clang
 	message(WARNING "You are using an unsupported compiler. Compilation has only been tested with MSVC, GCC and Clang.")
 	check_cxx_compiler_flag(-std=c++14 HAS_CXX14_FLAG)
 	if(HAS_CXX14_FLAG)
@@ -96,8 +100,10 @@ set(HEADERS
 	include/eos/fitting/linear_shape_fitting.hpp
 	include/eos/render/Mesh.hpp
 	include/eos/render/utils.hpp
+	include/eos/render/render.hpp
 	include/eos/render/render_affine.hpp
 	include/eos/render/detail/render_detail.hpp
+	include/eos/render/detail/render_affine_detail.hpp
 	include/eos/render/texture_extraction.hpp
 	include/eos/render/detail/texture_extraction_detail.hpp
 )

examples/fit-model.cpp

@@ -19,7 +19,7 @@
  */
 #include "eos/core/Landmark.hpp"
 #include "eos/core/LandmarkMapper.hpp"
-#include "eos/fitting/affine_camera_estimation.hpp"
+#include "eos/fitting/nonlinear_camera_estimation.hpp"
 #include "eos/fitting/linear_shape_fitting.hpp"
 #include "eos/render/utils.hpp"
 #include "eos/render/texture_extraction.hpp"
@@ -120,11 +120,11 @@ int main(int argc, char *argv[])
 		desc.add_options()
 			("help,h",
 				"display the help message")
-			("model,m", po::value<fs::path>(&modelfile)->required(),
+			("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<fs::path>(&imagefile)->required()->default_value("data/image_0001.png"),
+			("image,i", po::value<fs::path>(&imagefile)->required()->default_value("data/image_0010.png"),
 				"an input image")
-			("landmarks,l", po::value<fs::path>(&landmarksfile)->required()->default_value("data/image_0001.pts"),
+			("landmarks,l", po::value<fs::path>(&landmarksfile)->required()->default_value("data/image_0010.pts"),
 				"2D landmarks for the image, in ibug .pts format")
 			("mapping,p", po::value<fs::path>(&mappingsfile)->required()->default_value("../share/ibug2did.txt"),
 				"landmark identifier to model vertex number mapping")
@@ -178,7 +178,6 @@ int main(int argc, char *argv[])
 	vector<Vec2f> image_points; // the corresponding 2D landmark points
 
 	// Sub-select all the landmarks which we have a mapping for (i.e. that are defined in the 3DMM):
-	//std::transform(begin(landmarks), end(landmarks), begin(landmarks), [&landmark_mapper](const Landmark<Vec2f>& lm) {  });
 	for (int i = 0; i < landmarks.size(); ++i) {
 		auto converted_name = landmark_mapper.convert(landmarks[i].name);
 		if (!converted_name) { // no mapping defined for the current landmark
@@ -191,42 +190,32 @@ int main(int argc, char *argv[])
 		image_points.emplace_back(landmarks[i].coordinates);
 	}
 
-	// Estimate the camera from the 2D - 3D point correspondences
-	Mat affine_cam = fitting::estimate_affine_camera(image_points, model_points);
+	// Estimate the camera (pose) from the 2D - 3D point correspondences
+	fitting::OrthographicRenderingParameters rendering_params = fitting::estimate_orthographic_camera(image_points, model_points, image.cols, image.rows);
+	Mat affine_from_ortho = get_3x4_affine_camera_matrix(rendering_params, image.cols, image.rows);
 
-	// Draw the mean-face landmarks projected using the estimated camera:
-	for (auto&& vertex : model_points) {
-		Vec2f screen_point(Mat(affine_cam * Mat(vertex)).at<float>(0), Mat(affine_cam * Mat(vertex)).at<float>(1));
-		cv::circle(outimg, cv::Point2f(screen_point), 5, { 0.0f, 255.0f, 0.0f });
-	}
+	// The 3D head pose can be recovered as follows:
+	float yaw_angle = glm::degrees(rendering_params.r_y);
+	// and similarly for pitch (r_x) and roll (r_z).
 
 	// Estimate the shape coefficients by fitting the shape to the landmarks:
-	vector<float> fitted_coeffs = fitting::fit_shape_to_landmarks_linear(morphable_model, affine_cam, image_points, vertex_indices);
+	vector<float> fitted_coeffs = fitting::fit_shape_to_landmarks_linear(morphable_model, affine_from_ortho, image_points, vertex_indices);
 
-	// Obtain the full mesh and draw it using the estimated camera:
+	// Obtain the full mesh with the estimated coefficients:
 	render::Mesh mesh = morphable_model.draw_sample(fitted_coeffs, vector<float>());
-	outputfile += fs::path(".obj");
-	render::write_textured_obj(mesh, outputfile.string()); // save the mesh as obj
 
-	// Draw the projected points again, this time using the fitted model shape:
-	for (auto&& idx : vertex_indices) {
-		Vec4f model_point(mesh.vertices[idx][0], mesh.vertices[idx][1], mesh.vertices[idx][2], mesh.vertices[idx][3]);
-		Vec2f screen_point(Mat(affine_cam * Mat(model_point)).at<float>(0), Mat(affine_cam * Mat(model_point)).at<float>(1));
-		cv::circle(outimg, cv::Point2f(screen_point), 3, { 0.0f, 0.0f, 255.0f });
-	}
+	// Extract the texture from the image using given mesh and camera parameters:
+	Mat isomap = render::extract_texture(mesh, affine_from_ortho, image);
 
-	// Save an output image with the landmarks from the different stages:
-	//outputfile.replace_extension(".png");
-	//cv::imwrite(outputfile.string(), outimg);
-	outputfile.replace_extension(".png");
-	cv::imwrite(outputfile.string(), outimg);
+	// Save the mesh as textured obj:
+	outputfile += fs::path(".obj");
+	render::write_textured_obj(mesh, outputfile.string());
 
-	// Extract the texture and save the extracted texture map (isomap):
-	Mat isomap = render::extract_texture(mesh, affine_cam, image, render::TextureInterpolation::NearestNeighbour);
+	// And save the isomap:
 	outputfile.replace_extension(".isomap.png");
 	cv::imwrite(outputfile.string(), isomap);
 
-	cout << "Finished fitting and wrote result image and isomap " << outputfile.string() << "." << endl;
+	cout << "Finished fitting and wrote result mesh and isomap to files with basename " << outputfile.stem().stem() << "." << endl;
 
 	return EXIT_SUCCESS;
 }

include/eos/fitting/linear_shape_fitting.hpp

@@ -48,7 +48,7 @@ namespace eos {
  * Note: The standard deviations given should be a vector, i.e. different for each landmark. This is not implemented yet.
  *
  * @param[in] morphable_model The Morphable Model whose shape (coefficients) are estimated.
- * @param[in] affine_camera_matrix A 3x4 affine camera matrix from world to clip-space (should probably be of type CV_32FC1 as all our calculations are done with float).
+ * @param[in] 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, given in clip-coordinates.
  * @param[in] vertex_ids The vertex ids in the model that correspond to the 2D points.
  * @param[in] lambda The regularisation parameter (weight of the prior towards the mean).

include/eos/fitting/nonlinear_camera_estimation.hpp

@@ -64,7 +64,7 @@ struct Frustum
  * The rotation values are given in radians and estimated using the RPY convention.
  * Yaw is applied first to the model, then pitch, then roll (R * P * Y * vertex).
  */
-struct RenderingParameters
+struct OrthographicRenderingParameters
 {
 	float r_x; // Pitch.
 	float r_y; // Yaw. Positive means subject is looking left (we see her right cheek).
@@ -74,6 +74,92 @@ struct RenderingParameters
 	Frustum frustum;
 };
 
+/**
+ * @brief Converts a glm::mat4x4 to a cv::Mat.
+ *
+ * Note: move to render namespace
+ */
+cv::Mat to_mat(const glm::mat4x4& glm_matrix)
+{
+	// glm stores its matrices in col-major order in memory, OpenCV in row-major order.
+	// Hence we transpose the glm matrix to flip the memory layout, and then point OpenCV
+	// to that location.
+	auto glm_matrix_t = glm::transpose(glm_matrix);
+	cv::Mat opencv_mat(4, 4, CV_32FC1, &glm_matrix_t[0]);
+	// we need to clone because the underlying data of the original goes out of scope
+	return opencv_mat.clone();
+};
+
+/**
+ * @brief Creates a 4x4 model-view matrix from given fitting parameters.
+ *
+ * Together with the Frustum information, this describes the full
+ * orthographic rendering parameters of the OpenGL pipeline.
+ * Example:
+ *
+ * @code
+ * fitting::OrthographicRenderingParameters rendering_params = ...;
+ * glm::mat4x4 view_model = get_4x4_modelview_matrix(rendering_params);
+ * glm::mat4x4 ortho_projection = glm::ortho(rendering_params.frustum.l, rendering_params.frustum.r, rendering_params.frustum.b, rendering_params.frustum.t);
+ * glm::vec4 viewport(0, image.rows, image.cols, -image.rows); // flips y, origin top-left, like in OpenCV
+ *
+ * // project a point from 3D to 2D:
+ * glm::vec3 point_3d = ...; // from a mesh for example
+ * glm::vec3 point_2d = glm::project(point_3d, view_model, ortho_projection, viewport);
+ * @endcode
+ */
+glm::mat4x4 get_4x4_modelview_matrix(fitting::OrthographicRenderingParameters params)
+{
+	// rotation order: RPY * P
+	auto rot_mtx_x = glm::rotate(glm::mat4(1.0f), params.r_x, glm::vec3{ 1.0f, 0.0f, 0.0f });
+	auto rot_mtx_y = glm::rotate(glm::mat4(1.0f), params.r_y, glm::vec3{ 0.0f, 1.0f, 0.0f });
+	auto rot_mtx_z = glm::rotate(glm::mat4(1.0f), params.r_z, glm::vec3{ 0.0f, 0.0f, 1.0f });
+	auto t_mtx = glm::translate(glm::mat4(1.0f), glm::vec3{ params.t_x, params.t_y, 0.0f });
+	auto modelview = t_mtx * rot_mtx_z * rot_mtx_x * rot_mtx_y;
+	return modelview;
+};
+
+/**
+ * @brief Creates a 3x4 affine camera matrix from given fitting parameters. The
+ * matrix transforms points directly from model-space to screen-space.
+ *
+ * This function is mainly used since the linear shape fitting fitting::fit_shape_to_landmarks_linear
+ * expects one of these 3x4 affine camera matrices, as well as render::extract_texture.
+ */
+cv::Mat get_3x4_affine_camera_matrix(fitting::OrthographicRenderingParameters params, int width, int height)
+{
+	auto view_model = to_mat(get_4x4_modelview_matrix(params));
+	auto ortho_projection = to_mat(glm::ortho(params.frustum.l, params.frustum.r, params.frustum.b, params.frustum.t));
+	cv::Mat mvp = ortho_projection * view_model;
+
+	glm::vec4 viewport(0, height, width, -height); // flips y, origin top-left, like in OpenCV
+	// equivalent to what glm::project's viewport does, but we don't change z and w:
+	cv::Mat viewport_mat = (cv::Mat_<float>(4, 4) << viewport[2] / 2.0f, 0.0f,       0.0f, viewport[2] / 2.0f + viewport[0],
+												     0.0f,               viewport[3] / 2.0f, 0.0f, viewport[3] / 2.0f + viewport[1],
+													 0.0f,               0.0f,               1.0f, 0.0f,
+													 0.0f,               0.0f,               0.0f, 1.0f);
+
+	cv::Mat full_projection_4x4 = viewport_mat * mvp;
+	cv::Mat full_projection_3x4 = full_projection_4x4.rowRange(0, 3); // we take the first 3 rows, but then set the last one to [0 0 0 1]
+	full_projection_3x4.at<float>(2, 0) = 0.0f;
+	full_projection_3x4.at<float>(2, 1) = 0.0f;
+	full_projection_3x4.at<float>(2, 2) = 0.0f;
+	full_projection_3x4.at<float>(2, 3) = 1.0f;
+
+	return full_projection_3x4;
+};
+
+/**
+ * @brief Returns a glm/OpenGL compatible viewport vector that flips y and
+ * has the origin on the top-left, like in OpenCV.
+ *
+ * Note: Move to detail namespace / not used at the moment.
+ */
+glm::vec4 get_opencv_viewport(int width, int height)
+{
+	return glm::vec4(0, height, width, -height);
+};
+
 /**
  * @brief This algorithm estimates the rotation angles and translation of the model, as
  * well as the viewing frustum of the camera, given a set of corresponding 2D-3D points.
@@ -99,7 +185,7 @@ struct RenderingParameters
  * @param[in] height Height of the image (or viewport).
  * @return The estimated model and camera parameters.
  */
-RenderingParameters estimate_orthographic_camera(std::vector<cv::Vec2f> image_points, std::vector<cv::Vec4f> model_points, int width, int height)
+OrthographicRenderingParameters estimate_orthographic_camera(std::vector<cv::Vec2f> image_points, std::vector<cv::Vec4f> model_points, int width, int height)
 {
 	using cv::Mat;
 	assert(image_points.size() == model_points.size());
@@ -122,8 +208,8 @@ RenderingParameters estimate_orthographic_camera(std::vector<cv::Vec2f> image_po
 	auto info = lm.minimize(parameters); // we could or should use the return value
 	// 'parameters' contains the solution now.
 
-	Frustum camera_frustum{ -1.0f * aspect * parameters[5], 1.0f * aspect * parameters[5], -1.0f * parameters[5], 1.0f * parameters[5] };
-	return RenderingParameters{ static_cast<float>(parameters[0]), static_cast<float>(parameters[1]), static_cast<float>(parameters[2]), static_cast<float>(parameters[3]), static_cast<float>(parameters[4]), camera_frustum };
+	Frustum camera_frustum{ -1.0f * aspect * static_cast<float>(parameters[5]), 1.0f * aspect * static_cast<float>(parameters[5]), -1.0f * static_cast<float>(parameters[5]), 1.0f * static_cast<float>(parameters[5]) };
+	return OrthographicRenderingParameters{ static_cast<float>(parameters[0]), static_cast<float>(parameters[1]), static_cast<float>(parameters[2]), static_cast<float>(parameters[3]), static_cast<float>(parameters[4]), camera_frustum };
 };
 
 	} /* namespace fitting */

include/eos/render/Mesh.hpp

@@ -49,7 +49,7 @@ namespace eos {
 struct Mesh
 {
 	std::vector<cv::Vec4f> vertices; ///< 3D vertex positions.
-	std::vector<cv::Vec3f> colors; ///< Color information for each vertex. Expected to be in RGB order.
+	std::vector<cv::Vec3f> colors; ///< Colour information for each vertex. Expected to be in RGB order.
 	std::vector<cv::Vec2f> texcoords; ///< Texture coordinates for each vertex.
 
 	std::vector<std::array<int, 3>> tvi; ///< Triangle vertex indices
@@ -145,7 +145,7 @@ inline void write_textured_obj(Mesh mesh, std::string filename)
 	mtl_file << "map_Kd " << texture_filename.string() << std::endl;
 
 	return;
-}
+};
 
 	} /* namespace render */
 } /* namespace eos */

include/eos/render/detail/render_affine_detail.hpp

+/*
+ * Eos - A 3D Morphable Model fitting library written in modern C++11/14.
+ *
+ * File: include/eos/render/detail/render_affine_detail.hpp
+ *
+ * Copyright 2014, 2015 Patrik Huber
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#pragma once
+
+#ifndef RENDER_AFFINE_DETAIL_HPP_
+#define RENDER_AFFINE_DETAIL_HPP_
+
+#include "eos/render/detail/render_detail.hpp"
+
+#include "opencv2/core/core.hpp"
+
+/**
+ * Implementations of internal functions, not part of the
+ * API we expose and not meant to be used by a user.
+ *
+ * This file contains things specific to the affine rendering.
+ */
+namespace eos {
+	namespace render {
+		namespace detail {
+
+/**
+ * Takes a 3x4 affine camera matrix estimated with fitting::estimate_affine_camera
+ * and computes the cross product of the first two rows to create a third axis that
+ * is orthogonal to the first two.
+ * This allows us to produce z values and figure out correct depth ordering in the
+ * rendering and for texture extraction.
+ *
+ * @param[in] affine_camera_matrix A 3x4 affine camera matrix.
+ * @return The matrix with a third row inserted.
+ */
+cv::Mat calculate_affine_z_direction(cv::Mat affine_camera_matrix)
+{
+	using cv::Mat;
+	// Take the cross product of row 0 with row 1 to get the direction perpendicular to the viewing plane (= the viewing direction).
+	// Todo: We should check if we look/project into the right direction - the sign could be wrong?
+	Mat affine_cam_z_rotation = affine_camera_matrix.row(0).colRange(0, 3).cross(affine_camera_matrix.row(1).colRange(0, 3));
+	affine_cam_z_rotation /= cv::norm(affine_cam_z_rotation, cv::NORM_L2);
+
+	// The 4x4 affine camera matrix
+	Mat affine_cam_4x4 = Mat::zeros(4, 4, CV_32FC1);
+
+	// Replace the third row with the camera-direction (z)
+	Mat third_row_rotation_part = affine_cam_4x4.row(2).colRange(0, 3);
+	affine_cam_z_rotation.copyTo(third_row_rotation_part); // Set first 3 components. 4th component stays 0.
+
+	// Copy the first 2 rows from the input matrix
+	Mat first_two_rows_of_4x4 = affine_cam_4x4.rowRange(0, 2);
+	affine_camera_matrix.rowRange(0, 2).copyTo(first_two_rows_of_4x4);
+
+	// The 4th row is (0, 0, 0, 1):
+	affine_cam_4x4.at<float>(3, 3) = 1.0f;
+
+	return affine_cam_4x4;
+};
+
+/**
+ * Rasters a triangle into the given colour and depth buffer.
+ *
+ * In essence, loop through the pixels inside the triangle's bounding
+ * box, calculate the barycentric coordinates, and if inside the triangle
+ * and the z-test is passed, then draw the point using the barycentric
+ * coordinates for colour interpolation.
+
+ * Does not do perspective-correct weighting, and therefore only works
+ * with the affine rendering pipeline.
+ *
+ * No texturing at the moment.
+ *
+ * Note/Todo: See where and how this is used, and how similar it is to
+ * the "normal" raster_triangle. Maybe rename to raster_triangle_vertexcolour?
+ *
+ * @param[in] triangle A triangle.
+ * @param[in] colourbuffer The colour buffer to draw into.
+ * @param[in] depthbuffer The depth buffer to draw into and use for the depth test.
+ */
+void raster_triangle_affine(TriangleToRasterize triangle, cv::Mat colourbuffer, cv::Mat depthbuffer)
+{
+	for (int yi = triangle.min_y; yi <= triangle.max_y; ++yi)
+	{
+		for (int xi = triangle.min_x; xi <= triangle.max_x; ++xi)
+		{
+			// we want centers of pixels to be used in computations. Todo: Do we?
+			const float x = static_cast<float>(xi) + 0.5f;
+			const float y = static_cast<float>(yi) + 0.5f;
+
+			// these will be used for barycentric weights computation
+			const double one_over_v0ToLine12 = 1.0 / implicit_line(triangle.v0.position[0], triangle.v0.position[1], triangle.v1.position, triangle.v2.position);
+			const double one_over_v1ToLine20 = 1.0 / implicit_line(triangle.v1.position[0], triangle.v1.position[1], triangle.v2.position, triangle.v0.position);
+			const double one_over_v2ToLine01 = 1.0 / implicit_line(triangle.v2.position[0], triangle.v2.position[1], triangle.v0.position, triangle.v1.position);
+			// affine barycentric weights
+			const double alpha = implicit_line(x, y, triangle.v1.position, triangle.v2.position) * one_over_v0ToLine12;
+			const double beta = implicit_line(x, y, triangle.v2.position, triangle.v0.position) * one_over_v1ToLine20;
+			const double gamma = implicit_line(x, y, triangle.v0.position, triangle.v1.position) * one_over_v2ToLine01;
+
+			// if pixel (x, y) is inside the triangle or on one of its edges
+			if (alpha >= 0 && beta >= 0 && gamma >= 0)
+			{
+				const int pixel_index_row = yi;
+				const int pixel_index_col = xi;
+
+				const double z_affine = alpha*static_cast<double>(triangle.v0.position[2]) + beta*static_cast<double>(triangle.v1.position[2]) + gamma*static_cast<double>(triangle.v2.position[2]);
+				if (z_affine < depthbuffer.at<double>(pixel_index_row, pixel_index_col))
+				{
+					// attributes interpolation
+					// pixel_color is in RGB, v.color are RGB
+					cv::Vec3f pixel_color = alpha*triangle.v0.color + beta*triangle.v1.color + gamma*triangle.v2.color;
+
+					// clamp bytes to 255
+					const unsigned char red = static_cast<unsigned char>(255.0f * std::min(pixel_color[0], 1.0f)); // Todo: Proper casting (rounding?)
+					const unsigned char green = static_cast<unsigned char>(255.0f * std::min(pixel_color[1], 1.0f));
+					const unsigned char blue = static_cast<unsigned char>(255.0f * std::min(pixel_color[2], 1.0f));
+
+					// update buffers
+					colourbuffer.at<cv::Vec4b>(pixel_index_row, pixel_index_col)[0] = blue;
+					colourbuffer.at<cv::Vec4b>(pixel_index_row, pixel_index_col)[1] = green;
+					colourbuffer.at<cv::Vec4b>(pixel_index_row, pixel_index_col)[2] = red;
+					colourbuffer.at<cv::Vec4b>(pixel_index_row, pixel_index_col)[3] = 255; // alpha channel
+					depthbuffer.at<double>(pixel_index_row, pixel_index_col) = z_affine;
+				}
+			}
+		}
+	}
+};
+
+		} /* namespace detail */
+	} /* namespace render */
+} /* namespace eos */
+
+#endif /* RENDER_AFFINE_DETAIL_HPP_ */

include/eos/render/detail/render_detail.hpp

@@ -32,61 +32,95 @@ namespace eos {
 	namespace render {
 		namespace detail {
 
-/**
- * Takes a 3x4 affine camera matrix estimated with fitting::estimate_affine_camera
- * and computes the cross product of the first two rows to create a third axis that
- * is orthogonal to the first two.
- * This allows us to produce z values and figure out correct depth ordering in the
- * rendering and for texture extraction.
- *
- * @param[in] affine_camera_matrix A 3x4 affine camera matrix.
- * @return The matrix with a third row inserted.
- */
-cv::Mat calculate_affine_z_direction(cv::Mat affine_camera_matrix)
-{
-	using cv::Mat;
-	// Take the cross product of row 0 with row 1 to get the direction perpendicular to the viewing plane (= the viewing direction).
-	// Todo: We should check if we look/project into the right direction - the sign could be wrong?
-	Mat affine_cam_z_rotation = affine_camera_matrix.row(0).colRange(0, 3).cross(affine_camera_matrix.row(1).colRange(0, 3));
-	affine_cam_z_rotation /= cv::norm(affine_cam_z_rotation, cv::NORM_L2);
-
-	// The 4x4 affine camera matrix
-	Mat affine_cam_4x4 = Mat::zeros(4, 4, CV_32FC1);
-
-	// Replace the third row with the camera-direction (z)
-	Mat third_row_rotation_part = affine_cam_4x4.row(2).colRange(0, 3);
-	affine_cam_z_rotation.copyTo(third_row_rotation_part); // Set first 3 components. 4th component stays 0.
-
-	// Copy the first 2 rows from the input matrix
-	Mat first_two_rows_of_4x4 = affine_cam_4x4.rowRange(0, 2);
-	affine_camera_matrix.rowRange(0, 2).copyTo(first_two_rows_of_4x4);
-
-	// The 4th row is (0, 0, 0, 1):
-	affine_cam_4x4.at<float>(3, 3) = 1.0f;
-
-	return affine_cam_4x4;
-};
-
 /**
  * Just a representation for a vertex during rendering.
  *
  * Might consider getting rid of it.
+ * Used in render_affine and render.
  */
 class Vertex
 {
 public:
 	Vertex() {};
-	Vertex(const cv::Vec4f& position, const cv::Vec3f& color, const cv::Vec2f& texCoord) : position(position), color(color), texcrd(texCoord) {};
+	Vertex(const cv::Vec4f& position, const cv::Vec3f& color, const cv::Vec2f& texcoords) : position(position), color(color), texcoords(texcoords) {};
 
 	cv::Vec4f position;
-	cv::Vec3f color;
-	cv::Vec2f texcrd;
+	cv::Vec3f color; ///< in RGB order
+	cv::Vec2f texcoords;
+};
+
+class plane
+{
+public:
+	plane() {}
+
+	plane(float a, float b, float c, float d)
+	{
+		this->a = a;
+		this->b = b;
+		this->c = c;
+		this->d = d;
+	}
+
+	plane(const cv::Vec3f& normal, float d = 0.0f)
+	{
+		this->a = normal[0];
+		this->b = normal[1];
+		this->c = normal[2];
+		this->d = d;
+	}
+
+	plane(const cv::Vec3f& point, const cv::Vec3f& normal)
+	{
+		a = normal[0];
+		b = normal[1];
+		c = normal[2];
+		d = -(point.dot(normal));
+	}
+
+	plane(const cv::Vec3f& point1, const cv::Vec3f& point2, const cv::Vec3f& point3)
+	{
+		cv::Vec3f v1 = point2 - point1;
+		cv::Vec3f v2 = point3 - point1;
+		cv::Vec3f normal = (v1.cross(v2));
+		normal /= cv::norm(normal, cv::NORM_L2);
+
+		a = normal[0];
+		b = normal[1];
+		c = normal[2];
+		d = -(point1.dot(normal));
+	}
+
+	void normalize()
+	{
+		float length = sqrt(a*a + b*b + c*c);
+
+		a /= length;
+		b /= length;
+		c /= length;
+	}
+
+	float getSignedDistanceFromPoint(const cv::Vec3f& point) const
+	{
+		return a*point[0] + b*point[1] + c*point[2] + d;
+	}
+
+	float getSignedDistanceFromPoint(const cv::Vec4f& point) const
+	{
+		return a*point[0] + b*point[1] + c*point[2] + d;
+	}
+
+public:
+	float a, b, c;
+	float d;
 };
 
 /**
  * A representation for a triangle that is to be rasterised.
  * Stores the enclosing bounding box of the triangle that is
  * calculated during rendering and used during rasterisation.
+ *
+ * Used in render_affine and render.
  */
 struct TriangleToRasterize
 {
@@ -95,6 +129,26 @@ struct TriangleToRasterize
 	int max_x;
 	int min_y;
 	int max_y;
+	// Everything below is only used in the "normal" renderer, but not
+	// in render_affine.
+	double one_over_z0;
+	double one_over_z1;
+	double one_over_z2;
+	//double one_over_v0ToLine12;
+	//double one_over_v1ToLine20;
+	//double one_over_v2ToLine01;
+	plane alphaPlane;
+	plane betaPlane;
+	plane gammaPlane;
+	double one_over_alpha_c; // those are only used for texturing -> float
+	double one_over_beta_c;
+	double one_over_gamma_c;
+	float alpha_ffx;
+	float beta_ffx;
+	float gamma_ffx;
+	float alpha_ffy;
+	float beta_ffy;
+	float gamma_ffy;
 };
 
 /**
@@ -112,11 +166,17 @@ struct TriangleToRasterize
  */
 cv::Rect calculate_clipped_bounding_box(cv::Vec4f v0, cv::Vec4f v1, cv::Vec4f v2, int viewport_width, int viewport_height)
 {
+	/* Old, producing artifacts:
+	t.minX = max(min(t.v0.position[0], min(t.v1.position[0], t.v2.position[0])), 0.0f);
+	t.maxX = min(max(t.v0.position[0], max(t.v1.position[0], t.v2.position[0])), (float)(viewportWidth - 1));
+	t.minY = max(min(t.v0.position[1], min(t.v1.position[1], t.v2.position[1])), 0.0f);
+	t.maxY = min(max(t.v0.position[1], max(t.v1.position[1], t.v2.position[1])), (float)(viewportHeight - 1));*/
+
 	using std::min;
 	using std::max;
 	using std::floor;
 	using std::ceil;
-	int minX = max(min(floor(v0[0]), min(floor(v1[0]), floor(v2[0]))), 0.0f);
+	int minX = max(min(floor(v0[0]), min(floor(v1[0]), floor(v2[0]))), 0.0f); // Readded this comment after merge: What about rounding, or rather the conversion from double to int?
 	int maxX = min(max(ceil(v0[0]), max(ceil(v1[0]), ceil(v2[0]))), static_cast<float>(viewport_width - 1));
 	int minY = max(min(floor(v0[1]), min(floor(v1[1]), floor(v2[1]))), 0.0f);
 	int maxY = min(max(ceil(v0[1]), max(ceil(v1[1]), ceil(v2[1]))), static_cast<float>(viewport_height - 1));
@@ -148,30 +208,240 @@ double implicit_line(float x, float y, const cv::Vec4f& v1, const cv::Vec4f& v2)
 	return ((double)v1[1] - (double)v2[1])*(double)x + ((double)v2[0] - (double)v1[0])*(double)y + (double)v1[0] * (double)v2[1] - (double)v2[0] * (double)v1[1];
 };
 
-/**
- * Rasters a triangle into the given colour and depth buffer.
- *
- * In essence, loop through the pixels inside the triangle's bounding
- * box, calculate the barycentric coordinates, and if inside the triangle
- * and the z-test is passed, then draw the point using the barycentric
- * coordinates for colour interpolation.
+std::vector<Vertex> clip_polygon_to_plane_in_4d(const std::vector<Vertex>& vertices, const cv::Vec4f& plane_normal)
+{
+	std::vector<Vertex> clippedVertices;
 
- * Does not do perspective-correct weighting, and therefore only works
- * with the affine rendering pipeline.
- *
- * No texturing at the moment.
- *
- * @param[in] triangle A triangle.
- * @param[in] colourbuffer The colour buffer to draw into.
- * @param[in] depthbuffer The depth buffer to draw into and use for the depth test.
+	// We can have 2 cases:
+	//	* 1 vertex visible: we make 1 new triangle out of the visible vertex plus the 2 intersection points with the near-plane
+	//  * 2 vertices visible: we have a quad, so we have to make 2 new triangles out of it.
+
+	// See here for more info? http://math.stackexchange.com/questions/400268/equation-for-a-line-through-a-plane-in-homogeneous-coordinates
+
+	for (unsigned int i = 0; i < vertices.size(); i++)
+	{
+		int a = i; // the current vertex
+		int b = (i + 1) % vertices.size(); // the following vertex (wraps around 0)
+
+		float fa = vertices[a].position.dot(plane_normal); // Note: Shouldn't they be unit length?
+		float fb = vertices[b].position.dot(plane_normal); // < 0 means on visible side, > 0 means on invisible side?
+
+		if ((fa < 0 && fb > 0) || (fa > 0 && fb < 0)) // one vertex is on the visible side of the plane, one on the invisible? so we need to split?
+		{
+			cv::Vec4f direction = vertices[b].position - vertices[a].position;
+			float t = -(plane_normal.dot(vertices[a].position)) / (plane_normal.dot(direction)); // the parametric value on the line, where the line to draw intersects the plane?
+
+			// generate a new vertex at the line-plane intersection point
+			cv::Vec4f position = vertices[a].position + t*direction;
+			cv::Vec3f color = vertices[a].color + t*(vertices[b].color - vertices[a].color);
+			cv::Vec2f texCoord = vertices[a].texcoords + t*(vertices[b].texcoords - vertices[a].texcoords);	// We could omit that if we don't render with texture.
+
+			if (fa < 0) // we keep the original vertex plus the new one
+			{
+				clippedVertices.push_back(vertices[a]);
+				clippedVertices.push_back(Vertex(position, color, texCoord));
+			}
+			else if (fb < 0) // we use only the new vertex
+			{
+				clippedVertices.push_back(Vertex(position, color, texCoord));
+			}
+		}
+		else if (fa < 0 && fb < 0) // both are visible (on the "good" side of the plane), no splitting required, use the current vertex
+		{
+			clippedVertices.push_back(vertices[a]);
+		}
+		// else, both vertices are not visible, nothing to add and draw
+	}
+
+	return clippedVertices;
+};
+
+// used only in tex2D_linear_mipmap_linear
+// template?
+float clamp(float x, float a, float b)
+{
+	return std::max(std::min(x, b), a);
+};
+
+cv::Vec2f texcoord_wrap(const cv::Vec2f& texcoords)
+{
+	return cv::Vec2f(texcoords[0] - (int)texcoords[0], texcoords[1] - (int)texcoords[1]);
+};
+
+// forward decls
+cv::Vec3f tex2d_linear_mipmap_linear(const cv::Vec2f& texcoords, const Texture& texture, float dudx, float dudy, float dvdx, float dvdy);
+cv::Vec3f tex2d_linear(const cv::Vec2f& imageTexCoord, unsigned char mipmapIndex, const Texture& texture);
+
+cv::Vec3f tex2d(const cv::Vec2f& texcoords, const Texture& texture, float dudx, float dudy, float dvdx, float dvdy)
+{
+	return (1.0f / 255.0f) * tex2d_linear_mipmap_linear(texcoords, texture, dudx, dudy, dvdx, dvdy);
+};
+
+cv::Vec3f tex2d_linear_mipmap_linear(const cv::Vec2f& texcoords, const Texture& texture, float dudx, float dudy, float dvdx, float dvdy)
+{
+	using cv::Vec2f;
+	float px = std::sqrt(std::pow(dudx, 2) + std::pow(dvdx, 2));
+	float py = std::sqrt(std::pow(dudy, 2) + std::pow(dvdy, 2));
+	float lambda = std::log(std::max(px, py)) / CV_LOG2;
+	unsigned char mipmapIndex1 = detail::clamp((int)lambda, 0.0f, std::max(texture.widthLog, texture.heightLog) - 1);
+	unsigned char mipmapIndex2 = mipmapIndex1 + 1;
+
+	Vec2f imageTexCoord = detail::texcoord_wrap(texcoords);
+	Vec2f imageTexCoord1 = imageTexCoord;
+	imageTexCoord1[0] *= texture.mipmaps[mipmapIndex1].cols;
+	imageTexCoord1[1] *= texture.mipmaps[mipmapIndex1].rows;
+	Vec2f imageTexCoord2 = imageTexCoord;
+	imageTexCoord2[0] *= texture.mipmaps[mipmapIndex2].cols;
+	imageTexCoord2[1] *= texture.mipmaps[mipmapIndex2].rows;
+
+	cv::Vec3f color, color1, color2;
+	color1 = tex2d_linear(imageTexCoord1, mipmapIndex1, texture);
+	color2 = tex2d_linear(imageTexCoord2, mipmapIndex2, texture);
+	float lambdaFrac = std::max(lambda, 0.0f);
+	lambdaFrac = lambdaFrac - (int)lambdaFrac;
+	color = (1.0f - lambdaFrac)*color1 + lambdaFrac*color2;
+
+	return color;
+};
+
+cv::Vec3f tex2d_linear(const cv::Vec2f& imageTexCoord, unsigned char mipmap_index, const Texture& texture)
+{
+	int x = (int)imageTexCoord[0];
+	int y = (int)imageTexCoord[1];
+	float alpha = imageTexCoord[0] - x;
+	float beta = imageTexCoord[1] - y;
+	float oneMinusAlpha = 1.0f - alpha;
+	float oneMinusBeta = 1.0f - beta;
+	float a = oneMinusAlpha * oneMinusBeta;
+	float b = alpha * oneMinusBeta;
+	float c = oneMinusAlpha * beta;
+	float d = alpha * beta;
+	cv::Vec3f color;
+
+	using cv::Vec4b;
+	//int pixelIndex;
+	//pixelIndex = getPixelIndex_wrap(x, y, texture->mipmaps[mipmapIndex].cols, texture->mipmaps[mipmapIndex].rows);
+	int pixelIndexCol = x; if (pixelIndexCol == texture.mipmaps[mipmap_index].cols) { pixelIndexCol = 0; }
+	int pixelIndexRow = y; if (pixelIndexRow == texture.mipmaps[mipmap_index].rows) { pixelIndexRow = 0; }
+	//std::cout << texture.mipmaps[mipmapIndex].cols << " " << texture.mipmaps[mipmapIndex].rows << " " << texture.mipmaps[mipmapIndex].channels() << std::endl;
+	//cv::imwrite("mm.png", texture.mipmaps[mipmapIndex]);
+	color[0] = a * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[0];
+	color[1] = a * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[1];
+	color[2] = a * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[2];
+
+	//pixelIndex = getPixelIndex_wrap(x + 1, y, texture.mipmaps[mipmapIndex].cols, texture.mipmaps[mipmapIndex].rows);
+	pixelIndexCol = x + 1; if (pixelIndexCol == texture.mipmaps[mipmap_index].cols) { pixelIndexCol = 0; }
+	pixelIndexRow = y; if (pixelIndexRow == texture.mipmaps[mipmap_index].rows) { pixelIndexRow = 0; }
+	color[0] += b * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[0];
+	color[1] += b * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[1];
+	color[2] += b * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[2];
+
+	//pixelIndex = getPixelIndex_wrap(x, y + 1, texture.mipmaps[mipmapIndex].cols, texture.mipmaps[mipmapIndex].rows);
+	pixelIndexCol = x; if (pixelIndexCol == texture.mipmaps[mipmap_index].cols) { pixelIndexCol = 0; }
+	pixelIndexRow = y + 1; if (pixelIndexRow == texture.mipmaps[mipmap_index].rows) { pixelIndexRow = 0; }
+	color[0] += c * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[0];
+	color[1] += c * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[1];
+	color[2] += c * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[2];
+
+	//pixelIndex = getPixelIndex_wrap(x + 1, y + 1, texture.mipmaps[mipmapIndex].cols, texture.mipmaps[mipmapIndex].rows);
+	pixelIndexCol = x + 1; if (pixelIndexCol == texture.mipmaps[mipmap_index].cols) { pixelIndexCol = 0; }
+	pixelIndexRow = y + 1; if (pixelIndexRow == texture.mipmaps[mipmap_index].rows) { pixelIndexRow = 0; }
+	color[0] += d * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[0];
+	color[1] += d * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[1];
+	color[2] += d * texture.mipmaps[mipmap_index].at<Vec4b>(pixelIndexRow, pixelIndexCol)[2];
+
+	return color;
+};
+
+// Todo: Split this function into the general (core-part) and the texturing part.
+// Then, utils::extractTexture can re-use the core-part.
+// Note: Maybe a bit outdated "todo" above.
+boost::optional<TriangleToRasterize> process_prospective_tri(Vertex v0, Vertex v1, Vertex v2, int viewport_width, int viewport_height, bool enable_backface_culling)
+{
+	using cv::Vec2f;
+	using cv::Vec3f;
+	TriangleToRasterize t;
+	t.v0 = v0;	// no memcopy I think. the transformed vertices don't get copied and exist only once. They are a local variable in runVertexProcessor(), the ref is passed here, and if we need to rasterize it, it gets push_back'ed (=copied?) to trianglesToRasterize. Perfect I think. TODO: Not anymore, no ref here
+	t.v1 = v1;
+	t.v2 = v2;
+
+	// Only for texturing or perspective texturing:
+	//t.texture = _texture;
+	t.one_over_z0 = 1.0 / (double)t.v0.position[3];
+	t.one_over_z1 = 1.0 / (double)t.v1.position[3];
+	t.one_over_z2 = 1.0 / (double)t.v2.position[3];
+
+	// divide by w
+	// if ortho, we can do the divide as well, it will just be a / 1.0f.
+	t.v0.position = t.v0.position / t.v0.position[3];
+	t.v1.position = t.v1.position / t.v1.position[3];
+	t.v2.position = t.v2.position / t.v2.position[3];
+
+	// project from 4D to 2D window position with depth value in z coordinate
+	// Viewport transform:
+	/* (a possible optimisation might be to use matrix multiplication for this as well
+	   and do it for all triangles at once? See 'windowTransform' in:
+	   https://github.com/elador/FeatureDetection/blob/964f0b2107ce73ef2f06dc829e5084be421de5a5/libRender/src/render/RenderDevice.cpp)
 	*/
-void raster_triangle(TriangleToRasterize triangle, cv::Mat colourbuffer, cv::Mat depthbuffer)
+	Vec2f v0_screen = eos::render::clip_to_screen_space(Vec2f(t.v0.position[0], t.v0.position[1]), viewport_width, viewport_height);
+	t.v0.position[0] = v0_screen[0];
+	t.v0.position[1] = v0_screen[1];
+	Vec2f v1_screen = clip_to_screen_space(Vec2f(t.v1.position[0], t.v1.position[1]), viewport_width, viewport_height);
+	t.v1.position[0] = v1_screen[0];
+	t.v1.position[1] = v1_screen[1];
+	Vec2f v2_screen = clip_to_screen_space(Vec2f(t.v2.position[0], t.v2.position[1]), viewport_width, viewport_height);
+	t.v2.position[0] = v2_screen[0];
+	t.v2.position[1] = v2_screen[1];
+
+	if (enable_backface_culling) {
+		if (!are_vertices_ccw_in_screen_space(t.v0.position, t.v1.position, t.v2.position))
+			return boost::none;
+	}
+
+	// Get the bounding box of the triangle:
+	cv::Rect boundingBox = calculate_clipped_bounding_box(t.v0.position, t.v1.position, t.v2.position, viewport_width, viewport_height);
+	t.min_x = boundingBox.x;
+	t.max_x = boundingBox.x + boundingBox.width;
+	t.min_y = boundingBox.y;
+	t.max_y = boundingBox.y + boundingBox.height;
+
+	if (t.max_x <= t.min_x || t.max_y <= t.min_y) 	// Note: Can the width/height of the bbox be negative? Maybe we only need to check for equality here?
+		return boost::none;
+
+	// Which of these is for texturing, mipmapping, what for perspective?
+	// for partial derivatives computation
+	t.alphaPlane = plane(Vec3f(t.v0.position[0], t.v0.position[1], t.v0.texcoords[0] * t.one_over_z0),
+		Vec3f(t.v1.position[0], t.v1.position[1], t.v1.texcoords[0] * t.one_over_z1),
+		Vec3f(t.v2.position[0], t.v2.position[1], t.v2.texcoords[0] * t.one_over_z2));
+	t.betaPlane = plane(Vec3f(t.v0.position[0], t.v0.position[1], t.v0.texcoords[1] * t.one_over_z0),
+		Vec3f(t.v1.position[0], t.v1.position[1], t.v1.texcoords[1] * t.one_over_z1),
+		Vec3f(t.v2.position[0], t.v2.position[1], t.v2.texcoords[1] * t.one_over_z2));
+	t.gammaPlane = plane(Vec3f(t.v0.position[0], t.v0.position[1], t.one_over_z0),
+		Vec3f(t.v1.position[0], t.v1.position[1], t.one_over_z1),
+		Vec3f(t.v2.position[0], t.v2.position[1], t.one_over_z2));
+	t.one_over_alpha_c = 1.0f / t.alphaPlane.c;
+	t.one_over_beta_c = 1.0f / t.betaPlane.c;
+	t.one_over_gamma_c = 1.0f / t.gammaPlane.c;
+	t.alpha_ffx = -t.alphaPlane.a * t.one_over_alpha_c;
+	t.beta_ffx = -t.betaPlane.a * t.one_over_beta_c;
+	t.gamma_ffx = -t.gammaPlane.a * t.one_over_gamma_c;
+	t.alpha_ffy = -t.alphaPlane.b * t.one_over_alpha_c;
+	t.beta_ffy = -t.betaPlane.b * t.one_over_beta_c;
+	t.gamma_ffy = -t.gammaPlane.b * t.one_over_gamma_c;
+
+	// Use t
+	return boost::optional<TriangleToRasterize>(t);
+};
+
+void raster_triangle(TriangleToRasterize triangle, cv::Mat colourbuffer, cv::Mat depthbuffer, boost::optional<Texture> texture, bool enable_far_clipping)
 {
-	for (int yi = triangle.min_y; yi <= triangle.max_y; yi++)
+	using cv::Vec2f;
+	using cv::Vec3f;
+	for (int yi = triangle.min_y; yi <= triangle.max_y; ++yi)
 	{
-		for (int xi = triangle.min_x; xi <= triangle.max_x; xi++)
+		for (int xi = triangle.min_x; xi <= triangle.max_x; ++xi)
 		{
-			// we want centers of pixels to be used in computations. Do we?
+			// we want centers of pixels to be used in computations. Todo: Do we?
 			const float x = static_cast<float>(xi) + 0.5f;
 			const float y = static_cast<float>(yi) + 0.5f;
 
@@ -180,9 +450,9 @@ void raster_triangle(TriangleToRasterize triangle, cv::Mat colourbuffer, cv::Mat
 			const double one_over_v1ToLine20 = 1.0 / implicit_line(triangle.v1.position[0], triangle.v1.position[1], triangle.v2.position, triangle.v0.position);
 			const double one_over_v2ToLine01 = 1.0 / implicit_line(triangle.v2.position[0], triangle.v2.position[1], triangle.v0.position, triangle.v1.position);
 			// affine barycentric weights
-			const double alpha = implicit_line(x, y, triangle.v1.position, triangle.v2.position) * one_over_v0ToLine12;
-			const double beta = implicit_line(x, y, triangle.v2.position, triangle.v0.position) * one_over_v1ToLine20;
-			const double gamma = implicit_line(x, y, triangle.v0.position, triangle.v1.position) * one_over_v2ToLine01;
+			double alpha = implicit_line(x, y, triangle.v1.position, triangle.v2.position) * one_over_v0ToLine12;
+			double beta = implicit_line(x, y, triangle.v2.position, triangle.v0.position) * one_over_v1ToLine20;
+			double gamma = implicit_line(x, y, triangle.v0.position, triangle.v1.position) * one_over_v2ToLine01;
 
 			// if pixel (x, y) is inside the triangle or on one of its edges
 			if (alpha >= 0 && beta >= 0 && gamma >= 0)
@@ -191,11 +461,61 @@ void raster_triangle(TriangleToRasterize triangle, cv::Mat colourbuffer, cv::Mat
 				const int pixel_index_col = xi;
 
 				const double z_affine = alpha*static_cast<double>(triangle.v0.position[2]) + beta*static_cast<double>(triangle.v1.position[2]) + gamma*static_cast<double>(triangle.v2.position[2]);
-				if (z_affine < depthbuffer.at<double>(pixel_index_row, pixel_index_col))
+				// The '<= 1.0' clips against the far-plane in NDC. We clip against the near-plane earlier.
+				// TODO: Use enable_far_clipping here.
+				bool draw = true;
+				if (enable_far_clipping)
+				{
+					if (z_affine > 1.0)
+					{
+						draw = false;
+					}
+				}
+				//if (z_affine < depthbuffer.at<double>(pixelIndexRow, pixelIndexCol)/* && z_affine <= 1.0*/) // what to do in ortho case without n/f "squashing"? should we always squash? or a flag?
+				if (z_affine < depthbuffer.at<double>(pixel_index_row, pixel_index_col) && draw)
 				{
+					// perspective-correct barycentric weights
+					double d = alpha*triangle.one_over_z0 + beta*triangle.one_over_z1 + gamma*triangle.one_over_z2;
+					d = 1.0 / d;
+					alpha *= d*triangle.one_over_z0; // In case of affine cam matrix, everything is 1 and a/b/g don't get changed.
+					beta *= d*triangle.one_over_z1;
+					gamma *= d*triangle.one_over_z2;
+
 					// attributes interpolation
-					// pixel_color is in RGB, v.color are RGB
-					cv::Vec3f pixel_color = alpha*triangle.v0.color + beta*triangle.v1.color + gamma*triangle.v2.color;
+					Vec3f color_persp = alpha*triangle.v0.color + beta*triangle.v1.color + gamma*triangle.v2.color;
+					Vec2f texcoords_persp = alpha*triangle.v0.texcoords + beta*triangle.v1.texcoords + gamma*triangle.v2.texcoords;
+
+					Vec3f pixel_color;
+					// Pixel Shader:
+					if (texture) {	// We use texturing
+						// check if texture != NULL?
+						// partial derivatives (for mip-mapping)
+						const float u_over_z = -(triangle.alphaPlane.a*x + triangle.alphaPlane.b*y + triangle.alphaPlane.d) * triangle.one_over_alpha_c;
+						const float v_over_z = -(triangle.betaPlane.a*x + triangle.betaPlane.b*y + triangle.betaPlane.d) * triangle.one_over_beta_c;
+						const float one_over_z = -(triangle.gammaPlane.a*x + triangle.gammaPlane.b*y + triangle.gammaPlane.d) * triangle.one_over_gamma_c;
+						const float one_over_squared_one_over_z = 1.0f / std::pow(one_over_z, 2);
+
+						// partial derivatives of U/V coordinates with respect to X/Y pixel's screen coordinates
+						float dudx = one_over_squared_one_over_z * (triangle.alpha_ffx * one_over_z - u_over_z * triangle.gamma_ffx);
+						float dudy = one_over_squared_one_over_z * (triangle.beta_ffx * one_over_z - v_over_z * triangle.gamma_ffx);
+						float dvdx = one_over_squared_one_over_z * (triangle.alpha_ffy * one_over_z - u_over_z * triangle.gamma_ffy);
+						float dvdy = one_over_squared_one_over_z * (triangle.beta_ffy * one_over_z - v_over_z * triangle.gamma_ffy);
+
+						dudx *= texture.get().mipmaps[0].cols;
+						dudy *= texture.get().mipmaps[0].cols;
+						dvdx *= texture.get().mipmaps[0].rows;
+						dvdy *= texture.get().mipmaps[0].rows;
+
+						// The Texture is in BGR, thus tex2D returns BGR
+						Vec3f texture_color = detail::tex2d(texcoords_persp, texture.get(), dudx, dudy, dvdx, dvdy); // uses the current texture
+						pixel_color = Vec3f(texture_color[2], texture_color[1], texture_color[0]);
+						// other: color.mul(tex2D(texture, texCoord));
+						// Old note: for texturing, we load the texture as BGRA, so the colors get the wrong way in the next few lines...
+					}
+					else {	// We use vertex-coloring
+						// color_persp is in RGB
+						pixel_color = color_persp;
+					}
 
 					// clamp bytes to 255
 					const unsigned char red = static_cast<unsigned char>(255.0f * std::min(pixel_color[0], 1.0f)); // Todo: Proper casting (rounding?)

include/eos/render/detail/texture_extraction_detail.hpp

@@ -73,6 +73,12 @@ inline bool is_point_in_triangle(cv::Point2f point, cv::Point2f triV0, cv::Point
  * The vertices should be given in screen coordinates, but with their
  * z-values preserved, so they can be compared against the depthbuffer.
  *
+ * Obviously the depthbuffer given should have been created with the same projection
+ * matrix than the texture extraction is called with.
+ *
+ * Also, we don't do perspective-correct interpolation here I think, so only
+ * use it with affine and orthographic projection matrices.
+ *
  * @param[in] v0 First vertex, in screen coordinates (but still with their z-value).
  * @param[in] v1 Second vertex.
  * @param[in] v2 Third vertex.

include/eos/render/render.hpp

+/*
+ * Eos - A 3D Morphable Model fitting library written in modern C++11/14.
+ *
+ * File: include/eos/render/render.hpp
+ *
+ * Copyright 2014, 2015 Patrik Huber
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#pragma once
+
+#ifndef RENDER_HPP_
+#define RENDER_HPP_
+
+#include "eos/render/detail/render_detail.hpp"
+#include "eos/render/utils.hpp"
+
+#include "opencv2/core/core.hpp"
+
+#ifdef WIN32
+	#define BOOST_ALL_DYN_LINK	// Link against the dynamic boost lib. Seems to be necessary because we use /MD, i.e. link to the dynamic CRT.
+	#define BOOST_ALL_NO_LIB	// Don't use the automatic library linking by boost with VS2010 (#pragma ...). Instead, we specify everything in cmake.
+#endif
+#include "boost/optional.hpp"
+
+#include <array>
+#include <vector>
+#include <memory>
+
+namespace eos {
+	namespace render {
+
+/**
+ * This file implements a software renderer conforming to OpenGL conventions. The
+ * following are implementation notes, mostly for reference, or as a reminder of
+ * what exactly is going on. Don't try to understand them :-)
+ *
+ * The renderer was initially based on code by Wojciech Sterna
+ * (http://maxest.gct-game.net/content/vainmoinen/index.html), however, it has since
+ * then been completely rewritten. Still I'd like to thank him for making his code
+ * available and bravely answering my questions via email.
+ *
+ * Coordinate systems:
+ * When specifying the vertices: +x = right, +y = up, we look into -z.
+ * So z = 0.5 is in front of 0.0.
+ * Z-Buffer:
+ *
+ * Shirley: Specify n and f with negative values. which makes sense b/c the points
+ * are along the -z axis.
+ * Consequences: notably: orthogonal(2, 3): Shirley has denominator (n-f).
+ * In what space are the points in Shirley after this?
+ * OGL: We're in the orthographic viewing volume looking down -z.
+ * However, n and f are specified positive.
+
+ * B/c the 3D points in front of the cam obviously still have negative z values, the
+ * z-value is negated. So: n = 0.1, f = 100; With the given OpenGL ortho matrix,
+ * it means a point on the near-plane which will have z = -0.1 will land up
+ * on z_clip (which equals z_ndc with ortho because w=1) = -1, and a point on
+ * the far plane z = -100 will have z_ndc = +1.
+ *
+ * That's also why in the perspective case, w_clip is set to -z_eye because
+ * to project a point the formula is $x_p = (-n * x_e)/z_e$ (because our near is
+ * specified with positive values, but the near-plane is _really_ at -n); but now we
+ * just move the minus-sign to the denominator, $x_p = (n * x_e)/-z_e$, so in the projection matrix we can use
+ * the (positive) n and f values and afterwards we divide by w = -z_e.
+ *
+ * http://www.songho.ca/opengl/gl_projectionmatrix.html
+ *
+ * Random notes:
+ * clip-space: after applying the projection matrix.
+ * ndc: after division by w
+ * NDC cube: the range of x-coordinate from [l, r] to [-1, 1], the y-coordinate from [b, t] to [-1, 1] and the z-coordinate from [n, f] to [-1, 1].
+ *
+ * Note/Todo: I read that in screen space, OpenGL transform the z-values again to be between 0 and 1?
+ *
+ * In contrast to OGL, this renderer doesn't have state, it's just a function that gets called with all
+ * necessary parameters. It's easiest for our purposes.
+ *
+ * Here's the whole rendering pipeline:
+ * Model space
+ * -> model transforms
+ * World space
+ * -> camera (view/eye) transform
+ * View / eye / camera space ("truncated pyramid frustum". In case of ortho, it's already rectangular.)
+ * -> perspective/ortho projection
+ * Clip coords (x_c, y_c, z_c, w_c); the z-axis is flipped now. z [z=-n, z=-f] is mapped to [-1, +1] in case of ortho, but not yet in case of persp (it's also flipped though), but the not-[-1,1]-range is fine as we test against w_c. I.e. the larger the z-value, the further back we are.
+ * Do frustum culling (clipping) here. Test the clip-coords with w_c, and discard if a tri is completely outside.
+ * Of the partially visible tris, clip them against the near-plane and construct the visible part of the triangle.
+ * We only do this for the near-plane here. Clipping to the near plane must be done here because after w-division triangles crossing it would get distorted.
+ * "Then, OpenGL will reconstruct the edges of the polygon where clipping occurs."
+ * -> Then divide by the w component of the clip coordinates
+ * NDC. (now only 3D vectors: [x_ndc, y_ndc, z_ndc]). nearest points have z=-1, points on far plane have z=+1.
+ * -> window transform. (also, OGL does some more to the z-buffer?)
+ * Screen / window space
+ * Directly after window-transform (still processing triangles), do backface culling with areVerticesCCWInScreenSpace()
+ * Directly afterwards we calculate the triangle's bounding box and clip x/y (screen) against 0 and the viewport width/height.
+ * Rasterising: Clipping against the far plane here by only drawing those pixels with a z-value of <= 1.0f.
+ *
+ * OGL: "both clipping (frustum culling) and NDC transformations are integrated into GL_PROJECTION matrix"
+ *
+ * Note: In both the ortho and persp case, points at z=-n end up at -1, z=-f at +1. In case of persp proj., this happens only after the divide by w.
+ */
+
+/**
+ * Renders the given mesh onto a 2D image using 4x4 model-view and
+ * projection matrices. Conforms to OpenGL conventions.
+ *
+ * @param[in] mesh A 3D mesh.
+ * @param[in] model_view_matrix A 4x4 OpenGL model-view matrix.
+ * @param[in] projection_matrix A 4x4 orthographic or perspective OpenGL projection matrix.
+ * @param[in] viewport_width Screen width.
+ * @param[in] viewport_height Screen height.
+ * @param[in] texture An optional texture map (TODO: Not optional yet!).
+ * @param[in] enable_backface_culling Whether the renderer should perform backface culling. If true, only draw triangles with vertices ordered CCW in screen-space.
+ * @param[in] enable_near_clipping Screen height.
+ * @param[in] enable_far_clipping Screen height.
+ * @return A pair with the colourbuffer as its first element and the depthbuffer as the second element.
+ */
+std::pair<cv::Mat, cv::Mat> render(Mesh mesh, cv::Mat model_view_matrix, cv::Mat projection_matrix, int viewport_width, int viewport_height, const Texture& texture, bool enable_backface_culling = false, bool enable_near_clipping = true, bool enable_far_clipping = true)
+{
+	// Some internal documentation / old todos or notes:
+	// maybe change and pass depthBuffer as an optional arg (&?), because usually we never need it outside the renderer. Or maybe even a getDepthBuffer().
+	// modelViewMatrix goes to eye-space (camera space), projection does ortho or perspective proj.
+	// bool enable_texturing = false; Maybe re-add later, not sure
+	// take a cv::Mat texture instead and convert to Texture internally? no, we don't want to recreate mipmap levels on each render() call.
+
+	assert(mesh.vertices.size() == mesh.colors.size() || mesh.colors.empty()); // The number of vertices has to be equal for both shape and colour, or, alternatively, it has to be a shape-only model.
+	assert(mesh.vertices.size() == mesh.texcoords.size() || mesh.texcoords.empty()); // same for the texcoords
+	// another assert: If cv::Mat texture != empty, then we need texcoords?
+
+	using cv::Mat;
+	using std::vector;
+
+	Mat colourbuffer = Mat::zeros(viewport_height, viewport_width, CV_8UC4); // make sure it's CV_8UC4?
+	Mat depthbuffer = std::numeric_limits<float>::max() * Mat::ones(viewport_height, viewport_width, CV_64FC1);
+
+	// Vertex shader:
+	//processedVertex = shade(Vertex); // processedVertex : pos, col, tex, texweight
+	// Assemble the vertices, project to clip space, and store as detail::Vertex (the internal representation):
+	vector<detail::Vertex> clipspace_vertices;
+	clipspace_vertices.reserve(mesh.vertices.size());
+	for (int i = 0; i < mesh.vertices.size(); ++i) { // "previously": mesh.vertex
+		Mat clipspace_coords = projection_matrix * model_view_matrix * Mat(mesh.vertices[i]);
+		cv::Vec3f vertex_colour;
+		if (mesh.colors.empty()) {
+			vertex_colour = cv::Vec3f(0.5f, 0.5f, 0.5f);
+		}
+		else {
+			vertex_colour = mesh.colors[i];
+		}
+		clipspace_vertices.push_back(detail::Vertex(clipspace_coords, vertex_colour, mesh.texcoords[i]));
+	}
+
+	// All vertices are in clip-space now.
+	// Prepare the rasterisation stage.
+	// For every vertex/tri:
+	vector<detail::TriangleToRasterize> triangles_to_raster;
+	for (const auto& tri_indices : mesh.tvi) {
+		// Todo: Split this whole stuff up. Make a "clip" function, ... rename "processProspective..".. what is "process"... get rid of "continue;"-stuff by moving stuff inside process...
+		// classify vertices visibility with respect to the planes of the view frustum
+		// we're in clip-coords (NDC), so just check if outside [-1, 1] x ...
+		// Actually we're in clip-coords and it's not the same as NDC. We're only in NDC after the division by w.
+		// We should do the clipping in clip-coords though. See http://www.songho.ca/opengl/gl_projectionmatrix.html for more details.
+		// However, when comparing against w_c below, we might run into the trouble of the sign again in the affine case.
+		// 'w' is always positive, as it is -z_camspace, and all z_camspace are negative.
+		unsigned char visibility_bits[3];
+		for (unsigned char k = 0; k < 3; k++)
+		{
+			visibility_bits[k] = 0;
+			float x_cc = clipspace_vertices[tri_indices[k]].position[0];
+			float y_cc = clipspace_vertices[tri_indices[k]].position[1];
+			float z_cc = clipspace_vertices[tri_indices[k]].position[2];
+			float w_cc = clipspace_vertices[tri_indices[k]].position[3];
+			if (x_cc < -w_cc)			// true if outside of view frustum. False if on or inside the plane.
+				visibility_bits[k] |= 1;	// set bit if outside of frustum
+			if (x_cc > w_cc)
+				visibility_bits[k] |= 2;
+			if (y_cc < -w_cc)
+				visibility_bits[k] |= 4;
+			if (y_cc > w_cc)
+				visibility_bits[k] |= 8;
+			if (enable_near_clipping && z_cc < -w_cc) // near plane frustum clipping
+				visibility_bits[k] |= 16;
+			if (enable_far_clipping && z_cc > w_cc) // far plane frustum clipping
+				visibility_bits[k] |= 32;
+		} // if all bits are 0, then it's inside the frustum
+		// all vertices are not visible - reject the triangle.
+		if ((visibility_bits[0] & visibility_bits[1] & visibility_bits[2]) > 0)
+		{
+			continue;
+		}
+		// all vertices are visible - pass the whole triangle to the rasterizer. = All bits of all 3 triangles are 0.
+		if ((visibility_bits[0] | visibility_bits[1] | visibility_bits[2]) == 0)
+		{
+			boost::optional<detail::TriangleToRasterize> t = detail::process_prospective_tri(clipspace_vertices[tri_indices[0]], clipspace_vertices[tri_indices[1]], clipspace_vertices[tri_indices[2]], viewport_width, viewport_height, enable_backface_culling);
+			if (t) {
+				triangles_to_raster.push_back(*t);
+			}
+			continue;
+		}
+		// at this moment the triangle is known to be intersecting one of the view frustum's planes
+		std::vector<detail::Vertex> vertices;
+		vertices.push_back(clipspace_vertices[tri_indices[0]]);
+		vertices.push_back(clipspace_vertices[tri_indices[1]]);
+		vertices.push_back(clipspace_vertices[tri_indices[2]]);
+		// split the triangle if it intersects the near plane:
+		if (enable_near_clipping)
+		{
+			vertices = detail::clip_polygon_to_plane_in_4d(vertices, cv::Vec4f(0.0f, 0.0f, -1.0f, -1.0f)); // "Normal" (or "4D hyperplane") of the near-plane. I tested it and it works like this but I'm a little bit unsure because Songho says the normal of the near-plane is (0,0,-1,1) (maybe I have to switch around the < 0 checks in the function?)
+		}
+
+		// triangulation of the polygon formed of vertices array
+		if (vertices.size() >= 3)
+		{
+			for (unsigned char k = 0; k < vertices.size() - 2; k++)
+			{
+				boost::optional<detail::TriangleToRasterize> t = detail::process_prospective_tri(vertices[0], vertices[1 + k], vertices[2 + k], viewport_width, viewport_height, enable_backface_culling);
+				if (t) {
+					triangles_to_raster.push_back(*t);
+				}
+			}
+		}
+	}
+
+	// Fragment/pixel shader: Colour the pixel values
+	// for every tri:
+	for (const auto& tri : triangles_to_raster) {
+		detail::raster_triangle(tri, colourbuffer, depthbuffer, texture, enable_far_clipping);
+	}
+	return std::make_pair(colourbuffer, depthbuffer);
+};
+
+	} /* namespace render */
+} /* namespace eos */
+
+#endif /* RENDER_HPP_ */

include/eos/render/render_affine.hpp

@@ -23,6 +23,7 @@
 #define RENDER_AFFINE_HPP_
 
 #include "eos/render/detail/render_detail.hpp"
+#include "eos/render/detail/render_affine_detail.hpp"
 #include "eos/render/Mesh.hpp"
 
 #include "opencv2/core/core.hpp"
@@ -50,7 +51,7 @@ namespace eos {
  */
 std::pair<cv::Mat, cv::Mat> render_affine(Mesh mesh, cv::Mat affine_camera_matrix, int viewport_width, int viewport_height, bool do_backface_culling = true)
 {
-	assert(mesh.vertices.size() == mesh.colors.size() || mesh.colors.empty());// The number of vertices has to be equal for both shape and colour, or, alternatively, it has to be a shape-only model.
+	assert(mesh.vertices.size() == mesh.colors.size() || mesh.colors.empty()); // The number of vertices has to be equal for both shape and colour, or, alternatively, it has to be a shape-only model.
 	//assert(mesh.vertices.size() == mesh.texcoords.size() || mesh.texcoords.empty()); // same for the texcoords
 
 	using cv::Mat;
@@ -108,7 +109,7 @@ std::pair<cv::Mat, cv::Mat> render_affine(Mesh mesh, cv::Mat affine_camera_matri
 
 	// Raster all triangles, i.e. colour the pixel values and write the z-buffer
 	for (auto&& triangle : triangles_to_raster) {
-		detail::raster_triangle(triangle, colourbuffer, depthbuffer);
+		detail::raster_triangle_affine(triangle, colourbuffer, depthbuffer);
 	}
 	return std::make_pair(colourbuffer, depthbuffer);
 };

include/eos/render/texture_extraction.hpp

@@ -32,6 +32,7 @@
 
 #include <tuple>
 #include <cassert>
+#include <future>
 
 namespace eos {
 	namespace render {
@@ -47,20 +48,31 @@ enum class TextureInterpolation {
 };
 
 // Just a forward declaration
-inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, cv::Mat depthbuffer, TextureInterpolation mapping_type, int isomap_resolution);
+inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, cv::Mat depthbuffer, bool compute_view_angle, TextureInterpolation mapping_type, int isomap_resolution);
 
 /**
  * Extracts the texture of the face from the given image
  * and stores it as isomap (a rectangular texture map).
  *
+ * Note/#Todo: Only use TextureInterpolation::NearestNeighbour
+ * for the moment, the other methods don't have correct handling of
+ * the alpha channel (and will most likely throw an exception).
+ *
+ * Todo: These should be renamed to extract_texture_affine? Can we combine both cases somehow?
+ * Or an overload with RenderingParameters?
+ *
+ * For TextureInterpolation::NearestNeighbour, returns a 4-channel isomap
+ * with the visibility in the 4th channel (0=invis, 255=visible).
+ *
  * @param[in] mesh A mesh with texture coordinates.
  * @param[in] affine_camera_matrix An estimated 3x4 affine camera matrix.
- * @param[in] image The image to extract the texture from.
+ * @param[in] image The image to extract the texture from. Should be 8UC3, other types not supported yet.
+ * @param[in] compute_view_angle A flag whether the view angle of each vertex should be computed and returned. If set to true, the angle will be encoded into the alpha channel (0 meaning occluded or facing away 90°, 127 meaning facing a 45° angle and 255 meaning front-facing, and all values in between). If set to false, the alpha channel will only contain 0 for occluded vertices and 255 for visible vertices.
  * @param[in] mapping_type The interpolation type to be used for the extraction.
  * @param[in] isomap_resolution The resolution of the generated isomap. Defaults to 512x512.
  * @return The extracted texture as isomap (texture map).
  */
-inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, TextureInterpolation mapping_type, int isomap_resolution = 512)
+inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, bool compute_view_angle = false, TextureInterpolation mapping_type = TextureInterpolation::NearestNeighbour, int isomap_resolution = 512)
 {
 	// Render the model to get a depth buffer:
 	cv::Mat depthbuffer;
@@ -68,7 +80,7 @@ inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat
 	// Note: There's potential for optimisation here - we don't need to do everything that is done in render_affine to just get the depthbuffer.
 
 	// Now forward the call to the actual texture extraction function:
-	return extract_texture(mesh, affine_camera_matrix, image, depthbuffer, mapping_type, isomap_resolution);
+	return extract_texture(mesh, affine_camera_matrix, image, depthbuffer, compute_view_angle, mapping_type, isomap_resolution);
 };
 
 /**
@@ -78,22 +90,29 @@ inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat
  * To just run the texture extraction, see the overload
  * extract_texture(Mesh, cv::Mat, cv::Mat, TextureInterpolation, int).
  *
+ * It might be wise to remove this overload as it can get quite confusing
+ * with the zbuffer. Obviously the depthbuffer given should have been created
+ * with the same (affine or ortho) projection matrix than the texture extraction is called with.
+ *
  * @param[in] mesh A mesh with texture coordinates.
  * @param[in] affine_camera_matrix An estimated 3x4 affine camera matrix.
  * @param[in] image The image to extract the texture from.
  * @param[in] depthbuffer A pre-calculated depthbuffer image.
+ * @param[in] compute_view_angle A flag whether the view angle of each vertex should be computed and returned. If set to true, the angle will be encoded into the alpha channel (0 meaning occluded or facing away 90°, 127 meaning facing a 45° angle and 255 meaning front-facing, and all values in between). If set to false, the alpha channel will only contain 0 for occluded vertices and 255 for visible vertices.
  * @param[in] mapping_type The interpolation type to be used for the extraction.
  * @param[in] isomap_resolution The resolution of the generated isomap. Defaults to 512x512.
  * @return The extracted texture as isomap (texture map).
  */
-inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, cv::Mat depthbuffer, TextureInterpolation mapping_type, int isomap_resolution = 512)
+inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, cv::Mat depthbuffer, bool compute_view_angle = false, TextureInterpolation mapping_type = TextureInterpolation::NearestNeighbour, int isomap_resolution = 512)
 {
 	assert(mesh.vertices.size() == mesh.texcoords.size());
+	assert(image.type() == CV_8UC3); // the other cases are not yet supported
 
 	using cv::Mat;
 	using cv::Vec2f;
 	using cv::Vec3f;
 	using cv::Vec4f;
+	using cv::Vec3b;
 	using std::min;
 	using std::max;
 	using std::floor;
@@ -101,10 +120,15 @@ inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat
 
 	affine_camera_matrix = detail::calculate_affine_z_direction(affine_camera_matrix);
 
-	Mat isomap = Mat::zeros(isomap_resolution, isomap_resolution, CV_8UC3); // #Todo: We do want an alpha channel. Will be added soon-ish.
-	// #Todo: We should handle gray images, but output a 3-channel isomap nevertheless I think.
+	Mat isomap = Mat::zeros(isomap_resolution, isomap_resolution, CV_8UC4);
+	// #Todo: We should handle gray images, but output a 4-channel isomap nevertheless I think.
 
+	std::vector<std::future<void>> results;
 	for (const auto& triangle_indices : mesh.tvi) {
+
+		// Note: If there's a performance problem, there's no need to capture the whole mesh - we could capture only the three required vertices with their texcoords.
+		auto extract_triangle = [&mesh, &affine_camera_matrix, &triangle_indices, &depthbuffer, &isomap, &mapping_type, &image, &compute_view_angle]() {
+
 			// Find out if the current triangle is visible:
 			// We do a second rendering-pass here. We use the depth-buffer of the final image, and then, here,
 			// check if each pixel in a triangle is visible. If the whole triangle is visible, we use it to extract
@@ -116,13 +140,49 @@ inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat
 			// - We transform them later (below) a second time. Only do it once.
 
 			// Project the triangle vertices to screen coordinates, and use the depthbuffer to check whether the triangle is visible:
-		Vec4f v0 = Mat(affine_camera_matrix * Mat(mesh.vertices[triangle_indices[0]]));
-		Vec4f v1 = Mat(affine_camera_matrix * Mat(mesh.vertices[triangle_indices[1]]));
-		Vec4f v2 = Mat(affine_camera_matrix * Mat(mesh.vertices[triangle_indices[2]]));
+			const Vec4f v0 = Mat(affine_camera_matrix * Mat(mesh.vertices[triangle_indices[0]]));
+			const Vec4f v1 = Mat(affine_camera_matrix * Mat(mesh.vertices[triangle_indices[1]]));
+			const Vec4f v2 = Mat(affine_camera_matrix * Mat(mesh.vertices[triangle_indices[2]]));
 
 			if (!detail::is_triangle_visible(v0, v1, v2, depthbuffer))
 			{
-			continue;
+				//continue;
+				return;
+			}
+
+			float alpha_value;
+			if (compute_view_angle)
+			{
+				// Calculate how well visible the current triangle is:
+				// (in essence, the dot product of the viewing direction (0, 0, 1) and the face normal)
+				const Vec3f face_normal = calculate_face_normal(Vec3f(Mat(mesh.vertices[triangle_indices[0]]).rowRange(0, 3)), Vec3f(Mat(mesh.vertices[triangle_indices[1]]).rowRange(0, 3)), Vec3f(Mat(mesh.vertices[triangle_indices[2]]).rowRange(0, 3)));
+				// Transform the normal to "screen" (kind of "eye") space using the upper 3x3 part of the affine camera matrix (=the translation can be ignored):
+				Vec3f face_normal_transformed = Mat(affine_camera_matrix.rowRange(0, 3).colRange(0, 3) * Mat(face_normal));
+				face_normal_transformed /= cv::norm(face_normal_transformed, cv::NORM_L2); // normalise to unit length
+				// Implementation notes regarding the affine camera matrix and the sign:
+				// If the matrix given were the model_view matrix, the sign would be correct.
+				// However, affine_camera_matrix includes glm::ortho, which includes a z-flip.
+				// So we need to flip one of the two signs.
+				// * viewing_direction(0.0f, 0.0f, 1.0f) is correct if affine_camera_matrix were only a model_view matrix
+				// * affine_camera_matrix includes glm::ortho, which flips z, so we flip the sign of viewing_direction.
+				// We don't need the dot product since viewing_direction.xy are 0 and .z is 1:
+				const float angle = -face_normal_transformed[2]; // flip sign, see above
+				assert(angle >= -1.f && angle <= 1.f);
+				// angle is [-1, 1].
+				//  * +1 means   0° (same direction)
+				//  *  0 means  90°
+				//  * -1 means 180° (facing opposite directions)
+				// It's a linear relation, so +0.5 is 45° etc.
+				// An angle larger than 90° means the vertex won't be rendered anyway (because it's back-facing) so we encode 0° to 90°.
+				if (angle < 0.0f) {
+					alpha_value = 0.0f;
+				} else {
+					alpha_value = angle * 255.0f;
+				}
+			}
+			else {
+				// no visibility angle computation - if the triangle/pixel is visible, set the alpha chan to 255 (fully visible pixel).
+				alpha_value = 255.0f;
 			}
 
 			// Todo: Documentation
@@ -153,13 +213,13 @@ inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat
 			for (int x = min(dst_tri[0].x, min(dst_tri[1].x, dst_tri[2].x)); x < max(dst_tri[0].x, max(dst_tri[1].x, dst_tri[2].x)); ++x) {
 				for (int y = min(dst_tri[0].y, min(dst_tri[1].y, dst_tri[2].y)); y < max(dst_tri[0].y, max(dst_tri[1].y, dst_tri[2].y)); ++y) {
 					if (detail::is_point_in_triangle(cv::Point2f(x, y), dst_tri[0], dst_tri[1], dst_tri[2])) {
-					if (mapping_type == TextureInterpolation::Area){
+						if (mapping_type == TextureInterpolation::Area) {
 
 							// calculate positions of 4 corners of pixel in image (src)
-						Vec3f homogenous_dst_upper_left(x - 0.5, y - 0.5, 1.f);
-						Vec3f homogenous_dst_upper_right(x + 0.5, y - 0.5, 1.f);
-						Vec3f homogenous_dst_lower_left(x - 0.5, y + 0.5, 1.f);
-						Vec3f homogenous_dst_lower_right(x + 0.5, y + 0.5, 1.f);
+							Vec3f homogenous_dst_upper_left(x - 0.5f, y - 0.5f, 1.0f);
+							Vec3f homogenous_dst_upper_right(x + 0.5f, y - 0.5f, 1.0f);
+							Vec3f homogenous_dst_lower_left(x - 0.5f, y + 0.5f, 1.0f);
+							Vec3f homogenous_dst_lower_right(x + 0.5f, y + 0.5f, 1.0f);
 
 							Vec2f src_texel_upper_left = Mat(warp_mat_org_inv * Mat(homogenous_dst_upper_left));
 							Vec2f src_texel_upper_right = Mat(warp_mat_org_inv * Mat(homogenous_dst_upper_right));
@@ -181,35 +241,37 @@ inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat
 									if (detail::is_point_in_triangle(cv::Point2f(a, b), src_texel_upper_left, src_texel_lower_left, src_texel_upper_right) || detail::is_point_in_triangle(cv::Point2f(a, b), src_texel_lower_left, src_texel_upper_right, src_texel_lower_right)) {
 										if (a < image.cols && b < image.rows) { // if src_texel in triangle and in image
 											num_texels++;
-										color += image.at<cv::Vec3b>(b, a);
+											color += image.at<Vec3b>(b, a);
 										}
 									}
 								}
 							}
 							if (num_texels > 0)
 								color = color / num_texels;
-						else { // if no corresponding texel found, nearest neighbor interpolation
+							else { // if no corresponding texel found, nearest neighbour interpolation
 								// calculate corresponding position of dst_coord pixel center in image (src)
-							Vec3f homogenous_dst_coord = Vec3f(x, y, 1.f);
+								Vec3f homogenous_dst_coord = Vec3f(x, y, 1.0f);
 								Vec2f src_texel = Mat(warp_mat_org_inv * Mat(homogenous_dst_coord));
 
 								if ((cvRound(src_texel[1]) < image.rows) && cvRound(src_texel[0]) < image.cols) {
-								color = image.at<cv::Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]));
+									color = image.at<Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]));
 								}
 							}
-						isomap.at<cv::Vec3b>(y, x) = color;
+							isomap.at<Vec3b>(y, x) = color;
 						}
 						else if (mapping_type == TextureInterpolation::Bilinear) {
 
 							// calculate corresponding position of dst_coord pixel center in image (src)
-						Vec3f homogenous_dst_coord(x, y, 1.f);
+							Vec3f homogenous_dst_coord(x, y, 1.0f);
 							Vec2f src_texel = Mat(warp_mat_org_inv * Mat(homogenous_dst_coord));
 
 							// calculate distances to next 4 pixels
-						float distance_upper_left = sqrt(powf(src_texel[0] - floor(src_texel[0]), 2) + powf(src_texel[1] - floor(src_texel[1]), 2));
-						float distance_upper_right = sqrt(powf(src_texel[0] - floor(src_texel[0]), 2) + powf(src_texel[1] - ceil(src_texel[1]), 2));
-						float distance_lower_left = sqrt(powf(src_texel[0] - ceil(src_texel[0]), 2) + powf(src_texel[1] - floor(src_texel[1]), 2));
-						float distance_lower_right = sqrt(powf(src_texel[0] - ceil(src_texel[0]), 2) + powf(src_texel[1] - ceil(src_texel[1]), 2));
+							using std::sqrt;
+							using std::pow;
+							float distance_upper_left = sqrt(pow(src_texel[0] - floor(src_texel[0]), 2) + pow(src_texel[1] - floor(src_texel[1]), 2));
+							float distance_upper_right = sqrt(pow(src_texel[0] - floor(src_texel[0]), 2) + pow(src_texel[1] - ceil(src_texel[1]), 2));
+							float distance_lower_left = sqrt(pow(src_texel[0] - ceil(src_texel[0]), 2) + pow(src_texel[1] - floor(src_texel[1]), 2));
+							float distance_lower_right = sqrt(pow(src_texel[0] - ceil(src_texel[0]), 2) + pow(src_texel[1] - ceil(src_texel[1]), 2));
 
 							// normalise distances
 							float sum_distances = distance_lower_left + distance_lower_right + distance_upper_left + distance_upper_right;
@@ -219,27 +281,39 @@ inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat
 							distance_upper_right /= sum_distances;
 
 							// set color depending on distance from next 4 pixels
-						for (int color = 0; color < 3; color++){
-							float color_upper_left = image.at<cv::Vec3b>(floor(src_texel[1]), floor(src_texel[0]))[color] * distance_upper_left;
-							float color_upper_right = image.at<cv::Vec3b>(floor(src_texel[1]), ceil(src_texel[0]))[color] * distance_upper_right;
-							float color_lower_left = image.at<cv::Vec3b>(ceil(src_texel[1]), floor(src_texel[0]))[color] * distance_lower_left;
-							float color_lower_right = image.at<cv::Vec3b>(ceil(src_texel[1]), ceil(src_texel[0]))[color] * distance_lower_right;
+							for (int color = 0; color < 3; ++color) {
+								float color_upper_left = image.at<Vec3b>(floor(src_texel[1]), floor(src_texel[0]))[color] * distance_upper_left;
+								float color_upper_right = image.at<Vec3b>(floor(src_texel[1]), ceil(src_texel[0]))[color] * distance_upper_right;
+								float color_lower_left = image.at<Vec3b>(ceil(src_texel[1]), floor(src_texel[0]))[color] * distance_lower_left;
+								float color_lower_right = image.at<Vec3b>(ceil(src_texel[1]), ceil(src_texel[0]))[color] * distance_lower_right;
 
-							isomap.at<cv::Vec3b>(y, x)[color] = color_upper_left + color_upper_right + color_lower_left + color_lower_right;
+								isomap.at<Vec3b>(y, x)[color] = color_upper_left + color_upper_right + color_lower_left + color_lower_right;
 							}
 						}
 						else if (mapping_type == TextureInterpolation::NearestNeighbour) {
 
 							// calculate corresponding position of dst_coord pixel center in image (src)
-						Vec3f homogenous_dst_coord = Vec3f(x, y, 1.f);
-						Vec2f src_texel = Mat(warp_mat_org_inv * Mat(homogenous_dst_coord));
+							const Mat homogenous_dst_coord(Vec3f(x, y, 1.0f));
+							const Vec2f src_texel = Mat(warp_mat_org_inv * homogenous_dst_coord);
 
-						if ((cvRound(src_texel[1]) < image.rows) && (cvRound(src_texel[0]) < image.cols))
-							isomap.at<cv::Vec3b>(y, x) = image.at<cv::Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]));
+							if ((cvRound(src_texel[1]) < image.rows) && (cvRound(src_texel[0]) < image.cols) && cvRound(src_texel[0]) > 0 && cvRound(src_texel[1]) > 0)
+							{
+								cv::Vec4b isomap_pixel;
+								isomap.at<cv::Vec4b>(y, x)[0] = image.at<Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]))[0];
+								isomap.at<cv::Vec4b>(y, x)[1] = image.at<Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]))[1];
+								isomap.at<cv::Vec4b>(y, x)[2] = image.at<Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]))[2];
+								isomap.at<cv::Vec4b>(y, x)[3] = static_cast<uchar>(alpha_value); // pixel is visible
+							}
 						}
 					}
 				}
 			}
+		}; // end lambda auto extract_triangle();
+		results.emplace_back(std::async(extract_triangle));
+	} // end for all mesh.tvi
+	// Collect all the launched tasks:
+	for (auto&& r : results) {
+		r.get();
 	}
 	return isomap;
 };

include/eos/render/utils.hpp

@@ -22,7 +22,10 @@
 #ifndef RENDER_UTILS_HPP_
 #define RENDER_UTILS_HPP_
 
+#include "eos/render/Mesh.hpp"
+
 #include "opencv2/core/core.hpp"
+#include "opencv2/imgproc/imgproc.hpp"
 
 namespace eos {
 	namespace render {
@@ -78,6 +81,139 @@ inline cv::Vec2f screen_to_clip_space(const cv::Vec2f& screen_coordinates, int s
 	return cv::Vec2f(x_cs, y_cs);
 };
 
+/**
+ * Calculates the normal of a face (or triangle), i.e. the
+ * per-face normal. Return normal will be normalised.
+ * Assumes the triangle is given in CCW order, i.e. vertices
+ * in counterclockwise order on the screen are front-facing.
+ *
+ * @param[in] v0 First vertex.
+ * @param[in] v1 Second vertex.
+ * @param[in] v2 Third vertex.
+ * @return The unit-length normal of the given triangle.
+ */
+cv::Vec3f calculate_face_normal(const cv::Vec3f& v0, const cv::Vec3f& v1, const cv::Vec3f& v2)
+{
+	cv::Vec3f n = (v1 - v0).cross(v2 - v0); // v0-to-v1 x v0-to-v2
+	n /= cv::norm(n);
+	return n;
+};
+
+/**
+ * Draws the texture coordinates (uv-coords) of the given mesh
+ * into an image by looping over the triangles and drawing each
+ * triangle's texcoords.
+ *
+ * @param[in] mesh A mesh with texture coordinates.
+ * @param[in] image An optional image to draw onto.
+ * @return An image with the texture coordinate triangles drawn in it, 512x512 if no image is given.
+ */
+cv::Mat draw_texcoords(Mesh mesh, cv::Mat image = cv::Mat())
+{
+	using cv::Point2f;
+	using cv::Scalar;
+	if (image.empty())
+	{
+		image = cv::Mat(512, 512, CV_8UC4, Scalar(0.0f, 0.0f, 0.0f, 255.0f));
+	}
+
+	for (const auto& triIdx : mesh.tvi) {
+		cv::line(image, Point2f(mesh.texcoords[triIdx[0]][0] * image.cols, mesh.texcoords[triIdx[0]][1] * image.rows), Point2f(mesh.texcoords[triIdx[1]][0] * image.cols, mesh.texcoords[triIdx[1]][1] * image.rows), Scalar(255.0f, 0.0f, 0.0f));
+		cv::line(image, Point2f(mesh.texcoords[triIdx[1]][0] * image.cols, mesh.texcoords[triIdx[1]][1] * image.rows), Point2f(mesh.texcoords[triIdx[2]][0] * image.cols, mesh.texcoords[triIdx[2]][1] * image.rows), Scalar(255.0f, 0.0f, 0.0f));
+		cv::line(image, Point2f(mesh.texcoords[triIdx[2]][0] * image.cols, mesh.texcoords[triIdx[2]][1] * image.rows), Point2f(mesh.texcoords[triIdx[0]][0] * image.cols, mesh.texcoords[triIdx[0]][1] * image.rows), Scalar(255.0f, 0.0f, 0.0f));
+	}
+	return image;
+};
+
+// TODO: Should go to detail:: namespace, or texturing/utils or whatever.
+unsigned int get_max_possible_mipmaps_num(unsigned int width, unsigned int height)
+{
+	unsigned int mipmapsNum = 1;
+	unsigned int size = std::max(width, height);
+
+	if (size == 1)
+		return 1;
+
+	do {
+		size >>= 1;
+		mipmapsNum++;
+	} while (size != 1);
+
+	return mipmapsNum;
+};
+
+inline bool is_power_of_two(int x)
+{
+	return !(x & (x - 1));
+};
+
+class Texture
+{
+public:
+	// Todo: This whole class needs a major overhaul and documentation.
+
+	std::vector<cv::Mat> mipmaps;	// make Texture a friend class of renderer, then move this to private?
+	unsigned char widthLog, heightLog; // log2 of width and height of the base mip-level
+
+//private:
+	//std::string filename;
+	unsigned int mipmaps_num;
+};
+
+// throws: ocv exc,  runtime_ex
+Texture create_mipmapped_texture(cv::Mat image, unsigned int mipmapsNum = 0) {
+	assert(image.type() == CV_8UC3 || image.type() == CV_8UC4);
+
+	Texture texture;
+
+	texture.mipmaps_num = (mipmapsNum == 0 ? get_max_possible_mipmaps_num(image.cols, image.rows) : mipmapsNum);
+	/*if (mipmapsNum == 0)
+	{
+	uchar mmn = render::utils::MatrixUtils::getMaxPossibleMipmapsNum(image.cols, image.rows);
+	this->mipmapsNum = mmn;
+	} else
+	{
+	this->mipmapsNum = mipmapsNum;
+	}*/
+
+	if (texture.mipmaps_num > 1)
+	{
+		if (!is_power_of_two(image.cols) || !is_power_of_two(image.rows))
+		{
+			throw std::runtime_error("Error: Couldn't generate mipmaps, width or height not power of two.");
+		}
+	}
+	if (image.type() == CV_8UC3)
+	{
+		image.convertTo(image, CV_8UC4); // Most often, the input img is CV_8UC3. Img is BGR. Add an alpha channel
+		cv::cvtColor(image, image, CV_BGR2BGRA);
+	}
+
+	int currWidth = image.cols;
+	int currHeight = image.rows;
+	std::vector<cv::Mat> mipmaps;
+	for (int i = 0; i < texture.mipmaps_num; i++)
+	{
+		if (i == 0) {
+			mipmaps.push_back(image);
+		}
+		else {
+			cv::Mat currMipMap(currHeight, currWidth, CV_8UC4);
+			cv::resize(mipmaps[i - 1], currMipMap, currMipMap.size());
+			mipmaps.push_back(currMipMap);
+		}
+
+		if (currWidth > 1)
+			currWidth >>= 1;
+		if (currHeight > 1)
+			currHeight >>= 1;
+	}
+	texture.mipmaps = mipmaps;
+	texture.widthLog = (uchar)(std::log(mipmaps[0].cols) / CV_LOG2 + 0.0001f); // std::epsilon or something? or why 0.0001f here?
+	texture.heightLog = (uchar)(std::log(mipmaps[0].rows) / CV_LOG2 + 0.0001f); // Changed std::logf to std::log because it doesnt compile in linux (gcc 4.8). CHECK THAT
+	return texture;
+};
+
 	} /* namespace render */
 } /* namespace eos */