Fix for issue #21: Limit the inverse computation to epsilon to avoid …

…false positives in the node-ray intersection function

include/bvh/node_intersectors.hpp

@@ -89,7 +89,7 @@ struct FastNodeIntersector : public NodeIntersector<Bvh, FastNodeIntersector<Bvh
     FastNodeIntersector(const Ray<Scalar>& ray)
         : NodeIntersector<Bvh, FastNodeIntersector<Bvh>>(ray) 
     {
-        inverse_direction = ray.direction.inverse();
+        inverse_direction = ray.direction.safe_inverse();
         scaled_origin     = -ray.origin * inverse_direction;
     }
 

include/bvh/vector.hpp

@@ -5,6 +5,7 @@
 #include <algorithm>
 #include <array>
 #include <cmath>
+#include <limits>
 #include <type_traits>
 
 #include "bvh/platform.hpp"
@@ -65,6 +66,13 @@ struct Vector {
         return Vector([this] (size_t i) { return Scalar(1) / values[i]; });
     }
 
+    bvh__always_inline__ Vector safe_inverse() const {
+        static constexpr auto threshold = std::numeric_limits<Scalar>::epsilon();
+        return Vector([&] (size_t i) {
+            return Scalar(1) / (std::fabs(values[i]) < threshold ? threshold : values[i]);
+        });
+    }
+
     bvh__always_inline__ Vector& operator += (const Vector& other) {
         return *this = *this + other;
     }

test/benchmark.cpp

@@ -132,7 +132,7 @@ void render(
                 traversal_steps += statistics.traversal_steps;
                 intersections   += statistics.intersections;
             }
-            if(!hit) {
+            if (!hit) {
                 pixels[index] = pixels[index + 1] = pixels[index + 2] = 0;
             } else {
                 if (CollectStatistics) {