[Poj 2187]计算几何之凸包(三) {旋转卡壳初步}

{

上一节介绍了凸包的高效算法

和一个最远点对的应用

这一段将更好的解决最远点对问题

}

(若不做特殊说明 下文讨论的问题均是在欧氏空间

 若不做特殊说明 下文中距离均是指空间中欧氏距离)

==============================

一.简单枚举算法的不足

上一次介绍了一个基本的求平面最远点对的算法

即先求点集的凸包 然后枚举凸包上的点来求最远点集

这是利用了凸包上的点相比 点集中的点 一般是很少的 平均情况很好 并且我们也能AC这个问题

但是这是有局限性的 当凸包上的点达到O(N)的级别时 凸包的优化作用就不存在了

不过我们还要考虑到 凸包还起了对凸包上点集排序的作用

凸包有很多的优美的性质 我们可以加以利用 以得到更加高效的算法

旋转卡壳算法就是利用凸包特性的一类解决问题的方法

==============================

二.旋转卡壳算法

旋转卡(qiǎ)壳算法(Rotating Calipers Algorithm):

是解决一些与凸包有关问题的有效算法 就像一对卡壳卡住凸包旋转而得名

Every time one blade of the caliper lies flat against an edge of the polygon, it forms an antipodal pair with the point or edge touching the opposite blade. It turns out that the complete "rotation" of the caliper around the polygon detects all antipodal pairs and may be carried out in O(n) time.

http://en.wikipedia.org/wiki/Rotating_calipers

(图片来自:http://cgm.cs.mcgill.ca/~orm/rotcal.html)

被一对卡壳正好卡住的对应点对称为对踵点(Antipodal point)

http://en.wikipedia.org/wiki/Antipodal_point

可以证明对踵点的个数不超过3N/2个 也就是说对踵点的个数是O(N)的

对踵点的个数也是我们下面解决问题时间复杂度的保证

上第一个图是卡壳的一般情况 卡住两点 图二是卡住一条边和一个点

由于实现中 卡住两点的情况不好处理 我们通常关注第二种情况

在第二种情况中 我们可以看到 一个对踵点和对应边之间的距离比其他点要大

也就是一个对踵点和对应边所形成的三角形是最大的 下面我们会据此得到对踵点的简化求法

看一下官方的伪代码:

当时我看完了 就一个字 长... 我最讨厌冗长的程序了...

begin
     p0:=pn;
     q:=NEXT[p];
     while (Area(p,NEXT[p],NEXT[q]) > Area(p,NEXT[p],q)) do
          q:=NEXT[q];
          q0:=q;
          while (q != p0) do
               begin
                    p:=NEXT[p];
                    Print(p,q);
                    while (Area(p,NEXT[p],NEXT[q]) > Area(p,NEXT[p],q) do
                         begin
                              q:=NEXT[q];
                              if ((p,q) != (q0,p0)) then Print(p,q)
                              else return
                         end;
                    if (Area(p,NEXT[p],NEXT[q]) = Area(p,NEXT[p],q)) then
                      if ((p,q) != (q0,p0)) then Print(p,NEXT[q])
                      else Print(NEXT[p],q)
               end
end.

几经折腾 终于找到了一个不错的实现:http://www.cnblogs.com/DreamUp/archive/2010/09/16/1828131.html

不过不是很好理解 这里作一下说明

ch[m+1]:=ch[1]; j:=2;
 for i:=1 to m do
    begin
    while cross(ch[i],ch[j],ch[i+1])<cross(ch[i],ch[j+1],ch[i+1]) do
        begin inc(j); if j>m then j:=1; end;
    writeln(ch[i].x,' ',ch[i].y,' ',ch[j].x,' ',ch[j].y);
    end;

上面就是旋转卡壳寻找对踵点的过程

其中叉积函数Cross(A,B,C:Point):Real 返回AB到AC的二维定义下的叉积

这里主要用到了叉积求三角形面积的功能

我们对于一条对应边<CH i,CH Next[i]>求出距离这条边最远的点CH j

则由上面第二种情况可知 CH i 和 CH j 为一对对踵点 这样让 CH i 绕行凸包一周即可得到所有的对踵点

下面面这个图 由于本人的gif图制作水平拙劣 所以不好看

需要的可以下载几何画板察看原版GSP文件 点击这里下载GSP文件

接下来考虑 如何得到距离每条对应边的的最远点呢?

稍加分析 我们可以发现 凸包上的点依次与对应边产生的距离成单峰函数

具体证明可以从凸包定义入手 用反证法解决

这样我们再找到一个点 使下一个点的距离小于当前的点时就可以停止了

而且随着对应边的旋转 最远点也只会顺着这个方向旋转 我们可以从上一次的对踵点开始继续寻找这一次的

由于内层while循环的执行次数取决于j增加次数 j最多增加O(N)次

所以求出所有对踵点的时间复杂度为O(N)

还有有两点需要注意:

1.上面这段代码及代码的分析都是需要凸包上没有三点共线的

2.Next[i] 不需要手动求 在原代码中有很好的处理

最后指出网上很多文章的一个错误 一个点的对踵点并不是离这个点最远的点!

这样子的点对是根本不满足对踵点的性质的 即最为重要的单峰分布性质

下图是一个反例:

==============================

三.旋转卡壳算法的简单应用

至此我们终于可以更高效的解决平面最远点对问题了

有一个很重要的结论是 最远点对必然属于对踵点对集合

那么我们先求出凸包 然后求出对踵点对集合 然后选出距离最大的即可

用这个算法可以47ms AC这个问题 算上凸包的时间 总复杂度为O(Nlog2N)

代码如下:

{$inline on}
{$optimization on}
const    maxn=50000;
type    point=record x,y:longint; end;
var    n,i,x,m,ans,j:longint;
    ch,p:array[1..maxn+1]of point;
    s:array[1..maxn]of longint;
function cross(a,b,c:point):longint; inline;
begin
cross:=(b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
end;
function dist(a,b:point):longint; inline;
begin
dist:=sqr(a.x-b.x)+sqr(a.y-b.y);
end;
function cmp(a,b:point):boolean; inline;
begin
cmp:=(a.x<b.x)or(a.x=b.x)and(a.y<b.y);
end;
function max(a,b:longint):longint;
begin
if a>b then max:=a else max:=b;
end;
procedure swap(a,b:longint); inline;
var    x:point;
begin
x:=p[a]; p[a]:=p[b]; p[b]:=x;
end;
procedure hull(l,r:longint; a,b:point);
var    x,i,j,k:longint;
    y:point;
begin
x:=l; y:=p[l];
for k:=l to r do
    if (s[x]<s[k])or(s[x]=s[k])and(cmp(y,p[k]))
        then begin x:=k; y:=p[k]; end;
i:=l-1; j:=r+1;
for k:=l to r do
    begin
    inc(i); s[i]:=cross(p[k],a,y);
    if s[i]>0 then swap(i,k) else dec(i);
    end;
for k:=r downto l do
    begin
    dec(j); s[j]:=cross(p[k],y,b);
    if s[j]>0 then swap(j,k) else inc(j);
    end;
if l<=i then hull(l,i,a,y);
inc(m); ch[m]:=y;
if j<=r then hull(j,r,y,b);
end;
begin
assign(input,'Maxd.in'); reset(input);
assign(output,'Maxd.out'); rewrite(output);
readln(n);
for i:=1 to n do
    begin
    readln(p[i].x,p[i].y);
    if (x=0)or cmp(p[i],p[x]) then x:=i;
    end;
swap(1,x);
m:=1; ch[1]:=p[1]; hull(2,n,p[1],p[1]);
ch[m+1]:=ch[1]; j:=2; ans:=0; 
for i:=1 to m do
    begin
    while cross(ch[i],ch[j],ch[i+1])<cross(ch[i],ch[j+1],ch[i+1]) do
        begin inc(j); if j>m then j:=1; end;
    ans:=max(ans,dist(ch[i],ch[j]));
    end;
writeln(ans);
close(input); close(output);
end.

下一节介绍旋转卡壳的更多应用

之后开始介绍一点3D凸包