Graham Scan凸包算法

Posted 2020-09-20 czaoth

获得凸包的算法可以算是计算几何中最基础的算法之一了。寻找凸包的算法有很多种，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 1002
int 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;
}