Use rtree for large polygons (systemed#479)

include/coordinates.h

@@ -9,8 +9,10 @@
 
 #ifdef FAT_TILE_INDEX
 typedef uint32_t TileCoordinate;
+#define TILE_COORDINATE_MAX UINT32_MAX
 #else
 typedef uint16_t TileCoordinate;
+#define TILE_COORDINATE_MAX UINT16_MAX
 #endif
 class TileCoordinates_ {
 

include/shp_mem_tiles.h

@@ -14,7 +14,7 @@ class ShpMemTiles : public TileDataSource
 	void CreateNamedLayerIndex(const std::string &layerName);
 
 	// Used in shape file loading
-	OutputObjectRef AddObject(uint_least8_t layerNum,
+	OutputObjectRef StoreShapefileGeometry(uint_least8_t layerNum,
 		const std::string &layerName, 
 		enum OutputGeometryType geomType,
 		Geometry geometry, 
@@ -23,6 +23,10 @@ class ShpMemTiles : public TileDataSource
 	void AddObject(TileCoordinates const &index, OutputObjectRef const &oo) {
 		tileIndex[index].push_back(oo);
 	}
+	void AddObjectToLargeIndex(Box const &envelope, OutputObjectRef const &oo) {
+		std::lock_guard<std::mutex> lock(mutex);
+		box_rtree.insert(std::make_pair(envelope, oo));
+	}
 	std::vector<uint> QueryMatchingGeometries(const std::string &layerName, bool once, Box &box, 
 		std::function<std::vector<IndexValue>(const RTree &rtree)> indexQuery, 
 		std::function<bool(OutputObject const &oo)> checkQuery) const;

include/tile_data.h

@@ -19,6 +19,10 @@ class TileDataSource {
 	std::mutex mutex;
 	TileIndex tileIndex;
 	std::deque<OutputObject> objects;
+
+	// rtree index of large objects
+	using oo_rtree_param_type = boost::geometry::index::quadratic<128>;
+	boost::geometry::index::rtree< std::pair<Box,OutputObjectRef>, oo_rtree_param_type> box_rtree;
 
 	unsigned int baseZoom;
 
@@ -32,6 +36,8 @@ class TileDataSource {
 		MergeTileCoordsAtZoom(zoom, baseZoom, tileIndex, dstCoords);
 	}
 
+	void MergeLargeCoordsAtZoom(uint zoom, TileCoordinatesSet &dstCoords);
+
 	///This must be thread safe!
 	void MergeSingleTileDataAtZoom(TileCoordinates dstIndex, uint zoom, std::vector<OutputObjectRef> &dstTile) {
 		MergeSingleTileDataAtZoom(dstIndex, zoom, baseZoom, tileIndex, dstTile);
@@ -48,6 +54,13 @@ class TileDataSource {
 		tileIndex[index].push_back(oo);
 	}
 
+	void AddObjectToLargeIndex(Box const &envelope, OutputObjectRef const &oo) {
+		std::lock_guard<std::mutex> lock(mutex);
+		box_rtree.insert(std::make_pair(envelope, oo));
+	}
+
+	void MergeLargeObjects(TileCoordinates dstIndex, uint zoom, std::vector<OutputObjectRef> &dstTile);
+
 private:	
 	static void MergeTileCoordsAtZoom(uint zoom, uint baseZoom, const TileIndex &srcTiles, TileCoordinatesSet &dstCoords);
 	static void MergeSingleTileDataAtZoom(TileCoordinates dstIndex, uint zoom, uint baseZoom, const TileIndex &srcTiles, std::vector<OutputObjectRef> &dstTile);

src/osm_lua_processing.cpp

@@ -609,29 +609,45 @@ void OsmLuaProcessing::setWay(WayID wayId, LatpLonVec const &llVec, const tag_ma
 		// create a list of tiles this way passes through (tileSet)
 		unordered_set<TileCoordinates> tileSet;
 		try {
-			insertIntermediateTiles(osmStore.llListLinestring(llVecPtr->cbegin(),llVecPtr->cend()), this->config.baseZoom, tileSet);
+			Linestring ls = osmStore.llListLinestring(llVecPtr->cbegin(),llVecPtr->cend());
+			insertIntermediateTiles(ls, this->config.baseZoom, tileSet);
 
 			// then, for each tile, store the OutputObject for each layer
 			bool polygonExists = false;
+			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;
+				minTileX = std::min(index.x, minTileX);
+				minTileY = std::min(index.y, minTileY);
+				maxTileX = std::max(index.x, maxTileX);
+				maxTileY = std::max(index.y, maxTileY);
 				for (auto jt = this->outputs.begin(); jt != this->outputs.end(); ++jt) {
 					if (jt->first->geomType == POLYGON_) {
 						polygonExists = true;
 						continue;
 					}
-					osmMemTiles.AddObject(index, jt->first);
+					osmMemTiles.AddObject(index, jt->first); // not a polygon
 				}
 			}
 
 			// for polygon, fill inner tiles
 			if (polygonExists) {
-				fillCoveredTiles(tileSet);
-				for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {
-					TileCoordinates index = *it;
-					for (auto jt = this->outputs.begin(); jt != this->outputs.end(); ++jt) {
-						if (jt->first->geomType != POLYGON_) continue;
-						osmMemTiles.AddObject(index, jt->first);
+				bool tilesetFilled = false;
+				uint size = (maxTileX - minTileX + 1) * (maxTileY - minTileY + 1);
+				for (auto jt = this->outputs.begin(); jt != this->outputs.end(); ++jt) {
+					if (jt->first->geomType != POLYGON_) continue;
+					if (size>= 16) {
+						// Larger objects - add to rtree
+						Box box = Box(geom::make<Point>(minTileX, minTileY),
+						              geom::make<Point>(maxTileX, maxTileY));
+						osmMemTiles.AddObjectToLargeIndex(box, jt->first);
+					} else {
+						// Smaller objects - add to each individual tile index
+						if (!tilesetFilled) { fillCoveredTiles(tileSet); tilesetFilled = true; }
+						for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {
+							TileCoordinates index = *it;
+							osmMemTiles.AddObject(index, jt->first);
+						}
 					}
 				}
 			}
@@ -676,28 +692,46 @@ void OsmLuaProcessing::setRelation(int64_t relationId, WayVec const &outerWayVec
 			return;
 		}		
 
-		unordered_set<TileCoordinates> tileSet;
-		if (mp.size() == 1) {
-			insertIntermediateTiles(mp[0].outer(), this->config.baseZoom, tileSet);
-			fillCoveredTiles(tileSet);
-		} else {
-			for (Polygon poly: mp) {
-				unordered_set<TileCoordinates> tileSetTmp;
-				insertIntermediateTiles(poly.outer(), this->config.baseZoom, tileSetTmp);
-				fillCoveredTiles(tileSetTmp);
-				tileSet.insert(tileSetTmp.begin(), tileSetTmp.end());
-			}
-		}
-
 		for (auto jt = this->outputs.begin(); jt != this->outputs.end(); ++jt) {
 			// Store the attributes of the generated geometry
 			jt->first->setAttributeSet(attributeStore.store_set(jt->second));		
 		}
 
+		unordered_set<TileCoordinates> tileSet;
+		bool singleOuter = mp.size()==1;
+		for (Polygon poly: mp) {
+			unordered_set<TileCoordinates> tileSetTmp;
+			insertIntermediateTiles(poly.outer(), this->config.baseZoom, tileSetTmp);
+			fillCoveredTiles(tileSetTmp);
+			if (singleOuter) {
+				tileSet = std::move(tileSetTmp);
+			} else {
+				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;
-			for (auto jt = this->outputs.begin(); jt != this->outputs.end(); ++jt) {
-				osmMemTiles.AddObject(index, jt->first);
+			minTileX = std::min(index.x, minTileX);
+			minTileY = std::min(index.y, minTileY);
+			maxTileX = std::max(index.x, maxTileX);
+			maxTileY = std::max(index.y, maxTileY);
+		}
+		for (auto jt = this->outputs.begin(); jt != this->outputs.end(); ++jt) {
+			if (tileSet.size()>=16) {
+				// Larger objects - add to rtree
+				// note that the bbox is currently the envelope of the entire multipolygon,
+				// which is suboptimal in shapes like (_) ...... (_) where the outers are significantly disjoint
+				Box box = Box(geom::make<Point>(minTileX, minTileY),
+				              geom::make<Point>(maxTileX, maxTileY));
+				osmMemTiles.AddObjectToLargeIndex(box, jt->first);
+			} else {
+				// Smaller objects - add to each individual tile index
+				for (auto it = tileSet.begin(); it != tileSet.end(); ++it) {
+					TileCoordinates index = *it;
+					osmMemTiles.AddObject(index, jt->first);
+				}
 			}
 		}
 

src/read_shp.cpp

@@ -177,7 +177,7 @@ void readShapefile(const Box &clippingBox,
 				if (indexField>-1) { name=DBFReadStringAttribute(dbf, i, indexField); hasName = true;}
 
 				auto &attributeStore = osmLuaProcessing.getAttributeStore();
-				OutputObjectRef oo = shpMemTiles.AddObject(layerNum, layerName, POINT_, p, isIndexed, hasName, name, attributeStore.empty_set(), layer.minzoom);
+				OutputObjectRef oo = shpMemTiles.StoreShapefileGeometry(layerNum, layerName, POINT_, p, isIndexed, hasName, name, attributeStore.empty_set(), layer.minzoom);
 
 				addShapefileAttributes(dbf, oo, i, columnMap, columnTypeMap, layers, osmLuaProcessing);
 			}
@@ -199,7 +199,7 @@ void readShapefile(const Box &clippingBox,
 					if (indexField>-1) { name=DBFReadStringAttribute(dbf, i, indexField); hasName = true;}
 
 					auto &attributeStore = osmLuaProcessing.getAttributeStore();
-					OutputObjectRef oo = shpMemTiles.AddObject(layerNum, layerName, LINESTRING_, *it, isIndexed, hasName, name, attributeStore.empty_set(), layer.minzoom);
+					OutputObjectRef oo = shpMemTiles.StoreShapefileGeometry(layerNum, layerName, LINESTRING_, *it, isIndexed, hasName, name, attributeStore.empty_set(), layer.minzoom);
 
 					addShapefileAttributes(dbf, oo, i, columnMap, columnTypeMap, layers, osmLuaProcessing);
 				}
@@ -272,7 +272,7 @@ void readShapefile(const Box &clippingBox,
 
 				// create OutputObject
 				auto &attributeStore = osmLuaProcessing.getAttributeStore();
-				OutputObjectRef oo = shpMemTiles.AddObject(layerNum, layerName, POLYGON_, out, isIndexed, hasName, name, attributeStore.empty_set(), layer.minzoom);
+				OutputObjectRef oo = shpMemTiles.StoreShapefileGeometry(layerNum, layerName, POLYGON_, out, isIndexed, hasName, name, attributeStore.empty_set(), layer.minzoom);
 
 				addShapefileAttributes(dbf, oo, i, columnMap, columnTypeMap, layers, osmLuaProcessing);
 			}

src/shp_mem_tiles.cpp

@@ -51,7 +51,7 @@ void ShpMemTiles::CreateNamedLayerIndex(const std::string &layerName) {
 	indices[layerName]=RTree();
 }
 
-OutputObjectRef ShpMemTiles::AddObject(uint_least8_t layerNum,
+OutputObjectRef ShpMemTiles::StoreShapefileGeometry(uint_least8_t layerNum,
 	const std::string &layerName, enum OutputGeometryType geomType,
 	Geometry geometry, bool isIndexed, bool hasName, const std::string &name, AttributeStoreRef attributes, uint minzoom) {
 
@@ -116,15 +116,23 @@ OutputObjectRef ShpMemTiles::AddObject(uint_least8_t layerNum,
 }
 
 // Add an OutputObject to all tiles between min/max lat/lon
+// (only used for polygons)
 void ShpMemTiles::addToTileIndexByBbox(OutputObjectRef &oo, double minLon, double minLatp, double maxLon, double maxLatp) {
 	uint minTileX =  lon2tilex(minLon, baseZoom);
 	uint maxTileX =  lon2tilex(maxLon, baseZoom);
 	uint minTileY = latp2tiley(minLatp, baseZoom);
 	uint maxTileY = latp2tiley(maxLatp, baseZoom);
-	for (uint x=min(minTileX,maxTileX); x<=max(minTileX,maxTileX); x++) {
-		for (uint y=min(minTileY,maxTileY); y<=max(minTileY,maxTileY); y++) {
-			TileCoordinates index(x, y);
-			AddObject(index, oo);
+	uint size = (maxTileX - minTileX + 1) * (minTileY - maxTileY + 1);
+	if (size>=16) { 
+		// Larger objects - add to rtree
+		AddObjectToLargeIndex(Box(Point(minTileX, maxTileY), Point(maxTileX, minTileY)), oo);
+	} else {
+		// Smaller objects - add to each individual tile index
+		for (uint x=min(minTileX,maxTileX); x<=max(minTileX,maxTileX); x++) {
+			for (uint y=min(minTileY,maxTileY); y<=max(minTileY,maxTileY); y++) {
+				TileCoordinates index(x, y);
+				AddObject(index, oo);
+			}
 		}
 	}
 }

src/tile_data.cpp

@@ -29,6 +29,24 @@ void TileDataSource::MergeTileCoordsAtZoom(uint zoom, uint baseZoom, const TileI
 	}
 }
 
+// Find the tiles used by the "large objects" from the rtree index
+void TileDataSource::MergeLargeCoordsAtZoom(uint zoom, TileCoordinatesSet &dstCoords) {
+	for(auto const &result: box_rtree) {
+		int scale = pow(2, baseZoom-zoom);
+		TileCoordinate minx = result.first.min_corner().x() / scale;
+		TileCoordinate maxx = result.first.max_corner().x() / scale;
+		TileCoordinate miny = result.first.min_corner().y() / scale;
+		TileCoordinate maxy = result.first.max_corner().y() / scale;
+		for (int x=minx; x<=maxx; x++) {
+			for (int y=miny; y<=maxy; y++) {
+				TileCoordinates newIndex(x, y);
+				dstCoords.insert(newIndex);
+			}
+		}
+	}
+}
+
+// Copy objects from the tile at dstIndex (in the dataset srcTiles) into dstTile
 void TileDataSource::MergeSingleTileDataAtZoom(TileCoordinates dstIndex, uint zoom, uint baseZoom, const TileIndex &srcTiles, std::vector<OutputObjectRef> &dstTile) {
 	if (zoom==baseZoom) {
 		// at z14, we can just use tileIndex
@@ -59,22 +77,37 @@ void TileDataSource::MergeSingleTileDataAtZoom(TileCoordinates dstIndex, uint zo
 	}
 }
 
+// Copy objects from the large index into dstTile
+void TileDataSource::MergeLargeObjects(TileCoordinates dstIndex, uint zoom, std::vector<OutputObjectRef> &dstTile) {
+	int scale = pow(2, baseZoom - zoom);
+	TileCoordinates srcIndex1( dstIndex.x   *scale  ,  dstIndex.y   *scale  );
+	TileCoordinates srcIndex2((dstIndex.x+1)*scale-1, (dstIndex.y+1)*scale-1);
+	Box box = Box(geom::make<Point>(srcIndex1.x, srcIndex1.y),
+	              geom::make<Point>(srcIndex2.x, srcIndex2.y));
+	for(auto const &result: box_rtree | boost::geometry::index::adaptors::queried(boost::geometry::index::intersects(box)))
+		dstTile.push_back(result.second);
+}
+
 TileCoordinatesSet GetTileCoordinates(std::vector<class TileDataSource *> const &sources, unsigned int zoom) {
 	TileCoordinatesSet tileCoordinates;
 
 	// Create list of tiles
 	tileCoordinates.clear();
-	for(size_t i=0; i<sources.size(); i++)
+	for(size_t i=0; i<sources.size(); i++) {
 		sources[i]->MergeTileCoordsAtZoom(zoom, tileCoordinates);
+		sources[i]->MergeLargeCoordsAtZoom(zoom, tileCoordinates);
+	}
 
 	return tileCoordinates;
 }
 
 std::vector<OutputObjectRef> GetTileData(std::vector<class TileDataSource *> const &sources, TileCoordinates coordinates, unsigned int zoom)
 {
 	std::vector<OutputObjectRef> data;
-	for(size_t i=0; i<sources.size(); i++)
+	for(size_t i=0; i<sources.size(); i++) {
 		sources[i]->MergeSingleTileDataAtZoom(coordinates, zoom, data);
+		sources[i]->MergeLargeObjects(coordinates, zoom, data);
+	}
 
 	boost::sort::pdqsort(data.begin(), data.end());
 	data.erase(unique(data.begin(), data.end()), data.end());