计算几何,三点定圆,最小圆覆盖

问题引入

  题目描述

    给出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;
} 

 

 

 

 

 

 

 

版权声明:本文为wyy11原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://www.cnblogs.com/wyy11/p/13052410.html