我正在使用 Welzl 算法来查找点云的最小外接圆 (2d) 或最小外接球 (3d)。不幸的是，该算法具有非常高的递归深度，即输入点数。这个算法有迭代版本吗？我找不到任何并且不知道如何将递归更改为循环。
我发现了一些迭代最小外接圆/球算法，但它们的工作方式完全不同，并且没有 Welzl 算法预期的线性运行时间。
打乱点数组 P[0…n−1]。令 R = ∅ 且 i = 0。直到 i = n，
返回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()
Welzl 的 Python 3.10+ 算法。这使用显式堆栈而不是递归。下面只是
Circle @staticmethod
def from_points(points: list[Point]) -> 'Circle':
random.shuffle(points)
# Algorithm stages.
First = 0
Second = 1
stack = [(First, [], len(points))]
results = []
while stack:
params = stack.pop()
if params[0] is First:
r = params[1]
n = params[2]
if n == 0 or len(r) == 3:
match len(r):
case 3: c1 = Circle.from_3_points(r[0], r[1], r[2])
case 2: c1 = Circle.from_2_points(r[0], r[1])
case 1: c1 = Circle(r[0], 0)
case 0: c1 = Circle(Point(0, 0), 0)
case _: raise ValueError("Invalid number of points")
results.append(c1)
else:
idx = random.randint(0, n - 1)
p1 = points[idx]
points[idx], points[n - 1] = points[n - 1], points[idx]
stack.append((Second, r.copy(), n, p1))
stack.append((First, r, n - 1))
elif params[0] is Second:
r = params[1]
n = params[2]
p1 = params[3]
c1 = results.pop()
if c1.contains_point(p1):
results.append(c1)
else:
r.append(p1)
stack.append((First, r, n - 1))
else:
raise ValueError("Invalid stage")
return results.pop()
完整代码：
from dataclasses import dataclass
import random
@dataclass
class Point:
x: float
y: float
def distance(self, other: 'Point') -> float:
return ((self.x - other.x)**2 + (self.y - other.y)**2)**0.5
@dataclass
class Circle:
center: Point
radius: float
@staticmethod
def from_points(points: list[Point]) -> 'Circle':
"""Return the smallest circle that encloses the given points.
This function uses Welzl's algorithm to find the smallest circle that
encloses the given points.
Args:
points: list of points
Side effects:
The order of the points in the list is randomized.
Returns:
Circle: the smallest circle that encloses the given points.
"""
random.shuffle(points)
# Algorithm stages.
First = 0
Second = 1
stack = [(First, [], len(points))]
results = []
while stack:
params = stack.pop()
if params[0] is First:
r = params[1]
n = params[2]
if n == 0 or len(r) == 3:
match len(r):
case 3: c1 = Circle.from_3_points(r[0], r[1], r[2])
case 2: c1 = Circle.from_2_points(r[0], r[1])
case 1: c1 = Circle(r[0], 0)
case 0: c1 = Circle(Point(0, 0), 0)
case _: raise ValueError("Invalid number of points")
results.append(c1)
else:
idx = random.randint(0, n - 1)
p1 = points[idx]
points[idx], points[n - 1] = points[n - 1], points[idx]
stack.append((Second, r.copy(), n, p1))
stack.append((First, r, n - 1))
elif params[0] is Second:
r = params[1]
n = params[2]
p1 = params[3]
c1 = results.pop()
if c1.contains_point(p1):
results.append(c1)
else:
r.append(p1)
stack.append((First, r, n - 1))
else:
raise ValueError("Invalid stage")
return results.pop()
@staticmethod
def from_3_points(p1, p2, p3):
"""Calculate the circle passing through 3 points."""
p12 = Point(p2.x - p1.x, p2.y - p1.y)
p13 = Point(p3.x - p1.x, p3.y - p1.y)
c12 = p12.x * p12.x + p12.y * p12.y
c13 = p13.x * p13.x + p13.y * p13.y
c123 = p12.x * p13.y - p12.y * p13.x
cx = (p13.y * c12 - p12.y * c13) / (2. * c123)
cy = (p12.x * c13 - p13.x * c12) / (2. * c123)
center = Point(cx + p1.x, cy + p1.y)
return Circle(center, center.distance(p1))
@staticmethod
def from_2_points(p1: Point, p2: Point) -> 'Circle':
"""Return the circle passing through the two points."""
center_x = (p1.x + p2.x) / 2
center_y = (p1.y + p2.y) / 2
center = Point(center_x, center_y)
radius = center.distance(p1)
return Circle(center, radius)
def contains_point(self, point: Point) -> bool:
"""Return True if the circle contains the given point."""
return self.center.distance(point) <= self.radius
你可以这样理解迭代版本：
你可以用集合中的两个或三个点形成一个非平凡的最小外接圆；
发现位于该圆内部或之上的任何点都可以从进一步搜索中丢弃，因为它不能与这两个或三个形成更小的圆。
因此可以使用以下程序：
选择一个点；
A) 选择第二个点并构建包围这两个点的最小圆；
B）如果所有剩余点都是内部的，则完成；
取一个外点并构建包围这三个点的最小圆；如果这个圆经过两点，则转到 B
如果所有剩余点都是内部的，则完成；
取一个外部点并前往A。