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

作者: wjyyy 分类: 学习笔记,计算几何,递推 发布时间: 2019-03-12 21:22

点击量:193

前言

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

随机增量法

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

我们考虑递推的过程。

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

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

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

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

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

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

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

复杂度分析

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

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

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

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

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

总复杂度为 $latex 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;
}

3
说点什么

avatar
3 Comment threads
0 Thread replies
0 Followers
 
Most reacted comment
Hottest comment thread
0 Comment authors
Recent comment authors
  Subscribe  
最新 最旧 得票最多
提醒
trackback

… [Trackback]

[…] Find More on to that Topic: wjyyy.top/3347.html […]

trackback

… [Trackback]

[…] Information on that Topic: wjyyy.top/3347.html […]

trackback

… [Trackback]

[…] There you can find 94082 more Information to that Topic: wjyyy.top/3347.html […]

/* */