Logo_iSurf_org

Subscribe
to eAlerts

for Geom Site
Updates

CLICK HERE

 

Fast Convex Hull of a 2D Simple Polyline

by Dan Sunday

 

 

Algorithm 10 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(nlog-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. It is a very elegant and fast algorithm. It can be applied to the boundary of any connected 2D region, such as a geographic area (country, town, lake, etc).

How is this possible? Well, the general 2D hull algorithms first sort the vertex point set in O(nlog-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 15 about Tangents to and between Polygons, there are efficient algorithms for tangents to and between convex polygons. And, the solution for a polygon's convex hull is also the 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 was 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 that was difficult to implement. However, this was later simplified to a single stack method 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 the vertices sequentially once. Thus, it can be an on-line algorithm that reads a single input stream.
     
  3. Uses a double-ended queue (a deque) to store the 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(nlog-n) time because they 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 Hk–1 of the first k points OMEGA-k-1={P0,P1,...,Pk-1}, incrementally add the next point Pk, and then compute the next hull Hk. How does presorting facilitate this process? The answer is that at each stage, one knows that the next point Pk is exterior to the previous hull Hk–1, and thus one does not have to test for this. Otherwise, one would have to test Pk against all k edges of Hk–1, resulting in O(n2) such tests totaled over all stages of the algorithm. Instead, one immediately knows that that Pk is outside Hk–1, and can proceed to construct Hk by extending Hk–1. Thus, if presorting in O(nlog-n) time has been done, no such tests need to be made, and this results in a faster algorithm.

Additionally, most fast convex hull algorithms construct Hk from Pk and Hk–1 in a similar manner. Namely, they find the tangents from Pk to Hk–1 and use these as new edges for Hk. 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 Pk are at the stack top. So, one just needs to push Pk onto the stack after removing (pop off) any points at the stack top that get absorbed into the interior of the new hull Hk. This can be done by testing Pk against the line joining the top two points of the stack, and making sure that a left turn is being made for a counterclockwise (ccw) hull. The whole process is O(n) linear since each point is pushed onto the stack once and popped off at most once. For example, the Graham Scan algorithm radial sort computes this tangent on the stack as illustrated:

Pic_graham

The major advantage of presorting is to avoid expensive tests for whether each new incremental point Pk is included in the previous hull Hk–1.

However, if one knows that the initial ordered point set is a simple polyline, then inclusion can be tested by comparing Pk with only two edges of Hk–1. This results in O(n) inclusion edge tests in total over the life of the algorithm. Additionally, one can maintain Hk–1 as a stack that allows tangents from an exterior Pk 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 OMEGA-k={P0,P1,...,Pk} has its last vertex point Pk inside the hull Hk–1. The vertices of Hk–1 are a subset of the previous polyline vertices (but not necessarily in the same order). The geometry is illustrated in the diagram:

Pic_melkman1

The important observation to make is that the end of the polyline is now trapped inside a pocket formed by OMEGA-k, since it cannot cross over prior edges. Successive polyline points Pk+j (j > 0) are also trapped in this pocket until they escape to the exterior of Hk–1. However, a new vertex Pk+j can only escape across an edge of Hk–1 that contains Pk–1 since all other edges are blocked by OMEGA-k. Thus, each new Pk+j can be tested for inclusion with respect to at most two edges of Hk–1 until Pk+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 Hk–1 to form the extended hull Hk+j.

Pic_melkman2

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 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: OMEGA={P0,P1,...,Pn}. 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 Pk 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 PHI-cap={D_bot,...,D_top} where bot is the index at the bottom, and top is for the top of PHI-cap. The elements Di are vertices that form a polyline. When D_top=D_bot, then PHI-cap forms a polygon. In the Melkman hull algorithm, after processing vertex Pk, the deque PHI-cap_k satisfies:

  1. The polygon PHI-cap_k is the ccw convex hull Hk of the vertices OMEGA-k={P0,P1,...,Pk} already processed.
  2. D_top=D_bot is the most recent vertex processed that was added to PHI-cap_k.

If Pk is inside Hk–1, then PHI-cap_k=PHI-cap_k-1, and there is no associated processing. In this case, because the polyline is OMEGA-k simple, Pk is inside the subregion of Hk–1 bounded by the the vertices D_bot,D_bot+1,...,D_top-1,D_top, where "..." is the subpolyline of OMEGA-k that joins Dbot+1 and Dtop–1. And, Pk+j can only escape from this subregion by crossing over one of the edge segments D_bot D_bot+1 at the bottom or D_top-1 D_top at the top of PHI-cap_k. One can easily test for this by determining whether Pk+j is left of these edge segments. When it is left of both segments, this implies that Pk+j is still inside Hk–1; but as soon as Pk+j becomes right of either segment, then the vertex is exterior to Hk–1. Thus, at every incremental stage, inclusion of Pk in Hk–1 can be determined with only two isLeft() tests.

Pic_melkman3

When Pk is exterior to Hk–1, we must then change PHI-cap_k-1 to produce a new deque PHI-cap_k that satisfies the above two conditions. Clearly, by adding Pk to both the bottom and top of the deque, condition (2) is satisfied and Pk will be inside the new polygon defined by PHI-cap_k. However, other points already in PHI-cap_k-1 may get absorbed into the new hull Hk, and they need to be removed before Pk is added to the deque ends. As previously noted, vertices are removed from the two ends of PHI-cap_k-1 until the lines from Pk to the remaining deque endpoints form tangents to Hk–1. Again, this is easily determined by testing whether Pk 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 Hk, and we must remove that vertex from the deque. This continues until Pk is left of both the bottom and top edge segments of the deque. After that, we push Pk onto both ends, producing the required new deque PHI-cap_k.

Pic_melkman4

The speed of this algorithm is easily analyzed. Each vertex of OMEGA 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 DELTA-P0P1P2 is the final hull). Thus, the Melkman algorithm is a very efficient and elegant O(n) time and space algorithm.


Pseudo-Code: Melkman Algorithm

    Input: a simple polyline OMEGA with n vertices P[i]

    Put first 3 vertices onto deque D[] so that:
    a) 3rd vertex P[2] is at bottom and top of D[]
    b) on D[] they form a ccw triangle

    While there are more polyline vertices of OMEGA to process
    Get the next vertex P[i]
    {
        Note that:
        a) D[] is now the convex hull of already processed vertices
        b) D[bot] = D[top] = the last vertex added to D[]

        // Test if P[i] is inside D (as a polygon)
        If P[i] is left of both D[bot]D[bot+1] and D[top-1]D[top]
            Skip P[i] and Continue with the next vertex
        // Otherwise P[i] extends the hull and must be added to D[]

        // Get the tangent to the bottom
        While P[i] is right of D[bot]D[bot+1]
            Delete D[bot] from the bottom of D[]
        Insert P[i] at the bottom of D[]

        // Get the tangent to the top
        While P[i] is right of D[top-1]D[top]
            Pop D[top] from the top of D[]
        Push P[i] onto the top of D[]
    }

    Output: D[] = the ccw convex hull of OMEGA

 


Implementation

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

// Copyright 2001 softSurfer, 2012 Dan Sunday
// 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: 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;
}
 


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 2012 Dan Sunday, 2001 softSurfer