Calculate feature-dropping (except gamma) during geometry reordering

geojson.cpp

@@ -431,6 +431,7 @@ int serialize_geometry(json_object *geometry, json_object *properties, json_obje
 	index.end = *geompos;
 	index.segment = segment;
 	index.seq = *layer_seq;
+	index.t = sf.t;
 
 	// Calculate the center even if off the edge of the plane,
 	// and then mask to bring it back into the addressable area

main.cpp

@@ -24,6 +24,7 @@
 #include <string>
 #include <set>
 #include <map>
+#include <cmath>
 
 #ifdef __APPLE__
 #include <sys/types.h>
@@ -209,21 +210,62 @@ static void insert(struct mergelist *m, struct mergelist **head, unsigned char *
 	*head = m;
 }
 
-static void merge(struct mergelist *merges, int nmerges, unsigned char *map, FILE *f, int bytes, long long nrec, char *geom_map, FILE *geom_out, long long *geompos, long long *progress, long long *progress_max, long long *progress_reported) {
+static void merge(struct mergelist *merges, int nmerges, unsigned char *map, FILE *f, int bytes, long long nrec, char *geom_map, FILE *geom_out, long long *geompos, long long *progress, long long *progress_max, long long *progress_reported, int maxzoom, int basezoom, double droprate, double gamma) {
-	int i;
 	struct mergelist *head = NULL;
 
-	for (i = 0; i < nmerges; i++) {
+	for (size_t i = 0; i < nmerges; i++) {
 		if (merges[i].start < merges[i].end) {
 			insert(&(merges[i]), &head, map);
 		}
 	}
 
+	struct feature_drop {
+		double gap;
+		unsigned long long previndex;
+		double interval;
+		double scale;
+		double seq;
+	} fd[maxzoom + 1];
+
+	// Needs to be signed for interval calculation
+	for (ssize_t i = 0; i <= maxzoom; i++) {
+		fd[i].gap = 0;
+		fd[i].previndex = 0;
+		fd[i].interval = 0;
+
+		if (i < basezoom) {
+			fd[i].interval = std::exp(std::log(droprate) * (basezoom - i));
+		}
+
+		fd[i].scale = (double) (1LL << (64 - 2 * (i + 8)));
+		fd[i].seq = 0;
+	}
+
 	while (head != NULL) {
 		struct index *ix = (struct index *) (map + head->start);
 		fwrite_check(geom_map + ix->start, 1, ix->end - ix->start, geom_out, "merge geometry");
 		*geompos += ix->end - ix->start;
-		serialize_byte(geom_out, 0, geompos, "merge geometry");
+		double feature_minzoom = 0;
+
+		if (gamma >= 0 && (ix->t == VT_POINT ||
+				   (additional[A_LINE_DROP] && ix->t == VT_LINE) ||
+				   (additional[A_POLYGON_DROP] && ix->t == VT_POLYGON))) {
+			for (ssize_t i = maxzoom; i >= 0; i--) {
+				fd[i].seq++;
+			}
+			for (ssize_t i = maxzoom; i >= 0; i--) {
+				if (fd[i].seq >= 0) {
+					fd[i].seq -= fd[i].interval;
+				} else {
+					feature_minzoom = i + 1;
+					break;
+				}
+			}
+
+			// XXX manage_gap
+		}
+
+		serialize_byte(geom_out, feature_minzoom, geompos, "merge geometry");
 
 		// Count this as an 75%-accomplishment, since we already 25%-counted it
 		*progress += (ix->end - ix->start) * 3 / 4;
@@ -483,7 +525,7 @@ void start_parsing(int fd, FILE *fp, long long offset, long long len, volatile i
 	parser_created = true;
 }
 
-void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int splits, long long mem, const char *tmpdir, long long *availfiles, FILE *geomfile, FILE *indexfile, long long *geompos_out, long long *progress, long long *progress_max, long long *progress_reported) {
+void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int splits, long long mem, const char *tmpdir, long long *availfiles, FILE *geomfile, FILE *indexfile, long long *geompos_out, long long *progress, long long *progress_max, long long *progress_reported, int maxzoom, int basezoom, double droprate, double gamma) {
 	// Arranged as bits to facilitate subdividing again if a subdivided file is still huge
 	int splitbits = log(splits) / log(2);
 	splits = 1 << splitbits;
@@ -696,7 +738,7 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 				madvise(geommap, geomst.st_size, MADV_RANDOM);
 				madvise(geommap, geomst.st_size, MADV_WILLNEED);
 
-				merge(merges, nmerges, (unsigned char *) indexmap, indexfile, bytes, indexpos / bytes, geommap, geomfile, geompos_out, progress, progress_max, progress_reported);
+				merge(merges, nmerges, (unsigned char *) indexmap, indexfile, bytes, indexpos / bytes, geommap, geomfile, geompos_out, progress, progress_max, progress_reported, maxzoom, basezoom, droprate, gamma);
 
 				madvise(indexmap, indexst.st_size, MADV_DONTNEED);
 				if (munmap(indexmap, indexst.st_size) < 0) {
@@ -762,7 +804,7 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 				// counter backward but will be an honest estimate of the work remaining.
 				*progress_max += geomst.st_size / 4;
 
-				radix1(&geomfds[i], &indexfds[i], 1, prefix + splitbits, *availfiles / 4, mem, tmpdir, availfiles, geomfile, indexfile, geompos_out, progress, progress_max, progress_reported);
+				radix1(&geomfds[i], &indexfds[i], 1, prefix + splitbits, *availfiles / 4, mem, tmpdir, availfiles, geomfile, indexfile, geompos_out, progress, progress_max, progress_reported, maxzoom, basezoom, droprate, gamma);
 				already_closed = 1;
 			}
 		}
@@ -782,7 +824,7 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 	}
 }
 
-void radix(struct reader *reader, int nreaders, FILE *geomfile, int geomfd, FILE *indexfile, int indexfd, const char *tmpdir, long long *geompos) {
+void radix(struct reader *reader, int nreaders, FILE *geomfile, int geomfd, FILE *indexfile, int indexfd, const char *tmpdir, long long *geompos, int maxzoom, int basezoom, double droprate, double gamma) {
 	// Run through the index and geometry for each reader,
 	// splitting the contents out by index into as many
 	// sub-files as we can write to simultaneously.
@@ -849,7 +891,7 @@ void radix(struct reader *reader, int nreaders, FILE *geomfile, int geomfd, FILE
 
 	long long progress = 0, progress_max = geom_total, progress_reported = -1;
 	long long availfiles_before = availfiles;
-	radix1(geomfds, indexfds, nreaders, 0, splits, mem, tmpdir, &availfiles, geomfile, indexfile, geompos, &progress, &progress_max, &progress_reported);
+	radix1(geomfds, indexfds, nreaders, 0, splits, mem, tmpdir, &availfiles, geomfile, indexfile, geompos, &progress, &progress_max, &progress_reported, maxzoom, basezoom, droprate, gamma);
 
 	if (availfiles - 2 * nreaders != availfiles_before) {
 		fprintf(stderr, "Internal error: miscounted available file descriptors: %lld vs %lld\n", availfiles - 2 * nreaders, availfiles);
@@ -1384,7 +1426,7 @@ int read_input(std::vector<source> &sources, char *fname, const char *layername,
 	serialize_uint(geomfile, 0, &geompos, fname);
 	serialize_uint(geomfile, 0, &geompos, fname);
 
-	radix(reader, CPUS, geomfile, geomfd, indexfile, indexfd, tmpdir, &geompos);
+	radix(reader, CPUS, geomfile, geomfd, indexfile, indexfd, tmpdir, &geompos, maxzoom, basezoom, droprate, gamma);
 
 	/* end of tile */
 	serialize_byte(geomfile, -2, &geompos, fname);

main.hpp

@@ -3,7 +3,8 @@ struct index {
 	long long end;
 	unsigned long long index;
 	short segment;
-	unsigned long long seq : (64 - 16);  // pack with segment to stay in 32 bytes
+	short t : 2;
+	unsigned long long seq : (64 - 18);  // pack with segment and t to stay in 32 bytes
 };
 
 void checkdisk(struct reader *r, int nreader);

tile.cpp

@@ -792,11 +792,7 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 				continue;
 			}
 
-			if (t == VT_LINE && z + line_detail <= feature_minzoom) {
+			if (z < feature_minzoom) {
-				continue;
-			}
-
-			if (t == VT_POINT && z < feature_minzoom && gamma < 0) {
 				continue;
 			}
 
@@ -805,21 +801,10 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 				index = encode(bbox[0] / 2 + bbox[2] / 2, bbox[1] / 2 + bbox[3] / 2);
 			}
 
-			if (gamma >= 0 && (t == VT_POINT ||
+			if (gamma > 0) {
-					   (additional[A_LINE_DROP] && t == VT_LINE) ||
+				if (manage_gap(index, &previndex, scale, gamma, &gap)) {
-					   (additional[A_POLYGON_DROP] && t == VT_POLYGON))) {
-				seq++;
-				if (seq >= 0) {
-					seq -= interval;
-				} else {
 					continue;
 				}
-
-				if (gamma > 0) {
-					if (manage_gap(index, &previndex, scale, gamma, &gap)) {
-						continue;
-					}
-				}
 			}
 
 			if (additional[A_CALCULATE_FEATURE_DENSITY]) {