Faster than earcut: triangulating 3.7M map polygons in 489 ms

We have published the polygon triangulator of our map engine as a small library, earstrip, under the ISC license, the same one earcut uses. It is one C++ header, and a WebAssembly build on npm that takes the same arguments as earcut. In C++ it is 14% faster than earcut.hpp, in the browser 1.4 to 2.1 times as fast as earcut.js, and the comparison is near the end. This post is about how we got there.

A GPU fills triangles and nothing else. A map, on the other hand, is mostly polygons: buildings, parks, parking lots, the strip of grass between two tram tracks. Here are a few blocks of Śródmieście in Warsaw the way our renderer reads them from a map file, as outlines:

Before any of this gets a fill, every outline has to be cut into triangles. Hover over a polygon, or tap it, to see what it turns into. These blocks hold 331 polygons with 5,151 vertices, and they become 4,514 triangles.

Zoom level 13 of our offline map of the Masovian Voivodeship, the region around Warsaw, has 3,740,425 polygons. On one core of an M4 Pro, the constrained Delaunay triangulator in our engine needs 2.93 seconds to cut all of them. This October I replaced it with ear clipping. To write this post I rebuilt the ear clipper one idea at a time, one commit per idea, and timed every commit on those same 3.7 million polygons, so each section below ends with a number.

A triangulation of a set of points is Delaunay when no point lies inside the circle through the three corners of any triangle. Two triangles that share an edge are enough to see how the rule works. Drag the point D:

45°ABCDSmallest angle 45°; with BD instead it would be 37°

While D stays outside the dashed circle through A, B and C, the shared edge is AC. Push D into the circle and the edge flips to BD; pull it back out and it flips back. Each flip is a local repair, and a triangulation where no edge wants to flip is Delaunay. Of all the ways to triangulate the same points, it has the largest smallest angle. The orange arc marks the smallest angle of the two triangles, and the line under the figure shows what the other diagonal would give: wherever you put D, the edge the circle rule picks wins.

"Constrained" means that the polygon's own edges must stay edges of the triangulation, even where the circle rule would flip them away. Our triangulator inserts the vertices one at a time. For each new point it walks across the triangles built so far until it finds the one the point falls into, splits it, and flips edges around the new point until the circle rule holds again. It also forces the polygon's edges in, and at the end it keeps only the triangles that lie inside the polygon and outside its holes.

Well-shaped triangles matter when you interpolate something across them or run a simulation on them. A flat fill doesn't care: a GPU fills a needle-thin triangle as correctly as a fat one. And most of the work is spent on polygons that hardly need any. Almost half of the Masovian polygons, 1.8 million, have four vertices or fewer. A triangle needs no cutting at all and a convex quad is two triangles, yet the Delaunay code spends 419 ms on them, building each triangulation from scratch.

An ear of a polygon is a corner you can cut off with one straight cut. Take three consecutive vertices a, b and c. If the corner at b is convex and no other part of the polygon reaches into the triangle abc, that triangle is an ear. Cut it off, and what remains is a polygon with one vertex fewer. In 1975 Gary Meisters proved that every simple polygon with more than three vertices has at least two ears that don't overlap. So there is always another ear to cut, and a polygon with n vertices falls apart into n − 2 triangles.

Here is one of the buildings from the blocks above. Move the slider to cut its ears one at a time, or cut them all:

0 / 25

The green triangle is the next ear to go. The algorithm needs one test, "is this corner an ear?", and the simplest version of it checks the corner and then every other vertex of the boundary. cross is the 2×2 determinant found in every geometry codebase; it is positive when c lies to the left of the line from a to b, so for a boundary that runs counter-clockwise a positive cross(a, b, c) means a convex corner:

static double cross(double ax, double ay, double bx,

double by, double cx, double cy)

{

return (bx - ax) * (cy - ay) - (by - ay) * (cx - ax);

}

bool isEar(uint32_t v) const

{

const Node &b = _nodes[v];

const Node &a = _nodes[b.prev], &c = _nodes[b.next];

if (cross(a, b, c) <= 0)

return false; // turns the wrong way, or not at all

for (uint32_t w = c.next; w != b.prev;

w = _nodes[w].next) {

const Node &p = _nodes[w];

if (cross(a.x, a.y, b.x, b.y, p.x, p.y) < 0 ||

cross(b.x, b.y, c.x, c.y, p.x, p.y) < 0 ||

cross(c.x, c.y, a.x, a.y, p.x, p.y) < 0)

continue; // outside the triangle

if (!samePos(p, a) && !samePos(p, b) &&

!samePos(p, c))

return false;

}

return true;

}

A point on an edge of the triangle counts as inside: if the ear touched the rest of the boundary, cutting it off would leave a polygon that touches itself. The one exception is a point at exactly the position of a, b or c, which happens where an outline meets itself at a single vertex. Holes make that common.

Ear clipping works on one closed boundary, and 11,345 of the Masovian polygons have holes, such as courtyards and clearings in forests. The way to deal with them is to cut each hole open and stitch it to the outer ring with a bridge, an edge from a vertex of the hole to a vertex of the boundary, walked once in each direction. The boundary then runs around the outside, crosses the bridge, goes around the hole the other way and comes back. There is no area between the two copies of the bridge, so the polygon is the same, but now it is one ring.

The bridge must not cross anything. David Eberly's write-up of ear clipping has a neat way to find one. Take the rightmost vertex of the hole, m, and cast a ray to the right. The first edge the ray hits, at I, has an end P further to the right. If no boundary vertex lies inside the triangle m, I, P, the bridge m–P is safe. You can drag the hole around:

IPmThe tip of the dent lies in the triangle m, I, P. A bridge to P would cross the boundary, so it goes to the dent.

Higher up the triangle is empty and the bridge goes to P. Lower down the tip of the dent sits inside the triangle, and the bridge goes to the dent instead: of the dents in the triangle, Eberly takes the one closest in angle to the ray. Holes are bridged rightmost first, so a hole further left can bridge to one that is already part of the boundary.

Map data adds two cases of its own. Our coordinates are integers on a tile grid, so a ray along a grid line runs straight through a vertex all the time; the crossing test passes over such a vertex, and my first version got exactly this case wrong. And mappers often draw a courtyard that touches the outer wall at one point. Such a hole needs no bridge: the two rings are joined at the shared vertex.

With holes bridged, the first complete version walks around the boundary, tests corners with the isEar above, cuts every ear it finds and writes each triangle as three indices. On the Masovian polygons it takes 1.12 seconds, 2.6 times faster than Delaunay. The chart splits each time by polygon size; hover over a part of a bar to see its number, or click a row:

The naive row is shorter than Delaunay in every size class but the last. The 1,123 largest polygons, with 257 vertices or more, take 199 ms to clip and 90 ms with Delaunay. Each ear test there walks hundreds of vertices, and a polygon needs hundreds of tests, so the cost grows with the square of its size. The rest of this post is about making the ear test cheaper and needed less often.

Look at which vertices can actually stop a convex corner from being an ear. Say some vertices lie inside its triangle, and take the one closest to the corner itself. Its two edges can't cross the sides of the triangle, so both run back, away from the corner, and the polygon wraps around that vertex from the corner's side. That makes it a reflex vertex, a dent. So if anything is inside the triangle, a dent is inside it too, and the test only has to look at the dents.

The engine keeps a list of them. A dent leaves the list when cutting its neighbour makes it convex, and a convex corner never turns into a dent, so the list only shrinks. The test becomes:

bool isEar(uint32_t v) const

{

const Node &b = _nodes[v];

if (b.reflex)

return false;

const Node &a = _nodes[b.prev], &c = _nodes[b.next];

for (uint32_t i = 0; i < _live; ++i) {

const Item &p = _items[i]; // a dent

if (p.node == b.prev || p.node == b.next)

continue;

if (cross(a.x, a.y, b.x, b.y, p.x, p.y) < 0 ||

cross(b.x, b.y, c.x, c.y, p.x, p.y) < 0 ||

cross(c.x, c.y, a.x, a.y, p.x, p.y) < 0)

continue;

if ((p.x != a.x || p.y != a.y) &&

(p.x != b.x || p.y != b.y) &&

(p.x != c.x || p.y != c.y))

return false;

}

return true;

}

The 28.1 million ear tests now look at 185 million points instead of 353 million, and the total drops to 882 ms. That is less than you might expect, because building outlines have a lot of dents. The building with the arch above has 15 of them among its 27 vertices. And the smallest polygons get slower, from 34 to 54 ms: for a triangle or a quad, sorting the corners into dents and the rest costs more than the walk it saves.

A polygon without dents is convex, and every corner of a convex polygon is an ear, so it can be cut without a single test. Checking for that needs care. "Every corner turns left" is not enough:

PentagonStar

Corners turning left: 5 of 5. Changes of x direction: 4. Not convex.

The star turns left at all five corners, because it goes around its centre twice. The engine adds one more condition. Walking around a convex ring, the edges go right along the bottom and left along the top, so their x direction changes exactly twice, at the leftmost and the rightmost point. The pentagon passes; the star changes direction four times, at the circled corners, and goes to the ordinary ear tests. Both counts come out of the same loop that marks the dents, so the check costs almost nothing. Triangles and quads get their own shortcut before that loop.

The shortcut removes 4.4 million of the 28.1 million ear tests, but those were the cheapest ones: a ring without dents has nothing to scan. On its own it barely moves the total, 880 ms, though the polygons with four vertices or fewer win back part of what the dents cost them, from 54 to 48 ms. It matters more in the next step, where a convex ring becomes a single strip.

So far every triangle costs three indices. A GPU can also draw a triangle strip, where each new index adds a triangle made of it and the two before. A convex ring makes one strip if you zig-zag across it:

012345670 1 7 2 6 3 5 4

The ear clipper can write strips too. Cutting an ear makes the corners of its two neighbours sharper, so they are good candidates for the next ear, and cutting one of them extends the strip by a single index. The engine tries them first and keeps alternating between the two ends while the ears last. In its code the neighbours of the cut vertex are a and b:

// The next ear alternates between the ends of the cut.

while (true) {

if (sideA) {

const uint32_t z = n[a].prev;

_strips.push_back(n[z].index);

unlink(a);

a = z;

} else {

const uint32_t z = n[b].next;

_strips.push_back(n[z].index);

unlink(b);

b = z;

}

count--;

sideA = !sideA;

if (count < 3 || !isEar(sideA ? a : b))

break;

}

When the next corner is not an ear, the strip ends and the walk goes on from where it stopped. The 21.5 million triangles now take 35.8 million indices instead of 64.5 million. The ear tests don't get fewer: the naive walk already went on at the next corner after a cut, and now every strip also ends on a test that fails,