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 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;
}