Finding intersection points between 3 spheres

Question

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.

Answer 1

Andrew Wagner's solution is very good, and I needed to use it; but it can be improved: add some comments, clarify some variable names, and reuse variables better. Also, add a non-fatal option for the case where you just want a simple approximation in the case where there is no exact intersection.

def trilaterate(
    p: np.ndarray,      # 3*3
    norms: np.ndarray,  # *3
    fatal_radical: bool = False,
) -> np.array:  # 2*3
    """
    https://en.wikipedia.org/wiki/True-range_multilateration#Three_Cartesian_dimensions,_three_measured_slant_ranges

    :param p 3*3 array where the first axis is the observation and the second axis is xyz
    :param norms 3-array, the distance from each observation to an unknown centroid
    :param fatal_radical True to throw if no solution; false to assume a close-enough solution
    :return a 2x3 array of centroid solutions, or 1x3 if there is either one solution, or
            fatal_radical is false and there are no exact solutions

    Modified from https://stackoverflow.com/a/18654302/313768
    """
    r1, r2, r3 = norms
    p1, p2, p3 = p

    # Inter-point vectors
    p21 = p2 - p1
    p31 = p3 - p1
    d = np.linalg.norm(p21)

    # Basis vectors
    u = p21/d
    i = u.dot(p31)
    v = p31 - u*i
    v /= np.linalg.norm(v)
    j = v.dot(p31)
    w = np.cross(u, v)

    # Solution dimensions, still in projected space
    x = 0.5/d*(r1*r1 - r2*r2 + d*d)
    y = 0.5/j*(r1*r1 - r3*r3 - 2*i*x + i*i + j*j)
    radicand = r1*r1 - x*x - y*y

    # Solution vectors in original space
    sxy = p1 + u*x + v*y  # First two dimensions, missing 'z'

    if radicand < 0:
        if fatal_radical:
            raise ValueError(f'Negative radicand {radicand}: no exact solutions')
    if radicand <= 0:
        return sxy[np.newaxis, :]

    z = np.sqrt(radicand)
    sz = w*z              # 'z', prior to sign change
    sa = sxy + sz         # First solution
    sb = sxy - sz         # Second solution
    return np.stack((sa, sb))  # Both solutions: 2*3

Answer 2

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.

Answer 3

Probably easier than constructing 3D circles, because working mainly on lines and planes:

For each pair of spheres, get the equation of the plane containing their intersection circle, by subtracting the spheres equations (each of the form X^2+Y^2+Z^2+aX+bY+c*Z+d=0). Then you will have three planes P12 P23 P31.

These planes have a common line L, perpendicular to the plane Q by the three centers of the spheres. The two points you are looking for are on this line. The middle of the points is the intersection H between L and Q.

To implement this:

• compute the equations of P12 P23 P32 (difference of sphere equations)
• compute the equation of Q (solve a linear system, or compute a cross product)
• compute the coordinates of point H intersection of these four planes. (solve a linear system)
• get the normal vector U to Q from its equation (normalize a vector)
• compute the distance t between H and a solution X: t^2=R1^2-HC1^2, (C1,R1) are center and radius of the first sphere.
• solutions are H+tU and H-tU

A Cabri 3D construction showing the various planes and line L

Answer 4

Here is an answer in Python I just ported from the Wikipedia article. There is no need for an algorithm; there is a closed form solution.

import numpy                                             
from numpy import sqrt, dot, cross                       
from numpy.linalg import norm                            

# Find the intersection of three spheres                 
# P1,P2,P3 are the centers, r1,r2,r3 are the radii       
# Implementaton based on Wikipedia Trilateration article.                              
def trilaterate(P1,P2,P3,r1,r2,r3):                      
    temp1 = P2-P1                                        
    e_x = temp1/norm(temp1)                              
    temp2 = P3-P1                                        
    i = dot(e_x,temp2)                                   
    temp3 = temp2 - i*e_x                                
    e_y = temp3/norm(temp3)                              
    e_z = cross(e_x,e_y)                                 
    d = norm(P2-P1)                                      
    j = dot(e_y,temp2)                                   
    x = (r1*r1 - r2*r2 + d*d) / (2*d)                    
    y = (r1*r1 - r3*r3 -2*i*x + i*i + j*j) / (2*j)       
    temp4 = r1*r1 - x*x - y*y                            
    if temp4<0:                                          
        raise Exception("The three spheres do not intersect!");
    z = sqrt(temp4)                                      
    p_12_a = P1 + x*e_x + y*e_y + z*e_z                  
    p_12_b = P1 + x*e_x + y*e_y - z*e_z                  
    return p_12_a,p_12_b                       

Answer 5

after searching the web this is one of the first hits, so i am posting the most clean and easy solution i found after some hours of research here: Trilateration

This wiki site contains a full description of a fast and easy to understand vector approach, so one can code it with little effort.

Answer 6

Here is another interpretation of the picture which Eric posted above:

Let H be the plane spanned by the centers of the three spheres. Let C1,C2,C3 be the intersections of the spheres with H, then C1,C2,C3 are circles. Let Lij be the line connecting the two intersection points of Ci and Cj, then the three lines L12,L23,L13 intersect at one point P. Let M be the line orthogonal to H through P, then your two points of intersection lie on the line M; hence you just need to intersect M with either of the spheres.

Answer 7

Basically you need to do this in 3 steps. Let's say you've got three spheres, S1, S2, and S3.

1. C12 is the circle created by the intersection of S1 and S2.
2. C23 is the circle created by the intersection of S2 and S3.
3. P1, P2, are the intersection points of C12 and C13.

The only really hard part in here is the sphere intersection, and thankfully Mathworld has that solved pretty well. In fact, Mathworld also has the solution to the circle intersections.

From this information you should be able to create an algorithm.