Coordinates of the Intersection of Three Spheres (Trilateration)

三个球的交点坐标（三球交会定位） 苏剑林 · 2025-01-28 · Mathematics equation geometry

A few days ago, while thinking about a certain problem, I found myself reminded of the three-sphere intersection problem: given the center coordinates and radii of three spheres, find the coordinates of their intersection point(s). This seems like a clean, well-defined problem with an obvious practical application (satellite positioning, for instance), so surely someone must have already worked out a "standard solution." But after searching around, I found that neither English nor Chinese sources offered a standard derivation.

Of course, this isn't to say the problem is so hard that no one has been able to solve it — in fact it's a classic problem that was solved long ago. I was just surprised that no one seems to have written up the derivation online in a reasonably readable way, so this post aims to fill that gap.

A Special Case

First, suppose the equations of the three spheres are

\begin{align} &\text{ball 1:}\quad (\boldsymbol{x} - \boldsymbol{o}_1)^2 = r_1^2 \label{eq:s1} \\ &\text{ball 2:}\quad (\boldsymbol{x} - \boldsymbol{o}_2)^2 = r_2^2 \label{eq:s2} \\ &\text{ball 3:}\quad (\boldsymbol{x} - \boldsymbol{o}_3)^2 = r_3^2 \label{eq:s3} \\ \end{align}more

What we want to do is solve these three equations simultaneously for $\boldsymbol{x}$. Solving this system in general is fairly involved, but there's a simpler special case — when $\boldsymbol{o}_1=(0,0,0),\boldsymbol{o}_2=(a,0,0),\boldsymbol{o}_3=(b,c,0)$, the equations become

\begin{align} &\text{ball 1:}\quad x^2+y^2+z^2 = r_1^2 \label{eq:s4} \\ &\text{ball 2:}\quad (x-a)^2+y^2+z^2 = r_2^2 \label{eq:s5} \\ &\text{ball 3:}\quad (x-b)^2+(y-c)^2+z^2 = r_3^2 \label{eq:s6} \\ \end{align}

Subtracting equation $\eqref{eq:s1}$ from $\eqref{eq:s2}$ lets us solve for $x$, and then subtracting $\eqref{eq:s1}$ from $\eqref{eq:s3}$ lets us solve for $y$. Finally, we can solve for $z$:

\begin{align} &x = \frac{r_1^2 - r_2^2 + a^2}{2a} \label{eq:s7} \\[5pt] &y = \frac{r_1^2 - r_3^2 + b^2 - 2bx + c^2}{2c} \label{eq:s8} \\[5pt] &z = \pm \sqrt{r_1^2 - x^2 - y^2} \label{eq:s9} \end{align}

The final $\pm$ shows that, if an intersection exists, there will generally be two of them.

The General Case

The significance of the special case above is that $\boldsymbol{o}_1$ sits at the origin, $\boldsymbol{o}_2$ lies on the $x$-axis, and $\boldsymbol{o}_3$ lies in the $x,y$-plane. When the three spheres don't already satisfy this arrangement, we can simply choose a new coordinate system for them in which they do.

First, subtract $\boldsymbol{o}_1$ from all the coordinates, which places the origin at $\boldsymbol{o}_1$. Then let $\boldsymbol{o}_{ij} = \boldsymbol{o}_j - \boldsymbol{o}_i$, and take

\begin{equation}\boldsymbol{u}=\frac{\boldsymbol{o}_{12}}{\Vert\boldsymbol{o}_{12}\Vert}\end{equation}

as the $x$-axis; this puts $\boldsymbol{o}_2$ on the $x$-axis, with $a=\Vert\boldsymbol{o}_{12}\Vert=\boldsymbol{o}_{12}\cdot\boldsymbol{u}$. Next, take

\begin{equation}\boldsymbol{v}=\frac{\boldsymbol{o}_{13} - (\boldsymbol{o}_{13}\cdot \boldsymbol{u})\boldsymbol{u}}{\Vert\boldsymbol{o}_{13} - (\boldsymbol{o}_{13}\cdot \boldsymbol{u})\boldsymbol{u}\Vert}\end{equation}

as the $y$-axis; this puts $\boldsymbol{o}_3$ in the $x,y$-plane, with $b = \boldsymbol{o}_{13}\cdot \boldsymbol{u},c = \boldsymbol{o}_{13}\cdot \boldsymbol{v}$. Finally, complete the basis with $\boldsymbol{w}=\boldsymbol{u}\times \boldsymbol{v}$, so that $\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$ together form a new orthonormal basis.

Now that we have $a,b,c$, we can substitute into equation $\eqref{eq:s7},\eqref{eq:s8},\eqref{eq:s9}$ to compute the intersection coordinates $(x,y,z)$. But note that this intersection point is expressed in the new coordinate system $(\boldsymbol{u},\boldsymbol{v},\boldsymbol{w})$; converting back to the original coordinates, the intersection point is

\begin{equation}\boldsymbol{x}=x \boldsymbol{u} + y\boldsymbol{v} + z \boldsymbol{w} + \boldsymbol{o}_1\end{equation}

Putting It All Together

The derivation above is adapted from the Wikipedia entry "True-range multilateration", but that treatment wasn't entirely complete, so this post fills in the missing pieces.

Collecting all the steps together, we get

\begin{equation}\left\{\begin{aligned} &\,\boldsymbol{o}_{ij} = \boldsymbol{o}_j - \boldsymbol{o}_i \\[5pt] &\,\boldsymbol{u}=\frac{\boldsymbol{o}_{12}}{\Vert\boldsymbol{o}_{12}\Vert},\,\,\boldsymbol{v}=\frac{\boldsymbol{o}_{13} - (\boldsymbol{o}_{13}\cdot \boldsymbol{u})\boldsymbol{u}}{\Vert\boldsymbol{o}_{13} - (\boldsymbol{o}_{13}\cdot \boldsymbol{u})\boldsymbol{u}\Vert},\,\,\boldsymbol{w}=\boldsymbol{u}\times \boldsymbol{v}\\[5pt] &\,a=\boldsymbol{o}_{12}\cdot\boldsymbol{u},\,\,b = \boldsymbol{o}_{13}\cdot \boldsymbol{u},\,\,c = \boldsymbol{o}_{13}\cdot \boldsymbol{v} \\[5pt] &\,x = \frac{r_1^2 - r_2^2 + a^2}{2a} \\[5pt] &\,y = \frac{r_1^2 - r_3^2 + b^2 - 2bx + c^2}{2c} \\[5pt] &\,z = \pm \sqrt{r_1^2 - x^2 - y^2} \\[5pt] &\,\boldsymbol{x} = x \boldsymbol{u} + y\boldsymbol{v} + z \boldsymbol{w} + \boldsymbol{o}_1 \end{aligned}\right.\end{equation}

Reference implementation:

import numpy as np

def trilaterate(o1, o2, o3, r1, r2, r3):
    o12, o13 = o2 - o1, o3 - o1
    u = o12 / (o12**2).sum()**0.5
    v = (v := o13 -  o13.dot(u) * u) / (v**2).sum()**0.5
    w = np.cross(u, v)
    a, b, c = o12.dot(u), o13.dot(u), o13.dot(v)
    x = (r1**2 - r2**2 + a**2) / (2 * a)
    y = (r1**2 - r3**2 + b**2 - 2 * b * x + c**2) / (2 * c)
    z = (r1**2 - x**2 - y**2)**0.5
    p = x * u + y * v + o1
    return p + z * w, p - z * w

o1 = np.array([1, 2, -3])
o2 = np.array([2, 1, -1])
o3 = np.array([-3, 0, 2])
r1, r2, r3 = 4, 5, 6
trilaterate(o1, o2, o3, r1, r2, r3)

Summary

This post lays out a relatively concise solution to the three-sphere intersection problem.