ZOJ 1450 Minimal Circle_zoj1450-优快云博客

本文链接：https://blog.youkuaiyun.com/njcithy/article/details/117390044

本文解析了随机增量算法的核心思想，包括随机打乱点集、构造最小包围圆并动态调整的过程。重点介绍了核心代码实现和计算圆心算法，以及其在处理随机点集上的高效性。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

问题链接：https://zoj.pintia.cn/problem-sets/91827364500/problems/91827364949

随机增量算法主要思想：

随机打乱点集中所有点的顺序；
以序列中前2个点构造最小包围圆；
一次搜索其余的点，若某个点v在D外，则重新构造最小包围圆D'，使D'包含D中的所有点且以v为边界点。重复执行第三步，直到求出最小包围圆。

这个算法的关键在于第 3 步，以 v 作为边界点构造最小包围圆，从而增加了构造约束，减少了自由度。该算法的内存耗用少，速度快。在随机点集的条件下，期望运行时间为 O(n)。但是，若原始输入序列是有序的，则打乱点集的顺序也是需要耗费时间的。

核心代码：

void min_cover_circle(point &c, double &r, int n) {
	random_shuffle(p + 1, p + 1 + n);
	c = p[1];
	r = 0;
	for (int i = 2; i <= n; ++i) {
		if (dis(c, p[i] ) > r + eps) {
			c = p[i];
			r = 0;
			for (int j = 1; j < i; ++j) {
				if (dis(c, p[j]) > r + eps) {
					c.x = (p[i].x + p[j].x) / 2;
					c.y = (p[i].y + p[j].y) / 2;
					r = dis(c, p[i]);
					for (int k = 1; k < j; ++k) {
						if (dis(c, p[k]) > r + eps) {
							// 根据点p[i],p[j],p[k]计算包围圆圆心
							c = circumcenter(p[i], p[j], p[k]);
							r = dis(c, p[i]);
						}
					}
				}
			}
		}
	}
}

计算圆心算法：

设圆坐标为 $(x-a)^2+(y-b)^2=r^2$ ，且过点 $(x_1, y_1),(x_2,y_2),(x_3,y_3 )$ ，则圆心c坐标为

$c_x = \frac{C_1\cdot B_2-C_2\cdot B_1 }{A_1\cdot B_2 - A_2\cdot B_1}$ ， $c_y = \frac{A_1\cdot C_2-A_2\cdot C_1 }{A_1\cdot B_2 - A_2\cdot B_1}$

其中

$A_1=2\cdot (x_2-x_1)$ ， $B_1 = 2\cdot (y_2-y_1)$ ， $C_1 = x_2^2+y_2^2-x_1^2-y_2^2$

$A_2=2\cdot (x_3-x_2)$ ， $B_2 = 2\cdot (y_3-y_2)$ ， $C_2 = x_3^2+y_3^2-x_2^2-y_2^2$

point circumcenter(point A, point B, point C) {
	point ans;

	//double a1 = 2 * (B.x - A.x), b1 = 2 * (B.y - A.y), c1 = pow(B.x, 2) + pow(B.y, 2) - pow(A.x, 2) - pow(A.y, 2);
	//double a2 = 2 * (C.x - B.x), b2 = 2 * (C.y - B.y), c2 = pow(C.x, 2) + pow(C.y, 2) - pow(B.x, 2) - pow(B.y, 2);
	//double d = a1 * b2 - a2 * b1;
	//ans.x = (c1 * b2 - c2 * b1) / d;
	//ans.y = (a1 * c2 - a2 * c1) / d;

	double a1 = B.x - A.x, b1 = B.y - A.y, c1 = (pow(a1, 2) + pow(b1, 2)) / 2;
	double a2 = C.x - A.x, b2 = C.y - A.y, c2 = (pow(a2, 2) + pow(b2, 2)) / 2;
	double d = a1 * b2 - a2 * b1;
	ans.x = A.x + (c1 * b2 - c2 * b1) / d;
	ans.y = A.y + (a1 * c2 - a2 * c1) / d;
	
	return ans;
}

完整的代码：

#include <iostream>
#include <algorithm>
#include <math.h>

using namespace std;

const int maxn = 500;
const double eps = 1e-8;

struct point {
	double x, y;
	point(double _x = 0, double _y = 0) {
		x = _x; y = _y;
	}
	point operator + (const point &p) {
		return point(x + p.x, y + p.y);
	}
	point operator - (const point &p) {
		return point(x - p.x, y - p.y);
	}
};

double dis(point p, point q) {
	return pow(pow((p.x - q.x), 2) + pow((p.y - q.y), 2), 0.5);
}

point p[maxn];

point circumcenter(point A, point B, point C) {
	point ans;

	double a1 = B.x - A.x, b1 = B.y - A.y, c1 = (pow(a1, 2) + pow(b1, 2)) / 2;
	double a2 = C.x - A.x, b2 = C.y - A.y, c2 = (pow(a2, 2) + pow(b2, 2)) / 2;
	double d = a1 * b2 - a2 * b1;
	ans.x = A.x + (c1 * b2 - c2 * b1) / d;
	ans.y = A.y + (a1 * c2 - a2 * c1) / d;

	return ans;
}

void min_cover_circle(point &c, double &r, int n) {
	random_shuffle(p + 1, p + 1 + n);
	c = p[1];
	r = 0;
	for (int i = 2; i <= n; ++i) {
		if (dis(c, p[i] ) > r + eps) {
			c = p[i];
			r = 0;
			for (int j = 1; j < i; ++j) {
				if (dis(c, p[j]) > r + eps) {
					c.x = (p[i].x + p[j].x) / 2;
					c.y = (p[i].y + p[j].y) / 2;
					r = dis(c, p[i]);
					for (int k = 1; k < j; ++k) {
						if (dis(c, p[k]) > r + eps) {
							// 由p[i],p[j],p[k]三点计算包围圆圆心
							c = circumcenter(p[i], p[j], p[k]);
							r = dis(c, p[i]);
						}
					}
				}
			}
		}
	}
}

int main() {
	int n;
	point c;
	double r;

	while (scanf("%d", &n) != EOF && n) {
		for (int i = 1; i <= n; i++) {
			scanf("%lf %lf", &p[i].x, &p[i].y);
		}
		min_cover_circle(c, r, n);
		printf("%.2f %.2f %.2f\n", c.x, c.y, r);
	}
	return 0;
}

参考：《离散点集最小包围圆算法分析与改进》