nx_geo_s2.nx source

1// nx_geo_s2.nx -- LIB: GEO-005 S2-CLASS CELL INDEX via the HILBERT space-filling curve. The real 2// exceed over GEO-002's geohash (which interleaves lon/lat bits = the Z-ORDER / Morton curve). 3// 4// THE EXCEED ANGLE (measured, not asserted): a Hilbert curve is CONTINUOUS -- consecutive 1D cell 5// indices are ALWAYS 4-adjacent in 2D (Manhattan step exactly 1) -- so a 1D range scan never "tears" 6// across distant 2D regions. The Z-order curve that geohash uses TEARS: its consecutive step jumps by 7// up to the grid width N at every major quadrant boundary. Google's S2 chose Hilbert for exactly this 8// locality; here it is sovereign + integer. (cell index quality is what makes range queries / tile 9// coherence good; this is the measurable weakness of plain geohash.) 10// 11// Pure integer bit arithmetic (bit-tests via /,% to avoid the geohash const-index landmine), 12// deterministic. (x,y) are integer cell coords in [0, 2^k); a real lon/lat first quantizes to that 13// grid (GEO-002 already does the microdeg bisection). license_tier: ORIGINAL 14import "nx_syscalls.nx" 15 16func geo_s2_pow2(k: i64) -> i64 { var n: i64 = 1; var t: i64 = 0; while t < k { n = n * 2; t = t + 1 } return n } 17 18// (x,y) -> Hilbert distance d on the order-k curve (grid n=2^k). Standard rotate-and-fold algorithm. 19func geo_hilbert_xy2d(k: i64, x: i64, y: i64) -> i64 { 20 let n: i64 = geo_s2_pow2(k) 21 var d: i64 = 0 22 var xx: i64 = x 23 var yy: i64 = y 24 var s: i64 = n / 2 25 while s > 0 { 26 var rx: i64 = 0 27 if ((xx / s) % 2) == 1 { rx = 1 } 28 var ry: i64 = 0 29 if ((yy / s) % 2) == 1 { ry = 1 } 30 // q = (3*rx) XOR ry, computed without a bitwise xor. 31 var q: i64 = 3 * rx 32 if ry == 1 { if rx == 0 { q = 1 } else { q = 2 } } 33 d = d + s * s * q 34 // rotate the quadrant 35 if ry == 0 { 36 if rx == 1 { 37 xx = n - 1 - xx 38 yy = n - 1 - yy 39 } 40 let tmp: i64 = xx 41 xx = yy 42 yy = tmp 43 } 44 s = s / 2 45 } 46 return d 47} 48 49// inverse: Hilbert distance d -> (x,y), written to out2[0]=x, out2[1]=y. 50func geo_hilbert_d2xy(k: i64, d: i64, out2: *i64) -> i64 { 51 let n: i64 = geo_s2_pow2(k) 52 var rx: i64 = 0 53 var ry: i64 = 0 54 var s: i64 = 1 55 var t: i64 = d 56 var xx: i64 = 0 57 var yy: i64 = 0 58 while s < n { 59 rx = (t / 2) % 2 60 ry = ((t % 2) + rx) % 2 // = (t XOR rx) AND 1 61 if ry == 0 { 62 if rx == 1 { 63 xx = s - 1 - xx 64 yy = s - 1 - yy 65 } 66 let tmp: i64 = xx 67 xx = yy 68 yy = tmp 69 } 70 xx = xx + s * rx 71 yy = yy + s * ry 72 t = t / 4 73 s = s * 2 74 } 75 out2[0] = xx 76 out2[1] = yy 77 return 0 78} 79 80// (x,y) -> Z-order / Morton distance (the geohash-family curve) on the order-k grid, for comparison. 81func geo_zorder_xy2d(k: i64, x: i64, y: i64) -> i64 { 82 var d: i64 = 0 83 var i: i64 = 0 84 var bitval: i64 = 1 // 4^i = position weight of bit i of x; y's bit sits one above (bitval*2) 85 var xs: i64 = x 86 var ys: i64 = y 87 while i < k { 88 d = d + (xs % 2) * bitval + (ys % 2) * bitval * 2 89 xs = xs / 2 90 ys = ys / 2 91 bitval = bitval * 4 92 i = i + 1 93 } 94 return d 95}