Add --buffer-polygons-outward

README.md

@@ -496,6 +496,7 @@ the same layer, enclose them in an `all` expression so they will all be evaluate
  * `-pt` or `--no-tiny-polygon-reduction`: Don't combine the area of very small polygons into small squares that represent their combined area.
  * `-pT` or `--no-tiny-polygon-reduction-at-maximum-zoom`: Combine the area of very small polygons into small squares that represent their combined area only at zoom levels below the maximum.
  * `--tiny-polygon-size=`_size_: Use the specified _size_ for tiny polygons instead of the default 2. Anything above 6 or so will lead to visible artifacts with the default tile detail.
+ * `-aO` or `--buffer-polygons-outward`: Buffer polygons, except those that have been detected to be continous by `--no-simplification-of-shared-nodes`, outward slightly to prevent narrow polygon spindles fom collapsing away.
  * `-av` or `--visvalingam`: Use Visvalingam's simplification algorithm rather than Douglas-Peucker's.
 
 ### Attempts to improve shared polygon boundaries

geometry.cpp

@@ -197,6 +197,23 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *still_needs_sim
 
 			double area = get_area(geom, i, j);
 
+			long long minx = LLONG_MAX;
+			long long maxx = LLONG_MIN;
+			long long miny = LLONG_MAX;
+			long long maxy = LLONG_MIN;
+			for (auto const &d : geom) {
+				minx = std::min(minx, (long long) d.x);
+				maxx = std::max(maxx, (long long) d.x);
+				miny = std::min(miny, (long long) d.y);
+				maxy = std::max(maxy, (long long) d.y);
+			}
+			if (area > 0 && area <= pixel * pixel && area < (maxx - minx) * (maxy - miny) / 5) {
+				// if the polygon doesn't use most of its area,
+				// don't let it be dust, because the shape is
+				// probably something weird and interesting.
+				area = pixel * pixel * 2;
+			}
+
 			// XXX There is an ambiguity here: If the area of a ring is 0 and it is followed by holes,
 			// we don't know whether the area-0 ring was a hole too or whether it was the outer ring
 			// that these subsequent holes are somehow being subtracted from. I hope that if a polygon
@@ -485,6 +502,26 @@ drawvec impose_tile_boundaries(drawvec &geom, long long extent) {
 	return out;
 }
 
+bool is_shared_node(draw d, int z, int tx, int ty, struct node *shared_nodes_map, size_t nodepos) {
+	if (nodepos > 0) {
+		// offset to global
+		if (z != 0) {
+			d.x += tx * (1LL << (32 - z));
+			d.y += ty * (1LL << (32 - z));
+		}
+
+		// to quadkey
+		struct node n;
+		n.index = encode_quadkey((unsigned) d.x, (unsigned) d.y);
+
+		if (bsearch(&n, shared_nodes_map, nodepos / sizeof(node), sizeof(node), nodecmp) != NULL) {
+			return true;
+		}
+	}
+
+	return false;
+}
+
 drawvec simplify_lines(drawvec &geom, int z, int tx, int ty, int detail, bool mark_tile_bounds, double simplification, size_t retain, drawvec const &shared_nodes, struct node *shared_nodes_map, size_t nodepos) {
 	int res = 1 << (32 - detail - z);
 	long long area = 1LL << (32 - z);
@@ -514,21 +551,8 @@ drawvec simplify_lines(drawvec &geom, int z, int tx, int ty, int detail, bool ma
 				geom[i].necessary = true;
 			}
 
-			if (nodepos > 0) {
-				// offset to global
-				draw d = geom[i];
-				if (z != 0) {
-					d.x += tx * (1LL << (32 - z));
-					d.y += ty * (1LL << (32 - z));
-				}
-
-				// to quadkey
-				struct node n;
-				n.index = encode_quadkey((unsigned) d.x, (unsigned) d.y);
-
-				if (bsearch(&n, shared_nodes_map, nodepos / sizeof(node), sizeof(node), nodecmp) != NULL) {
-					geom[i].necessary = true;
-				}
+			if (is_shared_node(geom[i], z, tx, ty, shared_nodes_map, nodepos)) {
+				geom[i].necessary = true;
 			}
 		}
 	}
@@ -1399,33 +1423,24 @@ drawvec checkerboard_anchors(drawvec const &geom, int tx, int ty, int z, unsigne
 	return out;
 }
 
-static double anglediff(double dx1, double dy1, double dx2, double dy2) {
-	// https://gist.github.com/e-n-f/417cd85f6d9b00dab887276376e77f30
-	// https://twitter.com/enf/status/795798781186830341
+bool is_continous_poly(drawvec const &geom, size_t i, size_t j, int z, int tx, int ty, struct node *shared_nodes_map, size_t nodepos) {
+	size_t count = 0;
 
-	double cross = dx1 * dy2 - dy1 * dx2;
-	double sd = sqrt(dx1 * dx1 + dy1 * dy1);
-	double cd = sqrt(dx2 * dx2 + dy2 * dy2);
+	// j - 1 to avoid double-counting the duplicate last node
+	for (; i < j - 1; i++) {
+		if (is_shared_node(geom[i], z, tx, ty, shared_nodes_map, nodepos)) {
+			count++;
 
-	double dot = dx1 * dx2 + dy1 * dy2;
-	double ret;
-	if (dot < 0) {
-		if (cross < 0) {
-			ret = -M_PI - asin(cross / (sd * cd));
-		} else {
-			ret = M_PI - asin(cross / (sd * cd));
+			if (count >= 4) {
+				return true;
+			}
 		}
-	} else {
-		ret = asin(cross / (sd * cd));
 	}
-	if (ret < 0) {
-		// make it 0 to 360, not -180 to 180
-		ret += 2 * M_PI;
-	}
-	return ret;
+
+	return false;
 }
 
-drawvec buffer_poly(drawvec const &geom, double buffer) {
+drawvec buffer_poly(drawvec const &geom, double buffer, int z, int tx, int ty, struct node *shared_nodes_map, size_t nodepos) {
 	drawvec out = geom;
 
 	if (buffer != 0) {
@@ -1438,38 +1453,33 @@ drawvec buffer_poly(drawvec const &geom, double buffer) {
 					}
 				}
 
-				for (size_t k = i; k < j - 1; k++) {
-					draw p0 = geom[(k + 0 - i) % (j - i - 1) + i];
-					draw p1 = geom[(k + 1 - i) % (j - i - 1) + i];
-					draw p2 = geom[(k + 2 - i) % (j - i - 1) + i];
-
-					double adiff = 2 * M_PI - anglediff(p0.x - p1.x, p0.y - p1.y,
-									    p2.x - p1.x, p2.y - p1.y);
+				if (!is_continous_poly(geom, i, j, z, tx, ty, shared_nodes_map, nodepos)) {
+					for (size_t k = i; k < j - 1; k++) {
+						draw p0 = geom[(k + 0 - i) % (j - i - 1) + i];
+						draw p1 = geom[(k + 1 - i) % (j - i - 1) + i];
+						draw p2 = geom[(k + 2 - i) % (j - i - 1) + i];
 
-					double a1 = atan2(p1.y - p0.y, p1.x - p0.x);
+						double a10 = atan2(p1.y - p0.y, p1.x - p0.x);
+						double a21 = atan2(p2.y - p1.y, p2.x - p1.x);
 
-					printf("at %zu %zu %zu ",
-					       (k + 0 - i) % (j - i - 1) + i,
-					       (k + 1 - i) % (j - i - 1) + i,
-					       (k + 2 - i) % (j - i - 1) + i);
-					printf("%lld,%lld to %lld,%lld to %lld,%lld: ", p0.x, p0.y, p1.x, p1.y, p2.x, p2.y);
-					printf("angle %f then %f\n",
-					       atan2(p1.y - p0.y, p1.x - p0.x) * 180 / M_PI,
-					       atan2(p2.y - p1.y, p2.x - p1.x) * 180 / M_PI);
+						double dx = cos(a10 - 90 * M_PI / 180) + cos(a21 - 90 * M_PI / 180);
+						double dy = sin(a10 - 90 * M_PI / 180) + sin(a21 - 90 * M_PI / 180);
 
-					double a10 = atan2(p1.y - p0.y, p1.x - p0.x);
-					double a21 = atan2(p2.y - p1.y, p2.x - p1.x);
+						// the angle halfway between the angles
+						// perpendicular to a0->a1 (a10) and a1->a2 (a21)
+						double a2 = atan2(dy, dx);
 
-					double dx = cos(a10 - 90 * M_PI / 180) + cos(a21 - 90 * M_PI / 180);
-					double dy = sin(a10 - 90 * M_PI / 180) + sin(a21 - 90 * M_PI / 180);
-					double a2 = atan2(dy, dx);
+						out[(k + 1 - i) % (j - i - 1) + i].x = std::round(p1.x + buffer * cos(a2));
+						out[(k + 1 - i) % (j - i - 1) + i].y = std::round(p1.y + buffer * sin(a2));
+					}
 
-					out[(k + 1 - i) % (j - i - 1) + i].x = std::round(p1.x + buffer * cos(a2));
-					out[(k + 1 - i) % (j - i - 1) + i].y = std::round(p1.y + buffer * sin(a2));
+					out[j - 1].x = out[i].x;
+					out[j - 1].y = out[i].y;
+				} else {
+					printf("skipping continuous ring %zu to %zu\n", i, j);
 				}
 
-				out[j - 1].x = out[i].x;
-				out[j - 1].y = out[i].y;
+				i = j - 1;
 			}
 		}
 	}

geometry.hpp

@@ -104,6 +104,6 @@ std::string overzoom(mvt_tile tile, int oz, int ox, int oy, int nz, int nx, int
 std::string overzoom(std::string s, int oz, int ox, int oy, int nz, int nx, int ny,
 		     int detail, int buffer, std::set<std::string> const &keep, bool do_compress);
 
-drawvec buffer_poly(drawvec const &geom, double buffer);
+drawvec buffer_poly(drawvec const &geom, double buffer, int z, int tx, int ty, struct node *shared_nodes_map, size_t nodepos);
 
 #endif

main.cpp

@@ -3181,6 +3181,7 @@ int main(int argc, char **argv) {
 		{"no-tiny-polygon-reduction", no_argument, &prevent[P_TINY_POLYGON_REDUCTION], 1},
 		{"no-tiny-polygon-reduction-at-maximum-zoom", no_argument, &prevent[P_TINY_POLYGON_REDUCTION_AT_MAXZOOM], 1},
 		{"tiny-polygon-size", required_argument, 0, '~'},
+		{"buffer-polygons-outward", no_argument, &additional[A_BUFFER_POLYGONS_OUTWARD], 1},
 		{"no-simplification-of-shared-nodes", no_argument, &prevent[P_SIMPLIFY_SHARED_NODES], 1},
 		{"visvalingam", no_argument, &additional[A_VISVALINGAM], 1},
 

man/tippecanoe.1

@@ -637,6 +637,8 @@ the line or polygon within one tile unit of its proper location. You can probabl
 .IP \(bu 2
 \fB\fC\-\-tiny\-polygon\-size=\fR\fIsize\fP: Use the specified \fIsize\fP for tiny polygons instead of the default 2. Anything above 6 or so will lead to visible artifacts with the default tile detail.
 .IP \(bu 2
+\fB\fC\-aO\fR or \fB\fC\-\-buffer\-polygons\-outward\fR: Buffer polygons, except those that have been detected to be continous by \fB\fC\-\-no\-simplification\-of\-shared\-nodes\fR, outward slightly to prevent narrow polygon spindles fom collapsing away.
+.IP \(bu 2
 \fB\fC\-av\fR or \fB\fC\-\-visvalingam\fR: Use Visvalingam's simplification algorithm rather than Douglas\-Peucker's.
 .RE
 .SS Attempts to improve shared polygon boundaries

options.hpp

@@ -26,6 +26,7 @@
 #define A_HILBERT ((int) 'h')
 #define A_VISVALINGAM ((int) 'v')
 #define A_GENERATE_POLYGON_LABEL_POINTS ((int) 'P')
+#define A_BUFFER_POLYGONS_OUTWARD ((int) 'O')
 
 #define P_SIMPLIFY ((int) 's')
 #define P_SIMPLIFY_LOW ((int) 'S')

tile.cpp

@@ -612,7 +612,10 @@ void *partial_feature_worker(void *v) {
 		int out_detail = (*partials)[i].extra_detail;
 
 		drawvec geom = (*partials)[i].geoms[0];
-		drawvec before_scaling = geom;
+
+		if (additional[A_BUFFER_POLYGONS_OUTWARD] && t == VT_POLYGON) {
+			geom = buffer_poly(geom, (1LL << (32 - z - out_detail)) * sqrt(2) / 2, z, (*partials)[i].tx, (*partials)[i].ty, a->shared_nodes_map, a->nodepos);
+		}
 
 		to_tile_scale(geom, z, out_detail);
 
@@ -627,19 +630,10 @@ void *partial_feature_worker(void *v) {
 					check_polygon(geom);
 				}
 
-				if (true || geom.size() < 4) {
+				if (geom.size() < 4) {
 					if (area > 0) {
-						// Try reviving the polygon by buffering it outward a little
-
-						geom = buffer_poly(before_scaling, (1LL << (32 - z - out_detail)) * 2.0);
-						to_tile_scale(geom, z, out_detail);
-						geom = clean_or_clip_poly(geom, 0, 0, false, true);
-
-						if (geom.size() < 4) {
-							// OK, that didn't work, make a placeholder
-							// area is in world coordinates, calculated before scaling down
-							geom = revive_polygon(before, area, z, out_detail);
-						}
+						// area is in world coordinates, calculated before scaling down
+						geom = revive_polygon(before, area, z, out_detail);
 					} else {
 						geom.clear();
 					}