I am using vedo in Python, and I have two spheres. I need to find the circular surface, that is spanned by all points of intersection of the spheres. I couldn't find a good way to do it, and I wonder if there's a quick and easy solution to e.g. have this intersection surface plotted as a Disc.

I am using vedo in Python, and I have two spheres. I need to find the circular surface, that is spanned by all points of intersection of the spheres. I couldn't find a good way to do it, and I wonder if there's a quick and easy solution to e.g. have this intersection surface plotted as a Disc.

• python
• geometry
• vedo

2 Answers 2

Let the two spheres have centres a and b and radii A and B respectively. Let C be the distance between their centres, with c=b-a the vector from centre a to centre b.

If C > A+B then the spheres don't intersect and if C < abs(B-A) then one engulfs the other. Otherwise their intersection circle lies a (signed) distance d=(C²+A²-B²)/(2C) from a in the direction of b, and has a radius r=sqrt(A²-d²). (To see this, use Pythagoras's Theorem on both sides of the plane of intersection in the diagram below. It also works, if d is correctly signed, if the intersection is beyond the centre a.)

Take two unit vectors n₁ and n₁ that are perpendicular to the vector c that joins the centres. (A reasonable way of choosing these is to choose n₁ to be the (normalised) vector from the largest of the vectors c x i, c x j, c x k, where i, j, k are the unit cartesian vectors and x gives a cross product. The third vector n₂ is then the unit vector (c/C) x n₁.

Finally, the circle is generated as the set of points a+d(c/C)+n₁ r cos(t) + n₂ r sin(t) where t runs from 0 to 2.pi.

In code this might look like

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

a = np.array( [ 1.0, 2.0, 3.0 ] );   A = 4.0
b = np.array( [ 7.0, 6.0, 5.0 ] );   B = 6.0

c = b - a                                                  # vector between centres
C = np.linalg.norm( c )                                    # distance between centres
if C > A + B or C < abs( B - A ):
    print( "No intersection" )
    exit()

d = ( C ** 2 + A ** 2 - B ** 2 ) / ( 2 * C )               # distance from a (positive if in direction a to b)
r = math.sqrt( A ** 2 - d ** 2 )                           # radius of intersection circle
centre = a + d * c / C                                     # centre of intersection circle

n1a = np.array( [     0,  c[2], -c[1] ] )                  # candidates for in-plane vector n1
n1b = np.array( [ -c[2],     0, -c[0] ] )
n1c = np.array( [  c[1], -c[0],     0 ] )
n1 = n1a
if np.linalg.norm( n1b ) > np.linalg.norm( n1 ): n1 = n1b
if np.linalg.norm( n1c ) > np.linalg.norm( n1 ): n1 = n1c
n1 /= np.linalg.norm( n1 )                                 # normalise

n2 = np.cross( c, n1 ) / C                                 # second in-plane vector

numpts = 101
theta = np.linspace( 0, 2 * math.pi, numpts )
pts = np.zeros( ( 3, numpts ) )
for i in range( numpts ):
   pts[:,i] = centre + n1 * r * math.cos( theta[i] ) + n2 * r * math.sin( theta[i] )


def plotSphere( centre, radius ):
    u, v = np.mgrid[0:2*math.pi:50j, 0:math.pi:50j]
    x = centre[0] + radius * np.cos(u) * np.sin(v)
    y = centre[1] + radius * np.sin(u) * np.sin(v)
    z = centre[2] + radius * np.cos(v)
    ax.plot_surface( x, y, z, alpha=0.2 )


fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
plotSphere( a, A )
plotSphere( b, B )
ax.scatter( a[0], a[1], a[2] )
ax.scatter( b[0], b[1], b[2] )
ax.plot( pts[0,:], pts[1,:], pts[2,:], 'r-', linewidth=2 )
ax.set_aspect('equal')
plt.show()

With vedo you can do:

from vedo import Sphere, show
s1 = Sphere(r=1.2).pos(-1, 0, 0).c('red5').alpha(0.2)
s2 = Sphere(r=1).pos(0.1, 0.2, 0.3).c('blue5').alpha(0.2)
disc = s1.intersect_with(s2).triangulate()
disc.c('white').lw(1).lighting("off")
print(disc.coordinates.shape)
show(s1, s2, disc, axes=1)

This generalizes to any arbitrary polygonal surface.