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:
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.
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.
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.
Intersection of two spheres A sphere intersects the plane at infinity in a conic, which is called the absolute conic of the space. The intersection curve of two sphere always degenerates into the absolute conic and a circle. Therefore, the real intersection of two spheres is a circle.
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.
Mathematically, if two spheres touch, the place of contact is one point. If we consider two intersecting spheres, their intersection is a circle.
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With