Exact Bounding Boxes for Spheres/Ellipsoids

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.

4 thoughts on “Exact Bounding Boxes for Spheres/Ellipsoids”

  1. Hello,

    Very nice simplification. How does M relate to the ellipsoid's eigenvectors and eigenvalues? Or, put another way: I already know the ellipsoid's centre, its radii & eigenvalues (U, 3 × 3) and an eigenvector roation matrix (V, 3 × 3). The centre is M[1,2,3],4, but how can I get the other M elements from what I know? Or another way: is the top left 3×3 of the 4×4 affine transform simply UV (which I also already have)?

    http://en.wikipedia.org/wiki/Ellipsoid#Generalised_equations

    I'd like to implement your solution in my Ellipsoid class, to exclude ellipsoids from other calculations:

    https://github.com/mdoube/BoneJ/blob/master/src/org/doube/geometry/Ellipsoid.java

    Thanks,

    Michael

    1. If I'm not mistaken, you have \(\mathbf{A} = \mathbf{V} \, \mathbf{U} \, \mathbf{V}^\mathrm{T}\) (with \(\mathbf{U} = \mathrm{diag}(r_x^{-2}, r_y^{-2}, r_z^{-2})\)) such that your ellipsoid is determined by \((\vec{x} - \vec{v})^\mathrm{T} \, \mathbf{A} \, (\vec{x} - \vec{v}) = 1\). In this case, \(\mathbf{A}\) takes the place of \(\mathbf{Q}\) above, but the expressions are in terms of \(\mathbf{R} = \mathbf{Q}^{-1}\).

      A more convenient formulation would have \(\mathbf{U} = \mathrm{diag}(r_x^{-1}, r_y^{-1}, r_z^{-1})\). Then \(\mathbf{A} = \mathbf{M}^{-\mathrm{T}} \, \mathbf{M}^{-1} = \mathbf{V} \, \mathbf{U} \, (\mathbf{V} \, \mathbf{U})^\mathrm{T}\), so \(\mathbf{M}^{-1} = (\mathbf{V} \, \mathbf{U})^\mathrm{T} = \mathbf{U} \, \mathbf{V}^\mathrm{T}\) and \(\mathbf{M} = \mathbf{V}^{-\mathrm{T}} \, \mathbf{U}^{-1}\). Since \(\mathbf{V}\) is a rotation matrix, \(\mathbf{V}^{-\mathrm{T}} = \mathbf{V}\), so \(\mathbf{M} = \mathbf{V} \, \mathbf{U}^{-1}\). \(\mathrm{U}^{-1}\) is of course \(\mathrm{diag}(r_x, r_y, r_z)\). That means your bounding box will be given by

      $$
      \begin{align*}
      x & = v_x \pm \sqrt{r_x^2\,\mathbf{V}_{1,1}^2 + r_y^2\,\mathbf{V}_{1,2}^2 + r_z^2\,\mathbf{V}_{1,3}^2} \\
      y & = v_y \pm \sqrt{r_x^2\,\mathbf{V}_{2,1}^2 + r_y^2\,\mathbf{V}_{2,2}^2 + r_z^2\,\mathbf{V}_{2,3}^2} \\
      z & = v_z \pm \sqrt{r_x^2\,\mathbf{V}_{3,1}^2 + r_y^2\,\mathbf{V}_{3,2}^2 + r_z^2\,\mathbf{V}_{3,3}^2}.
      \end{align*}
      $$

      1. Hi Tavian,

        Thanks for the reply. I've worked out the affine transform, based on sequential scaling, rotation and transformation. The scaling matrix S is diag(ra, rb, rc), the rotation is V and the transform, M is given by M = VS. Then I use the values from M as you originally proposed. (seems to work and a massive speedup, thank you!)

        The result looks very similar to your comment above except that the off-diagonal terms are also multiplied, because the rotation matrix is postmultiplied by the diagonal scaling matrix. (Each column of the affine transform is multiplied by the corresponding element of the diagonal scaling transform)

        In Java:

        public double[] getXMinAndMax() {
        	final double m11 = ev[0][0] * ra;
        	final double m12 = ev[0][1] * rb;
        	final double m13 = ev[0][2] * rc;
        	final double d = Math.sqrt(m11 * m11 + m12 * m12 + m13 * m13);
        	double[] minMax = { cx - d, cx + d };
        	return minMax;
        }
         
        public double[] getYMinAndMax() {
        	final double m21 = ev[1][0] * ra;
        	final double m22 = ev[1][1] * rb;
        	final double m23 = ev[1][2] * rc;
        	final double d = Math.sqrt(m21 * m21 + m22 * m22 + m23 * m23);
        	double[] minMax = { cy - d, cy + d };
        	return minMax;
        }
         
        public double[] getZMinAndMax() {
        	final double m31 = ev[2][0] * ra;
        	final double m32 = ev[2][1] * rb;
        	final double m33 = ev[2][2] * rc;
        	final double d = Math.sqrt(m31 * m31 + m32 * m32 + m33 * m33);
        	double[] minMax = { cz - d, cz + d };
        	return minMax;
        }

Leave a Reply

Your email address will not be published. Required fields are marked *