Detect longitude wraparound only if it reduces self-intersections

geometry.cpp

@@ -1397,3 +1397,49 @@ drawvec checkerboard_anchors(drawvec const &geom, int tx, int ty, int z, unsigne
 
 	return out;
 }
+
+const std::pair<double, double> SAME_SLOPE = std::make_pair(-INT_MAX, INT_MAX);
+
+// https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
+//
+// beware of
+// https://stackoverflow.com/questions/9043805/test-if-two-lines-intersect-javascript-function/16725715#16725715
+// which does not seem to produce correct results.
+std::pair<double, double> get_line_intersection(long long p0_x, long long p0_y, long long p1_x, long long p1_y,
+						long long p2_x, long long p2_y, long long p3_x, long long p3_y) {
+	// bounding box reject, x
+	long long min01x = std::min(p0_x, p1_x);
+	long long max01x = std::max(p0_x, p1_x);
+	long long min23x = std::min(p2_x, p3_x);
+	long long max23x = std::max(p2_x, p3_x);
+	if (max01x < min23x || max23x < min01x) {
+		return std::make_pair(-1, -1);
+	}
+
+	// bounding box reject, y
+	long long min01y = std::min(p0_y, p1_y);
+	long long max01y = std::max(p0_y, p1_y);
+	long long min23y = std::min(p2_y, p3_y);
+	long long max23y = std::max(p2_y, p3_y);
+	if (max01y < min23y || max23y < min01y) {
+		return std::make_pair(-1, -1);
+	}
+
+	long long d01_x, d01_y, d23_x, d23_y;
+	d01_x = p1_x - p0_x;
+	d01_y = p1_y - p0_y;
+	d23_x = p3_x - p2_x;
+	d23_y = p3_y - p2_y;
+
+	long long det = (-d23_x * d01_y + d01_x * d23_y);
+
+	if (det != 0) {
+		double t, s;
+		t = (d23_x * (p0_y - p2_y) - d23_y * (p0_x - p2_x)) / (double) det;
+		s = (-d01_y * (p0_x - p2_x) + d01_x * (p0_y - p2_y)) / (double) det;
+
+		return std::make_pair(t, s);
+	}
+
+	return SAME_SLOPE;
+}

geometry.hpp

@@ -114,4 +114,7 @@ std::string overzoom(const std::string &s, int oz, int ox, int oy, int nz, int n
 		     std::unordered_map<std::string, attribute_op> const &attribute_accum,
 		     std::vector<std::string> const &unidecode_data);
 
+std::pair<double, double> get_line_intersection(long long p0_x, long long p0_y, long long p1_x, long long p1_y,
+						long long p2_x, long long p2_y, long long p3_x, long long p3_y);
+
 #endif

serial.cpp

@@ -298,9 +298,71 @@ serial_feature deserialize_feature(std::string const &geoms, unsigned z, unsigne
 }
 
 static long long scale_geometry(struct serialization_state *sst, long long *bbox, drawvec &geom) {
-	long long offset = 0;
+	if (additional[A_DETECT_WRAPAROUND]) {
-	long long prev = 0;
+		drawvec adjusted = geom;
-	bool has_prev = false;
+		std::vector<size_t> jumps;
+
+		long long offset = 0;
+		long long prev = 0;
+		bool has_prev = false;
+
+		for (size_t i = 0; i < geom.size(); i++) {
+			if (geom[i].op == VT_LINETO) {
+				long long x = geom[i].x;
+				x += offset;
+
+				if (has_prev) {
+					// jumps at least 180° in either direction
+
+					if (x - prev > (1LL << 31)) {
+						offset -= 1LL << 32;
+						x -= 1LL << 32;
+						jumps.push_back(i);
+					} else if (prev - x > (1LL << 31)) {
+						offset += 1LL << 32;
+						x += 1LL << 32;
+						jumps.push_back(i);
+					}
+				}
+
+				adjusted[i].x = x;
+				has_prev = true;
+				prev = x;
+			} else {
+				offset = 0;
+				prev = geom[i].x;
+			}
+		}
+
+		// count whether making these adjustments resulted in
+		// more or less self-intersections than in the original
+
+		size_t orig_bad = 0;
+		size_t adjusted_bad = 0;
+
+		for (auto i : jumps) {
+			for (size_t j = 1; j < geom.size(); j++) {
+				if (geom[j].op == VT_LINETO && j != i) {
+					auto before = get_line_intersection(geom[i - 1].x, geom[i - 1].y, geom[i].x, geom[i].y,
+									    geom[j - 1].x, geom[j - 1].y, geom[j].x, geom[j].y);
+					auto after = get_line_intersection(adjusted[i - 1].x, adjusted[i - 1].y, adjusted[i].x, adjusted[i].y,
+									   adjusted[j - 1].x, adjusted[j - 1].y, adjusted[j].x, adjusted[j].y);
+
+					if (before.first > 0 && before.first < 1) {
+						orig_bad++;
+					}
+					if (after.first > 0 && after.first < 1) {
+						adjusted_bad++;
+					}
+				}
+			}
+		}
+
+		if (orig_bad > adjusted_bad) {
+			geom = std::move(adjusted);
+		}
+	}
+
 	double scale = 1.0 / (1 << geometry_scale);
 
 	for (size_t i = 0; i < geom.size(); i++) {
@@ -308,30 +370,6 @@ static long long scale_geometry(struct serialization_state *sst, long long *bbox
 			long long x = geom[i].x;
 			long long y = geom[i].y;
 
-			if (additional[A_DETECT_WRAPAROUND]) {
-				if (geom[i].op == VT_LINETO) {
-					x += offset;
-					if (has_prev) {
-						// jumps at least 180° but not exactly 360°,
-						// which in some data sets is an intentional
-						// line across the world
-						if (x - prev > (1LL << 31) && x - prev != (1LL << 32)) {
-							offset -= 1LL << 32;
-							x -= 1LL << 32;
-						} else if (prev - x > (1LL << 31) && prev - x != (1LL << 32)) {
-							offset += 1LL << 32;
-							x += 1LL << 32;
-						}
-					}
-
-					has_prev = true;
-					prev = x;
-				} else {
-					offset = 0;
-					prev = x;
-				}
-			}
-
 			if (x < bbox[0]) {
 				bbox[0] = x;
 			}