15

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.

3
  • 1
    Can you confirm that you mean spherical surfaces rather than solids, and add something so that this is anything other than a math question and is somehow programming related? Sep 10, 2009 at 16:43
  • Well I want a C++ algorithm to do this, but I first need to understand the math behind it. As for the other part of your question, yes, just the surface of the sphere.
    – Adam
    Sep 10, 2009 at 16:49
  • just google for 'trilateration': en.wikipedia.org/wiki/Trilateration
    – thalm
    Feb 1, 2013 at 20:23

6 Answers 6

15

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                       
1
  • 6
    Andrew I love you. Note that P1, P2, P3 are created with numpy.array([x,y,z]). Mar 12, 2016 at 18:17
9

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

alt text

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

2
  • 3 planes need not have a a common line. 2 planes that intersect result in either a plane or a line. 3 planes that intersect result in either a plane, a line or a point.
    – ldog
    Sep 15, 2009 at 22:10
  • 1
    In this case, since the equations of the planes are not independent, the intersection is either a plane, or a line (possibly at infinity when the three centers are collinear). In the general case, the intersection will be a line. Sep 16, 2009 at 5:36
7

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

alt text
(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 alt text
(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.

9
  • The sphere-plane intersection, when it exists, is a circle, always. Sep 10, 2009 at 20:10
  • edited: The sphere-plane intersection, when it exists, is a circle, always[citation needed] Sep 10, 2009 at 20:15
  • Hmm are you sure? I remember in univ, I assumed this my prof marked it wrong... maybe I'm rememebring wrong let me consider it more thoroughly
    – ldog
    Sep 10, 2009 at 20:26
  • To see this easily (this is not a demonstration). Consider a sphere and a plane intersecting it. Get the perpendicular to the plane by the center of the sphere, and rotate everything around this line. Nothing changes, so the intersection must be invariant by rotation in the plane. Sep 10, 2009 at 20:29
  • 2
    Only true if you define a point as a circle, though, which may upset professors ;) Sep 11, 2009 at 8:34
5

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.

0
3

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.

1
  • This is a pretty neat way to look at it. I don't think it's less complicated than my proposed method, it simply moves the "complication" into algebraic manipulation. One thing to point out is it heavily uses the rotational symmetry of a sphere about every axis, which makes this method impossible to extend to any curves or surfaces without this property. Nice find though.
    – ldog
    May 21, 2013 at 18:34
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.