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 time. This month we present a different algorithm that improves this efficiency to linear time for a connected simple polyline with no abnormal selfintersections. 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 time, and then use a stack to compute the hull in 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 doubleended 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 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 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 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 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:
 Works for a simple polyline (it does not have to be closed).
 Does not do any preprocessing, and just processes the vertices sequentially once. Thus, it can be an online algorithm that reads a single input stream.
 Uses a doubleended queue (a deque) to store the incremental hull for the vertices already processed.
 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 convex hull algorithms for simple polylines and polygons? To understand this, recall that most convex hull algorithms for point sets take time because they initially sort the n points . After that, they generally only require 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 kth stage, they have constructed the hull H_{k–1} of the first k points , 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 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 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 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 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 (pop off) any points at the stack 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 (ccw) hull. The whole process is 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:
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 two edges of H_{k–1}. This results in 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 time.
But, how does one simplify the inclusion test? Suppose at the kth stage that the simple polyline has its 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 , 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 . 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 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 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: . At each stage, the algorithm determines and stores (on a doubleended 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 doubleended 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 where bot is the index at the bottom, and top is for the top of . The elements D_{i} are vertices that form a polyline. When , then forms a polygon. In the Melkman hull algorithm, after processing vertex P_{k}, the deque satisfies:
 The polygon is the ccw convex hull H_{k} of the vertices already processed.
 is the most recent vertex processed that was added to .
If P_{k} is inside H_{k–1}, then , and there is no associated processing. In this case, because the polyline is simple, P_{k} is inside the subregion of H_{k–1} bounded by the the vertices , where "..." is the subpolyline of that joins D_{bot+1} and D_{top–1}. And, P_{k+j} can only escape from this subregion by crossing over one of the edge segments at the bottom or at the top of . One can easily test for 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 only two isLeft() tests.
When P_{k} is exterior to H_{k–1}, we must then change to produce a new deque 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 . However, other points already in 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 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 push P_{k} onto both ends, producing the required new deque .
The speed of this algorithm is easily analyzed. Each vertex of 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 is the final hull). Thus, the Melkman algorithm is a very efficient and elegant time and space algorithm.
PseudoCode: Melkman Algorithm
Implementation
Here is a "C++" implementation of this 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(): test if a point is LeftOnRight 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 = n2, 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[top1], 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[top1], 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 <= (topbot); h++) H[h] = D[bot + h];
delete D; return h1; }
References
Greg Aloupis, A History of Lineartime 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, 296298 (1978)
Ronald Graham & F. Yao, "Finding the convex hull of a simple polygon", J. Algorithms 4(4), 32433 (1983)
D.T. Lee, "On finding the convex hull of a simple polygon", Int'l J. Comp. & Info. Sci. 12(2), 8798 (1983)
D. McCallum and D. Davis, "A linear algorithm for finding the convex hull of a simple polygon", Info. Proc. Letters 9, 201206 (1979)
A. Melkman, "Online construction of the convex hull of a simple polygon", Info. Proc. Letters 25, 1112 (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, p1355 (1972).
