Finding intersection points between 3 spheres

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                       

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

Tags:

Math