最小圆覆盖问题算法步骤与证明+代码模板

客官°小女子只卖身不卖艺 2022-11-07 13:10 192阅读 0赞

### 最小圆覆盖定义 ###

最小圆覆盖问题是什么呢？就是指在二维平面上有一堆点，然后我们要求一个最小半径的圆能够将所有点全部都包住，这就是最小圆覆盖问题。

--------------------

### 最小覆盖圆的性质 ###

#### 性质1：最小覆盖圆是唯一的 ####

证明：我们假设有两个圆O1，O2，他们半径都是r，都是最小覆盖圆，那么所有的点一定在两圆的交集部分。那我们以两圆交集部分的弦长为直径，做一个新圆，该圆依旧包含所有点，而且他的直径是原来圆的弦长，一定是小于原来圆的直径的，因此和原来的圆是最小覆盖圆相矛盾，所以最小覆盖圆是唯一的。  
![eg][]

#### 性质2：若圆O1是点集S的最小覆盖圆，则再新加一点p，他在集合S的外部，则新点集的最小覆盖圆一定过点p，即p一定在最小覆盖圆上。 ####

同理我们明显能看到，若新的最小覆盖圆将S点集包住也将p包住，过点p的圆半径一定最小。  
我们将这个大圆朝原先的小圆缩小，在缩小的过程中一定会过点p，p就是这个临界状态，因为在不断缩小，所以p一定在既包含S又包含p的最小覆盖圆上。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NzM1ODUx_size_16_color_FFFFFF_t_70]

#### 性质3：最终的最小覆盖圆一定只有两种情况：1是圆上至少有三个点，由三点限制一个圆（三点共圆）；2是圆上只有两个点，则该圆一定是以该两点的连线为直径的。 ####

如果最小覆盖圆上只有一个点的话，我们肯定可以通过平移加缩小找到一个更小的圆包住所有点，与最小不符。  
![eg][eg 1]

--------------------

### 最小圆覆盖算法步骤 ###

随机增量法：就是从所有点里随机选一个加入自己的集合中，然后维护这个集合，然后再随机选一个点加入集合，再维护… 直到所有点都加入这个集合，那么最终的这个集合就是满足我要求的情况了。  
然后我们找一个圆就是看能不能找到这个圆上的三个点。

首先我们**先将所有点随机化**。  
**三层循环**，**第一层**从2到n循环i，我们将圆初始化为以第一个点为圆心，半径为0的圆。然后我们要循环判断i点在圆内还是在圆外，若在圆内（在圆上也认为是圆内），那Ci和Ci-1是相同的，我们直接循环i+1就行了；如果i在圆外，那么Ci一定过点i（性质2），然后我们就进第二层循环。  
**第二层**从1到i-1循环j，我们将圆先初始化为以点i为圆心，0为半径的圆，然后不断去加点j，若j在圆内，就继续循环j+1；若j在圆外，就再进第三层循环。  
**第三层**从1到j-1循环k，将圆先初始化为以i和j为直径的圆，然后再循环判断k，若在圆内，就继续循环下一个k+1；若k在圆外，就将圆更新为i、j、k三点确定的圆，直到最后求出来的圆就是最小覆盖圆。

--------------------

### 证明 ###

如果只看步骤的话挺简单，可跳过证明，看后面的代码模板即可。

这个证明还是比较绕的，要用到前面的三条性质。我们的目的是为了找到三个点来确定这个圆。首先最外层循环1-n，圆最初是以第一个点为圆心，半径为0的圆。然后我们每次循环到i时，我们已经知道了前i-1个点的最小覆盖圆了，然后我们就判断点i，若i在圆内，那就不用管i了，直接循环下一个点i+1，因为Ci和Ci-1是同一个圆；若i在圆外，那么Ci一定过i点。  
现在我们知道了一个信息，**1-i的最小覆盖圆过i点**，但我们无法确定这个圆（因为我们只知道这个圆上的一个点），那么我们就再来一层循环：从1到i-1循环j，根据我们已知的信息（Ci必过i点），我们现在要求的就是**过i点的**、**包含1-j个点**的最小覆盖圆。所以我们将圆初始化为以i为圆心，以0为半径的圆，不断地去加j，判断j在圆内还是圆外，圆内就不用管；在圆外的话就说明我们求得的圆一定过点j（性质2依然适用）。  
那么现在我们就又知道了一个信息，**Ci必过i点和j点**，我们就利用这个信息，既然Ci必过i点和j点，那么我们就**只看过i和j的圆**。所以我们将圆初始化为以i、j为直径的圆（这个圆是包含i、j的最小覆盖圆）。然后我们再来一层循环：从1到j-1循环k，根据我们得到的信息，我们要求过i和j点且包含1-k个点的最小覆盖圆。然后我们循环判断k，若k在圆内，就下一个；若k在圆外，就说明新圆一定过k（性质2），我们就得到了一个过i、j、k三点的圆，因此也就确定了一个唯一的圆。  
当k循环到j-1时，就说明我们找到了一个包含1到j-1，并且i和j都在圆边上的最小覆盖圆；然后当j循环到i-1时，就说明我们找到了一个包含1到i-1且i点在圆边是的最小覆盖圆；然后当i循环到n时，就说明我们找到了一个包含1到i-1且i点在圆边上的最小覆盖圆，即包含前n个点的最小覆盖圆。

然后还有一个问题就是**已知三点，如何求圆**。就是求这三点组成的三角形的外接圆。求出三角形两个边的中垂线，再求出中垂线交点就是圆心，圆心到任意一点的距离就是半径。求中垂线的话得先求出两边的中点，再求出两个边绕端点旋转90°的向量，直到一点和这条直线的向量，就可以写出点向式了。关于向量旋转、求交点函数都在[基础知识][Link 1]中有。

然后来看看他的**时间复杂度**，虽然看着是三重循环，但是复杂度并不 O ( n 3 ) O(n^3) O(n3)，因为每次循环下一层前都有一个判断，平时该点是不是在圆上，而三点确定一个圆，所以我们只有 3 n \\frac 3 n n3的机会循环下一层，循环j层和k层的复杂度都是 3 n ∗ O ( n ) \\frac 3 n \* O(n) n3∗O(n)，最终的复杂度就是 O ( n ) ∗ 3 n ∗ O ( n ) ∗ 3 n ∗ O ( n ) O(n)\*\\frac 3 n \* O(n)\*\\frac 3 n \* O(n) O(n)∗n3∗O(n)∗n3∗O(n)，近似于 O ( n ) O(n) O(n)，这个复杂度还是可以接受的。因此，每次的随机化的过程是不可省略的，否则有可能会被出题人数据卡到 O ( n 3 ) O(n^3) O(n3)的复杂度去。

--------------------

### 代码模板 ###

#include<iostream>
    #include<cmath>
    #include<algorithm>
    using namespace std;
    const int N = 100010;
    const double eps = 1e-8,PI = acos(-1);	//会用到的常数
    #define x first //方便用pair
    #define y second
    typedef pair<double,double> PDD;
    struct Circle{ 		//存圆
    	PDD o;		//圆心o
    	double r;	//半径r
    }c;
    
    PDD p[N];	//p来存点
    int n;
    
    PDD operator+(PDD a,PDD b)	//重在运算符+，方便pair运算
    { 
    	return { a.x+b.x, a.y+b.y};
    }
    
    PDD operator-(PDD a,PDD b)
    { 
    	return { a.x-b.x, a.y-b.y};
    }
    
    double operator*(PDD a,PDD b)	//叉乘
    { 
    	return a.x*b.y - a.y*b.x;	//x1*y2 - x2*y1
    }
    
    PDD operator*(PDD a,double b)	//数乘
    { 
    	return { a.x*b, a.y*b};
    }
    
    PDD operator/(PDD a,double b)
    { 
    	return { a.x/b, a.y/b};
    }
    
    int judge(double a,double b)	//两个double数比大小，相等返回0，小于返回-1，大于返回1
    { 
    	if(fabs(a-b) < eps)
    		return 0;
    	if(a < b)
    		return -1;
    	return 1;
    }
    
    PDD rotate(PDD a,double b)	//向量a旋转b°后的向量，基础知识中有证明
    { 
    	return { a.x*cos(b)+a.y*sin(b), -a.x*sin(b)+a.y*cos(b)};
    }
    
    double lens(PDD a,PDD b)	//求两点距离
    { 
    	double dx = a.x - b.x;
    	double dy = a.y - b.y;
    	return sqrt(dx*dx + dy*dy);	//根号下x方加y方
    }
    
    PDD intersection(PDD p,PDD v,PDD q,PDD w)	//求两直线交点，详见基础知识
    { 
    	PDD u = p - q;
    	double t = w*u / (v*w);
    	return p + v*t;
    }
    
    pair<PDD,PDD> bisector(PDD a,PDD b)	//求两点中垂线
    { 
    	PDD p = (a+b) / 2;	//先求出a和b连线的中点，就是横纵坐标和的一半
    	PDD v = rotate(b-a,PI/2);	//然后求ab向量旋转90°的向量
    	return { p,v};		//就有了直线上一点p和一个向量v，根据点向式可以写出一条直线方程
    }
    
    Circle circle(PDD a,PDD b,PDD c)	//已知三点，求外接圆
    { 
    	auto n = bisector(a,b),m = bisector(a,c);	//先求出两个边的中垂线
    	PDD o = intersection(n.x,n.y,m.x,m.y);	//然后求两中垂线交点就是圆心o
    	double r = lens(o,a);	//求出oa长度就是半径
    	return { o,r};		//返回圆心和半径
    }
    
    int main()
    { 
    	cin >> n;
    	for(int i = 0;i < n;i++)
    		cin >> p[i].x >> p[i].y;	//输入个点
    	
    	random_shuffle(p,p+n);		//随机化p
    	c = { p[0],0};			//初始化圆为p[0]
    	for(int i = 1;i < n;i++)	//第一次循环
    	{ 
    		if(judge(c.r,lens(c.o,p[i])) == -1)	//如果p[i]在圆外
    		{ 
    			c = { p[i],0};			//初始化圆为p[i]
    			for(int j = 0;j < i;j++)	//第二层循环
    			{ 
    				if(judge(c.r,lens(c.o,p[j])) == -1)	//如果p[j]在圆外
    				{ 
    					c = { (p[i]+p[j])/2,lens(p[i],p[j])/2};	//初始化圆为p[i]和p[j]为直径的圆
    					for(int k = 0;k < j;k++)	//第三层循环
    					{ 
    						if(judge(c.r,lens(c.o,p[k])) == -1)	//如果p[k]在圆外
    							c = circle(p[i],p[j],p[k]);		//i、j、k三点求一个唯一的圆
    					}
    				}
    			}
    		}
    	}
    	printf("%.10lf\n",c.r);		//输出圆半径
    	printf("%.10lf %.10lf\n",c.o.x,c.o.y);	//输出圆心坐标
    
    	return 0;
    }

--------------------

模板题：[P1742：最小圆覆盖][P1742]、  
代码中几何知识[详见基础知识][Link 1]

除了证明比较绕，代码步骤还是比较好写的。

[eg]: /images/20221023/b9ae41577d5244479a76864cd3cab067.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NzM1ODUx_size_16_color_FFFFFF_t_70]: /images/20221023/99352b8dd998491c9fb76d675d929ff4.png
[eg 1]: /images/20221023/e9bbfb9dc4ef4f1d9a9562c81d19bf20.png
[Link 1]: https://blog.csdn.net/qq_45735851/article/details/114382634
[P1742]: https://www.luogu.com.cn/problem/P1742