I work in Python 3 on astrophysics projects. I need to compute the smallest circumscribed circle of a set of points in the sky (so described by Right Ascension and Declination). I have found a code here (real-time demo here).

It implements an algorithm in $\mathcal O(n)$ instead of the naive $\mathcal O(n^4)$ one, which is very interesting for me as I have large datasets. But it is written in the Euclidean plane. I have transformed it as much as possible to match spherical trigonometry on the sky, but there is one part I don't know how to convert.

My current code is then a mix of spherical geometry and plane approximation. Unfortunately, although my objects are very close in the sky (they are galaxies in compact groups), this approximation is not good enough and I need the fully-spherical computation.

Hereafter my modified code, with as many comments as possible (Original code, Spherical code, and so on) to describe what I have done and the part where I don't know what to do. It uses import numpy as np and import math.

# From https://www.nayuki.io/res/smallest-enclosing-circle/smallestenclosingcircle.py
# Demo on https://www.nayuki.io/page/smallest-enclosing-circle
def calc_sep1(a,b):

    return math.acos(min(1,math.sin(theta1)*math.sin(theta2) + 
                     math.cos(theta1)*math.cos(theta2)*math.cos(phi2 - phi1)) )

#Defines barycenter of a triangle
def calc_bary1(a,b,c):
    phi = [a[0],b[0],c[0]]
    theta = [a[1],b[1],c[1]]
    x = np.cos(theta) * np.cos(phi)
    y = np.cos(theta) * np.sin(phi)
    z = np.sin(theta)

    x_O, y_O, z_O = np.mean(x), np.mean(y), np.mean(z) 

    l = np.sqrt(x_O**2 + y_O**2 + z_O**2)

    (x_M, y_M, z_M) = (x_O, y_O, z_O)/l

    Dec_M = np.arcsin(z_M)
    if (y_M <0):
        RA_M = 2.0*np.pi - np.arccos(x_M / np.cos(Dec_M))
        RA_M = np.arccos(x_M / np.cos(Dec_M))

    return (RA_M,Dec_M)

def midpoint(a,b):

    Bx = math.cos(theta2)*math.cos(phi2-phi1)
    By = math.cos(theta2)*math.sin(phi2-phi1)
    theta_m = math.atan2(math.sin(theta1) + math.sin(theta2), math.sqrt((math.cos(theta1)+Bx)**2 + By**2) )
    phi_m = phi1 + math.atan2(By, math.cos(theta1)+Bx)

    return (phi_m,theta_m)

# Data conventions: A point is a pair of floats (x, y). 
#                   A circle is a triple of floats (center x, center y, radius).

# Returns the smallest circle that encloses all the given points. Runs in expected O(n) time, randomized.
# Input: A sequence of pairs of floats or ints, e.g. [(0,5), (3.1,-2.7)].
# Output: A triple of floats representing a circle.
# Note: If 0 points are given, None is returned. If 1 point is given, a circle of radius 0 is returned.
# Initially: No boundary points known
def make_circle(points):
    # Randomize order
    # Carfull: inpuot and outpu in degrees, inside: radiants
    shuffled = [(d2r(row['RA']), d2r(row['Dec'])) for index, row in points.iterrows()]
    # Progressively add points to circle or recompute circle
    c = None
    for (i, p) in enumerate(shuffled):
        if c is None or not is_in_circle(c, p):
            c = _make_circle_one_point(shuffled[ : i + 1], p)
    return (r2d(c[0]),r2d(c[1]),c[2])

# One boundary point known
def _make_circle_one_point(points, p):
    c = (p[0], p[1], 0.0)
    for (i, q) in enumerate(points):
        if not is_in_circle(c, q):
            if c[2] == 0.0:
                c = make_diameter(p, q)
                c = _make_circle_two_points(points[ : i + 1], p, q)
    return c

# Two boundary points known
def _make_circle_two_points(points, p, q):
    circ = make_diameter(p, q)
    left  = None
    right = None
    px, py = p
    qx, qy = q

    # For each point not in the two-point circle
    for r in points:
        if is_in_circle(circ, r):

        # Form a circumcircle and classify it on left or right side
        cross = _cross_product(px, py, qx, qy, r[0], r[1])
        c = make_circumcircle(p, q, r)
        if c is None:
        elif cross > 0.0 and (left is None or _cross_product(px, py, qx, qy, c[0], c[1]) > 
                              _cross_product(px, py, qx, qy, left[0], left[1])):
            left = c
        elif cross < 0.0 and (right is None or _cross_product(px, py, qx, qy, c[0], c[1]) < 
                              _cross_product(px, py, qx, qy, right[0], right[1])):
            right = c

    # Select which circle to return
    if left is None and right is None:
        return circ
    elif left is None:
        return right
    elif right is None:
        return left
        return left if (left[2] <= right[2]) else right

def make_diameter(a, b):
    # Original code:
#     cx = (a[0] + b[0]) / 2.0
#     cy = (a[1] + b[1]) / 2.0
#     r0 = math.hypot(cx - a[0], cy - a[1])
#     r1 = math.hypot(cx - b[0], cy - b[1])
#     return (cx, cy, max(r0, r1))

    # Sperical code:
    (phi_m,theta_m) = midpoint(a,b)
    return (phi_m,theta_m,calc_sep1(a,b)*.5)

def make_circumcircle(a, b, c):
    # Mathematical algorithm from Wikipedia: Circumscribed circle

    # I don't know how to sphericalize this part
    ox = (min(a[0], b[0], c[0]) + max(a[0], b[0], c[0])) / 2.0
    oy = (min(a[1], b[1], c[1]) + max(a[1], b[1], c[1])) / 2.0

    ax = a[0] - ox;  ay = a[1] - oy
    bx = b[0] - ox;  by = b[1] - oy
    cx = c[0] - ox;  cy = c[1] - oy

    d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2.0
    if d == 0.0:
        return None
    x = ox + ((ax * ax + ay * ay) * (by - cy) + 
              (bx * bx + by * by) * (cy - ay) + 
              (cx * cx + cy * cy) * (ay - by)) / d
    y = oy + ((ax * ax + ay * ay) * (cx - bx) + 
              (bx * bx + by * by) * (ax - cx) + 
              (cx * cx + cy * cy) * (bx - ax)) / d

    # Original code:
#     ra = math.hypot(x - a[0], y - a[1])
#     rb = math.hypot(x - b[0], y - b[1])
#     rc = math.hypot(x - c[0], y - c[1])

    # Spherical code:
    ra = calc_sep1((x,y),a)
    rb = calc_sep1((x,y),b)
    rc = calc_sep1((x,y),c)

    return (x, y, max(ra, rb, rc))


def is_in_circle(c, p):
#     return c is not None and math.hypot(p[0] - c[0], p[1] - c[1]) <= c[2] * _MULTIPLICATIVE_EPSILON
    return c is not None and calc_sep1(p,(c[0],c[1])) <= c[2] * _MULTIPLICATIVE_EPSILON

# Returns twice the signed area of the triangle defined by (x0, y0), (x1, y1), (x2, y2).
def _cross_product(x0, y0, x1, y1, x2, y2):
    # Original code:
    #return (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)

    # Spherical code:
    # (Uses http://mathforum.org/library/drmath/view/65316.html 
    # then adds sign accordingly to ordinary cross product)

    side1 = calc_sep1((x0, y0), (x1, y1))
    side2 = calc_sep1((x1, y1), (x2, y2))
    side3 = calc_sep1((x2, y2), (x0, y0))
    s = (side1+side2+side3)/2
    Area = 4* math.atan(math.sqrt(math.tan(s/2)*

    Eucl_cross_prod = (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)
    cross_prod_sign = 1 if (Eucl_cross_prod >=0) else -1

    return 2*cross_prod_sign*Area
  • 1
    $\begingroup$ Is your problem to find the smallest-circle that encloses a set of points that are constrained to a sphere? $\endgroup$
    – nicoguaro
    Commented Nov 20, 2018 at 14:56
  • 1
    $\begingroup$ My problem is to find the smallest circle that encloses a set of galaxies, for which I only have angles in spherical coordinates. It would be the same for cities with latitudes and longitudes on earth. So I suppose you can say they are constrained to a sphere, yes. $\endgroup$
    – Matt
    Commented Nov 21, 2018 at 10:30
  • $\begingroup$ Do you have any guarantees that all the points will be on the same hemisphere or not (not necessarily aligned north-south or east-west, just that the enclosing circle must be the same size as or smaller than a great circle on the sphere)? $\endgroup$ Commented Nov 24, 2018 at 5:24
  • $\begingroup$ The reason I ask is because if you can guarantee that all the points are on the same hemisphere, then you can modify the 2D algorithm to work on the surface of a sphere just fine. However, if you can't, then I think the easiest and possibly fastest solution would be to convert all points to 3D Cartesian coordinates on a unit sphere and use the 3D version of the smallest enclosing sphere. $\endgroup$ Commented Nov 24, 2018 at 5:48
  • 1
    $\begingroup$ @Markus Klyver I have a working algorithm in an old language (Super Mongo) which I have translated into Python but is bugged for the moment (I had found a way not to use my python routine until now). I have to debug it. $\endgroup$
    – Matt
    Commented Jan 12, 2021 at 15:38

1 Answer 1


I found this link which I translated into the following Python code (for a terrestrial application, for your application substitute calc_sep1 for geopy.distance.great_circle):

import numpy as np, math
from geopy.distance import great_circle

def makeCircumcircle(a, b, c):
    '''[Link](https://web.archive.org/web/20171023010630/http://mathforum.org/library/drmath/view/68373.html), intersection of
    sphere with plane containing all three points'''
    lats = np.radians([a[0], b[0], c[0]])
    lons = np.radians([a[1], b[1], c[1]])
    cosLats = np.cos(lats)
    xyz = np.transpose([cosLats * np.cos(lons), cosLats * np.sin(lons),
    N = np.cross(xyz[1] - xyz[0], xyz[2] - xyz[0])
    N /= np.linalg.norm(N)
    if np.dot(xyz[0], N) < 0.: N = - N
    centre = [math.degrees(math.asin(N[2])),
        math.degrees(math.atan2(N[1], N[0]))]
    return centre + [great_circle(a[:2], centre)]

>>> makeCircumcircle([0,0],[1,.5],[0,1])
[0.3750017848988929, 0.49999999999990025, Distance(69.4967288556)]
>>> makeCircumcircle([0,0],[1,1],[0,2])
[-1.8212037231487878e-13, 0.9999999999999474, Distance(111.195083724)]
>>> makeCircumcircle([60,-178],[60,178],[61,180])
[60.00003866151654, -180.0, Distance(111.190784754)]
>>> makeCircumcircle([89,-90],[89,90],[89,0])
[90.0, -180.0, Distance(111.195083724)]
>>> makeCircumcircle([89,-90],[89,90],[90,0])
[2.2633064952667684e-11, -180.0, Distance(10007.5575352)]

all seem to make sense to me.

  • 1
    $\begingroup$ The interesting part starts at more than 3 points. $\endgroup$
    – Matt
    Commented Jan 12, 2021 at 15:40

Your Answer

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

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