多边形快速凸包算法

前言
平面点集的凸包算法一文介绍了如何计算平面点集或者任意多边形的凸包。对于随机的平面点集，Graham scan和Andraw’s 单调链算法已经是最快的算法了。但是对于没有自相交的封闭的简单多边形，存在线性复杂度的算法。下面介绍这一优雅高效的算法。

一般的2D凸包算法，首先将点进行排序（时间复杂度）,然后利用栈操作在O(n)的时间复杂度内计算凸包。初始的排序决定了最终的时间复杂度。但是本文介绍的算法使用一个双端队列来进行操作，避免了排序。由于限定了多边形的简单性（平面单连通域），可以证明队列中的点构成凸包。该算法是由Melkman在1987年提出的。

一些多边形的特征算法可以通过其凸包来高效地求解，其凸包的解就是原来多边形的解。因此，对于简单多边形有一个快速凸包算法的话，可以加速相应算法的计算。例如多边形的直径、切线等算法。

背景
早在1972年，Sklansky就提出了一个O(n)时间复杂度的算法，并给出了实现。不幸的是，6年后Bykat证明他的算法是错误的。基于Sklansky的想法，Graham & Yao在1983年修正了这个算法，给出了一个使用栈操作的正确版本，但是算法的实现十分复杂。

最终，Melkman在1987年给出了一个简洁漂亮的O(n)算法：

1、适用于简单多段线（不自交）；
2、不需要任何预处理，直接按顶点顺序处理；
3、使用双端队列存储凸包结果；
4、算法的逻辑非常简单。
Melkman算法似乎不太可能被超越了。

简单多边形凸包算法
Melkman’s Algorithm
Melkman, 1987年设计了一种巧妙的方法来实现计算简单多段线的凸包，下面将对其进行详细描述。

Melkman Algorithm的策略非常直接，按原始顺序依次处理多段线上的每个点。假定输入多段线为S={P0,P1,…,Pn}。在每一步，算法将当前为止处理过的所有顶点形成的凸包存储在一个双端队列中。

接下来，考虑下一个顶点Pk。Pk有两种情况：（1）Pk在当前凸包内；（2）Pk在当前凸包外，并且形成一个新的凸包顶点。在case （2）中，原来凸包中的点可能变为在新凸包的内部，需要被丢弃，然后再将Pk加入队列。


Push P[i] onto the top of D[].

Output: D[]就是最终的凸包结果。

C++实现

// Assume that a class is already given for the object:
//    Point with coordinates {float x, y;}
//===================================================================
 
 
// isLeft(): test if a point is Left|On|Right of an infinite line.
//    Input:  three points P0, P1, and P2
//    Return: >0 for P2 left of the line through P0 and P1
//            =0 for P2 on the line
//            <0 for P2 right of the line
//    See: Algorithm 1 on Area of Triangles
inline float isLeft( Point P0, Point P1, Point P2 )
{
    return (P1.x - P0.x)*(P2.y - P0.y) - (P2.x - P0.x)*(P1.y - P0.y);
}
 
 
 
// simpleHull_2D(): Melkman's 2D simple polyline O(n) convex hull algorithm
//    Input:  P[] = array of 2D vertex points for a simple polyline
//            n   = the number of points in V[]
//    Output: H[] = output convex hull array of vertices (max is n)
//    Return: h   = the number of points in H[]
int simpleHull_2D( Point* P, int n, Point* H )
{
    // initialize a deque D[] from bottom to top so that the
    // 1st three vertices of P[] are a ccw triangle
    Point* D = new Point[2*n+1];
    int bot = n-2, top = bot+3;    // initial bottom and top deque indices
    D[bot] = D[top] = P[2];        // 3rd vertex is at both bot and top
    if (isLeft(P[0], P[1], P[2]) > 0) {
        D[bot+1] = P[0];
        D[bot+2] = P[1];           // ccw vertices are: 2,0,1,2
    }
    else {
        D[bot+1] = P[1];
        D[bot+2] = P[0];           // ccw vertices are: 2,1,0,2
    }
 
    // compute the hull on the deque D[]
    for (int i=3; i < n; i++) {   // process the rest of vertices
        // test if next vertex is inside the deque hull
        if ((isLeft(D[bot], D[bot+1], P[i]) > 0) &&
            (isLeft(D[top-1], D[top], P[i]) > 0) )
                 continue;         // skip an interior vertex
 
        // incrementally add an exterior vertex to the deque hull
        // get the rightmost tangent at the deque bot
        while (isLeft(D[bot], D[bot+1], P[i]) <= 0)
            ++bot;                 // remove bot of deque
        D[--bot] = P[i];           // insert P[i] at bot of deque
 
        // get the leftmost tangent at the deque top
        while (isLeft(D[top-1], D[top], P[i]) <= 0)
            --top;                 // pop top of deque
        D[++top] = P[i];           // push P[i] onto top of deque
    }
 
    // transcribe deque D[] to the output hull array H[]
    int h;        // hull vertex counter
    for (h=0; h <= (top-bot); h++)
        H[h] = D[bot + h];
 
    delete D;
    return h-1;
}