Menu
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.
558
Shuffle the point array P[0…n−1]. Let R = ∅ and i = 0. Until i = n,
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()
Welzl's algorithm for Python 3.10+. This uses an explicit stack instead of recursion. Below is just one method of the Circle class. The full code is in the next code view below it.
One benefit of converting recursive code to stack based iterative code is the execution stack won't be eaten up for large datasets that involve a lot of nested recursion to process. You can find a detailed description of the methodology for converting an algorithm from recursive to iterative here.
Welzl's algorithm from Circle class:
@staticmethod
def from_points(points: List[Point]) -> 'Circle':
"""Return the smallest circle that encloses the given points."""
# Algorithm stages.
FIRST = 0
SECOND = 1
stack = [(FIRST, ([], len(points)))]
results = []
while stack:
stage, params = stack.pop()
if stage is FIRST:
r, n = params
if n == 0 or len(r) == 3:
lr = len(r)
if lr == 3: c1 = Circle.from_3_points(*r)
elif lr == 2: c1 = Circle.from_2_points(*r)
elif lr == 1: c1 = Circle(r[0], 0)
elif lr == 0: c1 = Circle(Point(0, 0), 0)
else: 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 stage is SECOND:
r, n, p1 = params
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()
Circle class full code:
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."""
# Algorithm stages.
FIRST = 0
SECOND = 1
stack = [(FIRST, ([], len(points)))]
results = []
while stack:
stage, params = stack.pop()
if stage is FIRST:
r, n = params
if n == 0 or len(r) == 3:
lr = len(r)
if lr == 3: c1 = Circle.from_3_points(*r)
elif lr == 2: c1 = Circle.from_2_points(*r)
elif lr == 1: c1 = Circle(r[0], 0)
elif lr == 0: c1 = Circle(Point(0, 0), 0)
else: 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 stage is SECOND:
r, n, p1 = params
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: Point, p2: Point, p3: Point) -> 'Circle':
"""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
For reference, here's the recursive version of Welzl's:
@staticmethod
def from_points_recursive(points: List[Point]) -> 'Circle':
return Circle.welzl_recursive(points, [], len(points))
@staticmethod
def welzl_recursive(points: List[Point], r: List[Point], n: int) -> 'Circle':
if n == 0 or len(r) == 3:
lr = len(r)
if lr == 3: c1 = Circle.from_3_points(*r)
elif lr == 2: c1 = Circle.from_2_points(*r)
elif lr == 1: c1 = Circle(r[0], 0)
elif lr == 0: c1 = Circle(Point(0, 0), 0)
else: raise ValueError("Invalid number of points")
return c1
idx = random.randint(0, n - 1)
p1 = points[idx]
points[idx], points[n - 1] = points[n - 1], points[idx]
c1 = Circle.welzl_recursive(points, r.copy(), n - 1)
if c1.contains_point(p1):
return c1
r.append(p1)
return Circle.welzl_recursive(points, r, n - 1)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
