Finding the tightest axis-aligned bounding box for a sphere is trivial: the box extends from the center by the radius in all dimensions. But once the sphere is transformed, finding the minimal bounding box becomes trickier. Rotating a sphere, for example, shouldn't change its bounding box, but naïvely rotating the bounding box will expand it unnecessarily. Luckily there's a trick to computing minimal bounding boxes by representing the transformed sphere as a quadric surface.

A description of the technique can be found on Yining Karl Li's blog here, and in this Stack Overflow answer. Unfortunately, both descriptions gloss over the math, as well as some details that give rise to a simple and efficient implementation. The post on Li's blog also contains a small bug. In contrast, this post examines the technique in full detail. It relies on the projective representation of a quadric.

A sphere of radius 1, centred at the origin, can be defined as the set of points where \(x^2 + y^2 + z^2 = 1\). For the homogeneous point \(\vec{p} = \begin{bmatrix}x & y & z & 1\end{bmatrix}^\mathrm{T}\), the same sphere can be defined by

$$

\vec{p}^\mathrm{T} \, \mathbf{S} \, \vec{p} = 0,

$$ where $$

\mathbf{S} = \begin{bmatrix}

1 & 0 & 0 & 0 \\

0 & 1 & 0 & 0 \\

0 & 0 & 1 & 0 \\

0 & 0 & 0 & -1

\end{bmatrix}.

$$

If an arbitrary affine transformation \(\mathbf{M}\) is applied to the sphere, the equation becomes

$$

\begin{align*}

(\mathbf{M}^{-1}\,\vec{p})^\mathrm{T} \, \mathbf{S} \, (\mathbf{M}^{-1}\,\vec{p}) & = \vec{p}^\mathrm{T} \, (\mathbf{M}^{-\mathrm{T}}\,\mathbf{S}\,\mathbf{M}^{-1}) \, \vec{p} \\

{} & = 0.

\end{align*}

$$

We denote this transformed quadric matrix by \(\mathbf{Q} = \mathbf{M}^{-\mathrm{T}}\,\mathbf{S}\,\mathbf{M}^{-1}\). Note that \(\mathbf{Q}\) is symmetric.

## The Dual Form

The minimal bounding box for this surface will be tangent to it, so we need to find an expression for its tangent planes. In homogeneous coordinates, a plane can be represented by a row vector \(\vec{u}^\mathrm{T}\) such that \(\vec{u}\cdot\vec{p} = 0\). For a point \(\vec{p}\) on the surface, the plane \(\vec{u}^\mathrm{T} = \vec{p}^\mathrm{T} \, \mathbf{Q}\) touches the quadric at \(\vec{p}\), since

$$

\begin{align*}

\vec{u}\cdot\vec{p} & = \vec{u}^\mathrm{T} \, \vec{p} \\

{} & = \vec{p}^\mathrm{T} \, \mathbf{Q} \, \vec{p} \\

{} & = 0.

\end{align*}

$$

Furthermore, this must be the unique intersection point. To see why, consider a point \(\vec{q} = \vec{p} + r\,\hat{r}\) that falls on both the plane and the quadric.

$$

\begin{align*}

\vec{q}^\mathrm{T} \, \mathbf{Q} \, \vec{q} & = (\vec{p} + r\,\hat{r})^\mathrm{T} \, \mathbf{Q} \, \vec{q} \\

{} & = \vec{u}\cdot\vec{p} + r\,\hat{r}^\mathrm{T} \, \mathbf{Q} \, \vec{q} \\

{} & = r \, (\vec{q}^\mathrm{T} \, \mathbf{Q} \, \hat{r})^\mathrm{T} \\

{} & = r^2 \, (\hat{r}^\mathrm{T} \, \mathbf{Q} \, \hat{r}) \\

{} & = 0.

\end{align*}

$$

If \(\hat{r}^\mathrm{T} \, \mathbf{Q} \, \hat{r}\) is zero, then the whole line \(\vec{p} + t\,\hat{r}\) must lie on the quadric, which is impossible since its surface is curved. Therefore \(r = 0\) and \(\vec{q} = \vec{p}\). Since the plane has exactly one point in common with the quadric, \(\vec{u}^\mathrm{T}\) must be the tangent plane at that point. This gives a nice characterization of the tangent planes:

$$

\begin{align*}

\vec{u}^\mathrm{T} \, \mathbf{Q}^{-1} \, \vec{u} & = \vec{p}^\mathrm{T} \, \mathbf{Q} \, \mathbf{Q}^{-1} \, \mathbf{Q} \, \vec{p} \\

{} & = \vec{p}^\mathrm{T} \, \mathbf{Q} \, \vec{p} \\

{} & = 0.

\end{align*}

$$

## Solving for the Planes

Let \(\mathbf{R} = \mathbf{Q}^{-1} = \mathbf{M} \, \mathbf{S}^{-1} \, \mathbf{M}^\mathrm{T}\). Like \(\mathbf{Q}\), \(\mathbf{R}\) is symmetric. We want to find axis-aligned planes \(\vec{u}^\mathrm{T}\) such that \(\vec{u}^\mathrm{T} \, \mathbf{R} \, \vec{u} = 0\). The plane perpendicular to \(\hat{x}\) looks like \(\vec{u}^\mathrm{T} = \begin{bmatrix}1 & 0 & 0 & -x\end{bmatrix}\), so we have

$$

\begin{align*}

\vec{u}^\mathrm{T} \, \mathbf{R} \, \vec{u} & =

\begin{bmatrix}

1 & 0 & 0 & -x

\end{bmatrix}

\,

\mathbf{R}

\,

\begin{bmatrix}

1 \\

0 \\

0 \\

-x \\

\end{bmatrix} \\

{} & =

\begin{bmatrix}

1 & 0 & 0 & -x

\end{bmatrix}

\,

\begin{bmatrix}

\mathbf{R}_{1,1} - \mathbf{R}_{1,4}\,x \\

\mathbf{R}_{2,1} - \mathbf{R}_{2,4}\,x \\

\mathbf{R}_{3,1} - \mathbf{R}_{3,4}\,x \\

\mathbf{R}_{4,1} - \mathbf{R}_{4,4}\,x \\

\end{bmatrix} \\

{} & = \mathbf{R}_{4,4}\,x^2 - (\mathbf{R}_{1,4} + \mathbf{R}_{4,1})\,x + \mathbf{R}_{1,1} \\

{} & = \mathbf{R}_{4,4}\,x^2 - 2\,\mathbf{R}_{1,4}\,x + \mathbf{R}_{1,1} \\

{} & = 0. \\

x & = \frac{\mathbf{R}_{1,4} \pm \sqrt{\mathbf{R}_{1,4}^2 - \mathbf{R}_{1,1}\,\mathbf{R}_{4,4}}}{\mathbf{R}_{4,4}}. \\

\end{align*}

$$

Similarly, the \(\hat{y}\) and \(\hat{z}\) planes are given by

$$

\begin{align*}

y & = \frac{\mathbf{R}_{2,4} \pm \sqrt{\mathbf{R}_{2,4}^2 - \mathbf{R}_{2,2}\,\mathbf{R}_{4,4}}}{\mathbf{R}_{4,4}} \\

z & = \frac{\mathbf{R}_{3,4} \pm \sqrt{\mathbf{R}_{3,4}^2 - \mathbf{R}_{3,3}\,\mathbf{R}_{4,4}}}{\mathbf{R}_{4,4}}. \\

\end{align*}

$$

This is where the above linked posts stop, but a lot more simplification can be done.

## Implementation Details

Several details of the problem can make computing the planes more efficient. The first is that \(\mathbf{S}\) is involutory, meaning \(\mathbf{S}^{-1} = \mathbf{S}\). This means that the product \(\mathbf{M}\,\mathbf{S}^{-1}\) can be computed implicitly, it is simply \(\mathbf{M}\) with its last column negated. The last column of \(\mathbf{R} = \mathbf{M} \, \mathbf{S}^{-1} \, \mathbf{M}^\mathrm{T}\) is the same, because the last column of \(\mathbf{M}^\mathrm{T}\) is \(\begin{bmatrix}0 & 0 & 0 & 1\end{bmatrix}^\mathrm{T}\). In particular, \(\mathbf{R}_{4,4} = -1\).

Not all values of \(\mathbf{R}\) are used; in fact, only values from the last column and the diagonal appear in the formulae. We know the last column implicitly, and the diagonal has the formula

$$

\mathbf{R}_{i,i} = \left(\sum_{j=1}^3 \mathbf{M}_{i,j}^2\right) - \mathbf{M}_{i,4}^2.

$$

Plugging this identity into the plane equations simplifies them greatly:

$$

\begin{align*}

x & = \mathbf{M}_{1,4} \pm \sqrt{\mathbf{M}_{1,1}^2 + \mathbf{M}_{1,2}^2 + \mathbf{M}_{1,3}^2} \\

y & = \mathbf{M}_{2,4} \pm \sqrt{\mathbf{M}_{2,1}^2 + \mathbf{M}_{2,2}^2 + \mathbf{M}_{2,3}^2} \\

z & = \mathbf{M}_{3,4} \pm \sqrt{\mathbf{M}_{3,1}^2 + \mathbf{M}_{3,2}^2 + \mathbf{M}_{3,3}^2}.

\end{align*}

$$

This makes at least some intuitive sense, since the fourth column of \(\mathbf{M}\) determines the sphere's translation. See here for an implementation in my ray tracer Dimension.