Intersects: faster queries for negative case (#635)

systemed · Jan 8, 2024 · 3ff08a3 · 3ff08a3
1 parent fcddc4a
commit 3ff08a3
Show file tree

Hide file tree

Showing 3 changed files with 55 additions and 0 deletions.
diff --git a/include/shp_mem_tiles.h b/include/shp_mem_tiles.h
@@ -35,6 +35,7 @@ class ShpMemTiles : public TileDataSource
 		std::function<std::vector<IndexValue>(const RTree& rtree)> indexQuery, 
 		std::function<bool(const OutputObject& oo)> checkQuery
 	) const;
+	bool mayIntersect(const std::string& layerName, const Box& box) const;
 	std::vector<std::string> namesOfGeometries(const std::vector<uint>& ids) const;
 
 	template <typename GeometryT>
@@ -65,6 +66,7 @@ class ShpMemTiles : public TileDataSource
 	std::map<uint, std::string> indexedGeometryNames;			//  | optional names for each one
 	std::map<std::string, RTree> indices;			// Spatial indices, boost::geometry::index objects for shapefile indices
 	std::mutex indexMutex;
+	std::map<std::string, std::vector<bool>> bitIndices; // A bit it set if the z14 (or base zoom) tiles at x*width + y contains a shape. This lets us quickly reject negative Intersects queryes
 };
 
 #endif //_OSM_MEM_TILES

diff --git a/src/osm_lua_processing.cpp b/src/osm_lua_processing.cpp
@@ -178,6 +178,9 @@ double OsmLuaProcessing::AreaIntersecting(const string &layerName) {
 template <typename GeometryT>
 std::vector<uint> OsmLuaProcessing::intersectsQuery(const string &layerName, bool once, GeometryT &geom) const {
 	Box box; geom::envelope(geom, box);
+	if (!shpMemTiles.mayIntersect(layerName, box))
+		return std::vector<uint>();
+
 	std::vector<uint> ids = shpMemTiles.QueryMatchingGeometries(layerName, once, box,
 		[&](const RTree &rtree) { // indexQuery
 			vector<IndexValue> results;

diff --git a/src/shp_mem_tiles.cpp b/src/shp_mem_tiles.cpp
@@ -55,6 +55,36 @@ vector<string> ShpMemTiles::namesOfGeometries(const vector<uint>& ids) const {
 
 void ShpMemTiles::CreateNamedLayerIndex(const std::string& layerName) {
 	indices[layerName]=RTree();
+
+	bitIndices[layerName] = std::vector<bool>();
+	const uint8_t indexZoom = std::min(baseZoom, 14u);
+	bitIndices[layerName].resize((1 << indexZoom) * (1 << indexZoom));
+}
+
+bool ShpMemTiles::mayIntersect(const std::string& layerName, const Box& box) const {
+	// Check if any tiles in the bitmap might intersect this shape.
+	// If none, downstream code can skip querying the r-tree.
+	auto& bitvec = bitIndices.at(layerName);
+	const uint8_t indexZoom = std::min(baseZoom, 14u);
+
+	double lon1 = box.min_corner().x();
+	double lat1 = box.min_corner().y();
+	double lon2 = box.max_corner().x();
+	double lat2 = box.max_corner().y();
+
+	uint32_t x1 = lon2tilex(lon1, indexZoom);
+	uint32_t x2 = lon2tilex(lon2, indexZoom);
+	uint32_t y1 = lat2tiley(lat1, indexZoom);
+	uint32_t y2 = lat2tiley(lat2, indexZoom);
+
+	for (int x = std::min(x1, x2); x < std::max(x1, x2); x++) {
+		for (int y = std::min(y1, y2); y < std::max(y1, y2); y++) {
+			if (bitvec[x * (1 << indexZoom) + y])
+				return true;
+		}
+	}
+
+	return false;
 }
 
 void ShpMemTiles::StoreGeometry(
@@ -126,4 +156,24 @@ void ShpMemTiles::StoreGeometry(
 	indices.at(layerName).insert(std::make_pair(box, id));
 	if (hasName) { indexedGeometryNames[id] = name; }
 	indexedGeometries.push_back(*oo);
+
+	// Store a bitmap of which tiles at the basezoom might intersect
+	// this shape.
+	auto& bitvec = bitIndices.at(layerName);
+	const uint8_t indexZoom = std::min(baseZoom, 14u);
+	double lon1 = box.min_corner().x();
+	double lat1 = box.min_corner().y();
+	double lon2 = box.max_corner().x();
+	double lat2 = box.max_corner().y();
+
+	uint32_t x1 = lon2tilex(lon1, indexZoom);
+	uint32_t x2 = lon2tilex(lon2, indexZoom);
+	uint32_t y1 = lat2tiley(lat1, indexZoom);
+	uint32_t y2 = lat2tiley(lat2, indexZoom);
+
+	for (int x = std::min(x1, x2); x < std::max(x1, x2); x++) {
+		for (int y = std::min(y1, y2); y < std::max(y1, y2); y++) {
+			bitvec[x * (1 << indexZoom) + y] = true;
+		}
+	}
 }