多边形快速凸包算法

前言

平面点集的凸包算法一文介绍了如何计算平面点集或者任意多边形的凸包。对于随机的平面点集，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加入队列。

首先给双端队列两个标签:bot和top，在这中间的是当前凸包结果。在队列两端，都可以增加或删除元素。在顶部（top之上），我们称为push / pop；在底部（bot之下），我们将增删元素的操作称为insert / delete。不妨将队列记为，是原多段线中的点。当就形成了一个多边形。在Melkman算法中，处理顶点Pk后，满足：

1. 是多段线的逆时针方向的凸包；
2. ，是最近添加到中的点。

对于case（2）,我们需要改变，更新队列。在将Pk添加到队列两端之前，需要先将在新凸包内部的点删除。在队列的首尾，通过测试Pk是否在顶部的边的左侧，就可以判断此时top\bot处的点是否需要删除。继续这个检查，直到Pk在队列两端的边的左侧。最后，我们将Pk添加到队列两端。过程如下图：

根据上述过程，很容易分析算法的时间复杂度。每个顶点最多添加到队列中两次（top和bot各一次），队列中的点最多被移除一次，每添加/移除一个顶点，最多需要一次常数量级的isLeft判断。Melkman algorithm最多需要3n次isLeft测试和3n次队列操作。最佳性能是，2n次测试和4次队列操作（当最初的3个点构成最终的凸包结果时）。

因此，Melkman算法非常高效，时间复杂度和空间复杂度都是O(n).

算法伪代码如下：

Input: 有n个顶点的简单多段线P[i]

将初始3个点加入队列 D[] ，使得:
    a) P[2] 在 D[]的top和bot处
    b) 在D[]中，P0、P1、P2形成一个逆时针的三角形

依次处理i=2之后的每一个点，对于P[[i]，检查P[i]是否在D的内部：

     if  P[i]  在D[bot]D[bot+1] 和 D[top-1]D[top]的左侧 then

           跳过P[i],接着处理下一个点；

     while P[i] is right of D[bot]D[bot+1] do

          Delete D[bot] from the bottom of D[];

     Insert P[i] at the bottom of D[];

     while P[i] is right of D[top-1]D[top] do

          Pop D[top] from the top of D[]；

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