📜  最小围圈设置2 –韦尔兹(Welzl)的算法(1)

📅  最后修改于: 2023-12-03 15:26:27.528000             🧑  作者: Mango

最小围圈问题

最小围圈问题是一个计算几何问题,给定平面上的一些点,要求找到一个最小的圆,使得所有给定的点都在圆内或圆上。这个问题在计算机视觉、数据挖掘等领域得到了广泛应用,例如在图像处理中可以用来检测目标物体的位置和形状。

韦尔兹算法

韦尔兹算法是求解最小围圈问题的一种常用算法,其时间复杂度为 $O(n)$。该算法的思路是递归地缩小点集合,找出最小圆。

算法步骤
  1. 如果点集合为空,返回一个包含无穷远点的圆;
  2. 如果点集合只有一个点,则返回一个半径为零且圆心为该点的圆;
  3. 随机选择一个点 $p$,将其从点集中移除,递归求解点集 $S$ 的最小圆 $C$;
  4. 如果点 $p$ 在圆 $C$ 中,直接返回圆 $C$;
  5. 否则,将点 $p$ 加入点集 $S$ 中,递归求解点集 $S$ 的最小圆 $C'$;
  6. 如果点 $p$ 在圆 $C'$ 中,直接返回圆 $C'$;
  7. 否则,点 $p$ 必须在圆外。根据 Welzl 算法的关键,我们需要找到一个最小圆 $D$,使得圆 $D$ 包含点 $p$、圆 $C$、圆 $C'$ 中的所有点;
  8. 返回圆 $D$。
代码实现

下面是使用 Python 语言实现的代码:

import random
from typing import List, Tuple

def min_circle(points: List[Tuple[float, float]]) -> Tuple[Tuple[float, float], float]:
    def recursive_min_circle(points: List[Tuple[float, float]], circumference_points: List[Tuple[float, float]]):
        if not points or len(circumference_points) == 3:
            cx, cy, r = trivial_min_circle(circumference_points)
            return (cx, cy), r
        p = points.pop()
        circle = recursive_min_circle(points, circumference_points)
        if is_point_in_circle(circle, p):
            return circle
        circumference_points.append(p)
        return recursive_min_circle(points, circumference_points)

    def trivial_min_circle(points: List[Tuple[float, float]]) -> Tuple[float, float, float]:
        if len(points) == 0:
            return (0, 0), 0
        elif len(points) == 1:
            return points[0], 0
        elif len(points) == 2:
            (x0, y0), (x1, y1) = points
            cx, cy = (x0 + x1) / 2, (y0 + y1) / 2
            r = ((x0 - cx) ** 2 + (y0 - cy) ** 2) ** 0.5
            return (cx, cy), r
        else:
            x0, y0 = points[0]
            x1, y1 = points[1]
            x2, y2 = points[2]
            p = x1 - x0, y1 - y0
            q = x2 - x0, y2 - y0
            u = 2 * (p[0] * q[1] - p[1] * q[0])
            if u == 0:
                return (0, 0), 0
            cx = (q[1] * (p[0] ** 2 + p[1] ** 2) - p[1] * (q[0] ** 2 + q[1] ** 2)) / u
            cy = (p[0] * (q[0] ** 2 + q[1] ** 2) - q[0] * (p[0] ** 2 + p[1] ** 2)) / u
            r = ((x0 - cx) ** 2 + (y0 - cy) ** 2) ** 0.5
            return (cx, cy), r

    def is_point_in_circle(circle: Tuple[Tuple[float, float], float], point: Tuple[float, float]) -> bool:
        (cx, cy), r = circle
        x, y = point
        return ((x - cx) ** 2 + (y - cy) ** 2) <= r ** 2

    random.shuffle(points)
    return recursive_min_circle(points, [])

参考文献