Need a fast ray-box intersection that handles if a ray is parallel and in line with a plane of the box

Question

I've seen lots of spins on the ray-box intersection test, and a lot of them seem to boil down to this code: (omitting some details here)

  invRay.x = 1.0 / ray.x;

  tx1 = (box.min.x - rayOrigin.x) * invRay.x;
  tx2 = (box.max.x - rayOrigin.x) * invRay.x;

  tmin = min(tx1, tx2);
  tmax = max(tx1, tx2);

  // Repeat for Y and Z
  // Then check to ensure tmin < tmax

(I got mine specifically from this discussion and this source. Also note that if you use std::min and std::max, it's far from branchless)

However, this raises a significant problem.

What if the ray's origin in one or more dimensions is equal to the min or max points?

What I found is that I get this situation on an example ray and box:

• Ray origin: (3.225, 0.243, 2.221)
• Ray direction: (0, 0, -1)
• Box min: (3.223, 0.241, 2.119)
• Box max: (3.225, 0.243, 2.121)

Is that when tx1 is calculated, it correctly gets -inf. However, when tx2 is calculated, it gets (3.225 - 3.225) * inf which is coming out to -nan. This means that later, when we do:

float tmax = std::max(tx1, tx2);

We get max(-inf, -nan) which comes out to -inf and from that point on in the tests, tmax will always be -inf, making the overall ray-box test fail.

This is confirmed by some comments on the page where I got the code:

Something to consider here is that 0 * inf =nan which occurs when the ray starts exactly on the edge of a box

And

Having read both this and part 2, I'm unable to work out what changes I need to make in order to reliably detect intersections when a ray is exactly on the edge of an AABB.

Does anyone know a fast method for ray-box intersection with edge and corner hits?

I should note this is not being used with a BVH. I need an exact intersection test, I cannot increase the size of the box at all.

Answer 1

Here is how I found this can be solved.

(The below also happens at the start, just fyi)

  float tmin = std::numeric_limits<float>::lowest();
  float tmax = std::numeric_limits<float>::max();

Add an additional check:

if(rayDir.x == 0.0f)
{
  if (rayOrigin.x() < box.min.x() || rayOrigin.x() > box.max.x())
  {     
    return false;
  }
}
else 
{
  tx1 = (box.min.x - rayOrigin.x) * invRay.x;
  tx2 = (box.max.x - rayOrigin.x) * invRay.x;

  tmin = std::max(tmin, std::min(tx1, tx2));
  tmax = std::min(tmax, std::max(tx1, tx2));
}

This way the tx1/tx2 values never get calculated if the ray's direction vector is zero in that direction, aka if it's parallel to that plane.

It's not particularly compact or pretty, but it does neatly solve the issue.

Answer 2

This is a C++ implementation. Slab is the space between a pair of parallel planes.Test if ray R(t) = p + t * d intersect against bounding box bb.

bool IntersectRayBB( const double* p, const double* d, const BoundingBox& bb, double& tmin, double* q )
{
    double tolerance = 1e-9;
    tmin = 0.0;
    double tmax = FLOAT32_MAX; // Set to max distance ray can travel
    for( int i = 0; i < 3; i++ )
    {
        if( Abs( d[i] ) < tolerance )
        {
            // Ray is parallel to slab. No hit if origin not within slab.
            if( p[i] < bb.min[i] || p[i] > bb.max[i] )
                return false;
        }
        else
        {
            // Compute intersection t value of ray within near and far plane of slab
            double denom = 1.0 / d[i];
            double t1 = ( bb.min[i] - p[i] ) * denom;
            double t2 = ( bb.max[i] - p[i] ) * denom;
            // Make t1 be intersection with near plane, t2 with far plane
            if( t1 > t2 )
                std::Swap( t1, t2 );
            // Compute the intersection of slab intersection intervals
            if( t1 > tmin )
                tmin = t1;
            if( t2 > tmax )
                tmax = t2;
            // Exit with no collision as soon as slab intersection becomes empty
            if( tmin > tmax )
                return false;
        }
    }

    //Ray intersects all 3 slabs. Return point (q) and intersesction t value (tmin)
    q = p + d * tmin;
    return true;
}