Externalizing polygon shard detection (#146)

* Start of externalizing polygon shard detection

* Completely untested external quicksort

* Add unit test for external quicksort

* Remember to clean up temporary files

* Sort and scan the vertices

* Bring over more vertex logic

* Make nodes from vertices

* Checkpoint on switching over to global shared nodes

* Do the thing

* Revert unintended change to coalesced linestring behavior

* Take shared nodes into account in early simplification

* Let it do more sorting in memory

* Fix overnoding of collinear linestrings

* Remove duplicate nodes, since only duplicate vertices now matter

* Remember to delete temporary files

* Fix out of bounds memory access below, apparently

* Fix the actual undefined behavior

* Still running out of memory in one case. Find out where.

* Forgot the conditional

* Try again to make it not run out of memory

* Update version and changelog

CHANGELOG.md

@@ -1,3 +1,8 @@
+# 2.33.0
+
+* Further reduce memory usage of --no-simplification-of-shared-nodes by calculating the list of shared nodes globally using temporary files rather than in memory for each individual tile
+* Make --no-simplification-of-shared-nodes behave for LineStrings as it does for Polygons, preventing simplification only at crossings, convergences, and divergences, not at every point along collinear segments
+
 # 2.32.1
 
 * Reduce memory usage of --no-simplification-of-shared-nodes for polygons

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
+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
 	$(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread
 
 tippecanoe-enumerate: enumerate.o
@@ -73,7 +73,7 @@ tile-join: tile-join.o projection.o mbtiles.o mvt.o memfile.o dirtiles.o jsonpul
 tippecanoe-json-tool: jsontool.o jsonpull/jsonpull.o csv.o text.o geojson-loop.o
 	$(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread
 
-unit: unit.o text.o
+unit: unit.o text.o sort.o
 	$(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread
 
 tippecanoe-overzoom: overzoom.o mvt.o clip.o

README.md

@@ -492,7 +492,7 @@ the same layer, enclose them in an `all` expression so they will all be evaluate
  * `-ps` or `--no-line-simplification`: Don't simplify lines and polygons
  * `-pS` or `--simplify-only-low-zooms`: Don't simplify lines and polygons at maxzoom (but do simplify at lower zooms)
  * `--simplification-at-maximum-zoom=`_scale_: Use the specified _scale_ at maxzoom instead of the standard simplification scale (which still applies at lower zooms)
- * `-pn` or `--no-simplification-of-shared-nodes`: Don't simplify away nodes that appear in more than one feature or are used multiple times within the same feature, so that the intersection node will not be lost from intersecting roads. (This will not be effective if you also use `--coalesce`.) For polygons, the vertices where polygon rings meet will be detected, and the shared edges between the polygons will be simplified identically if possible. Use this instead of `--detect-shared-borders`.
+ * `-pn` or `--no-simplification-of-shared-nodes`: Don't simplify away nodes at which LineStrings or Polygon rings converge, diverge, or cross. (This will not be effective if you also use `--coalesce`.) In between intersection nodes, LineString segments or polygon edges will be simplified identically in each feature if possible. Use this instead of `--detect-shared-borders`.
  * `-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.

geometry.cpp

@@ -485,7 +485,7 @@ 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, size_t retain, drawvec const &shared_nodes) {
+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);
 
@@ -501,10 +501,35 @@ drawvec simplify_lines(drawvec &geom, int z, int detail, bool mark_tile_bounds,
 		}
 
 		if (prevent[P_SIMPLIFY_SHARED_NODES]) {
+			// This is kind of weird, because we have two lists of shared nodes to look through:
+			// * the drawvec, which is nodes that were introduced during clipping to the tile edge,
+			//   and which are in local tile coordinates
+			// * the shared_nodes_map, which was made globally before tiling began, and which
+			//   is in global quadkey coordinates.
+			// To look through the latter, we need to offset and encode the coordinates
+			// of the feature we are simplifying.
+
 			auto pt = std::lower_bound(shared_nodes.begin(), shared_nodes.end(), geom[i]);
 			if (pt != shared_nodes.end() && *pt == geom[i]) {
 				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;
+				}
+			}
 		}
 	}
 
@@ -661,6 +686,9 @@ drawvec fix_polygon(drawvec &geom) {
 				ring = tmp;
 			}
 
+			// Now we are rotating the ring to make the first/last point
+			// one that would be unlikely to be simplified away.
+
 			// calculate centroid
 			// a + 1 < size() because point 0 is duplicated at the end
 			long long xtotal = 0;

geometry.hpp

@@ -48,6 +48,10 @@ struct draw {
 		}
 	}
 
+	bool operator>(draw const &s) const {
+		return s < *this;
+	}
+
 	bool operator==(draw const &s) const {
 		return y == s.y && x == s.x;
 	}
@@ -73,7 +77,7 @@ drawvec clip_lines(drawvec &geom, int z, long long buffer);
 drawvec stairstep(drawvec &geom, int z, int detail);
 bool point_within_tile(long long x, long long y, int z);
 int quick_check(long long *bbox, int z, long long buffer);
-drawvec simplify_lines(drawvec &geom, int z, int detail, bool mark_tile_bounds, double simplification, size_t retain, drawvec const &shared_nodes);
+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);
 drawvec reorder_lines(drawvec &geom);
 drawvec fix_polygon(drawvec &geom);
 std::vector<drawvec> chop_polygon(std::vector<drawvec> &geoms);

main.cpp

@@ -65,6 +65,7 @@
 #include "text.hpp"
 #include "errors.hpp"
 #include "read_json.hpp"
+#include "sort.hpp"
 
 static int low_detail = 12;
 static int full_detail = -1;
@@ -121,7 +122,8 @@ void checkdisk(std::vector<struct reader> *r) {
 	for (size_t i = 0; i < r->size(); i++) {
 		// Pool and tree are used once.
 		// Geometry and index will be duplicated during sorting and tiling.
-		used += 2 * (*r)[i].geompos + 2 * (*r)[i].indexpos + (*r)[i].poolfile->off + (*r)[i].treefile->off;
+		used += 2 * (*r)[i].geompos + 2 * (*r)[i].indexpos + (*r)[i].poolfile->off + (*r)[i].treefile->off +
+			(*r)[i].vertexpos + (*r)[i].nodepos;
 	}
 
 	static int warned = 0;
@@ -1207,6 +1209,34 @@ void choose_first_zoom(long long *file_bbox, long long *file_bbox1, long long *f
 	}
 }
 
+int vertexcmp(const void *void1, const void *void2) {
+	vertex *v1 = (vertex *) void1;
+	vertex *v2 = (vertex *) void2;
+
+	if (v1->mid < v2->mid) {
+		return -1;
+	}
+	if (v1->mid > v2->mid) {
+		return 1;
+	}
+
+	if (v1->p1 < v2->p1) {
+		return -1;
+	}
+	if (v1->p1 > v2->p1) {
+		return 1;
+	}
+
+	if (v1->p2 < v2->p2) {
+		return -1;
+	}
+	if (v1->p2 > v2->p2) {
+		return 1;
+	}
+
+	return 0;
+}
+
 std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzoom, int basezoom, double basezoom_marker_width, sqlite3 *outdb, const char *outdir, std::set<std::string> *exclude, std::set<std::string> *include, int exclude_all, json_object *filter, double droprate, int buffer, const char *tmpdir, double gamma, int read_parallel, int forcetable, const char *attribution, bool uses_gamma, long long *file_bbox, long long *file_bbox1, long long *file_bbox2, const char *prefilter, const char *postfilter, const char *description, bool guess_maxzoom, bool guess_cluster_maxzoom, std::map<std::string, int> const *attribute_types, const char *pgm, std::map<std::string, attribute_op> const *attribute_accum, std::map<std::string, std::string> const &attribute_descriptions, std::string const &commandline, int minimum_maxzoom) {
 	int ret = EXIT_SUCCESS;
 
@@ -1219,11 +1249,15 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 		char treename[strlen(tmpdir) + strlen("/tree.XXXXXXXX") + 1];
 		char geomname[strlen(tmpdir) + strlen("/geom.XXXXXXXX") + 1];
 		char indexname[strlen(tmpdir) + strlen("/index.XXXXXXXX") + 1];
+		char vertexname[strlen(tmpdir) + strlen("/vertex.XXXXXXXX") + 1];
+		char nodename[strlen(tmpdir) + strlen("/node.XXXXXXXX") + 1];
 
 		snprintf(poolname, sizeof(poolname), "%s%s", tmpdir, "/pool.XXXXXXXX");
 		snprintf(treename, sizeof(treename), "%s%s", tmpdir, "/tree.XXXXXXXX");
 		snprintf(geomname, sizeof(geomname), "%s%s", tmpdir, "/geom.XXXXXXXX");
 		snprintf(indexname, sizeof(indexname), "%s%s", tmpdir, "/index.XXXXXXXX");
+		snprintf(vertexname, sizeof(vertexname), "%s%s", tmpdir, "/vertex.XXXXXXXX");
+		snprintf(nodename, sizeof(nodename), "%s%s", tmpdir, "/node.XXXXXXXX");
 
 		r->poolfd = mkstemp_cloexec(poolname);
 		if (r->poolfd < 0) {
@@ -1245,6 +1279,16 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 			perror(indexname);
 			exit(EXIT_OPEN);
 		}
+		r->vertexfd = mkstemp_cloexec(vertexname);
+		if (r->vertexfd < 0) {
+			perror(vertexname);
+			exit(EXIT_OPEN);
+		}
+		r->nodefd = mkstemp_cloexec(nodename);
+		if (r->nodefd < 0) {
+			perror(nodename);
+			exit(EXIT_OPEN);
+		}
 
 		r->poolfile = memfile_open(r->poolfd);
 		if (r->poolfile == NULL) {
@@ -1266,13 +1310,27 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 			perror(indexname);
 			exit(EXIT_OPEN);
 		}
+		r->vertexfile = fopen_oflag(vertexname, "w+b", O_RDWR | O_CLOEXEC);
+		if (r->vertexfile == NULL) {
+			perror(("open vertexfile " + std::string(vertexname)).c_str());
+			exit(EXIT_OPEN);
+		}
+		r->nodefile = fopen_oflag(nodename, "w+b", O_RDWR | O_CLOEXEC);
+		if (r->nodefile == NULL) {
+			perror(nodename);
+			exit(EXIT_OPEN);
+		}
 		r->geompos = 0;
 		r->indexpos = 0;
+		r->vertexpos = 0;
+		r->nodepos = 0;
 
 		unlink(poolname);
 		unlink(treename);
 		unlink(geomname);
 		unlink(indexname);
+		unlink(vertexname);
+		unlink(nodename);
 
 		// To distinguish a null value
 		{
@@ -1842,6 +1900,9 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 		fflush(stderr);
 	}
 
+	std::atomic<long long> vertexpos(0);
+	std::atomic<long long> nodepos(0);
+
 	for (size_t i = 0; i < CPUS; i++) {
 		if (fclose(readers[i].geomfile) != 0) {
 			perror("fclose geom");
@@ -1857,6 +1918,13 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 			perror("stat geom\n");
 			exit(EXIT_STAT);
 		}
+
+		vertexpos += readers[i].vertexpos;
+		nodepos += readers[i].nodepos;
+	}
+
+	if (!quiet) {
+		fprintf(stderr, "Merging string pool           \r");
 	}
 
 	// Create a combined string pool
@@ -1885,7 +1953,6 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 	}
 
 	unlink(poolname);
-
 	std::atomic<long long> poolpos(0);
 
 	for (size_t i = 0; i < CPUS; i++) {
@@ -1949,6 +2016,162 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 		madvise(stringpool, poolpos, MADV_RANDOM);
 	}
 
+	if (!quiet) {
+		fprintf(stderr, "Merging vertices              \r");
+	}
+
+	// Sort the vertices;
+	// find nodes where the same central point is part of two different vertices
+	{
+		std::string tmpname = std::string(tmpdir) + "/vertex2.XXXXXX";
+		int vertexfd = mkstemp((char *) tmpname.c_str());
+		if (vertexfd < 0) {
+			perror(("mkstemp vertexfile " + std::string(tmpname)).c_str());
+			exit(EXIT_OPEN);
+		}
+		unlink(tmpname.c_str());
+		FILE *vertex_out = fdopen(vertexfd, "w+b");
+		if (vertex_out == NULL) {
+			perror(tmpname.c_str());
+			exit(EXIT_OPEN);
+		}
+
+		std::vector<FILE *> vertex_readers;
+		for (size_t i = 0; i < CPUS; i++) {
+			vertex_readers.push_back(readers[i].vertexfile);
+			rewind(readers[i].vertexfile);
+		}
+		fqsort(vertex_readers, sizeof(vertex), vertexcmp, vertex_out, memsize / 10);
+
+		for (size_t i = 0; i < CPUS; i++) {
+			if (fclose(readers[i].vertexfile) != 0) {
+				perror("fclose vertex");
+				exit(EXIT_CLOSE);
+			}
+		}
+
+		rewind(vertex_out);
+
+		vertex prev(draw(VT_MOVETO, 0, 0), draw(VT_MOVETO, 0, 0), draw(VT_MOVETO, 0, 0));
+		vertex v(draw(VT_MOVETO, 0, 0), draw(VT_MOVETO, 0, 0), draw(VT_MOVETO, 0, 0));
+		while (fread((void *) &v, sizeof(vertex), 1, vertex_out)) {
+			if (v.mid == prev.mid && (v.p1 != prev.p1 || v.p2 != prev.p2)) {
+				long long x = v.mid.x * (1LL << geometry_scale);
+				long long y = v.mid.y * (1LL << geometry_scale);
+
+#if 0
+				double lon, lat;
+				tile2lonlat(x, y, 32, &lon, &lat);
+				printf("{\"type\":\"Feature\", \"properties\":{}, \"geometry\":{\"type\":\"Point\", \"coordinates\":[%f,%f]}}\n", lon, lat);
+#endif
+
+				struct node n;
+				n.index = encode_quadkey((unsigned) x, (unsigned) y);
+
+				fwrite_check((char *) &n, sizeof(struct node), 1, readers[0].nodefile, &readers[0].nodepos, "vertices");
+			}
+			prev = v;
+		}
+
+		fclose(vertex_out);
+	}
+
+	if (!quiet) {
+		fprintf(stderr, "Merging nodes                 \r");
+	}
+
+	// Sort nodes that can't be simplified away; scan the list to remove duplicates
+
+	FILE *shared_nodes;
+	node *shared_nodes_map = NULL;	// will be null if there are no shared nodes
+	{
+		// sort
+
+		std::string tmpname = std::string(tmpdir) + "/node2.XXXXXX";
+		int nodefd = mkstemp((char *) tmpname.c_str());
+		if (nodefd < 0) {
+			perror(("mkstemp nodefile " + std::string(tmpname)).c_str());
+			exit(EXIT_OPEN);
+		}
+		unlink(tmpname.c_str());
+		FILE *node_out;
+		node_out = fdopen(nodefd, "w+b");
+		if (node_out == NULL) {
+			perror(tmpname.c_str());
+			exit(EXIT_OPEN);
+		}
+
+		std::vector<FILE *> node_readers;
+		for (size_t i = 0; i < CPUS; i++) {
+			node_readers.push_back(readers[i].nodefile);
+			rewind(readers[i].nodefile);
+		}
+
+		fqsort(node_readers, sizeof(node), nodecmp, node_out, memsize / 10);
+
+		for (size_t i = 0; i < CPUS; i++) {
+			if (fclose(readers[i].nodefile) != 0) {
+				perror("fclose node");
+				exit(EXIT_CLOSE);
+			}
+		}
+
+		rewind(node_out);
+
+		// scan
+
+		tmpname = std::string(tmpdir) + "/node3.XXXXXX";
+		nodefd = mkstemp((char *) tmpname.c_str());
+		if (nodefd < 0) {
+			perror(("mkstemp nodefile " + std::string(tmpname)).c_str());
+			exit(EXIT_OPEN);
+		}
+		unlink(tmpname.c_str());
+		shared_nodes = fdopen(nodefd, "w+b");
+		if (shared_nodes == NULL) {
+			perror(tmpname.c_str());
+			exit(EXIT_OPEN);
+		}
+
+		// `written` is to see if this node has already been preserved
+		// and doesn't need to be preserved again
+		struct node written;
+		written.index = ULONG_MAX;
+
+		nodepos = 0;
+		struct node here;
+		while (fread((void *) &here, sizeof(here), 1, node_out)) {
+			if (nodecmp((void *) &here, (void *) &written) != 0) {
+				fwrite_check((void *) &here, sizeof(here), 1, shared_nodes, &nodepos, "shared nodes");
+				written = here;
+
+#if 0
+				unsigned wx, wy;
+				decode_quadkey(here.index, &wx, &wy);
+				double lon, lat;
+				tile2lonlat(wx, wy, 32, &lon, &lat);
+				printf("{\"type\":\"Feature\", \"properties\":{}, \"geometry\":{\"type\":\"Point\", \"coordinates\":[%f,%f]}}\n", lon, lat);
+#endif
+			}
+		}
+
+		fflush(shared_nodes);
+
+		if (nodepos > 0) {
+			shared_nodes_map = (node *) mmap(NULL, nodepos, PROT_READ, MAP_PRIVATE, nodefd, 0);
+			if (shared_nodes_map == (node *) MAP_FAILED) {
+				perror("mmap nodes");
+				exit(EXIT_MEMORY);
+			}
+		}
+
+		fclose(node_out);
+	}
+
+	if (!quiet) {
+		fprintf(stderr, "Merging index                 \r");
+	}
+
 	char indexname[strlen(tmpdir) + strlen("/index.XXXXXXXX") + 1];
 	snprintf(indexname, sizeof(indexname), "%s%s", tmpdir, "/index.XXXXXXXX");
 
@@ -2023,7 +2246,9 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 		long long s = progress_seq;
 		long long geompos_print = geompos;
 		long long poolpos_print = poolpos;
-		fprintf(stderr, "%lld features, %lld bytes of geometry, %lld bytes of string pool\n", s, geompos_print, poolpos_print);
+		long long vertexpos_print = vertexpos;
+		long long nodepos_print = nodepos;
+		fprintf(stderr, "%lld features, %lld bytes of geometry and attributes, %lld bytes of string pool, %lld bytes of vertices, %lld bytes of nodes\n", s, geompos_print, poolpos_print, vertexpos_print, nodepos_print);
 	}
 
 	if (indexpos == 0) {
@@ -2532,7 +2757,7 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 	std::atomic<unsigned> midx(0);
 	std::atomic<unsigned> midy(0);
 	std::vector<strategy> strategies;
-	int written = traverse_zooms(fd, size, stringpool, &midx, &midy, maxzoom, minzoom, outdb, outdir, buffer, fname, tmpdir, gamma, full_detail, low_detail, min_detail, pool_off, initial_x, initial_y, simplification, maxzoom_simplification, layermaps, prefilter, postfilter, attribute_accum, filter, strategies, iz);
+	int written = traverse_zooms(fd, size, stringpool, &midx, &midy, maxzoom, minzoom, outdb, outdir, buffer, fname, tmpdir, gamma, full_detail, low_detail, min_detail, pool_off, initial_x, initial_y, simplification, maxzoom_simplification, layermaps, prefilter, postfilter, attribute_accum, filter, strategies, iz, shared_nodes_map, nodepos);
 
 	if (maxzoom != written) {
 		if (written > minzoom) {
@@ -2555,6 +2780,8 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 		perror("close pool");
 	}
 
+	fclose(shared_nodes);
+
 	// mbtiles-style bounding box and center
 	double minlat = 0, minlon = 0, maxlat = 0, maxlon = 0, midlat = 0, midlon = 0;
 

man/tippecanoe.1

@@ -54,7 +54,7 @@ $ tippecanoe \-o file.mbtiles [options] [file.json file.json.gz file.fgb ...]
 .RE
 .PP
 If no files are specified, it reads GeoJSON from the standard input.
-If multiple files are specified, each is placed in its own layer.
+If multiple files are specified, each is placed in its own layer \[la]#input-files-and-layer-names\[ra]\&.
 .PP
 The GeoJSON features need not be wrapped in a FeatureCollection.
 You can concatenate multiple GeoJSON features or files together,
@@ -629,7 +629,7 @@ the line or polygon within one tile unit of its proper location. You can probabl
 .IP \(bu 2
 \fB\fC\-\-simplification\-at\-maximum\-zoom=\fR\fIscale\fP: Use the specified \fIscale\fP at maxzoom instead of the standard simplification scale (which still applies at lower zooms)
 .IP \(bu 2
-\fB\fC\-pn\fR or \fB\fC\-\-no\-simplification\-of\-shared\-nodes\fR: Don't simplify away nodes that appear in more than one feature or are used multiple times within the same feature, so that the intersection node will not be lost from intersecting roads. (This will not be effective if you also use \fB\fC\-\-coalesce\fR\&.) For polygons, the vertices where polygon rings meet will be detected, and the shared edges between the polygons will be simplified identically if possible. Use this instead of \fB\fC\-\-detect\-shared\-borders\fR\&.
+\fB\fC\-pn\fR or \fB\fC\-\-no\-simplification\-of\-shared\-nodes\fR: Don't simplify away nodes at which LineStrings or Polygon rings converge, diverge, or cross. (This will not be effective if you also use \fB\fC\-\-coalesce\fR\&.) In between intersection nodes, LineString segments or polygon edges will be simplified identically in each feature if possible. Use this instead of \fB\fC\-\-detect\-shared\-borders\fR\&.
 .IP \(bu 2
 \fB\fC\-pt\fR or \fB\fC\-\-no\-tiny\-polygon\-reduction\fR: Don't combine the area of very small polygons into small squares that represent their combined area.
 .IP \(bu 2

serial.cpp

@@ -393,6 +393,29 @@ static std::string strip_zeroes(std::string s) {
 	return s;
 }
 
+int nodecmp(const void *void1, const void *void2) {
+	node *n1 = (node *) void1;
+	node *n2 = (node *) void2;
+
+	if (n1->index < n2->index) {
+		return -1;
+	} else if (n1->index > n2->index) {
+		return 1;
+	}
+
+	return 0;
+}
+
+static void add_scaled_node(struct reader *r, serialization_state *sst, draw g) {
+	long long x = SHIFT_LEFT(g.x);
+	long long y = SHIFT_LEFT(g.y);
+
+	struct node n;
+	n.index = encode_quadkey((unsigned) x, (unsigned) y);
+
+	fwrite_check((char *) &n, sizeof(struct node), 1, r->nodefile, &r->nodepos, sst->fname);
+}
+
 // called from frontends
 int serialize_feature(struct serialization_state *sst, serial_feature &sf) {
 	struct reader *r = &(*sst->readers)[sst->segment];
@@ -446,11 +469,12 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf) {
 			scaled_geometry = simple_clip_poly(scaled_geometry, SHIFT_RIGHT(c.minx), SHIFT_RIGHT(c.miny), SHIFT_RIGHT(c.maxx), SHIFT_RIGHT(c.maxy), prevent[P_SIMPLIFY_SHARED_NODES]);
 		} else if (sf.t == VT_LINE) {
 			scaled_geometry = clip_lines(scaled_geometry, SHIFT_RIGHT(c.minx), SHIFT_RIGHT(c.miny), SHIFT_RIGHT(c.maxx), SHIFT_RIGHT(c.maxy));
-			scaled_geometry = remove_noop(scaled_geometry, sf.t, 0);
 		} else if (sf.t == VT_POINT) {
 			scaled_geometry = clip_point(scaled_geometry, SHIFT_RIGHT(c.minx), SHIFT_RIGHT(c.miny), SHIFT_RIGHT(c.maxx), SHIFT_RIGHT(c.maxy));
 		}
 
+		scaled_geometry = remove_noop(scaled_geometry, sf.t, 0);
+
 		sf.bbox[0] = LLONG_MAX;
 		sf.bbox[1] = LLONG_MAX;
 		sf.bbox[2] = LLONG_MIN;
@@ -480,6 +504,100 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf) {
 		return 1;
 	}
 
+	if (prevent[P_SIMPLIFY_SHARED_NODES]) {
+		scaled_geometry = remove_noop(scaled_geometry, sf.t, 0);
+
+		if (sf.t == VT_POLYGON || sf.t == VT_LINE) {
+			for (size_t i = 0; i < scaled_geometry.size(); i++) {
+				if (scaled_geometry[i].op == VT_MOVETO) {
+					size_t j;
+
+					for (j = i + 1; j < scaled_geometry.size(); j++) {
+						if (scaled_geometry[j].op != VT_LINETO) {
+							break;
+						}
+					}
+
+					if (sf.t == VT_POLYGON && j - i >= 4) {
+						for (size_t k = i; k < j - 1; k++) {
+							// % (j - i - 1) because we don't want the duplicate last point
+
+							struct vertex v(
+								scaled_geometry[(k - i + 0) % (j - i - 1) + i],
+								scaled_geometry[(k - i + 1) % (j - i - 1) + i],
+								scaled_geometry[(k - i + 2) % (j - i - 1) + i]);
+
+							fwrite_check((char *) &v, sizeof(struct vertex), 1, r->vertexfile, &r->vertexpos, sst->fname);
+						}
+					} else if (sf.t == VT_LINE && j - i >= 2) {
+						for (size_t k = i; k + 2 < j; k++) {
+							struct vertex v(
+								scaled_geometry[k + 0],
+								scaled_geometry[k + 1],
+								scaled_geometry[k + 2]);
+
+							fwrite_check((char *) &v, sizeof(struct vertex), 1, r->vertexfile, &r->vertexpos, sst->fname);
+						}
+					}
+
+					// since the starting point is never simplified away,
+					// don't let it be simplified away in any other polygons either.
+					// Needs to appear twice here so that the check below will see
+					// it as appearing in multiple features.
+					add_scaled_node(r, sst, scaled_geometry[i]);
+
+					if (sf.t == VT_LINE && j - i >= 2) {
+						// linestrings also need to preserve the last point
+
+						add_scaled_node(r, sst, scaled_geometry[j - 1]);
+					} else if (sf.t == VT_POLYGON && j - i >= 4) {
+						// To avoid letting polygons get simplified away to nothing,
+						// also keep the furthest-away point from the initial point
+						// (which Douglas-Peucker simplification would keep anyway,
+						// if its search weren't being split up by polygon side).
+
+						double far = 0;
+						size_t which = i;
+						for (size_t k = i + 1; k < j - 1; k++) {
+							double xd = scaled_geometry[k].x - scaled_geometry[i].x;
+							double yd = scaled_geometry[k].y - scaled_geometry[i].y;
+							double d = xd * xd + yd * yd;
+							if (d > far ||
+							    ((d == far) && (scaled_geometry[k] < scaled_geometry[which]))) {
+								far = d;
+								which = k;
+							}
+						}
+
+						add_scaled_node(r, sst, scaled_geometry[which]);
+
+						// And, likewise, the point most distant from those two points,
+						// which probably would also be the one that Douglas-Peucker
+						// would keep next.
+
+						far = 0;
+						size_t which2 = i;
+
+						for (size_t k = i + 1; k < j - 1; k++) {
+							double d = distance_from_line(scaled_geometry[k].x, scaled_geometry[k].y,
+										      scaled_geometry[i].x, scaled_geometry[i].y,
+										      scaled_geometry[which].x, scaled_geometry[which].y);
+							if ((d > far) ||
+							    ((d == far) && (scaled_geometry[k] < scaled_geometry[which2]))) {
+								far = d;
+								which2 = k;
+							}
+						}
+
+						add_scaled_node(r, sst, scaled_geometry[which2]);
+					}
+
+					i = j - 1;
+				}
+			}
+		}
+	}
+
 	if (!sf.has_id) {
 		if (additional[A_GENERATE_IDS]) {
 			sf.has_id = true;
@@ -583,7 +701,9 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf) {
 		for (size_t i = 0; i < scaled_geometry.size(); i++) {
 			ix += scaled_geometry[i].x + scaled_geometry[i].y;
 		}
-		ix = ix % scaled_geometry.size();
+		if (scaled_geometry.size() != 0) {
+			ix = ix % scaled_geometry.size();
+		}
 
 		// If off the edge of the plane, mask to bring it back into the addressable area
 		midx = SHIFT_LEFT(scaled_geometry[ix].x) & ((1LL << 32) - 1);

serial.hpp

@@ -78,14 +78,20 @@ struct reader {
 	int treefd = -1;
 	int geomfd = -1;
 	int indexfd = -1;
+	int vertexfd = -1;
+	int nodefd = -1;
 
 	struct memfile *poolfile = NULL;
 	struct memfile *treefile = NULL;
 	FILE *geomfile = NULL;
 	FILE *indexfile = NULL;
+	FILE *vertexfile = NULL;
+	FILE *nodefile = NULL;
 
 	std::atomic<long long> geompos;
 	std::atomic<long long> indexpos;
+	std::atomic<long long> vertexpos;
+	std::atomic<long long> nodepos;
 
 	long long file_bbox[4] = {0, 0, 0, 0};
 
@@ -97,7 +103,7 @@ struct reader {
 	char *geom_map = NULL;
 
 	reader()
-	    : geompos(0), indexpos(0) {
+	    : geompos(0), indexpos(0), vertexpos(0), nodepos(0) {
 	}
 
 	reader(reader const &r) {
@@ -105,11 +111,15 @@ struct reader {
 		treefd = r.treefd;
 		geomfd = r.geomfd;
 		indexfd = r.indexfd;
+		vertexfd = r.vertexfd;
+		nodefd = r.nodefd;
 
 		poolfile = r.poolfile;
 		treefile = r.treefile;
 		geomfile = r.geomfile;
 		indexfile = r.indexfile;
+		vertexfile = r.vertexfile;
+		nodefile = r.nodefile;
 
 		long long p = r.geompos;
 		geompos = p;
@@ -117,6 +127,12 @@ struct reader {
 		p = r.indexpos;
 		indexpos = p;
 
+		p = r.vertexpos;
+		vertexpos = p;
+
+		p = r.nodepos;
+		nodepos = p;
+
 		memcpy(file_bbox, r.file_bbox, sizeof(file_bbox));
 
 		geomst = r.geomst;
@@ -158,6 +174,52 @@ struct serialization_state {
 	int exclude_all = 0;
 };
 
+struct vertex {
+	// these are scaled geometry,
+	// but because scaling is disabled if P_SHARED_NODES is set,
+	// they are effectively also world coordinates
+	draw p1;
+	draw mid;
+	draw p2;
+
+	vertex(draw one, draw joint, draw two) {
+		if (one < two) {
+			p1 = one;
+			p2 = two;
+		} else {
+			p1 = two;
+			p2 = one;
+		}
+
+		mid = joint;
+	}
+
+	bool operator<(const vertex &v) const {
+		if (mid < v.mid) {
+			return true;
+		} else if (mid == v.mid) {
+			if (p1 < v.p1) {
+				return true;
+			} else if (p1 == v.p1) {
+				if (p2 < v.p2) {
+					return true;
+				}
+			}
+		}
+
+		return false;
+	}
+};
+
+struct node {
+	// this is in quadkey coordinates so that the nodes for each tile
+	// will be adjacent in memory, reducing potential thrashing during
+	// the binary search.
+	unsigned long long index;
+};
+
+int nodecmp(const void *void1, const void *void2);
+
 int serialize_feature(struct serialization_state *sst, serial_feature &sf);
 void coerce_value(std::string const &key, int &vt, std::string &val, std::map<std::string, int> const *attribute_types);
 

sort.cpp

@@ -0,0 +1,123 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <vector>
+#include <string>
+
+#define MAX_MEMORY (10 * 1024 * 1024)
+
+void fqsort(std::vector<FILE *> &inputs, size_t width, int (*cmp)(const void *, const void *), FILE *out, size_t mem) {
+	std::string pivot;
+	FILE *fp1, *fp2;
+
+	{
+		// read some elements into memory to choose a pivot from
+		//
+		// this is in its own scope so `buf` can go out of scope
+		// before trying to do any sub-sorts.
+
+		std::string buf;
+
+		bool read_everything = false;
+		for (size_t i = 0; i < inputs.size(); i++) {
+			if (buf.size() > mem) {
+				break;
+			}
+
+			while (true) {
+				std::string element;
+				element.resize(width);
+
+				size_t n = fread((void *) element.c_str(), width, 1, inputs[i]);
+				if (n == 0) {
+					if (i + 1 == inputs.size()) {
+						read_everything = true;
+					}
+					break;
+				}
+
+				buf.append(element);
+
+				if (buf.size() > mem) {
+					break;
+				}
+			}
+		}
+
+		qsort((void *) buf.c_str(), buf.size() / width, width, cmp);
+
+		// If that was everything we have to sort, we are done.
+
+		if (read_everything) {
+			fwrite((void *) buf.c_str(), buf.size() / width, width, out);
+			return;
+		}
+
+		// Otherwise, choose a pivot from it, make some temporary files,
+		// write what we have to those files, and then partition the rest
+		// of the input into them.
+
+		// This would be unstable if the pivot is one of several elements
+		// that compare equal. Does it matter?
+
+		size_t pivot_off = width * (buf.size() / width / 2);
+		pivot = std::string(buf, pivot_off, width);
+
+		std::string t1 = "/tmp/sort1.XXXXXX";
+		std::string t2 = "/tmp/sort2.XXXXXX";
+
+		int fd1 = mkstemp((char *) t1.c_str());
+		unlink(t1.c_str());
+		int fd2 = mkstemp((char *) t2.c_str());
+		unlink(t2.c_str());
+
+		fp1 = fdopen(fd1, "w+b");
+		if (fp1 == NULL) {
+			perror(t1.c_str());
+			exit(EXIT_FAILURE);
+		}
+		fp2 = fdopen(fd2, "w+b");
+		if (fp2 == NULL) {
+			perror(t2.c_str());
+			exit(EXIT_FAILURE);
+		}
+
+		fwrite((void *) buf.c_str(), sizeof(char), pivot_off, fp1);
+		fwrite((void *) ((char *) buf.c_str() + pivot_off), sizeof(char), buf.size() - pivot_off, fp2);
+	}
+
+	// read the remaining input into the temporary files
+
+	for (size_t i = 0; i < inputs.size(); i++) {
+		while (true) {
+			std::string element;
+			element.resize(width);
+
+			size_t n = fread((void *) element.c_str(), width, 1, inputs[i]);
+			if (n == 0) {
+				break;
+			}
+
+			if (cmp((void *) element.c_str(), (void *) pivot.c_str()) < 0) {
+				fwrite((void *) element.c_str(), width, 1, fp1);
+			} else {
+				fwrite((void *) element.c_str(), width, 1, fp2);
+			}
+		}
+	}
+
+	// Now sort the sub-ranges into the output.
+
+	rewind(fp1);
+	rewind(fp2);
+
+	std::vector<FILE *> v1;
+	v1.emplace_back(fp1);
+	fqsort(v1, width, cmp, out, mem);
+	fclose(fp1);
+
+	std::vector<FILE *> v2;
+	v2.emplace_back(fp2);
+	fqsort(v2, width, cmp, out, mem);
+	fclose(fp2);
+}

sort.hpp

@@ -0,0 +1,6 @@
+#ifndef SORT_HPP
+#define SORT_HPP
+
+void fqsort(std::vector<FILE *> &inputs, size_t width, int (*cmp)(const void *, const void *), FILE *out, size_t mem);
+
+#endif

tests/tl_2018_51685_roads/out/-Z11_-z11_--no-simplification-of-shared-nodes.json

@@ -130,21 +130,21 @@
 ,
 { "type": "Feature", "properties": { "LINEARID": "1103679538189", "FULLNAME": "Old Centreville Rd", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.458506, 38.787409 ], [ -77.458720, 38.785602 ], [ -77.458935, 38.783762 ], [ -77.459021, 38.783227 ], [ -77.459064, 38.782290 ], [ -77.459106, 38.781856 ], [ -77.459493, 38.778410 ], [ -77.459450, 38.778209 ], [ -77.454772, 38.775532 ], [ -77.453613, 38.774897 ], [ -77.452970, 38.774462 ], [ -77.452669, 38.774094 ], [ -77.452497, 38.773525 ], [ -77.452283, 38.773458 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "110410477217", "FULLNAME": "State Rte 28", "RTTYP": "S", "MTFCC": "S1200" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452111, 38.773425 ], [ -77.452283, 38.772622 ], [ -77.452326, 38.772387 ], [ -77.452455, 38.771886 ], [ -77.452497, 38.771685 ], [ -77.452583, 38.771317 ], [ -77.452712, 38.770681 ], [ -77.452798, 38.770514 ], [ -77.452884, 38.770246 ], [ -77.452970, 38.770012 ], [ -77.453098, 38.769811 ], [ -77.453270, 38.769510 ], [ -77.453356, 38.769343 ], [ -77.453570, 38.769042 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "110410477217", "FULLNAME": "State Rte 28", "RTTYP": "S", "MTFCC": "S1200" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452111, 38.773425 ], [ -77.452497, 38.771685 ], [ -77.452884, 38.770246 ], [ -77.453270, 38.769510 ], [ -77.453570, 38.769042 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "110410477454", "FULLNAME": "Centreville Rd", "RTTYP": "M", "MTFCC": "S1200" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452111, 38.773425 ], [ -77.452283, 38.772622 ], [ -77.452326, 38.772387 ], [ -77.452455, 38.771886 ], [ -77.452497, 38.771685 ], [ -77.452583, 38.771317 ], [ -77.452712, 38.770681 ], [ -77.452798, 38.770514 ], [ -77.452884, 38.770246 ], [ -77.452970, 38.770012 ], [ -77.453098, 38.769811 ], [ -77.453270, 38.769510 ], [ -77.453356, 38.769343 ], [ -77.453570, 38.769042 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "110410477454", "FULLNAME": "Centreville Rd", "RTTYP": "M", "MTFCC": "S1200" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452111, 38.773425 ], [ -77.452497, 38.771685 ], [ -77.452884, 38.770246 ], [ -77.453270, 38.769510 ], [ -77.453570, 38.769042 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1106073064085", "FULLNAME": "Blooms Quarry Ln", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452111, 38.773425 ], [ -77.451854, 38.773391 ], [ -77.451768, 38.773391 ], [ -77.451725, 38.773391 ], [ -77.451639, 38.773391 ], [ -77.451339, 38.773458 ], [ -77.451210, 38.773492 ], [ -77.450867, 38.773625 ], [ -77.450824, 38.773659 ], [ -77.450695, 38.773659 ], [ -77.450566, 38.773692 ], [ -77.450266, 38.773692 ], [ -77.449965, 38.773659 ], [ -77.449880, 38.773625 ], [ -77.449322, 38.773358 ], [ -77.448163, 38.772856 ], [ -77.447605, 38.772588 ], [ -77.447305, 38.772387 ], [ -77.447047, 38.772187 ], [ -77.446790, 38.771986 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1106073064085", "FULLNAME": "Blooms Quarry Ln", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452111, 38.773425 ], [ -77.451339, 38.773458 ], [ -77.450824, 38.773659 ], [ -77.450266, 38.773692 ], [ -77.449965, 38.773659 ], [ -77.449322, 38.773358 ], [ -77.447605, 38.772588 ], [ -77.446790, 38.771986 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1102214409896", "FULLNAME": "Centreville Rd", "RTTYP": "M", "MTFCC": "S1200" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452283, 38.773458 ], [ -77.452326, 38.773224 ], [ -77.452497, 38.772421 ], [ -77.452669, 38.771718 ], [ -77.452841, 38.770882 ], [ -77.452927, 38.770581 ], [ -77.453012, 38.770380 ], [ -77.453098, 38.770146 ], [ -77.453227, 38.769845 ], [ -77.453399, 38.769543 ], [ -77.453485, 38.769376 ], [ -77.453699, 38.769075 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1102214409896", "FULLNAME": "Centreville Rd", "RTTYP": "M", "MTFCC": "S1200" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452283, 38.773458 ], [ -77.452669, 38.771718 ], [ -77.452927, 38.770581 ], [ -77.453227, 38.769845 ], [ -77.453699, 38.769075 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1102214411210", "FULLNAME": "State Rte 28", "RTTYP": "S", "MTFCC": "S1200" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452283, 38.773458 ], [ -77.452326, 38.773224 ], [ -77.452497, 38.772421 ], [ -77.452669, 38.771718 ], [ -77.452841, 38.770882 ], [ -77.452927, 38.770581 ], [ -77.453012, 38.770380 ], [ -77.453098, 38.770146 ], [ -77.453227, 38.769845 ], [ -77.453399, 38.769543 ], [ -77.453485, 38.769376 ], [ -77.453699, 38.769075 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1102214411210", "FULLNAME": "State Rte 28", "RTTYP": "S", "MTFCC": "S1200" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452283, 38.773458 ], [ -77.452669, 38.771718 ], [ -77.452927, 38.770581 ], [ -77.453227, 38.769845 ], [ -77.453699, 38.769075 ] ] } }
 ,
 { "type": "Feature", "properties": { "LINEARID": "110410477418", "FULLNAME": "Manassas Dr", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.443743, 38.768640 ], [ -77.443914, 38.768841 ], [ -77.445374, 38.769711 ], [ -77.445545, 38.769778 ], [ -77.446747, 38.770581 ], [ -77.447648, 38.770982 ], [ -77.448807, 38.771250 ], [ -77.449880, 38.771417 ], [ -77.452283, 38.771752 ], [ -77.452497, 38.771685 ] ] } }
 ,
 { "type": "Feature", "properties": { "LINEARID": "110410477419", "FULLNAME": "Manassas Dr", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452497, 38.771685 ], [ -77.452326, 38.771584 ], [ -77.452025, 38.771518 ], [ -77.449880, 38.771250 ], [ -77.448764, 38.771083 ], [ -77.447777, 38.770815 ], [ -77.446790, 38.770380 ], [ -77.445674, 38.769644 ], [ -77.445545, 38.769543 ], [ -77.444601, 38.769042 ], [ -77.443743, 38.768640 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1106092774118", "FULLNAME": "Manassas Dr", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452283, 38.773458 ], [ -77.452111, 38.773425 ], [ -77.451854, 38.773391 ], [ -77.451768, 38.773391 ], [ -77.451725, 38.773391 ], [ -77.451639, 38.773391 ], [ -77.451339, 38.773458 ], [ -77.451210, 38.773492 ], [ -77.450867, 38.773625 ], [ -77.450824, 38.773659 ], [ -77.450695, 38.773659 ], [ -77.450566, 38.773692 ], [ -77.450266, 38.773692 ], [ -77.449965, 38.773659 ], [ -77.449880, 38.773625 ], [ -77.449322, 38.773358 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1106092774118", "FULLNAME": "Manassas Dr", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.452283, 38.773458 ], [ -77.452111, 38.773425 ], [ -77.451339, 38.773458 ], [ -77.450824, 38.773659 ], [ -77.450266, 38.773692 ], [ -77.449965, 38.773659 ], [ -77.449322, 38.773358 ] ] } }
 ,
 { "type": "Feature", "properties": { "LINEARID": "110410477551", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.450566, 38.772287 ], [ -77.450824, 38.772387 ], [ -77.450910, 38.772254 ], [ -77.450738, 38.772153 ], [ -77.450566, 38.772287 ] ] } }
 ,
@@ -154,7 +154,7 @@
 ,
 { "type": "Feature", "properties": { "LINEARID": "110410477440", "FULLNAME": "Sunnyside Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.449665, 38.772053 ], [ -77.449493, 38.771952 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1106087402596", "FULLNAME": "Blooms Quarry Rd", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.449322, 38.773358 ], [ -77.448163, 38.772856 ], [ -77.447605, 38.772588 ], [ -77.447305, 38.772387 ], [ -77.447047, 38.772187 ], [ -77.446790, 38.771986 ], [ -77.445030, 38.770614 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1106087402596", "FULLNAME": "Blooms Quarry Rd", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.449322, 38.773358 ], [ -77.447605, 38.772588 ], [ -77.446790, 38.771986 ], [ -77.445030, 38.770614 ] ] } }
 ,
 { "type": "Feature", "properties": { "LINEARID": "110410477452", "FULLNAME": "White Pine Dr", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.447348, 38.771417 ], [ -77.446961, 38.771250 ], [ -77.446876, 38.771350 ], [ -77.447176, 38.771584 ], [ -77.448249, 38.772321 ], [ -77.449064, 38.772689 ], [ -77.450180, 38.773157 ] ] } }
 ,
@@ -336,13 +336,13 @@
 ,
 { "type": "Feature", "properties": { "LINEARID": "110410477422", "FULLNAME": "Primrose Ln", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.434087, 38.763989 ], [ -77.433829, 38.764190 ], [ -77.433658, 38.764457 ], [ -77.433443, 38.764491 ], [ -77.433186, 38.764625 ], [ -77.433100, 38.764758 ], [ -77.433143, 38.765026 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1103678606081", "FULLNAME": "Silver Meteor Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.432799, 38.762650 ], [ -77.432885, 38.762851 ], [ -77.432971, 38.763119 ], [ -77.432971, 38.763186 ], [ -77.432971, 38.763253 ], [ -77.432842, 38.763420 ], [ -77.432628, 38.763822 ], [ -77.432585, 38.763855 ], [ -77.432499, 38.763888 ], [ -77.432284, 38.763955 ], [ -77.432241, 38.763955 ], [ -77.432156, 38.763955 ], [ -77.431855, 38.763822 ], [ -77.431812, 38.763822 ], [ -77.431812, 38.763788 ], [ -77.431812, 38.763721 ], [ -77.431941, 38.763353 ], [ -77.432241, 38.762952 ], [ -77.432284, 38.762918 ], [ -77.432456, 38.762784 ], [ -77.432542, 38.762717 ], [ -77.432628, 38.762684 ], [ -77.432714, 38.762684 ], [ -77.432799, 38.762650 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1103678606081", "FULLNAME": "Silver Meteor Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.432799, 38.762650 ], [ -77.432971, 38.763253 ], [ -77.432585, 38.763855 ], [ -77.432156, 38.763955 ], [ -77.431812, 38.763788 ], [ -77.431941, 38.763353 ], [ -77.432284, 38.762918 ], [ -77.432542, 38.762717 ], [ -77.432799, 38.762650 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1103679538423", "FULLNAME": "Black Hawk Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.432799, 38.762650 ], [ -77.432885, 38.762851 ], [ -77.432971, 38.763119 ], [ -77.432971, 38.763186 ], [ -77.432971, 38.763253 ], [ -77.432842, 38.763420 ], [ -77.432628, 38.763822 ], [ -77.432585, 38.763855 ], [ -77.432499, 38.763888 ], [ -77.432284, 38.763955 ], [ -77.432241, 38.763955 ], [ -77.432156, 38.763955 ], [ -77.431855, 38.763822 ], [ -77.431812, 38.763822 ], [ -77.431812, 38.763788 ], [ -77.431812, 38.763721 ], [ -77.431941, 38.763353 ], [ -77.432241, 38.762952 ], [ -77.432284, 38.762918 ], [ -77.432456, 38.762784 ], [ -77.432542, 38.762717 ], [ -77.432628, 38.762684 ], [ -77.432714, 38.762684 ], [ -77.432799, 38.762650 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1103679538423", "FULLNAME": "Black Hawk Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.432799, 38.762650 ], [ -77.432971, 38.763253 ], [ -77.432585, 38.763855 ], [ -77.432156, 38.763955 ], [ -77.431812, 38.763788 ], [ -77.431941, 38.763353 ], [ -77.432284, 38.762918 ], [ -77.432542, 38.762717 ], [ -77.432799, 38.762650 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1103678606039", "FULLNAME": "Silver Meteor Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.433400, 38.762583 ], [ -77.433100, 38.762650 ], [ -77.432928, 38.762684 ], [ -77.432799, 38.762650 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1103678606039", "FULLNAME": "Silver Meteor Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.433400, 38.762583 ], [ -77.432799, 38.762650 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1103679538422", "FULLNAME": "Black Hawk Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.433400, 38.762583 ], [ -77.433100, 38.762650 ], [ -77.432928, 38.762684 ], [ -77.432799, 38.762650 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1103679538422", "FULLNAME": "Black Hawk Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.433400, 38.762583 ], [ -77.432799, 38.762650 ] ] } }
 ,
 { "type": "Feature", "properties": { "LINEARID": "110410477423", "FULLNAME": "Railroad Dr", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.442241, 38.762617 ], [ -77.441897, 38.762650 ], [ -77.441554, 38.762851 ], [ -77.440696, 38.763688 ], [ -77.440567, 38.763989 ], [ -77.440567, 38.764089 ] ] } }
 ,
@@ -434,9 +434,9 @@
 ,
 { "type": "Feature", "properties": { "LINEARID": "110410477227", "FULLNAME": "Gary Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.429924, 38.767636 ], [ -77.430696, 38.767101 ], [ -77.431211, 38.766733 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "1106087387886", "FULLNAME": "Manassas Dr", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.427735, 38.765361 ], [ -77.427692, 38.765796 ], [ -77.427649, 38.767302 ], [ -77.427478, 38.768673 ], [ -77.427564, 38.769108 ], [ -77.427649, 38.769242 ], [ -77.427735, 38.769376 ], [ -77.427821, 38.769476 ], [ -77.427950, 38.769577 ], [ -77.428079, 38.769677 ], [ -77.428207, 38.769778 ], [ -77.428293, 38.769811 ], [ -77.428422, 38.769878 ], [ -77.428551, 38.769911 ], [ -77.428722, 38.769945 ], [ -77.428808, 38.769978 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "1106087387886", "FULLNAME": "Manassas Dr", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.427735, 38.765361 ], [ -77.427692, 38.765796 ], [ -77.427649, 38.767302 ], [ -77.427478, 38.768673 ], [ -77.427564, 38.769108 ], [ -77.427735, 38.769376 ], [ -77.428079, 38.769677 ], [ -77.428422, 38.769878 ], [ -77.428808, 38.769978 ] ] } }
 ,
-{ "type": "Feature", "properties": { "LINEARID": "110410477374", "FULLNAME": "Fairway Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.427564, 38.769108 ], [ -77.427649, 38.769242 ], [ -77.427735, 38.769376 ], [ -77.427821, 38.769476 ], [ -77.427950, 38.769577 ], [ -77.428079, 38.769677 ], [ -77.428207, 38.769778 ], [ -77.428293, 38.769811 ], [ -77.428422, 38.769878 ], [ -77.428551, 38.769911 ], [ -77.428722, 38.769945 ], [ -77.428808, 38.769978 ], [ -77.430611, 38.769911 ], [ -77.431340, 38.769711 ] ] } }
+{ "type": "Feature", "properties": { "LINEARID": "110410477374", "FULLNAME": "Fairway Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.427564, 38.769108 ], [ -77.427735, 38.769376 ], [ -77.428079, 38.769677 ], [ -77.428422, 38.769878 ], [ -77.428808, 38.769978 ], [ -77.430611, 38.769911 ], [ -77.431340, 38.769711 ] ] } }
 ,
 { "type": "Feature", "properties": { "LINEARID": "110410477279", "FULLNAME": "Payne Ct", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -77.428508, 38.767302 ], [ -77.428594, 38.766565 ], [ -77.428551, 38.766197 ] ] } }
 ,

tile-join.cpp

@@ -378,7 +378,7 @@ double max(double a, double b) {
 
 struct tilecmp {
 	bool operator()(std::pair<unsigned, unsigned> const &a, std::pair<unsigned, unsigned> const &b) {
-		// must match behavior of reader::operator<()
+		// must match behavior of tileset_reader::operator<()
 
 		if (a.first < b.first) {
 			return true;
@@ -394,13 +394,13 @@ struct tilecmp {
 	}
 } tilecmp;
 
-// The `reader` is an iterator through the tiles of a tileset,
+// The `tileset_reader` is an iterator through the tiles of a tileset,
 // in z/x/tms_y order.
 //
 // The basic idea is that it is used like this:
 //
 //  void blah(const char *fname) {
-//      reader r(fname);
+//      tileset_reader r(fname);
 //
 //      for (; !r.all_done(); r.advance()) {
 //          std::pair<zxy, std::string> tile = r.current();
@@ -411,13 +411,13 @@ struct tilecmp {
 //  }
 //
 // The complication is that you can actually keep calling current()
-// and advance() after the reader claims to be done, in which case
+// and advance() after the tileset_reader claims to be done, in which case
 // it will produce overzoomed tiles generated from the tiles in
 // the maxzoom tileset. The parent tiles for those overzoomed tiles
 // are retrieved internally using get_tile() rather than through the
 // main iteration query.
 
-struct reader {
+struct tileset_reader {
 	// z/x/y and data of the current tile
 	long long zoom = 0;
 	long long x = 0;
@@ -438,7 +438,7 @@ struct reader {
 	// for iterating mbtiles
 	sqlite3 *db = NULL;
 	sqlite3_stmt *stmt = NULL;
-	struct reader *next = NULL;
+	struct tileset_reader *next = NULL;
 
 	// for iterating dirtiles
 	std::vector<zxy> dirtiles;
@@ -449,7 +449,7 @@ struct reader {
 	char *pmtiles_map = NULL;
 	std::vector<pmtiles::entry_zxy> pmtiles_entries;
 
-	reader(const char *fname) {
+	tileset_reader(const char *fname) {
 		name = fname;
 		struct stat st;
 		if (stat(fname, &st) == 0 && (st.st_mode & S_IFDIR) != 0) {
@@ -499,7 +499,7 @@ struct reader {
 		}
 	}
 
-	// Checks the done status not only of this reader but also
+	// Checks the done status not only of this tileset_reader but also
 	// the others chained to it in the queue.
 	//
 	// Also claims not to be done if at least one overzoomed tile
@@ -513,7 +513,7 @@ struct reader {
 			return false;
 		}
 
-		for (struct reader *r = next; r != NULL; r = r->next) {
+		for (struct tileset_reader *r = next; r != NULL; r = r->next) {
 			if (!r->done) {
 				return false;
 			}
@@ -688,7 +688,7 @@ struct reader {
 
 	// Sort in z/x/tms_y order, because that is the order of the
 	// straightforward query of the mbtiles tiles table.
-	bool operator<(const struct reader &r) const {
+	bool operator<(const struct tileset_reader &r) const {
 		// must match behavior of tilecmp
 
 		if (zoom < r.zoom) {
@@ -752,10 +752,10 @@ struct reader {
 	}
 };
 
-struct reader *begin_reading(char *fname) {
+struct tileset_reader *begin_reading(char *fname) {
-	struct reader *r = new reader(fname);
+	struct tileset_reader *r = new tileset_reader(fname);
 
-	// The reason this prefetches is so the reader queue can be
+	// The reason this prefetches is so the tileset_reader queue can be
 	// priority-ordered, so the one with the next relevant tile
 	// is first in line.
 	r->advance();
@@ -776,7 +776,7 @@ struct arg {
 	std::set<std::string> *remove_layers = NULL;
 	int ifmatched = 0;
 	json_object *filter = NULL;
-	struct reader *readers = NULL;
+	struct tileset_reader *readers = NULL;
 };
 
 void *join_worker(void *v) {
@@ -821,7 +821,7 @@ void *join_worker(void *v) {
 	return NULL;
 }
 
-void dispatch_tasks(std::map<zxy, std::vector<std::string>> &tasks, std::vector<std::map<std::string, layermap_entry>> &layermaps, sqlite3 *outdb, const char *outdir, std::vector<std::string> &header, std::map<std::string, std::vector<std::string>> &mapping, std::set<std::string> &exclude, std::set<std::string> &include, int ifmatched, std::set<std::string> &keep_layers, std::set<std::string> &remove_layers, json_object *filter, struct reader *readers) {
+void dispatch_tasks(std::map<zxy, std::vector<std::string>> &tasks, std::vector<std::map<std::string, layermap_entry>> &layermaps, sqlite3 *outdb, const char *outdir, std::vector<std::string> &header, std::map<std::string, std::vector<std::string>> &mapping, std::set<std::string> &exclude, std::set<std::string> &include, int ifmatched, std::set<std::string> &keep_layers, std::set<std::string> &remove_layers, json_object *filter, struct tileset_reader *readers) {
 	pthread_t pthreads[CPUS];
 	std::vector<arg> args;
 
@@ -969,7 +969,7 @@ void handle_vector_layers(json_object *vector_layers, std::map<std::string, laye
 	}
 }
 
-void decode(struct reader *readers, std::map<std::string, layermap_entry> &layermap, sqlite3 *outdb, const char *outdir, struct stats *st, std::vector<std::string> &header, std::map<std::string, std::vector<std::string>> &mapping, std::set<std::string> &exclude, std::set<std::string> &include, int ifmatched, std::string &attribution, std::string &description, std::set<std::string> &keep_layers, std::set<std::string> &remove_layers, std::string &name, json_object *filter, std::map<std::string, std::string> &attribute_descriptions, std::string &generator_options, std::vector<strategy> *strategies) {
+void decode(struct tileset_reader *readers, std::map<std::string, layermap_entry> &layermap, sqlite3 *outdb, const char *outdir, struct stats *st, std::vector<std::string> &header, std::map<std::string, std::vector<std::string>> &mapping, std::set<std::string> &exclude, std::set<std::string> &include, int ifmatched, std::string &attribution, std::string &description, std::set<std::string> &keep_layers, std::set<std::string> &remove_layers, std::string &name, json_object *filter, std::map<std::string, std::string> &attribute_descriptions, std::string &generator_options, std::vector<strategy> *strategies) {
 	std::vector<std::map<std::string, layermap_entry>> layermaps;
 	for (size_t i = 0; i < CPUS; i++) {
 		layermaps.push_back(std::map<std::string, layermap_entry>());
@@ -1028,21 +1028,21 @@ void decode(struct reader *readers, std::map<std::string, layermap_entry> &layer
 			}
 		}
 
-		// The reason this prefetches is so the reader queue can be
+		// The reason this prefetches is so the tileset_reader queue can be
 		// priority-ordered, so the one with the next relevant tile
 		// is first in line.
 		readers->advance();
 
-		// pull the reader off the front of the queue for reordering
+		// pull the tileset_reader off the front of the queue for reordering
 
-		reader *r = readers;
+		tileset_reader *r = readers;
 		readers = readers->next;
 		r->next = NULL;
 
-		// put the reader back onto the queue,
+		// put the tileset_reader back onto the queue,
 		// in whatever sequence its next tile calls for
 
-		struct reader **rr;
+		struct tileset_reader **rr;
 		for (rr = &readers; *rr != NULL; rr = &((*rr)->next)) {
 			if (*r < **rr) {
 				break;
@@ -1066,8 +1066,8 @@ void decode(struct reader *readers, std::map<std::string, layermap_entry> &layer
 	dispatch_tasks(tasks, layermaps, outdb, outdir, header, mapping, exclude, include, ifmatched, keep_layers, remove_layers, filter, readers);
 	layermap = merge_layermaps(layermaps);
 
-	struct reader *next;
+	struct tileset_reader *next;
-	for (struct reader *r = readers; r != NULL; r = next) {
+	for (struct tileset_reader *r = readers; r != NULL; r = next) {
 		next = r->next;
 		r->close();
 
@@ -1218,7 +1218,7 @@ int main(int argc, char **argv) {
 	int filearg = 0;
 	json_object *filter = NULL;
 
-	struct reader *readers = NULL;
+	struct tileset_reader *readers = NULL;
 
 	CPUS = sysconf(_SC_NPROCESSORS_ONLN);
 
@@ -1417,10 +1417,10 @@ int main(int argc, char **argv) {
 				filearg = 1;
 				while (getline(read_file, sa)) {
 					char *c = const_cast<char *>(sa.c_str());
-					reader *r = begin_reading(c);
+					tileset_reader *r = begin_reading(c);
 
-					// put the new reader in priority order
+					// put the new tileset_reader in priority order
-					struct reader **rr;
+					struct tileset_reader **rr;
 					for (rr = &readers; *rr != NULL; rr = &((*rr)->next)) {
 						if (*r < **rr) {
 							break;
@@ -1509,10 +1509,10 @@ int main(int argc, char **argv) {
 
 	if (filearg == 0) {
 		for (i = optind; i < argc; i++) {
-			reader *r = begin_reading(argv[i]);
+			tileset_reader *r = begin_reading(argv[i]);
 
-			// put the new reader in priority order
+			// put the new tileset_reader in priority order
-			struct reader **rr;
+			struct tileset_reader **rr;
 			for (rr = &readers; *rr != NULL; rr = &((*rr)->next)) {
 				if (*r < **rr) {
 					break;

tile.cpp

@@ -480,7 +480,10 @@ struct partial_arg {
 	std::vector<struct partial> *partials = NULL;
 	int task = 0;
 	int tasks = 0;
+
 	drawvec *shared_nodes;
+	node *shared_nodes_map;
+	size_t nodepos;
 };
 
 drawvec revive_polygon(drawvec &geom, double area, int z, int detail) {
@@ -524,7 +527,7 @@ drawvec revive_polygon(drawvec &geom, double area, int z, int detail) {
 	}
 }
 
-double simplify_partial(partial *p, drawvec &shared_nodes) {
+double simplify_partial(partial *p, drawvec const &shared_nodes, node *shared_nodes_map, size_t nodepos) {
 	drawvec geom;
 
 	for (size_t j = 0; j < p->geoms.size(); j++) {
@@ -579,7 +582,7 @@ double simplify_partial(partial *p, drawvec &shared_nodes) {
 				}
 
 				// continues to simplify to line_detail even if we have extra detail
-				drawvec ngeom = simplify_lines(geom, z, line_detail, !(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), p->simplification, t == VT_POLYGON ? 4 : 0, shared_nodes);
+				drawvec ngeom = simplify_lines(geom, z, p->tx, p->ty, line_detail, !(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), p->simplification, t == VT_POLYGON ? 4 : 0, shared_nodes, shared_nodes_map, nodepos);
 
 				if (t != VT_POLYGON || ngeom.size() >= 3) {
 					geom = ngeom;
@@ -602,7 +605,7 @@ void *partial_feature_worker(void *v) {
 	std::vector<struct partial> *partials = a->partials;
 
 	for (size_t i = a->task; i < (*partials).size(); i += a->tasks) {
-		double area = simplify_partial(&((*partials)[i]), *(a->shared_nodes));
+		double area = simplify_partial(&((*partials)[i]), *(a->shared_nodes), a->shared_nodes_map, a->nodepos);
 
 		signed char t = (*partials)[i].t;
 		int z = (*partials)[i].z;
@@ -1031,7 +1034,8 @@ bool find_common_edges(std::vector<partial> &partials, int z, int line_detail, d
 			}
 		}
 		if (!(prevent[P_SIMPLIFY] || (z == maxzoom && prevent[P_SIMPLIFY_LOW]) || (z < maxzoom && additional[A_GRID_LOW_ZOOMS]))) {
-			simplified_arcs[ai->second] = simplify_lines(dv, z, line_detail, !(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), simplification, 4, drawvec());
+			// tx and ty are 0 here because we aren't trying to do anything with the shared_nodes_map
+			simplified_arcs[ai->second] = simplify_lines(dv, z, 0, 0, line_detail, !(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), simplification, 4, drawvec(), NULL, 0);
 		} else {
 			simplified_arcs[ai->second] = dv;
 		}
@@ -1294,48 +1298,6 @@ long long choose_minextent(std::vector<long long> &extents, double f) {
 	return extents[(extents.size() - 1) * (1 - f)];
 }
 
-struct joint {
-	long long x : 34;  // enough to wrap around the world either way
-	unsigned long long p1 : 64 - 34;
-
-	long long y : 33;  // enough to touch the top and bottom of the world
-	unsigned long long p2 : 64 - 34;
-
-	joint(draw one, draw hinge, draw two) {
-		if (one < two) {
-			std::swap(one, two);
-		}
-
-		long long coord1[2] = {one.x, one.y};
-		long long coord2[2] = {two.x, two.y};
-		p1 = fnv1a(std::string((const char *) &coord1, sizeof(coord1)));
-		p2 = fnv1a(std::string((const char *) &coord2, sizeof(coord2)));
-
-		x = hinge.x;
-		y = hinge.y;
-	}
-
-	bool operator<(const joint &o) const {
-		if (y < o.y) {
-			return true;
-		} else if (y == o.y) {
-			if (x < o.x) {
-				return true;
-			} else if (x == o.x) {
-				if (p1 < o.p1) {
-					return true;
-				} else if (p1 == o.p1) {
-					if (p2 < o.p2) {
-						return true;
-					}
-				}
-			}
-		}
-
-		return false;
-	}
-};
-
 struct write_tile_args {
 	struct task *tasks = NULL;
 	char *stringpool = NULL;
@@ -1387,6 +1349,8 @@ struct write_tile_args {
 	atomic_strategy *strategy = NULL;
 	int zoom = -1;
 	bool compressed;
+	struct node *shared_nodes_map;
+	size_t nodepos;
 };
 
 bool clip_to_tile(serial_feature &sf, int z, long long buffer) {
@@ -1934,7 +1898,7 @@ void add_sample_to(std::vector<T> &vals, T val, size_t &increment, size_t seq) {
 	}
 }
 
-long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, char *stringpool, int z, const unsigned tx, const unsigned ty, const int detail, int min_detail, sqlite3 *outdb, const char *outdir, int buffer, const char *fname, compressor **geomfile, int minzoom, int maxzoom, double todo, std::atomic<long long> *along, long long alongminus, double gamma, int child_shards, long long *pool_off, unsigned *initial_x, unsigned *initial_y, std::atomic<int> *running, double simplification, std::vector<std::map<std::string, layermap_entry>> *layermaps, std::vector<std::vector<std::string>> *layer_unmaps, size_t tiling_seg, size_t pass, unsigned long long mingap, long long minextent, double fraction, const char *prefilter, const char *postfilter, struct json_object *filter, write_tile_args *arg, atomic_strategy *strategy, bool compressed_input) {
+long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, char *stringpool, int z, const unsigned tx, const unsigned ty, const int detail, int min_detail, sqlite3 *outdb, const char *outdir, int buffer, const char *fname, compressor **geomfile, int minzoom, int maxzoom, double todo, std::atomic<long long> *along, long long alongminus, double gamma, int child_shards, long long *pool_off, unsigned *initial_x, unsigned *initial_y, std::atomic<int> *running, double simplification, std::vector<std::map<std::string, layermap_entry>> *layermaps, std::vector<std::vector<std::string>> *layer_unmaps, size_t tiling_seg, size_t pass, unsigned long long mingap, long long minextent, double fraction, const char *prefilter, const char *postfilter, struct json_object *filter, write_tile_args *arg, atomic_strategy *strategy, bool compressed_input, struct node *shared_nodes_map, size_t nodepos) {
 	double merge_fraction = 1;
 	double mingap_fraction = 1;
 	double minextent_fraction = 1;
@@ -1990,7 +1954,6 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 
 		double coalesced_area = 0;
 		drawvec shared_nodes;
-		std::vector<joint> shared_joints;
 
 		int tile_detail = line_detail;
 		size_t skipped = 0;
@@ -2272,89 +2235,8 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 				} else {
 					kept++;
 
-					if (prevent[P_SIMPLIFY_SHARED_NODES]) {
+					for (auto &p : sf.edge_nodes) {
-						if (sf.t == VT_LINE) {
+						shared_nodes.push_back(p);
-							for (auto &g : sf.geometry) {
-								shared_nodes.push_back(g);
-							}
-						} else if (sf.t == VT_POLYGON) {
-							sf.geometry = remove_noop(sf.geometry, sf.t, 0);
-
-							for (size_t i = 0; i < sf.geometry.size(); i++) {
-								if (sf.geometry[i].op == VT_MOVETO) {
-									size_t j;
-									for (j = i + 1; j < sf.geometry.size(); j++) {
-										if (sf.geometry[j].op != VT_LINETO) {
-											break;
-										}
-									}
-
-									// j - 1 because we don't want the duplicate last point
-									for (size_t k = i; k < j - 1; k++) {
-										shared_joints.emplace_back(
-											sf.geometry[k],
-											sf.geometry[(k + 1 - i) % (j - 1 - i) + i],
-											sf.geometry[(k + 2 - i) % (j - 1 - i) + i]);
-									}
-
-									// since the starting point is never simplified away,
-									// don't let it be simplified away in any other polygons either.
-									// Needs to appear twice here so that the check below will see
-									// it as appearing in multiple features.
-									shared_nodes.push_back(sf.geometry[i]);
-									shared_nodes.push_back(sf.geometry[i]);
-
-									// To avoid letting polygons get simplified away to nothing,
-									// also keep the furthest-away point from the initial point
-									// (which Douglas-Peucker simplification would keep anyway,
-									// if its search weren't being split up by polygon side).
-
-									double far = 0;
-									size_t which = i;
-									for (size_t k = i + 1; k < j - 1; k++) {
-										double xd = sf.geometry[k].x - sf.geometry[i].x;
-										double yd = sf.geometry[k].y - sf.geometry[i].y;
-										double d = xd * xd + yd * yd;
-										if (d > far ||
-										    ((d == far) && (sf.geometry[k] < sf.geometry[which]))) {
-											far = d;
-											which = k;
-										}
-									}
-
-									shared_nodes.push_back(sf.geometry[which]);
-									shared_nodes.push_back(sf.geometry[which]);
-
-									// And, likewise, the point most distant from those two points,
-									// which probably would also be the one that Douglas-Peucker
-									// would keep next.
-
-									far = 0;
-									size_t which2 = i;
-
-									for (size_t k = i + 1; k < j - 1; k++) {
-										double d = distance_from_line(sf.geometry[k].x, sf.geometry[k].y,
-													      sf.geometry[i].x, sf.geometry[i].y,
-													      sf.geometry[which].x, sf.geometry[which].y);
-										if ((d > far) ||
-										    ((d == far) && (sf.geometry[k] < sf.geometry[which2]))) {
-											far = d;
-											which2 = k;
-										}
-									}
-
-									shared_nodes.push_back(sf.geometry[which2]);
-									shared_nodes.push_back(sf.geometry[which2]);
-
-									i = j - 1;
-								}
-							}
-						}
-
-						for (auto &p : sf.edge_nodes) {
-							shared_nodes.push_back(p);
-							shared_nodes.push_back(p);
-						}
 					}
 
 					partial p;
@@ -2398,11 +2280,15 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 					partials.push_back(p);
 
 					unsimplified_geometry_size += sf.geometry.size() * sizeof(draw);
-					if (unsimplified_geometry_size > 10 * 1024 * 1024 && !additional[A_DETECT_SHARED_BORDERS] && !prevent[P_SIMPLIFY_SHARED_NODES]) {
+					if (unsimplified_geometry_size > 10 * 1024 * 1024 && !additional[A_DETECT_SHARED_BORDERS]) {
+						// we should be safe to simplify here with P_SIMPLIFY_SHARED_NODES, since they will
+						// have been assembled globally, although that also means that simplification
+						// may not be very effective for reducing memory usage.
+
 						drawvec dv;
 
 						for (; simplified_geometry_through < partials.size(); simplified_geometry_through++) {
-							simplify_partial(&partials[simplified_geometry_through], dv);
+							simplify_partial(&partials[simplified_geometry_through], dv, shared_nodes_map, nodepos);
 
 							for (auto &g : partials[simplified_geometry_through].geoms) {
 								if (partials[simplified_geometry_through].t == VT_POLYGON) {
@@ -2421,38 +2307,7 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 			coalesced_area = 0;
 		}
 
-		{
+		std::sort(shared_nodes.begin(), shared_nodes.end());
-			drawvec just_shared_nodes;
-			std::sort(shared_nodes.begin(), shared_nodes.end());
-
-			for (size_t i = 0; i + 1 < shared_nodes.size(); i++) {
-				if (shared_nodes[i] == shared_nodes[i + 1]) {
-					just_shared_nodes.push_back(shared_nodes[i]);
-
-					draw d = shared_nodes[i];
-					i++;  // consume the first, point at the duplicate
-					while (i + 1 < shared_nodes.size() && shared_nodes[i + 1] == d) {
-						i++;
-					}
-				}
-			}
-
-			shared_nodes = just_shared_nodes;
-
-			std::sort(shared_joints.begin(), shared_joints.end());
-			for (size_t i = 0; i + 1 < shared_joints.size(); i++) {
-				if (shared_joints[i].x == shared_joints[i + 1].x &&
-				    shared_joints[i].y == shared_joints[i + 1].y) {
-					if (shared_joints[i].p1 != shared_joints[i + 1].p1 ||
-					    shared_joints[i].p2 != shared_joints[i + 1].p2) {
-						shared_nodes.push_back(draw(VT_MOVETO, shared_joints[i].x, shared_joints[i].y));
-					}
-				}
-			}
-
-			std::sort(shared_nodes.begin(), shared_nodes.end());
-			shared_joints.clear();
-		}
 
 		for (size_t i = 0; i < partials.size(); i++) {
 			partial &p = partials[i];
@@ -2552,6 +2407,8 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 			args[i].tasks = tasks;
 			args[i].partials = &partials;
 			args[i].shared_nodes = &shared_nodes;
+			args[i].shared_nodes_map = shared_nodes_map;
+			args[i].nodepos = nodepos;
 
 			if (tasks > 1) {
 				if (pthread_create(&pthreads[i], NULL, partial_feature_worker, &args[i]) != 0) {
@@ -2667,8 +2524,9 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 				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);
 					if (!(prevent[P_SIMPLIFY] || (z == maxzoom && prevent[P_SIMPLIFY_LOW]))) {
-						layer_features[x].geom = simplify_lines(layer_features[x].geom, 32, 0,
+						// XXX revisit: why does this not take zoom into account?
-											!(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), simplification, layer_features[x].type == VT_POLYGON ? 4 : 0, shared_nodes);
+						layer_features[x].geom = simplify_lines(layer_features[x].geom, 32, 0, 0, 0,
+											!(prevent[P_CLIPPING] || prevent[P_DUPLICATION]), simplification, layer_features[x].type == VT_POLYGON ? 4 : 0, shared_nodes, NULL, 0);
 					}
 				}
 
@@ -3082,7 +2940,7 @@ void *run_thread(void *vargs) {
 
 			// fprintf(stderr, "%d/%u/%u\n", z, x, y);
 
-			long long len = write_tile(&dc, &geompos, arg->stringpool, z, x, y, z == arg->maxzoom ? arg->full_detail : arg->low_detail, arg->min_detail, arg->outdb, arg->outdir, arg->buffer, arg->fname, arg->geomfile, arg->minzoom, arg->maxzoom, arg->todo, arg->along, geompos, arg->gamma, arg->child_shards, arg->pool_off, arg->initial_x, arg->initial_y, arg->running, arg->simplification, arg->layermaps, arg->layer_unmaps, arg->tiling_seg, arg->pass, arg->mingap, arg->minextent, arg->fraction, arg->prefilter, arg->postfilter, arg->filter, arg, arg->strategy, arg->compressed);
+			long long len = write_tile(&dc, &geompos, arg->stringpool, z, x, y, z == arg->maxzoom ? arg->full_detail : arg->low_detail, arg->min_detail, arg->outdb, arg->outdir, arg->buffer, arg->fname, arg->geomfile, arg->minzoom, arg->maxzoom, arg->todo, arg->along, geompos, arg->gamma, arg->child_shards, arg->pool_off, arg->initial_x, arg->initial_y, arg->running, arg->simplification, arg->layermaps, arg->layer_unmaps, arg->tiling_seg, arg->pass, arg->mingap, arg->minextent, arg->fraction, arg->prefilter, arg->postfilter, arg->filter, arg, arg->strategy, arg->compressed, arg->shared_nodes_map, arg->nodepos);
 
 			if (len < 0) {
 				int *err = &arg->err;
@@ -3147,7 +3005,7 @@ void *run_thread(void *vargs) {
 	return NULL;
 }
 
-int traverse_zooms(int *geomfd, off_t *geom_size, char *stringpool, std::atomic<unsigned> *midx, std::atomic<unsigned> *midy, int &maxzoom, int minzoom, sqlite3 *outdb, const char *outdir, int buffer, const char *fname, const char *tmpdir, double gamma, int full_detail, int low_detail, int min_detail, long long *pool_off, unsigned *initial_x, unsigned *initial_y, double simplification, double maxzoom_simplification, std::vector<std::map<std::string, layermap_entry>> &layermaps, const char *prefilter, const char *postfilter, std::map<std::string, attribute_op> const *attribute_accum, struct json_object *filter, std::vector<strategy> &strategies, int iz) {
+int traverse_zooms(int *geomfd, off_t *geom_size, char *stringpool, std::atomic<unsigned> *midx, std::atomic<unsigned> *midy, int &maxzoom, int minzoom, sqlite3 *outdb, const char *outdir, int buffer, const char *fname, const char *tmpdir, double gamma, int full_detail, int low_detail, int min_detail, long long *pool_off, unsigned *initial_x, unsigned *initial_y, double simplification, double maxzoom_simplification, std::vector<std::map<std::string, layermap_entry>> &layermaps, const char *prefilter, const char *postfilter, std::map<std::string, attribute_op> const *attribute_accum, struct json_object *filter, std::vector<strategy> &strategies, int iz, struct node *shared_nodes_map, size_t nodepos) {
 	last_progress = 0;
 
 	// The existing layermaps are one table per input thread.
@@ -3352,6 +3210,8 @@ int traverse_zooms(int *geomfd, off_t *geom_size, char *stringpool, std::atomic<
 				args[thread].strategy = &strategy;
 				args[thread].zoom = z;
 				args[thread].compressed = (z != iz);
+				args[thread].shared_nodes_map = shared_nodes_map;
+				args[thread].nodepos = nodepos;
 
 				if (pthread_create(&pthreads[thread], NULL, run_thread, &args[thread]) != 0) {
 					perror("pthread_create");

tile.hpp

@@ -7,6 +7,7 @@
 #include <atomic>
 #include <map>
 #include "mbtiles.hpp"
+#include "serial.hpp"
 #include "jsonpull/jsonpull.h"
 
 enum attribute_op {
@@ -61,9 +62,9 @@ struct strategy {
 	strategy() = default;
 };
 
-long long write_tile(char **geom, char *stringpool, unsigned *file_bbox, int z, unsigned x, unsigned y, int detail, int min_detail, int basezoom, sqlite3 *outdb, const char *outdir, double droprate, int buffer, const char *fname, FILE **geomfile, int file_minzoom, int file_maxzoom, double todo, char *geomstart, long long along, double gamma, int nlayers, std::atomic<strategy> *strategy);
+// long long write_tile(char **geom, char *stringpool, unsigned *file_bbox, int z, unsigned x, unsigned y, int detail, int min_detail, int basezoom, sqlite3 *outdb, const char *outdir, double droprate, int buffer, const char *fname, FILE **geomfile, int file_minzoom, int file_maxzoom, double todo, char *geomstart, long long along, double gamma, int nlayers, std::atomic<strategy> *strategy);
 
-int traverse_zooms(int *geomfd, off_t *geom_size, char *stringpool, std::atomic<unsigned> *midx, std::atomic<unsigned> *midy, int &maxzoom, int minzoom, sqlite3 *outdb, const char *outdir, int buffer, const char *fname, const char *tmpdir, double gamma, int full_detail, int low_detail, int min_detail, long long *pool_off, unsigned *initial_x, unsigned *initial_y, double simplification, double maxzoom_simplification, std::vector<std::map<std::string, layermap_entry> > &layermap, const char *prefilter, const char *postfilter, std::map<std::string, attribute_op> const *attribute_accum, struct json_object *filter, std::vector<strategy> &strategies, int iz);
+int traverse_zooms(int *geomfd, off_t *geom_size, char *stringpool, std::atomic<unsigned> *midx, std::atomic<unsigned> *midy, int &maxzoom, int minzoom, sqlite3 *outdb, const char *outdir, int buffer, const char *fname, const char *tmpdir, double gamma, int full_detail, int low_detail, int min_detail, long long *pool_off, unsigned *initial_x, unsigned *initial_y, double simplification, double maxzoom_simplification, std::vector<std::map<std::string, layermap_entry> > &layermap, const char *prefilter, const char *postfilter, std::map<std::string, attribute_op> const *attribute_accum, struct json_object *filter, std::vector<strategy> &strategies, int iz, struct node *shared_nodes_map, size_t nodepos);
 
 int manage_gap(unsigned long long index, unsigned long long *previndex, double scale, double gamma, double *gap);
 

unit.cpp

@@ -1,6 +1,9 @@
 #define CATCH_CONFIG_MAIN
 #include "catch/catch.hpp"
 #include "text.hpp"
+#include "sort.hpp"
+#include <unistd.h>
+#include <limits.h>
 
 TEST_CASE("UTF-8 enforcement", "[utf8]") {
 	REQUIRE(check_utf8("") == std::string(""));
@@ -18,3 +21,47 @@ TEST_CASE("UTF-8 truncation", "[trunc]") {
 	REQUIRE(truncate16("0123456789😀😬😁😂😃😄😅😆", 17) == std::string("0123456789😀😬😁"));
 	REQUIRE(truncate16("0123456789あいうえおかきくけこさ", 16) == std::string("0123456789あいうえおか"));
 }
+
+int intcmp(const void *v1, const void *v2) {
+	return *((int *) v1) - *((int *) v2);
+}
+
+TEST_CASE("External quicksort", "fqsort") {
+	std::vector<FILE *> inputs;
+
+	size_t written = 0;
+	for (size_t i = 0; i < 5; i++) {
+		std::string tmpname = "/tmp/in.XXXXXXX";
+		int fd = mkstemp((char *) tmpname.c_str());
+		unlink(tmpname.c_str());
+		FILE *f = fdopen(fd, "w+b");
+		inputs.emplace_back(f);
+		size_t iterations = 2000 + rand() % 200;
+		for (size_t j = 0; j < iterations; j++) {
+			int n = rand();
+			fwrite((void *) &n, sizeof(int), 1, f);
+			written++;
+		}
+		rewind(f);
+	}
+
+	std::string tmpname = "/tmp/out.XXXXXX";
+	int fd = mkstemp((char *) tmpname.c_str());
+	unlink(tmpname.c_str());
+	FILE *f = fdopen(fd, "w+b");
+
+	fqsort(inputs, sizeof(int), intcmp, f, 256);
+	rewind(f);
+
+	int prev = INT_MIN;
+	int here;
+	size_t nread = 0;
+	while (fread((void *) &here, sizeof(int), 1, f)) {
+		REQUIRE(here >= prev);
+		prev = here;
+		nread++;
+	}
+
+	fclose(f);
+	REQUIRE(nread == written);
+}

version.hpp

@@ -1,6 +1,6 @@
 #ifndef VERSION_HPP
 #define VERSION_HPP
 
-#define VERSION "v2.32.1"
+#define VERSION "v2.33.0"
 
 #endif