本文介绍了一种计算给定点集最小覆盖圆的算法，包括详细步骤、实现思路及代码示例。该算法逐步增加约束条件，从两个点扩展到三个点，最终确定能够覆盖所有点的最小圆。
链接：
You are to write a program to find a circle which covers a set of points and has the minimal area. There will be no more than 100 points in one problem.
Input
The input contains several problems. The first line of each problem is a line containing only one integer N which indicates the number of points to be covered. The next N lines contain N points. Each point is represented by x and y coordinates separated by a space. After the last problem, there will be a line contains only a zero.
Output
For each input problem, you should give a one-line answer which contains three numbers separated by spaces. The first two numbers indicate the x and y coordinates of the result circle, and the third number is the radius of the circle. (use escape sequence %.2f)
Sample Input
2
0.0 0.0
3 0
5
0 0
0 1
1 0
1 1
2 2
0
Sample Output
1.50 0.00 1.50
1.00 1.00 1.41
Source: Asia 1997, Shanghai (Mainland China)
题意：
算法：计算几何
思路：
这里介绍的算法是，先任意选取两个点，以这两个点的连线为直径作圆。再以此判断剩余的点，看它们是否都在圆内（或圆上），如果都在，说明这个圆已经找到。如果没有都在：假设我们用的最开始的两个点为p1,p[2]，并且找到的第一个不在圆内（或圆上）的点为p[i]，于是我们用这个点p[i]去寻找覆盖p1到p[i-1]的最小覆盖圆。
那么，过确定点p[i]的从p1到p[i-1]的最小覆盖圆应该如何求呢？
我们先用p1和p[i]做圆，再从2到i-1判断是否有点不在这个圆上，如果都在，则说明已经找到覆盖1到i-1的圆。如果没有都在：假设我们找到第一个不在这个圆上的点为p[j]，于是我们用两个已知点p[j]与p[i]去找覆盖1到j-1的最小覆盖圆。
而对于两个已知点p[j]与p[i]求最小覆盖圆，只要从1到j-1中，第k个点求p[k],p[j],p[i]三个点的圆，再判断k+1到j-1是否都在圆上，若都在，说明找到圆；若有不在的，则再用新的点p[k]更新圆即可。
于是，这个问题就被转化为若干个子问题来求解了。
由于三个点确定一个圆，我们的过程大致上做的是从没有确定点，到有一个确定点，再到有两个确定点，再到有三个确定点来求圆的工作。
关于正确性的证明以及复杂度的计算这里就不介绍了，可以去看完整的算法介绍：ACM 计算几何 最小圆覆盖算法
关于细节方面：PS:细节部分推荐路过的 acmer 们搞自己的模板比较好。
a.通过三个点如何求圆？
先求叉积。
若叉积为0，即三个点在同一直线，那么找到距离最远的一对点，以它们的连线为直径做圆即可；
若叉积不为0，即三个点不共线，那么就是第二个问题，如何求三角形的外接圆？
b.如何求三角形外接圆？
假设三个点(x1,y1),(x2,y2),(x3,y3)；
设过(x1,y1),(x2,y2)的直线l1方程为Ax+By=C，它的中点为(midx,midy)=((x1+x2)/2,(y1+y2)/2)，l1中垂线方程为A1x+B1y=C1；则它的中垂线方程中A1=-B=x2-x1，B1=A=y2-y1,C1=-Bmidx+Amidy=((x2^2-x1^2)+(y2^2-y1^2))/2；
同理可以知道过(x1,y1),(x3,y3)的直线的中垂线的方程。
于是这两条中垂线的交点就是圆心。
c.如何求两条直线交点？
设两条直线为A1x+B1y=C1和A2x+B2y=C2。
设一个变量det=A1B2-A2B1；
如果det=0，说明两直线平行；若不等于0，则求交点：x=(B2C1 -B1C2)/det，y=(A1C2-A2C1)/det；
d.于是木有了。。
求三角形的外接圆：推荐资料
三角形垂心， 重心，外心坐标公式
Q = 1 / (1 / (P - a^2) + 1 / (P - b^2) + 1 / （P - c^2));
λ1=Q/(P-a^2),λ2=Q/(P-b^2),λ3=Q/(P-c^2)；(注:λ1+λ2+λ3=1)
垂心:H(x,y)=λ1*A(x,y)+λ2*B(x,y)+λ3*C(x,y),
重心:G(x,y)=1/3*A(x,y)+1/3*B(x,y)+1/3*C(x,y),
外心:O(x,y)=(1-λ1)/2*A(x,y)+(1-λ2)/2*B(x,y)+(1-λ3)/2*C(x,y);
//三角形的外接圆圆心,先前已经保证了 A,B,C 三点不共线
void TriangleCircumCircle(Point A, Point B, Point C)
{
Vector a = C-B, b = A-C, c = B-A;
double da = pow(dist(a), 2), db = pow(dist(b), 2), dc = pow(dist(c), 2);
double px = (da + db + dc) / 2.0;
double Q = 1 / (1/(px-da) + 1/(px-db) + 1/(px-dc));
//t1+t2+t3 = 1
double t1 = Q/(px-da), t2 = Q/(px-db), t3 = Q/(px-dc);
Point O = A*((1-t1)/2.0) + B*((1-t2)/2.0) + C*((1-t3)/2.0) ; //外心
//Point H = A*t1 + B*t2 + C*t3; //垂心
//Point G = A*(1/3) + B*(1/3) + C*(1/3); //重心
center = O;
}参考资料：
code：
| F | Accepted | 188 KB | 0 ms | C++ (g++ 4.4.5) | 4104 B | 2013-08-20 23:05:22 |
#include<stdio.h>
#include<math.h>
#include<algorithm>
#include<iostream>
using namespace std;
const int maxn = 110;
int n;
struct Point {
double x,y;
Point() {}
Point(double _x, double _y) {
x = _x; y = _y;
}
Point operator - (const Point &b) const { //确定向量
return Point(x-b.x, y-b.y);
}
Point operator + (const Point &b) const {
return Point(x+b.x, y+b.y);
}
Point operator * (const double &k) const {
return Point(k*x, k*y);
}
double operator * (const Point &b) const { //点积
return x*b.x + y*b.y;
}
double operator ^ (const Point &b) const { //叉积
return x*b.y - y*b.x;
}
}p[maxn];
typedef Point Vector;
double dist(Point a, Point b) { //向量 ab 的长度
return sqrt((a-b)*(a-b));
}
double dist(Point a) {
return sqrt(a*a);
}
Point MidPoint(Point a, Point b) { //求 ab 的中点
return Point( (a.x+b.x) / 2.0, (a.y+b.y) / 2.0);
}
double radius; //半径
Point center; //圆心
//三角形的外接圆圆心,先前已经保证了 A,B,C 三点不共线
void TriangleCircumCircle(Point A, Point B, Point C)
{
Vector a = C-B, b = A-C, c = B-A;
double da = pow(dist(a), 2), db = pow(dist(b), 2), dc = pow(dist(c), 2);
double px = (da + db + dc) / 2.0;
double Q = 1 / (1/(px-da) + 1/(px-db) + 1/(px-dc));
//t1+t2+t3 = 1
double t1 = Q/(px-da), t2 = Q/(px-db), t3 = Q/(px-dc);
Point O = A*((1-t1)/2.0) + B*((1-t2)/2.0) + C*((1-t3)/2.0) ; //外心
//Point H = A*t1 + B*t2 + C*t3; //垂心
//Point G = A*(1/3) + B*(1/3) + C*(1/3); //重心
center = O;
}
//如果第 j 个点 pj 不在第一个点和第 i 个点为直径的圆内
//注：m = j-1
//重新以第 i 个点和第 j 个点为直径继续判断
void MiniDiscWith2Point(Point pi, Point pj, int m)
{
center = MidPoint(pi, pj);
radius = dist(pi, pj) / 2.0;
for(int k = 1; k <= m; k++) //如果p[k]不在圆内,枚举三角形【pi, pj, p[k]】的外接圆
{
if(dist(center, p[k]) <= radius) continue;
//如果第 k 个点不在当前要判断的圆内
double cross = (pi-pj)^(p[k]-pj);
if(cross != 0) //如果三点不共线
{
TriangleCircumCircle(pi, pj, p[k]);
radius = dist(center, pi);
}
else //如果三点共线
{
double d1 = dist(pi, pj); //d1 肯定不是最大的
double d2 = dist(pi, p[k]);
double d3 = dist(pj, p[k]);
if(d2 >= d3)
{
center = MidPoint(pi, p[k]);
radius = dist(pi, p[k]) / 2.0;
}
else
{
center = MidPoint(pj, p[k]);
radius = dist(pj, p[k]) / 2.0;
}
}
}
}
//当前圆不能覆盖第 i 个点 pi, 以第一个点和第 i 个点为直径的圆再次判断
//注：m = i-1
void MiniDiscWithPoint(Point pi, int m)
{
center = MidPoint(pi, p[1]); //以 第一个点和第 i 个点为直径
radius = dist(pi, p[1]) / 2.0;
for(int j = 2; j <= m; j++) //判断第 2 个点到第 i-1 个点是否包含在当前园内
{
if(dist(center, p[j]) <= radius) continue;
else MiniDiscWith2Point(pi, p[j], j-1); //如果第 j 个点不在当前圆内
}
}
int main()
{
while(scanf("%d", &n) != EOF)
{
if(n == 0) break;
double x,y;
for(int i = 1; i <= n; i++)
{
scanf("%lf%lf", &x,&y);
p[i] = Point(x,y);
}
if(n == 1)
{
printf("%.2lf %.2lf 0.00\n", p[1].x, p[1].y);
continue;
}
//第一个圆任意枚举两点都可以, 个人意见是先以凸包直径上的点为圆直径,再判断
//不过考虑到求凸包直径稍微麻烦点,而总共才 100 个点
radius = dist(p[1], p[2]) / 2.0; //先枚举第一个点和第二个点为直径
center = MidPoint(p[1], p[2]);
for(int i = 3; i <= n; i++) // 前面必定保证了第一个点到第 i-1 个点都已经被覆盖
{
if(dist(center, p[i]) <= radius) continue;
//找到第一个不能被当前圆覆盖的点
else MiniDiscWithPoint(p[i], i-1); //当前的圆不能覆盖 p[i]
}
printf("%.2lf %.2lf %.2lf\n", center.x, center.y, radius);
}
return 0;
}
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
