Ramer–Douglas–Peucker line simplification

Above Almaty, in the foothills and gorges of the Trans-Ili Alatau, OpenStreetMap knows 1556 paths and tracks. Some were drawn from hikers’ GPS tracks (the data says so), some from satellite images. Together they hold 58,758 points: on average one per 15 metres, and on the most detailed trails one every two or three metres, with every turn of a switchback. To show them on a phone map or store them in an app you do not need that many: on screen, points a couple of metres apart merge into one. Move the tolerance slider: all the trails are simplified at once, and the selected one shows its corridor; no original point is farther from the simplified line than the tolerance.

At a tolerance of five metres, 17,621 of the 58,758 points remain, 3.3 times fewer, and the trail «3 брата — Кумбель» (Three Brothers – Kumbel) shrinks from 450 points to 33. If you take the same number of points by simply keeping every fifteenth, the line strays 14.5 m from the trail, while with this algorithm it strays at most 4.9 m. It chooses which points matter for the shape and spends them where the trail turns.

The idea

Keep the first and last points and draw a chord between them. Find the point farthest from the chord. If it is closer than the tolerance, all the points between the ends can be thrown out: the chord already passes near each of them. If it is farther, that point is needed; keep it and do the same for the two halves of the line, from the start to it and from it to the end.

This is the Ramer–Douglas–Peucker algorithm. The tolerance ε is its only parameter, and its meaning is clear: how far the simplified line may stray from the original. The distance here is from a point to the chord segment: to the perpendicular if it falls on the segment, otherwise to the nearer end. The classic description uses the distance to the line, but on switchbacks that is dangerous: a point far beyond the end of a chord can lie close to its extension and vanish. Below is the algorithm on one trail, chord by chord.

The grade 1B climbing route up Shkolnik Peak: 392 points over 679 metres, a point every 1.7 m. With a 5 m tolerance the first chord joins the start and end of the route, and then the switchback unfolds into ever shorter chords until each is within the tolerance of its points. 24 points remain.

Why the shape is not spoiled

The guarantee is one-sided, but it is enough: a corridor of width 2ε can be drawn around the simplified line, and the whole original line lies inside it. What the algorithm does not promise is the smallest number of points. It is greedy from the top down: once it keeps a point, it never reconsiders. The smallest number of points for the same tolerance can also be found exactly: build a graph with an edge i → j if the segment i–j passes within ε of all the points between them, and find the shortest path in it from the first point to the last. Hiroshi Imai and Masao Iri studied such exact methods in 1986, but even their best descendants take on the order of n² operations: all pairs of points have to be checked. Ramer–Douglas–Peucker is much faster on real lines but may keep more points than necessary. Nor does it promise the smallest deviation for a given number of points: with a very coarse tolerance, when only three or four points of a trail remain, uniform sampling sometimes lies closer; the first scene shows this at a 50 m tolerance.

Code

Recursion is natural here, but on a long line it may go too deep, so the listing uses a stack of chords instead. The example is a short switchback of thirteen points: the program prints what it decided on each chord. The second part is a pass that will come in handy below: a single run with zero tolerance gives each point its threshold.

And this is the first scene in numbers: all trails at once, tolerances from one metre to fifty, the length before and after, and a comparison with “every k-th point”.

The trap: a switchback gets shorter

The simplified line cuts corners. For the shape on a map that is what we want, but not for the length: a chord is always shorter than the arc it replaces. Over the whole trail network the loss is small: 895.5 km become 884.1 km at a 5 m tolerance, 1.3% less, and 835.1 km at 50 m. But on a switchback every turn is a few metres that a chord cuts off entirely. The 679 m route up Shkolnik is 564 m after simplification with a 5 m tolerance, 17% shorter, and 496 m with a 20 m tolerance. An app that computes the length of a hike from a simplified track will promise a quarter less walking than there will actually be.

Hence the rule: length, elevation gain and time are computed from the original track, and simplification is only for drawing and storing the shape. The second well-known trouble is topology: the simplified line may cross itself or a neighbouring line, although the originals did not cross. We checked the trails of this article: none of them crossed itself at either 5 or 50 metres; the distance to a segment is more careful here than the distance to a line. But for parcel boundaries, where a crossing is an error in the data, variants of the algorithm that check for it are used.

All zoom levels in one pass

A map shows trails at different zoom levels, and each needs its own tolerance, roughly a pixel. There is no need to rerun the algorithm for each level. The tree of chord splits does not depend on the tolerance: the point farthest from a chord is the same for any ε. Only where the splitting stops changes. So one pass with zero tolerance records for each point its distance to its chord, and at tolerance ε the points that remain are those whose distance exceeds ε for themselves and for all their ancestors in the tree.

The words about ancestors matter. Distances need not decrease down the tree: in the example from the listing point 5 is 11.70 from its chord, while its “parent”, point 2, is only 7.15 from its own. At a tolerance of 8 point 2 is thrown out, the chord 0–12 stays whole, and point 5 is never reached. So a point’s threshold is the smallest of the distances on the path from the root: min(11.70, 7.15) = 7.15. The listing checks that such a filter gives exactly the same as a full run, at every tolerance. The first scene works just this way: the slider only compares thresholds with the tolerance. TopoJSON has a similar trick: each point is given a weight in advance, the area of the triangle with its neighbours as in the Visvalingam–Whyatt algorithm, and simplification for a zoom level drops the points with a weight below the threshold.

What it costs

Each chord is a pass over the points between its ends. If the farthest point splits the line roughly in half, there are about log₂ n levels of splitting, with O(n) work in total at each: O(n log n) in all. Real trails behave this way: the second scene shows how quickly the chords get shorter. But if the farthest point always stands next to an end of the chord, each split nips off one point, and the total is 1 + 2 + … + n ≈ n²/2. Such a line exists: a zigzag converging to its end, where each next vertex is slightly closer to the axis than the previous one.

At a thousand points the zigzag costs tens of times more than the trail, and the gap grows with n. John Hershberger and Jack Snoeyink showed in 1992 how to make the worst case O(n log n) too: the point farthest from a chord always lies on the convex hull of the piece of line, and if the hulls of the pieces are maintained during splitting, it can be found quickly. Memory is O(n): the points themselves, “keep” flags and a stack of chords with at most n entries.

Where it works

• GIS and cartography. The function ST_Simplify in PostGIS, simplify in Shapely and DouglasPeuckerSimplifier in the JTS library simplify geometry with this algorithm.
• GPS tracks. Watches and phones often record a run or a hike as a point a second, while for a map and storage the points at the turns are enough.
• This catalogue. The street shapes on the Almaty map in the scenes with Dijkstra and A* were simplified by it with a 3 m tolerance; it shows in the data preparation script.
• Robots. A laser rangefinder gives a cloud of points on walls, and the same trick, “split at the farthest point”, extracts straight segments from it.

History

The algorithm was invented twice, almost at the same time. Urs Ramer published it in 1972 in the journal Computer Graphics and Image Processing as a way to approximate plane curves by polylines. A year later David Douglas and Thomas Peucker described the same trick in a cartographic journal for digitised lines on maps; their paper was called “Algorithms for the reduction of the number of points required to represent a digitized line or its caricature”. In the pattern recognition literature the same method is known as iterative end-point fit.

Problems

In the first problem the algorithm must be implemented exactly, in integers, without a single floating-point operation, and with the distance to a segment rather than a line. In the second you must say how many points remain on the whole trail map for thousands of different tolerances at once: running the algorithm for each tolerance takes too long.