I'm looking for an algorithm to find the common intersection points between 3 spheres.

Baring a complete algorithm, a thorough/detailed description of the math would be greatly helpful.

This is the only helpful resource I have found so far: http://mathforum.org/library/drmath/view/63138.html

But neither method described there is detailed enough for me to write an algorithm on.

I would prefer the purely algebraic method described in the second post, but what ever works.

UPDATE

An implementation of this answer in python complete with an example of usage can be found at this github repo.

It turns out the analytic solution is actually quite nice using this method and can tell you when a solution exists and when it doesn't (it is also possible to have exactly one solution.) There is no reason to use Newton's method.

IMHO, this is far easier to understand and simpler than trilateration given below. However, both techniques give correct answers in my testing.

ORIGINAL ANSWER

Consider the intersection of two spheres. To visualize it, consider the 3D line segment N connecting the two centers of the spheres. Consider this cross section

_{(source: googlepages.com)}

where the red-line is the cross section of the plane with normal N. By symmetry, you can rotate this cross-section from any angle, and the red line segments length can not change. This means that the resulting curve of the intersection of two spheres is a circle, and must lie in a plane with normal N.

That being said, lets get onto finding the intersection. First, we want to describe the resulting circle of the intersection of two spheres. You can not do this with 1 equation, a circle in 3D is essentially a curve in 3D and you cannot describe curves in 3D by 1 eq.

Consider the picture
_{(source: googlepages.com)}

let P be the point of intersection of the blue and red line. Let h be the length of the line segment along the red line from point P upwards. Let the distance between the two centers be denoted by d. Let x be the distance from the small circle center to P. Then we must have

x^2 +h^2 = r1^2
(d-x)^2 +h^2 = r2^2
==> h = sqrt(r1^2 - 1/d^2*(r1^2-r2^2+d^2)^2)

i.e. you can solve for h, which is the radius of the circle of intersection. You can find the center point C of the circle from x, along the line N that joins the 2 circle centers.

Then you can fully describe the circle as (X,C,U,V are all vector)

X = C + (h * cos t) U + (h * sin t) V for t in [0,2*PI)

where U and V are perpendicular vectors that lie in a plane with normal N.

The last part is the easiest. It remains only to find the intersection of this circle with the final sphere. This is simply a plug and chug of the equations (plug in for x,y,z in the last equation the parametric forms of x,y,z for the circle in terms of t and solve for t.)

edit ---

The equation that you will get is actually quite ugly, you will have a whole bunch of sine's and cosine's equal to something. To solve this you can do it 2 ways:

1. write the cosine's and sine's in terms of exponentials using the equality

   e^(it) = cos t + i sin t

   then group all the e^(it) terms and you should get a quadratic equations of e^(it)'s that you can solve for using the quadratic formula, then solve for t. This will give you the exact solution. This method will actually tell you exactly if a solution exists, two exist or one exist depending on how many of the points from the quadratic method are real.

2. use newton's method to solve for t, this method is not exact but its computationally much easier to understand, and it will work very well for this case.