A few days ago, while reflecting on a problem, I thought of the intersection of three spheres: given the center coordinates and radii of three spheres, find the coordinates of their intersection points. In theory, this is a clearly defined and concise problem with significant practical applications (such as satellite positioning), so a "standard solution" should have been established long ago. However, after searching around, I found that neither English nor Chinese resources seemed to provide a standard, well-documented 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 solution process online in a highly readable manner. Therefore, 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} Our goal is to solve for \boldsymbol{x} by simultaneously solving these three equations. Solving this system in its general form is relatively difficult, but there is a simpler case—when \boldsymbol{o}_1=(0,0,0), \boldsymbol{o}_2=(a,0,0), and \boldsymbol{o}_3=(b,c,0), the equations become: \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 Eq. [eq:s5] from Eq. [eq:s4], we can solve for x. Then, by subtracting Eq. [eq:s6] from Eq. [eq:s4], we can solve for y, and finally 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 \pm in the final step indicates that if an intersection exists, there are generally two points.
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 these conditions are not met, we can choose a new coordinate system for them.
First, subtract \boldsymbol{o}_1 from all coordinates so that \boldsymbol{o}_1 becomes the origin. Let \boldsymbol{o}_{ij} = \boldsymbol{o}_j - \boldsymbol{o}_i. We define the x-axis using the unit vector: \begin{equation} \boldsymbol{u}=\frac{\boldsymbol{o}_{12}}{\Vert\boldsymbol{o}_{12}\Vert} \end{equation} In this case, \boldsymbol{o}_2 lies on the x-axis, and we have a=\Vert\boldsymbol{o}_{12}\Vert=\boldsymbol{o}_{12}\cdot\boldsymbol{u}. Next, we define the y-axis using the unit vector: \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} This ensures 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, we complete the basis with \boldsymbol{w}=\boldsymbol{u}\times \boldsymbol{v}. Thus, \{\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}\} forms a new orthonormal basis.
With a, b, c determined, we can substitute them into equations [eq:s7], [eq:s8], and [eq:s9] to calculate the intersection coordinates (x,y,z). Note that these coordinates are relative to the new coordinate system (\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}). To restore the intersection point to the original coordinate system, we use: \begin{equation} \boldsymbol{x}=x \boldsymbol{u} + y\boldsymbol{v} + z \boldsymbol{w} + \boldsymbol{o}_1 \end{equation}
Overall Process
The derivation above is referenced from the Wikipedia entry for "True-range multilateration," but the process there was not entirely complete, so I have supplemented it here.
Summarizing all the steps: \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 problem of finding the intersection of three spheres.
Original URL: https://kexue.fm/archives/10684
For more details on reposting, please refer to: "Scientific Space FAQ"