Calculating Intersection of Three Spheres Step by Step

Question

How do I calculate the intersection of three spheres step by step?

Assume that the spheres are
$S_i(c_i, r_i)$ where $i = 1,2,3$, $c_i$ is the center coordinates of $S_i$ and $r_i$ is the radius of $S_i$.

First method comes to mind is finding the intersection of two spheres and caculate the plane-sphere intersection.

But again, I'm stuck at

$a^2 + b^2 + c^2 - d^2 - e^2 - f^2 - 2(ax+by+cz)+2(dx+ey+fz) = R^2 - Q^2$ where $S_1((a,b,c), R)$ and
$S_2((d,e,f), Q)$

Could you help me to go on?

Edit:
For center points $\{(a,b,c), (d,e,f), (e,f,g)\}$ and radii $\{r_0,r_1,r_2\}$ I have the matrix

$2\left[\begin{array}{ccc} d-a & e-b & f-c \\ g-d & h-e & i-f \\ a-g & b-h & c-i \end{array}\right]$

and I solve for
$\left[\begin{array}{ccc} P - Q \\ Q - R \\ R - P\end{array}\right]$

where
$P = ({r_0}^2 - a^2 - b^2 - c^2)$
$Q = ({r_1}^2 - d^2 - e^2 - f^2)$
$R = ({r_2}^2 - g^2 - h^2 - i^2)$

But I get $\left[\begin{array}{ccc} NaN \\ NaN \\ NaN\end{array}\right]$ for centers $\{(0,0,0), (3,0,0), (0,4,0)\}$ and radii $\{5,4,3\}$.

Have I built the matrices wrong?

Edit2:
I've made another trial with $(18,41,27),(34,12,12),(31,43,17)$ and $\{76.635500912, 98.493654618, 73.661387443 \}$ still, $NaN$

Answer 1

So by following the algorithm here, I was able to solve the second example you gave: $q = [64; 92; 61 ]$. The first example I do not think has a solution.

The code below is Matlab and is very poorly written. Hopefully, you will find it useful:

function q = trilateration( p, r )
% Inputs:
%   p       3x3 matrix where the columns are the xyz
%           locations of the sphere centers
%   r       3x1 vector of the sphere radii
%
% Outputs:
%   q       3x1 vector of the (positive) point of
%           intersection of the three spheres. If
%           imaginary, then no such point exists

p0 = p(:,1);
ex = unit( p(:,2) - p(:,1) );
i = ex' * ( p(:,3) - p(:,1) );
ey = unit( p(:,3) - p(:,1) - i * ex );
ez = cross( ex, ey );
T = [ ex, ey, ez ];

t = T' * ( p - repmat( p0, 1, 3 ) );

d = t(1,2);
j = t(2,3);

x = ( r(1)^2 - r(2)^2 + d^2 ) / (2*d);
y = ( r(1)^2 - r(3)^2 + i^2 + j^2 ) / (2*j) - i/j * x;
z = sqrt( r(1)^2 - x^2 - y^2 );

q = T * [ x; y; z ] + p0;

end

Here is what I get with your second example:

>> p = [ 18 41 27;
         34 12 12;
         31 43 17 ]';
r = [ 76.635500912; 98.493654618; 73.661387443 ];
>> q = trilateration( p, r )

q =

          64.0000000004705
          91.9999999997112
          61.0000000002842

Answer 2

Take the three equations in AnonSubmitter85's answer. Let's call these equations (1), (2), and (3). Now form three new equations by subtracting: (1) - (2), (2) - (3), and (3) - (1). The important thing is that the three new equations will not longer have any quadratic terms; they will just be linear in $x,y,z$, and therefore fairly easy to solve.

The three linear equations represent three planes, of course. What three planes are these? Well, if you intersect any two of the spheres, you will get a circle, as you said. Doing this two-at-a-time intersection will give you three circles. The three planes mentioned above are the planes containing these three circles. So, you were somewhat on the right track.

Edit: Despite the fact that it got two upvotes, what I wrote above is not really correct. It can't possibly be correct because intersecting the three planes will give only a single point (usually), whereas the intersection of three spheres will give two points (usually).

The three planes will actually intersect in a line, not a point. That's why you get NaN if you try to solve the system of three linear equations. You should take two of the planes and intersect them to get a line. Then intersect this line with a sphere to get the two desired points. My apologies. Moral: don't believe everything you read on the internet :-)

The tedioius details are worked out in this answer.

Also, good answers here, including some with code.