(Update: part 2)

Axis-aligned bounding boxes (AABBs) are universally used to bound finite objects in ray-tracing. Ray/AABB intersections are usually faster to calculate than exact ray/object intersections, and allow the construction of bounding volume hierarchies (BVHs) which reduce the number of objects that need to be considered for each ray. (More on BVHs in a later post.) This means that a ray-tracer spends a lot of its time calculating ray/AABB intersections, and therefore this code ought to be highly optimised.

The fastest method for performing ray/AABB intersections is the slab method. The idea is to treat the box as the space inside of three pairs of parallel planes. The ray is clipped by each pair of parallel planes, and if any portion of the ray remains, it intersected the box.

A simple implementation of this algorithm might look like this (in two dimensions for brevity):

bool intersection(box b, ray r) { double tmin = -INFINITY, tmax = INFINITY; if (ray.n.x != 0.0) { double tx1 = (b.min.x - r.x0.x)/r.n.x; double tx2 = (b.max.x - r.x0.x)/r.n.x; tmin = max(tmin, min(tx1, tx2)); tmax = min(tmax, max(tx1, tx2)); } if (ray.n.y != 0.0) { double ty1 = (b.min.y - r.x0.y)/r.n.y; double ty2 = (b.max.y - r.x0.y)/r.n.y; tmin = max(tmin, min(ty1, ty2)); tmax = min(tmax, max(ty1, ty2)); } return tmax >= tmin; } |

However, those divisions take quite a bit of time. Since when ray-tracing, the same ray is tested against many AABBs, it makes sense to pre-calculate the inverses of the direction components of the ray. If we can rely on the IEEE 754 floating-point properties, this also implicitly handles the edge case where a component of the direction is zero - the `tx1`

and `tx2`

values (for example) will be infinities of opposite sign if the ray is within the slabs, thus leaving `tmin`

and `tmax`

unchanged. If the ray is outside the slabs, `tx1`

and `tx2`

will be infinities with the same sign, thus making `tmin == +inf`

or `tmax == -inf`

, and causing the test to fail.

The final implementation would look like this:

bool intersection(box b, ray r) { double tx1 = (b.min.x - r.x0.x)*r.n_inv.x; double tx2 = (b.max.x - r.x0.x)*r.n_inv.x; double tmin = min(tx1, tx2); double tmax = max(tx1, tx2); double ty1 = (b.min.y - r.x0.y)*r.n_inv.y; double ty2 = (b.max.y - r.x0.y)*r.n_inv.y; tmin = max(tmin, min(ty1, ty2)); tmax = min(tmax, max(ty1, ty2)); return tmax >= tmin; } |

Since modern floating-point instruction sets can compute min and max without branches, this gives a ray/AABB intersection test with no branches or divisions.

My implementation of this in my ray-tracer Dimension can be seen here.

This doesn't seem to work out for me for negative direction vectors. Looking forward, I can see a box that's there, but looking backwards, I again see the same box mirrored, even though in this direction it's not "there". Following the algorithm step by step manually for both a positive and a negative z-dir (with x-dir and y-dir set to 0) gives the same near and far planes in both directions:

Right, because the test is only for whether the line intersects the box at all. The line extends both forwards and backwards. Just add a tmax >= 0 check. It's tmax, not tmin, since tmin will be < 0 if the ray originates inside the box.

Thanks. Works nicely and fast. I updated it a bit, to use SSE (though Vectormath), floats only.

https://gist.github.com/4412640#file-bbox-cpp-L14

You're welcome! I can't see that gist though (says "OAuth failure"). What kind of performance did the vectorisation give you?

To get the actual intersection, would I just use tmin, and multiply it with ray dir, adding ray origin? And what if I'm just interested in which face was intersected? x+,x-,y+,y-,z+ or z-?

Yes, that's what I'd do. Except if the ray origin is inside the box (

`tmin < 0`

), you need to use`tmax`

instead.To see what face was intersected, there's a few different ways. You can keep track of which slab is intersecting in the above algorithm at all times, but that slows it down.

For a cube centered at the origin, a neat trick is to take the component of the intersection point with the largest absolute value.

So if I return a vec3 with tmin, tmax and a float hit = step(tmin,tmax)*step(0,tmax), I basically know that if hit > 0 then ro+rd*tmin is my lower bound intersection point (entry point) and ro+rd*tmax is my higher bound intersection point (exit point), right? However, if tmin < 0, I'm inside the bounding box, which means I don't need the entry point and I can just use the ray origin as my starting point.

Sorry, not sure what you mean by

`step()`

. You can get the starting point as`ro + td*max(tmin, 0.0)`

.step(edge,A) returns 1 if A >= edge, 0 elsewise. So I know if it was a hit if tmax >= tmin AND tmax >= 0. Then we also know tmin is the closest intersection and tmax is the furthest.

I implemented this in C++ and it works perfectly for any ray direction. However, in GLSL it seems to have problems with negative directions.

Thanks Tavianator,

Yeah, the cube-trick requires sorting.

So for now I test proximity to face within an epsilon.

If close enough to face, I assume that face was hit.

I can live with the few false positives, as I shoot over 100M photons each frame anyway.

Oops... that should be 100K photons of course.

Haha I was

reallyimpressed for a second :)The cube trick does not require sorting, just selecting the max from three candidates.

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

True, if you want consistent handling of that case while staying branch-free, you have to do a little more work. This is worth another post actually, I'll write one up.

Thanks for posting this (so long ago)! I used this in my own code, and wondered if it is possible to save 2 of the subtractions by taking advantage of the fact that min(x+a, y+a) = min(x,y)+a (for some constant a):

I think that should work, but there's a couple bugs in your example code. And actually for the test as written you can omit

`r.offset`

entirely since it doesn't affect the`tmax >= tmin`

check at the end. But since you probably want a`tmax >= 0.0`

check too, this should work:I'll try it out and see how much faster it is, thanks!

What if I want to know the 't' of the intersection?

\(t = t_\min\), unless \(t_\min < 0\), in which case you're inside the box and \(t = t_\max\).

Very useful! faster the any method i tried. not because of the !(use of divisions), but because the implementation of the slab method is far simpler than others do. They use things(slow things) like square roots, DOT products and lots of checks with points. now my raycasting method has the speed i wanted. Thanks for that.

Hello, many thanks for this useful information!

I have a question about the surface normal. Is there an easy way to get the normal of the intersection point?

Thank you!

I would really like to know too!

Hi, great read!

One quick question.

I had no problem using the first implementation for my ray-tracing algorithm, but can't implement that faster version since I can't get the inverse of the direction vector(doesn't it have to be squared?).

Any help is appreciated, Thanks!

You can just compute (1/x, 1/y, 1/z) as the inverse. You don't have to square it.

Hey there,

considering I have box with min(1,1,1) and max(2,2,2) and a ray with origin(3,3,3) and direction (2,2,2). Oh. Wait a second. I think I just answered my question. It's ray direction, not target position.

If I want my ray to go from (3,3,3) towards the direction of point (2,2,2) I need a direction like (-1,-1,-1), correct? Oh my, this made me struggle way longer than it should have.

Nevermind me. Thank you for the write-up in both articles. :)

Does anyone know why this doesn't seem to work in GLSL for all angles?

INFINITY is defined as `const float INFINITY = 1.0 / 0.0;`

```

float rayAABBIntersection (Ray ray, AABB aabb)

{

float tx1 = (aabb.min.x - ray.origin.x)*ray.inverseDirection.x;

float tx2 = (aabb.max.x - ray.origin.x)*ray.inverseDirection.x;

float tmin = min(tx1, tx2);

float tmax = max(tx1, tx2);

float ty1 = (aabb.min.y - ray.origin.y)*ray.inverseDirection.y;

float ty2 = (aabb.max.y - ray.origin.y)*ray.inverseDirection.y;

tmin = max(tmin, min(ty1, ty2));

tmax = min(tmax, max(ty1, ty2));

float tz1 = (aabb.min.z - ray.origin.z)*ray.inverseDirection.z;

float tz2 = (aabb.max.z - ray.origin.z)*ray.inverseDirection.z;

tmin = max(tmin, min(tz1, tz2));

tmax = min(tmax, max(tz1, tz2));

if (tmin > tmax)

{

return INFINITY;

}

if (tmin < 0)

{

return tmax;

}

return tmin;

}

```