Welzl 算法的迭代版本 [英] Iterative version of Welzl's algorithm

问题描述

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

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.

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

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.

推荐答案

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

The Python code below is a conversion from working Rust code; but this converted code is untested, so consider it pseudocode. get_circ_2_pts() and get_circ_3_pts() are stubs for functions that take 2 and 3 Points respectively to get the intersecting circle. Formulae for these should be easy to find.

下面的代码模拟了已编译的本机代码在递归调用期间推送/恢复堆栈帧并将返回值传递回调用者等的行为.

The code below mimics how compiled native code behaves wrt pushing/restoring stack frames during recursive calls and passing return values back to the caller, etc.

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

This is not a good way to go about converting a recursive function to an iterative one, but it works. For native compiled code, it basically shifts the function's stack from stack memory to heap memory

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

There isn't going to be much payoff in terms of speed. And there also isn't any payoff in memory efficiency unless the Welzl algorithm implemented in an application is processing so many points that it's eating up stack space; in which case, this may address the problem.

有关这种递归到迭代转换方法的一般示例，您可以查看我关于此主题的其他帖子.

For a generalized example of this approach of recursive to iterative conversion, you can check out 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

这篇关于Welzl 算法的迭代版本的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！