Ray-Box Intersection

Ray-Box Intersection

Figure 1: equation of a line. The equation of a line can be written as y=mx+b. For the oblique line which equation is y=x-1 we have m=1 and b=-1.

In the following example we will assume that the box is aligned with the axis of our coordinate system. Such a box is called an axis-aligned box and since it will often be a bounding box (or bounding volume) which usually call them AABB (axis-aligned bounding box). Other boxes, those which are not oriented along the axis of a cartesian coordinate system, are called OBB (oriented bounding box).

Computing the intersection of a ray with a AABB is quite simple. All we need is to remember that a straight line can be defined using the following analytical form:

y = mx + b

In mathematics, the m term is called the slope (or gradient) and is actually responsible for the orientation of the line and b corresponds to the point where the line intersects the y-axis.

The ray can also be expressed with the following equation:

O + Rt

where O corresponds to the origin of the ray, R its direction and the parameter t can be any real value (negative, positive or zero). By changing t in this equation, we can create any point on the line defined by the ray's position and direction. This ray equation is very similar to the line equation. In fact they are the same if your replace O with b and R with m. To represent an axis-aligned bounding volume, all we need are two points representing the minimum and maximum extent of the box (called bounds in the code).

1
2
3
4
5
6
7
8
9
10
11
template<typename T> class Box3 { public: Box3(Vec3<T> vmin, Vec3<T> vmax) { bounds[0] = vmin; bounds[1] = vmax; } Vec3<T> bounds[2]; };

The bounds of the volume define a set of lines parallel to each axis of the coordinate system which we can also expressed using the line equation. For instance the line equation for the x component of the bounding volume's minimum extent can be written as:

y = B0x

To find where the ray intersects with this line we can write:

Ox + tRx = B0x (eq1)

Which can solved by re-ordering the terms:

t0x = ( B0x - Ox ) / Rx (eq2)

The x component of the bounding volume's maximum extent can be used in a similar way to compute t1x. Note than when the values for t are negative, the box is behind the ray.

Finally if we apply the same technique to the y and z components we will have at the end of this process a set of six values indicating where the ray intersects the box planes parallel to the x, y and z axis.

t0x = ( B0x - Ox ) / Rx
t1x = ( B1x - Ox ) / Rx
t0y = ( B0y - Oy ) / Ry
t1y = ( B1y - Oy ) / Ry
t0z = ( B0z - Oz ) / Rz
t1z = ( B1z - Oz ) / Rz

Figure 2: a ray intersecting a 2D box.

Note that if the ray is parallel to an axis it won't intersect with the bounding volume plane for this axis (in this case, the line equation for the ray is reduced to the constant b and there's no solution to equation 1). The problem at this stage is to find which of these six values correspond to an intersection of the ray with the box (if the ray intersects the box at all. So far, we have just found where the ray intersects the planes defined by each face of the cube). This can be easily found by operating a simple test on each of the computed t values. By looking at figure 2, the logic behind this test will become obvious (we will just be considering the 2D case for this example). The ray first intersects the planes defined by the minimum extent of the box in two places: t0x and t0y. However intersecting these planes doesn't necessarily mean that these intersecting points lie on the cube (if they don't lie on the cube, obviously the ray doesn't intersect the box). Mathematically we can find which one of these two points lie on the cube by comparing their values: we simply need to chose the point which value for t is the greatest. In pseudo-code we can write:

t0x = (B0x - Ox) / Rx
t0y = (B0y - Oy) / Ry
tmin = (t0x > t0y) ? t0x : t0y

The process to find the second point where the ray intersects the box is similar however in that case, we will use the t values computed using the planes defined by the maximum extent of the box and select the point which t value is the smallest.

t1x = (B1x - Ox) / Rx
t1y = (B1y - Oy) / Ry
tmax = (t1x < t1y) ? t1x : t1y

Figure 3: a ray intersecting a 2D box.

However the ray doesn't alway intersect the box. In figure 3, we are showing a couple of cases where the ray misses the cube. These cases can also be easily identified using simple comparisons between the t values. As you can see on the figure, the ray misses the box when t0x is greater than t1y and when t0y is greater than t1x.

Lets add this test to our code:

if (t0x > t1y || t0y > t1x)
return false

Finally we can extend the technique to the 3D case by computing t values for the z component and compare them to the t values we have computed so far for the x and y components:

t0z = (B0z - Oz) / Rz
t1z = (B1z - Oz) / Rz
if (tmin > t1z || t0z > tmax) return false
if (t0z > tmin) tmin = t0z
if (t1z < tmax) tmax = t1z

At this point, we know for sure that the ray intersect the box. However depending on the implementation you are using, you may have set a tmin and tmax values on the ray itself. You may need to check the tmin and tmax values you have computed in the ray-box intersection with the tmin and tmax values set on the ray before returning that the ray intersects the box:

if (tmin > ray.tmax || tmax < ray.tmin) return false

If we pass this test we can set the the ray tmin and tmax value with the tmin and tmax values from the function. These corresponds to the parametric values for t where the ray intersects the box:

if (tmin > ray.tmin) ray.tmin = tmin
if (tmax < ray.tmax) ray.tmax = tmax

Here is a full implementation of the method in C++ (min and max are the minimum and maximum extent of the bounding box):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
bool intersect(const Ray<T> &r)
{
T tmin = (min.x - r.orig.x) / r.dir.x;
T tmax = (max.x - r.orig.x) / r.dir.x;
if (tmin > tmax) swap(tmin, tmax);
T tymin = (min.y - r.orig.y) / r.dir.y;
T tymax = (max.y - r.orig.y) / r.dir.y;
if (tymin > tymax) swap(tymin, tymax);
if ((tmin > tymax) || (tymin > tmax))
return false;
if (tymin > tmin)
tmin = tymin;
if (tymax < tmax)
tmax = tymax;
T tzmin = (min.z - r.orig.z) / r.dir.z;
T tzmax = (max.z - r.orig.z) / r.dir.z;
if (tzmin > tzmax) swap(tzmin, tzmax);
if ((tmin > tzmax) || (tzmin > tmax))
return false;
if (tzmin > tmin)
tmin = tzmin;
if (tzmax < tmax)
tmax = tzmax;
if ((tmin > r.tmax) || (tmax < r.tmin)) return false;
if (r.tmin < tmin) r.tmin = tmin;
if (r.tmax > tmax) r.tmax = tmax;
return true;
}

Note that depending on the ray direction tmin might be greater than tmax. Considering that all the logic behind the tests we are doing later on relies on t0 being always smaller than t1 we need to swap the values if t1 is smaller than t0 (line 5, 8 and 17).

Optimising the Code

There is a couple of improvement we can make to the previous code to not only make it faster but also more robust. First we can replace the swap operation with the following code:

1
2
3
4
5
6
7
8
if (ray.dir.x >= 0) { tmin = (min.x - r.orig.x) / ray.dir.x; tmax = (max.x - r.orig.x) / ray.dir.x; } else { tmin = (max.x - r.orig.x) / ray.dir.x; tmax = (min.x - r.orig.x) / ray.dir.x; }

However there is a problem with this code. When the direction of ray is 0 it causes a division by zero. This though should be handled properly by the compiler which shall return +∞. The problem happens when the value for the ray direction is -0. When the ray direction is negative, we expect the second block of the if statement to be executed, but in the case where the ray direction is -0, the first block will be executed instead because in IEEE float point, the test -0 = 0 returns true. The values tmin/tmax will be set to +∞ and -∞ (instead of -∞/+∞) and the code will fail to detect an intersection. We can easily fix this problem by replacing the ray direction by the inverse of the ray direction:

1
2
3
4
5
6
7
8
9
Vec3<T> invdir = 1 / ray.dir; if (invdir.x >= 0) { tmin = (min.x - r.orig.x) * invdir.x; tmax = (max.x - r.orig.x) * invdir.x; } else { tmin = (max.x - r.orig.x) * invdir.x; tmax = (min.x - r.orig.x) * invdir.x; }

When the ray direction is -0 the inverse direction should be -∞ and the test on line 2 shall return false which is what we want. In the eventuality of testing the intersection of the ray against many boxes, we can save some time by pre-computing the inverse direction of the ray in the Ray's constructor and re-using it later in intersect. We can also pre-compute the sign of the ray direction which we can use to write the method in a much more compact form. Here is the final version of the ray-box intersection method:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
template<typename T> class Ray { public: Ray(Vec3<T> orig, Vec3<T> dir) : orig(orig), dir(dir), tmin(T(0)), tmax(std::numeric_limits<T>::max()) { invdir = T(1) / dir; sign[0] = (invdir.x < 0); sign[1] = (invdir.y < 0); sign[2] = (invdir.z < 0); } Vec3<T> orig, dir; /// ray orig and dir mutable T tmin, tmax; /// ray min and max distances Vec3<T> invdir; int sign[3]; }; bool intersect(const Ray<T> &r) const { T tmin, tmax, tymin, tymax, tzmin, tzmax; tmin = (bounds[r.sign[0]].x - r.orig.x) * r.invdir.x; tmax = (bounds[1-r.sign[0]].x - r.orig.x) * r.invdir.x; tymin = (bounds[r.sign[1]].y - r.orig.y) * r.invdir.y; tymax = (bounds[1-r.sign[1]].y - r.orig.y) * r.invdir.y; if ((tmin > tymax) || (tymin > tmax)) return false; if (tymin > tmin) tmin = tymin; if (tymax < tmax) tmax = tymax; tzmin = (bounds[r.sign[2]].z - r.orig.z) * r.invdir.z; tzmax = (bounds[1-r.sign[2]].z - r.orig.z) * r.invdir.z; if ((tmin > tzmax) || (tzmin > tmax)) return false; if (tzmin > tmin) tmin = tzmin; if (tzmax < tmax) tmax = tzmax; if (tmin > r.tmin) r.tmin = tmin; if (tmax < r.tmax) r.tmax = tmax; return true; }

The optimised version runs on our machine approximately 1.5 faster than the first version (the speedup might depend on the compiler).

The code from this lesson returns intersections with the box which are in front or behind the origin of the ray. For instance, if the ray's origin is inside the box (like in the image on the right), there will be two intersections: one in front of the ray and one behind. We know that an intersection is "behind" the origin of the ray when the value for t is negative. When t is positive, the intersection is in front of the origin of the ray. If your algorithm is not interested in intersections for values of t lower than 0, then you will have to carefully deal with these cases when you return from the ray-box intersection box (as it is often a source of bugs).

Reference

An Efficient and Robust Ray–Box Intersection Algorithm, Amy Williams et al. 2004.

Chapter 1 of 3