最小圆覆盖 学习笔记【计算几何】

点击量：432

前言

最小圆覆盖和最小矩形覆盖其实是有一定区别的。在矩形中有两组平行直线，在圆中则没有。

随机增量法

因为复杂度是在期望意义下计算的，所以无论题目给出什么样的数据，我们都要随机一下，这就是“随机”的含义。

我们考虑递推的过程。

假设我们现在已经求出来前 $ i-1$ 个点的最小覆盖圆。如果第 $ i$ 个点在这个圆内，那么这个点不产生影响。如果第 $ i$ 个点在圆外，我们考虑重构这个最小覆盖圆。

令 $ i$ 号点 $ P_i$ 为圆心，$ r=0$ 为半径，枚举 $ j:1\text{ to }i-1$，重新进行递推。

假设我们现在已经求出来 $ P_i$ 和第 $ 1\sim j-1$ 个点的最小覆盖圆。如果第 $ j$ 个点在这个圆内，那么这个点不产生影响。如果第 $ j$ 个点在圆外，那么同样重构这个最小覆盖圆。但是此时已经钦定了 $ i$ 号点在圆上，所以设 $ i$ 号点和 $ j$ 号点的中点为 $ M$。

令 $ M$ 点为圆心，$ r=\frac{\left|P_iP_j\right|}2$ 为半径，枚举 $ k:1\text{ to }j-1$，进一步递推。

现在我们已经知道了 $ P_i,P_j$ 一定在圆上。假设我们已经求出来 $ P_i,P_j$ 和 $ 1\sim k-1$ 号点的最小覆盖圆。如果 $ k$ 不在圆内，直接更新答案（此时是过 $ P_i$ 和 $ 1\sim j$ 号的点的最小覆盖圆）为 $ P_i,P_j,P_k$ 的外接圆。

过 $ P_i$ 和 $ 1\sim j$ 号点的最小覆盖圆找到了，我们下一步就可以做 $ 1\sim j+1$ 了，返回第 5 段（是自然段），进一步进行。

综上所述本题一共用了 3 层迭代，最终通过递推的方式找出了最小覆盖圆。

复杂度分析

但是这样的三重循环看样子是 $ O\left(n^3\right)$ 的啊。

我们现在回到一开始的随机操作，它的意义是让数据尽可能达到期望意义。

对于过 $ P_i,P_j$ 的圆的枚举，它全面地涉及到了 $ 1\sim j-1$ 中的每个点，至少循环了一遍，这一层复杂度是 $ O(j)$。

对于过 $ P_i$ 的圆的枚举，它要枚举 $ 1\sim i-1$ 中的那些点，考虑到三点定圆，覆盖到前 $ i$ 个点的最小覆盖圆也是由三个点定下来的，因此在前 $ i-1$ 号点中，期望有两个点做出贡献。因此能迭代进下一层的只有 $ O\left(\sum_{j=1}^i\frac 2j\right)$ 个。

同理，对于全局最小圆覆盖，在前 $ i$ 号点中，期望有三个点能做出贡献。能迭代进下一层的只有 $ O\left(\sum_{i=1}^n\frac 3i\right)$ 个。

总复杂度为 $ O\left(\sum_{i=1}^n\frac 3i\cdot \sum_{j=1}^i\frac 2j\cdot j\right)=O(n)$。

代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define db double
#define eps 1e-12
struct pnt
{
    db x,y;
    pnt(db x,db y)
    {
        this->x=x;
        this->y=y;
    }
    pnt(){}
    friend pnt operator +(pnt a,pnt b)
    {return pnt(a.x+b.x,a.y+b.y);}
    friend pnt operator -(pnt a,pnt b)
    {return pnt(a.x-b.x,a.y-b.y);}
    friend db operator *(pnt a,pnt b)
    {return a.x*b.y-a.y*b.x;}
    db dis()
    {return sqrt(x*x+y*y);}
    pnt times(db k)
    {return pnt(k*x,k*y);}
}p[100100];
pnt its(pnt a,pnt b,pnt c,pnt d)//ab和cd交点
{
    db u=(c-a)*(d-a);
    db D=u+(d-b)*(c-b);
    return a+(b-a).times(u/D);
}
int main()
{
    int n;
    scanf("%d",&n);
    for(int i=1;i<=n;++i)
        scanf("%lf%lf",&p[i].x,&p[i].y);
    std::random_shuffle(p+1,p+1+n);
    pnt c(0,0);
    db r=0.0;
    for(int i=1;i<=n;++i)
    {
        if((p[i]-c).dis()<r+eps)
            continue;
        c=p[i];
        r=0.0;
        for(int j=1;j<i;++j)
        {
            if((p[j]-c).dis()<r+eps)
                continue;
            r=(p[j]-p[i]).dis()/2;
            c=p[i]+(p[j]-p[i]).times(.5);
            for(int k=1;k<j;++k)
            {
                if((p[k]-c).dis()<r+eps)
                    continue;
                pnt a=p[j]-p[i];
                pnt f=p[i]+a.times(.5);
                pnt b=p[k]-p[j];
                pnt g=p[j]+b.times(.5);
                c=its(f,f+pnt(-a.y,a.x),g,g+pnt(-b.y,b.x));
                r=(c-p[i]).dis();
            }
        }
    }
    printf("%.10lf\n%.10lf %.10lf\n",r,c.x,c.y);
    return 0;
}