How to find the closest intersection between a line and several spheres without using a square root?

Question

I am doing some ray tracing and I have the equation of a line and several spheres with their coordinates and radius. I want to find the closest intersection with a sphere, which I know I can do with the distance formula on wikipedia. However, that requires the square root, which is computationally expensive, so I would like to know if there's some way to figure out which sphere intersects at the closest point without using square roots. Note that this is relative to each other so hopefully there's a way to do this without exact precision

Answer 1

If you have one or more spheres centered at $\vec{c}_i$ with radii $r_i$, and a ray in direction $\hat{n}$ from origin ($\hat{n}$ being an unit vector, $\hat{n}\cdot\hat{n} = 1$), then the distances at which the ray intersects (the side closest to origin of) each sphere $i$ is $$d_i = \hat{n} \cdot \vec{c}_i - \sqrt{\left(\hat{n} \cdot \vec{c}_i\right)^2 + r_i^2 - \vec{c}_i \cdot \vec{c}_i} \tag{1a}\label{NA1a}$$ if and only if an intersection does occur, $$\left(\hat{n} \cdot \vec{c}_i\right)^2 + r_i^2 - \vec{c}_i \cdot \vec{c}_i \ge 0 \tag{1b}\label{NA1b}$$ In the equal-to-zero case in $\eqref{NA1b}$, the ray is tangent to sphere $i$.

There is no way of omitting the square root in $\eqref{NA1a}$; even squaring $d_i$ still leaves you with a square root per sphere.

---

Note that since origin is at the camera focal point, you can precalculate a scalar constant $$C_i = \vec{c}_i \cdot \vec{c}_i - r_i^2$$ for each sphere; it is the same regardless of ray direction $\hat{n}$.

Then, for each ray $\hat{n}$ and for each sphere $i$, calculate $$D_i = \hat{n} \cdot \vec{c}_i, \quad S_i = D_i^2$$ If $S_i \lt C_i$, the ray does not intersect sphere $i$. Otherwise, the distance to the intersection is $$d_i = D_i - \sqrt{S_i - C_i}$$

---

In a computer program, when raycasting many rays towards many spheres, the main bottleneck is cache behaviour (and not the square root operation, if it is very slow compared to multiplication or addition).

A very effective optimization is to sort the spheres in increasing values of the minimum distance from sphere $i$ to origin, $L_i = \sqrt{\vec{c}_i \cdot \vec{c}_i} - r_i$.

For best cache behaviour, $L_i$ and $D_i$ would be in separate arrays ($S_i$ is better to calculate than to cache), and $r_i$ and $\vec{c}_i$ in a third array (where they can be used to calculate $D_i$ for each new ray direction, or $L_i$ after moving origin with all arrays sorted by $L_i$ afterwards). Note that these are in the same order regardless of ray direction $\hat{n}$.

This way, the program only needs to examine the spheres in that order, examining memory in sequential order, until the smallest intersection distance thus far is less than $L_i$ -- all intersections with spheres with a larger $L_i$ are further away.

Obviously, subdividing the space into cells (cuboids, cones, pyramids, or frustums) can speed up things even more, because the set of spheres needed to examine for each ray is smaller. For spatial subdivision, you'll then need to decide whether you'll duplicate spheres that intersect more than one cell, or whether each ray inspects more than one cell. This gets quite complicated quite fast, especially if the spheres have widely varying radii.