algorithm - Welzl 算法的迭代版本

Question

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

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

标签： algorithm

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| < 3，然后设置 i ← 0，否则设置 i ← i+1。

返回 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()