Smallest-circle problem

Run of Megiddo's algorithm phase, discarding from point set A,B,...,U needless points E, T.

Megiddo's algorithm[8] is based on the technique called prune and search reducing size of the problem by removal of n/16 of unnecessary points. That leads to the recurrence t(n)≤ t(15n/16)+cn giving t(n)=16cn.

The algorithm is rather complicated and it is reflected in big multiplicative constant. The reduction needs to solve twice the similar problem where center of the sought-after enclosing circle is constrained to lie on a given line. The solution of the subproblem is either solution of unconstrained problem or it is used to determine the half-plane where the unconstrained solution center is located.

The n/16 points to be discarded are found the following way: Points P_i are arranged to pairs what defines n/2 lines p_j as their bisectors. Median p_m of bisectors in order by their directions (oriented to the same half-plane determined by bisector p₁) is found and pairs from bisectors are made, such that in each pair one bisector has direction at most p_m and the other at least p_m (direction p₁ could be considered as -${\displaystyle \infty }$ or +${\displaystyle \infty }$ according our needs.) Let Q_k be intersection of bisectors of k-th pair.

Line q in the p₁ direction is placed to go through an intersection Q_x such that there is n/8 intersections in each half-plane defined by the line (median position). Constrained version of the enclosing problem is run on line q what determines half-plane where the center is located. Line q' in the p_m direction is placed to go through an intersection Q_x' such that there is n/16 intersections in each half of the half-plane not containing the solution. Constrained version of the enclosing problem is run on line q' what together with q determines quadrant, where the center is located. We consider the points Q_k in the quadrant not contained in a half-plane containing solution. One of the bisectors of pair defining Q_k has the direction ensuring which of points P_i defining the bisector is closer to each point in the quadrant containing the center of the enclosing circle. This point could be discarded.

The constrained version of the algorithm is also solved by the prune and search technique, but reducing problem size by removal of n/4 points leading to recurrence t(n)≤ t(3n/4)+cn giving t(n)=4cn.

The n/4 points to be discarded are found the following way: Points P_i are arranged to pairs. For each pair intersection Q_j of its bisector with the constraining line q is found. (If intersection does not exist we could remove one point from the pair immediately). Median M of points Q_j on the line q is found and in O(n) time is determined which halfline of q starting in M contains the solution of the constrained problem. We consider points Q_j from the other half. We know which of points P_i defining Q_j is closer to the each point of the halfline containing center of the enclosing circle of the constrained problem solution. This point could be discarded.

The half-plane where the unconstrained solution lies could be determined by the points P_i on the boundary of the constrained circle solution. (First and last point on the circle in each half-plane suffice. If the center belongs to their convex hull, it is unconstrained solution, otherwise direction to the nearest edge determines half-plane of the unconstrained solution.)

Welzl described the algorithm in a recursive form. For size of input set at most 3, trivial algorithm is created. It outputs empty circle, circle with radius 0, circle on diameter made by arc of input points for input sizes 0, 1 and 2 respectively. For input size 3 if the points are vertices of an acute triangle, the circumscribed circle is returned. Otherwise the circle having a diameter of the triangle's longest edge is returned. For bigger sizes the algorithm uses solution of 1 smaller size where the ignored point p is chosen randomly and uniformly. The returned circle D is checked and if it encloses p, it is returned as result. Otherwise we know point p is on border of the result. We rerun the algorithm again, but we call it with the information p is on the result circle boundary. Points which are known to be on the boundary are accumulated in the set R and they are excluded from the subset P from which random choices are made in the following calls.

The condition to use the trivial algorithm is |P|=0 or |R|=3, and the trivial algorithm is called for R. Equivalent condition would be |R|=3 or |P|+|R|≤3, and call of the trivial algorithm for R in the case |R|=3 and for union of P and R otherwise. This would prevent some recursive calls on bottom of the recursion.

algorithm welzl is[7]
    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 or |R| = 3 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 })

Welzl indicated implementation without need of recursion. Empty R and random permutation of P is used at start. Sets on bottom of recursion correspond to prefixes in the permutation. Trivial algorithm is called for R extended by prefix of P to size of 3, and following points are checked they are enclosed by its result D. Whenever check fails on point p, p is added to R and removed from P, and the prefix of the permutation up to p is reshufled. To make this implementation equivalent to the recursive one, we should maintain stack of prefix lengths and if longer prefix than the top would be inserted, last point from R should be removed and inserted to P, the top prefix on the stack should be removed at the same time. Therefore, the prefix length on stack would be nondecreasing from its top down. When |R|=3, the points up to prefix length on top position of the stack are not checked.

Typical implementation uses one list(array) of points of R ∪ P with R in its start, remembering the size of R, therefore sizeR could be used as stack pointer and its decrements while top prefix length was smaller than the current prefix length to be inserted does all required actions.

Welzl also proposes functionally different version of the algorithm using random permutation of points on input, where the points are not reshuffled when the point p is moved to front of the list. The algorithm is again presented in the recursive form.

Example of Welzl's algorithm run which returns wrong result

Matoušek, Sharir, Welzl [6] implicitly warned that the assumption point moved to R remains on enclosing circle's boundary for all subsets of R ∪ P does not hold, and the required additional condition is the radius of the enclosing circle does not decrease.

They proposed tiny change to the algorithm forcing the radius of the enclosing circle does not decrease. If p is not in enclosing circle D of R ∪ P - { p }, the radius r of enclosing circle of D must be smaller than radius r' of enclosing circle D' of R ∪ P. But radius of an enclosing circle of any subset of R ∪ P - { p } is at most r, so whenever r*>r is radius of an enclosing circle of a subset of R ∪ P, p must be in its boundary.

MSW algorithm chooses 3 random points of P as initial base R and they find initial D by running the trivial algorithm on them. The presented recursive version of the algorithm is very similar to Welzl's. Point p of P is chosen in random order and algorithm for the set P - { p } is called. If p is not enclosed by returned circle D, p replaces one of 3 points in the base (R ∪ P does not change, and |R| remains 3), such that the enclosing circle of the new base encloses the point removed from the base. The result recursion is restarted with changed P, R.

algorithm msw is
    input: Finite sets P and R of points in the plane |R|=3
    output: Minimal disk enclosing P ∪ R
    if P is empty then
        return trivial(R)
    choose p in P (randomly and uniformly)
    D := msw(P - { p }, R)
    if p is in D then
        return D
    q := nonbase(R ∪ { p }) (Welzl's algorithm for 4 points could be used to find what would not be in R)
    return msw(P - { p } ∪ { q }, R ∪ { p } - { q })

Use of base extended by one point is mentioned in the paper. If we put p to the extended base replacing extension, we could use Welzl's recurrence to rearrange the extended base to base and the extra point. You can look at it as allowing temporary decrease of radius r which would be restored when the extended base is finished. The move to front version of Welzl's algorithm works according to this schema except the randomization of choices. This is why it returns correct results, on the contrary to the Welzl's original algorithm.

Move to front version of welzl without recursion and one list for R ∪ P in view of msw algorithm could be simplified by not maintaining the size of R at all.

Experiments indicate the move to front version have performance slightly better (say 2% less checks p belongs to D) than the rerandomized one (of msw).