O(n)最小圆覆盖算法 (随机增量法)

为什么叫随机增量我也不知道
首先考虑一个点集A,设他的最小覆盖圆为g(A)。 g(A)是存在且唯一的。 (假设存在两个,则必定存在更小覆盖圆)
而且g(A)必定满足以下两者之一:
1. 圆上有三个(或以上)A中的点
2. 圆上有两个A中的点并且其连线为直径。
若上述条件不成立则显然有更小覆盖圆。(感性调整)

就像这样:
这里写图片描述

记K(A)为A中在g(A)上的点集合,这样的点下文叫做关键点。

下面给出一个关键的定理

若点b不在g(A)内,则b属于K(A+b)。

反证:
假设b不属于K(A+b),
则b在g(A+b)内(不是关键点),那么g(A)与g(A+b)是同圆(有没有他都一样)。又因为b不在g(A)内,矛盾。

设g(A,x)为,点集A的最小覆盖圆,并且满足x在圆上。 (x,y不是A中的点)
设g(A,x,y)为,点集A的最小覆盖圆,并且满足x,y在圆上
类似的,我们可以证明:
若点b不在g(A,x)内,则b是g(A,x)的关键点
若点b不在g(A,x,y)内,则b是g(A,x,y)的关键点

先对点随机排序。 这对我们的复杂度分析至关重要。
已知1..i-1的最小覆盖圆C的圆心与半径。
g(1)=圆心为点1,半径为0.

从2开始,我们顺序求g(1..i)。

FirstStep: 对于一个新的i

若i在C中则不需要额外操作。
否则,由上面证明的定理可得i在g(1..i) (点1..i的最小覆盖圆,下面都这么简写了)上。
只要求得g(1..i-1,i)即可得到g(1..i)。,相当于此时进入一个子过程变相求当前g(i)

转化1

初始圆=圆心为点i,半径为0圆。
顺序求g(1..i-1,i)

依次令j=从1到i-1,若点j不在g(1..j-1,i)中(若在则不需要额外操作),则同理得j在g(1..j,i)上。
只要求得g(1..j-1,i,j)即可得到g(1..j,i),相当于此时进入一个子过程变相求当前g(j,i)
(点集为1..j并且i在圆上的最小覆盖圆,表示有点不严谨理解一下)

转化2

初始圆=点i,j为直径所成的圆。
顺序求g(1..j-1,i,j)

依次令k=从1到j-1,若k不在g(1..k-1,i,j)中,则同理得k在g(1..k,i,j)上。

这下我们就知道g(1..k,i,j)上有三个点i,j,k了,不需要继续找在圆上的点了。 因此将圆设为找到的定圆,继续直到求出g(1..j-1,i,j)为止。

时间复杂度的证明

前提是随机序列!
关键分析,第一层来说: 对于一个新的点i,i在不当前圆中的概率为3/i.
证明: 因为g(1..i)只由至多三个关键点决定。这三个关键点是从1..i中选三个。选中i时才会进入下一层 (蜜汁证明….怪怪的但随机情况好像正确)

第二三层类似,然后分析一下复杂度就可以得到期望O(n)的结论。

实现

g是当前圆心,r是当前半径。

    for (int i=1; i<=10*n; i++) { //随机打乱。
        int x=rand()%n+1,y=rand()%n+1;
        swap(p[x],p[y]);
    }
    g=(P) {p[1].x,p[1].y}; r=0;
    for (int i=2; i<=n; i++) if (!zai(p[i])) {
        g=(P){p[i].x,p[i].y}; r=0;
        for (int j=1; j<i; j++) if (!zai(p[j])) {
            g=mid(p[i],p[j]);
            r=sqrt(dis2(p[i],p[j]))*0.5;
            for (int k=1; k<j; k++) if (!zai(p[k])) {
                make(p[i],p[j],p[k]);
            }
        }
    }

make是求出三点所得的圆。 长这样

struct P{db x,y;};
struct L{db k,b;};
L zx(P g,db k,db kx) {
    if (kx==0) return (L){G,g.x}; return (L){k/kx,g.y-k/kx*g.x};
}
P jiaodian(L u,L v) {
    if (v.k==G) swap(u,v);
    db x=0;
    if (u.k==G) x=u.b; else x=(v.b-u.b)/(u.k-v.k);
    return (P) {x,v.b+v.k*x};
}
void make(P a,P b,P c) {
    P yx=jiaodian(zx(mid(a,b),-(b.x-a.x),(b.y-a.y)),zx(mid(a,c),-(c.x-a.x),(c.y-a.y)));
    g=yx; r=sqrt(dis2(yx,a));
}

就是求中垂线交点。 注意特判与y轴平行的情况。 (本来想打向量的但是就这么一丢丢特判一下算了。。。)

猜你喜欢

转载自blog.csdn.net/jokerwyt/article/details/79221345