add ray range interection algorithm

src/Magnum/Math/Intersection.h

@@ -170,6 +170,22 @@ for plane normal @f$ \boldsymbol n @f$ and determinant @f$ w @f$.
 */
 template<class T> bool rangeFrustum(const Range3D<T>& range, const Frustum<T>& frustum);
 
+/**
+@brief Intersection of a ray with a range.
+@param rayOrigin    Origin of the ray
+@param invRayDir  Componentwise inverse of the ray direction
+@param range Range
+@return @cpp true @ce if the the ray intersect the range, @cpp false @ce otherwise.
+
+Note that you need to pass the inverse ray direction and not the ray direction.
+The purpose for this is to reduce the number of times you have to compute
+a ray inverse, when doing multiple ray range intersections (For example
+when traversing an aabb-tree).
+The algorithm implemented is a version of the classical slabs algorithm, see
+listing 1 in [Majercik et al.](http://jcgt.org/published/0007/03/04/).
+*/
+template<class T> bool rayRange(const Vector3<T>& rayOrigin, const Vector3<T>& invRayDir, const Range3D<T>& range);
+
 /**
 @brief Intersection of an axis-aligned box and a frustum
 @param aabbCenter   Center of the AABB
@@ -428,6 +444,14 @@ template<class T> bool rangeFrustum(const Range3D<T>& range, const Frustum<T>& f
     return true;
 }
 
+template<class T> bool rayRange(const Vector3<T>& rayOrigin, const Vector3<T>& invRayDir, const Range3D<T>& range) {
+    const Vector3<T> t0 = (range.min() - rayOrigin)*invRayDir;
+    const Vector3<T> t1 = (range.max() - rayOrigin)*invRayDir;
+    const std::pair<Vector3<T>, Vector3<T>> tminMax = minmax(t0,t1);
+
+    return tminMax.first.max() <= tminMax.second.min();
+}
+
 template<class T> bool aabbFrustum(const Vector3<T>& aabbCenter, const Vector3<T>& aabbExtents, const Frustum<T>& frustum) {
     for(const Vector4<T>& plane: frustum) {
         const Vector3<T> absPlaneNormal = Math::abs(plane.xyz());

src/Magnum/Math/Test/IntersectionTest.cpp

@@ -51,6 +51,7 @@ struct IntersectionTest: Corrade::TestSuite::Tester {
     void sphereConeView();
     void sphereConeViewNotRigid();
     void rangeCone();
+    void rayRange();
     void aabbCone();
 };
 
@@ -81,6 +82,7 @@ IntersectionTest::IntersectionTest() {
               &IntersectionTest::sphereConeView,
               &IntersectionTest::sphereConeViewNotRigid,
               &IntersectionTest::rangeCone,
+              &IntersectionTest::rayRange,
               &IntersectionTest::aabbCone});
 }
 
@@ -434,6 +436,55 @@ void IntersectionTest::rangeCone() {
         center, normal, angle));
 }
 
+
+void IntersectionTest::rayRange() {
+    const Vector3 origin{2,2,2};
+    const Range3D range{{-1, -1, -1},
+                        {1,  1,  1}};
+
+    const Vector3 center{0,0,1};
+    const Vector3 edge{0,-1,1};
+    const Vector3 corner{-1,-1,1};
+    const Float eps = 1e-6;
+
+    const Vector3 ones{1,1,1};
+
+    /* intersection at face center */
+    const Vector3 direction1 = center - origin;
+    const Vector3 invDir1 = ones/direction1;
+    CORRADE_VERIFY(Intersection::rayRange(origin, invDir1, range));
+
+    /* intersection close to edge */
+    const Vector3 direction2 = edge + Vector3{0,eps,0} - origin;
+    const Vector3 invDir2 = ones/direction2;
+    CORRADE_VERIFY(Intersection::rayRange(origin, invDir2, range));
+
+    /* no intersection close to edge */
+    const Vector3 direction3 = edge - Vector3{0,eps,0} - origin;
+    const Vector3 invDir3 = ones/direction3;
+    CORRADE_VERIFY(!Intersection::rayRange(origin, invDir3, range));
+
+    /* intersection close to corner */
+    const Vector3 direction4 = corner + Vector3{eps,eps,0} - origin;
+    const Vector3 invDir4 = ones/direction4;
+    CORRADE_VERIFY(Intersection::rayRange(origin, invDir4, range));
+
+    /* no intersection close to corner */
+    const Vector3 direction5 = corner - Vector3{eps,eps,0} - origin;
+    const Vector3 invDir5 = ones/direction5;
+    CORRADE_VERIFY(!Intersection::rayRange(origin, invDir5, range));
+
+    /* divide by zero test with intersection */
+    const Vector3 direction6{0,0,-1};
+    const Vector3 invDir6 = ones/direction6;
+    CORRADE_VERIFY(Intersection::rayRange({0,0,2}, invDir6, range));
+
+    /* divide by zero test without intersection */
+    const Vector3 direction7{0,0,1};
+    const Vector3 invDir7 = ones/direction7;
+    CORRADE_VERIFY(!Intersection::rayRange(origin, invDir7, range));
+}
+
 void IntersectionTest::aabbCone() {
     const Vector3 center{1.0f, -2.0f, 1.3f};
     const Vector3 normal{0.453154f, 0.422618f, 0.784886f};