Mapping reflections in 3D

Question

I am working on a program in python that takes a point $\begin{pmatrix} x_i \\y_i\\z_i \end{pmatrix} $, generates a displacement vector $\begin{pmatrix} x_d \\y_d \\z_d\end{pmatrix} $ and adds it to the initial pont to determine the next point$\begin{pmatrix} x_f \\y_f \\z_f\end{pmatrix} = \begin{pmatrix} x_i \\y_i\\z_i \end{pmatrix}+\begin{pmatrix} x_d \\y_d\\z_d \end{pmatrix} $

My trajectory is confined to a rectangular box $\begin{bmatrix} x_0, x_1 \end{bmatrix} \times \begin{bmatrix} y_0, y_1 \end{bmatrix} \times \begin{bmatrix} z_0, z_1 \end{bmatrix} $ So when the trajectory hits a wall it bounces off. My question is how can I construct a function that takes the initial point and the displacement vector (and the boundaries of the box) and returns the point where my trajectory will end up after bouncing off said wall?

Answer 1

There is a great way to do this, that works whenever the confinement volume is convex; and boxes and cuboids are definitely convex.

Define each of your walls as a plane (or half-space) using its unit normal vector $\hat{n}$ and its signed distance $d$ from origin. Have $\hat{n}$ point in the direction of the bounce.

For an axis-aligned box $(x_1, y_1, z_1) - (x_2, y_2, z_2)$, where $x_2 \gt x_1$, $y_2 \gt y_1$, $z_2 \gt z_1$, those planes are $$\begin{array}{ll} \hat{n}_1 = (1, 0, 0), & d_1 = x_1 \\ \hat{n}_2 = (-1, 0, 0), & d_2 = -x_2 \\ \hat{n}_3 = (0, 1, 0), & d_3 = y_1 \\ \hat{n}_4 = (0, -1, 0), & d_4 = -y_2 \\ \hat{n}_5 = (0, 0, 1), & d_5 = z_1 \\ \hat{n}_6 = (0, 0, -1), & d_6 = -z_2 \\ \end{array}$$

Point $\vec{p}$ signed distance $s_i$ from wall $i$ is $$s_i = \vec{p} \cdot \hat{n}_i - d_i$$ This is positive whenever the point is within the confined volume, and negative otherwise. Note that if your particle is a sphere of radius $r$, just use $$s_i = \vec{p} \cdot \hat{n}_i - d - r$$

To calculate the step for a point currently at $(x_c, y_c, z_c)$ with velocity $(x_\Delta, y_\Delta, z_\Delta)$, we start by calculating the next location for this point at $(x_n, y_n, z_n) = (x_c + t_\Delta x_\Delta, y_c + t_\Delta y_\Delta, z_c + t_\Delta z_\Delta)$. Note that if you use fixed time steps, then you can use $t_\Delta = 1$. However, if the particle velocities are not capped to the shortest distance between opposing walls, the particle may escape. One way to avoid this is to use (for each time step, for all particles in that time step) $$t_\Delta = \min\left(\tau_\max, ~ \frac{\Delta_\max}{\max\left(\lvert x_{\Delta, i} \rvert, \lvert y_{\Delta, i} \rvert, \lvert z_{\Delta, i} \rvert\right)}\right) \text{over all particles } i$$ which caps the maximum time step at $\tau_\max$, but will yield a smaller time step when necessary so that no particle is displaced by more than $\Delta_\max$ along an axis (so step distance is at most $\sqrt{3} \, \Delta_\max \approx 1.733 \, \Delta_\max$). If you choose $L$ based on the size of your confinement box, say $L = 0.5 \min(x_2 - x_1, y_2 - y_1, z_2 - z_1)$, then no matter what the velocities of your particles, they won't travel so fast they can escape (by colliding with the same wall twice in a time step - as this implementation can only bounce once from each wall per time step).

Then, we enter a loop that fixes $(x_n, y_n, z_n)$ for any bounces. This is because a single time step may involve more than one bounce, say near a corner. The loop maintains a list of planes we have already bounced from once (or say a flag per plane), and those planes are ignored from the following iterations.

Within the loop, you calculate $s_i$ for all planes we haven't yet bounced from in this timestep for this particle, and find the smallest signed distance, i.e $k$ so that $s_k = \min(s_i)$. If $s_k \gt 0$, there are no bounces, and we can exit the loop for the current particle, setting its location to $(x_n, y_n, z_n)$. Otherwise, we bounce the particle from wall $k$, and mark that wall for "bounced" so that we won't consider it in any following iterations (for this particle in this time step).

Calculate $c_k$, the signed distance of the particles current position $(x_c, y_c, z_c)$ to plane $k$. Note that $c_k \gt 0$, or it has already escaped; and $s_k \lt 0$ because we already checked. Because the trajectory during the time step is linear, we know that the wall divides the trajectory as $c_k:-s_k$, i.e. it encounters the wall at fraction $f$ of the distance (or time step), with the fraction $b$ reflected from the wall, $$f = \frac{c_k}{c_k - s_k}, \quad b = \frac{-s_k}{c_k - s_k}, \quad f + b = 1$$ Calculate the position $(x_w, y_w, z_w)$ exactly at the wall contact point, and the vector corresponding to the remaining displacement that we need to reflect $(x_d, y_d, z_d)$: $$\left\lbrace\begin{aligned} x_w &= b x_c + f x_n \\ y_w &= b y_c + f y_n \\ z_w &= b z_c + f z_n \\ \end{aligned}\right., \quad \left\lbrace\begin{aligned} x_d &= x_n - x_w \\ y_d &= y_n - y_w \\ z_d &= z_n - z_w \\ \end{aligned}\right.$$ Note that this displacement vector is the one we need to reflect with respect to plane $k$, to find where the particle (tentatively) should end up at the end of current time step. This is simple: to reflect vector $\vec{v}$ with respect to a plane having unit normal $\hat{n}$ (unit meaning $\lVert\hat{n}\lVert = 1$),

$$\vec{v}_r = \vec{v} - 2\left(\hat{n}\cdot\vec{v}\right)\hat{n}$$ Note that we need to reflect both the remaining displacement as well as the velocity vector for this particle. If we use $(x_k, y_k, z_k)$ for the normal $\hat{n}_k$ of this plane, we calculate $$\left\lbrace\begin{aligned} g_1 &= 2 ( x_k x_d + y_k y_d + z_k z_d ) \\ g_2 &= 2 ( x_k x_\Delta + y_k y_\Delta + z_k z_\Delta ) \\ \end{aligned}\right.$$ and then $$\left\lbrace\begin{aligned} x_n &= x_w + x_d - g_1 x_k \\ y_n &= y_w + y_d - g_1 y_k \\ z_n &= z_w + z_d - g_1 z_k \\ \end{aligned}\right., \quad \left\lbrace\begin{aligned} x_c &= x_w \\ y_c &= y_w \\ z_c &= z_w \\ \end{aligned}\right., \quad \left\lbrace\begin{aligned} x_\Delta &= x_\Delta - g_2 x_k \\ y_\Delta &= y_\Delta - g_2 y_k \\ z_\Delta &= z_\Delta - g_2 z_k \\ \end{aligned}\right.$$

At this point, we have reflected the particle from plane $k$, but it might have reflected from the other planes as well. That's why we need to iterate the loop for this particle during this time step, until there are either no more walls to check (in which case we could just reset to checking all walls again!), or the smallest signed distance we find is positive, $s_k \gt 0$.

Also note that during this iteration, we essentially split the time step into smaller parts, at each collision with a wall.

I wonder if I should restate the "rules" in basic Vector Algebra notation, with links to some details (like Wikipedia pages)? This answer is already rather long. Let me know in a comment, if you think it is necessary or useful.

Answer 2

Although I posted a description of the algorithm, the mathematical formulae used are easily lost in the wall of text. So, I decided to write a separate answer to list the basic vector algebra operations one needs here.

Define the walls as planes (properly, halfspaces), using an unit normal $\hat{n}$ ($\lVert\hat{n}\rVert = 1$) and signed distance from origin $d$, with the unit normal pointing away from the wall.

---

For an axis aligned box enclosure $(x_\min, y_\min, z_\min) - (x_\max, y_\max, z_\max)$ the six walls with unit normals pointing inside the box are $$\begin{array}{r|r} \hat{n} ~ ~ ~ ~ & d ~ ~ ~ ~ \\ \hline (1, 0, 0) & x_\min \\ (-1, 0, 0) & -x_\max \\ (0, 1, 0) & y_\min \\ (0, -1, 0) & -y_\max \\ (0, 0, 1) & z_\min \\ (0, 0, -1) & -z_\max \\ \end{array} \tag{1}\label{G1}$$

---

Distance from sphere of radius $r$ ($r = 0$ for a point) centered at $\vec{p}$, to halfspace $\hat{n}$, $d$, is $$L = \vec{p} \cdot \hat{n} - d - r\tag{2}\label{G2}$$ This will be positive outside the halfspace, and zero when the sphere touches the surface of the halfspace.

---

A sphere of radius $r$ (for a point, $r = 0$) located outside the halfspace, at $\vec{p}$ at time $t=0$, having a velocity vector $\vec{v}$, and therefore its location as a function of time being $$\vec{p}(t) = \vec{p} + t \vec{v}$$ will touch halfspace $\hat{n}$, $d$, at time $t \gt 0$ if and only if $$\hat{n} \cdot \vec{v} \lt 0$$ where $$t = \frac{d + r - \hat{n} \cdot \vec{p}}{\hat{n} \cdot \vec{v}} \tag{3}\label{G3}$$ When dealing with multiple spheres and/or walls, solve all positive $t$, and handle them in order of increasing $t$. In a simulation, this lets you skip ahead in the simulation to the next bounce event, without having to calculate any intermediate steps. If you store $\vec{p}_0$, $\vec{v}$, $t_0$ (so $\vec{p}(t_0) = \vec{p}_0$), and the time $t_1$ for its next bounce (and the wall it bounces from) for each sphere, you only need to update that one sphere for each bounce. The position of the sphere at time $t_0 \le t \le t_1$ is then $$\vec{p}(t) = \vec{p}_0 + (t - t_0) \vec{v}, \quad t_0 \le t \le t_1$$

---

To reflect vector $\vec{v}$ (outside the halfspace) with respect to halfspace $\hat{n}$, $d$, $$\begin{aligned} \vec{v}_r &= \vec{v} - 2 \left( \vec{v} \cdot \hat{n} \right) \hat{n} \\ ~ &= \vec{v} - \hat{n} \bigl( 2 \hat{n} \cdot \vec{v} \bigr) \\ \end{aligned} \tag{4}\label{G4}$$ Note that the part in parentheses evaluates to a scalar (a real number, not a vector).

---

A sphere of radius $r$ ($r = 0$ for a point) traveling between $\vec{p}_o$ outside the halfspace and $\vec{p}_i$ inside the halfspace $\hat{n}$, $d$, will touch the halfspace at $\vec{p}_b$, $$\begin{aligned} L_o &= \vec{p}_o \cdot \hat{n} - d - r \\ L_i &= \vec{p}_i \cdot \hat{n} - d - r \\ \vec{p}_b &= \displaystyle \vec{p}_o + \frac{L_o}{L_o - L_i} \left( \vec{p}_i - \vec{p}_o \right) = \displaystyle \frac{L_o \vec{p}_i - L_i \vec{p}_o}{L_o - L_i} \\ \end{aligned} \tag{5a} \label{G5a}$$ noting that $L_o \gt 0$ because $\vec{p}_o$ is outside the halfspace, and $L_i \lt 0$ because $\vec{p}_i$ is inside the halfspace.

The "untraveled" part that is reflected off the wall is $$\begin{aligned} \vec{p}_r &= \frac{-L_i}{L_o - L_i} \left( \vec{p}_i - \vec{p}_o \right) = \vec{p}_i - \vec{p}_b \\ \end{aligned}$$ where $\vec{p}_b$ is the above calculated position where the sphere touches the halfspace.

In other words, if the sphere does not bounce off any other walls, it will end up at $\vec{p}_e$, $$\begin{aligned} \vec{p}_e &= \vec{p}_b - \hat{n} \bigl( 2 \hat{n} \cdot ( \vec{p}_i - \vec{p}_b ) \bigr) \\ \end{aligned} \tag{5b}\label{G5b}$$ where the part in larger parentheses evaluates to a scalar (a real number, not a vector).