最小圆覆盖(随机增量法)
计算几何,三点定圆,最小圆覆盖
问题引入
题目描述
给出N个点,让你画一个最小的包含所有点的圆。
输入格式
先给出点的个数N,2<=N<=100000,再给出坐标Xi,Yi.(-10000.0<=xi,yi<=10000.0)
输出格式
输出圆的半径,及圆心的坐标,保留10位小数
在解这道题时,先要弄懂另一个知识点–三点定圆
给出三个不共线的点,求这三个点构成的三角形的外接圆
三角形的外心:三边中垂线交点,到三角形三个顶点距离相同
根据三角形外心的性质可知,三角形的外心即为三角形外接圆的圆心
//外心,三角形外接圆的圆心 Point out(Point a,Point b,Point c){ Point s; double a1=2*(b.x-a.x),b1=2*(b.y-a.y),a2=2*(c.x-a.x),b2=2*(c.y-a.y); double c1=b.x*b.x-a.x*a.x+b.y*b.y-a.y*a.y,c2=c.x*c.x-a.x*a.x+c.y*c.y-a.y*a.y; s.x=(b1*c2-b2*c1)/(a2*b1-a1*b2); s.y=(a2*c1-a1*c2)/(a2*b1-a1*b2); return s; }
求解最小圆覆盖:
若当前已考虑了前 i 个点的最小圆覆盖,则当考虑第 i+1个点时
(1)若该点在圆内,则跳过该点,继续判断下一个点
(2)若该点不在圆内,则该点应该由前 i 个点构成的最小圆内
具体步骤
1、首先要将点集随机化(防止毒瘤数据)
random_shuffle(p,p+n); //p为点集
2、设置一个半径为0的初始圆,开始遍历所有点
3、若第 i 个点不在当前确定的最小圆内,则将最小圆变成由第 i 个点为圆心,半径为0的圆,执行步骤4
4、遍历前 i-1 个点记为点 j,若 点 j 不在由点 i确定的最小圆内,则将最小圆变成以点 i 和点 j的中点为圆心,半径为两点距离的一半的圆,执行步骤5
5、遍历前 j-1 个点记为点 k,若点k不在由点 i 和点 j 确定的最小圆上,则由点i, j ,k三点,根据三点定圆,确定出最小圆。
Circle min_circular(){ random_shuffle(p,p+n); //随机化 Circle C=Circle(Point(0,0),0); for(int i=0;i<n;i++) if(sgn(Length(C.c-p[i])-C.r)>0){ //不在圆上 C.c=p[i],C.r=0; for(int j=0;j<i;j++) if(sgn(Length(C.c-p[j])-C.r)>0){ C.c=(p[i]+p[j])/2.0; C.r=Length(C.c-p[i]); for(int k=0;k<j;k++) if(sgn(Length(C.c-p[k])-C.r)>0){ C.c=out(p[i],p[j],p[k]); C.r=Length(C.c-p[i]); } } } return C; }
时间复杂度:乍一看貌似是O(n3),但实际上均摊下来,期望的复杂度是O(n)
例题:https://www.luogu.com.cn/problem/P1742
AC代码:
#include<iostream> #include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #include<vector> #include<queue> #include<stack> #include<map> #define INF 0x3f3f3f3f using namespace std; typedef long long LL; typedef unsigned long long ULL; const int maxn=1e6+6; const int mod=1e9+7; const double pi=acos(-1.0); const double inf=1e100; const double eps=1e-8; int sgn(double d){ //精度控制 if(fabs(d)<eps) return 0; if(d>0) return 1; return -1; } struct Point{ //二维平面点 double x,y; Point(double x=0, double y=0 ):x(x),y(y) { } }; typedef Point Vector ; struct Circle{ Point c; double r; Circle(Point c, double r):c(c), r(r) {}; //通过圆心角求坐标 Point point(double a){ return Point(c.x+cos(a)*r, c.y+sin(a)*r); } }; int n; Point p[maxn]; //向量加法 Vector operator + (Vector A, Vector B){ return Vector(A.x+B.x, A.y+B.y); } //向量减法 Vector operator - (Vector A, Vector B){ return Vector(A.x-B.x, A.y-B.y); } //向量与常数的乘积 Vector operator * (Vector A, double p){ return Vector(A.x*p, A.y*p); } //向量与常数的除法 Vector operator / (Vector A, double p){ return Vector(A.x/p, A.y/p); } //向量的点积(|A|*|B|*cosa) double Dot(Vector A,Vector B){ //锐角为正,钝角为负,直角为0 return A.x*B.x + A.y*B.y; } //向量的叉积(|A|*|B|*sina) double Cross(Vector A, Vector B){ return A.x*B.y - A.y*B.x; //B在A的逆时针为正,否则为负 } //取模(长度) double Length(Vector A){ return sqrt(Dot(A,A)); } //外心,三角形外接圆的圆心 Point out(Point a,Point b,Point c){ Point s; double a1=2*(b.x-a.x),b1=2*(b.y-a.y),a2=2*(c.x-a.x),b2=2*(c.y-a.y); double c1=b.x*b.x-a.x*a.x+b.y*b.y-a.y*a.y,c2=c.x*c.x-a.x*a.x+c.y*c.y-a.y*a.y; s.x=(b1*c2-b2*c1)/(a2*b1-a1*b2); s.y=(a2*c1-a1*c2)/(a2*b1-a1*b2); return s; } //求最小圆覆盖 Circle min_circular(){ random_shuffle(p,p+n); //随机化 Circle C=Circle(Point(0,0),0); for(int i=0;i<n;i++) if(sgn(Length(C.c-p[i])-C.r)>0){ //不在圆上 C.c=p[i],C.r=0; for(int j=0;j<i;j++) if(sgn(Length(C.c-p[j])-C.r)>0){ C.c=(p[i]+p[j])/2.0; C.r=Length(C.c-p[i]); for(int k=0;k<j;k++) if(sgn(Length(C.c-p[k])-C.r)>0){ C.c=out(p[i],p[j],p[k]); C.r=Length(C.c-p[i]); } } } return C; } int main(){ scanf("%d",&n); for(int i=0;i<n;i++) scanf("%lf%lf",&p[i].x,&p[i].y); Circle ans=min_circular(); printf("%.10lf\n",ans.r); printf("%.10lf %.10lf\n",ans.c.x, ans.c.y); return 0; }