Welzl 算法的迭代版本

Question

我正在使用 Welzl 的算法来查找点云的最小封闭圆 (2d) 或最小封闭球 (3d)。不幸的是，该算法具有非常高的递归深度，即输入点的数量。这个算法有迭代版本吗？我找不到任何东西，也不知道如何将递归更改为循环。

我发现了一些迭代的最小封闭圆/球算法，但它们的工作方式完全不同，并且没有 Welzl 算法的预期线性运行时间。

Answer 1

打乱点数组 P[0…n−1]。令 R = ∅ 且 i = 0。直到 i = n，

• 如果 P[i] ∈ R 或 R 定义了一个包含 P[i] 的圆，则设 i ← i+1。
• 否则，设置 R ← {P[i]} ∪ (R ∩ {P[i+1], ..., P[n−1]})。如果 |R|

返回 R。

在 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

下面的 Python 代码是从可用的 Rust 代码转换而来的；但是这个转换后的代码是未经测试的，所以认为它是伪代码。 get_circ_2_pts() 和 get_circ_3_pts() 是函数的存根，分别需要 2 和 3 个 Points 来获得相交的圆。这些公式应该很容易找到。

下面的代码模拟了编译后的本机代码在递归调用期间如何推送/恢复堆栈帧并将返回值传回给调用者等。

这不是将递归函数转换为迭代函数的好方法，但它确实有效。对于原生编译代码，它基本上是将函数的堆栈从堆栈内存转移到堆内存

在速度方面不会有太多回报。而且内存效率也没有任何回报，除非应用程序中实现的 Welzl 算法处理的点太多以至于占用了堆栈空间；在这种情况下，这可能会解决问题。

有关这种递归到迭代转换方法的通用示例，您可以查看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