Intersection Coordinates of Three Spheres (Trilateration)

By 苏剑林 | January 28, 2025

A few days ago, while pondering a problem, I thought of the three-sphere intersection problem—that is, given the coordinates of the centers and the radii of three spheres, find the coordinates of their intersection points. Logically, this is a well-defined and concise problem with clear application backgrounds (such as satellite positioning), so a "standard answer" should have been provided long ago. However, after searching around, I found that neither English nor Chinese resources seem to provide a readable, standard derivation process.

Of course, this is not to say that the problem is so difficult that no one can solve it; in fact, it is a classic problem that was solved long ago. I was simply surprised that no one seemed to have written down the derivation process on the internet in a highly readable way, so this article attempts to fill that gap.

Special Case

First, let the equations of the three spheres be: \begin{align} &\text{Sphere 1:}\quad (\boldsymbol{x} - \boldsymbol{o}_1)^2 = r_1^2 \label{eq:s1} \\ &\text{Sphere 2:}\quad (\boldsymbol{x} - \boldsymbol{o}_2)^2 = r_2^2 \label{eq:s2} \\ &\text{Sphere 3:}\quad (\boldsymbol{x} - \boldsymbol{o}_3)^2 = r_3^2 \label{eq:s3} \end{align} What we need to do is solve for $\boldsymbol{x}$ by simultaneously solving these three equations. Solving the general equations is relatively difficult, but there is a simpler example—when $\boldsymbol{o}_1=(0,0,0), \boldsymbol{o}_2=(a,0,0), \boldsymbol{o}_3=(b,c,0)$, the equations are: \begin{align} &\text{Sphere 1:}\quad x^2+y^2+z^2 = r_1^2 \label{eq:s4} \\ &\text{Sphere 2:}\quad (x-a)^2+y^2+z^2 = r_2^2 \label{eq:s5} \\ &\text{Sphere 3:}\quad (x-b)^2+(y-c)^2+z^2 = r_3^2 \label{eq:s6} \end{align} By subtracting equation $\eqref{eq:s4}$ from $\eqref{eq:s5}$, we can solve for $x$. Then, by subtracting equation $\eqref{eq:s4}$ from $\eqref{eq:s6}$, we can solve for $y$, and finally, we can determine $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$ indicates that if the intersection points exist, there will generally be two of them.

General Case

The significance of the above example is that $\boldsymbol{o}_1$ is at the origin, $\boldsymbol{o}_2$ is on the $x$-axis, and $\boldsymbol{o}_3$ is in the $xy$-plane. When they do not satisfy these conditions, we can choose a new coordinate system for them.

First, subtract $\boldsymbol{o}_1$ from all coordinates, achieving a system where $\boldsymbol{o}_1$ is the origin. Then we denote $\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. Thus $\boldsymbol{o}_2$ lies on the $x$-axis, and we have $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, so that $\boldsymbol{o}_3$ lies in the $xy$-plane, with $b = \boldsymbol{o}_{13}\cdot \boldsymbol{u}$ and $c = \boldsymbol{o}_{13}\cdot \boldsymbol{v}$. Finally, append $\boldsymbol{w}=\boldsymbol{u}\times \boldsymbol{v}$, so that $\boldsymbol{u}, \boldsymbol{v}, \boldsymbol{w}$ form a new orthonormal basis.

Now that we have $a, b, c$, we can substitute them into equations $\eqref{eq:s7}, \eqref{eq:s8}, \eqref{eq:s9}$ to calculate the intersection coordinates $(x,y,z)$. However, note that these coordinates are obtained in the new coordinate system $(\boldsymbol{u}, \boldsymbol{v}, \boldsymbol{w})$. The intersection points in the original coordinates are: \begin{equation}\boldsymbol{x} = x \boldsymbol{u} + y\boldsymbol{v} + z \boldsymbol{w} + \boldsymbol{o}_1\end{equation}

Total Process

The above solution process is referenced from the Wikipedia entry "True-range multilateration," but the process there was not quite complete, so this article has supplemented it to be full.

Combining all the steps together: \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 article has organized a relatively concise solution for the three-sphere intersection problem.