Let a rectangular cuboid of size $(L_x, L_y, L_z)$ be located with one corner at the origin and aligned with the $(x,y,z)$ axes. Let $\overrightarrow{sr}$ be a vector from point $s$ to point $r$. $s$ is known to be outside the cube, while $r$ could be inside or outside the cube but neither is on the faces. The goal is to find if the line goes through the cube and which side it intersects first. If the line is on one of the planes, it falls under the definition of not going through the cube, namely, we are only interested on a single point crossing.
This can be done easily enough by parameterizing the line by $s+(r-s)t=p(t)$. The trivial calculation can be time-consuming. We need to intersect the line with 6 planes, constraint the results to the limits of the square on each plane, and finally determine the first encounter.
The thing is, due to the convenient location of the cube, this computation may contain many shortcuts. for example, if we define 6 normals directed outwards of the rectangular we can eliminate the last part by noticing the sign of the dot product between the line and each of the normals. A positive dot product indicates it is a first encounter while a negative one means it is not. Due to the relatively easy 6 normals, this dot multiplication is reduced to analyzing the sign of a single component in the direction vector of the line.
I wish to implement this in a program with a minimum amount of calculations. I am looking for the shortest, closed solution to such a problem under these assumptions.
I am looking for mathematical assumptions\trics\accelerations that can shorten the calculation and not programming optimization techniques.
Best Answer
Note: A previous version of this answer considered only the planes of the box sides, and not whether the intersection points were within the corresponding box face. This version includes a complete Python example/verification of the corrected approach.
Let's first examine the intersection of the line segment from $\vec{r} = (x_0, y_0, z_0)$ to $\vec{s} = (x_1, y_1, z_1)$ with the planes of the axis-aligned box with one vertex at origin, and the diagonally opposite vertex at $(L_x, L_y, L_z)$.
Parametrising the line as a vector-valued function using $0 \le t \le 1$, we have $$\vec{p}(t) = (1 - t) \vec{r} + t \vec{s} = \vec{r} + t (\vec{s} - \vec{r})$$ and the intersections with the six planes of the box faces are at $$\begin{array}{c|c|c|c|c} \text{Face} & t & x & y & z \\ \hline 1 & \frac{x_0}{x_0 - x_1} & 0 & \frac{x_0 y_1 - x_1 y_0}{x_0-x_1} & \frac{x_0 z_1 - x_1 z_0}{x_0 - x_1} \\ \hline 2 & \frac{x_0 - L_x}{x_0 - x_1} & L_x & \frac{(L_x - x_1) y_0 - (L_x - x_0) y_1}{x_0 - x_1} & \frac{(L_x - x_1) z_0 - (L_x - x_0) z_1}{x_0 - x_1} \\ \hline 3 & \frac{y_0}{y_0 - y_1} & \frac{x_1 y_0 - x_0 y_1}{y_0 - y_1} & 0 & \frac{y_0 z_1 - y_1 z_0}{y_0 - y_1} \\ \hline 4 & \frac{y_0 - L_y}{y_0 - y_1} & \frac{(L_y - y_1) x_0 - (L_y - y_0) x_1}{y_0 - y_1} & L_y & \frac{(L_y - y_1) z_0 - (L_y - y_0) z_1}{y_0 - y_1} \\ \hline 5 & \frac{z_0}{z_0 - z_1} & \frac{x_1 z_0 - x_0 z_1}{z_0 - z_1} & \frac{y_1 z_0 - y_0 z_1}{z_0 - z_1} & 0 \\ \hline 6 & \frac{z_0 - L_z}{z_0 - z_1} & \frac{ (L_z - z_1) x_0 - (L_z - z_0) x_1}{z_0 - z_1} & \frac{ (L_z - z_1) y_0 - (L_z - z_0) y_1}{z_0 - z_1} & L_z \\ \end{array}$$
To avoid the divisions, we can use $$\tau = t \lvert (x_0 - x_1)(y_0 - y_1)(z_0 - z_1) \rvert$$ Note that if we verify the start and end coordinates cross the desired face in the correct direction, then $0 \le t \le 1$.
If we use thirteen temporary variables, $$\begin{array}{lcl} ~ & \quad & a_{xyz} = \lvert (x_0 - x_1)(y_0 - y_1)(z_0 - z_1) \rvert \\ s_{xy} = ( x_0 - x_1 ) L_y & \quad & a_{xy} = \lvert (x_0 - x_1)(y_0 - y_1) \rvert \\ s_{xz} = ( x_0 - x_1 ) L_z & \quad & a_{xz} = \lvert (x_0 - x_1)(z_0 - z_1) \rvert \\ s_{yx} = ( y_0 - y_1 ) L_x & \quad & a_{yz} = \lvert (y_0 - y_1)(z_0 - z_1) \rvert \\ s_{yz} = ( y_0 - y_1 ) L_z & \quad & c_{xy} = x_1 y_0 - x_0 y_1 \\ s_{zx} = ( z_0 - z_1 ) L_x & \quad & c_{xz} = x_1 z_0 - x_0 z_1 \\ s_{zy} = ( z_0 - z_1 ) L_y & \quad & c_{yz} = y_1 z_0 - z_0 y_1 \\ \end{array}$$ we can calculate the $\tau$ and the required conditions for intersection on that face via $$\begin{array}{c|c|c|c|c} \text{Face} & ~ & ~ & ~ & \tau \\ \hline 1 & x_0 \lt 0 \lt x_1 & 0 \le c_{xy} \le -s_{xy} & 0 \le c_{xz} \le -s_{xz} & -x_0 a_{yz} \\ 2 & x_1 \lt L_x \lt x_0 & s_{yx}-s_{xy} \le c_{xy} \le s_{yx} & s_{zx}-s_{xz} \le c_{xz} \le s_{zx} & (x_0 - L_x) a_{yz} \\ 3 & y_0 \lt 0 \lt y_1 & s_{yx} \le c_{xy} \le 0 & 0 \le c_{yz} \le -s_{yz} & -y_0 a_{xz} \\ 4 & y_1 \lt L_y \lt y_0 & -s_{xy} \le c_{xy} \le s_{yx}-s_{xy} & s_{zy}-z_{yz} \le c_{yz} \le s_{zy} & (y_0 - L_y) a_{xz} \\ 5 & z_0 \lt 0 \lt z_1 & s_{zx} \le c_{xz} \le 0 & s_{zy} \le c_{xy} \lt 0 & -z_0 a_{xy} \\ 6 & z_1 \lt L_z \lt z_0 & -s_{xz} \le c_{xz} \le s_{zx}-s_{xz} & -s_{yz} \le c_{yz} \le s_{zy}-s_{yz} & (z_0 - L_z) a_{xy} \\ \end{array}$$ Note that the first condition cannot be true at the same time for faces $1$ and $2$, faces $3$ and $4$, or faces $5$ and $6$.
Precalculating the variables (using three temporaries) requires $16$ multiplications, $6$ subtractions, $4$ absolute values, and $16$ assignments (including the three temporaries).
In the worst case, we do all six primary requirement checks (12 comparisons), plus four additional comparisons per axis, for a worst case total of $24$ comparisons; also requiring $3$ multiplications, $9$ subtractions, and $4$ negations.
To find the first intersecting face, we need an additional $3$ comparisons and up to $7$ assignments.
Thus, the worst-case overall cost is $19$ multiplications, $15$ subtractions, $8$ negations or absolute values, and up to $27$ comparisons (but only $9$ of those are conditional jumps).
Here is a Python3 implementation of the above approach:
To verify, append
where
verify_box_ray()uses the simpler notation in this answer to ease verification of its correctness.If we form an integer-valued function $N(\vec{p})$ via $$\begin{aligned} c &= 1 \quad \text{ if } x \lt 0 \\ ~ &+ 2 \quad \text{ if } x \gt L_x \\ ~ &+ 3 \quad \text{ if } y \lt 0 \\ ~ &+ 6 \quad \text{ if } y \gt L_y \\ ~ &+ 9 \quad \text{ if } z \lt 0 \\ ~ &+ 18 \quad \text{ if } z \gt L_z \\ \end{aligned}$$ we can classify both $\vec{r}$ and $\vec{s}$ into one of $27$ classes; for a total of 729 cases. A large number of these cases yield an immediate answer, but some need one, two, or three face tests to yield a definitive answer.
Essentially, this cell-based approach only eliminates at most six comparisons (and six conditional jumps), so this is unlikely to be more efficient than the above.
In some cases, we may want to change to a coordinate system where we are looking for the intersections between integer coordinate planes; or, 3D lattice wall intersections.
For this particular box, we only need to divide all $x$ coordinates by $L_x$, all $y$ coordinates by $L_y$, and all $z$ coordinates by $L_z$.
(Note that this is a linear transform that preserves $t$. Also note that with floating-point numbers, multiplying by a reciprocal does not produce identical results. True division yields more precise results. Multiplying by a reciprocal is less precise, because the reciprocal is rounded to within the floating point range and precision used, as a middle step.)
Note that this does not find out which lattice cells the line segment or ray intersects, but which lattice cell walls the line segment or ray intersects. One can be derived from the other rather trivially, but they are not the same thing exactly.
To do this, one must first calculate the lengths in units of $t$ between successive parallel unit lattice cell walls. If we assume $(x_0, y_0, z_0)$ and $(x_1, y_1, z_1)$ are already in the scaled coordinates, then these are just the reciprocals: $$\begin{aligned} \Delta_x &= \displaystyle \frac{1}{x_1 - x_0} \\ \Delta_y &= \displaystyle \frac{1}{y_1 - y_0} \\ \Delta_z &= \displaystyle \frac{1}{z_1 - z_0} \\ \end{aligned}$$ The sign indicates which direction (positive or negative) the ray progresses.
The first $t$ where the ray intersects the integer coordinate planes, are $$\begin{aligned} t_{x \min} &= \begin{cases} \Delta_x (\lfloor x_0 \rfloor - x_0), & \Delta_x \lt 0 \\ \infty, & \Delta_x = 0 \\ \Delta_x (\lceil x_0 \rceil - x_0), & \Delta_x \gt 0 \\ \end{cases} \\ t_{y \min} &= \begin{cases} \Delta_y (\lfloor y_0 \rfloor - y_o), & \Delta_y \lt 0 \\ \infty, & \Delta_y = 0 \\ \Delta_y (\lceil y_0 \rceil - y_0), & \Delta_y \gt 0 \\ \end{cases} \\ t_{z \min} &= \begin{cases} \Delta_z (\lfloor z_0 \rfloor - z_0), & \Delta_z \lt 0 \\ \infty, & \Delta_z = 0 \\ \Delta_z (\lceil z_0 \rceil - z_0), & \Delta_z \gt 0 \\ \end{cases} \\ \end{aligned}$$ Similarly for the two other axes.
The idea in lattice wall sequence finding is for first calculate the first intersections along each axis, say $t_x$, $t_y$, and $t_z$. The very first intersection is the smallest of the three. You use that for $t$, replacing it with the corresponding next intersection: $t_x \gets t_x + \lvert\Delta_x\rvert$, $t_y \gets t_y + \lvert\Delta_y\rvert$, or $t_z \gets t_z + \lvert\Delta_z\rvert$, respectively, depending on which one was the smallest this iteration.
This does mean that each lattice cell or cell wall iteration is just a couple of comparisons and a couple of additions (one for $t$, the other for the lattice cell integer coordinate that changes). For this reason, this is very commonly used in voxel raycasting.
Note that to apply this for OP's stated problem, we might have to iterate through a lot of unit cell walls, if $\lvert x_0 \rvert \gg L_x$, $\lvert x_1 \rvert \gg L_x$, $\lvert y_0 \rvert \gg L_y$, $\lvert y_1 \rvert \gg L_y$, $\lvert z_0 \rvert \gg L_z$, or $\lvert z_1 \rvert \gg L_z$. If it is known that the line segment start and end are always within the origin cell or a neighboring cell, then this approach might yield an even more efficient solution (noting the cost of six divisions), but I doubt it, and therefore did not bother to work it out completely.