Introduction and Problem Definition
Given n points in the plane, the Convex Hull can be defined as the smallest convex set containing them. This structure is important in computer graphics given its differents applications to shape analysis, collision detection, bounding box computing, among others. It is natural then, to question how efficiently we can perform the calculations in order to find the convex hull on a given set of points. It is now know that for an arbitrary input set of points, the solution can be found in a complexity time of O(nlogn). For this a sort algorithm of order O(logn) is needed as a first step. The remaining steps will require a time of order O(n) implying our claim of an algorithm of order O(nlogn). However, as McCallum and Avis first showed, it is possible to compute the convex hull with a complexity time of order O(n) whenever the input set is a simple polygon (polygon without self-intersections).
In this document, I going to explore a particular algorithm due to Avraham A. Melkman to find the convex hull of a simple polyline with complexity time O(n). A simple polyline is a curve made of lines without self-intersections.
The Melkman Algorithm
The Melkman algorithm was introduced by Avraham A. Melkman in 1985 in his paper Online Construction of the convex hull of a simple polyline. He presented and incremental algorithm computing with a complexity time of O(n). I'm going to use the same notations of his paper. The code presented here is based on an implmentation of the algorithm in C++ using the OpenGL API for the graphics.
A polyline of n points is defined as an array of points P=(v1,v2,...vn). The simplicity of the Melkman algorithm resides in the use of a deque as a main structure to store the points in the current convex hull. A deque is a double-ended queue, where you can pop and push from both top and bottom. It is noted as D=(Db,...,Dt) where Db means the bottom element and Dt the top element. The magic of the algorithm is a consequence of the polyline being a simple polyline. Let's ilustrate this with a simple example. Suposse we have just 3 points in the plane defining a simple polyline. It is obvious that the convex hull its just the triangle with vertices the 3 points.
Here the green lines indicates edges of the convex hull that coincide with edges of the polyline, the red line indicate a edge of the convex hull that is not part of the polyline. Now suppose we want to add a new vertex v4 and compute the convex hull of the new set of 4 vertices. Let's see the possibilities for this new vertex. Because the polyline is simple, we know that the new vertex cannot be to the left of the oriented line v1-v2, and to right of the oriented line v3-v1 at the same time ("v1-v2" means: From v1 to v2, and so on).
Also, if the point is inside the triangle, the convex hull remain the same. The thing here is that we know the point is inside the triangle just because the new vertex is to the right of the oriented line v3-v1 and also to the left side of the oriented line v3-v2 (notice that we don't need to check if the point is on the right side of the oriented line v1-v2, because we have a simple polyline!). Thus when inserting a 4th vertex, we know the convex hull will remain the same except when the point is in the gray region of the figure below.
Then, whe need to update the convex hull just in the cases II and III, i.e: when the point to the right of the oriented line v3-v2 or when the point is to the left of the oriented line v3-v1. Because checking this is of order O(1), doing this process incrementally along a simple polyline of n points will result in the desired algorithm! of O(n). We will see how the deque structure will help the update step of the convex hull. Let's see the code now. The following routine takes 3 points a,b,c as input and returns 1 or 0 depending on whether the point c is to the right or left, of the oriented line from a to b (we are assuming like in the paper, that there no exist three collinear sucesive points):
Observe that we are just computing the vector product between the vectors (c-b) and (b-a), and querying about the sign of the z-component. Below is an excerpt of the Melkman algorithm. Here the simple polyline is represented by the array points[length] of planar points. To download the complete source code in C++ with OpenGL (on linux) click here. Notice that we are assuming that the points in the final convex hull will appear in clockwise order, as in the original paper.
Correctness of the Algorithm
It is clear that the melkman algorithm in linear in the number of points of the input polyline, because each vertex is processed (pushed-poped) at most once. But why this algorithm works?. To answer this, first remember the next characterization of the convex hull H of a simple polyline P: (i) H is convex, (ii) P is contained in H, (iii) the set of vertices of H is a subset of the set of verices of P. In his paper, Melkman proved by induction that the following hypothesis is true each time a new vertex appear in the first while loop:
First observe that all vertices rejected in the first while loop are inside, or on the current convex hull contained in D. In fact if a vertex v is rejected there it is because both v_cross_sgn(db,db+1,v)>0 and v_cross_sgn(dt-1,dt,v)>0, i.e v is to the right of both oriented edges db-(db+1) and (dt-1)-dt. Moreover the polyline P connects dt with v, and because P is simple, the edge v-dt do not intersect the part of P between db+1 and dt-1. Hence v must lie inside the polygon inscribed by D.
Let's proceed by induction. The part of the algorithm before the for loop ensures that the hypothesis HP holds when three vertices have been read (assuming they are not collinear). By induction suppose that HP holds when k vertices have been read. Because of HP and the previous paragraph we can assume that v satisfies either (db,db+1,v)=-1 or (dt-1,dt,v)=-1 (or both) and than we are going to step over the next two while loops. Now suppose that we have the deque D=(db,...dt) and that after v is read and processed the new deque is D'=(dk,...,dm) with dk=dm=v (from our previuos assumption). Let's proof (i),(ii) and (iii) for D'.
It is easly seen by induction that P is contained in D', because it is contained in D and each time a vertex dj is popped; the polygon (v,di,...dj,v) is contained within the polygon (v,di,...,dj+1,v). Similarly for the removal of a vertex. To proof D' to be convex, we show that it describes a closed simple curve and that (di,di+1,di+2)>0 for i=k,...m-2, and (dm-1,dm,dk+1)>0. Because D is simple, D' can only be nonsimple if the edges (v,dk+1) and (dm-1,v) intersect edges of D, and this cannot be the case since P is simple. Also (di,di+1,di+2)>0 for i=k,...m-2, is either already true for D or else it is the result of the two final while loops. For (dm-1,dm,dk+1)>0 suppose by contradiction that 0>=(dm-1,dm,dk+1). Then v=dm is on or to the right of the oriented edge (dm-1)-(dk+1) and at the same time it must be to the right of the polyline (dk+1,....,dm-1) implying that v is inside D, a contradiction. This conludes the proof of HP.