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

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.

• Nuts, I have this bug in some of my code... I think, in this case, you want 0 * inf to be 0, since it happens at the origin of the Ray. That should result in a ray lying on one face of the box to be a hit, the opposite face a miss. But you need to check for the zero. Feb 1 '20 at 15:27
• I think the pragmatic fix for BVH use would be to inflate the boxes by 1 ulps. I suspect that would beat any fix I can think of as all of them involve several extra instructions. I will keep thinking about it though, it is an interesting problem. Feb 2 '20 at 0:30
• Making the boxes physically larger in my case would make the results totally worthless. I can check if the result is NAN but that doesn't guarantee a hit, only that the origin started on one of those dimensions.Also, making them larger would just mean it would be the same problem, just at a different point in space, no? Or do you mean, if I detect a nan, extend the box and then try again? Feb 3 '20 at 2:50
• @TylerShellberg I meant always make them larger. This will give a tiny % of false positives, which for a BVH is of no importance. If you require an exact intersection test, you might want to edit your question to specify that and also if you care about fixing the problem on both the min/max side of each axis or only on one side. Feb 4 '20 at 1:13
• @Oliver I am not using this for a BVH. Feb 4 '20 at 18:25

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();

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.

• You can specialize the intersection function based on the zero components (with 7 different versions), and avoid the runtime cost of the checks as you traverse for a given ray. Feb 3 '20 at 21:00
• @DanielMGessel Could you clarify what you mean by that? Feb 3 '20 at 21:09
• Yeah, guess that was a bit sloppy. If you’re traversing a BVH for a single ray (on a scalar processor, not a gpu) consider having a version of the traversal code for the general case where no component of the ray direction is zero (most rays), eliminating your added checks, and a version with the checks. Select traversal code based on ray direction, Might speed things up on a cpu. Feb 4 '20 at 0:25
• Just grokked the pbr book’s solution. I’m slow today. Check it out: pbr-book.org/3ed-2018/Shapes/Basic_Shape_Interface.html Feb 4 '20 at 0:44
• @DanielMGessel How is checking what type of ray, then choosing a function, rather than the function determining it any faster? Either way, the exact same checks need to happen on the ray. It doesn't seem like less work is being done. Feb 4 '20 at 15:48

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;
}