Welzl算法的迭代版本

Question

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

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

Answer 1

对点数组P0…进行混洗n−1.设R=∅，I= 0。直到i= n，

如果Pi集合R或R定义了包含Pi的圆，则set I i+1.

• Otherwise，∈R←{Pi}←(R∪{∩，…，Pn−1})。如果|R| < 3，则set I←0，否则set I←。

返回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个Point来获得相交圆的函数的存根。这些公式应该很容易找到。

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

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

在速度方面不会有太大的回报。而且在内存效率方面也没有任何回报，除非在应用程序中实现的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