perf(precheck): locate holes with a band index, validate shapefile polygons once

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<uint> QueryMatchingGeometries(

src/geom.cpp

@@ -177,13 +177,17 @@ struct PrecheckBox {
 };
 struct PrecheckScratch {
 	std::vector<PrecheckSeg> segs;
-	std::vector<PrecheckBox> boxes, holeBoxes;
+	std::vector<PrecheckBox> boxes;
 	std::vector<std::pair<uint32_t, uint32_t>> pairs;
+	std::vector<uint32_t> bandStart, bandFill, bandSegs, crossed;
 	void trim() {
 		if (segs.capacity() > kPrecheckRetained) std::vector<PrecheckSeg>().swap(segs);
 		if (boxes.capacity() > kPrecheckRetained) std::vector<PrecheckBox>().swap(boxes);
-		if (holeBoxes.capacity() > kPrecheckRetained) std::vector<PrecheckBox>().swap(holeBoxes);
 		if (pairs.capacity() > kPrecheckRetained) std::vector<std::pair<uint32_t, uint32_t>>().swap(pairs);
+		if (bandStart.capacity() > kPrecheckRetained) std::vector<uint32_t>().swap(bandStart);
+		if (bandFill.capacity() > kPrecheckRetained) std::vector<uint32_t>().swap(bandFill);
+		if (bandSegs.capacity() > kPrecheckRetained) std::vector<uint32_t>().swap(bandSegs);
+		if (crossed.capacity() > kPrecheckRetained) std::vector<uint32_t>().swap(crossed);
 	}
 };
 thread_local PrecheckScratch precheckScratch;
@@ -330,20 +334,27 @@ bool segmentsClear(std::vector<PrecheckSeg> &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<PrecheckSeg> &segs = scratch.segs;
+	if (inners.empty()) return true;
+	if (segs.size() >= std::numeric_limits<uint32_t>::max()) return false;
+
+	double minY = std::numeric_limits<double>::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<std::size_t>(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<std::size_t>(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<uint32_t>::max()) return false;
+
+	// Segments of band b: bandSegs[bandStart[b] .. bandStart[b+1])
+	std::vector<uint32_t> &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<uint32_t>(k);
+
+	// 1 = h inside the outer ring and no other hole, -1 = not, 0 = unsure
+	std::vector<uint32_t> &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<std::size_t>(m, 8);
+		int loc = 0;
+		for (std::size_t k = 0; k < tries && loc == 0; k++)
+			loc = locate(hole[k * m / tries], static_cast<uint32_t>(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<uint32_t>::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<PrecheckSeg> &segs = scratch.segs;
 	segs.clear();
-	std::vector<PrecheckBox> &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<uint32_t>(i + 1), segs, holeBoxes[i])) return false;
-		holeBoxes[i].id = static_cast<uint32_t>(i);
-	}
+	PrecheckBox holeBox;
+	for (std::size_t i = 0; i < inners.size(); i++)
+		if (!ringClear(inners[i], true, tol, static_cast<uint32_t>(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:

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<Point> 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<MultiPolygon>(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<OutputObject>(geomType, layerNum, oid, attrIdx, minzoom);
 			std::vector<OutputObject> oolist { *oo };
 			addGeometryToIndex(boost::get<MultiPolygon>(geometry), oolist, 0);

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 {

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() {