计算几何之二维凸包：卷包裹算法、Graham Scan Algorithm、旋转卡壳算法
2016-08-24 14:27
501 查看
转自：http://www.cnblogs.com/Booble/
原文链接：
http://www.cnblogs.com/Booble/archive/2011/02/28/1967179.html#3229390
http://www.cnblogs.com/Booble/archive/2011/03/10/1980089.html
http://www.cnblogs.com/Booble/archive/2011/04/03/2004865.html
第一篇：
一.凸集&凸包
(下文中所有的集合 若不作特殊说明 都是指欧氏空间上的集合)
凸集(Convex Set):任意两点的连线都在这个集合内的集合就是一个凸集.
A set in Euclideanspace is convex set if it contains all the line segments connecting anypair of its points.
http://mathworld.wolfram.com/Convex.html
凸包(Convex Hull):包含集合S的所有凸集的交集就是集合S的凸包.
The convex hull ofa set of points S in N dimensions is the intersection of all convex setscontaining S.
我们经常关注一个点集的凸包 这也是计算几何学的一个基本问题
我们现在已经有成熟的算法可以求出平面点集的凸包和空间点集的凸包
甚至有的算法可以方便的求出任意维度欧氏空间内的一个点集的凸包
凸包有着优美而实用的性质 我们可以利用凸包把一个杂乱的点集所包含的信息进行有效的概括 梳理
====================================================================
二.平面点集的凸包的算法
(下文中所有凸包 若不作特殊说明 都是指平面点集的凸包)
有两种直观的理解凸包的方式
在木板上钉钉子 用一个有弹性的橡皮筋去框住所有钉子
橡皮筋形成的图形就是这个钉子所构成的点集的凸包
还有一种理解我们用一根麻绳绑住一个外面的钉子 然后拉着麻绳绕所有钉子一圈
这个麻绳最后也构成了点集的凸包
其中 第二种理解是我们一个经典算法 卷包裹法(Gift Wrapping)的思路
卷包裹算法从一个必然在凸包上的点开始向着一个方向依次选择最外侧的点
当回到最初的点时 所选出的点集就是所要求的凸包
这里还有两个问题不是很清楚:
1.怎么确定一个肯定在凸包上的点?
这个问题很好解决 取一个最左边的也就是横坐标最小的点
如果有多个这样的点 就取这些点里 纵坐标最小的
这样可以很好的处理共线的情况
2.如何确定下一个点(即最外侧的点)?
我们需要利用向量的叉积来解决这个问题
-----------------------------------------------------------------------
向量的叉积(Cross Product)原本是三维空间中的问题 在二维中也有巧妙的应用
http://mathworld.wolfram.com/CrossProduct.html
(下文中所有的叉积 若不作特殊说明 都是指二维中新定义的叉积
下文中所有的向量乘法 若不作特殊说明 都是指向量的叉积)
我们定义二维向量<x1,y1>和<x2,y2>的叉积为一个实数Cp=x1*y2-x2*y1
叉积有两条性质很常用 也很好用
1.叉积的一半是一个三角形的有向面积
这个公式可以避免面积计算的误差 如果点是整点 那么所有运算都是整数
2.向量的叉积的符号代表着向量旋转的方向
向量的叉积是不满足交换律的
向量A乘以向量B 如果为正则为A逆时针旋转向B 否则为顺时针
当然这里A转向B的角总是考虑一个小于180度以内的角 否则就会出错
-----------------------------------------------------------------------
有了向量 我们就可以选取一个最外侧的点了
比如现在我们卷包裹卷到J点我们要选取一个最外侧的点
当然比较红色的到角可以直接得到最外侧的点 不过不方便
我们可以考虑那个蓝色的到角
利用向量 我们可以比较哪个点"更外侧"
比如点K和点I 我们利用向量JK乘以向量JI得到一个数 这个数应该是负数 说明I比K更外侧
两个向量的比较具有传递性 所以我们可以像N个数里取最大的数一样取出最外侧的
遍历所有点 每个点都和现有最外侧的点比较 得到新的最外侧的点
至此两个问题都得以解决 我们可以写出满足一般要求的卷包裹算法了
不过还遗留有一个问题 就是处理共线的问题
有时候我们需要凸包边上的点也考虑到 有时候却需要去掉这些点
我们通常称在凸包顶点处的点为极点
如果我们只要求保留极点而去除在边上的点
我们只需在取外侧的点的时候 碰到共线的点取最远的
相反 如果我们要保留所有在边上的点我们只需要在共线的点中取最近的
这样整个卷包裹法终于完成了
给出完整的代码:
GiftWrapping
GiftWrapping
{$inline on}
const zero=1e-6; maxn=100000;
type point=record x,y:extended; end;
var p:array[1..maxn]of point;
ch:array[1..maxn]of longint;
temp,n,m,i,j,k:longint;
function sgn(x:extended):longint; inline;
begin
if abs(x)<zero then exit(0);
if x<0 then sgn:=-1 else sgn:=1;
end;
function cross(a,b,c:point):extended; 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):extended; inline;
begin
dist:=sqr(a.x-b.x)+sqr(a.y-b.y);
end;
function cmp(a,b,c:point):boolean; inline;
var temp:longint;
begin
temp:=sgn(cross(a,b,c));
if temp<>0 then exit(temp<0); {*B}
cmp:=dist(a,b)<dist(a,c); {*A}
end;
begin
assign(input,'Hull.in'); reset(input);
assign(output,'Hull1.out'); rewrite(output);
readln(n);
for i:=1 to n do
begin
readln(p[i].x,p[i].y);
if (j=0)or(p[i].x<p[j].x)or
(sgn(p[i].x-p[j].x)=0)and(p[i].y<p[j].y)
then j:=i;
end;
temp:=j;
while true do
begin
k:=0;
inc(m); ch[m]:=j;
for i:=1 to n do
if (i<>j)and((k=0)or
cmp(p[j],p[k],p[i]))
then k:=i;
if k=temp then break;
j:=k;
end;
for i:=1 to m do
writeln(p[ch[i]].x:0:2,' ',p[ch[i]].y:0:2);
close(input); close(output);
end.
*A Change Direction
*B Remove Colinear Points
还有两点要补充说明一下:
1.卷包裹算法的复杂度是O(NH)
N是全部的点数 H是最终在凸包上的点数
所以卷包裹算法很适合凸包上的点很少的时候 通常随机数据很快
但是构造出的凸包上的点很多的数据 这个算法就会很慢
比如所有的点都在一个圆周上
2.卷包裹算法输出的点是有序的
这也是对二维凸包算法的一个基本要求
通常只有保证有序才能进行进一步的计算
通过改变CMP函数可以改变上文中提到的共线(取/不取)以及这里的序(顺时针/逆时针)的问题
====================================================================
三.凸包的面积和周长
凸包的一个简单算法介绍完了
来看一个具体的问题 Poj 1113 http://poj.org/problem?id=1113
问题是求到一个简单多边形距离最少为L的点的轨迹的周长
我们可以先求凸包 然后得到计算公式
实际上 圆的部分可以拼接起来变成一个整圆
这样就好算多了 注意把共线的都去掉 不然误差会很大的
代码就不贴了 很简单的
下一篇介绍更快速的平面凸包算法
Graham Scan 和 QuickHull
强烈推荐QuickHull
第二篇：
一.卷包裹算法(Gift Wrapping Algorithm)的特性
前面提到过卷包裹算法的复杂度问题
由于卷包裹算法是两重循环实现的 因此很好分析它的复杂度
1 while true
do
2 begin
3 k:=0;
4 inc(m); ch[m]:=j;
5 for i:=1
to n do
6 if (i<>j)and((k=0)or
7 cmp(p[j],p[k],p[i]))
8 then k:=i;
9 if k=temp
then break;
10 j:=k;
11 end;
内部循环N次 外循环的次数决定于凸包上的点数H
所以卷包裹算法复杂度为O(HN) 这个复杂度很有特点
考虑一个随机点集的凸包上的点数往往很少 但是最坏情况下是O(N)级别的
比如所有的点都在一个圆上的时候
这就决定了卷包裹算法很适合随机的点集 但是如果是刻意构造的数据或是比较特殊的数据
就会达到O(N^2)的最坏复杂度 这是我们不愿意看到的
下面介绍最坏情况下复杂度更好的Graham扫描算法(Graham ScanAlgorithm)
====================================================================
二.Graham扫描算法(Graham ScanAlgorithm)
Graham扫描算法维护一个凸壳 通过不断在凸壳中加入新的点和去除影响凸性的点 最后形成凸包
The Graham scan isa method of computing the convex hull of a finite set of points in the planewith time complexity O(n log n). It is named after Ronald Graham, who publishedthe original algorithm in 1972.[1] The algorithm finds all vertices
of theconvex hull ordered along its boundary.
http://en.wikipedia.org/wiki/Graham_scan
算法主体由两部分组成 先是排序 后是扫描 分块讲解一下
----------------------------------------------------------------------------------------------------------
1.点集排序
为了得到加入新点的顺序 Graham扫描法的第一步是对点集排序
排序是对杂乱的点集进行了梳理 这也是这种算法能够得到更高效率的根本原因
排序的方法也有两种 极角坐标排序(极角序) 和 直角坐标排序(水平序)
前者好理解一些 但是在实现的时候 后者更方便
先说极角序 为了极角排序 我们先得得到一个参考点
一般的 我们取最左边(横坐标最小)的点作为参考点 如果有多个这样的点就取最下面的(纵坐标最小)
看这样一个例子 这是一个任意给出的平面点集:
参考点的定义:在横坐标最小的情况下取纵坐标最小的点
所以所有的点只能在这个黄色的半平面中 而且正上方为闭(可取得) 正下方为开(不可取)
这就决定了参考点的性质:点集中任意两点和参考点所成的到角为锐角
这样我们取得参考点 然后再考虑极角排序
极角排序以参考点为极角坐标系原点 各个点的极角为关键字
由于上面我们得到的参考点的性质 我们可以设所有点的极角均在(-90,90]之间
排序完成后应该是这样的:
比较极角我们仍然可以利用向量的叉积
叉积在这里已经介绍了 http://www.cnblogs.com/Booble/archive/2011/02/28/1967179.html
同样由于参考点的性质 所有向量之间的到角都是在180度以内 不会产生错误
----------------------------------------------------------------------------------------------------------
2.Graham的栈扫描
Graham的扫描是一个很优美的过程 用到的数据结构也很简单 仅仅是一个栈而已
核心的思想是按照排好的序 依次加入新点得到新的边
如果和上一条边成左转关系就压栈继续 如果右转就弹栈直到和栈顶两点的边成左转关系 压栈继续
实现的时候我们不用存边 只需要含顺序在栈里存点 相邻两点就是一条边
由于我们时时刻刻都保证栈内是一个凸壳 所以最后扫描完毕 就得到了一个凸包
下面还是继续上面的那个样例 演示一下栈扫描的过程
这样Graham扫描算法基本完成
复杂度是排序O(Nlog2N) 扫描O(N) {每个点仅仅出入栈一次}
合起来是一个O(Nlog2N)的算法 很优秀
----------------------------------------------------------------------------------------------------------
3.双重共线点难题
和卷包裹算法一样 我们同样还要考虑共线点问题
而且Graham扫描算法的共线问题更复杂 所以需要仔细考虑
i).排序时的共线问题
如果极角相同 我们应该怎么定先后呢?
我们得加上第二关键字距离 比如极角相同 距离参考点近的先
不过不管是近的先还是 远的先 开始和结束的两条边总是矛盾的
我们必须对其中一条特殊处理 除了结束边外距离近的先 结束边上距离远的先
这就是为什么极角排序不是很好实现的原因了 下面会介绍一下水平序
ii).扫描时的共线问题
这个和卷包裹算法的处理放法如出一辙
如果需要保留共线的点就在到角相同时取距离最近的
如果仅仅需要凸包极点就取距离最远的
----------------------------------------------------------------------------------------------------------
4.直角坐标排序(水平序)
直角坐标排序方法没有了极角排序的不足
以横坐标为第一关键字 纵坐标为第二关键字 排序点集
然后从第一个点开始 分别利用Graham扫描生成左链和右链
需要注意以下两点:
i).左链和右链的旋转方向是相反的
ii).注意生成第二条链要忽略第一条链上的点
由于水平序的CMP函数比较简单 代码也更短
还有一件有意思的事 Graham扫描算法是1972年提出的 卷包裹算法是1973年提出的
其实不奇怪 这两个算法本来也没有优劣所用的思想不同各自善于处理不同的情况
Graham扫描所用时间在点集随机时 还不如卷包裹算法快
正如第一节所说 卷包裹算法已经可以处理大部分情况 也是一个可取的算法
实际应用中点集更趋向于均匀 而不是集中在凸包上
----------------------------------------------------------------------------------------------------------
最后给一下比较难实现的极角排序代码
GrahamScan
{$inline on}
const zero=1e-6; maxn=100000;
type point=record x,y:extended; end;
var p,ch:array[1..maxn]of point;
i,j,n,m,temp:longint;
function sgn(x:extended):longint; inline;
begin
if abs(x)<zero then exit(0);
if x<0 then sgn:=-1 else sgn:=1;
end;
function cross(a,b,c,d:point):extended; inline;
begin
cross:=(b.x-a.x)*(d.y-c.y)-(d.x-c.x)*(b.y-a.y);
end;
function dist(a,b:point):extended; inline;
begin
dist:=sqr(a.x-b.x)+sqr(a.y-b.y);
end;
function cmp1(a,b,c:point):boolean; inline;
var temp:longint;
begin
temp:=sgn(cross(a,b,a,c));
if temp<>0 then exit(temp>0); {*B}
cmp1:=dist(a,b)<dist(a,c);
end;
function cmp2(a,b,c:point):boolean; inline;
var temp:longint;
begin
temp:=sgn(cross(a,b,b,c));
if temp<>0 then exit(temp<0); {*B}
cmp2:=dist(a,b)<dist(a,c); {*A}
end;
procedure swap(var a,b:point); inline;
var temp:point;
begin
temp:=a; a:=b; b:=temp;
end;
procedure sort(l,r:longint);
var i,j:longint;
x:point;
begin
i:=l; j:=r;
x:=p[(l+r)>>1];
repeat
while cmp1(p[1],p[i],x) do inc(i);
while cmp1(p[1],x,p[j]) do dec(j);
if not(i>j)
then begin
swap(p[i],p[j]);
inc(i); dec(j);
end;
until i>j;
if i<r then sort(i,r);
if l<j then sort(l,j);
end;
begin
assign(input,'Hull.in'); reset(input);
assign(output,'Hull2.out'); rewrite(output);
readln(n);
for i:=1 to n do
begin
readln(p[i].x,p[i].y);
if (p[i].x<p[1].x)or
(sgn(p[i].x-p[1].x)=0)and(p[i].y<p[1].y)
then swap(p[1],p[i]);
end;
sort(2,n);
while i>2 do
begin
if sgn(cross(p[1],p[i],p[1],p[i-1]))<>0
then break; dec(i);
end;
temp:=(i+n+1)>>1;
for j:=n downto temp do
swap(p[j],p[i+n-j]);
m:=1; ch[1]:=p[1];
inc(n); p
:=p[1];
for i:=2 to n do
begin
while m>1 do
begin
if not cmp2(ch[m-1],ch[m],p[i])
then break; dec(m);
end;
inc(m); ch[m]:=p[i];
end;
for i:=1 to m-1 do
writeln(ch[i].x:0:2,' ',ch[i].y:0:2);
close(input); close(output);
end.
*A Change Direction
*B Remove Colinear Points
====================================================================
三.快速凸包算法(Quickhull Algorithm)
对比Graham扫描算法和卷包裹算法
我们发现 Graham扫描算法在凸包上的点很密集时仍然适用
卷包裹算法在凸包上点集随机分布时是很高效的
那么有没有两个优点都具备的算法呢?
是有的! 快速凸包算法(Quickhull Algorithm)就是这样的一个算法
快速凸包算法是一个和快速排序(Quicksort Algorithm)神似的算法
尽管快速排序的最坏复杂度可以达到O(N^2)
但是有着极小的常数 实现方便 思路优美 绝大多数情况特别高效的快速排序 还是赢得了更多人的青睐
快速凸包算法也是这样 尽管可以构造一个数据使之达到O(N^2)的复杂度
但这需要刻意针对程序经过分析 才能做到 是实际应用中根本不会碰到的情况
在点集均匀分布时 快速凸包的复杂度更是达到了O(N) 是上面两种算法难以企及的
在绝大多数情况下 平均复杂度是O(Nlog2N) 也很高效
快速凸包继承了快速排序分治的思想 这是一个递归的过程
------------------------------------------------------------------------------------
伪代码如下:
1 void 快速凸包(P:点集 , S:向量
/*S.p,S.q:点)*/ ){
2 /* P 在 S
左侧*/
3 选取 P 中距离 S 最远的 点 Y ;
4 向量 A <- { S.p , Y }; 向量 B <- { Y , S.q } ;
5 点集 Q <- 在 P 中 且在 A 左侧的点 ;
6 点集 R <- 在 P 中 且在 B 左侧的点 ;
/* 划分 */
7 快速凸包 ( Q , A ) ; /*
分治 */
8 输出 (点 Y) ; /*
按中序输出 保证顺序*/
9 快速凸包 ( P , B ) ; /*
分治 */
10 }
初始化就选取 最左下和最右上的点 划分好 然后调用两次快速凸包函数分别求上下凸包
其中 选取和划分都需要用到向量的叉乘 注意方向
另外 存储点集可以用数组里的连续一段 传参的时候就传左右下标
划分的时候就是给数组交换元素 使新的点集 变成连续的一段
这里有一个很好的演示
http://www.cs.princeton.edu/courses/archive/spr10/cos226/demo/ah/QuickHull.html
还要补充说明一下
快速凸包在所有点都在圆周上的时候还是O(Nlog2N) 不过会比Graham扫描算法慢一些
可以说 这种数据是Graham扫描法的最好情况 一遍走完就行了
构造快速凸包的最坏情况就是使划分不均等 和构造快速排序最坏情况一样
贴以下我很辛苦写出来的代码吧 为了好看些 还心血来潮用cpp写了...
QuickHull
#include <iostream>
#include <math.h>
#define maxn 100000
#define zero 1e-12
#define sgn(x) (fabs(x)<zero?0:(x>0?1:-1))
#define cross(a,b,c) ((b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x))
#define cmp(a,b) (a.x<b.x || sgn(a.x-b.x)==0 && a.y<b.y)
using namespace std;
struct point{
double x,y;
}p[maxn];
double s[maxn];
void hull(int l,int r,point a,point b){
int x=l,i=l-1,j=r+1,k;
for (k=l;k<=r;k++){
double temp=sgn(s[x]-s[k]);
if (temp<0 || temp==0 && cmp(p[x],p[k])) x=k;
}
point y=p[x];
for (k=l;k<=r;k++){
s[++i]=cross(p[k],a,y);
if (sgn(s[i])>0) swap(p[i],p[k]); else i--;
}
for (k=r;k>=l;k--){
s[--j]=cross(p[k],y,b);
if (sgn(s[j])>0) swap(p[j],p[k]); else j++;
}
if (l<=i) hull(l,i,a,y);
printf("%.4lf %.4lf\n",y.x,y.y);
if (j<=r) hull(j,r,y,b);
}
int main(){
freopen("CH2D.in","r",stdin);
freopen("CH2D1.out","w",stdout);
int n,i,x=0;
scanf("%d",&n);
for (i=1;i<=n;i++){
scanf("%lf %lf",&p[i].x,&p[i].y);
if (x==0 || cmp(p[i],p[x])) x=i;
}
swap(p[1],p[x]);
printf("%.4lf %.4lf\n",p[1].x,p[1].y);
hull(2,n,p[1],p[1]);
return 0;
}
====================================================================
四.凸包算法复杂度下界
(引自<算法艺术与信息学竞赛>)
====================================================================
五.高效算法的应用
这里有一个平面最远点对的问题 可以利用凸包解决
Poj 2187 http://poj.org/problem?id=2187
由于最远点对必然在凸包上
我们先求凸包 然后枚举凸包上的点 不过这个复杂度最坏是O(N^2)的
不过凸包在很多情况下会改善问题的平均复杂度
凸包上的点通常很少 所以这个问题也可以过
篇幅关系 在下一节会介绍降低最坏复杂度的旋转卡壳算法
====================================================================
文中部分图片来源 http://westhoffswelt.de/blog/0040_quickhull_introduction_and_php_implementation.html http://en.wikipedia.org/wiki/Graham_scan
第三篇：
{
上一节介绍了凸包的高效算法
和一个最远点对的应用
这一段将更好的解决最远点对问题
}
(若不做特殊说明 下文讨论的问题均是在欧氏空间
若不做特殊说明 下文中距离均是指空间中欧氏距离)
==============================
一.简单枚举算法的不足
上一次介绍了一个基本的求平面最远点对的算法
即先求点集的凸包 然后枚举凸包上的点来求最远点集
这是利用了凸包上的点相比 点集中的点 一般是很少的 平均情况很好 并且我们也能AC这个问题
但是这是有局限性的 当凸包上的点达到O(N)的级别时 凸包的优化作用就不存在了
不过我们还要考虑到 凸包还起了对凸包上点集排序的作用
凸包有很多的优美的性质 我们可以加以利用 以得到更加高效的算法
旋转卡壳算法就是利用凸包特性的一类解决问题的方法
==============================
二.旋转卡壳算法
旋转卡(qiǎ)壳算法(Rotating CalipersAlgorithm):
是解决一些与凸包有关问题的有效算法 就像一对卡壳卡住凸包旋转而得名
Every time oneblade of the caliper lies flat against an edge of the polygon, it forms anantipodal pair with the point or edge touching the opposite blade. It turns outthat the complete "rotation" of the caliper around the polygondetects 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
不过不是很好理解 这里作一下说明
1 ch[m+1]:=ch[1]; j:=2;
2 for i:=1
to m do
3 begin
4 whilecross(ch[i],ch[j],ch[i+1])<cross(ch[i],ch[j+1],ch[i+1])
do
5 begininc(j);
if j>m then j:=1;
end;
6 writeln(ch[i].x,' ',ch[i].y,' ',ch[j].x,' ',ch[j].y);
7 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)
代码如下:
Maxd
{$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凸包
原文链接：
http://www.cnblogs.com/Booble/archive/2011/02/28/1967179.html#3229390
http://www.cnblogs.com/Booble/archive/2011/03/10/1980089.html
http://www.cnblogs.com/Booble/archive/2011/04/03/2004865.html
第一篇：
一.凸集&凸包
(下文中所有的集合 若不作特殊说明 都是指欧氏空间上的集合)
凸集(Convex Set):任意两点的连线都在这个集合内的集合就是一个凸集.
A set in Euclideanspace is convex set if it contains all the line segments connecting anypair of its points.
http://mathworld.wolfram.com/Convex.html
凸包(Convex Hull):包含集合S的所有凸集的交集就是集合S的凸包.
The convex hull ofa set of points S in N dimensions is the intersection of all convex setscontaining S.
我们经常关注一个点集的凸包 这也是计算几何学的一个基本问题
我们现在已经有成熟的算法可以求出平面点集的凸包和空间点集的凸包
甚至有的算法可以方便的求出任意维度欧氏空间内的一个点集的凸包
凸包有着优美而实用的性质 我们可以利用凸包把一个杂乱的点集所包含的信息进行有效的概括 梳理
====================================================================
二.平面点集的凸包的算法
(下文中所有凸包 若不作特殊说明 都是指平面点集的凸包)
有两种直观的理解凸包的方式
在木板上钉钉子 用一个有弹性的橡皮筋去框住所有钉子
橡皮筋形成的图形就是这个钉子所构成的点集的凸包
还有一种理解我们用一根麻绳绑住一个外面的钉子 然后拉着麻绳绕所有钉子一圈
这个麻绳最后也构成了点集的凸包
其中 第二种理解是我们一个经典算法 卷包裹法(Gift Wrapping)的思路
卷包裹算法从一个必然在凸包上的点开始向着一个方向依次选择最外侧的点
当回到最初的点时 所选出的点集就是所要求的凸包
这里还有两个问题不是很清楚:
1.怎么确定一个肯定在凸包上的点?
这个问题很好解决 取一个最左边的也就是横坐标最小的点
如果有多个这样的点 就取这些点里 纵坐标最小的
这样可以很好的处理共线的情况
2.如何确定下一个点(即最外侧的点)?
我们需要利用向量的叉积来解决这个问题
-----------------------------------------------------------------------
向量的叉积(Cross Product)原本是三维空间中的问题 在二维中也有巧妙的应用
http://mathworld.wolfram.com/CrossProduct.html
(下文中所有的叉积 若不作特殊说明 都是指二维中新定义的叉积
下文中所有的向量乘法 若不作特殊说明 都是指向量的叉积)
我们定义二维向量<x1,y1>和<x2,y2>的叉积为一个实数Cp=x1*y2-x2*y1
叉积有两条性质很常用 也很好用
1.叉积的一半是一个三角形的有向面积
这个公式可以避免面积计算的误差 如果点是整点 那么所有运算都是整数
2.向量的叉积的符号代表着向量旋转的方向
向量的叉积是不满足交换律的
向量A乘以向量B 如果为正则为A逆时针旋转向B 否则为顺时针
当然这里A转向B的角总是考虑一个小于180度以内的角 否则就会出错
-----------------------------------------------------------------------
有了向量 我们就可以选取一个最外侧的点了
比如现在我们卷包裹卷到J点我们要选取一个最外侧的点
当然比较红色的到角可以直接得到最外侧的点 不过不方便
我们可以考虑那个蓝色的到角
利用向量 我们可以比较哪个点"更外侧"
比如点K和点I 我们利用向量JK乘以向量JI得到一个数 这个数应该是负数 说明I比K更外侧
两个向量的比较具有传递性 所以我们可以像N个数里取最大的数一样取出最外侧的
遍历所有点 每个点都和现有最外侧的点比较 得到新的最外侧的点
至此两个问题都得以解决 我们可以写出满足一般要求的卷包裹算法了
不过还遗留有一个问题 就是处理共线的问题
有时候我们需要凸包边上的点也考虑到 有时候却需要去掉这些点
我们通常称在凸包顶点处的点为极点
如果我们只要求保留极点而去除在边上的点
我们只需在取外侧的点的时候 碰到共线的点取最远的
相反 如果我们要保留所有在边上的点我们只需要在共线的点中取最近的
这样整个卷包裹法终于完成了
给出完整的代码:
GiftWrapping
GiftWrapping
{$inline on}
const zero=1e-6; maxn=100000;
type point=record x,y:extended; end;
var p:array[1..maxn]of point;
ch:array[1..maxn]of longint;
temp,n,m,i,j,k:longint;
function sgn(x:extended):longint; inline;
begin
if abs(x)<zero then exit(0);
if x<0 then sgn:=-1 else sgn:=1;
end;
function cross(a,b,c:point):extended; 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):extended; inline;
begin
dist:=sqr(a.x-b.x)+sqr(a.y-b.y);
end;
function cmp(a,b,c:point):boolean; inline;
var temp:longint;
begin
temp:=sgn(cross(a,b,c));
if temp<>0 then exit(temp<0); {*B}
cmp:=dist(a,b)<dist(a,c); {*A}
end;
begin
assign(input,'Hull.in'); reset(input);
assign(output,'Hull1.out'); rewrite(output);
readln(n);
for i:=1 to n do
begin
readln(p[i].x,p[i].y);
if (j=0)or(p[i].x<p[j].x)or
(sgn(p[i].x-p[j].x)=0)and(p[i].y<p[j].y)
then j:=i;
end;
temp:=j;
while true do
begin
k:=0;
inc(m); ch[m]:=j;
for i:=1 to n do
if (i<>j)and((k=0)or
cmp(p[j],p[k],p[i]))
then k:=i;
if k=temp then break;
j:=k;
end;
for i:=1 to m do
writeln(p[ch[i]].x:0:2,' ',p[ch[i]].y:0:2);
close(input); close(output);
end.
*A Change Direction
*B Remove Colinear Points
还有两点要补充说明一下:
1.卷包裹算法的复杂度是O(NH)
N是全部的点数 H是最终在凸包上的点数
所以卷包裹算法很适合凸包上的点很少的时候 通常随机数据很快
但是构造出的凸包上的点很多的数据 这个算法就会很慢
比如所有的点都在一个圆周上
2.卷包裹算法输出的点是有序的
这也是对二维凸包算法的一个基本要求
通常只有保证有序才能进行进一步的计算
通过改变CMP函数可以改变上文中提到的共线(取/不取)以及这里的序(顺时针/逆时针)的问题
====================================================================
三.凸包的面积和周长
凸包的一个简单算法介绍完了
来看一个具体的问题 Poj 1113 http://poj.org/problem?id=1113
问题是求到一个简单多边形距离最少为L的点的轨迹的周长
我们可以先求凸包 然后得到计算公式
实际上 圆的部分可以拼接起来变成一个整圆
这样就好算多了 注意把共线的都去掉 不然误差会很大的
代码就不贴了 很简单的
下一篇介绍更快速的平面凸包算法
Graham Scan 和 QuickHull
强烈推荐QuickHull
第二篇：
一.卷包裹算法(Gift Wrapping Algorithm)的特性
前面提到过卷包裹算法的复杂度问题
由于卷包裹算法是两重循环实现的 因此很好分析它的复杂度
1 while true
do
2 begin
3 k:=0;
4 inc(m); ch[m]:=j;
5 for i:=1
to n do
6 if (i<>j)and((k=0)or
7 cmp(p[j],p[k],p[i]))
8 then k:=i;
9 if k=temp
then break;
10 j:=k;
11 end;
内部循环N次 外循环的次数决定于凸包上的点数H
所以卷包裹算法复杂度为O(HN) 这个复杂度很有特点
考虑一个随机点集的凸包上的点数往往很少 但是最坏情况下是O(N)级别的
比如所有的点都在一个圆上的时候
这就决定了卷包裹算法很适合随机的点集 但是如果是刻意构造的数据或是比较特殊的数据
就会达到O(N^2)的最坏复杂度 这是我们不愿意看到的
下面介绍最坏情况下复杂度更好的Graham扫描算法(Graham ScanAlgorithm)
====================================================================
二.Graham扫描算法(Graham ScanAlgorithm)
Graham扫描算法维护一个凸壳 通过不断在凸壳中加入新的点和去除影响凸性的点 最后形成凸包
The Graham scan isa method of computing the convex hull of a finite set of points in the planewith time complexity O(n log n). It is named after Ronald Graham, who publishedthe original algorithm in 1972.[1] The algorithm finds all vertices
of theconvex hull ordered along its boundary.
http://en.wikipedia.org/wiki/Graham_scan
算法主体由两部分组成 先是排序 后是扫描 分块讲解一下
----------------------------------------------------------------------------------------------------------
1.点集排序
为了得到加入新点的顺序 Graham扫描法的第一步是对点集排序
排序是对杂乱的点集进行了梳理 这也是这种算法能够得到更高效率的根本原因
排序的方法也有两种 极角坐标排序(极角序) 和 直角坐标排序(水平序)
前者好理解一些 但是在实现的时候 后者更方便
先说极角序 为了极角排序 我们先得得到一个参考点
一般的 我们取最左边(横坐标最小)的点作为参考点 如果有多个这样的点就取最下面的(纵坐标最小)
看这样一个例子 这是一个任意给出的平面点集:
参考点的定义:在横坐标最小的情况下取纵坐标最小的点
所以所有的点只能在这个黄色的半平面中 而且正上方为闭(可取得) 正下方为开(不可取)
这就决定了参考点的性质:点集中任意两点和参考点所成的到角为锐角
这样我们取得参考点 然后再考虑极角排序
极角排序以参考点为极角坐标系原点 各个点的极角为关键字
由于上面我们得到的参考点的性质 我们可以设所有点的极角均在(-90,90]之间
排序完成后应该是这样的:
比较极角我们仍然可以利用向量的叉积
叉积在这里已经介绍了 http://www.cnblogs.com/Booble/archive/2011/02/28/1967179.html
同样由于参考点的性质 所有向量之间的到角都是在180度以内 不会产生错误
----------------------------------------------------------------------------------------------------------
2.Graham的栈扫描
Graham的扫描是一个很优美的过程 用到的数据结构也很简单 仅仅是一个栈而已
核心的思想是按照排好的序 依次加入新点得到新的边
如果和上一条边成左转关系就压栈继续 如果右转就弹栈直到和栈顶两点的边成左转关系 压栈继续
实现的时候我们不用存边 只需要含顺序在栈里存点 相邻两点就是一条边
由于我们时时刻刻都保证栈内是一个凸壳 所以最后扫描完毕 就得到了一个凸包
下面还是继续上面的那个样例 演示一下栈扫描的过程
这样Graham扫描算法基本完成
复杂度是排序O(Nlog2N) 扫描O(N) {每个点仅仅出入栈一次}
合起来是一个O(Nlog2N)的算法 很优秀
----------------------------------------------------------------------------------------------------------
3.双重共线点难题
和卷包裹算法一样 我们同样还要考虑共线点问题
而且Graham扫描算法的共线问题更复杂 所以需要仔细考虑
i).排序时的共线问题
如果极角相同 我们应该怎么定先后呢?
我们得加上第二关键字距离 比如极角相同 距离参考点近的先
不过不管是近的先还是 远的先 开始和结束的两条边总是矛盾的
我们必须对其中一条特殊处理 除了结束边外距离近的先 结束边上距离远的先
这就是为什么极角排序不是很好实现的原因了 下面会介绍一下水平序
ii).扫描时的共线问题
这个和卷包裹算法的处理放法如出一辙
如果需要保留共线的点就在到角相同时取距离最近的
如果仅仅需要凸包极点就取距离最远的
----------------------------------------------------------------------------------------------------------
4.直角坐标排序(水平序)
直角坐标排序方法没有了极角排序的不足
以横坐标为第一关键字 纵坐标为第二关键字 排序点集
然后从第一个点开始 分别利用Graham扫描生成左链和右链
需要注意以下两点:
i).左链和右链的旋转方向是相反的
ii).注意生成第二条链要忽略第一条链上的点
由于水平序的CMP函数比较简单 代码也更短
还有一件有意思的事 Graham扫描算法是1972年提出的 卷包裹算法是1973年提出的
其实不奇怪 这两个算法本来也没有优劣所用的思想不同各自善于处理不同的情况
Graham扫描所用时间在点集随机时 还不如卷包裹算法快
正如第一节所说 卷包裹算法已经可以处理大部分情况 也是一个可取的算法
实际应用中点集更趋向于均匀 而不是集中在凸包上
----------------------------------------------------------------------------------------------------------
最后给一下比较难实现的极角排序代码
GrahamScan
{$inline on}
const zero=1e-6; maxn=100000;
type point=record x,y:extended; end;
var p,ch:array[1..maxn]of point;
i,j,n,m,temp:longint;
function sgn(x:extended):longint; inline;
begin
if abs(x)<zero then exit(0);
if x<0 then sgn:=-1 else sgn:=1;
end;
function cross(a,b,c,d:point):extended; inline;
begin
cross:=(b.x-a.x)*(d.y-c.y)-(d.x-c.x)*(b.y-a.y);
end;
function dist(a,b:point):extended; inline;
begin
dist:=sqr(a.x-b.x)+sqr(a.y-b.y);
end;
function cmp1(a,b,c:point):boolean; inline;
var temp:longint;
begin
temp:=sgn(cross(a,b,a,c));
if temp<>0 then exit(temp>0); {*B}
cmp1:=dist(a,b)<dist(a,c);
end;
function cmp2(a,b,c:point):boolean; inline;
var temp:longint;
begin
temp:=sgn(cross(a,b,b,c));
if temp<>0 then exit(temp<0); {*B}
cmp2:=dist(a,b)<dist(a,c); {*A}
end;
procedure swap(var a,b:point); inline;
var temp:point;
begin
temp:=a; a:=b; b:=temp;
end;
procedure sort(l,r:longint);
var i,j:longint;
x:point;
begin
i:=l; j:=r;
x:=p[(l+r)>>1];
repeat
while cmp1(p[1],p[i],x) do inc(i);
while cmp1(p[1],x,p[j]) do dec(j);
if not(i>j)
then begin
swap(p[i],p[j]);
inc(i); dec(j);
end;
until i>j;
if i<r then sort(i,r);
if l<j then sort(l,j);
end;
begin
assign(input,'Hull.in'); reset(input);
assign(output,'Hull2.out'); rewrite(output);
readln(n);
for i:=1 to n do
begin
readln(p[i].x,p[i].y);
if (p[i].x<p[1].x)or
(sgn(p[i].x-p[1].x)=0)and(p[i].y<p[1].y)
then swap(p[1],p[i]);
end;
sort(2,n);
while i>2 do
begin
if sgn(cross(p[1],p[i],p[1],p[i-1]))<>0
then break; dec(i);
end;
temp:=(i+n+1)>>1;
for j:=n downto temp do
swap(p[j],p[i+n-j]);
m:=1; ch[1]:=p[1];
inc(n); p
:=p[1];
for i:=2 to n do
begin
while m>1 do
begin
if not cmp2(ch[m-1],ch[m],p[i])
then break; dec(m);
end;
inc(m); ch[m]:=p[i];
end;
for i:=1 to m-1 do
writeln(ch[i].x:0:2,' ',ch[i].y:0:2);
close(input); close(output);
end.
*A Change Direction
*B Remove Colinear Points
====================================================================
三.快速凸包算法(Quickhull Algorithm)
对比Graham扫描算法和卷包裹算法
我们发现 Graham扫描算法在凸包上的点很密集时仍然适用
卷包裹算法在凸包上点集随机分布时是很高效的
那么有没有两个优点都具备的算法呢?
是有的! 快速凸包算法(Quickhull Algorithm)就是这样的一个算法
快速凸包算法是一个和快速排序(Quicksort Algorithm)神似的算法
尽管快速排序的最坏复杂度可以达到O(N^2)
但是有着极小的常数 实现方便 思路优美 绝大多数情况特别高效的快速排序 还是赢得了更多人的青睐
快速凸包算法也是这样 尽管可以构造一个数据使之达到O(N^2)的复杂度
但这需要刻意针对程序经过分析 才能做到 是实际应用中根本不会碰到的情况
在点集均匀分布时 快速凸包的复杂度更是达到了O(N) 是上面两种算法难以企及的
在绝大多数情况下 平均复杂度是O(Nlog2N) 也很高效
快速凸包继承了快速排序分治的思想 这是一个递归的过程
------------------------------------------------------------------------------------
伪代码如下:
1 void 快速凸包(P:点集 , S:向量
/*S.p,S.q:点)*/ ){
2 /* P 在 S
左侧*/
3 选取 P 中距离 S 最远的 点 Y ;
4 向量 A <- { S.p , Y }; 向量 B <- { Y , S.q } ;
5 点集 Q <- 在 P 中 且在 A 左侧的点 ;
6 点集 R <- 在 P 中 且在 B 左侧的点 ;
/* 划分 */
7 快速凸包 ( Q , A ) ; /*
分治 */
8 输出 (点 Y) ; /*
按中序输出 保证顺序*/
9 快速凸包 ( P , B ) ; /*
分治 */
10 }
初始化就选取 最左下和最右上的点 划分好 然后调用两次快速凸包函数分别求上下凸包
其中 选取和划分都需要用到向量的叉乘 注意方向
另外 存储点集可以用数组里的连续一段 传参的时候就传左右下标
划分的时候就是给数组交换元素 使新的点集 变成连续的一段
这里有一个很好的演示
http://www.cs.princeton.edu/courses/archive/spr10/cos226/demo/ah/QuickHull.html
还要补充说明一下
快速凸包在所有点都在圆周上的时候还是O(Nlog2N) 不过会比Graham扫描算法慢一些
可以说 这种数据是Graham扫描法的最好情况 一遍走完就行了
构造快速凸包的最坏情况就是使划分不均等 和构造快速排序最坏情况一样
贴以下我很辛苦写出来的代码吧 为了好看些 还心血来潮用cpp写了...
QuickHull
#include <iostream>
#include <math.h>
#define maxn 100000
#define zero 1e-12
#define sgn(x) (fabs(x)<zero?0:(x>0?1:-1))
#define cross(a,b,c) ((b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x))
#define cmp(a,b) (a.x<b.x || sgn(a.x-b.x)==0 && a.y<b.y)
using namespace std;
struct point{
double x,y;
}p[maxn];
double s[maxn];
void hull(int l,int r,point a,point b){
int x=l,i=l-1,j=r+1,k;
for (k=l;k<=r;k++){
double temp=sgn(s[x]-s[k]);
if (temp<0 || temp==0 && cmp(p[x],p[k])) x=k;
}
point y=p[x];
for (k=l;k<=r;k++){
s[++i]=cross(p[k],a,y);
if (sgn(s[i])>0) swap(p[i],p[k]); else i--;
}
for (k=r;k>=l;k--){
s[--j]=cross(p[k],y,b);
if (sgn(s[j])>0) swap(p[j],p[k]); else j++;
}
if (l<=i) hull(l,i,a,y);
printf("%.4lf %.4lf\n",y.x,y.y);
if (j<=r) hull(j,r,y,b);
}
int main(){
freopen("CH2D.in","r",stdin);
freopen("CH2D1.out","w",stdout);
int n,i,x=0;
scanf("%d",&n);
for (i=1;i<=n;i++){
scanf("%lf %lf",&p[i].x,&p[i].y);
if (x==0 || cmp(p[i],p[x])) x=i;
}
swap(p[1],p[x]);
printf("%.4lf %.4lf\n",p[1].x,p[1].y);
hull(2,n,p[1],p[1]);
return 0;
}
====================================================================
四.凸包算法复杂度下界
(引自<算法艺术与信息学竞赛>)
====================================================================
五.高效算法的应用
这里有一个平面最远点对的问题 可以利用凸包解决
Poj 2187 http://poj.org/problem?id=2187
由于最远点对必然在凸包上
我们先求凸包 然后枚举凸包上的点 不过这个复杂度最坏是O(N^2)的
不过凸包在很多情况下会改善问题的平均复杂度
凸包上的点通常很少 所以这个问题也可以过
篇幅关系 在下一节会介绍降低最坏复杂度的旋转卡壳算法
====================================================================
文中部分图片来源 http://westhoffswelt.de/blog/0040_quickhull_introduction_and_php_implementation.html http://en.wikipedia.org/wiki/Graham_scan
第三篇：
{
上一节介绍了凸包的高效算法
和一个最远点对的应用
这一段将更好的解决最远点对问题
}
(若不做特殊说明 下文讨论的问题均是在欧氏空间
若不做特殊说明 下文中距离均是指空间中欧氏距离)
==============================
一.简单枚举算法的不足
上一次介绍了一个基本的求平面最远点对的算法
即先求点集的凸包 然后枚举凸包上的点来求最远点集
这是利用了凸包上的点相比 点集中的点 一般是很少的 平均情况很好 并且我们也能AC这个问题
但是这是有局限性的 当凸包上的点达到O(N)的级别时 凸包的优化作用就不存在了
不过我们还要考虑到 凸包还起了对凸包上点集排序的作用
凸包有很多的优美的性质 我们可以加以利用 以得到更加高效的算法
旋转卡壳算法就是利用凸包特性的一类解决问题的方法
==============================
二.旋转卡壳算法
旋转卡(qiǎ)壳算法(Rotating CalipersAlgorithm):
是解决一些与凸包有关问题的有效算法 就像一对卡壳卡住凸包旋转而得名
Every time oneblade of the caliper lies flat against an edge of the polygon, it forms anantipodal pair with the point or edge touching the opposite blade. It turns outthat the complete "rotation" of the caliper around the polygondetects 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
不过不是很好理解 这里作一下说明
1 ch[m+1]:=ch[1]; j:=2;
2 for i:=1
to m do
3 begin
4 whilecross(ch[i],ch[j],ch[i+1])<cross(ch[i],ch[j+1],ch[i+1])
do
5 begininc(j);
if j>m then j:=1;
end;
6 writeln(ch[i].x,' ',ch[i].y,' ',ch[j].x,' ',ch[j].y);
7 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)
代码如下:
Maxd
{$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凸包
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- 计算几何之凸包_卷包裹算法
- 计算几何 Graham_scan凸包 pku 1113
- 算法模板——计算几何2（二维凸包——Andrew算法）
- 计算几何--求凸包模板--Graham算法--poj 1113
- [Poj 1113] 计算几何之凸包(一) {卷包裹算法}
- 计算几何---凸包问题（Graham/Andrew Scan ）
- 计算几何之凸包_卷包裹算法
- 二维计算几何模板--多边形/凸包
- 【Codeforces Round 335 (Div 2)E】【计算几何-凸包 线性规划 三分凸包上最优点】Freelancer's Dreams 二维属性 充最少的钱变得满足要求 [计算几何-凸包模
- 计算几何-经典算法-凸包
- <转> 最小凸包算法(Convex Hull)(1)-Graham扫描法 -计算几何-算法导论
- 凸包GiftWrapping GrahamScan 算法实现
- [计算几何]凸包的旋转卡壳算法
- 计算几何之凸包----Graham扫描法
- [Poj 2187] 计算几何之凸包(二) {更高效的算法}
- HDU 3692 Shade of Hallelujah Mountain （计算几何，三维空间点的旋转，二维凸包）
- 【计算几何】 Andrew凸包算法 + 旋转卡壳(以求点对最长距离为例) -- 以 POJ 2187 Beauty Contest 为例
- 【计算几何初步-凸包-Graham扫描法-极角序】【HDU1348】 WALL
- HDU 4449 Building Design 第37届ACM/ICPC 金华赛区H题 （计算几何，三维凸包+空间坐标旋转+二维凸包）
- 计算几何 ： 凸包学习笔记 --- Graham 扫描法