Cutting map polygons into triangles
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, :
Before any of this gets a fill, every outline has to be cut into . 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 3.76 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 and timed every step on those same 3.7 million polygons, so each section below ends with a number.
Delaunay
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:
While D stays outside the dashed circle through A, B and C, the shared edge is AC. and the edge flips to BD; 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.
Ears
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 :
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.
Holes
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:
the triangle is empty and the bridge goes to P. 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 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 at one point. Such a hole needs no bridge: the two rings are joined at the shared vertex.
The naive clipper
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.36
seconds, 2.8 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:
is shorter than in every size class but the last. The 1,123 largest polygons, with 257 vertices or more, take 251 ms to clip and 157 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.
Only dents can block
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 33.9 million ear tests now look at 231 million points instead of 378 million, and the total drops to . 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.
Polygons that need no tests
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:
The 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 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 33.9 million ear tests and brings the total to . The polygons with four vertices or fewer go from 63 to 56 ms.
Strips
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:
0 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. And because the neighbours of a cut ear are usually ears themselves, the clipper needs 25.0 million tests instead of 29.5 million. Together that is , and the largest polygons drop from 189 to 98 ms.
Buckets
A forest or a lake can have thousands of vertices and hundreds of dents, and each ear test still reads every one of them. When a polygon has more than 32 dents, the engine sorts them into buckets along the longer side of their bounding box, about two dents per bucket. An ear test looks only at the buckets that its triangle spans. As dents turn convex and leave, the index is rebuilt once half of it is gone.
Points scanned fall from 214 million to 122 million, almost all of it on the large polygons: the polygons with 257 vertices or more go from 98 to 48 ms. Total: .
Counting without branches
The loop in the dent test exits at the first dent inside the triangle. Most of the time there is none: the test reads all the dents and says yes. An early exit then buys nothing, and the three comparisons per point are branches the CPU has to predict. The engine first runs a pass that only counts dents inside the triangle, with no branches at all:
const double abx = b.x - a.x, aby = b.y - a.y;
const double bcx = c.x - b.x, bcy = c.y - b.y;
const double cax = a.x - c.x, cay = a.y - c.y;
uint32_t hits = 0;
for (uint32_t i = begin; i < end; ++i) {
const Item &p = items[i];
const double side = std::min(
std::min(abx * (p.y - a.y) - aby * (p.x - a.x),
bcx * (p.y - b.y) - bcy * (p.x - b.x)),
cax * (p.y - c.y) - cay * (p.x - c.x));
hits += (uint32_t)(side >= 0) & (uint32_t)(p.node != ia)
& (uint32_t)(p.node != ic);
}
if (hits == 0)
return true;
The smallest of the three determinants is non-negative exactly when the point is inside the closed triangle. Only when the count is not zero does the exact loop from before run, to rule out dents that sit on a corner. The bucket version does the same with a bounding-box test, 64 points at a time into a bit mask. This is the last step: , 6.8 times less than Delaunay on the same polygons. The middle sizes gain the most: polygons with 33 to 256 vertices go from 312 to 224 ms.
When there are no ears
Meisters' theorem is about simple polygons. Map data is rounded to a tile grid, and rounding can break a simple polygon. Here is the outline of Warszawa Śródmieście, a railway station less than a kilometre from the blocks above, rounded to grids of different size. The orange rings mark vertices that land on the same point as their neighbour:
per tile side, the grid our map files use, a grid step is 18 cm. Even there 31 of the 154 vertices merge with a neighbour, because mappers often place points a few centimetres apart, and the outline survives. nearly half of the vertices merge, and the clipper still cuts it. the rounded outline touches itself in a way that leaves no ear at all, and at 10 bits it happens to work again.
When the clipper walks all the way around the boundary without finding an ear, it drops the corners that have no area and walks once more. If that fails too, it gives up, and the engine triangulates that polygon with the Delaunay code, which still exists for this. In the Masovian data 76 polygons end up there. They come from a map built before the last change on the other end: our map generator now rounds every polygon to the final grid with GEOS's snap rounding and checks that ear clipping can cut it. In a fresh build of the Masovian map that repaired 273 polygons, and ear clipping can cut every polygon in the file.
Delaunay has one more job. Every hole needs its own bridge, and finding one scans the whole boundary, so a polygon with thousands of holes pays for thousands of scans. Our world map has a polygon with 11,675 holes: ear clipping takes 606 ms on it, Delaunay 17 ms. The engine now sends every polygon with 128 holes or more to Delaunay. In the Masovian data that is two polygons.
In the browser
Our web renderer runs the same C++ engine compiled to WebAssembly, so I also compared our clipper with earcut, the ear clipper behind MapLibre and Mapbox GL, in its C++ version earcut.hpp 3.2.4 and in JavaScript as earcut.js 3.2.4. This is a separate harness: it reads the same polygons from a flat file and feeds them to native, WebAssembly and JavaScript builds, which is why its native time is lower than in the chart above. Every build gives exactly the same triangles as its native version on all 4.4 million polygons of zoom levels 13 and 9. Milliseconds for the 3.74 million polygons of zoom 13, best of 7, one thread, V8 in Node and JavaScriptCore in Bun:
| Our clipper | earcut | |
|---|---|---|
| Native C++ | 510 | 561 |
| WebAssembly, loop inside the module, V8 / JSC | 599 / 586 | 964 / 850 |
| WebAssembly, one call from JS per polygon, V8 / JSC | 690 / 702 | 1,050 / 955 |
| JavaScript, V8 / JSC | 1,126 / 930 | 1,037 / 1,104 |
Natively earcut.hpp is 10% slower than our code. Compiled to WebAssembly, our clipper loses 17% and earcut.hpp loses 72%, almost all of it in its main loop that searches for ears. It isn't the allocator: builds with mimalloc and emmalloc give the same time. I haven't dug further. In V8, earcut.hpp in WebAssembly ends up only 7% faster than plain earcut.js, and once each polygon is a separate call from JavaScript, it is slower.
The JavaScript row for our clipper is a straight port of the C++ class. Its first version kept the nodes in parallel typed arrays, the way you would write it in C, and took 1,396 ms in V8: bounds checks on every access, and an ear test too big for the inliner. Plain objects reused from a pool brought it to about 1,100 ms.
What changed in the app
Loading the 4,260 zoom-13 tiles of the Masovian map for the first time took 5,681 ms before and takes 1,531 ms
now. The .cache file is gone. It held the triangles from the first load so that later loads could read them
back, and for this map it was 117 MB, 59% of the size of the map itself. Reading those triangles was still
faster than cutting them: a repeated load took 1,071 ms with the file and takes 1,517 ms without it. I took that
trade for the disk space, and the half-second may come back if triangulation becomes lazy and only covers the
polygons that are actually on screen.