凸包算法——Divide and Conquer

Divide and Conquer

• 步骤:

  1. 分割点集: 将点集分为两个子集。
  2. 递归计算: 对每个子集递归计算凸包。
  3. 合并结果: 合并两个子集的凸包，得到最终的凸包。

• 时间复杂度: O(n log n)

初始代码：

from __future__ import division
from numpy import *

link = lambda a, b: concatenate((a, b[1:]))
edge = lambda a, b: concatenate(([a], [b]))

def qhull2D(sample):
    def dome(sample, base):
        h, t = base
        dists = dot(sample-h, dot(((0, -1), (1, 0)), (t-h)))
        outer = sample[dists > 0]
        if len(outer):
            pivot = sample[argmax(dists)]
            return link(dome(outer, edge(h, pivot)),
                        dome(outer, edge(pivot, t)))
        else:
            return base

    if len(sample) > 2:
        axis = sample[:, 0]
        base = take(sample, [argmin(axis), argmax(axis)], 0)
        return link(dome(sample, base), dome(sample, base[::-1]))
    else:
        return sample

但是在进行运算的时候，会出现无限递归的情况，排查后，增加递归深度检查，防止过深的递归调用；在 dome 函数中添加了处理 dists 值过小的情况，避免在精度不足时进入死循环;确保每次递归处理的数据都是有意义的，避免重复计算相同的数据；重构了代码如下：

from __future__ import division
from numpy import *

def link(a, b):
    return concatenate((a, b[1:]))

def edge(a, b):
    return concatenate(([a], [b]))

def qhull2D(sample):
    def dome(sample, base, depth=0, max_depth=1000):
        h, t = base
        dists = dot(sample - h, dot(((0, -1), (1, 0)), (t - h)))
        
        if len(dists) == 0:
            return base

        # Handle cases where all distances are very small
        if all(abs(dists) < 1e-10):
            return base

        outer = sample[dists > 0]

        # Handle empty outer case
        if len(outer) == 0:
            return base

        pivot_index = argmax(dists)
        pivot = sample[pivot_index]

        # Ensure depth does not exceed maximum allowed
        if depth > max_depth:
            return base

        # Recursive case
        left_dome = dome(outer, edge(h, pivot), depth + 1, max_depth)
        right_dome = dome(outer, edge(pivot, t), depth + 1, max_depth)

        return link(left_dome, right_dome)

    if len(sample) > 2:
        axis = sample[:, 0]
        base = take(sample, [argmin(axis), argmax(axis)], axis=0)
        left_dome = dome(sample, base)
        right_dome = dome(sample, base[::-1])
        
        return link(left_dome, right_dome)
    else:
        return sample