smallest_circle.py

"""===============================================================================

        FILE: /Users/nailbiter/Documents/datawise/alex_python_toolbox/alex_python_toolbox/smallest_circle.py

       USAGE: (not intended to be directly executed)

 DESCRIPTION: implementation of Welzl's algorithm (https://en.wikipedia.org/wiki/Smallest-circle_problem#Welzl's_algorithm) on sphere

     OPTIONS: ---
REQUIREMENTS: ---
        BUGS: ---
       NOTES: ---
      AUTHOR: Alex Leontiev (nailbiter@dtws-work.in)
ORGANIZATION: Datawise Inc.
     VERSION: ---
     CREATED: 2020-11-20T20:08:01.309305
    REVISION: ---

==============================================================================="""

import numpy as np
import logging
import random


def _circumscribed_triangle(a, b, c):
    """
    input : a,b,c points in $\\mathbb{R}^3$ space, each as 3-tuple of coordinates
    """
    a, b, c = (np.array(x) for x in (a, b, c))
    normal_vector = np.cross(a-b, a-c)
    d1, d2 = (np.cross(normal_vector, a-x) for x in (b, c))
    p1, p2 = ((a+x)/2 for x in (b, c))
    m = np.array([d1, -d2])
    if np.linalg.det(m[:, :-1]) == 0:
        m = m[:, 1:]
        rhs = (p2-p1)[1:]
    else:
        m = m[:, :-1]
        rhs = (p2-p1)[:-1]
    if np.linalg.det(m) == 0:
        min_pt, max_pt = [np.array(
            [f([x[i] for x in [a, b, c]]) for i in range(3)]) for f in [np.min, np.max]]
        return list((min_pt+max_pt)/2)
    else:
        xy = np.linalg.solve(m.transpose(), rhs)
        return list(xy[0]*d1+p1)


def _s2_to_r3(theta, phi, r=1, phi_theta_in_degrees=False):
    if phi_theta_in_degrees:
        phi, theta = [x*np.pi/180.0 for x in [phi, theta]]
    assert 0 <= theta <= np.pi
    assert 0 <= phi <= 2*np.pi
    assert r >= 0
    # formula taken from https://en.wikipedia.org/wiki/Spherical_coordinate_system#Cartesian_coordinates
    return (r*np.sin(theta)*np.cos(phi), r*np.sin(theta)*np.sin(phi), r*np.cos(theta))


def _r3_to_s2(x, y, z, phi_theta_in_degrees=False):
    """
    return (theta,phi,r)
    """
    r = np.linalg.norm(np.array([x, y, z]))
    assert r > 0
    theta, phi = (np.arccos(z/r), np.arctan2(y, x))
    if phi_theta_in_degrees:
        theta, phi = [x*180.0/np.pi for x in [theta, phi]]
    logger = logging.getLogger("_r3_to_s2")
    logger.info(f"{(x,y,z)} => {theta,phi,r}")
    return (theta, phi, r)


def _circumscribed_triangle_on_s2(a, b, c, phi_theta_in_degrees=False):
    """
    input: a,b,c points on $\\mathbb{S}^2$, each as a 2-tuple (theta,phi) in radians or degrees
    return (theta,phi)
    """
    if phi_theta_in_degrees:
        a, b, c = [np.deg2rad(x) for x in [a, b, c]]
    res = _circumscribed_triangle(*[_s2_to_r3(*x) for x in [a, b, c]])
    theta, phi, _ = _r3_to_s2(*res)
    if phi_theta_in_degrees:
        theta, phi = [np.rad2deg(x) for x in [theta, phi]]
    return (theta, phi)


def __pseudo_s2_distance(a, b):
    p1, p2 = [np.array(_s2_to_r3(*x)) for x in [a, b]]
    return np.linalg.norm(p1-p2)


def _welzl(P, R):
    if len(P) == 0:
        if len(R) == 0:
            return None, -1
        elif len(R) == 1:
            return R[0], 0
        elif len(R) == 2:
            p1, p2 = [np.array(_s2_to_r3(theta, phi)) for theta, phi in R]
            res = _r3_to_s2(*((p1+p2)/2))[:2]
            return res, __pseudo_s2_distance(res, R[0])
        else:
            res = _circumscribed_triangle_on_s2(*R[:3])
            # FIXME: need to check other points?
            return res, __pseudo_s2_distance(res, R[0])
    pt_i = random.choice(range(len(P)))
    P_ = [pt for i, pt in enumerate(P) if i != pt_i]
    c, r = _welzl(P_, R)
    pt = P[pt_i]
    if c is not None and __pseudo_s2_distance(pt, c) <= r:
        return c, r
    else:
        return _welzl(P_, R+[pt])


def minimal_circle_on_s2(pts, phi_theta_in_degrees=False):
    """
    input: pts=[pt1,pt2,...ptN], where ptI=(theta,phi)
    return: (theta,phi)

    remark: lat = 90-theta, lon = phi
    """
    if phi_theta_in_degrees:
        pts = [np.deg2rad(pt) for pt in pts]
    (theta, phi), _ = _welzl(pts, [])
    if phi_theta_in_degrees:
        theta, phi = [np.rad2deg(x) for x in [theta, phi]]
    return theta, phi