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

geometry.h

@@ -21,37 +21,89 @@ struct Triangle {
     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)
-        : coordinates(coordinates), triangle(triangle) {}
+        : pos(Vector(coordinates)), triangle(triangle) {}
 
-    Point operator+(const Point &other) const {
+    Vector pos;
-        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) {}
+    Ray(Vector origin, Vector direction) : origin(origin), direction(direction) {}
 
-    float *origin;
+    Vector origin;
 
-    float *direction;
+    Vector direction;
 
     bool intersects_triangle(Triangle *triangle) {
-        // TODO
+        // Ray-triangle-intersection with the Möller–Trumbore algorithm
-        return false;
+        // 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

@@ -28,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];
         };
     }
 
@@ -47,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) {
@@ -80,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 *near =
-                                                                                   : node->left;
+            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] << ", "
+        std::cout << "Checking " << node->point->pos[0] << ", " << node->point->pos[1] << ", "
-                  << node->point->coordinates[1] << ", " << node->point->coordinates[2] << ", "
+                  << node->point->pos[2] << ", " << std::endl;
-                  << std::endl;
 
         // Check for collisions in this order (stopping if an intersection is found):
         // 1. In the nearer section
@@ -96,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);
@@ -106,7 +105,7 @@ class KDTree {
         if (near_result != nullptr) { return near_result; }
 
         // No collision in the nearer side, so check for a collision directly here
-        if (ray.intersects_triangle(node->point->triangle)) {
+        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;
         }
@@ -115,14 +114,10 @@ class KDTree {
         // 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 ...
+            // For this, calculate a new ray origin and continue towards that direction, but with
-            float new_origin[3]{ray.origin[0] + t * ray.direction[0],
+            // the new origin (we can leave behind what we already checked)
-                                ray.origin[1] + t * ray.direction[1],
+            return intersect_ray_recurse(Ray(ray.origin + ray.direction * t, ray.direction), far,
-                                ray.origin[2] + t * ray.direction[2]};
+                                         max_distance - t);
-
-            // ... 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);
         }
 
         // If nothing worked, return a nullptr
@@ -134,10 +129,9 @@ class KDTree {
 
         Point *point = node->point;
 
-        str += std::string(depth, '-') + std::to_string(point->coordinates[0]) + ", " +
+        str += std::string(depth, '-') + std::to_string(point->pos[0]) + ", " +
-               std::to_string(point->coordinates[1]) + ", " +
+               std::to_string(point->pos[1]) + ", " + std::to_string(point->pos[2]) +
-               std::to_string(point->coordinates[2]) + " with axis " + std::to_string(node->axis) +
+               " with axis " + std::to_string(node->axis) + "\n";
-               "\n";
 
         to_string_recurse(str, node->left, depth + 1);
         to_string_recurse(str, node->right, depth + 1);