1
$\begingroup$

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.

# Circumscribed circle - PARTIALLY ADAPTED TO SPHERICAL GEOMETRY! 
# 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):
    (phi1,theta1)=(a[0],a[1])
    (phi2,theta2)=(b[0],b[1])

    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))
    else:
        RA_M = np.arccos(x_M / np.cos(Dec_M))

    return (RA_M,Dec_M)

def midpoint(a,b):
    (phi1,theta1)=(a[0],a[1])
    (phi2,theta2)=(b[0],b[1])

    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()]
    random.shuffle(shuffled)   
    # 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)
            else:
                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):
            continue

        # 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:
            continue
        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
    else:
        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))


_MULTIPLICATIVE_EPSILON = 1 + 1e-14

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)*
                                  math.tan((s-side1)/2)*
                                  math.tan((s-side2)/2)*
                                  math.tan((s-side3)/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
$\endgroup$
  • $\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 acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.