Optimize Hilbert tile ID <-> XYZ conversion (#208)

pmtiles/tile_id.go

@@ -1,5 +1,9 @@
 package pmtiles
 
+import (
+	"math/bits"
+)
+
 func rotate(n uint64, x *uint64, y *uint64, rx uint64, ry uint64) {
 	if ry == 0 {
 		if rx == 1 {
@@ -10,82 +14,42 @@ func rotate(n uint64, x *uint64, y *uint64, rx uint64, ry uint64) {
 	}
 }
 
-func tOnLevel(z uint8, pos uint64) (uint8, uint32, uint32) {
-	var n uint64 = 1 << z
-	rx, ry, t := pos, pos, pos
-	var tx uint64
-	var ty uint64
-	var s uint64
-	for s = 1; s < n; s *= 2 {
-		rx = 1 & (t / 2)
-		ry = 1 & (t ^ rx)
-		rotate(s, &tx, &ty, rx, ry)
-		tx += s * rx
-		ty += s * ry
-		t /= 4
-	}
-	return uint8(z), uint32(tx), uint32(ty)
-}
-
 // ZxyToID converts (Z,X,Y) tile coordinates to a Hilbert TileID.
 func ZxyToID(z uint8, x uint32, y uint32) uint64 {
-	var acc uint64
-	var tz uint8
-	for ; tz < z; tz++ {
-		acc += (0x1 << tz) * (0x1 << tz)
-	}
-	var n uint64 = 1 << z
-	var rx uint64
-	var ry uint64
-	var d uint64
-	tx := uint64(x)
-	ty := uint64(y)
-	for s := n / 2; s > 0; s /= 2 {
-		if tx&s > 0 {
-			rx = 1
-		} else {
-			rx = 0
-		}
-		if ty&s > 0 {
-			ry = 1
-		} else {
-			ry = 0
-		}
-		d += s * s * ((3 * rx) ^ ry)
+	var acc uint64 = ((1 << (z * 2)) - 1) / 3
+	var tx, ty uint64 = uint64(x), uint64(y)
+	for a := int32(z - 1); a >= 0; a-- {
+		var rx uint64 = (tx >> a) & 1
+		var ry uint64 = (ty >> a) & 1
+		var s uint64 = (1 << a)
 		rotate(s, &tx, &ty, rx, ry)
+		acc += s * s * ((3 * rx) ^ ry)
 	}
-	return acc + d
+	return acc
 }
 
 // IDToZxy converts a Hilbert TileID to (Z,X,Y) tile coordinates.
 func IDToZxy(i uint64) (uint8, uint32, uint32) {
-	var acc uint64
-	var z uint8
-	for {
-		var numTiles uint64
-		numTiles = (1 << z) * (1 << z)
-		if acc+numTiles > i {
-			return tOnLevel(z, i-acc)
-		}
-		acc += numTiles
-		z++
+	var z uint8 = uint8((64 - bits.LeadingZeros64(3*i+1) - 1) / 2)
+	var acc uint64 = (1<<(z*2) - 1) / 3
+	var pos uint64 = i - acc
+	var tx, ty uint64 = 0, 0
+	for a := uint8(0); a < z; a++ {
+		var rx uint64 = (pos / 2) & 1
+		var ry uint64 = (pos ^ rx) & 1
+		var s uint64 = 1 << a
+		rotate(s, &tx, &ty, rx, ry)
+		tx += s * rx
+		ty += s * ry
+		pos /= 4
 	}
+	return z, uint32(tx), uint32(ty)
 }
 
 // ParentID efficiently finds a parent Hilbert TileID without converting to (Z,X,Y).
 func ParentID(i uint64) uint64 {
-	var acc uint64
-	var lastAcc uint64
-	var z uint8
-	for {
-		var numTiles uint64
-		numTiles = (1 << z) * (1 << z)
-		if acc+numTiles > i {
-			return lastAcc + (i-acc)/4
-		}
-		lastAcc = acc
-		acc += numTiles
-		z++
-	}
-
+	var z uint8 = uint8((64 - bits.LeadingZeros64(3*i+1) - 1) / 2)
+	var acc uint64 = (1<<(z*2) - 1) / 3
+	var parentAcc uint64 = (1<<((z-1)*2) - 1) / 3
+	return parentAcc + (i-acc)/4
 }

pmtiles/tile_id_test.go

@@ -1,8 +1,9 @@
 package pmtiles
 
 import (
-	"github.com/stretchr/testify/assert"
 	"testing"
+
+	"github.com/stretchr/testify/assert"
 )
 
 func TestZxyToId(t *testing.T) {
@@ -49,8 +50,7 @@ func TestManyTileIds(t *testing.T) {
 func TestExtremes(t *testing.T) {
 	var tz uint8
 	for tz = 0; tz < 32; tz++ {
-		var dim uint32
-		dim = (1 << tz) - 1
+		var dim uint32 = (1 << tz) - 1
 		z, x, y := IDToZxy(ZxyToID(tz, 0, 0))
 		assert.Equal(t, tz, z)
 		assert.Equal(t, uint32(0), x)
@@ -92,4 +92,9 @@ func TestParent(t *testing.T) {
 	assert.Equal(t, ZxyToID(1, 1, 1), ParentID(ZxyToID(2, 2, 3)))
 	assert.Equal(t, ZxyToID(1, 1, 1), ParentID(ZxyToID(2, 3, 2)))
 	assert.Equal(t, ZxyToID(1, 1, 1), ParentID(ZxyToID(2, 3, 3)))
+
+	assert.Equal(t, ZxyToID(18, 500, 1), ParentID(ZxyToID(19, 1000, 3)))
+	assert.Equal(t, ZxyToID(18, 500, 2), ParentID(ZxyToID(19, 1000, 4)))
+	assert.Equal(t, ZxyToID(18, 1, 500), ParentID(ZxyToID(19, 3, 1000)))
+	assert.Equal(t, ZxyToID(18, 2, 500), ParentID(ZxyToID(19, 4, 1000)))
 }