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, as 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
  • $\begingroup$ Is your problem to find the smallest-circle that encloses a set of points that are constrained to a sphere? $\endgroup$ – nicoguaro Nov 20 '18 at 14:56
  • $\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 Nov 21 '18 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$ – helloworld922 Nov 24 '18 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$ – helloworld922 Nov 24 '18 at 5:48
  • $\begingroup$ If your star are close enough to be able to assume that the arc distance is almost the same of the Euclidean ($\sin x= x$ approximation), you can convert ascension and declination to Cartesian coordinates and use the original code. You just need to scale the declination according to the cosine of the ascension (if I correctly remember their definition). $\endgroup$ – N74 Nov 24 '18 at 17:26

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.