Add Hilbert Curve feature sequencing as an option

CHANGELOG.md

@@ -1,6 +1,7 @@
 ## 1.32.7
 
 * Point features may now be coalesced into MultiPoint features with --coalesce.
+* Add --hilbert option to put features in Hilbert Curve sequence
 
 ## 1.32.6
 

README.md

@@ -425,7 +425,7 @@ the same layer, enclose them in an `all` expression so they will all be evaluate
  * `-as` or `--drop-densest-as-needed`: If a tile is too large, try to reduce it to under 500K by increasing the minimum spacing between features. The discovered spacing applies to the entire zoom level.
  * `-ad` or `--drop-fraction-as-needed`: Dynamically drop some fraction of features from each zoom level to keep large tiles under the 500K size limit. (This is like `-pd` but applies to the entire zoom level, not to each tile.)
  * `-an` or `--drop-smallest-as-needed`: Dynamically drop the smallest features (physically smallest: the shortest lines or the smallest polygons) from each zoom level to keep large tiles under the 500K size limit. This option will not work for point features.
- * `-aN` or `--coalesce-smallest-as-needed`: Dynamically combine the smallest features (physically smallest: the shortest lines or the smallest polygons) from each zoom level into other nearby features to keep large tiles under the 500K size limit. This option will not work for point features, and will probably not help very much with LineStrings. It is mostly intended for polygons, to maintain the full original area covered by polygons while still reducing the feature count somehow. The attributes of the small polygons are *not* preserved into the combined features, only their geometry. Furthermore, the polygons that nested polygons are coalesced to may not necessarily be the immediately enclosing features.
+ * `-aN` or `--coalesce-smallest-as-needed`: Dynamically combine the smallest features (physically smallest: the shortest lines or the smallest polygons) from each zoom level into other nearby features to keep large tiles under the 500K size limit. This option will not work for point features, and will probably not help very much with LineStrings. It is mostly intended for polygons, to maintain the full original area covered by polygons while still reducing the feature count somehow. The attributes of the small polygons are *not* preserved into the combined features (except through `--accumulate-attribute`), only their geometry. Furthermore, the polygons to which nested polygons are coalesced may not necessarily be the immediately enclosing features.
  * `-aD` or `--coalesce-densest-as-needed`: Dynamically combine the densest features from each zoom level into other nearby features to keep large tiles under the 500K size limit. (Again, mostly useful for polygons.)
  * `-aS` or `--coalesce-fraction-as-needed`: Dynamically combine a fraction of features from each zoom level into other nearby features to keep large tiles under the 500K size limit. (Again, mostly useful for polygons.)
  * `-pd` or `--force-feature-limit`: Dynamically drop some fraction of features from large tiles to keep them under the 500K size limit. It will probably look ugly at the tile boundaries. (This is like `-ad` but applies to each tile individually, not to the entire zoom level.) You probably don't want to use this.
@@ -458,9 +458,10 @@ the same layer, enclose them in an `all` expression so they will all be evaluate
 ### Reordering features within each tile
 
  * `-pi` or `--preserve-input-order`: Preserve the original input order of features as the drawing order instead of ordering geographically. (This is implemented as a restoration of the original order at the end, so that dot-dropping is still geographic, which means it also undoes `-ao`).
- * `-ao` or `--reorder`: Reorder features to put ones with the same attributes in sequence, to try to get them to coalesce. You probably want to use this if you use `--coalesce`.
- * `-ac` or `--coalesce`: Coalesce adjacent features that have the same attributes. This can be useful if you have lots of small polygons with identical attributes and you would like to merge them together.
+ * `-ac` or `--coalesce`: Coalesce consecutive features that have the same attributes. This can be useful if you have lots of small polygons with identical attributes and you would like to merge them together.
+ * `-ao` or `--reorder`: Reorder features to put ones with the same attributes in sequence (instead of ones that are approximately spatially adjacent), to try to get them to coalesce. You probably want to use this if you use `--coalesce`.
  * `-ar` or `--reverse`: Try reversing the directions of lines to make them coalesce and compress better. You probably don't want to use this.
+ * `-ah` or `--hilbert`: Put features in Hilbert Curve order instead of the usual Z-Order. This improves the odds that spatially adjacent features will be sequentially adjacent, and should improve density calculations and spatial coalescing. It should be the default eventually.
 
 ### Adding calculated attributes
 
@@ -469,7 +470,7 @@ the same layer, enclose them in an `all` expression so they will all be evaluate
 
 ### Trying to correct bad source geometry
 
- * `-aw` or `--detect-longitude-wraparound`: Detect when adjacent points within a feature jump to the other side of the world, and try to fix the geometry.
+ * `-aw` or `--detect-longitude-wraparound`: Detect when consecutive points within a feature jump to the other side of the world, and try to fix the geometry.
  * `-pw` or `--use-source-polygon-winding`: Instead of respecting GeoJSON polygon ring order, use the original polygon winding in the source data to distinguish inner (clockwise) and outer (counterclockwise) polygon rings.
  * `-pW` or `--reverse-source-polygon-winding`: Instead of respecting GeoJSON polygon ring order, use the opposite of the original polygon winding in the source data to distinguish inner (counterclockwise) and outer (clockwise) polygon rings.
  * `--clip-bounding-box=`*minlon*`,`*minlat*`,`*maxlon*`,`*maxlat*: Clip all features to the specified bounding box.

geometry.cpp

@@ -874,8 +874,8 @@ drawvec reorder_lines(drawvec &geom) {
 	// instead of down and to the right
 	// so that it will coalesce better
 
-	unsigned long long l1 = encode(geom[0].x, geom[0].y);
-	unsigned long long l2 = encode(geom[geom.size() - 1].x, geom[geom.size() - 1].y);
+	unsigned long long l1 = encode_index(geom[0].x, geom[0].y);
+	unsigned long long l2 = encode_index(geom[geom.size() - 1].x, geom[geom.size() - 1].y);
 
 	if (l1 > l2) {
 		drawvec out;

main.cpp

@@ -269,7 +269,7 @@ struct drop_state {
 int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, double gamma) {
 	int feature_minzoom = 0;
 	unsigned xx, yy;
-	decode(ix->ix, &xx, &yy);
+	decode_index(ix->ix, &xx, &yy);
 
 	if (gamma >= 0 && (ix->t == VT_POINT ||
 			   (additional[A_LINE_DROP] && ix->t == VT_LINE) ||
@@ -2062,7 +2062,7 @@ int read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzo
 		long long ip;
 		for (ip = 0; ip < indices; ip++) {
 			unsigned xx, yy;
-			decode(map[ip].ix, &xx, &yy);
+			decode_index(map[ip].ix, &xx, &yy);
 
 			long long nprogress = 100 * ip / indices;
 			if (nprogress != progress) {
@@ -2590,6 +2590,7 @@ int main(int argc, char **argv) {
 		{"reorder", no_argument, &additional[A_REORDER], 1},
 		{"coalesce", no_argument, &additional[A_COALESCE], 1},
 		{"reverse", no_argument, &additional[A_REVERSE], 1},
+		{"hilbert", no_argument, &additional[A_HILBERT], 1},
 
 		{"Adding calculated attributes", 0, 0, 0},
 		{"calculate-feature-density", no_argument, &additional[A_CALCULATE_FEATURE_DENSITY], 1},
@@ -3022,6 +3023,14 @@ int main(int argc, char **argv) {
 		}
 	}
 
+	if (additional[A_HILBERT]) {
+		encode_index = encode_hilbert;
+		decode_index = decode_hilbert;
+	} else {
+		encode_index = encode_quadkey;
+		decode_index = decode_quadkey;
+	}
+
 	// Wait until here to project the bounding box, so that the behavior is
 	// the same no matter what order the projection and bounding box are
 	// specified in

man/tippecanoe.1

@@ -388,7 +388,7 @@ the single specified tile at that zoom level.
 .PP
 If you know the precision to which you want your data to be represented,
 this table shows the approximate precision corresponding to various
-\fB\fC\-z\fR options\fB\fCif you use the default\fR\-d` detail of 12:
+\fB\fC\-z\fR options if you use the default \fB\fC\-d\fR detail of 12:
 .TS
 allbox;
 cb cb cb
@@ -535,7 +535,7 @@ compensate for the larger marker, or \fB\fC\-Bf\fR\fInumber\fP to allow at most
 .IP \(bu 2
 \fB\fC\-an\fR or \fB\fC\-\-drop\-smallest\-as\-needed\fR: Dynamically drop the smallest features (physically smallest: the shortest lines or the smallest polygons) from each zoom level to keep large tiles under the 500K size limit. This option will not work for point features.
 .IP \(bu 2
-\fB\fC\-aN\fR or \fB\fC\-\-coalesce\-smallest\-as\-needed\fR: Dynamically combine the smallest features (physically smallest: the shortest lines or the smallest polygons) from each zoom level into other nearby features to keep large tiles under the 500K size limit. This option will not work for point features, and will probably not help very much with LineStrings. It is mostly intended for polygons, to maintain the full original area covered by polygons while still reducing the feature count somehow. The attributes of the small polygons are \fInot\fP preserved into the combined features, only their geometry. Furthermore, the polygons that nested polygons are coalesced to may not necessarily be the immediately enclosing features.
+\fB\fC\-aN\fR or \fB\fC\-\-coalesce\-smallest\-as\-needed\fR: Dynamically combine the smallest features (physically smallest: the shortest lines or the smallest polygons) from each zoom level into other nearby features to keep large tiles under the 500K size limit. This option will not work for point features, and will probably not help very much with LineStrings. It is mostly intended for polygons, to maintain the full original area covered by polygons while still reducing the feature count somehow. The attributes of the small polygons are \fInot\fP preserved into the combined features (except through \fB\fC\-\-accumulate\-attribute\fR), only their geometry. Furthermore, the polygons to which nested polygons are coalesced may not necessarily be the immediately enclosing features.
 .IP \(bu 2
 \fB\fC\-aD\fR or \fB\fC\-\-coalesce\-densest\-as\-needed\fR: Dynamically combine the densest features from each zoom level into other nearby features to keep large tiles under the 500K size limit. (Again, mostly useful for polygons.)
 .IP \(bu 2
@@ -585,11 +585,13 @@ the line or polygon within one tile unit of its proper location. You can probabl
 .IP \(bu 2
 \fB\fC\-pi\fR or \fB\fC\-\-preserve\-input\-order\fR: Preserve the original input order of features as the drawing order instead of ordering geographically. (This is implemented as a restoration of the original order at the end, so that dot\-dropping is still geographic, which means it also undoes \fB\fC\-ao\fR).
 .IP \(bu 2
-\fB\fC\-ao\fR or \fB\fC\-\-reorder\fR: Reorder features to put ones with the same attributes in sequence, to try to get them to coalesce. You probably want to use this if you use \fB\fC\-\-coalesce\fR\&.
+\fB\fC\-ac\fR or \fB\fC\-\-coalesce\fR: Coalesce consecutive features that have the same attributes. This can be useful if you have lots of small polygons with identical attributes and you would like to merge them together.
 .IP \(bu 2
-\fB\fC\-ac\fR or \fB\fC\-\-coalesce\fR: Coalesce adjacent features that have the same attributes. This can be useful if you have lots of small polygons with identical attributes and you would like to merge them together.
+\fB\fC\-ao\fR or \fB\fC\-\-reorder\fR: Reorder features to put ones with the same attributes in sequence (instead of ones that are approximately spatially adjacent), to try to get them to coalesce. You probably want to use this if you use \fB\fC\-\-coalesce\fR\&.
 .IP \(bu 2
 \fB\fC\-ar\fR or \fB\fC\-\-reverse\fR: Try reversing the directions of lines to make them coalesce and compress better. You probably don't want to use this.
+.IP \(bu 2
+\fB\fC\-ah\fR or \fB\fC\-\-hilbert\fR: Put features in Hilbert Curve order instead of the usual Z\-Order. This improves the odds that spatially adjacent features will be sequentially adjacent, and should improve density calculations and spatial coalescing. It should be the default eventually.
 .RE
 .SS Adding calculated attributes
 .RS
@@ -601,7 +603,7 @@ the line or polygon within one tile unit of its proper location. You can probabl
 .SS Trying to correct bad source geometry
 .RS
 .IP \(bu 2
-\fB\fC\-aw\fR or \fB\fC\-\-detect\-longitude\-wraparound\fR: Detect when adjacent points within a feature jump to the other side of the world, and try to fix the geometry.
+\fB\fC\-aw\fR or \fB\fC\-\-detect\-longitude\-wraparound\fR: Detect when consecutive points within a feature jump to the other side of the world, and try to fix the geometry.
 .IP \(bu 2
 \fB\fC\-pw\fR or \fB\fC\-\-use\-source\-polygon\-winding\fR: Instead of respecting GeoJSON polygon ring order, use the original polygon winding in the source data to distinguish inner (clockwise) and outer (counterclockwise) polygon rings.
 .IP \(bu 2

options.hpp

@@ -24,6 +24,7 @@
 #define A_CLUSTER_DENSEST_AS_NEEDED ((int) 'C')
 #define A_GENERATE_IDS ((int) 'i')
 #define A_CONVERT_NUMERIC_IDS ((int) 'I')
+#define A_HILBERT ((int) 'h')
 
 #define P_SIMPLIFY ((int) 's')
 #define P_SIMPLIFY_LOW ((int) 'S')

projection.cpp

@@ -5,6 +5,9 @@
 #include <atomic>
 #include "projection.hpp"
 
+unsigned long long (*encode_index)(unsigned int wx, unsigned int wy) = NULL;
+void (*decode_index)(unsigned long long index, unsigned *wx, unsigned *wy) = NULL;
+
 struct projection projections[] = {
 	{"EPSG:4326", lonlat2tile, tile2lonlat, "urn:ogc:def:crs:OGC:1.3:CRS84"},
 	{"EPSG:3857", epsg3857totile, tiletoepsg3857, "urn:ogc:def:crs:EPSG::3857"},
@@ -93,7 +96,60 @@ void tiletoepsg3857(long long ix, long long iy, int zoom, double *ox, double *oy
 	*oy = ((1LL << 32) - 1 - iy - (1LL << 31)) * M_PI * 6378137.0 / (1LL << 31);
 }
 
-unsigned long long encode(unsigned int wx, unsigned int wy) {
+// https://en.wikipedia.org/wiki/Hilbert_curve
+
+void hilbert_rot(unsigned long long n, unsigned *x, unsigned *y, unsigned long long rx, unsigned long long ry) {
+	if (ry == 0) {
+		if (rx == 1) {
+			*x = n - 1 - *x;
+			*y = n - 1 - *y;
+		}
+
+		unsigned t = *x;
+		*x = *y;
+		*y = t;
+	}
+}
+
+unsigned long long hilbert_xy2d(unsigned long long n, unsigned x, unsigned y) {
+	unsigned long long d = 0;
+	unsigned long long rx, ry;
+
+	for (unsigned long long s = n / 2; s > 0; s /= 2) {
+		rx = (x & s) != 0;
+		ry = (y & s) != 0;
+
+		d += s * s * ((3 * rx) ^ ry);
+		hilbert_rot(s, &x, &y, rx, ry);
+	}
+
+	return d;
+}
+
+void hilbert_d2xy(unsigned long long n, unsigned long long d, unsigned *x, unsigned *y) {
+	unsigned long long rx, ry;
+	unsigned long long t = d;
+
+	*x = *y = 0;
+	for (unsigned long long s = 1; s < n; s *= 2) {
+		rx = 1 & (t / 2);
+		ry = 1 & (t ^ rx);
+		hilbert_rot(s, x, y, rx, ry);
+		*x += s * rx;
+		*y += s * ry;
+		t /= 4;
+	}
+}
+
+unsigned long long encode_hilbert(unsigned int wx, unsigned int wy) {
+	return hilbert_xy2d(1LL << 32, wx, wy);
+}
+
+void decode_hilbert(unsigned long long index, unsigned *wx, unsigned *wy) {
+	hilbert_d2xy(1LL << 32, index, wx, wy);
+}
+
+unsigned long long encode_quadkey(unsigned int wx, unsigned int wy) {
 	unsigned long long out = 0;
 
 	int i;
@@ -111,7 +167,7 @@ unsigned long long encode(unsigned int wx, unsigned int wy) {
 static std::atomic<unsigned char> decodex[256];
 static std::atomic<unsigned char> decodey[256];
 
-void decode(unsigned long long index, unsigned *wx, unsigned *wy) {
+void decode_quadkey(unsigned long long index, unsigned *wx, unsigned *wy) {
 	static std::atomic<int> initialized(0);
 	if (!initialized) {
 		for (size_t ix = 0; ix < 256; ix++) {

projection.hpp

@@ -5,8 +5,6 @@ void lonlat2tile(double lon, double lat, int zoom, long long *x, long long *y);
 void epsg3857totile(double ix, double iy, int zoom, long long *x, long long *y);
 void tile2lonlat(long long x, long long y, int zoom, double *lon, double *lat);
 void tiletoepsg3857(long long x, long long y, int zoom, double *ox, double *oy);
-unsigned long long encode(unsigned int wx, unsigned int wy);
-void decode(unsigned long long index, unsigned *wx, unsigned *wy);
 void set_projection_or_exit(const char *optarg);
 
 struct projection {
@@ -19,4 +17,13 @@ struct projection {
 extern struct projection *projection;
 extern struct projection projections[];
 
+extern unsigned long long (*encode_index)(unsigned int wx, unsigned int wy);
+extern void (*decode_index)(unsigned long long index, unsigned *wx, unsigned *wy);
+
+unsigned long long encode_quadkey(unsigned int wx, unsigned int wy);
+void decode_quadkey(unsigned long long index, unsigned *wx, unsigned *wy);
+
+unsigned long long encode_hilbert(unsigned int wx, unsigned int wy);
+void decode_hilbert(unsigned long long index, unsigned *wx, unsigned *wy);
+
 #endif

serial.cpp

@@ -460,7 +460,7 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf) {
 		std::vector<unsigned long long> locs;
 		for (size_t i = 0; i < sf.geometry.size(); i++) {
 			if (sf.geometry[i].op == VT_MOVETO || sf.geometry[i].op == VT_LINETO) {
-				locs.push_back(encode(sf.geometry[i].x << geometry_scale, sf.geometry[i].y << geometry_scale));
+				locs.push_back(encode_index(sf.geometry[i].x << geometry_scale, sf.geometry[i].y << geometry_scale));
 			}
 		}
 		std::sort(locs.begin(), locs.end());
@@ -552,7 +552,7 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf) {
 	// and then mask to bring it back into the addressable area
 	long long midx = (sf.bbox[0] / 2 + sf.bbox[2] / 2) & ((1LL << 32) - 1);
 	long long midy = (sf.bbox[1] / 2 + sf.bbox[3] / 2) & ((1LL << 32) - 1);
-	bbox_index = encode(midx, midy);
+	bbox_index = encode_index(midx, midy);
 
 	if (additional[A_DROP_DENSEST_AS_NEEDED] || additional[A_COALESCE_DENSEST_AS_NEEDED] || additional[A_CLUSTER_DENSEST_AS_NEEDED] || additional[A_CALCULATE_FEATURE_DENSITY] || additional[A_INCREASE_GAMMA_AS_NEEDED] || sst->uses_gamma || cluster_distance != 0) {
 		sf.index = bbox_index;
@@ -633,7 +633,6 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf) {
 				}
 			}
 		}
-
 	}
 
 	if (!sst->filters) {