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

前言

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

随机增量法

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

我们考虑递推的过程。

假设我们现在已经求出来前 \(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;
}