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 &lt; 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.
Here is how I found this can be solved.

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

  float tmin = std::numeric_limits&lt;float&gt;::lowest();
  float tmax = std::numeric_limits&lt;float&gt;::max();


Add an additional check:

if(rayDir.x == 0.0f)
{
  if (rayOrigin.x() &lt; box.min.x() || rayOrigin.x() &gt; 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.