Clipping monster polygons

include/tile_data.h

@@ -248,7 +248,7 @@ class TileDataSource {
 	uint16_t z6OffsetDivisor;
 
 	// Guards objects, objectsWithIds.
-	std::vector<std::mutex> objectsMutex;
+	mutable std::vector<std::mutex> objectsMutex;
 	// The top-level vector has 1 entry per z6 tile, indexed by x*64 + y
 	// The inner vector contains the output objects that are contained in that z6 tile
 	//
@@ -291,6 +291,33 @@ class TileDataSource {
 	std::mutex multi_linestring_store_mutex;
 	std::unique_ptr<multi_linestring_store_t> multi_linestring_store;
 
+	// Some polygons are very large, e.g. Hudson Bay covers hundreds of thousands of
+	// z14 tiles.
+	//
+	// Because it is very large, it is expensive to clip it once, let alone hundreds
+	// of thousands of times.
+	//
+	// Fortunately, we can avoid redundant clipping. If we clip a polygon and observe
+	// that it fills its entire tile, we know that it must fill all its descendant
+	// tiles at higher zooms.
+	//
+	// We divide the map into 4,096 z6 tiles, then track polygons within that z6
+	// tile if they'd cover at least an entire z9 tile.
+	//
+	// Because z9 is 3 zooms higher than z6, we can use a uint64_t as a bitset
+	// to track the tiles that would be entirely filled.
+	//
+	// To minimize lock contention, we shard by NodeID.
+	std::vector<std::map<std::pair<uint32_t, NodeID>, uint64_t>> jumboCoveringPolygons;
+	// Like above, but for the case where we're in the hole of the polygon,
+	// and so render nothing.
+	std::vector<std::map<std::pair<uint32_t, NodeID>, uint64_t>> jumboExcludedPolygons;
+
+	// Like above, but cover z10..z13, rather than z6..z9
+	std::vector<std::map<std::pair<uint32_t, NodeID>, uint64_t>> largeCoveringPolygons;
+	std::vector<std::map<std::pair<uint32_t, NodeID>, uint64_t>> largeExcludedPolygons;
+
+
 public:
 	TileDataSource(size_t threadNum, unsigned int baseZoom, bool includeID);
 
@@ -335,7 +362,7 @@ class TileDataSource {
 		TileCoordinates coordinates
 	);
 
-	Geometry buildWayGeometry(OutputGeometryType const geomType, NodeID const objectID, const TileBbox &bbox) const;
+	Geometry buildWayGeometry(OutputGeometryType const geomType, NodeID const objectID, const TileBbox &bbox);
 	LatpLon buildNodeGeometry(OutputGeometryType const geomType, NodeID const objectID, const TileBbox &bbox) const;
 
 	void open() {
@@ -427,6 +454,7 @@ class TileDataSource {
 
 
 private:	
+	void updateLargePolygonMaps(NodeID objectID, const TileBbox& bbox, const MultiPolygon& mp);
 };
 
 TileCoordinatesSet getTilesAtZoom(

src/tile_data.cpp

@@ -2,6 +2,7 @@
 #include <iostream>
 #include "tile_data.h"
 
+#include <bitset>
 #include <ciso646>
 
 using namespace std;
@@ -14,7 +15,11 @@ TileDataSource::TileDataSource(size_t threadNum, unsigned int baseZoom, bool inc
 	objectsMutex(threadNum * 4),
 	objects(CLUSTER_ZOOM_AREA),
 	objectsWithIds(CLUSTER_ZOOM_AREA),
-	baseZoom(baseZoom)
+	baseZoom(baseZoom),
+	jumboCoveringPolygons(threadNum * 4),
+	jumboExcludedPolygons(threadNum * 4),
+	largeCoveringPolygons(threadNum * 4),
+	largeExcludedPolygons(threadNum * 4)
 {
 }
 
@@ -127,9 +132,48 @@ void TileDataSource::collectLargeObjectsForTile(
 	}
 }
 
+bool bboxInLargePolygonBitset(uint16_t rootZoom, const TileBbox& bbox, uint64_t bits) {
+	std::bitset<64> bitset(bits);
+
+	const uint16_t leafZoom = rootZoom + 3;
+	const size_t rootX = bbox.index.x / (1 << (bbox.zoom - rootZoom));
+	const size_t rootY = bbox.index.y / (1 << (bbox.zoom - rootZoom));
+
+	// If zoom >= leafZoom: can short-circuit if our tile is set.
+	if (bbox.zoom >= leafZoom) {
+		size_t leafX = bbox.index.x / (1 << (bbox.zoom - leafZoom));
+		size_t leafY = bbox.index.y / (1 << (bbox.zoom - leafZoom));
+		return bitset.test((leafX - (rootX * 8)) * 8 + (leafY - (rootY * 8)));
+	}
+
+	// If zoom is 7 or 8, it's a bit more involved -- we need our entire
+	// range of z9 tiles to be set.
+	if (bbox.zoom == rootZoom + 1 || bbox.zoom == rootZoom + 2) {
+		const size_t stride = 1 << (leafZoom - bbox.zoom);
+
+		const size_t leafX = rootX * 8;
+		const size_t leafY = rootY * 8;
+
+		// Figure out where we'll start reading.
+		const size_t startX = bbox.index.x * (1 << (leafZoom - bbox.zoom)) - leafX;
+		const size_t startY = bbox.index.y * (1 << (leafZoom - bbox.zoom)) - leafY;
+
+		bool success = true;
+		for (int x = startX; x < startX + stride; x++) {
+			for (int y = startY; y < startY + stride; y++) {
+				success = success && bitset.test(x * 8 + y);
+			}
+		}
+		return success;
+	}
+
+	// You shouldn't call this for zooms <= rootZoom.
+	return false;
+}
+
 // Build node and way geometries
 Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType, 
-                                          NodeID const objectID, const TileBbox &bbox) const {
+                                          NodeID const objectID, const TileBbox &bbox) {
 	switch(geomType) {
 		case POINT_: {
 			auto p = retrieve_point(objectID);
@@ -175,6 +219,76 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType,
 		}
 
 		case POLYGON_: {
+			if (bbox.zoom > 6) {
+				// If we're at a high enough zoom, see if we can short-circuit clipping.
+				//
+				// This is possible for very large polygons that cover (or fail to cover)
+				// an entire z6/z7/z8/z9 tile.
+
+				const size_t rootXz6 = bbox.index.x / (1 << (bbox.zoom - 6));
+				const size_t rootYz6 = bbox.index.y / (1 << (bbox.zoom - 6));
+				const uint16_t z6Tile = rootXz6 * (1 << 6) + rootYz6;
+
+				const size_t rootXz10 = bbox.index.x / (1 << (bbox.zoom - 10));
+				const size_t rootYz10 = bbox.index.y / (1 << (bbox.zoom - 10));
+				const uint16_t z10Tile = rootXz10 * (1 << 10) + rootYz10;
+
+				const uint16_t coveringShard = objectID % jumboCoveringPolygons.size();
+				const uint16_t excludedShard = objectID % jumboExcludedPolygons.size();
+				const uint16_t lockShard = objectID % objectsMutex.size();
+				uint64_t jumboCoveringBits = 0;
+				uint64_t jumboExcludedBits = 0;
+				uint64_t largeCoveringBits = 0;
+				uint64_t largeExcludedBits = 0;
+				{
+					std::lock_guard<std::mutex> lock(objectsMutex[lockShard]);
+					{
+						const auto& rv = jumboCoveringPolygons[coveringShard].find(std::make_pair(z6Tile, objectID));
+						if (rv != jumboCoveringPolygons[coveringShard].end()) {
+							jumboCoveringBits = rv->second;
+						}
+					}
+
+					{
+						const auto& rv = jumboExcludedPolygons[excludedShard].find(std::make_pair(z6Tile, objectID));
+						if (rv != jumboExcludedPolygons[excludedShard].end()) {
+							jumboExcludedBits = rv->second;
+						}
+					}
+
+					if (bbox.zoom > 10) {
+						{
+							const auto& rv = largeCoveringPolygons[coveringShard].find(std::make_pair(z10Tile, objectID));
+							if (rv != largeCoveringPolygons[coveringShard].end()) {
+								largeCoveringBits = rv->second;
+							}
+						}
+
+						{
+							const auto& rv = largeExcludedPolygons[excludedShard].find(std::make_pair(z10Tile, objectID));
+							if (rv != largeExcludedPolygons[excludedShard].end()) {
+								largeExcludedBits = rv->second;
+							}
+						}
+					}
+				}
+
+				if ((jumboCoveringBits != 0 && bbox.zoom > 6 && bboxInLargePolygonBitset(6, bbox, jumboCoveringBits)) ||
+						(largeCoveringBits != 0 && bbox.zoom > 10 && bboxInLargePolygonBitset(10, bbox, largeCoveringBits))) {
+					MultiPolygon mp;
+					Polygon p;
+					boost::geometry::assign(p, bbox.clippingBox);
+					mp.push_back(p);
+					return mp;
+				}
+
+				if ((jumboExcludedBits != 0 && bbox.zoom > 6 && bboxInLargePolygonBitset(6, bbox, jumboExcludedBits)) ||
+						(largeExcludedBits != 0 && bbox.zoom > 10 && bboxInLargePolygonBitset(10, bbox, largeExcludedBits))) {
+					MultiPolygon mp;
+					return mp;
+				}
+			}
+
 			auto const &input = retrieve_multi_polygon(objectID);
 
 			Box box = bbox.clippingBox;
@@ -235,11 +349,15 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType,
 					MultiPolygon output;
 					geom::intersection(input, box, output);
 					geom::correct(output);
+
+					updateLargePolygonMaps(objectID, bbox, output);
 					return output;
 				} else {
 					// occasionally also wrong_topological_dimension, disconnected_interior
 				}
 			}
+
+			updateLargePolygonMaps(objectID, bbox, mp);
 			return mp;
 		}
 
@@ -248,6 +366,73 @@ Geometry TileDataSource::buildWayGeometry(OutputGeometryType const geomType,
 	}
 }
 
+void TileDataSource::updateLargePolygonMaps(NodeID objectID, const TileBbox& bbox, const MultiPolygon& mp) {
+	// Did this polygon clip to the full extent of the tile, or clip to nothing?
+	//
+	// If yes, flag it so future tiles at higher zooms can avoid doing an
+	// expensive clip.
+	if (bbox.zoom < 6 || bbox.zoom > 13)
+		return;
+
+	bool isExcluded = false;
+	bool isCovered = false;
+
+	if (mp.size() == 0) {
+		isExcluded = true;
+	} else if (mp.size() == 1 && mp[0].outer().size() == 5) {
+		Polygon bboxAsPolygon;
+		boost::geometry::assign(bboxAsPolygon, bbox.clippingBox);
+
+		isCovered = boost::geometry::equals(mp[0], bboxAsPolygon);
+	}
+
+	if (!isExcluded && !isCovered)
+		return;
+
+	const bool isJumbo = bbox.zoom <= 9;
+	const uint16_t rootZoom = isJumbo ? 6 : 10;
+	const uint16_t leafZoom = rootZoom + 3;
+
+	// Truncate to the root x/y tile
+	const size_t rootX = bbox.index.x / (1 << (bbox.zoom - rootZoom));
+	const size_t rootY = bbox.index.y / (1 << (bbox.zoom - rootZoom));
+
+	const uint32_t rootTile = rootX * (1 << rootZoom) + rootY;
+
+	auto& store = isJumbo ? (isCovered ? jumboCoveringPolygons : jumboExcludedPolygons)
+		: (isCovered ? largeCoveringPolygons : largeExcludedPolygons);
+	const uint16_t shard = objectID % store.size();
+	const uint16_t lockShard = objectID % objectsMutex.size();
+
+	uint64_t bits = 0;
+	const std::pair<uint32_t, NodeID> key = std::make_pair(rootTile, objectID);
+	std::lock_guard<std::mutex> lock(objectsMutex[lockShard]);
+	const auto& rv = store[shard].find(key);
+	if (rv != store[shard].end())
+		bits = rv->second;
+
+	std::bitset<64> covered(bits);
+
+	// We're at zoom 6, 7, 8, or 9, so we'll cover an 8x8, 4x4, 2x2 or 1x1 area
+	// of z9 tiles, and similarly for z10, .., z13.
+	const size_t stride = 1 << (leafZoom - bbox.zoom);
+
+	const size_t leafX = rootX * 8;
+	const size_t leafY = rootY * 8;
+
+	// Figure out where we'll start writing.
+	const size_t startX = bbox.index.x * (1 << (leafZoom - bbox.zoom)) - leafX;
+	const size_t startY = bbox.index.y * (1 << (leafZoom - bbox.zoom)) - leafY;
+
+	for (int x = startX; x < startX + stride; x++) {
+		for (int y = startY; y < startY + stride; y++) {
+			covered.set(x * 8 + y);
+		}
+	}
+
+	store[shard][key] = covered.to_ullong();
+}
+
 LatpLon TileDataSource::buildNodeGeometry(OutputGeometryType const geomType, 
                                           NodeID const objectID, const TileBbox &bbox) const {
 	switch(geomType) {
@@ -406,7 +591,7 @@ void TileDataSource::addGeometryToIndex(
 			tileSet.insert(tileSetTmp.begin(), tileSetTmp.end());
 		}
 	}
-	
+
 	TileCoordinate minTileX = TILE_COORDINATE_MAX, maxTileX = 0, minTileY = TILE_COORDINATE_MAX, maxTileY = 0;
 	for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {
 		TileCoordinates index = *it;

src/tile_worker.cpp

@@ -96,7 +96,7 @@ void MergeIntersecting(MultiPolygon &input, MultiPolygon &to_merge) {
 
 template <typename T>
 void CheckNextObjectAndMerge(
-	const TileDataSource* source,
+	TileDataSource* source,
 	OutputObjectsConstIt& jt,
 	OutputObjectsConstIt ooSameLayerEnd, 
 	const TileBbox& bbox,
@@ -147,7 +147,7 @@ void RemoveInnersBelowSize(MultiPolygon &g, double filterArea) {
 }
 
 void ProcessObjects(
-	const TileDataSource* source,
+	TileDataSource* source,
 	const AttributeStore& attributeStore,
 	OutputObjectsConstIt ooSameLayerBegin,
 	OutputObjectsConstIt ooSameLayerEnd, 

src/tilemaker.cpp

@@ -434,6 +434,7 @@ int main(int argc, char* argv[]) {
 
 		// Launch the pool with threadNum threads
 		boost::asio::thread_pool pool(threadNum);
+		//boost::asio::thread_pool pool(2);
 
 		// Mutex is hold when IO is performed
 		std::mutex io_mutex;
@@ -543,6 +544,37 @@ int main(int argc, char* argv[]) {
 					unsigned int zoom = tileCoordinates[i].first;
 					TileCoordinates coords = tileCoordinates[i].second;
 					std::vector<std::vector<OutputObjectID>> data;
+
+					if (zoom >= 6) {
+						const size_t rootX = coords.x / (1 << (zoom - 6));
+						const size_t rootY = coords.y / (1 << (zoom - 6));
+
+						// drill into 6/17/20, it's a good test case.
+						//if (!(rootX == 17 && rootY == 20)) continue;
+					}
+					/*
+					if (!(
+							(zoom == 9 && 
+								((coords.x == 143 && (coords.y == 159 || coords.y == 160)) ||
+								(coords.x == 144 && (coords.y == 159 || coords.y == 160)))) ||
+							(zoom == 8 && 
+							 ((coords.x == 71 && coords.y == 79) ||
+								(coords.x == 71 && coords.y == 80) ||
+								(coords.x == 72 && coords.y == 79) ||
+								(coords.x == 72 && coords.y == 80)
+														)) ||
+							(zoom == 7 &&
+							 ((coords.x == 35 && coords.y == 39) ||
+								(coords.x == 35 && coords.y == 40) ||
+								(coords.x == 36 && coords.y == 39) ||
+								(coords.x == 36 && coords.y == 40))) ||
+							(zoom == 6 &&
+							 ((coords.x == 17 && coords.y == 19) ||
+								(coords.x == 17 && coords.y == 20) ||
+								(coords.x == 18 && coords.y == 19) ||
+								(coords.x == 18 && coords.y == 20)))
+							)) continue;
+					*/
 					for (auto source : sources) {
 						data.emplace_back(source->getObjectsForTile(sortOrders, zoom, coords));
 					}