Iterative version of Welzl's algorithm

Question

I am using Welzl's algorithm to find the smallest enclosing circle (2d) or smallest enclosing sphere (3d) of a point cloud. Unfortunately, the algorithm has a very high recursion depth, namely the number of input points. Is there an iterative version of this algorithm? I could not find any and have no idea how the recursion can be changed to a loop.

I found some iterative smallest enclosing circle/sphere algorithms, but they work completely different and do not have the expected linear runtime of Welzl's algorithm.

Answer 1

Shuffle the point array P[0…n−1]. Let R = ∅ and i = 0. Until i = n,

• If P[i] ∈ R or R defines a circle that contains P[i], then set i ← i+1.
• Otherwise, set R ← {P[i]} ∪ (R ∩ {P[i+1], …, P[n−1]}). If |R| < 3, then set i ← 0, else set i ← i+1.

Return R.

Tested implementation in Python 3:

from itertools import combinations
from random import randrange, shuffle


def mag_squared(v):
    return v.real ** 2 + v.imag ** 2


def point_in_circle(p, circle):
    if not circle:
        return False
    if len(circle) == 1:
        (q,) = circle
        return q == p
    if len(circle) == 2:
        q, r = circle
        return mag_squared((p - q) + (p - r)) <= mag_squared(q - r)
    if len(circle) == 3:
        q, r, s = circle
        # Translate p to the origin.
        q -= p
        r -= p
        s -= p
        # Orient the triangle counterclockwise.
        qr = r - q
        qs = s - q
        a, b = qr.real, qr.imag
        c, d = qs.real, qs.imag
        determinant = a * d - b * c
        assert determinant != 0
        if determinant < 0:
            r, s = s, r
        # Test whether the origin lies in the circle.
        a, b, c = q.real, q.imag, mag_squared(q)
        d, e, f = r.real, r.imag, mag_squared(r)
        g, h, i = s.real, s.imag, mag_squared(s)
        determinant = a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g)
        return determinant >= 0
    assert False, str(circle)


def brute_force(points):
    for k in range(4):
        for circle in combinations(points, k):
            if any(
                point_in_circle(p, circle[:i] + circle[i + 1 :])
                for (i, p) in enumerate(circle)
            ):
                continue
            if all(point_in_circle(p, circle) for p in points):
                return circle
    assert False, str(points)


def welzl_recursive(points, required=()):
    points = list(points)
    if not points or len(required) == 3:
        return required
    p = points.pop(randrange(len(points)))
    circle = welzl_recursive(points, required)
    if point_in_circle(p, circle):
        return circle
    return welzl_recursive(points, required + (p,))


def welzl(points):
    points = list(points)
    shuffle(points)
    circle_indexes = {}
    i = 0
    while i < len(points):
        if i in circle_indexes or point_in_circle(
            points[i], [points[j] for j in circle_indexes]
        ):
            i += 1
        else:
            circle_indexes = {j for j in circle_indexes if j > i}
            circle_indexes.add(i)
            i = 0 if len(circle_indexes) < 3 else i + 1
    return [points[j] for j in circle_indexes]


def random_binomial():
    return sum(2 * randrange(2) - 1 for i in range(100))


def random_point():
    return complex(random_binomial(), random_binomial())


def test(implementation):
    for rep in range(1000):
        points = [random_point() for i in range(randrange(20))]
        expected = set(brute_force(points))
        for j in range(10):
            shuffle(points)
            got = set(implementation(points))
            assert expected == got or (
                all(point_in_circle(p, expected) for p in got)
                and all(point_in_circle(p, got) for p in expected)
            )


def graphics():
    points = [random_point() for i in range(100)]
    circle = welzl(points)
    print("%!PS")
    print(0, "setlinewidth")
    inch = 72
    print(8.5 * inch / 2, 11 * inch / 2, "translate")
    print(inch / 16, inch / 16, "scale")
    for p in points:
        print(p.real, p.imag, 1 / 4, 0, 360, "arc", "stroke")
    for p in circle:
        print(p.real, p.imag, 1 / 4, 0, 360, "arc", "fill")
    print("showpage")


test(welzl_recursive)
test(welzl)
graphics()

Answer 2

The Python code below is a conversion from working Rust code; but this converted code is untested, so consider it pseudocode. get_circ_2_pts() and get_circ_3_pts() are stubs for functions that take 2 and 3 Points respectively to get the intersecting circle. Formulae for these should be easy to find.

The code below mimics how compiled native code behaves wrt pushing/restoring stack frames during recursive calls and passing return values back to the caller, etc.

This is not a good way to go about converting a recursive function to an iterative one, but it works. For native compiled code, it basically shifts the function's stack from stack memory to heap memory

There isn't going to be much payoff in terms of speed. And there also isn't any payoff in memory efficiency unless the Welzl algorithm implemented in an application is processing so many points that it's eating up stack space; in which case, this may address the problem.

For a generalized example of this approach of recursive to iterative conversion, you can check out my other post on this topic.

from random import randint
from copy   import deepcopy

class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y
        
    def distance(self, p):
        ((p.x - self.x)**2 + (p.y - self.y)**2)**.5
        
class Circle:
    def __init__(self, center, radius):
        self.center = center
        self.radius = radius

    def contains_point(self, p):
        return self.center.distance(p) <= self.radius


class Frame:
    def __init__(self, size, stage='FIRST'):
        self.stage = stage
        self.r     = []
        self.n     = size
        self.p     = Point(0, 0)
        self.ret   = None
        
    def clone(self):
        return deepcopy(self)

def welzl(points):
    circle = None
    stack  = []
    stack.append(Frame(len(points)))
    
    while len(stack):
        frame = stack.pop()
        
        if frame.n == 0 or len(frame.r) == 3:
            r_len = len(frame.r)
            ret   = None
            
            if   r_len == 0: ret = Circle(Point(0, 0), 0)
            elif r_len == 1: ret = Circle(frame.r[0], 0)
            elif r_len == 2: ret = get_circ_2_pts(*frame.r[:2])
            elif r_len == 3: ret = get_circ_3_pts(*frame.r[:3])
                
            if len(stack):
                stack[-1].ret = ret
            else:
                circle = ret
        else:
            if frame.stage == 'FIRST':
                n           = frame.n - 1
                idx         = randint(0, n)
                frame.p     = points[idx]
                points[idx] = points[n]
                points[n]   = frame.p
                call        = frame.clone()
                                
                frame.stage = 'SECOND'
                stack.append(frame)
                call.n      = n
                call.stage  = 'FIRST'
                stack.append(call)
                
            elif frame.stage == 'SECOND':
                if frame.ret.contains_point(frame.p):
                    if len(stack):
                        stack[-1].ret = frame.ret
                    else:
                        circle = frame.ret
                else:
                    frame.stage = 'THIRD'
                    stack.append(frame)                    
                    call        = frame.clone()
                    call.n     -= 1
                    call.stage  = 'FIRST'
                    call.r.append(frame.p)
                    stack.append(call)
                
            elif frame.stage == 'THIRD':
                if len(stack):
                    stack[-1].ret = frame.ret
                else:
                    circle = frame.ret
                
    return circle