Thứ Bảy, 26 tháng 9, 2015

The Convex Hull of a Planar Point Set

Computing a convex hull (or just "hull") is one of the first sophisticated geometry algorithms, and there are many variations of it. The most common form of this algorithm involves determining the smallest convex set (called the "convex hull") containing a discrete set of points. This algorithm also applies to a polygon, or just any set of line segments, whose hull is the same as the hull of its vertex point set. There are numerous applications for convex hulls: collision avoidance, hidden object determination, and shape analysis to name a few. And they are a minimal linear bounding container.
The most popular hull algorithms are the "Graham scan" algorithm [Graham, 1972] and the "divide-and-conquer" algorithm [Preparata & Hong, 1977]. Implementations of both these algorithms are readily available (see [O'Rourke, 1998]). Both are O(nlog-n) time algorithms, but the Graham has a low runtime constant in 2D and runs very fast there. However, the Graham algorithm does not generalize to 3D and higher dimensions whereas the divide-and-conquer algorithm has a natural extension. We do not consider 3D algorithms here (see [O'Rourke, 1998] for more information).
Here is a list of some well-known 2D hull algorithms. Let n = # points in the input set, and h = # vertices on the output hull. Note that h.le.n, so nh.le.n2. The list is ordered by date of first publication.
Algorithm
Speed
Discovered By
Brute Force
O(n4)
[Anon, the dark ages]
Gift Wrapping
O(nh)
[Chand & Kapur, 1970]
Graham Scan
O(nlog-n)
[Graham, 1972]
Jarvis March
O(nh)
[Jarvis, 1973]
QuickHull
O(nh)
[Eddy, 1977], [Bykat, 1978]
Divide-and-Conquer
O(nlog-n)
[Preparata & Hong, 1977]
Monotone Chain
O(nlog-n)
[Andrew, 1979]
Incremental
O(nlog-n)
[Kallay, 1984]
Marriage-before-Conquest
O(nlog-h)
[Kirkpatrick & Seidel, 1986]


Convex Hulls

The convex hull of a geometric object (such as a point set or a polygon) is the smallest convex set containing that object. There are many equivalent definitions for a convex set S. The most basic of these is:
Def 1. A set S is convex if whenever two points P and Q are inside S, then the whole line segment PQ is also in S.
Pic_convex1
But this definition does not readily lead to algorithms for constructing convex sets. A more useful definition states:
Def 2. A set S is convex if it is exactly equal to the intersection of all the half planes containing it.
It can be shown that these two definitions are equivalent. However, the second one gives us a better computational handle, especially when the set S is the intersection of a finite number of half planes. In this case, the boundary of S is polygon in 2D, and polyhedron in 3D, with which it can be identified.
Def 3. The convex hull of a finite point set S = {P} is the smallest 2D convex polygon OMEGA (or polyhedron in 3D) that contains S. That is, there is no other convex polygon (or polyhedron) LAMBDA with S.in.LAMBDA.in.OMEGA.
Also, this convex hull has the smallest area and the smallest perimeter of all convex polygons that contain S.


2D Hull Algorithms

For this algorithm we will cover two similar fast 2D hull algorithms: the Graham scan, and Andrew's Monotone Chain scan. They both use a similar idea, and are implemented as a stack. In practice, they are both very fast, but Andrew's algorithm will execute slightly faster since its sort comparisons and rejection tests are more efficient. An implementation of Andrew's algorithm is given below in our chainHull_2D() routine.


Andrew's Monotone Chain Algorithm

[Andrew, 1979] discovered an alternative to the Graham scan that uses a linear lexographic sort of the point set by the x and y-coordinates. This is an advantage if this ordering is already known for a set, which is sometimes the case. But even if sorting is required, this is a faster sort than the angular Graham-scan sort with its more complicated comparison function. The "Monotone Chain" algorithm computes the upper and lower hulls of a monotone chain of points, which is why we refer to it as the "Monotone Chain" algorithm. Like the Graham scan, it runs in O(nlog-n) time due to the sort time. After that, it only takes O(n) time to compute the hull. This algorithm also uses a stack in a manner very similar to Graham's algorithm.
First the algorithm sorts the point set S={P0,P1,...,Pn-1} by increasing x and then y coordinate values. Let the minimum and maximum x-coordinates be xmin and xmax. Clearly, P0.x=x-min, but there may be other points with this minimum x-coordinate. Let Pminmin be a point in S with P.x=x-min first and then min y among all those points. Also, let Pminmax be the point with P.x=x-min first and then max y second. Note that Pminmin=Pminmax when there is a unique x-minimum point. Similarly define Pmaxmin and Pmaxmax as the points with P.x=x-max first, and then y min or max second. Again note that Pmaxmin=Pmaxmax when there is a uniquex-maximum point. Next, join the lower two points, Pminmin and Pmaxmin to define a lower line L-min. Also, join the upper two points, Pminmax and Pmaxmax to define an upper line L-max. These points and lines are shown in the following example diagram.
Pic_chain
The algorithm now proceeds to construct a lower convex vertex chain OMEGA-min below L-min and joining the two lower points Pminmin and Pmaxmin; and also an upper convex vertex chain OMEGA-max above L-max and joining the two upper points Pmaxmax and Pminmax. Then the convex hull OMEGA of S is constructed by joining OMEGA-min and OMEGA-max together.
The lower or upper convex chain is constructed using a stack algorithm almost identical to the one used for the Graham scan. For the lower chain, start with Pminmin on the stack. Then process the points of S in sequence. For OMEGA-min, only consider points strictly below the lower line L-min. Suppose that at any stage, the points on the stack are the convex hull of points below L-min that have already been processed. Now consider the next pointPk that is below L-min. If the stack contains only the one point Pminmin then put Pk onto the stack and proceed to the next stage. Otherwise, determine whether Pk is strictly left of the line between the top two points on the stack. If it is, put Pk onto the stack and proceed. If it is not, pop the top point off the stack, and test Pkagainst the stack again. Continue until Pk gets pushed onto the stack. After this stage, the stack again contains the vertices of the lower hull for the points already considered. The geometric rationale is exactly the same as for the Graham scan. After all points have been processed, push Pmaxmin onto the stack to complete the lower convex chain.
The upper convex chain OMEGA-max is constructed in an analogous manner. But, process S in decreasing order {Pn-1,Pn-2,...,P0}, starting at Pmaxmax, and only considering points above L-max. Once the two hull chains have been found, it is easy to join them together (but be careful to avoid duplicating the endpoints).

Pseudo-Code: Andrew's Monotone Chain Algorithm
    Input: a set S = {P = (P.x,P.y)} of N points

    Sort S by increasing x and then y-coordinate.
    Let P[] be the sorted array of N points.

    Get the points with 1st x min or max and 2nd y min or max
        minmin = index of P with min x first and min y second
        minmax = index of P with min x first and max y second
        maxmin = index of P with max x first and min y second
        maxmax = index of P with max x first and max y second

    Compute the lower hull stack as follows:
    (1) Let L_min be the lower line joining P[minmin] with  P[maxmin].
    (2) Push P[minmin] onto the stack.
    (3) for i = minmax+1 to maxmin-1 (the points between xmin and xmax)
        {
            if (P[i] is above or on L_min)
                 Ignore it and continue.
            while (there are at least 2 points on the stack)
            {
                 Let PT1 = the top point on the stack.
                 Let PT2 = the second point on the stack.
                 if (P[i] is strictly left of the line from PT2 to PT1)
                     break out of this while loop.
                 Pop the top point PT1 off the stack.
            }
            Push P[i] onto the stack.
        }
    (4) Push P[maxmin] onto the stack.

    Similarly, compute the upper hull stack.

    Let OMEGA = the join of the lower and upper hulls.

    Output: OMEGA = the convex hull of S.
 
Implementation
Here is a "C++" implementation of the Chain Hull algorithm.
// Copyright 2001 softSurfer, 2012 Dan Sunday
// This code may be freely used, distributed and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
 
// Assume that a class is already given for the object:
//    Point with coordinates {float x, y;}
//===================================================================
 
// isLeft(): tests 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);
}
//===================================================================

// chainHull_2D(): Andrew's monotone chain 2D convex hull algorithm
//     Input:  P[] = an array of 2D points
//                  presorted by increasing x and y-coordinates
//             n =  the number of points in P[]
//     Output: H[] = an array of the convex hull vertices (max is n)
//     Return: the number of points in H[]
int
chainHull_2D( Point* P, int n, Point* H )
{
    // the output array H[] will be used as the stack
    int    bot=0, top=(-1);   // indices for bottom and top of the stack
    int    i;                 // array scan index

    // Get the indices of points with min x-coord and min|max y-coord
    int minmin = 0, minmax;
    float xmin = P[0].x;
    for (i=1; i<n; i++)
        if (P[i].x != xmin) break;
    minmax = i-1;
    if (minmax == n-1) {       // degenerate case: all x-coords == xmin
        H[++top] = P[minmin];
        if (P[minmax].y != P[minmin].y) // a  nontrivial segment
            H[++top] =  P[minmax];
        H[++top] = P[minmin];          // add polygon endpoint
        return top+1;
    }

    // Get the indices of points with max x-coord and min|max y-coord
    int maxmin, maxmax = n-1;
    float xmax = P[n-1].x;
    for (i=n-2; i>=0; i--)
        if (P[i].x != xmax) break;
    maxmin = i+1;

    // Compute the lower hull on the stack H
    H[++top] = P[minmin];    // push  minmin point onto stack
    i = minmax;
    while (++i <= maxmin)
    {
        // the lower line joins P[minmin]  with P[maxmin]
        if (isLeft( P[minmin], P[maxmin], P[i])  >= 0 && i < maxmin)
            continue;  // ignore P[i] above or on the lower line

        while (top > 0)  
// there are at least 2 points on the stack
        {
            // test if  P[i] is left of the line at the stack top
            if (isLeft(  H[top-1], H[top], P[i]) > 0)
                 break;         // P[i] is a new hull  vertex
            else
                 top--;         // pop top point off  stack
        }
        H[++top] = P[i];        // push P[i] onto stack
    }

    // Next, compute the upper hull on the stack H above  the bottom hull
    if (maxmax != maxmin)      // if  distinct xmax points
         H[++top] = P[maxmax];  // push maxmax point onto stack
    bot = top;
// the bottom point of the upper hull stack
    i = maxmin;
    while (--i >= minmax)
    {
        // the upper line joins P[maxmax]  with P[minmax]
        if (isLeft( P[maxmax],P[minmax],P[i])>= 0 && i > minmax)
            continue;  // ignore P[i] below or on the upper line

        while (top > bot) // at least 2 points on the upper stack
        {
            // test if  P[i] is left of the line at the stack top
            if (isLeft(  H[top-1], H[top], P[i]) > 0)
                 break;     // P[i] is a new hull  vertex
            else
                 top--;     // pop top point off  stack
        }
        H[++top] = P[i];   // push P[i] onto stack
    }
    if (minmax != minmin)
        H[++top] = P[minmin];
// push  joining endpoint onto stack
    return top+1;
}

                                                        From: http://geomalgorithms.com/