平面上圆的凸包算法
我们之前探讨过这个有趣的问题:
平面上有若干圆,求包含这些圆的所有凸集的交。
根据之前讨论的结果,直接按圆心排序做扫描线的方法是错误的。我们需要考虑圆上的每一个点是否可以作为凸包上的一部分。
然而圆上的点的数目是无限多的。我们需要借助离散化的思想:因为我们发现凸包一定是由若干圆弧和线段构成的。而其中线段一定是切线段。由此,我们只需要将每两两圆的切点取出来求一个凸包。显然,在圆凸包上出现的切线段一定会在切点凸包中出现;而切点凸包中其余的线段则一定是弧弦。
但是,这个算法需要枚举每两个圆之间的切线。
但实际上,真的有那么多点会用到吗?
显然不是。
完整的切点信息太冗余;尝试用圆来描述凸包。
如果存在半平面包含所有圆且和圆A圆B相切,并且方向从圆A指向圆B,我们则称AB是一条凸包的边;这样以来,我们便可以用一个圆的序列表示出整个凸包的轮廓。
我们设\(H\)为这个最终序列,则对于任何小于序列长度的\(i\)有\(H_iH_{i+1}\)是凸包的边;同理,夹在切线\(H_{i-1}H_i\)和切线\(H_iH_{i+1}\)之间在\(H_i\)上的弧也是凸包的一部分。
接下来,我们考虑序列长度的上界,这同时也和凸包中切线以及圆弧的个数的上界;
显然序列\(H\)有以下性质:
- 每相邻两个元素不相等;
- 不存在\(\cdots A\cdots B \cdots A\cdots B\cdots\) 这样的子序列
我们可以证明,任何满足上面两个性质的一个值域为1~n的整数序列长度都不会超过\(2n-1\)
证明:
设一个存在长度为2n+d的序列,满足性质一二。其中d为非负整数。
如果该序列未出现1~n中某一个元素,则将2n个元素中的重复元素换成未出现的元素,直到每一种元素都出现。经过了这一过程,序列一定仍然合法。
之后,序列中一定有n个互不相同的元素。剩下的另外n+d的个中,假设最多有m个互相不同的。
我们将两两相同的元素之间都连一条线,显然这些线之间可以有包含关系,但是不能交错。对于这些只有包含和不相交关系的区间,我们可以插入n-m个元素,来保证其相同的元素都被隔开。但是实际上我们会有n+d-m+1个相同元素相邻。显然是不可能合法的。这与原设矛盾。
综上所述,2n-1是序列长度的一个上界。
我们保证了对于任意N个圆,都有其凸包序列长度小于2N,那么怎么算出整个集合S的凸包呢?
考虑如果我们已经计算出了P集合与Q集合的两个凸包序列,要想合并两个序列,我们只需要像旋转卡壳一样,用两个方向相同的且与分别和P,Q相切的半平面来判断任意时刻外侧的圆,这样来构造新的序列。
顺着这样的思路,我们考虑分治:
//构造集合S的凸包序列;
序列 分治( 圆集合 S ) {
如果( |S| == 1 ) {
返回 S中唯一元素单独组成的序列;
}
将S集合按任意的顺序分成大小最多相差1的不相交的两个子集P,Q;
序列 A = 分治(P), B = 分治(Q);
序列 Final = 合并(A,B);
返回 Final;
}
what?这是啥?易语言??
其中 合并函数 的伪代码如下:
//合并两个序列的函数;
序列 合并( 序列 P, 序列 Q ) {
序列 S = 空;
有向直线 L = 随意, Lp = L平行并且和P序列相切的有向直线, Lq = L平行并且和Q序列相切的有向直线;
序列元素 p = Lp相切的圆,q = Lq相切的圆;
循环(P,Q中存在没有访问过的元素) {
如果(Lp在Lq外侧) {
S序列后面添加p;
更新(L,p,q);
} 否则 {
S序列后面添加q;
更新(L,q,p);
}
Lp = 和L平行的和p相切的有向直线;
Lq = 和L平行的和q相切的有向直线;
}
返回 S;
}
其中 更新函数 的内容如下:
角度 夹角( 有向直线 a, 有向直线 b ) {
如果(b不存在) 返回 无穷大;
返回 有向直线a逆时针旋转到有向线段b的角度;
}
有向直线 切线( 序列元素 a, 序列元素 b ) {
返回 圆a到圆b的半平面切线;
}
//将L以及p,q跟新到合适的位置;
void 更新( 有向直线& L, 序列元素& x, 序列元素& y ) {
角度 a1 = 夹角(L, 切线(x,y)), a2 = 夹角(L, 切线(x,序列中x的下一个元素));
角度 a3 = 夹角(L, 切线(y,x)), a4 = 夹角(L, 切线(y,序列中y的下一个元素));
如果( a1是a1,a2,a3中最小的 ) {
S序列后面添加x;
如果( a3是a2,a3,a4中最小的 ) {
S序列后面添加y;
}
}
如果(a2 < a4) {
L = 切线(x,序列中x的下一个元素);
x = 下一个元素;
} 否则 {
L = 切线(y,序列中y的下一个元素);
y = 下一个元素;
}
}
这里还有一些关于算法正确性和复杂度的证说明:
- 关于更新函数,因为我们将p,q中的外侧圆作为序列参数中的前一个,所以在第一大块的判断可以处理完p和q的连边关系;
- 关于复杂度,显然,由于S序列长度不超过2n-1;则分治的复杂度为\(T(n) = 2T(n/2)+O(n)=O(n\log n)\)
- 因为点集的凸包算法复杂度下界是NlogN,而点集凸包算法实际上又是这个问题半径无穷小的特殊情况,所以NlogN是该问题时间复杂度的一个下届。
后记
我们之所以不采用增量法,原因是圆的大小不一,因而不方便组织插入圆的顺序;及其容易出现删点不是连续区间之内的问题而难以组织数据来做到快速删除并删除;而分治算法则充分利用了n圆序列长度不超过2n-1的性质,提高了旋转卡壳的信息利用。
引用
[1] David Rappaport, A convex hull algorithm for discs, and applications , 1991:173-177
[2] [数据已删除]的超级大脑, 2000: memory access failed-memory access failed