I'm not sure how this ever worked before

drop.cpp

@@ -1,36 +1,39 @@
+#include <cmath>
 #include "drop.hpp"
 #include "options.hpp"
 #include "geometry.hpp"
 
 unsigned long long preserve_point_density_threshold = 0;
 
-int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, double gamma) {
+int calc_feature_minzoom(struct index *ix, struct drop_state ds[], int maxzoom, double gamma) {
 	int 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--) {
+		for (ssize_t i = 0; i <= maxzoom; i++) {
-			ds[i].seq++;
+			// This zoom level is now lighter on features than it should be.
+			ds[i].error -= 1.0 / ds[i].interval;
+			// printf("z%zd: error %f with interval %f\n", i, ds[i].error, ds[i].interval);
 		}
+
 		ssize_t chosen = maxzoom + 1;
-		for (ssize_t i = maxzoom; i >= 0; i--) {
+		for (ssize_t i = 0; i <= maxzoom; i++) {
-			if (ds[i].seq < 0) {
+			if (ds[i].error < 0) {
-				feature_minzoom = i + 1;
+				// this zoom level is too light, so it is time to emit a feature.
+				feature_minzoom = i;
 
-				// The feature we are pushing out
+				// this feature now appears in this zoom level and all higher zoom levels,
-				// appears in zooms i + 1 through maxzoom,
+				// so each of them has this feature as its last feature, and each of them
-				// so track where that was so we can make sure
+				// is now one feature heavier than before.
-				// not to cluster something else that is *too*
+				for (ssize_t j = i; j <= maxzoom; j++) {
-				// far away into it.
-				for (ssize_t j = i + 1; j <= maxzoom; j++) {
 					ds[j].previndex = ix->ix;
+					ds[j].error += ds[j].interval / ds[j].interval;
+					// printf("z%zd: now error %f\n", j, ds[j].error);
 				}
 
-				chosen = i + 1;
+				chosen = i;
 				break;
-			} else {
-				ds[i].seq -= ds[i].interval;
 			}
 		}
 
@@ -44,8 +47,12 @@ int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, d
 				if (ix->ix - ds[i].previndex > ((1LL << (32 - i)) / preserve_point_density_threshold) * ((1LL << (32 - i)) / preserve_point_density_threshold)) {
 					feature_minzoom = i;
 
-					for (ssize_t j = i; j <= maxzoom; j++) {
+					// this feature now appears in this zoom level and all higher zoom levels below `chosen`,
+					// so each of them has this feature as its last feature, and each of them
+					// is now one feature heavier than before.
+					for (ssize_t j = i; j < chosen; j++) {
 						ds[j].previndex = ix->ix;
+						ds[j].error += ds[j].interval / ds[j].interval;
 					}
 
 					break;
@@ -56,3 +63,32 @@ int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, d
 
 	return feature_minzoom;
 }
+
+void prep_drop_states(struct drop_state ds[], int maxzoom, int basezoom, double droprate) {
+	if (basezoom < 0) {
+		basezoom = maxzoom;
+	}
+
+	// Needs to be signed for interval calculation
+	// printf("prep! max %d, base %d, rate %f\n", maxzoom, basezoom, droprate);
+	for (ssize_t i = 0; i <= maxzoom; i++) {
+		ds[i].previndex = 0;
+		ds[i].interval = 1;  // every feature appears in every zoom level at or above the basezoom
+
+		if (i < basezoom) {
+			// at zoom levels below the basezoom, the fraction of points that are dropped is
+			// the drop rate to the power of the number of zooms this zoom is below the basezoom
+			//
+			// for example:
+			// basezoom:       1    (droprate ^ 0)
+			// basezoom - 1:   2.5  (droprate ^ 1)
+			// basezoom - 2:   6.25 (droprate ^ 2)
+			// ...
+			// basezoom - n:        (droprate ^ n)
+			ds[i].interval = std::exp(std::log(droprate) * (basezoom - i));
+			// printf("%zd: interval %f\n", i, ds[i].interval);
+		}
+
+		ds[i].error = 0;
+	}
+}

drop.hpp

@@ -8,7 +8,7 @@
 // Note that the fields are in a specific order so that `segment` and `t` will
 // packed together with `seq` so that the total structure size will be only 32 bytes
 // instead of 40. (Could we save a few more, perhaps, by tracking `len` instead of
-// `end` and limiting the size of individual features to 32 bits?)
+// `end` and limiting the size of individual features to 2^32 bytes?)
 
 struct index {
 	// first and last+1 byte of the feature in the geometry temp file
@@ -25,7 +25,7 @@ struct index {
 	unsigned short t : 2;
 
 	// sequence number (sometimes with gaps in numbering) of the feature in the original input file
-	unsigned long long seq : (64 - 18);  // pack with segment and t to stay in 32 bytes
+	unsigned long long seq : (64 - 16 - 2);	 // pack with segment and t to stay in 32 bytes
 
 	index()
 	    : t(0),
@@ -33,15 +33,33 @@ struct index {
 	}
 };
 
+// Each zoom level has a drop_state that is used to account for the fraction of
+// point features that are supposed to be dropped in that zoom level. As it goes
+// through the spatially-sorted features, it is basically doing a diffusion dither
+// to keep the density of features in each vicinity at each zoom level
+// approximately correct by including or excluding individual features
+// to maintain the balance.
+
 struct drop_state {
-	double gap;
+	// the z-index or hilbert index of the last feature that was placed in this zoom level
 	unsigned long long previndex;
+
+	// the preservation rate (1 or more) for features in this zoom level.
+	// 1 would be to keep all the features; 2 would drop every other feature;
+	// 4 every fourth feature, and so on.
 	double interval;
-	double seq;  // floating point because interval is
+
+	// the current accumulated error in this zoom level:
+	// positive if too many features have been dropped;
+	// negative if not enough features have been dropped.
+	//
+	// this is floating-point because the interval is.
+	double error;
 };
 
 extern unsigned long long preserve_point_density_threshold;
 
-int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, double gamma);
+int calc_feature_minzoom(struct index *ix, struct drop_state ds[], int maxzoom, double gamma);
+void prep_drop_states(struct drop_state ds[], int maxzoom, int basezoom, double droprate);
 
 #endif

main.cpp

@@ -985,21 +985,6 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 	}
 }
 
-void prep_drop_states(struct drop_state *ds, int maxzoom, int basezoom, double droprate) {
-	// Needs to be signed for interval calculation
-	for (ssize_t i = 0; i <= maxzoom; i++) {
-		ds[i].gap = 0;
-		ds[i].previndex = 0;
-		ds[i].interval = 0;
-
-		if (i < basezoom) {
-			ds[i].interval = std::exp(std::log(droprate) * (basezoom - i));
-		}
-
-		ds[i].seq = 0;
-	}
-}
-
 static size_t calc_memsize() {
 	size_t mem;
 

unit.cpp

@@ -3,6 +3,8 @@
 #include "text.hpp"
 #include "drop.hpp"
 
+unsigned int additional[256] = {0};
+
 TEST_CASE("UTF-8 enforcement", "[utf8]") {
 	REQUIRE(check_utf8("") == std::string(""));
 	REQUIRE(check_utf8("hello world") == std::string(""));
@@ -24,4 +26,16 @@ TEST_CASE("index structure packing", "[index]") {
 	REQUIRE(sizeof(struct index) == 32);
 }
 
-unsigned int additional[256] = {0};
+TEST_CASE("prep drop states", "[prep_drop_state]") {
+	struct drop_state ds[25];
+
+	prep_drop_states(ds, 24, 16, 2);
+	REQUIRE(ds[24].interval == 1);
+	REQUIRE(ds[17].interval == 1);
+	REQUIRE(ds[16].interval == 1);
+	REQUIRE(ds[15].interval == 2);
+	REQUIRE(ds[14].interval == 4);
+
+	// to fix: because of floating point error this is not quite true
+	// REQUIRE(ds[0].interval == 65536);
+}