Algorithm 15
[Home] [Overview] [History] [Algorithms] [Books] [Web Sites] [Gift Shop]

Convex Hull of a 2D Simple Polyline

by Dan Sunday

The algorithm about The Convex Hull of a Planar Point Set or Polygon showed how to compute the convex hull of any 2D point set or polygon with no restrictions.  The polygon could have been simple or not, connected or not.  It could even have been just a random set of segments or points.  The algorithms given, the "Graham Scan" and the "Andrew Chain", computed the hull in O(n log n) time.  This month we present a different algorithm that improves this efficiency to O(n) linear time for a connected simple polyline with no abnormal self-intersections.

How is this possible?  Well, the general 2D hull algorithms first sort the vertex point set in O(n log n) time, and then use a stack to compute the hull in O(n) time.  All the nonlinear overhead is in the initial sorting.  However, this month's algorithm uses the given sequential ordering of a simple polygon's edges along with a similar algorithm using a "deque" (a double-ended queue).  It turns out that the property of simplicity is enough to guarantee that the vertices on the deque are the convex hull.  In fact, the algorithm by [Melkman, 1987] that we present just assumes that the vertices form a simple polyline which is more general than earlier O(n) algorithms for simple polygons.

Before proceeding, we note that some polygon constructions can be performed more efficiently using the convex hull of the polygon.  For example, in Algorithm 14 about Tangents to and between Polygons, we saw there are efficient algorithms for tangents to and between convex polygons.  Further, the solution for a polygon's hull is also a solution for the polygon itself.  Thus, having a fast O(n) simple polygon hull algorithm also speeds up computing tangents for arbitrary simple polygons.

Background

The evolution of this algorithm has an interesting history that is given in detail at the outstanding web site of [Aloupis, 2000].  Aloupis describes both correct and incorrect algorithms and ideas discovered along the way, culminating in the [Melkman, 1987] algorithm, as well as attempts to make further improvements.

Initially, the O(n) improvement for simple polygons was proposed, and an algorithm implementation given, by [Sklansky, 1972].  Unfortunately, six years later, [Bykat, 1978] showed that the algorithm was flawed.  Nevertheless, the Sklansky conjecture was still valid, and his central idea was used by others to develop a correct algorithm.  The first valid algorithm was by [McCallum and Davis, 1979], but they had a complicated two stack method.  However, this was later simplified to a single stack by [Graham & Yao, 1983] and [Lee, 1983] (see also [Preparata & Shamos, 1985]).  At this point, correct algorithms worked for simple polygons, and did O(n) preprocessing to initialize a stack with a known (extreme) point on the hull.

Then, [Melkman, 1987] made a significant breakthrough by describing an algorithm that:

1. Works for a simple polyline (it does not have to be closed).

2. Does not do any preprocessing, and just processes vertices sequentially once.  Thus, this is an on-line algorithm that reads a single input stream.

3. Uses a double-ended queue (a deque) to store an incremental hull for the vertices already processed.

4. Greatly simplifies the logic of the algorithm.

It seems unlikely that the Melkman algorithm will be surpassed.

Simple Polyline Hull Algorithms

The Basic Incremental Strategy

Why should there be a faster O(n) convex hull algorithms for simple polylines and polygons?  To understand this, recall that most convex hull algorithms for point sets take O(n log n) time because initially sort the n points.  After that, they generally only require O(n) time.  So, one needs to ask why sorting is needed; that is, what does it accomplish?  Consider how other algorithms proceed after sorting is done.

Most 2D convex hull algorithms (see: The Convex Hull of a Planar Point Set) use a basic incremental strategy.  At the k-th stage, they have constructed the hull H_k-1 of the first k points {P₀, P₁, ..., P_k-1}, incrementally add the next point P_k, and then compute the next hull H_k.  How does presorting facilitate this process?  The answer is that at each stage, one knows that the next point P_k is exterior to the previous hull H_k-1, and thus one does not have to test for this.  Otherwise, one would have to test P_k against all k edges of H_k-1, resulting in O(n²) such tests totaled over all stages of the algorithm.  Instead, one immediately knows that that P_k is outside H_k-1, and can proceed to construct H_k by extending H_k-1.  Thus, if presorting in O(n log n) time has been done, no such tests need to be made, and this results in a faster algorithm. 

Additionally, most convex hull algorithms construct H_k from P_k and H_k-1 in a similar manner.  Namely, they find the tangents from P_k to H_k-1 and use these as new edges for H_k.  These tangents could be computed from scratch at each stage in O(log k) time (see: Tangents to Polygons), but many algorithms maintain a stack that simplifies the process.  Because of the presorting, the points on the stack are also sorted, and the ones closest to P_k are at the stack top.  So, one just needs to push P_k onto the stack after removing (popping) any points at the top that get absorbed into the interior of the new hull H_k.  This can be done by testing P_k against the line joining the top two points of the stack, and making sure that a left turn is being made (for a counterclockwise hull).  The whole process is O(n) linear since any given point is pushed once and popped at most once from the stack.  For example, the Graham Scan algorithm radial sort computes this tangent on the stack as illustrated:

Nevertheless, the major advantage of presorting is to avoid expensive tests for whether each new incremental point P_k is included in the previous hull H_k-1.  However, if one knows that the initial ordered point set is a simple polyline, then inclusion can be tested by comparing P_k with only one edge of H_k-1.  This results in only O(n) inclusion edge tests in total over the life of the algorithm.  Additionally, one can maintain H_k-1 as a stack that allows tangents from an exterior P_k to be computed by popping elements from the stack.  But, most importantly, the presorting is skipped, and the whole algorithm takes only O(n) time.

But, how does one simplify the inclusion test?  Suppose at the k-th stage that the simple polyline W_k = {P₀, P₁, ..., P_k} has the last vertex point P_k inside the hull H_k-1.  The vertices of H_k-1 are a subset of the previous polyline vertices (but not necessarily in the same order).  The geometry is illustrated in the diagram:

The important observation to make is that the end of the polyline is now trapped inside a pocket formed by W_k since it cannot cross over prior edges.  Successive polyline points P_k+j (j>0) are also trapped in this pocket until they escape to the exterior of H_k-1.  However, a new vertex P_k+j can only escape across an edge of H_k-1 that contains P_k-1 since all other edges are blocked by W_k.  Thus, each new P_k+j can be tested for inclusion with respect to at most two edges of H_k-1 until P_k+j is found to be outside one of those edges.  After that, tangents have to be found from this external vertex to the prior hull H_k-1 to form the next extended hull H_k+j. 

All simple polygon or polyline convex hull algorithms implement this strategy in one form or another.

The Melkman Algorithm

[Melkman, 1987] devised an ingenious method for organizing and implementing the operations to compute the hull of a simple polyline that we will now describe [Note: We have modified Melkman's algorithm to produce a convex hull with a counterclockwise (ccw) orientation, whereas the published algorithm gives a clockwise hull.  We have also changed his notation to match ours].

The strategy of the Melkman algorithm is straightforward.  It sequentially processes each of the polyline vertices in order.  Let the input polyline be given by the ordered vertex set: W = {P₀, P₁, ..., P_n}.  At each stage, the algorithm determines and stores (on a double-ended queue) those vertices that form the ordered hull for all polyline vertices considered so far.  Then, the next vertex P_k is considered.  It satisfies one of two conditions: either (1) it is inside the currently constructed hull, and can be ignored; or (2) it is outside the current hull, and becomes a new hull vertex extending the old hull.  However, in case (2), vertices that are on the list for the old hull, may become interior to the new hull, and need to be discarded before adding the new vertex to the new list.

The double-ended queue (called a "deque") has both a top and a bottom.  At both ends of the deque, elements can be either added or removed.  At the top, we say an element is pushed or popped; while at the bottom, we say an element is inserted or deleted.  The deque is given by an ordered list D = { d_bot, ..., d_top }where bot is the index at the bottom, and top is for the top of D.  The elements d_i are vertices that form a polyline.  When d_top = d_bot, then D forms a polygon.  In the Melkman hull algorithm, after processing vertex P_k, the deque D_k satisfies:

1. The polygon D_k is the ccw convex hull H_k of the vertices W_k = {P₀, ..., P_k} already processed.

2. d_top = d_bot is the most recent vertex processed that was added to D_k.

If P_k is inside H_k−1, then D_k = D_k−1, and there is no associated processing.  In this case, P_k is inside the subregion of H_k−1 bounded by the the vertices: d_bot, d_bot+1,..., d_top−1, d_top where "..." is the subpolyline of W_k that joins d_bot+1 and d_top−1.  P_k+j can only escape from this subregion  by crossing over one of the edge segments d_botd_bot+1 or d_top−1d_top at the bottom or the top of D_k−1.  One can easily test this by determining whether P_k+j is left of these edge segments.  When it is left of both segments, this implies that P_k+j is still inside H_k−1; but as soon as P_k+j becomes right of either segment, then the vertex is exterior to H_k−1.  Thus, at every incremental stage, inclusion of  P_k in H_k−1 can be determined with just two isLeft() tests.

When P_k is exterior to H_k−1, we must then change D_k−1 to produce a new deque D_k that satisfies the above two conditions.  Clearly, by adding P_k to both the bottom and top of the deque, condition (2) is satisfied and P_k will be inside the new polygon defined by D_k.  However, other points already in D_k−1 may get absorbed into the new hull H_k, and they need to be removed before P_k is added to the deque ends.  As previously noted, vertices are removed from the two ends of D_k−1 until the lines from P_k to the remaining deque endpoints form tangents to H_k−1.  Again, this is easily determined by testing whether P_k is left of the edge segments at the bottom and top of the deque.  If it is not to the left of a segment, then the current vertex at the end of the deque would be inside H_k, and we must remove that vertex from the deque.  This continues until P_k is left of both the bottom and top edge segments of the deque.  After that, we add P_k to both ends, producing the required new deque D_k.

The speed of this algorithm is easily analyzed.  Each vertex of W can get put on the deque at most twice (once at each end).  Thereafter, elements on the deque can be removed at most once.  Each of these events, to potentially add or remove a vertex to/from the deque, is associated with exactly one (constant) isLeft() test.  Thus, the worst case behavior of the Melkman algorithm is bounded by 3n tests and queue operations.  The best case behavior would have 2n tests and only 4 queue operations (when the initial triangle P₀P₁P₂ is the final hull).  Thus, the Melkman algorithm is an efficient and beautiful O(n) time and space algorithm.

Pseudo-Code: Melkman Algorithm

Implementation

Here is a "C++" implementation of this algorithm.

// Copyright 2002, softSurfer (www.softsurfer.com)
// This code may be freely used 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(): 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: the January 2001 Algorithm 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():
//    Input:  V[] = polyline array of 2D vertex points
//            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* V, int n, Point* H )
{
    // initialize a deque D[] from bottom to top so that the
    // 1st three vertices of V[] are a counterclockwise 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] = V[2];       // 3rd vertex is at both bot and top
    if (isLeft(V[0], V[1], V[2]) > 0) {
        D[bot+1] = V[0];
        D[bot+2] = V[1];          // ccw vertices are: 2,0,1,2
    }
    else {
        D[bot+1] = V[1];
        D[bot+2] = V[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], V[i]) > 0) &&
            (isLeft(D[top-1], D[top], V[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], V[i]) <= 0)
            ++bot;                // remove bot of deque
        D[--bot] = V[i];          // insert V[i] at bot of deque

        // get the leftmost tangent at the deque top
        while (isLeft(D[top-1], D[top], V[i]) <= 0)
            --top;                // pop top of deque
        D[++top] = V[i];          // push V[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;
}

References

Greg Aloupis, A History of Linear-time Convex Hull Algorithms for Simple Polygons, McGill Univ website (2000)

A. Bykat, "Convex hull of a finite set of points in two dimensions", Info. Proc. Letters 7, 296-298 (1978)

Ronald Graham & F. Yao, "Finding the convex hull of a simple polygon", J. Algorithms 4(4), 324-33 (1983)

D.T. Lee, "On finding the convex hull of a simple polygon", Int'l J. Comp. & Info. Sci. 12(2), 87-98 (1983)

D. McCallum and D. Davis, "A linear algorithm for finding the convex hull of a simple polygon", Info. Proc. Letters 9, 201-206 (1979)

A. Melkman, "On-line construction of the convex hull of a simple polygon", Info. Proc. Letters 25, 11-12 (1987)

Franco Preparata & Michael Shamos, Computational Geometry: An Introduction, Section 4.1.4 "Convex hull of a simple polygon" (1985)

Sklansky J., "Measuring Concavity on a Rectangular Mosaic", IEEE Transactions on Computing 21, p-1355 (1972).

Copyright © 2001-2006 softSurfer.  All rights reserved.
Email comments and suggestions to feedback@softsurfer.com