Thinking Inside the Box

There's a non-trivial algorithm in Computer Graphics, the so-called Kay-Kajiya test, that is actually trivial to derive by carefully thinking about the problem it solves. A axis-aligned box is the 3-dimensional generalization of a closed interval. Intuitively, it's a rectangular box whose six sides are parallel to the three axes; formally, given points P = (p_x, p_y, p_z) and Q = (q_x, q_y, q_z), it is the set of points B = { X = (x, y, z) | p_x ≤ x ≤ q_x ∧ p_y ≤ y ≤ q_y ∧ p_z ≤ z ≤ q_z }. As you can see, it is the "product" of three real intervals; as such, it can be generalized for any number of dimensions, but this would require an uncountable set of subscripts, which is not really practical. Pragmatically, I'll resort to naturally-indexed sets of coordinates—arrays:

type point = float array
type vect  = float array

The algorithm is in fact generic in the number of dimensions:

let dim = 3

With that, a box is just the left-bottom-front ("low") and the right-top-back ("high") points:

type bbox  = { l : point; h : point; }

A parametric line is given by a direction vector and a point belonging to the line:

type line  = { p : point; v : vect ; }

Now, the line l intersects the box B if and only if they have a point in common. This can have a set-based expression l ∩ B ≠ ∅, or a predicate-based existential expression ⟨∃X :: X ∈ l ∧ X ∈ B⟩ which is more amenable to formal manipulation:

   ⟨∃X :: X ∈ l ∧ X ∈ B⟩
≡ { parametricity of line, definition of box }
   ⟨∃X :: ⟨∃t : t ∈ R : X = l.p + t⋅l.v ⟩ ∧ ⟨⋀i : 0 ≤ i < D : B.l.i ≤ X.i ≤ B.h.i ⟩⟩
≡ { distributivity of ∧ over ∃ }
   ⟨∃X :: ⟨∃t : t ∈ R : X = l.p + t⋅l.v ∧ ⟨⋀i : 0 ≤ i < D : B.l.i ≤ X.i ≤ B.h.i ⟩⟩⟩
≡ { interchange }
   ⟨∃t : t ∈ R : ⟨∃X :: X = l.p + t⋅l.v ∧ ⟨⋀i : 0 ≤ i < D : B.l.i ≤ X.i ≤ B.h.i ⟩⟩⟩
≡ { one-point rule, defn. vector, point }
   ⟨∃t : t ∈ R : ⟨⋀i : 0 ≤ i < D : B.l.i ≤ l.p.i + t⋅l.v.i ≤ B.h.i ⟩⟩
≡ { algebra }
   ⟨∃t : t ∈ R : ⟨⋀i : 0 ≤ i < D : (B.l.i - l.p.i)/l.v.i ≤ t ≤ (B.h.i - l.p.i)/l.v.i⟩⟩

This is a set of D simultaneous linear inequalities. Any solution t must then satisfy the most restrictive of the inequalities, that is, their extreme bounds:

⇐ { upper and lower bounds }
   ⟨∃t : t ∈ R : ⟨max i : 0 ≤ i < D : (B.l.i - l.p.i)/l.v.i⟩ ≤ t ≤ ⟨min i : 0 ≤ i < D : (B.h.i - l.p.i)/l.v.i⟩⟩
≡ { real interval }
   ⟨max i : 0 ≤ i < D : (B.l.i - l.p.i)/l.v.i⟩ ≤ ⟨min i : 0 ≤ i < D : (B.h.i - l.p.i)/l.v.i⟩

There are a number of complications, though. The use of the parametric definition of the line makes implicit use of the fact that it is not really a line but an oriented ray; this means that if the box is "behind" the ray it won't be hit by the test. The solution for this is to sort the extremes to be tested. Also, if a vector component is zero the test must explicitly decide (that is, test to avoid dividing by zero) if a ray belonging to the side but not to the volume is actually inside or not the box. With these caveats, the following code is self-explanatory:

let intersects l b =
  let rec loop i hither thither =
    if i == dim then hither <= thither else
    let v, p, l, u = l.v.(i), l.p.(i), b.l.(i), b.h.(i) in
    if v <> 0. then
      let t0, t1 =
        let t0, t1 = (l -. p) /. v, (u -. p) /. v in
        if t0 < t1 then t0, t1 else t1, t0 in
      let hither  = max hither  t0
      and thither = min thither t1 in
      if hither > thither then false else
      loop (i + 1) hither thither
    else if l > p || p > u then false else
    loop (i + 1) hither thither
  in loop 0 neg_infinity infinity

The usual explanation of the algorithm is both motivated by and restricted to an appeal to geometric handwaving. A simple manipulation of the definition makes the derivation automatic and forced.