Convex hull

Below are all earthquakes of magnitude five and above for 2015–2025, more than nineteen thousand points. Nobody drew a map of the continents, but it emerges by itself: the epicentres lie along the boundaries of tectonic plates, and the Pacific is outlined by the ring of fire. Choose a region, and the algorithm builds a fence around its epicentres, the shortest convex polygon with all the points inside. See how many of them the fence rests on.

Of almost two hundred epicentres of the Tian Shan and the Pamirs, ten hold the fence; the rest are inside. Such a polygon is called the convex hull. If you drive nails into the points and stretch a rubber band around them, it takes the shape of the convex hull.

The idea

Sort the points from left to right and go through them, collecting the lower half of the fence on a stack. Going round a convex polygon counterclockwise, you keep turning left. So if the last two points of the stack and the new point make a right turn (or go straight), the middle point is not a vertex: throw it out and check again. Then do the same from right to left for the upper half.

This is Andrew's monotone chain. Each point gets onto the stack once and comes off it at most once, so the pass is linear, and all the time goes into sorting.

Left or right

All the algorithm needs to know about geometry is which way the path O → A → B turns. The answer is given by the sign of the cross product of the vectors OA and OB:

Greater than zero is a left turn, less a right one, zero means the three points are on one line. No angles, no trigonometry, no division: only multiplications and subtractions. For integer coordinates the answer is exact, and that matters: comparing angles with atan2 lets floating-point error into the algorithm. More on vectors and orientation in vectors and orientation.

Why the result is the hull

The leftmost and rightmost points (after sorting, the first and the last) certainly lie on the hull: no point is further left or right than they are. The lower chain joins them, turning left all the time, so it is convex. No point thrown out can be a vertex: it was thrown out when the segment from the previous point of the chain to the new one passed below it or through it, that is, it turned out to be above that segment, inside. And all points not thrown out lie no lower than the chain, because any point that turned out to be below would itself, on arrival, have thrown out those in the way. The upper chain is the same, mirrored, and together they give the whole fence.

The condition “≤ 0” in the loop also throws out points lying exactly on an edge: the answer has only real corners. If a problem needs all points on the boundary, “< 0” is used, but then coinciding points need more care.

Code

Both halves in one function, almost word for word in Python and Go. The area is by the shoelace formula: the sum of cross products of neighbouring vertices.

And this is the first scene in numbers: the epicentres of the Tian Shan and the Pamirs from the USGS catalogue, converted from degrees to kilometres. At a latitude of forty degrees a degree of longitude is almost a quarter shorter than a degree of latitude, and without the correction the fence would come out squashed.

What it costs

Sorting is O(n log n), the pass O(n): each point is pushed onto the stack at most once and popped at most once. In total O(n log n).

Nothing faster is possible in general, and the proof is neat: if the hull could be built faster, sorting could be done faster than n log n. Put the numbers x₁, …, x_n on a parabola, as points (x_i, x_i²). All of them end up on the hull, and walking the hull lists them in order. And comparison sorting cannot be done faster than n log n.

You can win, though, when the hull has few points. Jarvis's “gift wrapping” algorithm (1973) runs in O(n·h), where h is the number of hull vertices: for the epicentres of the Tian Shan, where h = 10, this is better. And Chan's algorithm (1996) runs in O(n log h), better than both. Press “All on a circle” in the sandbox: there h = n, the worst case for Jarvis.

Where it comes in handy

• Geodata. The outline of a cloud of GPS points, the coverage area of stations, the area where a phenomenon was observed: the first thing computed on receiving a set of coordinates.
• Collisions. Physics engines first check whether the convex hulls of objects intersect: with convex shapes everything is much simpler to compute.
• The farthest pair. The two farthest points of a set always lie on the hull, and they can be found with rotating calipers in linear time after the hull is built. That is the second problem.
• Statistics. Peeling hulls off layer by layer, like an onion, gives the “depth” of a point in a cloud, an outlier-resistant analogue of the median for the plane.

History

In the late sixties the best known methods built a hull in time on the order of n². Ronald Graham was then at Bell Labs, and, as the story is usually told, hulls of sets of many thousands of points were needed there, for which a quadratic algorithm was too slow. Graham sorted the points by angle around the lowest one and went through them with a stack; the 1972 paper was the first to give an O(n log n) bound. Graham recalled considering it a minor piece of work, and it became one of the most cited of a long career. In 1979 Andrew replaced the sort by angle with a sort by coordinate: without angles the code is simpler and does not suffer from rounding errors.

Problems

In the first, build the hull and compute its area exactly, without fractions. In the second, find the farthest pair of points among two hundred thousand, where trying pairs will not fit in the time.