From 72df366761b09f164747c1edcb34c484caad5e8d Mon Sep 17 00:00:00 2001 From: karl Date: Sat, 16 Jan 2021 22:21:08 +0100 Subject: [PATCH] Points hold multiple Triangles; fix recursive intersection Also comments and other minor fixes. --- geometry.h | 26 +++++++++++--- kdtree.h | 99 +++++++++++++++++++++++++++++++++--------------------- main.cpp | 5 +-- 3 files changed, 85 insertions(+), 45 deletions(-) diff --git a/geometry.h b/geometry.h index 506298f..847122b 100644 --- a/geometry.h +++ b/geometry.h @@ -1,5 +1,6 @@ #pragma once +#include #include // 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 triangles) : pos(pos), triangles(triangles) {} Vector pos; - Triangle *triangle; + std::list triangles; }; struct Triangle { Triangle(Vector p1, Vector p2, Vector p3) : p1(p1), p2(p2), p3(p3) {} std::vector create_point_objects() { - return std::vector{new Point(p1, this), new Point(p2, this), new Point(p3, this)}; + return std::vector{new Point(p1, std::list{this}), + new Point(p2, std::list{this}), + new Point(p3, std::list{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. diff --git a/kdtree.h b/kdtree.h index 2765e01..c2ef7c3 100644 --- a/kdtree.h +++ b/kdtree.h @@ -2,6 +2,7 @@ #include "geometry.h" #include +#include #include #include #include @@ -9,11 +10,32 @@ class KDTree { public: - KDTree(std::vector points) { root = build(points, 0); } + KDTree(std::vector 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(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(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) { - // Exit condition: There was no collision - if (node == nullptr) { return nullptr; } + void intersect_ray_recurse(const Triangle *&nearest_triangle, Vector &result, Ray ray, + Node *node, int depth, float &nearest) { + // 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): - // 1. In the nearer section - // 2. With the point in this current node - // 3. In the further section + if (ray.direction[node->axis] == 0.0) { + // The ray is parallel to the splitting axis, so we only need to check within this box. + intersect_ray_recurse(nearest_triangle, result, ray, near, depth + 1, nearest); + } 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 - // some area. `t` represents the distance from this node to the splitting plane. - float t = ray.direction[node->axis] != 0.0 - ? (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); + // Check this side for intersections up to the distance of the currently best + // intersection + intersect_ray_recurse(nearest_triangle, result, ray, near, depth + 1, nearest); - // If the nearer segment had a collision, we're done! We're only interested in the closest - // collision. - 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)) { - // We do have a collision here, so we're done and can return this point! - return node->point->triangle; + // If the far side is closer than the distance to the best current intersection, check + // that side too + if (t < nearest) { + intersect_ray_recurse(nearest_triangle, result, ray, far, depth + 1, nearest); + } } - - // 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) { diff --git a/main.cpp b/main.cpp index a7a8b80..bfb6d3b 100644 --- a/main.cpp +++ b/main.cpp @@ -16,8 +16,9 @@ int main() { std::cout << tree.to_string(); // Intersection check - Triangle *intersection = - tree.intersect_ray(Ray(new float[3]{0.5, 0.5, -1.0}, new float[3]{0.0, 0.1, 1.0})); + Vector result(0, 0, 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; }