Smallest circle problem on $\mathbb{S}^2$

Original problem

Smallest circle problem is a well-known problem in planar geometry, which consists of finding the circle surrounding the given
finite set of points in place having the smallest radius.

While this problem can be further generalized to $N$-dimensions (with circles replaced by $\mathbb{S}^{N-1}$ spheres), for now on we will solely
restrict ourselves to the planar case.

Some simple reasoning suggests that the solution exists and unique for any finite set of points on a plane. Now, how to find it?
Again, some reasoning and perhaps trying out some examples would suggest one that the minimal bounding circle should contain at least 2 points on its boundary, otherwise it could be minified further
(the black circle on the picture below contains only one point on its edge, hence it can be minified, so to become the blue circle).

Now, if minimal bounding circle contains exactly two points, these have to be located on the diameter of the circle, otherwise it can be minified again, as blue circle can be minified to red circle for
two points A and B on the picture below.

However, as we know from the planar geometry, three points uniquely determine the circumscribed circle. Hence, by iterating through all
2- and 3-point subsets of the original finite set of points, we can try out all of the candidates and thus select the smallest bounding circle. This naive algorithm has the complexity of $O(n^3)$ ($n$ being the
number of the points we want to bound) and can be further improved.

In fact, while it was conjectured that the best time complexity for this problem would be $O(n \log n)$, in 1983 Nimrod Megiddo has proposed the linear-time algorithm. In 1991, another linear-time algorithm
was proposed by Emo Welzl, this time the recursive one. Since the latter has particularly simple structure, I will describe it here

Welzl's algorithm starts with seemingly more difficult problem: for two finite sets $P$ and $R$ of points on plane it attempts to find a minimal circle $O=welzl(P,R)$
which contains $P$ and has $R$ on its boundary.
If $P$ is empty, the problem reduces to the search of circumscribed circle for a polygon, which is well-defined and well-known, since such circumscribed circle is basically determined by
arbitrary three points.

Otherwise, the recursion step is done by selecting random point $p\in P$. Then, exactly one of two things can happen:

• $p$ lies on the boundary of $O$, in which case $O=welzl(P-{p},R \cup {p})$
• $p$ lies inside $O$, in which case $O=welzl(P-{p},R)$

and the problem can be solved recursively. The termination of recursion is guaranteed by the fact that size of $P$ reduces on every recursion step

We reproduce the pseudocode below (taken from Wikipedia with slight modifications):

algorithm welzl is
    input: Finite sets P and R of points in the plane |R|≤ 3.
    output: Minimal disk enclosing P with R on the boundary.

    if P is empty then
        return trivial(R)
    choose p in P (randomly and uniformly)
    D := welzl(P - { p }, R)
    if p is in D then
        return D

    return welzl(P - { p }, R ∪ { p })

Here trivial(R) depends on a size of $R$:

• if $R$ is empty, it returns the imaginary "degenerate circle" of radius -1, which contains no points;
• if $R$ consists of one point $p$, the result is the circle with center in $p$ and radius 0;
• if $R$ consists of two points $a$ and $b$, the result is the circle with center located in the midpoint of $a$ and $b$ and with the radius equal half of the distance between $a$ and $b$ (so to have
  them on diameter)
• if $R$ consists of more than 3 points, random 3 of them are chosen and corresponding circumscribed circle $O$ is generated (in case 3 points chosen form degenerated triangle i.e. all three lie on one line,
  circumscribed circle is generated according to one of the previous two cases); if $O$ contains all points from $R$ on its boundary, it is returned as is, otherwise the degenerate circle is returned.

Setting on $\mathbb{S}^2$

Now, suppose that we are not on a plane, but rather on two-dimensional sphere $\mathbb{S}^2$. The smallest circle problem can be defined here as well, since circles and their radii can live on a sphere
as well (circle on $\mathbb{S}^2$ is basically an intersection of sphere and plane). The only concept which would require some rethinking is the concept of "lying inside the circle".
When we say "outside" in application to some closed curve (e.g. circle) we intuitively take without proof
the fact that "well-behaved" closed curve does split the $\mathbb{R}^2$ plane into two connected components, of which one is bounded (we call it "inside") and one is not (we call it "outside").
By the way, this is very nontrivial to prove strictly
and is known as Jordan curve theorem.

However, since sphere is
compact, the "inside" and "outside" here are technically both "inside". Nevertheless, we can rule out this difficulty by restricting ourselves to sets of points contained in only one hemisphere. Since
hemisphere without the equator has the same topology as $\mathbb{R}^2$ plane (i.e. they are homeomorphic, even diffeomorphic). In future we will always implicitly assume this condition.

Now, one would naturally want to adapt Welzl's algorithm to the current setting. Essentially, the only minor technical complication is the computation of circumscribed triangle on a sphere.

Luckily, however,
three points in space necessary lie in one plane. Moreover, three distinct points on a sphere will necessary lie on the unique plane (since no line can cross
a sphere in three distinct points).
And since as we mentioned above, the circle on a sphere is determined by intersection with the plane, the problem
essentially reduces to the planar case. It other words, we first compute the plane passing through the three points, and then simply compute the intersection of this plane with the sphere, which
gives us the circumscribed triangle.

The source code is available as a single python file: https://gist.github.com/nailbiter/f418a1952ed1af9ac35313febb0b0ba2.

One can call the function minimal_circle_on_s2 with the pts argument being the array of spherical coordinates of points on a sphere. By default, coordinates are assumed to be in radians, but if
degrees are more preferable, one can set the flag phi_theta_in_degrees to True. The output is the tuple $(theta,phi)$ of spherical coordinates of the bounding circle.