Simplify and optimize Hilbert tile ID <-> XYZ conversion (#527)

* Optimize Hilbert tile ID <-> XYZ conversion

* biome

* tile_id: more improvements

* tile_id(cpp): uint32

cpp/pmtiles.hpp

@@ -321,35 +321,18 @@ uint64_t decode_varint(const char **data, const char *end) {
 	return decode_varint_impl(data, end);
 }
 
-void rotate(int64_t n, int64_t &x, int64_t &y, int64_t rx, int64_t ry) {
+void rotate(int64_t n, uint32_t &x, uint32_t &y, uint32_t rx, uint32_t ry) {
 	if (ry == 0) {
-		if (rx == 1) {
+		if (rx != 0) {
 			x = n - 1 - x;
 			y = n - 1 - y;
 		}
-		int64_t t = x;
+		uint32_t t = x;
 		x = y;
 		y = t;
 	}
 }
 
-zxy t_on_level(uint8_t z, uint64_t pos) {
-	int64_t n = 1LL << z;
-	int64_t rx, ry, s, t = pos;
-	int64_t tx = 0;
-	int64_t ty = 0;
-
-	for (s = 1; s < n; s *= 2) {
-		rx = 1LL & (t / 2);
-		ry = 1LL & (t ^ rx);
-		rotate(s, tx, ty, rx, ry);
-		tx += s * rx;
-		ty += s * ry;
-		t /= 4;
-	}
-	return zxy(z, static_cast<int>(tx), static_cast<int>(ty));
-}
-
 int write_varint(std::back_insert_iterator<std::string> data, uint64_t value) {
 	int n = 1;
 
@@ -405,17 +388,32 @@ entryv3 find_tile(const std::vector<entryv3> &entries, uint64_t tile_id) {
 
 }  // end anonymous namespace
 
+// Note: std::bit_width is available in C++20 and later.
+static uint8_t bit_width(uint64_t n) {
+    uint8_t count = 0;
+    while (n > 0) { count++; n >>= 1; }
+    return count;
+}
+
 inline zxy tileid_to_zxy(uint64_t tileid) {
-	uint64_t acc = 0;
+	if (tileid > 6148914691236517204) {
-	for (uint8_t t_z = 0; t_z < 32; t_z++) {
-		uint64_t num_tiles = (1LL << t_z) * (1LL << t_z);
-		if (acc + num_tiles > tileid) {
-			return t_on_level(t_z, tileid - acc);
-		}
-		acc += num_tiles;
-	}
 		throw std::overflow_error("tile zoom exceeds 64-bit limit");
 	}
+	uint8_t z = (bit_width(3 * tileid + 1) - 1) / 2;
+	uint64_t acc = ((1L << (z * 2)) - 1) / 3;
+	uint64_t pos = tileid - acc;
+	uint32_t x = 0, y = 0;
+	for (uint8_t a = 0; a < z; a++) {
+        uint64_t s = 1 << a;
+        uint32_t rx = s & (pos / 2);
+        uint32_t ry = s & (pos ^ rx);
+		rotate(s, x, y, rx, ry);
+        pos >>= 1;
+        x += rx;
+        y += ry;
+	}
+	return zxy(z, static_cast<int>(x), static_cast<int>(y));
+}
 
 inline uint64_t zxy_to_tileid(uint8_t z, uint32_t x, uint32_t y) {
 	if (z > 31) {
@@ -424,19 +422,17 @@ inline uint64_t zxy_to_tileid(uint8_t z, uint32_t x, uint32_t y) {
 	if (x > (1U << z) - 1U || y > (1U << z) - 1U) {
 		throw std::overflow_error("tile x/y outside zoom level bounds");
 	}
-	uint64_t acc = 0;
+	uint64_t acc = ((1LL << (z * 2U)) - 1) / 3;
-	for (uint8_t t_z = 0; t_z < z; t_z++) acc += (1LL << t_z) * (1LL << t_z);
+	uint32_t tx = x, ty = y;
-	int64_t n = 1LL << z;
+	int a = z - 1;
-	int64_t rx, ry, s, d = 0;
+	for (uint32_t s = 1LL << a; s > 0; s >>= 1) {
-	int64_t tx = x;
+		uint32_t rx = s & tx;
-	int64_t ty = y;
+		uint32_t ry = s & ty;
-	for (s = n / 2; s > 0; s /= 2) {
-		rx = (tx & s) > 0;
-		ry = (ty & s) > 0;
-		d += s * s * ((3LL * rx) ^ ry);
 		rotate(s, tx, ty, rx, ry);
+		acc += ((3LL * rx) ^ ry) << a;
+		a--;
 	}
-	return acc + d;
+	return acc;
 }
 
 // returns an uncompressed byte buffer

cpp/test.cpp

@@ -72,7 +72,6 @@ MU_TEST(test_roundtrip) {
 	} catch (const std::runtime_error &e) {
 		caught = true;
 	}
-
 	mu_check(caught);
 
 	caught = false;
@@ -81,7 +80,6 @@ MU_TEST(test_roundtrip) {
 	} catch (const std::runtime_error &e) {
 		caught = true;
 	}
-
 	mu_check(caught);
 
 	caught = false;
@@ -90,7 +88,6 @@ MU_TEST(test_roundtrip) {
 	} catch (const std::runtime_error &e) {
 		caught = true;
 	}
-
 	mu_check(caught);
 
 	caught = false;

js/src/index.ts

@@ -55,44 +55,22 @@ export function readVarint(p: BufferPosition): number {
   return readVarintRemainder(val, p);
 }
 
-function rotate(n: number, xy: number[], rx: number, ry: number): void {
+function rotate(
+  n: number,
+  x: number,
+  y: number,
+  rx: number,
+  ry: number
+): [number, number] {
   if (ry === 0) {
-    if (rx === 1) {
+    if (rx !== 0) {
-      xy[0] = n - 1 - xy[0];
+      return [n - 1 - y, n - 1 - x];
-      xy[1] = n - 1 - xy[1];
     }
-    const t = xy[0];
+    return [y, x];
-    xy[0] = xy[1];
-    xy[1] = t;
   }
+  return [x, y];
 }
 
-function idOnLevel(z: number, pos: number): [number, number, number] {
-  const n = 2 ** z;
-  let rx = pos;
-  let ry = pos;
-  let t = pos;
-  const xy = [0, 0];
-  let s = 1;
-  while (s < n) {
-    rx = 1 & (t / 2);
-    ry = 1 & (t ^ rx);
-    rotate(s, xy, rx, ry);
-    xy[0] += s * rx;
-    xy[1] += s * ry;
-    t = t / 4;
-    s *= 2;
-  }
-  return [z, xy[0], xy[1]];
-}
-
-const tzValues: number[] = [
-  0, 1, 5, 21, 85, 341, 1365, 5461, 21845, 87381, 349525, 1398101, 5592405,
-  22369621, 89478485, 357913941, 1431655765, 5726623061, 22906492245,
-  91625968981, 366503875925, 1466015503701, 5864062014805, 23456248059221,
-  93824992236885, 375299968947541, 1501199875790165,
-];
-
 /**
  * Convert Z,X,Y to a Hilbert TileID.
  */
@@ -100,43 +78,52 @@ export function zxyToTileId(z: number, x: number, y: number): number {
   if (z > 26) {
     throw new Error("Tile zoom level exceeds max safe number limit (26)");
   }
-  if (x > 2 ** z - 1 || y > 2 ** z - 1) {
+  if (x >= 1 << z || y >= 1 << z) {
     throw new Error("tile x/y outside zoom level bounds");
   }
-
+  let acc = ((1 << z) * (1 << z) - 1) / 3;
-  const acc = tzValues[z];
+  let a = z - 1;
-  const n = 2 ** z;
+  let [tx, ty] = [x, y];
-  let rx = 0;
+  for (let s = 1 << a; s > 0; s >>= 1) {
-  let ry = 0;
+    const rx = tx & s;
-  let d = 0;
+    const ry = ty & s;
-  const xy = [x, y];
+    acc += ((3 * rx) ^ ry) * (1 << a);
-  let s = n / 2;
+    [tx, ty] = rotate(s, tx, ty, rx, ry);
-  while (s > 0) {
+    a--;
-    rx = (xy[0] & s) > 0 ? 1 : 0;
-    ry = (xy[1] & s) > 0 ? 1 : 0;
-    d += s * s * ((3 * rx) ^ ry);
-    rotate(s, xy, rx, ry);
-    s = s / 2;
   }
-  return acc + d;
+  return acc;
+}
+
+function tileIdToZ(i: number): number {
+  const c = 3 * i + 1;
+  if (c < 0x100000000) {
+    return 31 - Math.clz32(c);
+  }
+  return 63 - Math.clz32(c / 0x100000000);
 }
 
 /**
  * Convert a Hilbert TileID to Z,X,Y.
  */
 export function tileIdToZxy(i: number): [number, number, number] {
-  let acc = 0;
+  const z = tileIdToZ(i) >> 1;
-  const z = 0;
+  if (z > 26)
-
-  for (let z = 0; z < 27; z++) {
-    const numTiles = (0x1 << z) * (0x1 << z);
-    if (acc + numTiles > i) {
-      return idOnLevel(z, i - acc);
-    }
-    acc += numTiles;
-  }
-
     throw new Error("Tile zoom level exceeds max safe number limit (26)");
+  const acc = ((1 << z) * (1 << z) - 1) / 3;
+
+  let t = i - acc;
+  let x = 0;
+  let y = 0;
+  const n = 1 << z;
+  for (let s = 1; s < n; s <<= 1) {
+    const rx = s & (t / 2);
+    const ry = s & (t ^ rx);
+    [x, y] = rotate(s, x, y, rx, ry);
+    t = t / 2;
+    x += rx;
+    y += ry;
+  }
+  return [z, x, y];
 }
 
 /**

python/pmtiles/tile.py

@@ -1,6 +1,6 @@
-from enum import Enum
-import io
 import gzip
+import io
+from enum import Enum
 
 
 class Entry:
@@ -16,28 +16,13 @@ class Entry:
         return f"id={self.tile_id} offset={self.offset} length={self.length} runlength={self.run_length}"
 
 
-def rotate(n, xy, rx, ry):
+def rotate(n, x, y, rx, ry):
     if ry == 0:
-        if rx == 1:
+        if rx != 0:
-            xy[0] = n - 1 - xy[0]
+            x = n - 1 - x
-            xy[1] = n - 1 - xy[1]
+            y = n - 1 - y
-        xy[0], xy[1] = xy[1], xy[0]
+        x, y = y, x
-
+    return x, y
-
-def t_on_level(z, pos):
-    n = 1 << z
-    rx, ry, t = pos, pos, pos
-    xy = [0, 0]
-    s = 1
-    while s < n:
-        rx = 1 & (t // 2)
-        ry = 1 & (t ^ rx)
-        rotate(s, xy, rx, ry)
-        xy[0] += s * rx
-        xy[1] += s * ry
-        t //= 4
-        s *= 2
-    return z, xy[0], xy[1]
 
 
 def zxy_to_tileid(z, x, y):
@@ -45,41 +30,38 @@ def zxy_to_tileid(z, x, y):
         raise OverflowError("tile zoom exceeds 64-bit limit")
     if x > (1 << z) - 1 or y > (1 << z) - 1:
         raise ValueError("tile x/y outside zoom level bounds")
-    acc = 0
+
-    tz = 0
+    acc = ((1 << (z * 2)) - 1) // 3
-    while tz < z:
+    a = z - 1
-        acc += (0x1 << tz) * (0x1 << tz)
+    while a >= 0:
-        tz += 1
+        s = 1 << a
-    n = 1 << z
+        rx = s & x
-    rx = 0
+        ry = s & y
-    ry = 0
+        acc += ((3 * rx) ^ ry) << a
-    d = 0
+        (x, y) = rotate(s, x, y, rx, ry)
-    xy = [x, y]
+        a -= 1
-    s = n // 2
+    return acc
-    while s > 0:
-        if (xy[0] & s) > 0:
-            rx = 1
-        else:
-            rx = 0
-        if (xy[1] & s) > 0:
-            ry = 1
-        else:
-            ry = 0
-        d += s * s * ((3 * rx) ^ ry)
-        rotate(s, xy, rx, ry)
-        s //= 2
-    return acc + d
 
 
 def tileid_to_zxy(tile_id):
-    num_tiles = 0
+    z = ((3 * tile_id + 1).bit_length() - 1) // 2
-    acc = 0
+    if z >= 32:
-    for z in range(0,32):
-        num_tiles = (1 << z) * (1 << z)
-        if acc + num_tiles > tile_id:
-            return t_on_level(z, tile_id - acc)
-        acc += num_tiles
         raise OverflowError("tile zoom exceeds 64-bit limit")
+    acc = ((1 << (z * 2)) - 1) // 3
+    pos = tile_id - acc
+    x = 0
+    y = 0
+    s = 1
+    n = 1 << z
+    while s < n:
+        rx = (pos // 2) & s
+        ry = (pos ^ rx) & s
+        (x, y) = rotate(s, x, y, rx, ry)
+        x += rx
+        y += ry
+        pos >>= 1
+        s <<= 1
+    return (z, x, y)
 
 
 def find_tile(entries, tile_id):
@@ -195,12 +177,15 @@ def serialize_directory(entries):
 
     return gzip.compress(b_io.getvalue())
 
+
 class SpecVersionUnsupported(Exception):
     pass
 
+
 class MagicNumberNotFound(Exception):
     pass
 
+
 def deserialize_header(buf):
     if buf[0:7].decode() != "PMTiles":
         raise MagicNumberNotFound()