Scratchapixel 2.0
Sign in
A Minimal Ray-Tracer: Rendering Simple Shapes (Sphere, Cube, Disk, Plane, etc.)

Your are subscriber and logged in!

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 we also designate by the acronym AABB (axis-aligned bounding box). Boxes 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 equation:

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. If this 1D line is parallel to the x-axis as shown in figure 1 (with the line whose equation is y=2), then \(m\) is equal to 0.

The ray can also be expressed with the following equation:

$$O + D * t$$

where \(O\) corresponds to the origin of the ray, \(D\) its direction and the parameter \(t\) can be any real value (negative, positive or zero). By changing \(t\) in this equation, we can define 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 you replace \(O\) with \(b\) and \(D\) 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).

class Box3 { public: Box3(cont Vec3f &vmin, const Vec3f &vmax) { bounds[0] = vmin; bounds[1] = vmax; } Vec3f 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

\(B0x\) in this example, corresponds to bounds[0].x in the code above. To find where the ray intersects this line we can write:

Ox + tDx = B0x (eq1)

Which can solved by reordering the terms:

t0x = (B0x - Ox) / Dx (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) / Dx t1x = (B1x - Ox) / Dx t0y = (B0y - Oy) / Dy t1y = (B1y - Oy) / Dy t0z = (B0z - Oz) / Dz t1z = (B1z - Oz) / Dz

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) / Dx t0y = (B0y - Oy) / Dy 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) / Dx t1y = (B1y - Oy) / Dy tmax = (t1x < t1y) ? t1x : t1y

Figure 3: a ray intersecting a 2D box.

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\).

Let's 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) / Dz t1z = (B1z - Oz) / Dz if (tmin > t1z || t0z > tmax) return false if (t0z > tmin) tmin = t0z if (t1z < tmax) tmax = t1z

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

bool intersect(const Ray &r) { float tmin = (min.x - r.orig.x) / r.dir.x; float tmax = (max.x - r.orig.x) / r.dir.x; if (tmin > tmax) swap(tmin, tmax); float tymin = (min.y - r.orig.y) / r.dir.y; float 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; float tzmin = (min.z - r.orig.z) / r.dir.z; float 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; 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\).

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:

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:

Vec3f 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:

class Ray { public: Ray(const Vec3f &orig, const Vec3f &dir) : orig(orig), dir(dir) { invdir = 1 / dir; sign[0] = (invdir.x < 0); sign[1] = (invdir.y < 0); sign[2] = (invdir.z < 0); } Vec3 orig, dir; // ray orig and dir Vec3 invdir; int sign[3]; }; bool intersect(const Ray &r) const { float 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; 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 below), 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 deal with these values when you return from the ray-box intersection function.

Check the file raybox.cpp (which you can find in the last chapter of this lesson) for a complete example.

Reference

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