From 4423a791fa59e42e8c05d9529c8d3b3382060bea Mon Sep 17 00:00:00 2001 From: "daniel.korn" Date: Thu, 29 Oct 2026 09:17:21 +0100 Subject: [PATCH 1/2] perf: exact rectangle clip for valid polygons, track multipolygon validity (planet -25% wall, -31% CPU) --- Makefile | 6 + include/geom.h | 10 +- include/osm_mem_tiles.h | 3 +- include/tile_data.h | 12 +- src/geom.cpp | 285 ++++++++++++++++++++++++++++++++++++- src/osm_lua_processing.cpp | 11 +- src/tile_data.cpp | 31 +++- test/geom.test.cpp | 73 ++++++++++ 8 files changed, 415 insertions(+), 16 deletions(-) create mode 100644 test/geom.test.cpp diff --git a/Makefile b/Makefile index b651ec64..57227d4c 100644 --- a/Makefile +++ b/Makefile @@ -157,6 +157,7 @@ test: \ test_attribute_store \ test_config_validator \ test_deque_map \ + test_geom \ test_helpers \ test_options_parser \ test_pbf_reader \ @@ -196,6 +197,11 @@ test_deque_map: \ test/deque_map.test.o $(CXX) $(CXXFLAGS) -o test.deque_map $^ $(INC) $(LIB) $(LDFLAGS) && ./test.deque_map +test_geom: \ + src/geom.o \ + test/geom.test.o + $(CXX) $(CXXFLAGS) -o test.geom $^ $(INC) $(LIB) $(LDFLAGS) && ./test.geom + test_helpers: \ src/helpers.o \ src/external/libdeflate/lib/adler32.o \ diff --git a/include/geom.h b/include/geom.h index a599da42..0a35df52 100644 --- a/include/geom.h +++ b/include/geom.h @@ -81,8 +81,10 @@ 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, bool knownInvalid = false); +// best-effort input so callers never regress. If wholeValid is given, it is set +// to true only if mp as a whole is known to pass geom::is_valid afterwards (a true +// result alone can leave individually valid polygons that overlap each other). +bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid = false, bool *wholeValid = nullptr); void union_many(std::vector &mps); @@ -91,6 +93,10 @@ char bit_code(Point const &p, Box const &bbox); // fast_clip returns true if the geometry was modified (false means identity). bool fast_clip(Ring &points, Box const &bbox); bool fast_clip(MultiPolygon &mp, Box const &bbox); +// Exact clip of a valid multipolygon to bbox that keeps it valid (fast_clip can +// leave zero-width bridges along the box edge). Returns false if it could not +// clip robustly; `out` is then unspecified and the caller must use geom::intersection. +bool rect_clip(const MultiPolygon &mp, Box const &bbox, MultiPolygon &out); // Cheap conservative validity tests, run before geom::is_valid. They return // true only when the geometry is certainly valid (is_valid would return true diff --git a/include/osm_mem_tiles.h b/include/osm_mem_tiles.h index e8d7ed9b..1966bc4b 100644 --- a/include/osm_mem_tiles.h +++ b/include/osm_mem_tiles.h @@ -61,7 +61,8 @@ class OsmMemTiles : public TileDataSource { // CorrectGeometry found them Valid (or Reversed: valid after reversing the ring) // on the same lon/latp coordinates. bool polygonKnownValid(NodeID objectID) const override { - return objectID >= OSM_THRESHOLD && (IS_WAY(objectID) || IS_WAY_REV(objectID)); + if (objectID < OSM_THRESHOLD) return TileDataSource::polygonKnownValid(objectID); + return IS_WAY(objectID) || IS_WAY_REV(objectID); } const NodeStore& nodeStore; diff --git a/include/tile_data.h b/include/tile_data.h index 5bc9fa67..3e8d20e9 100644 --- a/include/tile_data.h +++ b/include/tile_data.h @@ -298,6 +298,8 @@ class TileDataSource { std::vector linestringStores; std::vector multilinestringStores; std::vector multipolygonStores; + // Per multipolygonStores entry: stored geometry is known to pass geom::is_valid + std::vector> multipolygonValid; ClipCache multiPolygonClipCache; ClipCache multiLinestringClipCache; @@ -359,6 +361,7 @@ class TileDataSource { pointStores[0].push_back(Point(0,0)); linestringStores[0].push_back(linestring_t()); multipolygonStores[0].push_back(multi_polygon_t()); + multipolygonValid[0].push_back(false); multilinestringStores[0].push_back(multi_linestring_t()); } void reportSize() const; @@ -380,7 +383,11 @@ class TileDataSource { virtual void populateMultiPolygon(MultiPolygon& dst, NodeID objectID); // True if the polygon populateMultiPolygon returns for objectID is known to // have passed geom::is_valid on exactly the same coordinates. - virtual bool polygonKnownValid(NodeID objectID) const { return false; } + virtual bool polygonKnownValid(NodeID objectID) const { + const auto& valid = multipolygonValid[getShard(objectID)]; + const auto offset = getId(objectID); + return offset < valid.size() && valid[offset]; + } unsigned int getIndexZoom() const { return indexZoom; } inline size_t getId(NodeID id) const { @@ -415,7 +422,8 @@ class TileDataSource { return shard.at(offset); } - NodeID storeMultiPolygon(const MultiPolygon& src); + // knownValid: src has passed geom::is_valid (see polygonKnownValid) + NodeID storeMultiPolygon(const MultiPolygon& src, bool knownValid = false); multi_polygon_t const &retrieveMultiPolygon(NodeID id) const { const auto& shardId = getShard(id); diff --git a/src/geom.cpp b/src/geom.cpp index 50f54705..24c008f0 100644 --- a/src/geom.cpp +++ b/src/geom.cpp @@ -484,11 +484,15 @@ static bool repair_one_polygon(const Polygon &p, double minArea, MultiPolygon &o return false; } -bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid) +bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid, bool *wholeValid) { + if (wholeValid) *wholeValid = false; if (!knownInvalid) { const bool pre = multiPolygonDefinitelyValid(mp); - if (pre || geom::is_valid(mp)) return true; + if (pre || geom::is_valid(mp)) { + if (wholeValid) *wholeValid = true; + return true; + } } // Repair PER POLYGON, area-preserving. Running make_valid/buffer on the whole @@ -527,6 +531,8 @@ bool repair_multi_polygon(MultiPolygon &mp, bool knownInvalid) } mp = std::move(out); + // A single valid polygon is a valid multipolygon + if (wholeValid) *wholeValid = allValid && mp.size() == 1; return allValid; } @@ -699,3 +705,278 @@ bool fast_clip(MultiPolygon &mp, Box const &bbox) { mp.end()); return modified || mp.size() != sizeBefore; } + +// --------------- +// Exact clipping of a valid polygon to an axis-aligned box. +// +// Sutherland-Hodgman (fast_clip) is exact for convex input but leaves zero-width +// "bridges" along the box edge wherever a concave polygon leaves and re-enters +// the box; Boost reports those as self-intersections. This clipper instead cuts +// every ring into the chains that lie inside the box and reconnects them along +// the box boundary (Weiler-Atherton specialised to a rectangle), so a valid input +// gives a valid output. Whenever the input does not fit its assumptions (numerical +// trouble, a hole it cannot place) it gives up and returns false, and the caller +// falls back to Boost. + +namespace { + +struct RectClipper { + double x0, y0, x1, y1, w, h, perimeter; + + explicit RectClipper(const Box &box) : + x0(box.min_corner().x()), y0(box.min_corner().y()), + x1(box.max_corner().x()), y1(box.max_corner().y()), + w(x1 - x0), h(y1 - y0), perimeter(2 * (x1 - x0) + 2 * (y1 - y0)) {} + + bool inside(const Point &p) const { + return p.x() >= x0 && p.x() <= x1 && p.y() >= y0 && p.y() <= y1; + } + + // Position along the box boundary, clockwise from the top-left corner. + // Returns false if p is not exactly on the boundary. + bool boundaryPos(const Point &p, double &s) const { + if (p.y() == y1) s = p.x() - x0; + else if (p.x() == x1) s = w + (y1 - p.y()); + else if (p.y() == y0) s = w + h + (x1 - p.x()); + else if (p.x() == x0) s = 2 * w + h + (p.y() - y0); + else return false; + return true; + } + + // Liang-Barsky: the part [t0,t1] of segment a->b inside the closed box, with + // its end points snapped exactly onto the box edge they were clipped against. + bool clipSegment(const Point &a, const Point &b, double &t0, double &t1, Point &pa, Point &pb) const { + const double dx = b.x() - a.x(), dy = b.y() - a.y(); + const double p[4] = { -dx, dx, -dy, dy }; + const double q[4] = { a.x() - x0, x1 - a.x(), a.y() - y0, y1 - a.y() }; + int e0 = -1, e1 = -1; + t0 = 0; t1 = 1; + for (int k = 0; k < 4; k++) { + if (p[k] == 0) { + if (q[k] < 0) return false; + continue; + } + const double r = q[k] / p[k]; + if (p[k] < 0) { if (r > t0) { t0 = r; e0 = k; } } + else { if (r < t1) { t1 = r; e1 = k; } } + } + if (t0 > t1) return false; + pa = e0 < 0 ? a : onEdge(a, dx, dy, t0, e0); + pb = e1 < 0 ? b : onEdge(a, dx, dy, t1, e1); + return true; + } + + Point onEdge(const Point &a, double dx, double dy, double t, int edge) const { + switch (edge) { + case 0: return Point(x0, std::min(std::max(a.y() + t * dy, y0), y1)); + case 1: return Point(x1, std::min(std::max(a.y() + t * dy, y0), y1)); + case 2: return Point(std::min(std::max(a.x() + t * dx, x0), x1), y0); + default: return Point(std::min(std::max(a.x() + t * dx, x0), x1), y1); + } + } + + // Corners in clockwise boundary order, starting top-left + Point corner(int i) const { + switch (i) { + case 0: return Point(x0, y1); + case 1: return Point(x1, y1); + case 2: return Point(x1, y0); + default: return Point(x0, y0); + } + } + double cornerPos(int i) const { + switch (i) { + case 0: return 0; + case 1: return w; + case 2: return w + h; + default: return 2 * w + h; + } + } +}; + +enum class RingClip { Inside, Outside, Cut, Fail }; + +struct RectChain { + std::vector points; + double entry, exit; + bool used = false; +}; + +// Cut a closed ring into the chains inside the box. Inside: every vertex is +// inside (so is the ring). Outside: no part of the ring enters the box. +RingClip clipRing(const RectClipper &rc, const Ring &ring, std::vector &chains) { + if (ring.size() < 4) return RingClip::Fail; + const std::size_t n = ring.size() - 1; // closed: ring[n] == ring[0] + std::size_t start = n; + for (std::size_t i = 0; i < n; i++) + if (!rc.inside(ring[i])) { start = i; break; } + if (start == n) return RingClip::Inside; + + // Start at a vertex outside the box, so every chain begins with an entry. + std::vector current; + bool building = false, any = false; + auto finish = [&]() -> bool { + RectChain c; + if (current.size() < 2 || !rc.boundaryPos(current.front(), c.entry) || !rc.boundaryPos(current.back(), c.exit)) return false; + c.points.swap(current); + chains.push_back(std::move(c)); + current.clear(); + building = false; + any = true; + return true; + }; + for (std::size_t k = 0; k < n; k++) { + const Point &a = ring[(start + k) % n]; + const Point &b = ring[(start + k + 1) % n]; + double t0, t1; + Point pa, pb; + const bool hit = rc.clipSegment(a, b, t0, t1, pa, pb) && t1 > t0; + if (!hit || (building && t0 > 0)) { + if (building && !finish()) return RingClip::Fail; + if (!hit) continue; + } + if (!building) { + current.push_back(pa); + building = true; + } + if (!geom::equals(pb, current.back())) current.push_back(pb); + if (t1 < 1 && !finish()) return RingClip::Fail; + } + // The walk ends at the outside start vertex, so the last chain must have exited + if (building) return RingClip::Fail; + return any ? RingClip::Cut : RingClip::Outside; +} + +// Crossing-number test: 1 inside, -1 outside, 0 on the boundary. +int pointInRing(const Point &p, const Ring &ring) { + bool in = false; + for (std::size_t i = 0, j = ring.size() - 1; i < ring.size(); j = i++) { + const Point &a = ring[i], &b = ring[j]; + if (geom::equals(a, p)) return 0; + if ((a.y() > p.y()) != (b.y() > p.y())) { + const double x = a.x() + (p.y() - a.y()) * (b.x() - a.x()) / (b.y() - a.y()); + if (x == p.x()) return 0; + if (x > p.x()) in = !in; + } + } + return in ? 1 : -1; +} + +// Where a ring that does not cross `outer` lies relative to it (0: unsure). +int ringInRing(const Ring &ring, const Ring &outer) { + for (const Point &p : ring) { + const int r = pointInRing(p, outer); + if (r != 0) return r; + } + return 0; +} + +bool clipPolygon(const RectClipper &rc, const Polygon &poly, MultiPolygon &out) { + std::vector chains; + const RingClip outerClip = clipRing(rc, poly.outer(), chains); + if (outerClip == RingClip::Fail) return false; + if (outerClip == RingClip::Inside) { + out.push_back(poly); + return true; + } + + std::vector insideHoles, outsideHoles; + for (const Ring &inner : poly.inners()) { + switch (clipRing(rc, inner, chains)) { + case RingClip::Fail: return false; + case RingClip::Inside: insideHoles.push_back(&inner); break; + case RingClip::Outside: outsideHoles.push_back(&inner); break; + case RingClip::Cut: break; + } + } + + const std::size_t firstOuter = out.size(); + if (chains.empty()) { + // No ring crosses the box: it is either inside the polygon or disjoint from it + const Point centre((rc.x0 + rc.x1) / 2, (rc.y0 + rc.y1) / 2); + const int inOuter = pointInRing(centre, poly.outer()); + if (inOuter == 0) return false; + bool inside = inOuter > 0; + for (const Ring *hole : outsideHoles) { + const int inHole = pointInRing(centre, *hole); + if (inHole == 0) return false; + if (inHole > 0) inside = false; + } + if (!inside) return insideHoles.empty(); + Polygon boxPoly; + for (int i = 0; i < 4; i++) boxPoly.outer().push_back(rc.corner(i)); + boxPoly.outer().push_back(rc.corner(0)); + out.push_back(std::move(boxPoly)); + } else { + // Polygon interiors are on the right of every ring (clockwise outers, + // anticlockwise holes), so from each exit the boundary walk goes clockwise + // to the next entry. + std::vector> entries; + entries.reserve(chains.size()); + for (std::size_t i = 0; i < chains.size(); i++) entries.emplace_back(chains[i].entry, i); + std::sort(entries.begin(), entries.end()); + + for (std::size_t first = 0; first < chains.size(); first++) { + if (chains[first].used) continue; + Ring ring; + std::size_t c = first; + while (true) { + RectChain &chain = chains[c]; + chain.used = true; + for (const Point &p : chain.points) + if (ring.empty() || !geom::equals(p, ring.back())) ring.push_back(p); + + auto it = std::lower_bound(entries.begin(), entries.end(), std::make_pair(chain.exit, std::size_t(0))); + if (it == entries.end()) it = entries.begin(); + const std::size_t next = it->second; + if (chains[next].used && next != first) return false; + + double span = chains[next].entry - chain.exit; + if (span < 0) span += rc.perimeter; + std::pair passed[4]; + int numPassed = 0; + for (int k = 0; k < 4; k++) { + double d = rc.cornerPos(k) - chain.exit; + if (d < 0) d += rc.perimeter; + if (d > 0 && d < span) passed[numPassed++] = { d, k }; + } + std::sort(passed, passed + numPassed); + for (int k = 0; k < numPassed; k++) ring.push_back(rc.corner(passed[k].second)); + + if (next == first) break; + c = next; + } + if (!geom::equals(ring.front(), ring.back())) ring.push_back(ring.front()); + if (ring.size() < 4 || geom::area(ring) <= 0) return false; + Polygon p; + p.outer() = std::move(ring); + out.push_back(std::move(p)); + } + } + + // Holes wholly inside the box go to whichever new outer contains them + for (const Ring *hole : insideHoles) { + bool placed = false; + for (std::size_t i = firstOuter; i < out.size() && !placed; i++) { + const int r = ringInRing(*hole, out[i].outer()); + if (r == 0) return false; + if (r > 0) { + out[i].inners().push_back(*hole); + placed = true; + } + } + if (!placed) return false; + } + return true; +} + +} // namespace + +bool rect_clip(const MultiPolygon &mp, Box const &bbox, MultiPolygon &out) { + out.clear(); + const RectClipper rc(bbox); + if (!(rc.w > 0 && rc.h > 0)) return false; + for (const Polygon &poly : mp) + if (!clipPolygon(rc, poly, out)) return false; + return true; +} diff --git a/src/osm_lua_processing.cpp b/src/osm_lua_processing.cpp index c719eaac..ba34304e 100644 --- a/src/osm_lua_processing.cpp +++ b/src/osm_lua_processing.cpp @@ -864,8 +864,13 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) { if (isRelation) { try { MultiPolygon mp = multiPolygonCached(); - if(CorrectGeometry(mp) == CorrectGeometryResult::Invalid) return; - NodeID id = osmMemTiles.storeMultiPolygon(mp); + auto correctionResult = CorrectGeometry(mp); + if(correctionResult == CorrectGeometryResult::Invalid) return; + // make_valid usually succeeds; checking that once here saves + // re-validating the relation in every tile it is clipped to. + const bool knownValid = correctionResult == CorrectGeometryResult::Valid || + PrecheckValid(mp) || geom::is_valid(mp); + NodeID id = osmMemTiles.storeMultiPolygon(mp, knownValid); OutputObject oo(geomType, layers.layerMap[layerName], id, 0, layerMinZoom); outputs.push_back(std::make_pair(std::move(oo), attributes)); } catch(std::out_of_range &err) { @@ -893,7 +898,7 @@ void OsmLuaProcessing::Layer(const string &layerName, bool area) { id = USE_WAY_STORE_REVERSED | originalOsmID; wayEmitted = true; } else - id = osmMemTiles.storeMultiPolygon(mp); + id = osmMemTiles.storeMultiPolygon(mp, correctionResult != CorrectGeometryResult::Corrected); OutputObject oo(geomType, layers.layerMap[layerName], id, 0, layerMinZoom); outputs.push_back(std::make_pair(std::move(oo), attributes)); } diff --git a/src/tile_data.cpp b/src/tile_data.cpp index 3fff009a..89b367e0 100644 --- a/src/tile_data.cpp +++ b/src/tile_data.cpp @@ -34,8 +34,8 @@ static const char* validityFailureName(geom::validity_failure_type failure) { // 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); +static bool repairMultiPolygon(MultiPolygon &mp, NodeID objectID, bool *wholeValid = nullptr) { + bool ok = repair_multi_polygon(mp, false, wholeValid); if (!ok && verbose) std::cerr << ("multipolygon repair failed for object " + std::to_string(objectID) + "\n"); return ok; @@ -60,6 +60,7 @@ TileDataSource::TileDataSource(size_t threadNum, unsigned int indexZoom, bool in linestringStores(threadNum), multilinestringStores(threadNum), multipolygonStores(threadNum), + multipolygonValid(threadNum), multiPolygonClipCache(ClipCache(threadNum, indexZoom)), multiLinestringClipCache(ClipCache(threadNum, indexZoom)) { @@ -398,6 +399,20 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType, geom::correct(mp); geom::validity_failure_type failure = geom::validity_failure_type::no_failure; bool valid = multiPolygonDefinitelyValid(mp); + if (!valid && inputKnownValid) { + // fast_clip leaves zero-width bridges along the box edge where a concave + // polygon leaves and re-enters the box. Clip the valid input exactly + // instead, which is much cheaper than the is_valid/intersection/repair + // path below; rect_clip bails out on anything it can't handle. + MultiPolygon original; + if (cachedClip == nullptr) populateMultiPolygon(original, objectID); + MultiPolygon clipped; + if (rect_clip(cachedClip == nullptr ? original : input, box, clipped) && + (multiPolygonDefinitelyValid(clipped) || geom::is_valid(clipped))) { + multiPolygonClipCache.add(bbox, objectID, clipped, true); + return clipped; + } + } if (!valid) valid = geom::is_valid(mp,failure); if (!valid) { @@ -432,14 +447,16 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType, // 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); + bool outputValid; + repairMultiPolygon(output, objectID, &outputValid); + // Children of this tile can then take the fast paths for valid input + multiPolygonClipCache.add(bbox, objectID, output, outputValid); return output; } else if (!valid) { // 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); + repairMultiPolygon(mp, objectID, &valid); } } @@ -658,7 +675,7 @@ NodeID TileDataSource::storeLinestring(const Linestring& src) { return rv; } -NodeID TileDataSource::storeMultiPolygon(const MultiPolygon& src) { +NodeID TileDataSource::storeMultiPolygon(const MultiPolygon& src, bool knownValid) { const auto& store = multipolygonStore.get(this); multi_polygon_t dst; @@ -676,6 +693,8 @@ NodeID TileDataSource::storeMultiPolygon(const MultiPolygon& src) { NodeID offset = store.second->size(); store.second->emplace_back(std::move(dst)); + // The store is leased to this thread, so its validity flags are too + multipolygonValid[store.first].push_back(knownValid); NodeID rv = (store.first << (TILE_DATA_ID_SIZE - shardBits)) + offset; return rv; } diff --git a/test/geom.test.cpp b/test/geom.test.cpp new file mode 100644 index 00000000..cdabafa1 --- /dev/null +++ b/test/geom.test.cpp @@ -0,0 +1,73 @@ +#include +#include +#include "external/minunit.h" +#include "geom.h" + +static MultiPolygon fromWkt(const char *wkt) { + MultiPolygon mp; + geom::read_wkt(wkt, mp); + geom::correct(mp); + return mp; +} + +// rect_clip must give a valid result of the expected area +static void checkClip(const char *wkt, const Box &box, std::size_t expectedPolygons, double expectedArea) { + MultiPolygon mp = fromWkt(wkt); + mu_check(geom::is_valid(mp)); + + MultiPolygon clipped; + mu_check(rect_clip(mp, box, clipped)); + mu_check(geom::is_valid(clipped)); + mu_check(clipped.size() == expectedPolygons); + mu_check(std::abs(geom::area(clipped) - expectedArea) < 1e-9); +} + +MU_TEST(test_rect_clip_concave) { + // A U shape whose arms leave the box through its top: Sutherland-Hodgman + // joins the arms with a zero-width bridge along the box edge + const char *u = "MULTIPOLYGON(((0 0,0 10,2 10,2 2,8 2,8 10,10 10,10 0,0 0)))"; + checkClip(u, Box(Point(-1, -1), Point(11, 5)), 1, 32); + // Cutting the base off leaves the two arms as separate polygons + checkClip(u, Box(Point(-1, 4), Point(11, 12)), 2, 24); +} + +MU_TEST(test_rect_clip_holes) { + const char *holes = "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0),(1 1,3 1,3 3,1 3,1 1),(6 4,9 4,9 6,6 6,6 4)))"; + // One hole inside the box, one cut by its right edge + checkClip(holes, Box(Point(-1, -1), Point(7, 7)), 1, 43); + // Box wholly inside the polygon: the box itself + checkClip(holes, Box(Point(4, 7), Point(5, 8)), 1, 1); + // Box wholly inside a hole: nothing + checkClip(holes, Box(Point(1.5, 1.5), Point(2.5, 2.5)), 0, 0); +} + +MU_TEST(test_rect_clip_trivial) { + const char *square = "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))"; + // Disjoint + checkClip(square, Box(Point(20, 20), Point(30, 30)), 0, 0); + // Polygon wholly inside the box is returned as it is + MultiPolygon mp = fromWkt(square), clipped; + mu_check(rect_clip(mp, Box(Point(-1, -1), Point(11, 11)), clipped)); + mu_check(geom::equals(mp, clipped)); +} + +MU_TEST(test_rect_clip_vertex_on_edge) { + // The ring crosses the box edge at a (duplicated) vertex exactly on it. + // (Boost's intersection returns nothing for the first box.) + const char *poly = "MULTIPOLYGON(((0 0,0 4,0 4,0 10,10 10,10 6,12 4,10 2,10 0,0 0)))"; + checkClip(poly, Box(Point(-1, 4), Point(11, 11)), 1, 61.5); + checkClip(poly, Box(Point(-1, -1), Point(11, 11)), 1, 103); +} + +MU_TEST_SUITE(test_suite_geom) { + MU_RUN_TEST(test_rect_clip_concave); + MU_RUN_TEST(test_rect_clip_holes); + MU_RUN_TEST(test_rect_clip_trivial); + MU_RUN_TEST(test_rect_clip_vertex_on_edge); +} + +int main() { + MU_RUN_SUITE(test_suite_geom); + MU_REPORT(); + return MU_EXIT_CODE; +} From 5afb9d9e59239df68b477ded8d5d9524e754bba6 Mon Sep 17 00:00:00 2001 From: "daniel.korn" Date: Thu, 29 Oct 2026 14:23:37 +0100 Subject: [PATCH 2/2] fix(rect_clip): split rings that touch the clip box from inside A ring can touch the box edge at a single vertex, with the polygon on both sides of that vertex along the edge. rect_clip kept it as one chain, so the walk along the box edge ran through the vertex and the result was a ring touching itself. Boost rejects that as invalid, and the Boost intersection fallback then lost most or all of the polygon. Split the chain at such a vertex, so it is both an exit and an entry and the result is two polygons meeting there. Touches where the polygon lies between the two arms (a wedge) are left as they were. Co-Authored-By: Claude Opus 5.5 (1M context) --- src/geom.cpp | 51 +++++++++++++++++++++++++++++++++++++++++++--- test/geom.test.cpp | 11 ++++++++++ 2 files changed, 59 insertions(+), 3 deletions(-) diff --git a/src/geom.cpp b/src/geom.cpp index 24c008f0..867a84f9 100644 --- a/src/geom.cpp +++ b/src/geom.cpp @@ -731,6 +731,32 @@ struct RectClipper { bool inside(const Point &p) const { return p.x() >= x0 && p.x() <= x1 && p.y() >= y0 && p.y() <= y1; } + bool strictlyInside(const Point &p) const { + return p.x() > x0 && p.x() < x1 && p.y() > y0 && p.y() < y1; + } + // Segment a->b, both inside the box, runs through its interior (not along an edge) + bool throughInterior(const Point &a, const Point &b) const { + return strictlyInside(Point((a.x() + b.x()) / 2, (a.y() + b.y()) / 2)); + } + // Segment b->c, b on the boundary, starts off into the interior + bool entersInterior(const Point &b, const Point &c) const { + double t0, t1; + Point pa, pb; + return clipSegment(b, c, t0, t1, pa, pb) && t0 == 0 && t1 > 0 && throughInterior(b, pb); + } + // A ring p->t->q touches the boundary at t (not a corner) from inside. True if the + // polygon interior (right of the ring) covers the boundary on both sides of t, + // i.e. the edge direction lies in the clockwise sweep from t->q to t->p. + bool pinchAt(const Point &p, const Point &t, const Point &q) const { + if ((t.x() == x0 || t.x() == x1) && (t.y() == y0 || t.y() == y1)) return false; + const double ex = (t.y() == y0 || t.y() == y1) ? 1 : 0, ey = 1 - ex; + const double qx = q.x() - t.x(), qy = q.y() - t.y(); + auto clockwiseFromQ = [&](double vx, double vy) { + const double a = -std::atan2(qx * vy - qy * vx, qx * vx + qy * vy); + return a < 0 ? a + 2 * M_PI : a; + }; + return clockwiseFromQ(ex, ey) < clockwiseFromQ(p.x() - t.x(), p.y() - t.y()); + } // Position along the box boundary, clockwise from the top-left corner. // Returns false if p is not exactly on the boundary. @@ -800,6 +826,7 @@ struct RectChain { std::vector points; double entry, exit; bool used = false; + bool touchExit = false; // ends at a vertex touching the boundary from inside; the next chain starts there }; // Cut a closed ring into the chains inside the box. Inside: every vertex is @@ -840,7 +867,21 @@ RingClip clipRing(const RectClipper &rc, const Ring &ring, std::vector= 2 && !rc.strictlyInside(b) && rc.throughInterior(current[current.size() - 2], b) && + rc.entersInterior(b, ring[(start + k + 2) % n]) && + rc.pinchAt(current[current.size() - 2], b, ring[(start + k + 2) % n])) { + // b touches the boundary from inside (reached through the interior, left + // into the interior) and the polygon lies along the boundary on both sides + // of it. Kept in one chain, the boundary walk would run through b and make + // the ring touch itself there: split the chain at b, so b is an exit and + // an entry, with the polygon pinched in two there. + if (!finish()) return RingClip::Fail; + chains.back().touchExit = true; + current.push_back(b); + building = true; + } } // The walk ends at the outside start vertex, so the last chain must have exited if (building) return RingClip::Fail; @@ -926,13 +967,17 @@ bool clipPolygon(const RectClipper &rc, const Polygon &poly, MultiPolygon &out) for (const Point &p : chain.points) if (ring.empty() || !geom::equals(p, ring.back())) ring.push_back(p); - auto it = std::lower_bound(entries.begin(), entries.end(), std::make_pair(chain.exit, std::size_t(0))); + // From a touch exit, the entry at the same point is the far side of the + // touch: walk on to the next one + auto it = chain.touchExit + ? std::upper_bound(entries.begin(), entries.end(), std::make_pair(chain.exit, std::numeric_limits::max())) + : std::lower_bound(entries.begin(), entries.end(), std::make_pair(chain.exit, std::size_t(0))); if (it == entries.end()) it = entries.begin(); const std::size_t next = it->second; if (chains[next].used && next != first) return false; double span = chains[next].entry - chain.exit; - if (span < 0) span += rc.perimeter; + if (span < 0 || (span == 0 && chain.touchExit)) span += rc.perimeter; std::pair passed[4]; int numPassed = 0; for (int k = 0; k < 4; k++) { diff --git a/test/geom.test.cpp b/test/geom.test.cpp index cdabafa1..8464f04a 100644 --- a/test/geom.test.cpp +++ b/test/geom.test.cpp @@ -59,11 +59,22 @@ MU_TEST(test_rect_clip_vertex_on_edge) { checkClip(poly, Box(Point(-1, -1), Point(11, 11)), 1, 103); } +MU_TEST(test_rect_clip_touch_from_inside) { + // A notch from the top comes down to touch the box's bottom edge at (5 0), then + // goes back up. The box walk must not run through the touch: two polygons. + const char *notch = "MULTIPOLYGON(((-5 -5,-5 8,4 8,5 0,6 8,15 8,15 -5,-5 -5)))"; + checkClip(notch, Box(Point(0, 0), Point(10, 10)), 2, 72); + // The same touching the left edge + const char *sideNotch = "MULTIPOLYGON(((-5 -5,-5 15,8 15,8 6,0 5,8 4,8 -5,-5 -5)))"; + checkClip(sideNotch, Box(Point(0, 0), Point(10, 10)), 2, 72); +} + MU_TEST_SUITE(test_suite_geom) { MU_RUN_TEST(test_rect_clip_concave); MU_RUN_TEST(test_rect_clip_holes); MU_RUN_TEST(test_rect_clip_trivial); MU_RUN_TEST(test_rect_clip_vertex_on_edge); + MU_RUN_TEST(test_rect_clip_touch_from_inside); } int main() {