repair corrupt (multi)polygons as far as possible. Fixes many clipping (#908)

Author: geoneutrino

include/geom.h

@@ -77,6 +77,12 @@ void make_valid(GeometryT &geom) { }
 
 void make_valid(MultiPolygon &mp);
 
+// Attempt to repair an invalid areal geometry in place: dissolve-based
+// make_valid first (preserves area), then a zero-width buffer as a last
+// resort. Returns true if mp is valid afterwards; on failure mp is left as the
+// best-effort input so callers never regress.
+bool repair_multi_polygon(MultiPolygon &mp);
+
 void union_many(std::vector<MultiPolygon> &mps);
 
 Point intersect_edge(Point const &a, Point const &b, char edge, Box const &bbox);

src/geom.cpp

@@ -3,6 +3,8 @@
 
 #include <boost/geometry/geometries/segment.hpp>
 #include <boost/geometry/index/rtree.hpp>
+#include <boost/geometry/algorithms/buffer.hpp>
+#include <boost/geometry/strategies/buffer.hpp>
 
 #include "geometry/correct.hpp"
 
@@ -144,6 +146,83 @@ void make_valid(MultiPolygon &mp)
 	mp = result;
 }
 
+// Repair a single (possibly invalid) polygon in an area-preserving way and
+// append the resulting valid polygon(s) to `out`. Returns true on success.
+// `minArea` is the lower bound on the repaired area we are willing to accept.
+static bool repair_one_polygon(const Polygon &p, double minArea, MultiPolygon &out)
+{
+	// 1) Dissolve (resolves self-intersections of this single polygon).
+	try {
+		MultiPolygon fixed;
+		geometry::correct(p, fixed, 1E-12);
+		if (geom::is_valid(fixed) && std::abs(geom::area(fixed)) >= minArea) {
+			for (auto &fp : fixed) out.push_back(std::move(fp));
+			return true;
+		}
+	} catch (const std::exception &) {
+		// fall through to the buffer attempt
+	}
+
+	// 2) Zero-width buffer as a last resort.
+	try {
+		MultiPolygon buffered;
+		geom::strategy::buffer::distance_symmetric<double> distanceStrategy(0.0);
+		geom::strategy::buffer::side_straight sideStrategy;
+		geom::strategy::buffer::join_miter joinStrategy;
+		geom::strategy::buffer::end_flat endStrategy;
+		geom::strategy::buffer::point_square pointStrategy;
+
+		geom::buffer(p, buffered, distanceStrategy, sideStrategy, joinStrategy, endStrategy, pointStrategy);
+		geom::correct(buffered);
+		if (geom::is_valid(buffered) && std::abs(geom::area(buffered)) >= minArea) {
+			for (auto &bp : buffered) out.push_back(std::move(bp));
+			return true;
+		}
+	} catch (const std::exception &) {
+		// keep best-effort polygon
+	}
+
+	return false;
+}
+
+bool repair_multi_polygon(MultiPolygon &mp)
+{
+	if (geom::is_valid(mp)) return true;
+
+	// Repair PER POLYGON, area-preserving. Running make_valid/buffer on the whole
+	// multipolygon can catastrophically COLLAPSE large/complex inputs: a clipped
+	// reservoir with >1000 rings dropped ~99% of its area, leaving missing lake
+	// tiles at low zoom. Conversely, simply keeping the whole invalid geometry
+	// lets a single self-intersecting ring render as a spurious "spike".
+	//
+	// Fixing each polygon independently gets the best of both: a self-touching
+	// ring is cleaned (no spike) while the rest stays intact, and we avoid the
+	// O(n^2) cross-polygon union that caused the collapse. A polygon whose repair
+	// would not preserve its area is kept as-is (invalid but complete renders;
+	// only a tiny local artefact, never a dropped area).
+	MultiPolygon out;
+	bool allValid = true;
+	for (const auto &p : mp) {
+		if (geom::is_valid(p)) {
+			out.push_back(p);
+			continue;
+		}
+		// Lenient threshold: resolving a self-intersection legitimately changes a
+		// single polygon's (shoelace) area, so anything down to half the original
+		// is accepted. Per-polygon repair cannot trigger the cross-polygon union
+		// that previously caused the catastrophic ~99% collapse, so this only
+		// rejects a genuine local collapse.
+		const double minArea = 0.5 * std::abs(geom::area(p));
+		if (!repair_one_polygon(p, minArea, out)) {
+			out.push_back(p);
+			allValid = false;
+		}
+	}
+
+	mp = std::move(out);
+	return allValid;
+}
+
 // ---------------
 // Union multipolygons
 // from https://github.com/boostorg/geometry/discussions/947

src/tile_data.cpp

@@ -1,5 +1,6 @@
 #include <algorithm>
 #include <iostream>
+#include <sstream>
 #include "tile_data.h"
 #include "coordinates_geom.h"
 #include "leased_store.h"
@@ -8,6 +9,36 @@
 using namespace std;
 extern bool verbose;
 
+// Human-readable name for a Boost.Geometry validity failure, used for diagnostics.
+static const char* validityFailureName(geom::validity_failure_type failure) {
+	switch (failure) {
+		case geom::no_failure:                        return "no_failure";
+		case geom::failure_few_points:                return "few_points";
+		case geom::failure_wrong_topological_dimension: return "wrong_topological_dimension";
+		case geom::failure_spikes:                    return "spikes";
+		case geom::failure_duplicate_points:          return "duplicate_points";
+		case geom::failure_not_closed:                return "not_closed";
+		case geom::failure_self_intersections:        return "self_intersections";
+		case geom::failure_wrong_orientation:         return "wrong_orientation";
+		case geom::failure_interior_rings_outside:    return "interior_rings_outside";
+		case geom::failure_nested_interior_rings:     return "nested_interior_rings";
+		case geom::failure_disconnected_interior:     return "disconnected_interior";
+		case geom::failure_intersecting_interiors:    return "intersecting_interiors";
+		case geom::failure_wrong_corner_order:        return "wrong_corner_order";
+		case geom::failure_invalid_coordinate:        return "invalid_coordinate";
+		default:                                      return "unknown";
+	}
+}
+
+// Thin wrapper around the shared repair_multi_polygon() that adds per-object
+// verbose diagnostics. See geom.cpp for the dissolve + zero-width buffer logic.
+static bool repairMultiPolygon(MultiPolygon &mp, NodeID objectID) {
+	bool ok = repair_multi_polygon(mp);
+	if (!ok && verbose)
+		std::cerr << ("multipolygon repair failed for object " + std::to_string(objectID) + "\n");
+	return ok;
+}
+
 thread_local LeasedStore<TileDataSource::point_store_t> pointStore;
 thread_local LeasedStore<TileDataSource::linestring_store_t> linestringStore;
 thread_local LeasedStore<TileDataSource::multi_linestring_store_t> multilinestringStore;
@@ -332,21 +363,39 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType,
 			geom::validity_failure_type failure = geom::validity_failure_type::no_failure;
 			bool valid = geom::is_valid(mp,failure);
 			if (!valid) {
+				if (verbose) {
+					// Build the whole line first and emit it with a single stream write:
+					// tilemaker runs multi-threaded and chained operator<< calls are not
+					// atomic, so per-token writes interleave into unreadable output.
+					std::ostringstream msg;
+					msg << "invalid multipolygon for object " << objectID
+					    << " at z" << bbox.zoom << " " << bbox.index.x << "/" << bbox.index.y
+					    << ": " << validityFailureName(failure) << "\n";
+					std::cerr << msg.str();
+				}
 				if (failure==geom::failure_spikes) {
 					geom::remove_spikes(mp);
 					failure = geom::validity_failure_type::no_failure;
 					valid = geom::is_valid(mp,failure);
 				}
 				if (!valid && (failure==geom::failure_self_intersections || failure==geom::failure_intersecting_interiors)) {
+					// fast_clip can introduce self-intersections; redo the clip with the
+					// slower but robust Boost intersection against the original geometry.
 					MultiPolygon output;
 					geom::intersection(input, box, output);
 					geom::correct(output);
 
-					// retry with Boost intersection if fast_clip has caused self-intersections
+					// The intersection result can itself still be invalid for very complex
+					// multipolygons (e.g. large reservoirs), which previously produced
+					// dropped or holey tiles. Repair it before returning.
+					repairMultiPolygon(output, objectID);
 					multiPolygonClipCache.add(bbox, objectID, output);
 					return output;
 				} else if (!valid) {
-					// occasionally also wrong_topological_dimension, disconnected_interior
+					// occasionally also wrong_topological_dimension, disconnected_interior:
+					// defects geom::correct cannot mend. Repair mp in place; on failure it
+					// is left unchanged so behaviour never regresses.
+					repairMultiPolygon(mp, objectID);
 				}
 			}
 

src/tile_worker.cpp

@@ -235,13 +235,31 @@ void writeMultiPolygon(
 	geom::correct(current);
 
 	geom::validity_failure_type failure;
-	if (verbose && !geom::is_valid(current, failure)) { 
-		cout << "output multipolygon has " << boost_validity_error(failure) << endl; 
-
-		if (!geom::is_valid(mp, failure)) 
-			cout << "input multipolygon has " << boost_validity_error(failure) << endl; 
-		else
-			cout << "input multipolygon valid" << endl;
+	if (!geom::is_valid(current, failure)) {
+		if (verbose) {
+			cout << "output multipolygon has " << boost_validity_error(failure) << endl;
+
+			if (!geom::is_valid(mp, failure))
+				cout << "input multipolygon has " << boost_validity_error(failure) << endl;
+			else
+				cout << "input multipolygon valid" << endl;
+		}
+		
+		if (simplifyLevel > 0) {
+			// Simplification can turn a valid input into a self-intersecting/spiky
+			// one; such polygons are silently dropped by many renderers (missing
+			// features). Repair (dissolve, then zero-width buffer) before writing.
+			bool repaired = repair_multi_polygon(current);
+
+			if (geom::is_empty(current))
+				return;
+
+			if (verbose && !repaired) {
+				geom::validity_failure_type postFailure;
+				if (!geom::is_valid(current, postFailure))
+					cout << "output multipolygon STILL invalid after repair: " << boost_validity_error(postFailure) << endl;
+			}
+		}
 	}
 
 	vtzero::polygon_feature_builder fbuilder{vtLayer};

src/visvalingam.cpp

@@ -255,11 +255,17 @@ Polygon simplifyVis(const Polygon &p, double max_distance) {
 	}
 	return output;
 }
-MultiPolygon simplifyVis(const MultiPolygon &mp, double max_distance) { 
+MultiPolygon simplifyVis(const MultiPolygon &mp, double max_distance) {
 	MultiPolygon output;
 	for (const auto &p : mp) {
 		output.emplace_back(simplifyVis(p, max_distance));
 	}
-	make_valid(output);
+	// Per-ring simplification can leave the multipolygon invalid. Use the
+	// area-preserving repair instead of a bare make_valid: on large/complex
+	// geometries an unconditional dissolve can collapse the polygon and drop
+	// almost all of its area (which showed up as missing lake tiles at low
+	// zoom). repair_multi_polygon keeps the simplified geometry untouched if a
+	// repair would not preserve the covered area.
+	repair_multi_polygon(output);
 	return output;
 }