Welzl算法求最小圆覆盖

给定平面上的$n$个点，求一个最小的圆使得所有的点都落在圆内。

求解上述问题便是$Welzl$算法的作用。其期望复杂度可以达到$O(n)$。

$Welzl$算法的正确性基于一个定理：
定理：对于平面上任意$n$个点，在其最小覆盖圆外取第$n+1$个点。那么第$n+1$个点一定在这$n+1$个点的最小覆盖圆上。

有了这个定理，就可以进行$Welzl$算法了。先取3个点建立一个圆(不共线的三个点能确定一个圆，如果共线就取距离最远的两个点作为直径建立圆)，然后遍历剩下的所有点，对于一个遍历到的点$k$：

• 如果该点在圆内，那么最小覆盖圆不变。
• 如果该点在圆外，根据上述定理，该点一定在想要求得的最小覆盖圆上，又因为三个点才能确定一个圆，所以需要枚举$k$之前的点来找剩下的两个点。直到找到两个点，使得与$k$确定的圆能够将$k$之前的所有点包在里面(或者在圆上)。

写成代码就是：

// a[] 存放需要求最小圆覆盖的所有点，n 为点的数量
// radius 为圆半径，center为圆心，都是全局变量
// point_in 函数判断一个点是否在最小覆盖圆内
// circle_center 函数有两个重载，一个用两个点确定圆，一个用三个点确定圆
// dist 函数计算距离
void min_circle_cover(point a[], int n)
{
	radius = 0;
	center = a[0];
	for (int i = 1; i < n; i++) {
		if (!point_in(a[i])) {
			center = a[i];
			radius = 0;
			for (int j = 0; j < i; j++) {
				if (!point_in(a[j])) {
					circle_center(a[i], a[j], center);
					radius = dist(a[j], center);
					for (int k = 0; k < j; k++) {
						if (!point_in(a[k])) {
							circle_center(a[i], a[j], a[k], center);
							radius = dist(a[k], center);
						}
					}
				}
			}
		}
	}
}

关于这个期望线性复杂度，我并不会证明。

只要调用$min\_circle\_cover$之前调用$random\_shuffle$随机排列一下，就可以保证期望线性复杂度。

这里给一个完整代码（洛谷p1742的AC代码）：

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

const int maxn = 1e5 + 10;

const double eps = 1e-8;
int cmp(double x) {
	if (fabs(x) < eps) return 0;
	if (x > 0) return 1;
	return -1;
}

struct point {
	double x, y;
	point() {}
	point(double a, double b) : x(a), y(b) {}
	friend point operator-(const point& a, const point& b) {
		return point(a.x - b.x, a.y - b.y);
	}
	double norm() {
		return sqrt(pow(x, 2) + pow(y, 2));
	}
}p[maxn];
double dist(const point& a, const point& b) {
	return (a - b).norm();
}

void circle_center(point p0, point p1, point& cp) {
	cp.x = (p0.x + p1.x) / 2;
	cp.y = (p0.y + p1.y) / 2;
}

void circle_center(point p0, point p1, point p2, point& cp) {
	double a1 = p1.x - p0.x, b1 = p1.y - p0.y, c1 = (a1 * a1 + b1 * b1) / 2;
	double a2 = p2.x - p0.x, b2 = p2.y - p0.y, c2 = (a2 * a2 + b2 * b2) / 2;
	double d = (a1 * b2 - a2 * b1);
	cp.x = p0.x + (c1 * b2 - c2 * b1) / d;
	cp.y = p0.y + (a1 * c2 - a2 * c1) / d;
}

double radius; point center;
bool point_in(const point& p) {
	return cmp(dist(p, center) - radius) < 0;
}

void min_circle_cover(point a[], int n)
{
	radius = 0;
	center = a[0];
	for (int i = 1; i < n; i++) if (!point_in(a[i])) {
		center = a[i], radius = 0;
		for (int j = 0; j < i; j++) if (!point_in(a[j])) {
			circle_center(a[i], a[j], center);
			radius = dist(a[j], center);
			for (int k = 0; k < j; k++) if (!point_in(a[k])) {
				circle_center(a[i], a[j], a[k], center);
				radius = dist(a[k], center);
			}
		}
	}
}

int main()
{
	int n;
	scanf("%d", &n);
	for (int i = 0; i < n; i++)
	{
		double x, y;
		scanf("%lf%lf", &x, &y);
		p[i] = point(x, y);
	}
	random_shuffle(p, p + n);
	min_circle_cover(p, n);
	printf("%.10lf\n%.10lf %.10lf", radius, center.x, center.y);

	return 0;
}

例题：
洛谷P1742 最小圆覆盖 [传送门]
洛谷P2533 [AHOI2012]信号塔 [传送门]