Add a critical point where the set of rings using a polygon edge changes

geometry.cpp

@@ -1042,17 +1042,19 @@ drawvec impose_tile_boundaries(drawvec &geom, long long extent) {
 	return out;
 }
 
-drawvec simplify_lines(drawvec &geom, int z, int detail, bool mark_tile_bounds, double simplification) {
+drawvec simplify_lines(drawvec &geom, int z, int detail, bool mark_tile_bounds, double simplification, bool already_marked) {
 	int res = 1 << (32 - detail - z);
 	long long area = 1LL << (32 - z);
 
-	for (size_t i = 0; i < geom.size(); i++) {
-		if (geom[i].op == VT_MOVETO) {
-			geom[i].necessary = 1;
-		} else if (geom[i].op == VT_LINETO) {
-			geom[i].necessary = 0;
-		} else {
-			geom[i].necessary = 1;
+	if (!already_marked) {
+		for (size_t i = 0; i < geom.size(); i++) {
+			if (geom[i].op == VT_MOVETO) {
+				geom[i].necessary = 1;
+			} else if (geom[i].op == VT_LINETO) {
+				geom[i].necessary = 0;
+			} else {
+				geom[i].necessary = 1;
+			}
 		}
 	}
 
@@ -1073,7 +1075,19 @@ drawvec simplify_lines(drawvec &geom, int z, int detail, bool mark_tile_bounds,
 			geom[j - 1].necessary = 1;
 
 			if (j - i > 1) {
-				douglas_peucker(geom, i, j - i, res * simplification);
+				if (already_marked && geom[j - 1] < geom[i]) {
+					drawvec dv;
+					for (size_t k = j; k > i; k--) {
+						dv.push_back(geom[k - 1]);
+					}
+					douglas_peucker(dv, 0, j - i, res * simplification);
+					size_t l = 0;
+					for (size_t k = j; k > i; k--) {
+						geom[k - 1] = dv[l++];
+					}
+				} else {
+					douglas_peucker(geom, i, j - i, res * simplification);
+				}
 			}
 			i = j - 1;
 		}

geometry.hpp

@@ -28,6 +28,22 @@ struct draw {
 
 	draw() {
 	}
+
+	bool operator<(draw const &s) const {
+		if (y < s.y || (y == s.y && x < s.x)) {
+			return true;
+		} else {
+			return false;
+		}
+	}
+
+	bool operator==(draw const &s) const {
+		return y == s.y && x == s.x;
+	}
+
+	bool operator!=(draw const &s) const {
+		return y != s.y || x != s.x;
+	}
 };
 
 typedef std::vector<draw> drawvec;
@@ -43,7 +59,7 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double
 drawvec clip_lines(drawvec &geom, int z, int detail, long long buffer);
 bool point_within_tile(long long x, long long y, int z, int detail, long long buffer);
 int quick_check(long long *bbox, int z, int detail, long long buffer);
-drawvec simplify_lines(drawvec &geom, int z, int detail, bool mark_tile_bounds, double simplification);
+drawvec simplify_lines(drawvec &geom, int z, int detail, bool mark_tile_bounds, double simplification, bool already_marked);
 drawvec reorder_lines(drawvec &geom);
 drawvec fix_polygon(drawvec &geom);
 std::vector<drawvec> chop_polygon(std::vector<drawvec> &geoms);

main.cpp

@@ -1801,6 +1801,7 @@ int main(int argc, char **argv) {
 		{"drop-polygons", no_argument, &additional[A_POLYGON_DROP], 1},
 		{"prefer-radix-sort", no_argument, &additional[A_PREFER_RADIX_SORT], 1},
 		{"calculate-feature-density", no_argument, &additional[A_CALCULATE_FEATURE_DENSITY], 1},
+		{"simplify-polygons-together", no_argument, &additional[A_SIMPLIFY_TOGETHER], 1},
 
 		{"no-line-simplification", no_argument, &prevent[P_SIMPLIFY], 1},
 		{"simplify-only-low-zooms", no_argument, &prevent[P_SIMPLIFY_LOW], 1},

options.hpp

@@ -4,6 +4,7 @@
 #define A_LINE_DROP ((int) 'l')
 #define A_DEBUG_POLYGON ((int) 'd')
 #define A_POLYGON_DROP ((int) 'p')
+#define A_SIMPLIFY_TOGETHER ((int) 't')
 #define A_PREFER_RADIX_SORT ((int) 'R')
 #define A_CALCULATE_FEATURE_DENSITY ((int) 'g')
 

tile.cpp

@@ -460,7 +460,11 @@ void *partial_feature_worker(void *v) {
 					geom = remove_noop(geom, t, 32 - z - line_detail);
 				}
 
-				drawvec ngeom = simplify_lines(geom, z, line_detail, !(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), (*partials)[i].simplification);
+				bool already_marked = false;
+				if (additional[A_SIMPLIFY_TOGETHER] && t == VT_POLYGON) {
+					already_marked = true;
+				}
+				drawvec ngeom = simplify_lines(geom, z, line_detail, !(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), (*partials)[i].simplification, already_marked);
 
 				if (t != VT_POLYGON || ngeom.size() >= 3) {
 					geom = ngeom;
@@ -561,6 +565,117 @@ int manage_gap(unsigned long long index, unsigned long long *previndex, double s
 	return 0;
 }
 
+void find_common_edges(std::vector<partial> &partials) {
+	std::map<drawvec, std::set<size_t>> edges;
+
+	// Construct a mapping from all polygon edges to the set of rings
+	// that each edge appears in. (The ring number is across all polygons;
+	// we don't need to look it back up, just to tell where it changes.)
+
+	size_t ring = 0;
+	for (size_t i = 0; i < partials.size(); i++) {
+		if (partials[i].t == VT_POLYGON) {
+			for (size_t j = 0; j < partials[i].geoms.size(); j++) {
+				for (size_t k = 0; k + 1 < partials[i].geoms[j].size(); k++) {
+					if (partials[i].geoms[j][k].op == VT_MOVETO) {
+						ring++;
+					}
+
+					if (partials[i].geoms[j][k + 1].op == VT_LINETO) {
+						drawvec dv;
+						if (partials[i].geoms[j][k] < partials[i].geoms[j][k + 1]) {
+							dv.push_back(partials[i].geoms[j][k]);
+							dv.push_back(partials[i].geoms[j][k + 1]);
+						} else {
+							dv.push_back(partials[i].geoms[j][k + 1]);
+							dv.push_back(partials[i].geoms[j][k]);
+						}
+
+						auto e = edges.find(dv);
+						if (e != edges.end()) {
+							e->second.insert(ring);
+						} else {
+							std::set<size_t> s;
+							s.insert(ring);
+							edges.insert(std::pair<drawvec, std::set<size_t>>(dv, s));
+						}
+					}
+				}
+			}
+		}
+	}
+
+	// Now mark all the points where the set of rings using the edge on one side
+	// is not the same as the set of rings using the edge on the other side.
+
+	for (size_t i = 0; i < partials.size(); i++) {
+		if (partials[i].t == VT_POLYGON) {
+			for (size_t j = 0; j < partials[i].geoms.size(); j++) {
+				// following simplify_lines()
+				for (size_t k = 0; k < partials[i].geoms[j].size(); k++) {
+					if (partials[i].geoms[j][k].op == VT_MOVETO) {
+						partials[i].geoms[j][k].necessary = 1;
+					} else if (partials[i].geoms[j][k].op == VT_LINETO) {
+						partials[i].geoms[j][k].necessary = 0;
+					} else {
+						partials[i].geoms[j][k].necessary = 1;
+					}
+				}
+
+				for (size_t k = 1; k + 1 < partials[i].geoms[j].size(); k++) {
+					if (partials[i].geoms[j][k].op == VT_LINETO && partials[i].geoms[j][k + 1].op == VT_LINETO) {
+						drawvec left, right;
+						if (partials[i].geoms[j][k - 1] < partials[i].geoms[j][k]) {
+							left.push_back(partials[i].geoms[j][k - 1]);
+							left.push_back(partials[i].geoms[j][k]);
+						} else {
+							left.push_back(partials[i].geoms[j][k]);
+							left.push_back(partials[i].geoms[j][k - 1]);
+						}
+
+						if (partials[i].geoms[j][k] < partials[i].geoms[j][k + 1]) {
+							right.push_back(partials[i].geoms[j][k]);
+							right.push_back(partials[i].geoms[j][k + 1]);
+						} else {
+							right.push_back(partials[i].geoms[j][k + 1]);
+							right.push_back(partials[i].geoms[j][k]);
+						}
+
+						auto e1 = edges.find(left);
+						auto e2 = edges.find(right);
+
+						if (left[1] < left[0]) {
+							fprintf(stderr, "left misordered\n");
+						}
+						if (right[1] < right[0]) {
+							fprintf(stderr, "left misordered\n");
+						}
+
+						if (e1 == edges.end() || e2 == edges.end()) {
+							fprintf(stderr, "Internal error: polygon edge lookup failed for %lld,%lld to %lld,%lld or %lld,%lld to %lld,%lld\n", left[0].x, left[0].y, left[1].x, left[1].y, right[0].x, right[0].y, right[1].x, right[1].y);
+
+							for (auto ei = edges.begin(); ei != edges.end(); ++ei) {
+								if (ei->first[1] < ei->first[0]) {
+									fprintf(stderr, "%lld,%lld to %lld,%lld %lu\n",
+										ei->first[0].x, ei->first[0].y,
+										ei->first[1].x, ei->first[1].y,
+										ei->second.size());
+								}
+							}
+
+							exit(EXIT_FAILURE);
+						}
+
+						if (e1->second != e2->second) {
+							partials[i].geoms[j][k].necessary = 1;
+						}
+					}
+				}
+			}
+		}
+	}
+}
+
 long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *stringpool, int z, unsigned tx, unsigned ty, int detail, int min_detail, int basezoom, sqlite3 *outdb, double droprate, int buffer, const char *fname, FILE **geomfile, int minzoom, int maxzoom, double todo, volatile long long *along, long long alongminus, double gamma, int child_shards, long long *meta_off, long long *pool_off, unsigned *initial_x, unsigned *initial_y, volatile int *running, double simplification, std::vector<std::map<std::string, layermap_entry>> *layermaps, std::vector<std::vector<std::string>> *layer_unmaps) {
 	int line_detail;
 	double fraction = 1;
@@ -868,6 +983,10 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 			}
 		}
 
+		if (additional[A_SIMPLIFY_TOGETHER]) {
+			find_common_edges(partials);
+		}
+
 		int tasks = ceil((double) CPUS / *running);
 		if (tasks < 1) {
 			tasks = 1;
@@ -987,7 +1106,7 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 				if (layer_features[x].coalesced && layer_features[x].type == VT_LINE) {
 					layer_features[x].geom = remove_noop(layer_features[x].geom, layer_features[x].type, 0);
 					layer_features[x].geom = simplify_lines(layer_features[x].geom, 32, 0,
-										!(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), simplification);
+										!(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), simplification, false);
 				}
 
 				if (layer_features[x].type == VT_POLYGON) {