最小包围圆|设置 2 – Welzl 算法

简介

Welzl算法是用来寻找二维平面上一组点的最小圆覆盖的算法。它基于一个简单的观察结果：最小圆覆盖一定是这些点子集的某个Delaunay 三角形的外接圆或者是这个三角形的某条边上的圆。

算法流程

Welzl算法基于递归地构建最小圆覆盖。具体流程如下：

1. 从点集P中随机选出一个点,作为最小圆覆盖的中心点c
2. 遍历P，如果所有点都在以c为圆心，以其中一个点为半径的圆内，则已找到最小圆，结束
3. 随机选出一个点q不在以c为圆心以其中一个点为半径的圆内，则以cq为直径或者qt为直径（t为P中不在c和q所在圆内的点）求出新的最小覆盖圆
4. 对P中不在c和q所在圆内的点进行递归搜索，直到最后所有点都被包括在圆内

实现代码

Welzl算法的代码实现如下：

#include <bits/stdc++.h>
using namespace std;

const int N = 1e5 + 5;

struct Point {
    double x, y;
};

double dist(Point a, Point b) {
    return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}

Point center(Point a, Point b, Point c) {
    double x1 = c.x - a.x, y1 = c.y - a.y;
    double x2 = b.x - a.x, y2 = b.y - a.y;
    double x3 = b.x - c.x, y3 = b.y - c.y;
    double A = 2*(x1*y2 - x2*y1);
    double C = 2*(x2*y3 - x3*y2);
    double x = (y3*(x1*x1+y1*y1) + y1*(x2*x2+y2*y2) + y2*(x3*x3+y3*y3))/A/3.0;
    double y = (x1*(y2*y2+x2*x2) + x2*(y3*y3+x3*x3) + x3*(y1*y1+x1*x1))/C/3.0;
    return (Point){x, y};
}

double get_radius(Point a, Point b, Point c) {
    Point p = center(a, b, c);
    return dist(p, a);
}

bool is_in_circle(Point p, Point a, Point b, Point c) {
    return dist(p, a) <= get_radius(a, b, c);
}

void welzl_algorithm(Point* p, int n, Point* ans, int& num) {
    if (n == 0 || num == 3) {
        if (num == 1 || num == 2 && !is_in_circle(p[2], ans[0], ans[1], ans[1])) {
            ans[0] = p[0];
        } else if (num == 2 || num == 3 && !(is_in_circle(p[3], ans[0], ans[1], ans[2]) && is_in_circle(p[3], ans[1], ans[0], ans[2]) && is_in_circle(p[3], ans[2], ans[0], ans[1]))) {
            ans[0] = p[0];
            ans[1] = p[1];
        }
        num = min(num, 2);
        return ;
    }
    welzl_algorithm(p + 1, n - 1, ans, num);
    int t = num;
    for (int i = 0; i < t; i++) {
        if (dist(ans[i], p[0]) > get_radius(ans[i], p[0], ans[(i+1)%t])) {
            swap(ans[i], ans[--num]);
            i--;
            t--;
        }
    }
    if (!is_in_circle(p[0], ans[0], ans[1], ans[2])) {
        ans[num++] = p[0];
        assert(1 <= num && num <= 3);
        welzl_algorithm(p + 1, n - 1, ans, num);
        swap(ans[num-1], ans[0]);
    }
}

总结

Welzl算法是一种十分优美的算法，虽然涉及到的几何知识有些深，但是算法实现上并不难，利用递归的方式，可以在$O(n)$的时间复杂度内求出最小包围圆。