首页 > 代码库 > Graham Scan凸包算法
Graham Scan凸包算法
获得凸包的算法可以算是计算几何中最基础的算法之一了。寻找凸包的算法有很多种,Graham Scan算法是一种十分简单高效的二维凸包算法,能够在O(nlogn)的时间内找到凸包。
首先介绍一下二维向量的叉积(这里和真正的叉积还是不同的):对于二维向量a=(x1,y2)和b=(x2,y2),a×b定义为x1*y2-y1*x2。而它的几何意义就是|a||b|sin<a,b>。如果a与b夹角小于180度(逆时针),那么这个值就是正值,大于180度就是负值。需要注意的是,左乘和右乘是不同的。如图所示:
Graham Scan算法的做法是先定下一个起点,一般是最左边的点和最右边的点,然后一个个点扫过去,如果新加入的点和之前已经找到的点所构成的“壳”凸性没有变化,就继续扫,否则就把已经找到的最后一个点删去,再比较凸性,直到凸性不发生变化。分别扫描上下两个“壳”,合并在一起,凸包就找到了。这么说很抽象,我们看图来解释:
我们找下“壳”,上下其实是一样的。首先加入两个点A和C:
然后插入第三个点G,并计算AC×CG的叉积,却发现叉积小于0,也就是说逆时针方向上∠ACG大于180度,于是删去C点,加入G点:
然后就是依照这个步骤便能加入D点。在AD上方是以D为起点。就能够找到AGD和DFEA两个凸壳。合并就得到了凸包。
关于扫描的顺序,有坐标序和极角序两种。坐标序是比较两个点的x坐标,如果小的先被扫描(扫描上凸壳的时候反过来);如果两个点x坐标相同,那么就比较y坐标,小的先被扫描(扫描上凸壳的时候也是反过来)。极角序使用arctan2函数的返回值进行比较,我没写过所以也不是很清楚。
程序可以写得很精简,以下是我用C++写得凸包程序
/*d[]是一个Point的数组,Point有两个两个属性x和y,同时支持减法操作和det(叉积)。convex数组保存被选中的凸包的点的编号,cTotal是凸包中点的个数*/bool cmpPoint(const Point &a, const Point &b) //比较坐标序所用的比较函数{ if (a.x!=b.x) return a.x<b.x; return a.y<b.y;}void get_convex_hull(){ sort(d,d+N,cmpPoint); int Total=0,tmp; for (int i=0;i<N;++i) //扫描下凸壳 { while ( (Total>1) && ((d[convex[Total-1]]-d[convex[Total-2]]).det( //获得凸包中最后两个点的向量 d[i]-d[convex[Total-1]])<=0) ) Total--; //获得准备插入的点和凸包中最后一点的向量,计算叉积 convex[Total++]=i; } tmp=Total; for (int i=N-2;i>=0;--i) //扫描上凸壳 { while ( (Total>tmp) && ((d[convex[Total-1]]-d[convex[Total-2]]).det( d[i]-d[convex[Total-1]])<=0) ) Total--; convex[Total++]=i; } cTotal=Total;}
我们来看一道题:POJ1113 Wall,题意是给一些点,找一个闭合曲线C,使C能包住所有的点,并且给定的点到C的距离最小为L,问C的周长。稍微画一画就知道这个C的周长是这些点所构成的凸包的周长加上以L为半径的圆的周长。于是求一个凸包再加上2πL就可以了。我的程序如下:
#include <cstdio>#include <cstring>#include <cstdlib>#include <algorithm>#include <cmath>using std::sort;#define MAXN 1002int N,L;double sqr(double a){ return a*a; }struct Point{ double x,y; inline Point operator- (const Point &t) { Point ret; ret.x=x-t.x; ret.y=y-t.y; return ret; } inline Point operator+ (const Point &t) { Point ret; ret.x=x+t.x; ret.y=y+t.y; return ret; } inline int det(const Point &t) { return x*t.y-t.x*y; } inline double dist(Point &t) { return sqrt(sqr(x-t.x)+sqr(y-t.y)); }}d[MAXN];bool cmpPoint(const Point &a, const Point &b){ if (a.x!=b.x) return a.x<b.x; return a.y<b.y;}int convex[MAXN],cTotal;void get_convex_hull(){ sort(d,d+N,cmpPoint); int Total=0,tmp; for (int i=0;i<N;++i) { while ( (Total>1) && ((d[convex[Total-1]]-d[convex[Total-2]]).det( d[i]-d[convex[Total-1]])<=0) ) Total--; convex[Total++]=i; } tmp=Total; for (int i=N-2;i>=0;--i) { while ( (Total>tmp) && ((d[convex[Total-1]]-d[convex[Total-2]]).det( d[i]-d[convex[Total-1]])<=0) ) Total--; convex[Total++]=i; } cTotal=Total;}int main(){ scanf("%d%d",&N,&L); for (int i=0;i<N;++i) { scanf("%lf%lf",&d[i].x,&d[i].y); } get_convex_hull(); double Ans=0; for (int i=0;i<cTotal-1;++i) { Ans+=d[convex[i]].dist(d[convex[i+1]]); } Ans+=d[convex[0]].dist(d[convex[cTotal-1]]); Ans+=3.1415926*2*L; printf("%.0lf\n",Ans); return 0;}
Graham Scan凸包算法