最小圆覆盖

本来不想学的…于是今天就碰到一道大裸题…

例题：bzoj2823 求最小圆覆盖n个点。

伪代码如下：

把所有点随机化，设为(x[1],y[1])...(x[n],y[n])
开始把圆心设为x[1]，半径设为0
for i=2 to n
    如果i号点在当前圆内则跳过
    //那么i号点就在圆周上
    把1号点和i号点作为直径作一个圆
    for j=1 to i-1
        如果j号点在当前圆内则跳过
        考虑以j号点和i号点作为直径作一个圆
        for k=1 to j-1
            如果k号点在当前圆内则跳过
            以i,j,k三点组成的三角形的外心作为新圆心
            如果i,j,k三点共线就取ij、ik、jk连线最长的那条作为直径作为新圆

据说期望复杂度是O(n)的。

细节还是比较多的。

首先共线叉积是比较容易判的。如果叉积为0那么就共线。

那么三角形的外心怎么求呢？

如果有3个点，A(x1,y1)，B(x2,y2)，C(x3,y3)。

AB解析式：(y2-y1)x+(x1-x2)y+x2*y1-x1*y2=0

AB中点：((x1+x2)/2,(y1+y2)/2)

那么AB中垂线解析式就是(y1-y2)x+(x1-x2)y+一些常数，把AB中点的坐标带进去减一下就行了。

然后我们求一下AB和AC中垂线的交点。

两条直线

ax+by=c dx+ey=f

aex+bey=ce dbx+eby=fb

x=(fb-ce)/(ae-db)

同理

adx+bdy=cd adx+aey=af

y=(af-cd)/(ae-bd)

当然如果ae=bd就共线。

所以就解决啦。

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <string.h>
#include <math.h>
#include <vector>
#include <limits>
#include <set>
#include <map>
using namespace std;
int n;
#define SZ 2333333
typedef double db;
struct pnt
{
    db x,y;
    pnt() {}
    pnt(db a,db b) {x=a; y=b;}
}ps[SZ];
pnt operator - (pnt a,pnt b) {return pnt(a.x-b.x,a.y-b.y);}
db operator * (pnt a,pnt b) {return a.x*b.y-a.y*b.x;}
db pf(db x) {return x*x;}
db dis(pnt a,pnt b) {return sqrt(pf(a.x-b.x)+pf(a.y-b.y));}
pnt mid(pnt a,pnt b) {return pnt((a.x+b.x)/2,(a.y+b.y)/2);}
struct lne
{
    db a,b,c;
    lne() {a=b=c=0;}
    lne(db x,db y,db z) {a=x;b=y;c=z;}
};
pnt operator * (lne a,lne b)
{
    return pnt((b.c*a.b-a.c*b.b)/(a.a*b.b-a.b*b.a),
                (a.c*b.a-b.c*a.a)/(a.a*b.b-a.b*b.a));
}
lne zcx(pnt a,pnt b)
{
    db mx=(a.x+b.x)/2,my=(a.y+b.y)/2;
    db la=a.x-b.x,lb=a.y-b.y,lc=-(mx*la+my*lb);
    return lne(la,lb,lc);
}
pnt wx(pnt a,pnt b,pnt c) {return zcx(a,b)*zcx(a,c);}
db eps=1e-6;
pnt np(pnt a,pnt b,pnt c)
{
    if(fabs((b-a)*(c-a))<eps)
    {
        double abl=dis(a,b),acl=dis(a,c),bcl=dis(b,c),ans=-1;
        pnt ap;
        if(abl>ans) ans=abl, ap=mid(a,b);
        if(acl>ans) ans=acl, ap=mid(a,c);
        if(bcl>ans) ans=bcl, ap=mid(b,c);
        return ap;
    }
    return wx(a,b,c);
}
db gr(pnt o,pnt a,pnt b,pnt c) {return max(max(dis(o,a),dis(o,b)),dis(o,c));}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%lf%lf",&ps[i].x,&ps[i].y);
    random_shuffle(ps+1,ps+1+n);
    pnt yx=ps[1]; double r=0;
    for(int i=2;i<=n;i++)
    {
        if(dis(yx,ps[i])<r+eps) continue;
        yx=mid(ps[1],ps[i]);
        r=dis(yx,ps[i]);
        for(int j=1;j<i;j++)
        {
            if(dis(yx,ps[j])<r+eps) continue;
            yx=mid(ps[i],ps[j]);
            r=dis(yx,ps[i]);
            for(int k=1;k<j;k++)
            {
                if(dis(yx,ps[k])<r+eps) continue;
                yx=np(ps[i],ps[j],ps[k]);
                r=gr(yx,ps[i],ps[j],ps[k]);
            }
        }
    }
    printf("%.2lf %.2lf %.2lf\n",yx.x,yx.y,r);
}