Ray-box Intersection Theory

Question

Hello:

I wish to determine the intersection point between a ray and a box. The box is defined by its min 3D coordinate and max 3D coordinate and the ray is defined by its origin and the direction to which it points.

Currently, I am forming a plane for each face of the box and I'm intersecting the ray with the plane. If the ray intersects the plane, then I check whether or not the intersection point is actually on the surface of the box. If so, I check whether it is the closest intersection for this ray and I return the closest intersection.

The way I check whether the plane-intersection point is on the box surface itself is through a function

bool PointOnBoxFace(R3Point point, R3Point corner1, R3Point corner2)
{
  double min_x = min(corner1.X(), corner2.X());
  double max_x = max(corner1.X(), corner2.X());
  double min_y = min(corner1.Y(), corner2.Y());
  double max_y = max(corner1.Y(), corner2.Y());
  double min_z = min(corner1.Z(), corner2.Z());
  double max_z = max(corner1.Z(), corner2.Z());
  if(point.X() >= min_x && point.X() <= max_x && 
     point.Y() >= min_y && point.Y() <= max_y &&
     point.Z() >= min_z && point.Z() <= max_z)
     return true;

  return false;
}

where corner1 is one corner of the rectangle for that box face and corner2 is the opposite corner. My implementation works most of the time but sometimes it gives me the wrong intersection. Please see image:

The image shows rays coming from the camera's eye and hitting the box surface. The other rays are the normals to the box surface. It can be seen that the one ray in particular (it's actually the normal that is seen) comes out from the "back" of the box, whereas the normal should be coming up from the top of the box. This seems to be strange since there are multiple other rays that hit the top of the box correctly.

I was wondering if the way I'm checking whether the intersection point is on the box is correct or if I should use some other algorithm.

Thanks.

Answer 1

Code looks fine. Try to find this particular ray and debug it.

Answer 2

EDIT: Ignore this answer (see comments below where I am quite convincingly shown the error of my ways).

You are testing whether the point is inside the volume, but the point is on the periphery of that volume, so you may find that it is an "infinitesimal" distance outside the volume. Try growing the box by some small epsilon, e.g.:

double epsilon = 1e-10; // Depends the scale of things in your code.
double min_x = min(corner1.X(), corner2.X()) - epsilon;
double max_x = max(corner1.X(), corner2.X()) + epsilon;
double min_y = min(corner1.Y(), corner2.Y()) - epsilon;
...

Technically, the correct way to compare almost-equal numbers is to cast their bit representations to ints and compare the integers for some small offset, e.g., in C:

#define EPSILON 10 /* some small int; tune to suit */
int almostequal(double a, double b) {
    return llabs(*(long long*)&a - *(long long*)&b) < EPSILON;
}

Of course, this isn't so easy in C#, but perhaps unsafe semantics can achieve the same effect. (Thanks to @Novox for his comments, which lead me to a nice page explaining this technique in detail.)

Answer 3

Could it be that that ray ends up passing exactly through the edge of the box? Floating point roundoff errors might cause it to be missed by both the right and the back face.

Answer 4

Increasing things by epsilon is not actually a great way to do this, as you now have a border of size epsilon at the edge of your box through which rays can pass. So you'll get rid of this (relatively common) weird set of errors, and end up with another (rarer) set of weird errors.

I assume that you're already envisioning that your ray is traveling at some speed along its vector and find the time of intersection with each plane. So, for example, if you are intersecting the plane at x=x0, and your ray is going in direction (rx,ry,rz) from (0,0,0), then the time of intersection is t = x0/rx. If t is negative, ignore it--you're going the other way. If t is zero, you have to decide how to handle that special case--if you're in a plane already, do you bounce off it, or pass through it? You may also want to handle rx==0 as a special case (so that you can hit the edge of the box).

Anyway, now you have exactly the coordinates where you struck that plane: they are (t*rx , t*ry , t*rz). Now you can just read off whether t*ry and t*rz are within the rectangle they need to be in (i.e. between the min and max for the cube along those axes). You don't test the x coordinate because you already know that you hit it Again, you have to decide whether/how to handle hitting corners as a special case. Furthermore, now you can order your collisions with the various surfaces by time and pick the first one as your collision point.

This allows you to compute, without resorting to arbitrary epsilon-factors, whether and where your ray intersects your cube, to the accuracy possible with floating point arithmetic.

So you just need three functions like the one you've already got: one for testing whether you hit within yz assuming you hit x, and the corresponding ones for xz and xy assuming that you hit y and z respectively.

---

Edit: code added to (verbosely) show how to do the tests differently for each axis:

#define X_FACE 0
#define Y_FACE 1
#define Z_FACE 2
#define MAX_FACE 4

// true if we hit a box face, false otherwise
bool hit_face(double uhit,double vhit,
                 double umin,double umax,double vmin,double vmax)
{
  return (umin <= uhit && uhit <= umax && vmin <= vhit && vhit <= vmax);
}

// 0.0 if we missed, the time of impact otherwise
double hit_box(double rx,double ry, double rz,
                double min_x,double min_y,double min_z,
                double max_x,double max_y,double max_z)
{
  double times[6];
  bool hits[6];
  int faces[6];
  double t;
  if (rx==0) { times[0] = times[1] = 0.0; }
  else {
    t = min_x/rx;
    times[0] = t; faces[0] = X_FACE; 
    hits[0] = hit_box(t*ry , t*rz , min_y , max_y , min_z , max_z);
    t = max_x/rx;
    times[1] = t; faces[1] = X_FACE + MAX_FACE;
    hits[1] = hit_box(t*ry , t*rz , min_y , max_y , min_z , max_z);
  }
  if (ry==0) { times[2] = times[3] = 0.0; }
  else {
    t = min_y/ry;
    times[2] = t; faces[2] = Y_FACE;
    hits[2] = hit_box(t*rx , t*rz , min_x , max_x , min_z , max_z);
    t = max_y/ry;
    times[3] = t; faces[3] = Y_FACE + MAX_FACE;
    hits[3] = hit_box(t*rx , t*rz , min_x , max_x , min_z , max_z);
  }
  if (rz==0) { times[4] = times[5] = 0.0; }
  else {
    t = min_z/rz;
    times[4] = t; faces[4] = Z_FACE;
    hits[4] = hit_box(t*rx , t*ry , min_x , max_x , min_y , max_y);
    t = max_z/rz;
    times[5] = t; faces[5] = Z_FACE + MAX_FACE;
    hits[5] = hit_box(t*rx , t*ry , min_x , max_x , min_y , max_y);
  }
  int first = 6;
  t = 0.0;
  for (int i=0 ; i<6 ; i++) {
    if (times[i] > 0.0 && (times[i]<t || t==0.0)) {
      first = i;
      t = times[i];
    }
  }
  if (first>5) return 0.0;  // Found nothing
  else return times[first];  // Probably want hits[first] and faces[first] also....
}

(I just typed this, didn't compile it, so beware of bugs.) (Edit: just corrected an i -> first.)

Anyway, the point is that you treat the three directions separately, and test to see whether the impact has occurred within the right box in (u,v) coordinates, where (u,v) are either (x,y), (x,z), or (y,z) depending on which plane you hit.

Answer 5

PointOnBoxFace should be a two-dimensional check instead of three-dimensional. For example, if you're testing against the z = z_min plane, then you should only need to compare x and y to their respective boundaries. You've already figured out that z coordinate is correct. Floating point precision is likely tripping you up as you "re-check" the third coordinate.

For example, if z_min is 1.0, you first test against that plane. You find an intersection point of (x, y, 0.999999999). Now, even though x and y are within bounds, the z isn't quite right.