Python: Intersection of spheres

Question

I am extremely new to programming but I decided to take on an interesting project as I recently learnt how to represent a sphere in parametric form. When intersecting three spheres, there are two points of intersections that are distinct unless they only overlap at a singular point.

Parametric representation of a sphere:

The code I have is modified from the answer from Python/matplotlib : plotting a 3d cube, a sphere and a vector?, adding the ability to dictate the x, y and z origin and the radius of the sphere. Many similar questions were written in C++, Java, and C#, which I cannot understand at all (I barely know what I am doing so go easy on me).

My Code:

import numpy as np

def make_sphere_x(x, radius):
  u, v = np.mgrid[0:2 * np.pi:5000j, 0:np.pi:2500j]
  x += radius * np.cos(u) * np.sin(v)
  return x

def make_sphere_y(y, radius):
  u, v = np.mgrid[0:2 * np.pi:5000j, 0:np.pi:2500j]
  y += radius * np.sin(u) * np.sin(v)
  return y

def make_sphere_z(z, radius):
  u, v = np.mgrid[0:2 * np.pi:5000j, 0:np.pi:2500j]
  z += radius * np.cos(v)
  return z

#x values
sphere_1_x = make_sphere_x(0, 2)
sphere_2_x = make_sphere_x(1, 3)
sphere_3_x = make_sphere_x(-1, 4)
#y values
sphere_1_y = make_sphere_y(0, 2)
sphere_2_y = make_sphere_y(1, 3)
sphere_3_y = make_sphere_y(0, 4)
#z values
sphere_1_z = make_sphere_z(0, 2)
sphere_2_z = make_sphere_z(1, 3)
sphere_3_z = make_sphere_z(-2, 4)

#intercept of x-values
intercept_x = list(filter(lambda x: x in sphere_1_x, sphere_2_x))
intercept_x = list(filter(lambda x: x in intercept_x, sphere_3_x))
print(intercept_x)

Problems:

1. Clearly there must be a better way of finding the intercepts. Right now, the code generates points at equal intervals, with the number of intervals I specify under the imaginary number in np.mgrid. If this is increased, the chances of an intersection should increase (I think) but when I try to increase it to 10000j or above, it just spits a memory error.

2. There are obvious gaps in the array and this method would most likely be erroneous even if I have access to a super computer and can crank up the value to an obscene value. Right now the code results in a null set.

3. The code is extremely inefficient, not that this is a priority but people like things in threes right?

Feel free to flame me for rookie mistakes in coding or asking questions on Stack Overflow. Your help is greatly valued.

Answer 1

Using scipy.optimize.fsolve you can find the root of a given function, given an initial guess that is somewhere in the range of your solution. I used this approach to solve your problem and it seems to work for me. The only downside is that it only provides you one intersection. To find the second one you would have to tinker with the initial conditions until fsolve finds the second root.

First we define our spheres by defining (arbitrary) radii and centers for each sphere:

a1 = np.array([0,0,0])
r1 = .4
a2 = np.array([.3,0,0])
r2 = .5
a3 = np.array([0,.3,0])
r3 = .5

We then define how to transform back into cartesian coordinates, given angles u,v

def position(a,r,u,v):
    return a + r*np.array([np.cos(u)*np.sin(v),np.sin(u)*np.sin(v),np.cos(v)])

Now we think about what equation we need to find the root of. For any intersection point, it holds that for perfect u1,v1,u2,v2,u3,v3 the positions position(a1,r1,u1,v1) = position(a2,r2,u2,v2) = position(a3,r3,u3,v3) are equal. We thus find three equations which must be zeros, namely the differences of two position vectors. In fact, as every vector has 3 components, we have 9 equations which is more than enough to determine our 6 variables.

We find the function to minimize as:

def f(args):
    u1,v1,u2,v2,u3,v3,_,_,_ = args
    pos1 = position(a1,r1,u1,v1)
    pos2 = position(a2,r2,u2,v2)
    pos3 = position(a3,r3,u3,v3)
    return np.array([pos1 - pos2, pos1 - pos3, pos2 - pos3]).flatten()

fsolve needs the same amount of input and output arguments. As we have 9 equations but only 6 variables I simply used 3 dummy variables so the dimensions match. Flattening the array in the last line is necessary as fsolve only accepts 1D-Arrays.

Now the intersection can be found using fsolve and a (pretty random) guess:

guess = np.array([np.pi/4,np.pi/4,np.pi/4,np.pi/4,np.pi/4,np.pi/4,0,0,0])

x0 = fsolve(f,guess)
u1,v1,u2,v2,u3,v3,_,_,_ = x0

You can check that the result is correct by plugging the angles you received into the position function.

Answer 2

The problem would be better tackled using trigonometry.

Reducing the problem into 2D circles, we could do:

import math
import numpy

class Circle():
    def __init__(self, cx, cy, r):
        """initialise Circle and set main properties"""
        self.centre = numpy.array([cx, cy])
        self.radius = r

    def find_intercept(self, c2):
        """find the intercepts between the current Circle and a second c2"""
        #Find the distance between the circles
        s = c2.centre - self.centre
        self.dx, self.dy = s
        self.d = math.sqrt(numpy.sum(s**2))
        #Test if there is an overlap.  Note: this won't detect if one circle completly surrounds the other.
        if self.d > (self.radius + c2.radius):
            print("no interaction")
        else:
            #trigonometry
            self.theta = math.atan2(self.dy,self.dx)
            #cosine rule
            self.cosA = (c2.radius**2 - self.radius**2 + self.d**2)/(2*c2.radius*self.d)
            self.A = math.acos(self.cosA)

            self.Ia = c2.centre - [math.cos(self.A+self.theta)*c2.radius, math.sin(self.A+self.theta)*c2.radius]
            self.Ib = c2.centre - [math.cos(self.A-self.theta)*c2.radius,-math.sin(self.A-self.theta)*c2.radius]

            print("Interaction points are : ", self.Ia, " and: ", self.Ib)


#define two arbitrary circles 
c1 = Circle(2,5,5)
c2 = Circle(1,6,4)


#find the intercepts
c1.find_intercept(c2)

#test results by reversing the operation
c2.find_intercept(c1)