Points hold multiple Triangles; fix recursive intersection

Also comments and other minor fixes.

geometry.h

@@ -1,5 +1,6 @@
 #pragma once
 
+#include <list>
 #include <vector>
 
 // Forward declarations
@@ -21,6 +22,16 @@ struct Vector {
         return Vector(c[0] - other.c[0], c[1] - other.c[1], c[2] - other.c[2]);
     }
 
+    // Arbitrary but functional definition of '<', primarily for use in an std::map.
+    bool operator<(const Vector &other) const {
+        if ((c[2] < other.c[2])) { return true; }
+        if ((c[2] == other.c[2]) && (c[1] < other.c[1])) { return true; }
+        if ((c[2] == other.c[2]) && (c[1] == other.c[1]) && (c[0] < other.c[0])) { return true; }
+
+        return false;
+    }
+
+    // Component-wise multiplication with a scalar
     Vector operator*(float scalar) const {
         return Vector(c[0] * scalar, c[1] * scalar, c[2] * scalar);
     }
@@ -36,18 +47,22 @@ struct Vector {
 };
 
 struct Point {
-    Point(Vector pos, Triangle *triangle) : pos(pos), triangle(triangle) {}
+    Point(Vector pos) : pos(pos) {}
+
+    Point(Vector pos, std::list<Triangle *> triangles) : pos(pos), triangles(triangles) {}
 
     Vector pos;
 
-    Triangle *triangle;
+    std::list<Triangle *> triangles;
 };
 
 struct Triangle {
     Triangle(Vector p1, Vector p2, Vector p3) : p1(p1), p2(p2), p3(p3) {}
 
     std::vector<Point *> create_point_objects() {
-        return std::vector<Point *>{new Point(p1, this), new Point(p2, this), new Point(p3, this)};
+        return std::vector<Point *>{new Point(p1, std::list<Triangle *>{this}),
+                                    new Point(p2, std::list<Triangle *>{this}),
+                                    new Point(p3, std::list<Triangle *>{this})};
     }
 
     Vector p1;
@@ -74,7 +89,7 @@ struct Ray {
 
     Vector direction;
 
-    bool intersects_triangle(Triangle *triangle) {
+    bool intersects_triangle(const Triangle *triangle, Vector &result, float &t) {
         // 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;
@@ -103,8 +118,9 @@ struct Ray {
 
         // At this stage we can compute t to find out where the intersection point is on the
         // line.
-        float t = f * edge2.dot(q);
+        t = f * edge2.dot(q);
         if (t > EPSILON) {
+            result = origin + direction * t;
             return true;
         } else {
             // This means that there is a line intersection but not a ray intersection.

kdtree.h

@@ -2,6 +2,7 @@
 
 #include "geometry.h"
 #include <algorithm>
+#include <chrono>
 #include <glm/glm.hpp>
 #include <iostream>
 #include <string>
@@ -9,11 +10,32 @@
 
 class KDTree {
   public:
-    KDTree(std::vector<Point *> points) { root = build(points, 0); }
+    KDTree(std::vector<Point *> points) {
+        std::cout << "Starting to build tree with " << points.size() << " points took ";
+        auto start = std::chrono::high_resolution_clock::now();
+        root = build(points, 0);
+        auto end = std::chrono::high_resolution_clock::now();
+
+        std::cout << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
+                  << " microseconds" << std::endl;
+    }
 
     ~KDTree() = default; // TODO: Delete all allocated Nodes
 
-    Triangle *intersect_ray(Ray ray) { return intersect_ray_recurse(ray, root, 1000.0, 0); }
+    const Triangle *intersect_ray(const Ray ray, Vector &result) {
+        auto start = std::chrono::high_resolution_clock::now();
+
+        float nearest = 1000.0; // Initial max distance
+        const Triangle *nearest_triangle = nullptr;
+        intersect_ray_recurse(nearest_triangle, result, ray, root, 0, nearest);
+
+        auto end = std::chrono::high_resolution_clock::now();
+
+        std::cout << "Intersection took "
+                  << std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()
+                  << " microseconds" << std::endl;
+        return nearest_triangle;
+    }
 
     std::string to_string() {
         std::string str = "";
@@ -87,50 +109,51 @@ class KDTree {
                         build(right_of_median, depth + 1));
     }
 
-    Triangle *intersect_ray_recurse(Ray ray, Node *node, float max_distance, int depth) {
+    void intersect_ray_recurse(const Triangle *&nearest_triangle, Vector &result, Ray ray,
-        // Exit condition: There was no collision
+                               Node *node, int depth, float &nearest) {
-        if (node == nullptr) { return nullptr; }
+        // Exit condition: This node was iterated towards, but does not exist
+        if (node == nullptr) { return; }
 
-        // Is the left or right child node closer to this point?
+        // Check for a collision here
+        // Iterate over all Triangles which this Point is involved in
+        for (const Triangle *triangle : node->point->triangles) {
+            Vector current_result(0, 0, 0);
+            float current_distance;
+
+            // If we have a collision, and it is closer to the ray origin than the closest previous
+            // collision, remember it
+            if (ray.intersects_triangle(triangle, current_result, current_distance)) {
+                if (current_distance < nearest) {
+                    nearest = current_distance;
+                    nearest_triangle = triangle;
+                    result = current_result;
+                }
+            }
+        }
+
+        // Is the ray origin within the left or right child node bounding box?
         Node *near =
             ray.origin[node->axis] > node->point->pos[node->axis] ? node->right : node->left;
         Node *far = near == node->right ? node->left : node->right;
 
-        // Check for collisions in this order (stopping if an intersection is found):
+        if (ray.direction[node->axis] == 0.0) {
-        // 1. In the nearer section
+            // The ray is parallel to the splitting axis, so we only need to check within this box.
-        // 2. With the point in this current node
+            intersect_ray_recurse(nearest_triangle, result, ray, near, depth + 1, nearest);
-        // 3. In the further section
+        } else {
+            // Calculate the distance from the ray origin to the splitting axis
+            float t =
+                (node->point->pos[node->axis] - ray.origin[node->axis]) / ray.direction[node->axis];
 
-        // If the axes are not parallel, our max_distance decreases, since we've already covered
+            // Check this side for intersections up to the distance of the currently best
-        // some area. `t` represents the distance from this node to the splitting plane.
+            // intersection
-        float t = ray.direction[node->axis] != 0.0
+            intersect_ray_recurse(nearest_triangle, result, ray, near, depth + 1, nearest);
-                      ? (node->point->pos[node->axis] - ray.origin[node->axis]) /
-                            ray.direction[node->axis]
-                      : max_distance;
-        Triangle *near_result = intersect_ray_recurse(ray, near, t, depth + 1);
 
-        // If the nearer segment had a collision, we're done! We're only interested in the closest
+            // If the far side is closer than the distance to the best current intersection, check
-        // collision.
+            // that side too
-        if (near_result != nullptr) { return near_result; }
+            if (t < nearest) {
-
+                intersect_ray_recurse(nearest_triangle, result, ray, far, depth + 1, nearest);
-        // No collision in the nearer side, so check for a collision directly here
+            }
-        if (ray.intersects_triangle(node->point->triangle)) {
-            // We do have a collision here, so we're done and can return this point!
-            return node->point->triangle;
         }
-
-        // 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 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, depth + 1);
-        }
-
-        // If nothing worked, return a nullptr
-        return nullptr;
     }
 
     void to_string_recurse(std::string &str, Node *node, int depth) {

main.cpp

@@ -16,8 +16,9 @@ int main() {
     std::cout << tree.to_string();
 
     // Intersection check
-    Triangle *intersection =
+    Vector result(0, 0, 0);
-        tree.intersect_ray(Ray(new float[3]{0.5, 0.5, -1.0}, new float[3]{0.0, 0.1, 1.0}));
+    const Triangle *intersection =
+        tree.intersect_ray(Ray(new float[3]{0.5, 0.5, -1.0}, new float[3]{0.0, 0.1, 1.0}), result);
 
     if (intersection != nullptr) { std::cout << "Hit!" << std::endl; }