Added visibility-testing by ray-casting to edge-fitting

include/eos/fitting/closest_edge_fitting.hpp

@@ -33,25 +33,84 @@
 #include "glm/vec4.hpp"
 #include "glm/mat4x4.hpp"
 
+#include "boost/optional.hpp"
+
 #include <vector>
 #include <algorithm>
+#include <utility>
 
 namespace eos {
 	namespace fitting {
 
+/**
+ * @brief Computes the intersection of the given ray with the given triangle.
+ *
+ * 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.
+ * When culling is on, rays intersecting triangles from the back will be discarded -
+ * otherwise, the triangles normal direction w.r.t. the ray direction is just ignored.
+ *
+ * Note: The use of optional might turn out as a performance problem, as this
+ * function is called loads of time - how costly is it to construct a boost::none optional?
+ *
+ * @param[in] ray_origin Ray origin.
+ * @param[in] ray_direction Ray direction.
+ * @param[in] v0 First vertex of a triangle.
+ * @param[in] v1 Second vertex of a triangle.
+ * @param[in] v2 Third vertex of a triangle.
+ * @param[in] enable_backculling When culling is on, rays intersecting triangles from the back will be discarded.
+ * @return Whether the ray intersects the triangle, and if yes, including the distance.
+ */
+std::pair<bool, boost::optional<float>> ray_triangle_intersect(const glm::vec3& ray_origin, const glm::vec3& ray_direction, const glm::vec3& v0, const glm::vec3& v1, const glm::vec3& v2, bool enable_backculling)
+{
+	using glm::vec3;
+	const float epsilon = 1e-6f;
+
+	vec3 v0v1 = v1 - v0;
+	vec3 v0v2 = v2 - v0;
+
+	vec3 pvec = glm::cross(ray_direction, v0v2);
+
+	float det = glm::dot(v0v1, pvec);
+	if (enable_backculling)
+	{
+		// If det is negative, the triangle is back-facing.
+		// If det is close to 0, the ray misses the triangle.
+		if (det < epsilon)
+			return { false, boost::none };
+	}
+	else {
+		// If det is close to 0, the ray and triangle are parallel.
+		if (std::abs(det) < epsilon)
+			return { false, boost::none };
+	}
+	float inv_det = 1 / det;
+
+	vec3 tvec = ray_origin - v0;
+	auto u = glm::dot(tvec, pvec) * inv_det;
+	if (u < 0 || u > 1)
+		return { false, boost::none };
+
+	vec3 qvec = glm::cross(tvec, v0v1);
+	auto v = glm::dot(ray_direction, qvec) * inv_det;
+	if (v < 0 || u + v > 1)
+		return { false, boost::none };
+
+	auto t = glm::dot(v0v2, qvec) * inv_det;
+
+	return { true, t };
+};
+
 /**
  * @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 returns a list of the (unique) vertices that make the boundary edges.
+ * 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.
  *
- * Notes:
- *   Actually the function only needs the vertices part of the Mesh (i.e. a
- *   std::vector<glm::vec4>).
- *   But it could use of the face normals, if they're already computed. But our
- *   Meshes don't contain normals anyway right now.
- *
  * @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.
@@ -59,8 +118,6 @@ namespace eos {
  */
 std::vector<int> occluding_boundary_vertices(const eos::render::Mesh& mesh, const morphablemodel::EdgeTopology& edge_topology, glm::mat4x4 R)
 {
-	std::vector<int> occluding_vertices; // The model's contour vertices
-
 	// 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); });
@@ -94,6 +151,7 @@ std::vector<int> occluding_boundary_vertices(const eos::render::Mesh& mesh, cons
 	}
 	// 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!
@@ -104,10 +162,48 @@ std::vector<int> occluding_boundary_vertices(const eos::render::Mesh& mesh, cons
 	std::sort(begin(occluding_vertices), end(occluding_vertices));
 	occluding_vertices.erase(std::unique(begin(occluding_vertices), end(occluding_vertices)), end(occluding_vertices));
 
-	// % Remove vertices from occluding boundary list that are not visible:
-	// Todo! Not doing yet. Have to render or ray-cast.
+	// Perform ray-casting to find out which vertices are not visible (i.e. self-occluded):
+	std::vector<bool> visibility;
+	for (const auto& vertex_idx : occluding_vertices)
+	{
+		bool visible = true;
+		// For every tri of the rotated mesh:
+		for (auto&& tri : mesh.tvi)
+		{
+			auto& v0 = rotated_vertices[tri[0]];
+			auto& v1 = rotated_vertices[tri[1]];
+			auto& v2 = rotated_vertices[tri[2]];
+
+			glm::vec3 ray_origin = rotated_vertices[vertex_idx];
+			glm::vec3 ray_direction(0.0f, 0.0f, 1.0f); // we shoot the ray from the vertex towards the camera
+			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);
+	}
 
-	return occluding_vertices;
+	// 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;
 };
 
 /** A simple vector-of-vectors adaptor for nanoflann, without duplicating the storage.