Simplify buildings algorithm (#918)

Author: systemed

CMakeLists.txt

@@ -132,6 +132,7 @@ file(GLOB tilemaker_src_files
 	src/shp_mem_tiles.cpp
 	src/shp_processor.cpp
 	src/significant_tags.cpp
+	src/simplify_buildings.cpp
 	src/sorted_node_store.cpp
 	src/sorted_way_store.cpp
 	src/tag_map.cpp

Makefile

@@ -125,6 +125,7 @@ tilemaker: \
 	src/shp_mem_tiles.o \
 	src/shp_processor.o \
 	src/significant_tags.o \
+	src/simplify_buildings.o \
 	src/sorted_node_store.o \
 	src/sorted_way_store.o \
 	src/tag_map.o \

docs/CONFIGURATION.md

@@ -77,7 +77,7 @@ You can add optional parameters to layers:
 * `simplify_level` - how much to simplify features (in degrees of longitude) on the zoom level `simplify_below-1`
 * `simplify_length` - how much to simplify features (in kilometers) on the zoom level `simplify_below-1`, preceding `simplify_level`
 * `simplify_ratio` - (optional: the default value is 2.0) the actual simplify level will be `simplify_level * pow(simplify_ratio, (simplify_below-1) - <current zoom>)`
-* `simplify_algorithm` - which simplification algorithm to use (defaults to Douglas-Peucker; you can also specify `"visvalingam"`, which can be better for landuse and similar polygons)
+* `simplify_algorithm` - which simplification algorithm to use (defaults to Douglas-Peucker; you can also specify `"visvalingam"`, which can be better for landuse and similar polygons, or `"buildings"`, which preserves rectilinear shapes)
 * `filter_below` - filter areas by minimum size below this zoom level
 * `filter_area` - minimum size (in square degrees of longitude) for the zoom level `filter_below-1`
 * `feature_limit` - restrict the number of features written to each tile

include/geom.h

@@ -34,6 +34,7 @@ typedef boost::geometry::model::polygon<Point> Polygon;
 typedef boost::geometry::model::multi_polygon<Polygon> MultiPolygon;
 typedef boost::geometry::model::multi_linestring<Linestring> MultiLinestring;
 typedef boost::geometry::model::box<Point> Box;
+typedef boost::geometry::model::segment<Point> Segment;
 typedef boost::geometry::ring_type<Polygon>::type Ring;
 typedef boost::geometry::interior_type<Polygon>::type InteriorRing;
 typedef boost::variant<Point,Linestring,MultiLinestring,MultiPolygon> Geometry;

include/shared_data.h

@@ -46,6 +46,7 @@ struct LayerDef {
 
 	static const uint DOUGLAS_PEUCKER = 0;
 	static const uint VISVALINGAM = 1;
+	static const uint BUILDINGS = 2;
 };
 
 ///\brief Defines layers used in map rendering

include/simplify_buildings.h

@@ -0,0 +1,38 @@
+/*! \file */ 
+#ifndef _SIMPLIFY_BUILDINGS_H
+#define _SIMPLIFY_BUILDINGS_H
+
+struct BuildingSimplifyConfig {
+	// "Short" segment: the short side of the stub being removed (in coordinate units).
+	// For lat/lon coordinates use e.g. 5.0/100000.0; scale accordingly for projected coords.
+	double distance_filter = 5.0 / 100000.0;
+
+	// Area guard: don't remove stubs where d1*d2 (product of the two stub sides) exceeds this.
+	// Prevents collapsing significant right-angle jogs.  Set very large to disable.
+	double area_filter = 50.0;
+
+	// Narrow-edge exception: bypass the area guard when either stub side is shorter than this.
+	double area_narrow = 2.0;
+
+	// Minimum sine of the angle between the two outer lines (ring[km1]→ring[k] and
+	// ring[kp2]→ring[kp3]).  When these lines are nearly parallel their intersection is at
+	// infinity; the fallback midpoint passes the length check via the triangle inequality
+	// and destroys a significant corner.  0.1 ≈ sin(5.7°) — reject anything more parallel.
+	static constexpr double PARALLEL_TOL = 0.1;
+
+	// Degrees tolerance for removing collinear points (straight-line removal).
+	static constexpr double COLLINEAR_TOL = 8.0;
+
+	// Degrees tolerance for snapping the intersection to a right angle.
+	static constexpr double SNAP_TOL = 3.0;
+
+	// Parametric tolerance for the proper-intersection test: values in [tol, 1-tol] count.
+	// Keeps endpoint-touching segments from being treated as intersecting.
+	static constexpr double INTERSECT_TOL = 0.001;
+};
+
+void simplify_building(Polygon &poly, const BuildingSimplifyConfig &cfg);
+void simplify_building(MultiPolygon &mp, const BuildingSimplifyConfig &cfg);
+void simplifyBuildings(MultiPolygon &mp, double max_distance);
+
+#endif //_SIMPLIFY_BUILDINGS_H

src/geom.cpp

@@ -225,26 +225,68 @@ bool repair_multi_polygon(MultiPolygon &mp)
 
 // ---------------
 // Union multipolygons
-// from https://github.com/boostorg/geometry/discussions/947
+// Groups polygons into connected components by bbox intersection (via R-tree +
+// union-find), then runs the binary reduction only within each component.
+// Disjoint polygons are concatenated directly, skipping the expensive union_()
+// call that Boost still runs even for non-overlapping geometry.
 void union_many(std::vector<MultiPolygon> &to_unify) {
-	if (to_unify.size()<2) return;
-	size_t step = 1;
-	size_t half_step;
+	if (to_unify.size() < 2) return;
+
+	namespace bgi = boost::geometry::index;
+	typedef std::pair<Box, size_t> BoxIdx;
+
+	std::vector<Box> boxes(to_unify.size());
+	for (size_t i = 0; i < to_unify.size(); i++)
+		boost::geometry::envelope(to_unify[i], boxes[i]);
+
+	// Union-find with path compression
+	std::vector<size_t> parent(to_unify.size());
+	for (size_t i = 0; i < to_unify.size(); i++) parent[i] = i;
+	std::function<size_t(size_t)> find = [&](size_t x) -> size_t {
+		return parent[x] == x ? x : (parent[x] = find(parent[x]));
+	};
+
+	// Incrementally build an R-tree; unite each polygon with all prior ones
+	// whose bboxes intersect (transitivity is handled by union-find)
+	bgi::rtree<BoxIdx, bgi::quadratic<16>> rtree;
+	for (size_t i = 0; i < to_unify.size(); i++) {
+		for (auto const &v : rtree | bgi::adaptors::queried(bgi::intersects(boxes[i]))) {
+			size_t ri = find(i), rj = find(v.second);
+			if (ri != rj) parent[ri] = rj;
+		}
+		rtree.insert({boxes[i], i});
+	}
 
-	// the outer loop doubles the distance between two polygons to be merged at every iteration
-	do {
-		half_step = step;
-		step *= 2;
-		size_t i = 0;
+	// Bucket indices by component root (roots are in 0..n-1)
+	std::vector<std::vector<size_t>> components(to_unify.size());
+	for (size_t i = 0; i < to_unify.size(); i++)
+		components[find(i)].push_back(i);
 
-		// the inner loop merges polygons at i and i+half_step storing the result at i
+	// Union within each component; concatenate all results into to_unify[0]
+	MultiPolygon result;
+	for (auto &comp : components) {
+		if (comp.empty()) continue;
+		if (comp.size() == 1) {
+			for (auto &p : to_unify[comp[0]]) result.push_back(std::move(p));
+			continue;
+		}
+		std::vector<MultiPolygon> sub;
+		sub.reserve(comp.size());
+		for (size_t idx : comp) sub.push_back(std::move(to_unify[idx]));
+
+		size_t step = 1, half_step;
 		do {
-			MultiPolygon unified;
-			boost::geometry::union_(to_unify.at(i), to_unify.at(i + half_step), unified);
-			to_unify.at(i) = std::move(unified);
-			i += step;
-		} while (i + half_step < to_unify.size());
-	} while (step < to_unify.size());
+			half_step = step; step *= 2;
+			for (size_t i = 0; i + half_step < sub.size(); i += step) {
+				MultiPolygon unified;
+				boost::geometry::union_(sub[i], sub[i + half_step], unified);
+				sub[i] = std::move(unified);
+			}
+		} while (step < sub.size());
+
+		for (auto &p : sub[0]) result.push_back(std::move(p));
+	}
+	to_unify[0] = std::move(result);
 }
 
 // ---------------

src/shared_data.cpp

@@ -324,7 +324,8 @@ void Config::readConfig(rapidjson::Document &jsonConfig, bool &hasClippingBox, B
 		int    combinePolyBelow=it->value.HasMember("combine_polygons_below") ? it->value["combine_polygons_below"].GetInt() : 0;
 		bool sortZOrderAscending = it->value.HasMember("z_order_ascending") ? it->value["z_order_ascending"].GetBool() : (featureLimit==0);
 		string algo           = it->value.HasMember("simplify_algorithm") ? it->value["simplify_algorithm"].GetString() : "";
-		uint simplifyAlgo = algo=="visvalingam" ? LayerDef::VISVALINGAM : LayerDef::DOUGLAS_PEUCKER;
+		uint simplifyAlgo = algo=="visvalingam" ? LayerDef::VISVALINGAM : 
+		                    algo=="buildings"   ? LayerDef::BUILDINGS : LayerDef::DOUGLAS_PEUCKER;
 		string source = it->value.HasMember("source") ? it->value["source"].GetString() : "";
 		vector<string> sourceColumns;
 		bool allSourceColumns = false;