u02: Implement barycentric coordinates

u02/include/util.h

@@ -30,3 +30,9 @@ std::pair<int, int> transform_inv(Transformation transformation, int x, int y);
 // to make the given endpoints of a line conform
 // to the standard case for rasterization.
 Transformation transformation_to_standard_case(int x0, int y0, int x1, int y1);
+
+// Returns the barycentric coordinates of the point given by (xp, yp)
+// in the triangle given by the three points (x0, y0), (x1, y1), (x2, y2).
+std::tuple<float, float, float> barycentric_coordinates(int x0, int y0, int x1,
+                                                        int y1, int x2, int y2,
+                                                        int xp, int yp);

u02/src/tests.cpp

@@ -4,6 +4,7 @@
 #include <catch2/generators/catch_generators.hpp>
 #include <catch2/generators/catch_generators_adapters.hpp>
 #include <catch2/generators/catch_generators_random.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
 
 #include "bresenham_circle_tool.h"
 #include "bresenham_line_tool.h"
@@ -14,6 +15,8 @@
 #include "recursive_fill_tool.h"
 #include "util.h"
 
+using Catch::Matchers::WithinRel;
+
 TEST_CASE("Transform Mirror") {
   // elementary operations
   REQUIRE(transform(TRANSFORM_MIRROR_X, 10, 20) == std::make_pair(10, -20));
@@ -370,3 +373,57 @@ TEST_CASE("Bresenham circle (prop: √(x²+y²)-r<ε)") {
   }
   REQUIRE(pass);
 }
+
+TEST_CASE("Barycentric coordinates: Edge cases") {
+  int x0 = 0, y0 = 0, x1 = 0, y1 = 10, x2 = 10, y2 = 0;
+  float b1, b2, b3;
+
+  // point on vertex
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 0, 0);
+  REQUIRE(b1 == 1);
+  REQUIRE(b2 == 0);
+  REQUIRE(b3 == 0);
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 0, 10);
+  REQUIRE(b1 == 0);
+  REQUIRE(b2 == 1);
+  REQUIRE(b3 == 0);
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 10, 0);
+  REQUIRE(b1 == 0);
+  REQUIRE(b2 == 0);
+  REQUIRE(b3 == 1);
+
+  // point on edge
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 0, 5);
+  REQUIRE(b1 == 0.5);
+  REQUIRE(b2 == 0.5);
+  REQUIRE(b3 == 0);
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 5, 0);
+  REQUIRE(b1 == 0.5);
+  REQUIRE(b2 == 0);
+  REQUIRE(b3 == 0.5);
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, 5, 5);
+  REQUIRE(b1 == 0);
+  REQUIRE(b2 == 0.5);
+  REQUIRE(b3 == 0.5);
+}
+
+TEST_CASE("Barycentric coordinates (prop: Σ = 1)") {
+  int x0 = GENERATE(take(2, random(-100, 100)));
+  int y0 = GENERATE(take(2, random(-100, 100)));
+  int x1 = GENERATE(take(2, random(-100, 100)));
+  int y1 = GENERATE(take(2, random(-100, 100)));
+  int x2 = GENERATE(take(2, random(-100, 100)));
+  int y2 = GENERATE(take(2, random(-100, 100)));
+
+  int x = GENERATE(take(5, random(-100, 100)));
+  int y = GENERATE(take(5, random(-100, 100)));
+  float b1, b2, b3;
+
+  std::tie(b1, b2, b3) = barycentric_coordinates(x0, y0, x1, y1, x2, y2, x, y);
+
+  REQUIRE_THAT(b1 + b2 + b3, WithinRel(1.0, 1e-4));
+}

u02/src/util.cpp

@@ -139,3 +139,23 @@ Transformation transformation_to_standard_case(int x0, int y0, int x1, int y1) {
 
   return transformation;
 }
+
+std::tuple<float, float, float> barycentric_coordinates(int x0, int y0, int x1,
+                                                        int y1, int x2, int y2,
+                                                        int xp, int yp) {
+  // Source:
+  // https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Vertex_approach
+  float b1 = x1 * y2 - x2 * y1 + xp * (y1 - y2) + yp * (x2 - x1);
+  float b2 = x2 * y0 - x0 * y2 + xp * (y2 - y0) + yp * (x0 - x2);
+  float b3 = x0 * y1 - x1 * y0 + xp * (y0 - y1) + yp * (x1 - x0);
+
+  // reciprocal computed directly for performance
+  float area_factor =
+      1 / static_cast<float>(x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1));
+
+  b1 *= area_factor;
+  b2 *= area_factor;
+  b3 *= area_factor;
+
+  return {b1, b2, b3};
+}