diff --git a/include/geom.h b/include/geom.h index d99acc1e..0897b63c 100644 --- a/include/geom.h +++ b/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 &mps); Point intersect_edge(Point const &a, Point const &b, char edge, Box const &bbox); diff --git a/src/geom.cpp b/src/geom.cpp index 69633372..07d70c41 100644 --- a/src/geom.cpp +++ b/src/geom.cpp @@ -3,6 +3,8 @@ #include #include +#include +#include #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 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 diff --git a/src/tile_data.cpp b/src/tile_data.cpp index 10c2f11b..800e774a 100644 --- a/src/tile_data.cpp +++ b/src/tile_data.cpp @@ -1,5 +1,6 @@ #include #include +#include #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 pointStore; thread_local LeasedStore linestringStore; thread_local LeasedStore 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); } } diff --git a/src/tile_worker.cpp b/src/tile_worker.cpp index 7e894332..1375410f 100644 --- a/src/tile_worker.cpp +++ b/src/tile_worker.cpp @@ -235,13 +235,32 @@ 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(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 (!geom::is_valid(mp, failure)) - cout << "input multipolygon has " << boost_validity_error(failure) << endl; - else - cout << "input multipolygon valid" << endl; + // Simplification (and the subsequent spike removal) can turn a valid + // input polygon into a self-intersecting or spiky one. Such invalid + // polygons are silently dropped by many vector-tile renderers, which + // shows up as missing features in individual tiles (e.g. holes in a + // lake at low zoom). Repair (dissolve, then zero-width buffer) before + // writing so only valid geometry is emitted. + 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}; diff --git a/src/visvalingam.cpp b/src/visvalingam.cpp index 7ae43e8f..b393345d 100644 --- a/src/visvalingam.cpp +++ b/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; }