Shadow ray intersection test optimizations

src/Objects/Surface/Mesh.cpp

@@ -30,4 +30,15 @@ Mesh::intersect(const Ray& ray,
     }
     return minT < std::numeric_limits<FloatT>::infinity() ? minT : -1;
 }
+
+bool
+Mesh::intersectsBefore(const Ray& ray, FloatT maxT, FloatT epsilon) const
+{
+    for (auto& triangle : triangles) {
+        if (triangle.intersectsBefore(ray, maxT, epsilon))
+            return true;
+    }
+
+    return false;
+}
 }

src/Objects/Surface/Mesh.hpp

@@ -63,6 +63,22 @@ class Mesh : public Surface
                      LinearAlgebra::Vec3& normalOut,
                      FloatT epsilon = 0) const override;
 
+    /**
+     * @brief Finds if ray intersects this mesh with a t value lower than
+     * maxT
+     *
+     * Looks for an intersection with t in range (0, maxT)
+     *
+     * @param ray
+     * @param maxT Upper limit of t value for intersection test to succeed
+     * (exclusive)
+     * @return true There is a triangle closer than t
+     * @return false There is no such triangle
+     */
+    bool intersectsBefore(const Ray& ray,
+                          FloatT maxT,
+                          FloatT epsilon = 0) const override;
+
 protected:
     /**
      * @brief Triangles in this mesh

src/Objects/Surface/Sphere.cpp

@@ -47,4 +47,25 @@ Sphere::intersect(const Ray& ray,
     normalOut = (ray.origin + ray.direction * t - center).normalize();
     return t;
 }
+
+bool
+Sphere::intersectsBefore(const Ray& ray, FloatT maxT, FloatT epsilon) const
+{
+    // see intersect() for an explanation of this calculation
+
+    auto a = ray.direction.squaredNorm();
+    auto b = 2 * ray.direction.dot(ray.origin - center);
+    auto c = (ray.origin - center).squaredNorm() - radius * radius;
+
+    auto delta = b * b - 4 * a * c;
+
+    if (delta < 0)
+        return false;
+
+    auto t = (-b - sqrt(delta)) / (2 * a);
+    if (t <= 0)
+        t = (-b + sqrt(delta)) / (2 * a);
+
+    return t > 0 && t < maxT;
+}
 }

src/Objects/Surface/Sphere.hpp

@@ -54,6 +54,22 @@ class Sphere : public Surface
                      LinearAlgebra::Vec3& normalOut,
                      FloatT epsilon = 0) const override;
 
+    /**
+     * @brief Finds if ray intersects this sphere with a t value lower than
+     * maxT
+     *
+     * Looks for an intersection with t in range (0, maxT)
+     *
+     * @param ray
+     * @param maxT Upper limit of t value for intersection test to succeed
+     * (exclusive)
+     * @return true There is an intersection closer than t
+     * @return false There is no such intersection
+     */
+    bool intersectsBefore(const Ray& ray,
+                          FloatT maxT,
+                          FloatT epsilon = 0) const override;
+
 protected:
     /**
      * @brief Center

src/Objects/Surface/Surface.hpp

@@ -45,6 +45,22 @@ class Surface
                              LinearAlgebra::Vec3& normalOut,
                              FloatT epsilon = 0) const = 0;
 
+    /**
+     * @brief Finds if ray intersects this surface with a t value lower than
+     * maxT
+     *
+     * Looks for an intersection with t in range (0, maxT)
+     *
+     * @param ray
+     * @param maxT Upper limit of t value for intersection test to succeed
+     * (exclusive)
+     * @return true There is an intersection closer than t
+     * @return false There is no such intersection
+     */
+    virtual bool intersectsBefore(const Ray& ray,
+                                  FloatT maxT,
+                                  FloatT epsilon = 0) const = 0;
+
     /**
      * @brief Material of this Surface
      *

src/Objects/Surface/Triangle.cpp

@@ -54,6 +54,39 @@ Triangle::intersect(const Ray& ray, FloatT epsilon) const
     return t;
 }
 
+bool
+Triangle::intersectsBefore(const Ray& ray, FloatT maxT, FloatT epsilon) const
+{
+    // see intersect() for an explanation of this calculation
+
+    auto edge21 = v1 - v2;
+    auto edge31 = v1 - v3;
+    auto rhs = v1 - ray.origin;
+
+    auto detA =
+      LinearAlgebra::Mat3(ray.direction, edge21, edge31).determinant();
+
+    auto detT = LinearAlgebra::Mat3(rhs, edge21, edge31).determinant();
+    auto t = detT / detA;
+
+    if (t <= 0 || t >= maxT)
+        return false;
+
+    auto detB = LinearAlgebra::Mat3(ray.direction, rhs, edge31).determinant();
+    auto beta = detB / detA;
+
+    if (beta <= -epsilon || beta >= 1 + epsilon)
+        return false;
+
+    auto detC = LinearAlgebra::Mat3(ray.direction, edge21, rhs).determinant();
+    auto gamma = detC / detA;
+
+    if (gamma <= -epsilon || beta + gamma >= 1 + epsilon)
+        return false;
+
+    return true;
+}
+
 LinearAlgebra::Vec3
 Triangle::getNormal() const
 {

src/Objects/Surface/Triangle.hpp

@@ -52,6 +52,24 @@ class Triangle
      */
     FloatT intersect(const Ray& ray, FloatT epsilon = 0) const;
 
+    /**
+     * @brief Finds if ray intersects this triangle with a t value lower than
+     * maxT
+     *
+     * Looks for an intersection with t in range (0, maxT)
+     *
+     * @param ray
+     * @param maxT Upper limit of t value for intersection test to succeed
+     * (exclusive)
+     * @param epsilon Baricentric coordinates are allowed to be in the range
+     * (-epsilon, 1 + epsilon) instead of (0, 1)
+     * @return true This triangle is closer than t
+     * @return false There is no intersection in this range
+     */
+    bool intersectsBefore(const Ray& ray,
+                          FloatT maxT,
+                          FloatT epsilon = 0) const;
+
     /**
      * @brief Normal vector
      *

src/PathTracer/PathTracer.cpp

@@ -125,9 +125,7 @@ PathTracer::lightVisible(const LinearAlgebra::Vec3& point,
     // than the light, t > 1 means the surface is behind the light
 
     for (auto surface : scene.surfaces) {
-        LinearAlgebra::Vec3 normal;
-        auto t = surface->intersect(ray, normal, scene.intersectionTestEpsilon);
-        if (t != -1 && t < 1)
+        if (surface->intersectsBefore(ray, 1))
             return false;
     }
     return true;