Thicken spindly polygons to keep them from collapsing

Makefile

@@ -58,7 +58,7 @@ C = $(wildcard *.c) $(wildcard *.cpp)
 INCLUDES = -I/usr/local/include -I.
 LIBS = -L/usr/local/lib
 
-tippecanoe: geojson.o jsonpull/jsonpull.o tile.o pool.o mbtiles.o geometry.o projection.o memfile.o mvt.o serial.o main.o text.o dirtiles.o pmtiles_file.o plugin.o read_json.o write_json.o geobuf.o flatgeobuf.o evaluator.o geocsv.o csv.o geojson-loop.o json_logger.o visvalingam.o compression.o clip.o sort.o
+tippecanoe: geojson.o jsonpull/jsonpull.o tile.o pool.o mbtiles.o geometry.o projection.o memfile.o mvt.o serial.o main.o text.o dirtiles.o pmtiles_file.o plugin.o read_json.o write_json.o geobuf.o flatgeobuf.o evaluator.o geocsv.o csv.o geojson-loop.o json_logger.o visvalingam.o compression.o clip.o sort.o earcut.o
 	$(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread
 
 tippecanoe-enumerate: enumerate.o

earcut.cpp

@@ -0,0 +1,98 @@
+#include "geometry.hpp"
+#include "mapbox/geometry/earcut.hpp"
+
+using Coord = long long;
+using N = size_t;
+using Point = std::array<Coord, 2>;
+
+drawvec reinforce(drawvec const &pts, std::vector<std::vector<Point>> polygon, double scale) {
+	std::vector<N> indices = mapbox::earcut<N>(polygon);
+
+	drawvec out2;
+	for (size_t i = 0; i + 2 < indices.size(); i += 3) {
+		std::vector<double> lengths;
+
+		for (size_t j = 0; j < 3; j++) {
+			size_t v1 = i + j;
+			size_t v2 = i + ((j + 1) % 3);
+			size_t v3 = i + ((j + 2) % 3);
+
+#if 0
+			out2.push_back(draw(VT_MOVETO, pts[indices[v1]].x, pts[indices[v1]].y));
+			out2.push_back(draw(VT_LINETO, pts[indices[v2]].x, pts[indices[v2]].y));
+			out2.push_back(draw(VT_LINETO, pts[indices[v3]].x, pts[indices[v3]].y));
+			out2.push_back(draw(VT_LINETO, pts[indices[v1]].x, pts[indices[v1]].y));
+#endif
+
+			double px, py;
+			double d = distance_from_line_noclamp(pts[indices[v1]].x, pts[indices[v1]].y,
+							      pts[indices[v2]].x, pts[indices[v2]].y,
+							      pts[indices[v3]].x, pts[indices[v3]].y,
+							      &px, &py);
+
+			if (d < scale) {
+				double ang = atan2(py - pts[indices[v1]].y, px - pts[indices[v1]].x);
+				double dist = scale - d;
+
+				out2.push_back(draw(VT_MOVETO, pts[indices[v2]].x, pts[indices[v2]].y));
+				out2.push_back(draw(VT_LINETO, pts[indices[v2]].x + dist * cos(ang), pts[indices[v2]].y + dist * sin(ang)));
+				out2.push_back(draw(VT_LINETO, pts[indices[v3]].x + dist * cos(ang), pts[indices[v3]].y + dist * sin(ang)));
+				out2.push_back(draw(VT_LINETO, pts[indices[v3]].x, pts[indices[v3]].y));
+				out2.push_back(draw(VT_LINETO, pts[indices[v2]].x, pts[indices[v2]].y));
+			}
+		}
+	}
+
+	return out2;
+}
+
+drawvec fix_by_triangulation(drawvec const &dv, int z, int detail) {
+	std::vector<std::vector<Point>> polygon;
+	drawvec out, out2;
+	double scale = 1LL << (32 - z - detail);
+
+	for (size_t i = 0; i < dv.size(); i++) {
+		if (dv[i].op == VT_MOVETO) {
+			size_t j;
+			for (j = i + 1; j < dv.size(); j++) {
+				if (dv[j].op != VT_LINETO) {
+					break;
+				}
+			}
+
+			if (get_area(dv, i, j) > 0) {
+				// outer ring, so process whatever we have so far
+				drawvec additional = reinforce(out, polygon, scale);
+				for (auto const &d : additional) {
+					out2.push_back(d);
+				}
+
+				polygon.clear();
+				out.clear();
+			}
+
+			std::vector<Point> ring;
+			// j - 1 because earcut docs indicate that it doesn't expect
+			// a duplicate last point in each ring
+			for (size_t k = i; k < j - 1; k++) {
+				Point p = {(long long) dv[k].x, (long long) dv[k].y};
+				ring.push_back(p);
+				out.push_back(dv[k]);
+			}
+			polygon.push_back(ring);
+
+			i = j - 1;
+		}
+	}
+
+	drawvec additional = reinforce(out, polygon, scale);
+	for (auto const &d : additional) {
+		out2.push_back(d);
+	}
+
+	for (auto const &d : dv) {
+		out2.push_back(d);
+	}
+
+	return out2;
+}

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) / 3) {
+				// 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
@@ -334,7 +351,7 @@ bool point_within_tile(long long x, long long y, int z) {
 	return x >= 0 && y >= 0 && x < area && y < area;
 }
 
-double distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y) {
+double distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y, double *px, double *py) {
 	long long p2x = segB_x - segA_x;
 	long long p2y = segB_y - segA_y;
 	double something = p2x * p2x + p2y * p2y;
@@ -349,6 +366,36 @@ double distance_from_line(long long point_x, long long point_y, long long segA_x
 	double x = segA_x + u * p2x;
 	double y = segA_y + u * p2y;
 
+	if (px != NULL) {
+		*px = std::round(x);
+	}
+	if (py != NULL) {
+		*py = std::round(y);
+	}
+
+	double dx = x - point_x;
+	double dy = y - point_y;
+
+	double out = std::round(sqrt(dx * dx + dy * dy) * 16.0) / 16.0;
+	return out;
+}
+
+double distance_from_line_noclamp(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y, double *px, double *py) {
+	long long p2x = segB_x - segA_x;
+	long long p2y = segB_y - segA_y;
+	double something = p2x * p2x + p2y * p2y;
+	double u = (0 == something) ? 0 : ((point_x - segA_x) * p2x + (point_y - segA_y) * p2y) / (something);
+
+	double x = segA_x + u * p2x;
+	double y = segA_y + u * p2y;
+
+	if (px != NULL) {
+		*px = std::round(x);
+	}
+	if (py != NULL) {
+		*py = std::round(y);
+	}
+
 	double dx = x - point_x;
 	double dy = y - point_y;
 
@@ -400,7 +447,7 @@ static void douglas_peucker(drawvec &geom, int start, int n, double e, size_t ke
 		if (geom[start + first] < geom[start + second]) {
 			farthest_element_index = first;
 			for (i = first + 1; i < second; i++) {
-				double temp_dist = distance_from_line(geom[start + i].x, geom[start + i].y, geom[start + first].x, geom[start + first].y, geom[start + second].x, geom[start + second].y);
+				double temp_dist = distance_from_line(geom[start + i].x, geom[start + i].y, geom[start + first].x, geom[start + first].y, geom[start + second].x, geom[start + second].y, NULL, NULL);
 
 				double distance = std::fabs(temp_dist);
 
@@ -412,7 +459,7 @@ static void douglas_peucker(drawvec &geom, int start, int n, double e, size_t ke
 		} else {
 			farthest_element_index = second;
 			for (i = second - 1; i > first; i--) {
-				double temp_dist = distance_from_line(geom[start + i].x, geom[start + i].y, geom[start + second].x, geom[start + second].y, geom[start + first].x, geom[start + first].y);
+				double temp_dist = distance_from_line(geom[start + i].x, geom[start + i].y, geom[start + second].x, geom[start + second].y, geom[start + first].x, geom[start + first].y, NULL, NULL);
 
 				double distance = std::fabs(temp_dist);
 
@@ -1046,7 +1093,7 @@ double label_goodness(const drawvec &dv, long long x, long long y) {
 		}
 
 		if (i > 0 && dv[i].op == VT_LINETO) {
-			dist = distance_from_line(x, y, dv[i - 1].x, dv[i - 1].y, dv[i].x, dv[i].y);
+			dist = distance_from_line(x, y, dv[i - 1].x, dv[i - 1].y, dv[i].x, dv[i].y, NULL, NULL);
 			if (dist < closest) {
 				closest = dist;
 			}

geometry.hpp

@@ -96,7 +96,8 @@ drawvec clip_lines(drawvec &geom, long long x1, long long y1, long long x2, long
 drawvec clip_point(drawvec &geom, long long x1, long long y1, long long x2, long long y2);
 void visvalingam(drawvec &ls, size_t start, size_t end, double threshold, size_t retain);
 int pnpoly(const drawvec &vert, size_t start, size_t nvert, long long testx, long long testy);
-double distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y);
+double distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y, double *px, double *py);
+double distance_from_line_noclamp(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y, double *px, double *py);
 
 std::string overzoom(mvt_tile tile, 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,
@@ -106,4 +107,6 @@ std::string overzoom(std::string s, int oz, int ox, int oy, int nz, int nx, int
 		     int detail, int buffer, std::set<std::string> const &keep, bool do_compress,
 		     std::vector<std::pair<unsigned, unsigned>> *next_overzoomed_tiles);
 
+drawvec fix_by_triangulation(drawvec const &dv, int z, int detail);
+
 #endif