最小圆覆盖（经典算法【三点定圆】）

问题描述

给定n个点，用一个最小的圆把这些点全部覆盖，求这个圆的圆心半径

 

于是，这个问题就被转化为若干个子问题来求解了

由于三个点确定一个圆，我们的过程大致上做的是从没有确定点，到有一个确定点，再到有两个确定点，再到有三个确定点来求圆的工作

时间复杂度：O(N)

空间复杂度：O(N)

 

小细节

Q1.

过三点如何求圆？

A1.

先求叉积

若叉积为0，即三个点在同一直线，那么找到距离最远的一对点，以它们的连线为直径做圆即可；

若叉积不为0，即三个点不共线，那么就求三角形的外接圆

Q2.

如何求三角形外接圆？

node jiao(node p,node v,node q,node w)
//p+tv
//q+tw
{
    node u=p-q;
    double t=Cross(w,u)/Cross(v,w);
    return p+v*t;
}

 

给出板子（几何算法）：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<ctime>
using namespace std;
typedef long long lt;
#define eps 1e-6
#define sqr(x) ((x)*(x))
 
const int maxn=1000010;
int n;
struct point{
	double x,y;
}p[maxn],O;
double R;//半径 
 
double getd(point a,point b){ //求直径 
	return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}
point getO(point p1,point p2,point p3) { //求圆心 
    point res;
    double a=p2.x-p1.x;
	double b=p2.y-p1.y;
	double c=p3.x-p2.x;
	double d=p3.y-p2.y;
    double e=sqr(p2.x)+sqr(p2.y)-sqr(p1.x)-sqr(p1.y);
    double f=sqr(p3.x)+sqr(p3.y)-sqr(p2.x)-sqr(p2.y);
    res.x=(f*b-e*d)/(c*b-a*d)/2.0; 
    res.y=(a*f-e*c)/(a*d-b*c)/2.0;
    return res;
}
void mincir() {
    O=p[1]; 
	R=0;
    for(int i=1;i<=n;++i){
	    if(getd(p[i],O)-R>eps) { //不在圆内 
	        O=p[i]; 
			R=0;
	        for(int j=1;j<i;j++) { 
		        if(getd(p[j],O)-R>eps) {//不在圆内 
		            O=(point){(p[i].x+p[j].x)/2.0,(p[i].y+p[j].y)/2.0};
		            R=getd(p[i],p[j])/2.0;
		            for(int k=1;k<j;++k)
		            if(getd(p[k],O)-R>eps) {//不在圆内 
		                O=getO(p[i],p[j],p[k]);  //外接圆 
		                R=getd(p[i],O);
		            }
		        }
		    } 
	    }
	} 
} 
 
int main()
{
    cin>>n;
    for(int i=1;i<=n;++i)
    scanf("%lf%lf",&p[i].x,&p[i].y);
    random_shuffle(p+1,p+1+n);// random_shuffle()随机打乱函数 首指针 尾指针 
    mincir();
    printf("%.3f",R);
    return 0;
}

给出板子（模拟退火法）：

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdio>
#include<cmath>
using namespace std;
const double eps=1e-8;
struct POINT{
    double x,y,z;
}p[510];
POINT op;//最小圆的圆心
int n;
inline double dist(POINT &a,POINT &b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
void solve(){
    double ret,delta=100.0;
    double maxDis,tempDis;
    int i,id;
    while(delta>eps){
        id=0;
        maxDis=dist(op,p[id]);
        for(i=1;i<n;i++){
            tempDis=dist(op,p[i]);
            if(tempDis>maxDis){
                maxDis=tempDis;
                id=i;
            }
        }
        ret=maxDis;
        op.x+=(p[id].x-op.x)/maxDis*delta;
        op.y+=(p[id].y-op.y)/maxDis*delta;
        delta*=0.98;
    }
    printf("%.3f\n",ret);
}
int main(){
	scanf("%d",&n);
    op.x=op.y=0;
    for(int i=0;i<n;i++){
        scanf("%lf%lf",&p[i].x,&p[i].y);
        op.x+=p[i].x;
        op.y+=p[i].y;
    }
    op.x/=n;
    op.y/=n;
    solve();
    
    return 0;
}

 

这里给出一个模板题：P1742

　最小圆覆盖问题：给定n个点的平面坐标，求一个半径最小的圆，把n个点全部包围，部分点在圆上。（两种算法：几何算法和模拟退火算法）

　　几何算法：（1）加第1个点P1。C1的圆心就是P1，半径为0。

　　　　　　　（2）加第二个点P2。新的C2的圆心是线段P1P2的中心，半径为两点距离的一半。这一步操作是两点定圆。

　　　　　　　（3）加第三个点P3。若P3在圆内或圆上，忽略；若不在，则以P3为圆心，重复（1）和（2），若还是不行则用三点定圆。

　　　　　　　（4）加第四个点P4。若P4在圆内或圆上，忽略；若不在，则以P4为圆心，从前三个点中选择一个点重复（1）和（2）即两点定圆，若还是不行则选三个点进行三点定圆(一定有)。

　　　　　　　（5）继续加入新的点。

　　复杂度分析：3层for循环，貌似是O(n3)，但是当点的分布是随机的时候，可以通过概论计算得到实际复杂度接近O(n)，代码中使用random_shuffle()函数实现。

#include <bits/stdc++.h>
using namespace std;
#define eps 1e-8
const int maxn = 100000+5;
int sgn(double x)
{
    if (fabs(x)<eps) return 0;
    else return x<0? -1:1;
}
struct Point
{
    double x,y;
};
double Distance(Point A, Point B){return hypot(A.x-B.x,A.y-B.y);}
//求三角形abc的外接圆圆心
Point circle_center(const Point a, const Point b, const Point c)
{
    Point center;
    double a1 = b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2;
    double a2 = c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2;
    double d=a1*b2-a2*b1;
    center.x=a.x+(c1*b2-c2*b1)/d;
    center.y=a.y+(a1*c2-a2*c1)/d;
    return center;
}
//求最小圆覆盖，返回圆心c和半径r:
void min_cover_circle(Point *p, int n,Point &c, double &r)
{
    random_shuffle(p,p+n);             //打乱所有点
    c = p[0];r = 0;                    //第一个点
    for (int i = 1; i < n; ++i)        //剩下所有点
    {
        if (sgn(Distance(p[i],c)-r)>0) //pi在圆外部
        {                              
            c=p[i];r=0;                //将圆心设为pi半径为0
            for (int j = 0; j < i; ++j) //重新检查前面的点
            {
                if (sgn(Distance(p[j],c)-r)>0)//两点定圆
                {
                    c.x=(p[i].x+p[j].x)/2;
                    c.y=(p[i].y+p[j].y)/2;
                    r=Distance(p[j],c);
                    for (int k = 0; k < j; ++k)
                    {
                        if (sgn(Distance(p[k],c)-r)>0)
                        {
                            c=circle_center(p[i],p[j],p[k]);
                            r=Distance(p[i],c);
                        }
                    }
                }
            }
        }
    }
}
int main(int argc, char const *argv[])
{
    int n;
    Point p[maxn];
    Point c;double r;
    while(~scanf("%d",&n)&&n)
    {
        for (int i = 0; i < n; ++i)
        {
            scanf("%lf%lf",&p[i].x,&p[i].y);
        }
        min_cover_circle(p,n,c,r);
        printf("%.10lf\n%.10lf %.10lf\n",r,c.x,c.y);
    }
    return 0;
}

模拟退火算法：

void min_cover_circle(Point *p, int n,Point &c, double &r)
{
    double T = 100.0;                 //初始温度
    double delta = 0.98;             //降温系数
    c = p[0];                        
    int pos;
    while(T>eps)                     //eps终止温度
    {
        pos = 0; r = 0;
        for (int i = 0; i <= n-1; ++i) //初始:p[0]是圆心，半径是0
        {
            if (Distance(c,p[i])>r)   
            {
                r = Distance(c,p[i]); //距圆心最远的点看到
                pos = i;
            }
            c.x += (p[pos].x-c.x)/r * T; //逼近最后的解
            c.y += (p[pos].y-c.y)/r * T;
            T *= delta;
        }
    }
}