凸包最小外切圆(最小覆盖圆)

定义: 给出平面上的一些点,求覆盖这些点的最小圆。

具体问题可见hdu 2215

方法一:O(n3)

参考

解决:
1.求这些离散点的凸包;
2.枚举凸包上任意3点形成的三角形,求能覆盖这个三角形的最小圆(不是外切圆),找出所有最小圆中半径最大的一个。

细节:
1.若枚举的三个点构成钝角三角形,则最小圆半径为最长边的一半。否则,半径 r=a*b*c/(4*S) ,其中S为三角形面积,用叉积很好求。
2.钝角三角形判断:若最长边为c,则a*a+b*b<c*c(类比直角三角)

时间复杂度O(n3)

//n个点,每次输入一个点的x,y值,每个点有0.5大小的半径,求最小覆盖圆
#include<bits/stdc++.h>
#define ll long long
#define lf double
#define IOS std::ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define Rep(i,l,r) for(int i=(l);i<=(r);i++)
using namespace std;

int n;

struct Point
{
    double x,y,z;
};
typedef Point vector_t;

double Dot(vector_t A, vector_t B){
    return A.x*B.x + A.y*B.y;
}
double Length(vector_t A){
    return sqrt(Dot(A, A));
}

double Cross(const vector_t& x, const vector_t& y)
{
	double z = x.x * y.y - x.y * y.x;
	return z;
}

Point operator-(const Point&a,const Point&b)
{
    Point c;
    c.x=a.x-b.x;
    c.y=a.y-b.y;
    c.z=0;
    return c;
}

Point p[500];
int cmp(const Point&a,const Point&b){
    double d=Cross(a-p[0],b-p[0]);
    if(d>0)
        return 1;
    if(d==0 && Length(a-p[0])<=Length(b-p[0]))
        return 1;
    return 0;
}
int build(int n,Point *p,Point *s)
{
	//*p为散点,*s为凸包的点
    sort(p,p+n,cmp);
    int m = 0;
    for(int i = 0;i < n;i++)
    {
        while(m>1&&Cross((s[m-1]-s[m-2]),(p[i]-s[m-2]))<=0)m--;
        s[m++] = p[i];
    }
    int lim = m;
    for(int i = n-2;i >= 0;i--)
    {
        while(m>lim&&Cross((s[m-1]-s[m-2]),(p[i]-s[m-2]))<=0)m--;
        s[m++] = p[i];
    }
    if(m!=1)m--;
    return m;
}

Point s[500];

int main()
{
    while(cin>>n)
    {
        if(n==0)break;
        Rep(i,0,n-1){
            double x,y;
            cin>>x>>y;
            p[i].x=x;p[i].y=y;
        }
        if(n==1){
            printf("0.50\n");
            continue;
        }
        
        int m=build(n,p,s);
        double R=0;
        if(m<=2){
            double r=(Length(p[1]-p[0])+1)/2;
            R=max(R,r);
        }
        else{
            for(int i=0;i<m;i++){
                for(int j=i+1;j<m;j++){
                    for(int k=j+1;k<m;k++){
                        double len1=Length(s[j]-s[i]),len2=Length(s[k]-s[i]),
                        	len3=Length(s[k]-s[j]);
                        if(len1>len2)swap(len1,len2);
                        if(len2>len3)swap(len2,len3);
                        double r=1e9;
                        if((len1*len1)+(len2*len2)<(len3*len3)){
                            r=len3/2;
                        }
                        else{
                            double S=abs(Cross(s[j]-s[i],s[k]-s[i]));
                            r=len1*len2*len3/(2*S);
                        }
                        R=max(R,r+0.5);
                    }
                }
            }
        }
        printf("%.2lf\n",R);
    }
    return 0;
}

方法二:O(n)

参考:最小圆覆盖算法

更推荐的方法。均摊复杂度O(n),因为有随机化,最差复杂度也是O(n3)。

例题:P1742 最小圆覆盖

#include<bits/stdc++.h>
#define ll long long
#define lf double
#define IOS std::ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define Rep(i,l,r) for(int i=(l);i<=(r);i++)
using namespace std;

struct Point{double x,y;};
struct Circle{Point o;double r;};
typedef Point vector_t;

double Dot(vector_t A, vector_t B){return A.x*B.x + A.y*B.y;}
double Length(vector_t A){return sqrt(Dot(A, A));}

double Cross(const vector_t& x, const vector_t& y) {
	return x.x * y.y - x.y * y.x;
}  //叉积
vector_t operator*(const vector_t& v, double t) {
	return (vector_t){v.x * t, v.y * t};
}
vector_t operator+(const vector_t& a, const vector_t& b) {return (vector_t){a.x + b.x, a.y + b.y};}
Point operator-(const Point& a, const Point& b) {return (Point){a.x - b.x, a.y - b.y};}

Circle GetCircumcircleTriangle(const Point& p1, const Point& p2, const Point&p3)
{
    double a,b,c,d,e,f;
    a=p1.x-p2.x; b=p1.y-p2.y; c=(Dot(p2,p2)-Dot(p1,p1))/2;
    d=p1.x-p3.x; e=p1.y-p3.y; f=(Dot(p3,p3)-Dot(p1,p1))/2;
    double x=(f*b-c*e) / (a*e-b*d);
    double y=(f*a-c*d) / (b*d-e*a);
    Point o={x,y};double r=Length(p1-o);
    return (Circle){o,r};
}

Circle GetMinCoverageCircle(Point* p,int n)
{
    double r = 0;
    Point o = {0,0};
    random_shuffle(p,p+n);
    for(int i=0;i<n;i++){
        if(Length(p[i]-o)<=r)continue;
        o=p[i];r=0;
        for(int j=0;j<i;j++){
            if(Length(p[j]-o)<=r)continue;
            o.x=(p[i].x+p[j].x)/2;
            o.y=(p[i].y+p[j].y)/2;
            r=Length(p[i]-o);
            for(int k=0;k<j;k++){
                if(Length(p[k]-o)<=r)continue;
                Circle tmp = GetCircumcircleTriangle(p[i],p[j],p[k]);
                o=tmp.o;r=tmp.r;
            }
        }
    }
    return (Circle){o,r};
}
int n;
Point p[100005];
int main()
{
    cin>>n;
    Rep(i,0,n-1){
        double x,y;
        cin>>x>>y;
        p[i].x=x;p[i].y=y;
    }
    
    Circle ans=GetMinCoverageCircle(p,n);
    printf("%.10lf\n%.10lf %.10lf\n",ans.r,ans.o.x,ans.o.y);
    return 0;
}
  • 1
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值