From ffd911dd639bd541f140baa7d68eddec8d6b1a64 Mon Sep 17 00:00:00 2001 From: "daniel.korn" Date: Sat, 31 Oct 2026 18:18:34 +0100 Subject: [PATCH] perf(precheck): locate holes with a band index, validate shapefile polygons once The validity precheck located every hole by walking the whole outer ring and tested hole nesting pairwise, both capped at 2^26 vertex visits. Huge ocean polygons with thousands of islands exceeded the caps and fell back to boost::geometry::is_valid, twice (processShapeGeometry and again in ShpMemTiles::StoreGeometry), which made the planet shapefile phase last as long as its slowest polygon (91 s). - holesClearlyPlaced: bucket all ring segments into horizontal bands and cast one ray per hole, counting crossing parity per ring; no work cap - rayClearlyCrosses: exact half-open y comparisons, so rays across large polygons are no longer "unsure" at nearby vertices - shapefile polygons that are valid and inside the clip box skip the clip and pass knownValid to StoreGeometry Planet shapefile phase 95 s -> 21 s (-78%), total wall -3.9%; precheck re-verified against is_valid on planet shapefiles and Germany with 0 mismatches. --- include/shp_mem_tiles.h | 3 +- src/geom.cpp | 139 +++++++++++++++++++++++++++++++--------- src/shp_mem_tiles.cpp | 7 +- src/shp_processor.cpp | 20 ++++-- test/geom.test.cpp | 39 +++++++++++ 5 files changed, 166 insertions(+), 42 deletions(-) diff --git a/include/shp_mem_tiles.h b/include/shp_mem_tiles.h index a5cf4f3f..0bab561c 100644 --- a/include/shp_mem_tiles.h +++ b/include/shp_mem_tiles.h @@ -27,7 +27,8 @@ class ShpMemTiles : public TileDataSource const std::string& name, uint minzoom, int32_t score, - AttributeIndex attrIdx + AttributeIndex attrIdx, + bool knownValid = false // polygon: caller has checked it is valid ); std::vector QueryMatchingGeometries( diff --git a/src/geom.cpp b/src/geom.cpp index ed290645..b2998c68 100644 --- a/src/geom.cpp +++ b/src/geom.cpp @@ -177,13 +177,17 @@ struct PrecheckBox { }; struct PrecheckScratch { std::vector segs; - std::vector boxes, holeBoxes; + std::vector boxes; std::vector> pairs; + std::vector bandStart, bandFill, bandSegs, crossed; void trim() { if (segs.capacity() > kPrecheckRetained) std::vector().swap(segs); if (boxes.capacity() > kPrecheckRetained) std::vector().swap(boxes); - if (holeBoxes.capacity() > kPrecheckRetained) std::vector().swap(holeBoxes); if (pairs.capacity() > kPrecheckRetained) std::vector>().swap(pairs); + if (bandStart.capacity() > kPrecheckRetained) std::vector().swap(bandStart); + if (bandFill.capacity() > kPrecheckRetained) std::vector().swap(bandFill); + if (bandSegs.capacity() > kPrecheckRetained) std::vector().swap(bandSegs); + if (crossed.capacity() > kPrecheckRetained) std::vector().swap(crossed); } }; thread_local PrecheckScratch precheckScratch; @@ -330,20 +334,27 @@ bool segmentsClear(std::vector &segs, double tol) { return true; } +// Whether the ray from h to +x crosses segment a->b: 1 = yes, -1 = no, 0 = unsure. +// Vertices are classed as above (y > hy) or not with exact comparisons (the +// half-open rule), so a ray through a vertex is counted consistently and the +// crossing parity is exact; only the side of h is decided with a tolerance. +inline int rayClearlyCrosses(const Point &h, const Point &a, const Point &b, double tol) { + const double hx = h.x(), hy = h.y(); + if ((a.y() > hy) == (b.y() > hy)) return -1; + if (std::max(a.x(), b.x()) < hx) return -1; + if (std::min(a.x(), b.x()) > hx) return 1; + const int o = clearOrient(a, b, h, tol); + if (o == 0) return 0; + return (b.y() > a.y()) == (o > 0) ? 1 : -1; +} + // Crossing-number test: 1 = h clearly inside r, -1 = clearly outside, 0 = unsure. int pointClearlyInRing(const Point &h, const Ring &r, double tol) { bool inside = false; - const double hx = h.x(), hy = h.y(); for (std::size_t i = 0; i + 1 < r.size(); i++) { - const Point &a = r[i], &b = r[i + 1]; - if ((a.y() > hy + tol && b.y() > hy + tol) || (a.y() < hy - tol && b.y() < hy - tol)) continue; - if (std::fabs(a.y() - hy) <= tol || std::fabs(b.y() - hy) <= tol) return 0; - // a and b are clearly on opposite sides of the horizontal line through h - if (std::max(a.x(), b.x()) < hx - tol) continue; - if (std::min(a.x(), b.x()) > hx + tol) { inside = !inside; continue; } - const int o = clearOrient(a, b, h, tol); - if (o == 0) return 0; - if ((b.y() > a.y()) == (o > 0)) inside = !inside; // the ray to +x crosses a->b + const int c = rayClearlyCrosses(h, r[i], r[i + 1], tol); + if (c == 0) return 0; + if (c > 0) inside = !inside; } return inside ? 1 : -1; } @@ -359,12 +370,91 @@ int ringClearlyInRing(const Ring &inner, const Ring &r, double tol) { return 0; } +// Holes lie inside the outer ring (ring 0) and outside each other. `segs` holds the +// segments of all rings of p, known not to touch. They are bucketed into horizontal +// bands, so that a ray from a vertex of each hole is tested against the segments of +// that vertex's band only, counting crossings per ring (odd = inside that ring). +bool holesClearlyPlaced(const Polygon &p, PrecheckScratch &scratch, double tol) { + const auto &inners = p.inners(); + const std::vector &segs = scratch.segs; + if (inners.empty()) return true; + if (segs.size() >= std::numeric_limits::max()) return false; + + double minY = std::numeric_limits::infinity(), maxY = -minY; + for (const PrecheckSeg &s : segs) { minY = std::min(minY, s.minY); maxY = std::max(maxY, s.maxY); } + minY -= tol; maxY += tol; + // About 4 segments per band; fewer bands while tall segments would make the index much bigger than segs + std::size_t bands = std::max(1, segs.size() / 4); + double scale = 0; + auto bandOf = [&](double y) { + const double b = (y - minY) * scale; + return b <= 0 ? std::size_t(0) : std::min(bands - 1, static_cast(b)); + }; + std::size_t entries; + for (;;) { + scale = bands / (maxY - minY); + entries = 0; + for (const PrecheckSeg &s : segs) entries += bandOf(s.maxY + tol) - bandOf(s.minY - tol) + 1; + if (entries <= 4 * segs.size() || bands == 1) break; + bands = (bands + 1) / 2; + } + if (entries >= std::numeric_limits::max()) return false; + + // Segments of band b: bandSegs[bandStart[b] .. bandStart[b+1]) + std::vector &start = scratch.bandStart, &fill = scratch.bandFill, &items = scratch.bandSegs; + start.assign(bands + 1, 0); + for (const PrecheckSeg &s : segs) + for (std::size_t b = bandOf(s.minY - tol), hi = bandOf(s.maxY + tol); b <= hi; b++) start[b + 1]++; + for (std::size_t b = 0; b < bands; b++) start[b + 1] += start[b]; + fill.assign(start.begin(), start.end() - 1); + items.resize(entries); + for (std::size_t k = 0; k < segs.size(); k++) + for (std::size_t b = bandOf(segs[k].minY - tol), hi = bandOf(segs[k].maxY + tol); b <= hi; b++) + items[fill[b]++] = static_cast(k); + + // 1 = h inside the outer ring and no other hole, -1 = not, 0 = unsure + std::vector &crossed = scratch.crossed; + auto locate = [&](const Point &h, uint32_t ownRing) { + crossed.clear(); + const std::size_t b = bandOf(h.y()); + for (uint32_t k = start[b]; k < start[b + 1]; k++) { + const PrecheckSeg &s = segs[items[k]]; + if (s.ring == ownRing) continue; + const int c = rayClearlyCrosses(h, s.a[0], s.a[1], tol); + if (c == 0) return 0; + if (c > 0) crossed.push_back(s.ring); + } + std::sort(crossed.begin(), crossed.end()); + bool inOuter = false; + for (std::size_t i = 0; i < crossed.size(); ) { + std::size_t j = i; + while (j < crossed.size() && crossed[j] == crossed[i]) j++; + if ((j - i) % 2 == 1) { + if (crossed[i] != 0) return -1; // inside another hole + inOuter = true; + } + i = j; + } + return inOuter ? 1 : -1; + }; + + // A ring known not to touch the others lies wherever any of its vertices does + for (std::size_t i = 0; i < inners.size(); i++) { + const Ring &hole = inners[i]; + const std::size_t m = hole.size() - 1; + const std::size_t tries = std::min(m, 8); + int loc = 0; + for (std::size_t k = 0; k < tries && loc == 0; k++) + loc = locate(hole[k * m / tries], static_cast(i + 1)); + if (loc != 1) return false; + } + return true; +} + bool polygonClear(const Polygon &p, PrecheckScratch &scratch, PrecheckBox &outerBox) { const Ring &outer = p.outer(); const auto &inners = p.inners(); if (inners.size() >= std::numeric_limits::max()) return false; - // Every hole is located by walking the outer ring - if (inners.size() * outer.size() > kPrecheckMaxWork) return false; double maxAbs = 0; if (!ringFinite(outer, maxAbs)) return false; for (const auto &inner : inners) @@ -373,29 +463,16 @@ bool polygonClear(const Polygon &p, PrecheckScratch &scratch, PrecheckBox &outer std::vector &segs = scratch.segs; segs.clear(); - std::vector &holeBoxes = scratch.holeBoxes; - holeBoxes.resize(inners.size()); if (!ringClear(outer, false, tol, 0, segs, outerBox)) return false; - for (std::size_t i = 0; i < inners.size(); i++) { - if (!ringClear(inners[i], true, tol, static_cast(i + 1), segs, holeBoxes[i])) return false; - holeBoxes[i].id = static_cast(i); - } + PrecheckBox holeBox; + for (std::size_t i = 0; i < inners.size(); i++) + if (!ringClear(inners[i], true, tol, static_cast(i + 1), segs, holeBox)) return false; // No two rings touch, and no ring touches itself. if (!segmentsClear(segs, tol)) return false; // Holes lie inside the outer ring and are not nested in each other. - for (std::size_t i = 0; i < inners.size(); i++) { - if (boxesApart(holeBoxes[i], outerBox, tol)) return false; - if (ringClearlyInRing(inners[i], outer, tol) != 1) return false; - } - std::size_t work = 0; - return forOverlappingBoxes(holeBoxes, tol, [&](uint32_t i, uint32_t j) { - work += inners[i].size() + inners[j].size(); - return work <= kPrecheckMaxWork && - ringClearlyInRing(inners[i], inners[j], tol) == -1 && - ringClearlyInRing(inners[j], inners[i], tol) == -1; - }); + return holesClearlyPlaced(p, scratch, tol); } // a's outer ring, known not to touch any ring of b, lies outside b's interior: diff --git a/src/shp_mem_tiles.cpp b/src/shp_mem_tiles.cpp index e5c69e52..fc9bee81 100644 --- a/src/shp_mem_tiles.cpp +++ b/src/shp_mem_tiles.cpp @@ -137,7 +137,8 @@ void ShpMemTiles::StoreGeometry( const std::string& name, uint minzoom, int32_t score, - AttributeIndex attrIdx + AttributeIndex attrIdx, + bool knownValid ) { geom::model::box box; @@ -175,9 +176,9 @@ void ShpMemTiles::StoreGeometry( case POLYGON_: { - // Checked once here, so tiles can clip valid polygons on the fast path + // Checked once here (unless the caller has), so tiles can clip valid polygons on the fast path const MultiPolygon &mp = boost::get(geometry); - NodeID oid = storeMultiPolygon(mp, multiPolygonDefinitelyValid(mp) || geom::is_valid(mp)); + NodeID oid = storeMultiPolygon(mp, knownValid || multiPolygonDefinitelyValid(mp) || geom::is_valid(mp)); oo = std::make_shared(geomType, layerNum, oid, attrIdx, minzoom); std::vector oolist { *oo }; addGeometryToIndex(boost::get(geometry), oolist, 0); diff --git a/src/shp_processor.cpp b/src/shp_processor.cpp index f3b3c4f2..b32291b2 100644 --- a/src/shp_processor.cpp +++ b/src/shp_processor.cpp @@ -366,8 +366,8 @@ void ShpProcessor::processShapeGeometry(SHPObject* shape, AttributeIndex attrIdx // Make valid if needs be string reason; - const bool precheckValid = multiPolygonDefinitelyValid(multi); - if (!precheckValid && !geom::is_valid(multi, reason)) { + const bool valid = multiPolygonDefinitelyValid(multi) || geom::is_valid(multi, reason); + if (!valid) { if (verbose) cerr << "Shapefile entity " << shape->nShapeId << " type " << shapeType << " is invalid. Parts:" << shape->nParts << ". Reason:" << reason; make_valid(multi); if (verbose) { @@ -376,11 +376,17 @@ void ShpProcessor::processShapeGeometry(SHPObject* shape, AttributeIndex attrIdx cerr << endl; } } - // clip to bounding box - MultiPolygon out; - geom::intersection(multi, clippingBox, out); - if (boost::size(out)>0) { - shpMemTiles.StoreGeometry(layerNum, layer.name, POLYGON_, out, layer.indexed, hasName, name, minzoom, score, attrIdx); + // clip to bounding box (a valid polygon wholly inside it is stored as it is) + Box envelope; + geom::envelope(multi, envelope); + if (valid && geom::covered_by(envelope, clippingBox)) { + shpMemTiles.StoreGeometry(layerNum, layer.name, POLYGON_, multi, layer.indexed, hasName, name, minzoom, score, attrIdx, true); + } else { + MultiPolygon out; + geom::intersection(multi, clippingBox, out); + if (boost::size(out)>0) { + shpMemTiles.StoreGeometry(layerNum, layer.name, POLYGON_, out, layer.indexed, hasName, name, minzoom, score, attrIdx); + } } } else { diff --git a/test/geom.test.cpp b/test/geom.test.cpp index e33b5990..2a31d97d 100644 --- a/test/geom.test.cpp +++ b/test/geom.test.cpp @@ -116,6 +116,44 @@ MU_TEST(test_precheck_many_holes) { checkPrecheck(mp, false); } +static Ring squareRing(double x, double y, double size) { + Ring r; + geom::append(r, Point(x, y)); geom::append(r, Point(x + size, y)); + geom::append(r, Point(x + size, y + size)); geom::append(r, Point(x, y + size)); + geom::append(r, Point(x, y)); + return r; +} + +MU_TEST(test_precheck_large_polygon) { + // A 40000-vertex outer ring (a circle) with 3600 holes: far more than walking the + // outer ring for every hole would allow (the old 2^26 work cap) + Polygon p; + const int n = 40000; + for (int i = 0; i <= n; i++) { + const double a = 2 * M_PI * (i % n) / n; + geom::append(p.outer(), Point(100 + 100 * std::cos(a), 100 + 100 * std::sin(a))); + } + for (int i = 0; i < 60; i++) + for (int j = 0; j < 60; j++) + p.inners().push_back(squareRing(40 + i * 2, 40 + j * 2, 1)); + MultiPolygon mp; + mp.push_back(p); + geom::correct(mp); + checkPrecheck(mp, true); + + // A hole outside the outer ring + MultiPolygon outside = mp; + outside[0].inners().push_back(squareRing(300, 100, 1)); + geom::correct(outside); + checkPrecheck(outside, false); + + // A hole nested in a hole on the other side of the polygon + MultiPolygon nested = mp; + nested[0].inners().insert(nested[0].inners().begin(), squareRing(40.25, 40.25, 0.5)); + geom::correct(nested); + checkPrecheck(nested, false); +} + MU_TEST_SUITE(test_suite_geom) { MU_RUN_TEST(test_rect_clip_concave); MU_RUN_TEST(test_rect_clip_holes); @@ -124,6 +162,7 @@ MU_TEST_SUITE(test_suite_geom) { MU_RUN_TEST(test_rect_clip_touch_from_inside); MU_RUN_TEST(test_precheck_overlapping_bboxes); MU_RUN_TEST(test_precheck_many_holes); + MU_RUN_TEST(test_precheck_large_polygon); } int main() {