From 9b69ab4b64910359cbb0576ad5065ee271c93154 Mon Sep 17 00:00:00 2001 From: Taku Fukada Date: Fri, 21 Feb 2025 16:37:49 +0900 Subject: [PATCH] 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 | 74 ++++++++++++++--------------- cpp/test.cpp | 3 -- js/src/index.ts | 105 ++++++++++++++++++----------------------- python/pmtiles/tile.py | 97 ++++++++++++++++--------------------- 4 files changed, 122 insertions(+), 157 deletions(-) diff --git a/cpp/pmtiles.hpp b/cpp/pmtiles.hpp index e8dd71f..2819f33 100644 --- a/cpp/pmtiles.hpp +++ b/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(tx), static_cast(ty)); -} - int write_varint(std::back_insert_iterator data, uint64_t value) { int n = 1; @@ -405,16 +388,31 @@ entryv3 find_tile(const std::vector &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; - 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; + if (tileid > 6148914691236517204) { + throw std::overflow_error("tile zoom exceeds 64-bit limit"); } - 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(x), static_cast(y)); } inline uint64_t zxy_to_tileid(uint8_t z, uint32_t x, uint32_t y) { @@ -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; - for (uint8_t t_z = 0; t_z < z; t_z++) acc += (1LL << t_z) * (1LL << t_z); - int64_t n = 1LL << z; - int64_t rx, ry, s, d = 0; - int64_t tx = x; - int64_t ty = y; - for (s = n / 2; s > 0; s /= 2) { - rx = (tx & s) > 0; - ry = (ty & s) > 0; - d += s * s * ((3LL * rx) ^ ry); + uint64_t acc = ((1LL << (z * 2U)) - 1) / 3; + uint32_t tx = x, ty = y; + int a = z - 1; + for (uint32_t s = 1LL << a; s > 0; s >>= 1) { + uint32_t rx = s & tx; + uint32_t ry = s & ty; rotate(s, tx, ty, rx, ry); + acc += ((3LL * rx) ^ ry) << a; + a--; } - return acc + d; + return acc; } // returns an uncompressed byte buffer diff --git a/cpp/test.cpp b/cpp/test.cpp index 0d9088c..c0d3b7c 100644 --- a/cpp/test.cpp +++ b/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; diff --git a/js/src/index.ts b/js/src/index.ts index 3c117af..437ed2f 100644 --- a/js/src/index.ts +++ b/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) { - xy[0] = n - 1 - xy[0]; - xy[1] = n - 1 - xy[1]; + if (rx !== 0) { + return [n - 1 - y, n - 1 - x]; } - const t = xy[0]; - xy[0] = xy[1]; - xy[1] = t; + return [y, x]; } + 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"); } - - const acc = tzValues[z]; - const n = 2 ** z; - let rx = 0; - let ry = 0; - let d = 0; - const xy = [x, y]; - let s = n / 2; - while (s > 0) { - 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; + let acc = ((1 << z) * (1 << z) - 1) / 3; + let a = z - 1; + let [tx, ty] = [x, y]; + for (let s = 1 << a; s > 0; s >>= 1) { + const rx = tx & s; + const ry = ty & s; + acc += ((3 * rx) ^ ry) * (1 << a); + [tx, ty] = rotate(s, tx, ty, rx, ry); + a--; } - 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 = 0; + const z = tileIdToZ(i) >> 1; + if (z > 26) + throw new Error("Tile zoom level exceeds max safe number limit (26)"); + const acc = ((1 << z) * (1 << z) - 1) / 3; - for (let z = 0; z < 27; z++) { - const numTiles = (0x1 << z) * (0x1 << z); - if (acc + numTiles > i) { - return idOnLevel(z, i - acc); - } - acc += numTiles; + 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; } - - throw new Error("Tile zoom level exceeds max safe number limit (26)"); + return [z, x, y]; } /** diff --git a/python/pmtiles/tile.py b/python/pmtiles/tile.py index 0143256..6e0429f 100644 --- a/python/pmtiles/tile.py +++ b/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: - xy[0] = n - 1 - xy[0] - xy[1] = n - 1 - xy[1] - xy[0], xy[1] = xy[1], xy[0] - - -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] + if rx != 0: + x = n - 1 - x + y = n - 1 - y + x, y = y, x + return x, y 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 - while tz < z: - acc += (0x1 << tz) * (0x1 << tz) - tz += 1 - n = 1 << z - rx = 0 - ry = 0 - d = 0 - xy = [x, y] - s = n // 2 - 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 + + acc = ((1 << (z * 2)) - 1) // 3 + a = z - 1 + while a >= 0: + s = 1 << a + rx = s & x + ry = s & y + acc += ((3 * rx) ^ ry) << a + (x, y) = rotate(s, x, y, rx, ry) + a -= 1 + return acc def tileid_to_zxy(tile_id): - num_tiles = 0 - acc = 0 - 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") + z = ((3 * tile_id + 1).bit_length() - 1) // 2 + if z >= 32: + 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() @@ -282,7 +267,7 @@ def serialize_header(h): write_int32(max_lon_e7) max_lat_e7 = h.get("max_lat_e7", int(90 * 10000000)) write_int32(max_lat_e7) - write_uint8(h.get("center_zoom",h["min_zoom"])) + write_uint8(h.get("center_zoom", h["min_zoom"])) write_int32(h.get("center_lon_e7", round((min_lon_e7 + max_lon_e7) / 2))) write_int32(h.get("center_lat_e7", round((min_lat_e7 + max_lat_e7) / 2)))