Divide-and-conquer polygon cleaning

clip.cpp

@@ -571,8 +571,7 @@ drawvec remove_noop(drawvec geom, int type, int shift) {
 				}
 
 				if (geom[i + 1].op == VT_CLOSEPATH) {
-					// followed by closepath: not possible
-					fprintf(stderr, "Shouldn't happen\n");
+					// followed by closepath: only possible after close_poly()
 					i++;  // also remove unused closepath
 					continue;
 				}
@@ -2469,3 +2468,87 @@ drawvec fix_polygon(const drawvec &geom, bool use_winding, bool reverse_winding)
 
 	return out;
 }
+
+bool line_is_too_small(drawvec const &geometry, int z, int detail) {
+	if (geometry.size() == 0) {
+		return true;
+	}
+
+	long long x = 0, y = 0;
+	for (auto &g : geometry) {
+		if (g.op == VT_MOVETO) {
+			x = std::llround((double) g.x / (1LL << (32 - detail - z)));
+			y = std::llround((double) g.y / (1LL << (32 - detail - z)));
+		} else {
+			long long xx = std::llround((double) g.x / (1LL << (32 - detail - z)));
+			long long yy = std::llround((double) g.y / (1LL << (32 - detail - z)));
+
+			if (xx != x || yy != y) {
+				return false;
+			}
+		}
+	}
+
+	return true;
+}
+
+void coalesce_polygon(drawvec &geom, bool scale_up) {
+	// wagyu should be able to straightforwardly handle
+	// anything under a few hundred thousand vertices
+	if (geom.size() < 100000) {
+		geom = clean_or_clip_poly(geom, 0, 0, false, scale_up);
+		return;
+	}
+
+	// These geometries were assembled in geometric order,
+	// so sub-batches of them should hopefully union into
+	// reasonable sets.
+	//
+	// Find the first outer ring after halfway point.
+
+	for (size_t i = geom.size() / 2; i < geom.size(); i++) {
+		if (geom[i].op == VT_MOVETO) {
+			size_t j;
+			for (j = i + 1; j < geom.size(); j++) {
+				if (geom[j].op != VT_LINETO) {
+					break;
+				}
+			}
+
+			if (get_area(geom, i, j) > 0) {
+				// If we have an outer ring, split there
+				// and coalesce the two halves
+
+				// Copy second half to new vector
+				std::vector<draw> geom2;
+				geom2.resize(geom.size() - i);
+				for (size_t k = i; k < geom.size(); k++) {
+					geom2[k - i] = geom[k];
+				}
+
+				// Resize vector to include only first half
+				geom.resize(i);
+
+				// Clean each half individually
+				coalesce_polygon(geom, scale_up);
+				coalesce_polygon(geom2, scale_up);
+
+				// Copy second half back with first
+				size_t brk = geom.size();
+				geom.resize(brk + geom2.size());
+				for (size_t k = 0; k < geom2.size(); k++) {
+					geom[brk + k] = geom2[k];
+				}
+				geom2.clear();
+
+				// Clean the combined geometry
+				geom = clean_or_clip_poly(geom, 0, 0, false, scale_up);
+			}
+
+			i = j - 1;
+		}
+	}
+
+	// Can't find a breakpoint; take what we can get.
+	geom = clean_or_clip_poly(geom, 0, 0, false, scale_up);
+}

geometry.hpp

@@ -173,4 +173,7 @@ void get_quadkey_bounds(long long xmin, long long ymin, long long xmax, long lon
 
 clipbbox parse_clip_poly(std::string arg);
 
+bool line_is_too_small(drawvec const &geometry, int z, int detail);
+void coalesce_polygon(drawvec &geom, bool scale_up);
+
 #endif

tile.cpp

@@ -639,7 +639,7 @@ static double simplify_feature(serial_feature *p, drawvec const &shared_nodes, n
 					// unioned exactly
 					//
 					// don't try to scale up because these are still world coordinates
-					geom = clean_or_clip_poly(geom, 0, 0, false, false);
+					coalesce_polygon(geom, false);
 				}
 
 				// continues to simplify to line_detail even if we have extra detail
@@ -695,7 +695,7 @@ static void *simplification_worker(void *v) {
 
 				if (!a->trying_to_stop_early) {
 					// we can try scaling up because this is now tile scale
-					geom = clean_or_clip_poly(geom, 0, 0, false, true);
+					coalesce_polygon(geom, true);
 					if (additional[A_DEBUG_POLYGON]) {
 						check_polygon(geom);
 					}
@@ -1551,26 +1551,6 @@ bool find_feature_to_accumulate_onto(std::vector<std::shared_ptr<serial_feature>
 	return false;
 }
 
-static bool line_is_too_small(drawvec const &geometry, int z, int detail) {
-	if (geometry.size() == 0) {
-		return true;
-	}
-
-	long long x = std::round((double) geometry[0].x / (1LL << (32 - detail - z)));
-	long long y = std::round((double) geometry[0].y / (1LL << (32 - detail - z)));
-
-	for (auto &g : geometry) {
-		long long xx = std::round((double) g.x / (1LL << (32 - detail - z)));
-		long long yy = std::round((double) g.y / (1LL << (32 - detail - z)));
-
-		if (xx != x || yy != y) {
-			return false;
-		}
-	}
-
-	return true;
-}
-
 // Keep only a sample of 100K extents for feature dropping,
 // to avoid spending lots of memory on a complete list when there are
 // hundreds of millions of features.
@@ -2201,7 +2181,7 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 								drawvec to_clean = features[simplified_geometry_through]->geometry;
 
 								// don't scale up because this is still world coordinates
-								to_clean = clean_or_clip_poly(to_clean, 0, 0, false, false);
+								coalesce_polygon(to_clean, false);
 								features[simplified_geometry_through]->geometry = std::move(to_clean);
 							}
 						}
@@ -2470,7 +2450,7 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 					if (layer_features[x]->t == VT_POLYGON) {
 						if (layer_features[x]->coalesced) {
 							// we can try scaling up because this is tile coordinates
-							layer_features[x]->geometry = clean_or_clip_poly(layer_features[x]->geometry, 0, 0, false, true);
+							coalesce_polygon(layer_features[x]->geometry, true);
 						}
 
 						layer_features[x]->geometry = close_poly(layer_features[x]->geometry);

unit.cpp

@@ -152,3 +152,12 @@ TEST_CASE("mvt_geometry bbox") {
 	REQUIRE(start == 0x1c84fc0000000000);
 	REQUIRE(end == 0x1c84ffffffffffff);
 }
+
+TEST_CASE("line_is_too_small") {
+	drawvec dv;
+	dv.emplace_back(VT_MOVETO, 4243099709, 2683872952);
+	dv.emplace_back(VT_LINETO, 4243102487, 2683873977);
+	dv.emplace_back(VT_MOVETO, -51867587, 2683872952);
+	dv.emplace_back(VT_LINETO, -51864809, 2683873977);
+	REQUIRE(line_is_too_small(dv, 0, 10));
+}