Implement ray-triangle-intersection; separate Vector and Point for this

Also, the Ray got an unimplemented `intersects_triangle` method, so the kdtree intersection function should be complete.

geometry.h

@@ -0,0 +1,109 @@
+// Forward declarations
+struct Point;
+
+struct Node {
+    Node(int axis, Point *point, Node *left, Node *right)
+        : axis(axis), point(point), left(left), right(right) {}
+
+    int axis;
+
+    Point *point;
+
+    Node *left;
+    Node *right;
+};
+
+struct Triangle {
+    Triangle(Point *p1, Point *p2, Point *p3) : p1(p1), p2(p2), p3(p3) {}
+
+    Point *p1;
+    Point *p2;
+    Point *p3;
+};
+
+struct Vector {
+    Vector(float coordinates[3]) : c(coordinates) {}
+
+    // Avoid having to write vector.c[index], instead allow vector[index]
+    float operator[](int i) const { return c[i]; }
+    float &operator[](int i) { return c[i]; }
+
+    Vector operator+(const Vector &other) const {
+        return Vector(new float[3]{c[0] + other.c[0], c[1] + other.c[1], c[2] + other.c[2]});
+    }
+
+    Vector operator-(const Vector &other) const {
+        return Vector(new float[3]{c[0] - other.c[0], c[1] - other.c[1], c[2] - other.c[2]});
+    }
+
+    Vector operator*(float scalar) const {
+        return Vector(new float[3]{c[0] * scalar, c[1] * scalar, c[2] * scalar});
+    }
+
+    Vector cross(const Vector &other) {
+        // TODO
+        return other;
+    }
+
+    float dot(const Vector &other) {
+        // TODO
+        return 0.0;
+    }
+
+    float *c;
+};
+
+struct Point {
+    Point(float coordinates[3], Triangle *triangle)
+        : pos(Vector(coordinates)), triangle(triangle) {}
+
+    Vector pos;
+
+    Triangle *triangle;
+};
+
+struct Ray {
+    Ray(Vector origin, Vector direction) : origin(origin), direction(direction) {}
+
+    Vector origin;
+
+    Vector direction;
+
+    bool intersects_triangle(Triangle *triangle) {
+        // Ray-triangle-intersection with the Möller–Trumbore algorithm
+        // https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
+        const float EPSILON = 0.0000001;
+
+        Vector p1 = triangle->p1->pos;
+        Vector p2 = triangle->p2->pos;
+        Vector p3 = triangle->p3->pos;
+
+        Vector edge1 = p2 - p1;
+        Vector edge2 = p3 - p1;
+
+        Vector h = direction.cross(edge2);
+        float a = edge1.dot(h);
+
+        if (a > -EPSILON && a < EPSILON) return false; // This ray is parallel to this triangle.
+
+        float f = 1.0 / a;
+        Vector s = origin - p1;
+        float u = f * s.dot(h);
+
+        if (u < 0.0 || u > 1.0) return false;
+
+        Vector q = s.cross(edge1);
+        float v = f * direction.dot(q);
+        if (v < 0.0 || u + v > 1.0) return false;
+
+        // At this stage we can compute t to find out where the intersection point is on the
+        // line.
+        float t = f * edge2.dot(q);
+        if (t > EPSILON) {
+            return true;
+        } else {
+            // This means that there is a line intersection but not a ray intersection.
+            return false;
+        }
+    }
+};

kdtree.h

@@ -1,64 +1,10 @@
+#include "geometry.h"
 #include <algorithm>
 #include <glm/glm.hpp>
 #include <iostream>
 #include <string>
 #include <vector>
 
-// Forward declarations
-struct Node;
-struct Point;
-struct Triangle;
-
-struct Node {
-    Node(int axis, Point *point, Node *left, Node *right)
-        : axis(axis), point(point), left(left), right(right) {}
-
-    int axis;
-
-    Point *point;
-
-    Node *left;
-    Node *right;
-};
-
-struct Triangle {
-    Triangle(Point *p1, Point *p2, Point *p3) : p1(p1), p2(p2), p3(p3) {}
-
-    Point *p1;
-    Point *p2;
-    Point *p3;
-};
-
-struct Point {
-    Point(float coordinates[3], Triangle *triangle)
-        : coordinates(coordinates), triangle(triangle) {}
-
-    Point operator+(const Point &other) const {
-        return Point(new float[3]{coordinates[0] + other.coordinates[0],
-                                  coordinates[1] + other.coordinates[1],
-                                  coordinates[2] + other.coordinates[2]},
-                     nullptr);
-    }
-
-    Point operator*(float scalar) const {
-        return Point(
-            new float[3]{coordinates[0] * scalar, coordinates[1] * scalar, coordinates[2] * scalar},
-            nullptr);
-    }
-
-    float *coordinates;
-
-    Triangle *triangle;
-};
-
-struct Ray {
-    Ray(float origin[3], float direction[3]) : origin(origin), direction(direction) {}
-
-    float *origin;
-
-    float *direction;
-};
-
 class KDTree {
   public:
     KDTree(std::vector<Point *> points) { root = build(points, 0); }
@@ -82,7 +28,7 @@ class KDTree {
     // greater coordinate in the given axis.
     auto get_point_comparator(int axis) {
         return [axis](Point *p1, Point *p2) {
-            return p1->coordinates[axis] < p2->coordinates[axis];
+            return p1->pos[axis] < p2->pos[axis];
         };
     }
 
@@ -101,7 +47,7 @@ class KDTree {
             Point *min = *std::min_element(points.begin(), points.end(), comparator);
             Point *max = *std::max_element(points.begin(), points.end(), comparator);
 
-            float extent = max->coordinates[it_axis] - min->coordinates[it_axis];
+            float extent = max->pos[it_axis] - min->pos[it_axis];
 
             // Is it greater than max_extent?
             if (extent > max_extent) {
@@ -134,13 +80,12 @@ class KDTree {
         if (node == nullptr) { return nullptr; }
 
         // Is the left or right child node closer to this point?
-        Node *near = ray.origin[node->axis] > node->point->coordinates[node->axis] ? node->right
-                                                                                   : node->left;
+        Node *near =
+            ray.origin[node->axis] > node->point->pos[node->axis] ? node->right : node->left;
         Node *far = near == node->right ? node->left : node->right;
 
-        std::cout << "Checking " << node->point->coordinates[0] << ", "
-                  << node->point->coordinates[1] << ", " << node->point->coordinates[2] << ", "
-                  << std::endl;
+        std::cout << "Checking " << node->point->pos[0] << ", " << node->point->pos[1] << ", "
+                  << node->point->pos[2] << ", " << std::endl;
 
         // Check for collisions in this order (stopping if an intersection is found):
         // 1. In the nearer section
@@ -150,7 +95,7 @@ class KDTree {
         // If the axes are not parallel, our max_distance decreases, since we've already covered
         // some area. `t` represents the distance from this node to the splitting plane.
         float t = ray.direction[node->axis] != 0.0
-                      ? (node->point->coordinates[node->axis] - ray.origin[node->axis]) /
+                      ? (node->point->pos[node->axis] - ray.origin[node->axis]) /
                             ray.direction[node->axis]
                       : max_distance;
         Point *near_result = intersect_ray_recurse(ray, near, t);
@@ -160,21 +105,19 @@ class KDTree {
         if (near_result != nullptr) { return near_result; }
 
         // No collision in the nearer side, so check for a collision directly here
-        Point *collision_here = nullptr;
-        // TODO: Ray-triangle-intersection here
+        if (node->point->triangle && ray.intersects_triangle(node->point->triangle)) {
+            // We do have a collision here, so we're done and can return this point!
+            return node->point;
+        }
 
         // No collision here either. Does it make sense to also check the far node?
         // Only if the axes are not parallel and if that area is not behind us
         if (ray.direction[node->axis] != 0.0 && t >= 0.0) {
             // It does make sense to check the far node.
-            // For this, calculate a new ray origin ...
-            float new_origin[3]{ray.origin[0] + t * ray.direction[0],
-                                ray.origin[1] + t * ray.direction[1],
-                                ray.origin[2] + t * ray.direction[2]};
-
-            // ... and continue towards that direction, but with the new origin (we can
-            // leave behind what we already checked)
-            return intersect_ray_recurse(Ray(new_origin, ray.direction), far, max_distance - t);
+            // For this, calculate a new ray origin and continue towards that direction, but with
+            // the new origin (we can leave behind what we already checked)
+            return intersect_ray_recurse(Ray(ray.origin + ray.direction * t, ray.direction), far,
+                                         max_distance - t);
         }
 
         // If nothing worked, return a nullptr
@@ -186,10 +129,9 @@ class KDTree {
 
         Point *point = node->point;
 
-        str += std::string(depth, '-') + std::to_string(point->coordinates[0]) + ", " +
-               std::to_string(point->coordinates[1]) + ", " +
-               std::to_string(point->coordinates[2]) + " with axis " + std::to_string(node->axis) +
-               "\n";
+        str += std::string(depth, '-') + std::to_string(point->pos[0]) + ", " +
+               std::to_string(point->pos[1]) + ", " + std::to_string(point->pos[2]) +
+               " with axis " + std::to_string(node->axis) + "\n";
 
         to_string_recurse(str, node->left, depth + 1);
         to_string_recurse(str, node->right, depth + 1);