最小圆覆盖

我们有一些点,需要用一个圆将所有点覆盖,要求圆的半径最小,这样的问题我们称为最小圆覆盖
如何求最小圆呢?
首先我们要先了解2最小圆覆盖的两条性质:
1.最小覆盖圆是唯一的
2.若点P不在S的最小覆盖圆内部,则P一定在 P ⋂ S P \bigcap S PS
的最小覆盖圆边上
通过这两个性质,我们不难得出,如果点i不在 1 - i-1 的最小覆盖圆上,那么i一定在1 - i的最小覆盖圆的边上,这样,我们就得到了构成圆的一个点,对于圆来说,我们至少需要找到三个点才能构成一个圆,
所以我们要在1 - i-1 里面找到第二个点j,对于点j,在找的时候需要满足两个条件,不在1 - j - 1内部,和i必须在圆边上,最后我们再在1 - j-1里面找一个点k,条件同上,这样我们就确定了一个圆。

#include <cstdio>
#include <iostream>
#include <cmath>
#include <utility>
#include <climits>
#include <cstring>
#include <ctime>
#include <algorithm>
#define eps 1e-12
#define N_MAX 100005
#define PI acos(-1)
using std::pair;
int n;
int sign(double x){ //符号
    if(fabs(x) < eps)return 0;
    if(x < 0)return -1;
    return 1;
}
struct Point{
    double x,y;
    Point(double x = 0,double y = 0):x(x),y(y){}
    Point operator+(Point a){
        return Point(x + a.x,y + a.y);
    }
    Point operator-(Point a){
        return Point(x - a.x,y - a.y);
    }
    Point operator*(double t){
        return Point(x * t,y * t);
    }
    Point operator/(double t){
        return Point(x / t,y / t);
    }
    bool operator==(Point a){
        return x == a.x && y == a.y;
    }
    double lenth(){ //向量的模
        return sqrt(x * x + y * y);
    }
}point[N_MAX];
typedef Point Vector;
struct Line{
    Point a,b;
    Line(){}
    Line(Point a,Point b){
        this->a = a;
        this->b = b;
    }
};
struct Circle{
    Point o;
    double r;
    Circle(){}
    Circle(Point o,double r){
        this->o = o;
        this->r = r;
    }
}c;
double cross(Vector A,Vector B){  // 向量外积
    return A.x * B.y - B.x * A.y;
}
double dot(Vector A,Vector B){ //向量内积
        return B.x * A.x + B.y * A.y;
}
bool isclock(Vector v1,Vector v2){ // 判断v2是否在v1的顺时针方向
    if(cross(v1,v2) < 0)return true;
    return false;
}
double GetAngle(Vector v1,Vector v2){ // 两向量夹角
    return acos(dot(v1,v2) / v1.lenth() / v2.lenth());
}
Vector rotate(Vector v,double angle){ //向量v偏转angle后的位置
    return Vector(v.x * cos(angle) + v.y * sin(angle),-v.x * sin(angle) + v.y * cos(angle));
}
Point GetPointForLine(Line X,Line Y){ // 算直线交点
    Vector v = X.b;
    Vector w = Y.b;
    Vector u = X.a - Y.a;
    double base = (cross(w,u)) / (cross(v,w));
    return X.a + v * base;
}
double DistanceForLine(Point P,Line X){ // 算点到直线的距离
    Vector vp = X.b - X.a,vq = P - X.a;
    return fabs(cross(vp,vq))/vp.lenth();
}
double DistanceForSegment(Point P,Line X){ //点到线段的距离
    if(X.a == X.b)return (P - X.a).lenth();
    Vector v1 = X.b - X.a,v2 = P - X.a,v3 = P - X.b;
    if(sign(dot(v1,v2)) < 0) return v2.lenth();
    if(sign(dot(v1,v3)) > 0) return v3.lenth();
    return DistanceForLine(P,X);
}
Point GetProjectionInLine(Point P,Line X){ //点到直线上的投影
    Vector v = X.b - X.a;
    return X.a + v * (dot(v,P - X.a) / dot(v,v));
}
bool PointOnSegment(Point P,Line X){ // 点是否在线段上
    return sign(cross(P - X.a,P - X.b)) == 0 && sign(dot(P - X.a,P - X.b)) <= 0;
}
bool SegmentIntersection(Line X,Line Y){ //两线段是否相交
    double c1 = cross(X.b - X.a,Y.a - Y.b),c2 = cross(X.b - X.a,Y.b - X.a);
    double c3 = cross(Y.b - Y.a,X.b - Y.a),c4 = cross(Y.b - Y.a,X.a - Y.a);
    return sign(c1) * sign(c2) <= 0 && sign(c3) * sign(c4) <= 0;
}
pair<Point,Point> PointWithCircleAndLine(Point o,Line X,double r){ // 求圆与直线的交点
    Point o2 = o;
    o2.x += X.a.y - X.b.y;
    o2.y += X.b.x - X.a.x;
    o2 = GetPointForLine(Line(o,o2),X);
    double base = sqrt(r * r - (o2 - o).lenth() * (o2 - o).lenth());
    Vector e = (X.b - X.a) / (X.b - X.a).lenth();
    return std::make_pair(o2 - e * base,o2 + e * base);
}
Line GetPerpendicularBisector(Point P,Point Q){
    return Line((P + Q) / 2,rotate(Q - P,PI / 2));
}
Circle GetCircle(Point A,Point B,Point C){
    Line AB = GetPerpendicularBisector(A,B),AC = GetPerpendicularBisector(A,C);
    Point o = GetPointForLine(AB,AC);
    return Circle(o,(o - A).lenth());
}
int main(int argc,char *argv[]){
    srand(time(NULL));
    scanf("%d",&n);
    for(int i = 0;i < n;i++)scanf("%lf %lf",&point[i].x,&point[i].y);
    std::random_shuffle(point,point + n);
    c = Circle(point[0],0);
    for(int i = 1;i < n;i++){
        if(sign(c.r - (c.o - point[i]).lenth()) < 0){
            c = Circle(point[i],0);
            for(int j = 0;j < i;j++){
                if(sign(c.r - (c.o - point[j]).lenth()) < 0){
                    c = Circle((point[i] + point[j]) / 2,(point[i] - point[j]).lenth() / 2);
                    for(int k = 0;k < j;k++){
                        if(sign(c.r - (c.o - point[k]).lenth()) < 0){
                            c = GetCircle(point[i],point[j],point[k]);
                        }
                    }
                }
            }
        }
    }
    printf("%.10lf\n",c.r);
    printf("%.10lf %.10lf\n",c.o.x,c.o.y);
    return 0;
}
  • 1
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值