From d43526ee88f9a673d914037153f46e507347c4ac Mon Sep 17 00:00:00 2001
From: systemed <richard@systemeD.net>
Date: Fri, 24 Mar 2023 16:00:10 +0000
Subject: [PATCH] Use faster bbox-specialised clipping algorithm

---
 include/geom.h        |  5 +++
 src/geom.cpp          | 77 +++++++++++++++++++++++++++++++++++++++++++
 src/output_object.cpp | 22 +++----------
 3 files changed, 87 insertions(+), 17 deletions(-)

diff --git a/include/geom.h b/include/geom.h
index 41bfa9e5..f946927f 100644
--- a/include/geom.h
+++ b/include/geom.h
@@ -89,5 +89,10 @@ void make_valid(GeometryT &geom) { }
 
 void make_valid(MultiPolygon &mp);
 
+Point intersect_edge(Point const &a, Point const &b, char edge, Box const &bbox);
+char bit_code(Point const &p, Box const &bbox);
+void fast_clip(Ring &points, Box const &bbox);
+void fast_clip(MultiPolygon &mp, Box const &bbox);
+
 #endif //_GEOM_TYPES_H
 
diff --git a/src/geom.cpp b/src/geom.cpp
index eb921830..25112f4a 100644
--- a/src/geom.cpp
+++ b/src/geom.cpp
@@ -142,3 +142,80 @@ void make_valid(MultiPolygon &mp)
 	}
 	mp = result;
 }
+
+// ---------------
+// Sutherland-Hodgeman clipper
+// ported from Volodymyr Agafonkin's https://github.com/mapbox/lineclip
+
+// intersect a segment against one of the 4 lines that make up the bbox
+Point intersect_edge(Point const &a, Point const &b, char edge, Box const &bbox) {
+	if (edge & 8) return Point(a.x() + (b.x() - a.x()) * (bbox.max_corner().y() - a.y()) / (b.y() - a.y()), bbox.max_corner().y()); // top
+	if (edge & 4) return Point(a.x() + (b.x() - a.x()) * (bbox.min_corner().y() - a.y()) / (b.y() - a.y()), bbox.min_corner().y()); // bottom
+	if (edge & 2) return Point(bbox.max_corner().x(), a.y() + (b.y() - a.y()) * (bbox.max_corner().x() - a.x()) / (b.x() - a.x())); // right
+	if (edge & 1) return Point(bbox.min_corner().x(), a.y() + (b.y() - a.y()) * (bbox.min_corner().x() - a.x()) / (b.x() - a.x())); // left
+	throw std::runtime_error("intersect called on non-intersection");
+}
+
+// bit code reflects the point position relative to the bbox:
+//         left  mid  right
+//    top  1001  1000  1010
+//    mid  0001  0000  0010
+// bottom  0101  0100  0110
+
+char bit_code(Point const &p, Box const &bbox) {
+	char code = 0;
+	if      (p.x() < bbox.min_corner().x()) code |= 1; // left
+	else if (p.x() > bbox.max_corner().x()) code |= 2; // right
+	if      (p.y() < bbox.min_corner().y()) code |= 4; // bottom
+	else if (p.y() > bbox.max_corner().y()) code |= 8; // top
+	return code;
+}
+
+// Sutherland-Hodgeman polygon clipping algorithm
+void fast_clip(Ring &points, Box const &bbox) {
+	// clip against each side of the clip rectangle
+	for (char edge = 1; edge <= 8; edge *= 2) {
+		Ring result;
+		Point prev = points[points.size() - 1];
+		bool prevInside = (bit_code(prev, bbox) & edge)==0;
+
+		for (unsigned int i = 0; i<points.size(); i++) {
+			Point p = points[i];
+			bool inside = (bit_code(p, bbox) & edge)==0;
+
+			// if segment goes through the clip window, add an intersection
+			if (inside!=prevInside) result.emplace_back(intersect_edge(prev, p, edge, bbox));
+			if (inside) result.emplace_back(p); // add a point if it's inside
+			prev = p;
+			prevInside = inside;
+		}
+		points = std::move(result);
+		if (points.size()==0) break;
+	}
+}
+
+// Wrappers for polygon/multipolygon
+void fast_clip(Polygon &polygon, Box const &bbox) {
+	fast_clip(polygon.outer(), bbox);
+	if (polygon.outer().empty()) {
+		polygon.inners().resize(0);
+		return;
+	}
+	for (auto &inner: polygon.inners()) {
+		fast_clip(inner, bbox);
+	}
+	polygon.inners().erase(std::remove_if(
+		polygon.inners().begin(), polygon.inners().end(), 
+		[](const Ring &inner) -> bool { return inner.empty(); }),
+		polygon.inners().end());
+}
+
+void fast_clip(MultiPolygon &mp, Box const &bbox) {
+	for (auto &polygon: mp) {
+		fast_clip(polygon, bbox);
+	}
+	mp.erase(std::remove_if(
+		mp.begin(), mp.end(), 
+		[](const Polygon &poly) -> bool { return poly.outer().empty(); }),
+		mp.end());
+}
diff --git a/src/output_object.cpp b/src/output_object.cpp
index 275d60d2..f9feade7 100644
--- a/src/output_object.cpp
+++ b/src/output_object.cpp
@@ -167,23 +167,11 @@ Geometry buildWayGeometry(OSMStore &osmStore, OutputObject const &oo, const Tile
 					std::min(box.max_corner().y(), extBox.max_corner().y()));
 			}
 
-			Polygon clippingPolygon;
-			geom::convert(box, clippingPolygon);
-	
-			try {
-				MultiPolygon mp;
-				if (geom::within(input, clippingPolygon)) { geom::assign(mp, input); return mp; }
-				// Work around boost::geometry intersection issues by doing each constituent polygon at a time
-				for (auto &p : input) {
-					MultiPolygon out;
-					geom::intersection(p, clippingPolygon, out);
-					for (auto &o : out) mp.emplace_back(move(o));
-				}
-				return mp;
-			} catch (geom::overlay_invalid_input_exception &err) {
-				std::cout << "Couldn't clip polygon (self-intersection)" << std::endl;
-				return MultiPolygon(); // blank
-			}
+			MultiPolygon mp;
+			geom::assign(mp, input);
+			fast_clip(mp, box);
+			geom::correct(mp);
+			return mp;
 		}
 
 		default: