1 // Copyright 2014 Google Inc. All rights reserved. 2 // 3 // Licensed under the Apache License, Version 2.0 (the "License"); 4 // you may not use this file except in compliance with the License. 5 // You may obtain a copy of the License at 6 // 7 // http://www.apache.org/licenses/LICENSE-2.0 8 // 9 // Unless required by applicable law or agreed to in writing, software 10 // distributed under the License is distributed on an "AS IS" BASIS, 11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 // See the License for the specific language governing permissions and 13 // limitations under the License. 14 15 package s2 16 17 import ( 18 "bytes" 19 "fmt" 20 "io" 21 "math" 22 "sort" 23 "strconv" 24 "strings" 25 26 "github.com/golang/geo/r1" 27 "github.com/golang/geo/r2" 28 "github.com/golang/geo/r3" 29 "github.com/golang/geo/s1" 30 ) 31 32 // CellID uniquely identifies a cell in the S2 cell decomposition. 33 // The most significant 3 bits encode the face number (0-5). The 34 // remaining 61 bits encode the position of the center of this cell 35 // along the Hilbert curve on that face. The zero value and the value 36 // (1<<64)-1 are invalid cell IDs. The first compares less than any 37 // valid cell ID, the second as greater than any valid cell ID. 38 // 39 // Sequentially increasing cell IDs follow a continuous space-filling curve 40 // over the entire sphere. They have the following properties: 41 // 42 // - The ID of a cell at level k consists of a 3-bit face number followed 43 // by k bit pairs that recursively select one of the four children of 44 // each cell. The next bit is always 1, and all other bits are 0. 45 // Therefore, the level of a cell is determined by the position of its 46 // lowest-numbered bit that is turned on (for a cell at level k, this 47 // position is 2 * (MaxLevel - k)). 48 // 49 // - The ID of a parent cell is at the midpoint of the range of IDs spanned 50 // by its children (or by its descendants at any level). 51 // 52 // Leaf cells are often used to represent points on the unit sphere, and 53 // this type provides methods for converting directly between these two 54 // representations. For cells that represent 2D regions rather than 55 // discrete point, it is better to use Cells. 56 type CellID uint64 57 58 // SentinelCellID is an invalid cell ID guaranteed to be larger than any 59 // valid cell ID. It is used primarily by ShapeIndex. The value is also used 60 // by some S2 types when encoding data. 61 // Note that the sentinel's RangeMin == RangeMax == itself. 62 const SentinelCellID = CellID(^uint64(0)) 63 64 // sortCellIDs sorts the slice of CellIDs in place. 65 func sortCellIDs(ci []CellID) { 66 sort.Sort(cellIDs(ci)) 67 } 68 69 // cellIDs implements the Sort interface for slices of CellIDs. 70 type cellIDs []CellID 71 72 func (c cellIDs) Len() int { return len(c) } 73 func (c cellIDs) Swap(i, j int) { c[i], c[j] = c[j], c[i] } 74 func (c cellIDs) Less(i, j int) bool { return c[i] < c[j] } 75 76 const ( 77 // FaceBits is the number of bits used to encode the face number. 78 FaceBits = 3 79 // NumFaces is the number of faces. 80 NumFaces = 6 81 82 // MaxLevel is the number of levels needed to specify a leaf cell. 83 MaxLevel = 30 84 85 // PosBits is the total number of position bits. The extra bit (61 rather 86 // than 60) lets us encode each cell as its Hilbert curve position at the 87 // cell center (which is halfway along the portion of the Hilbert curve that 88 // fills that cell). 89 PosBits = 2*MaxLevel + 1 90 91 // MaxSize is the maximum index of a valid leaf cell plus one. The range of 92 // valid leaf cell indices is [0..MaxSize-1]. 93 MaxSize = 1 << MaxLevel 94 95 wrapOffset = uint64(NumFaces) << PosBits 96 ) 97 98 // CellIDFromFacePosLevel returns a cell given its face in the range 99 // [0,5], the 61-bit Hilbert curve position pos within that face, and 100 // the level in the range [0,MaxLevel]. The position in the cell ID 101 // will be truncated to correspond to the Hilbert curve position at 102 // the center of the returned cell. 103 func CellIDFromFacePosLevel(face int, pos uint64, level int) CellID { 104 return CellID(uint64(face)<<PosBits + pos | 1).Parent(level) 105 } 106 107 // CellIDFromFace returns the cell corresponding to a given S2 cube face. 108 func CellIDFromFace(face int) CellID { 109 return CellID((uint64(face) << PosBits) + lsbForLevel(0)) 110 } 111 112 // CellIDFromLatLng returns the leaf cell containing ll. 113 func CellIDFromLatLng(ll LatLng) CellID { 114 return cellIDFromPoint(PointFromLatLng(ll)) 115 } 116 117 // CellIDFromToken returns a cell given a hex-encoded string of its uint64 ID. 118 func CellIDFromToken(s string) CellID { 119 if len(s) > 16 { 120 return CellID(0) 121 } 122 n, err := strconv.ParseUint(s, 16, 64) 123 if err != nil { 124 return CellID(0) 125 } 126 // Equivalent to right-padding string with zeros to 16 characters. 127 if len(s) < 16 { 128 n = n << (4 * uint(16-len(s))) 129 } 130 return CellID(n) 131 } 132 133 // ToToken returns a hex-encoded string of the uint64 cell id, with leading 134 // zeros included but trailing zeros stripped. 135 func (ci CellID) ToToken() string { 136 s := strings.TrimRight(fmt.Sprintf("%016x", uint64(ci)), "0") 137 if len(s) == 0 { 138 return "X" 139 } 140 return s 141 } 142 143 // IsValid reports whether ci represents a valid cell. 144 func (ci CellID) IsValid() bool { 145 return ci.Face() < NumFaces && (ci.lsb()&0x1555555555555555 != 0) 146 } 147 148 // Face returns the cube face for this cell ID, in the range [0,5]. 149 func (ci CellID) Face() int { return int(uint64(ci) >> PosBits) } 150 151 // Pos returns the position along the Hilbert curve of this cell ID, in the range [0,2^PosBits-1]. 152 func (ci CellID) Pos() uint64 { return uint64(ci) & (^uint64(0) >> FaceBits) } 153 154 // Level returns the subdivision level of this cell ID, in the range [0, MaxLevel]. 155 func (ci CellID) Level() int { 156 return MaxLevel - findLSBSetNonZero64(uint64(ci))>>1 157 } 158 159 // IsLeaf returns whether this cell ID is at the deepest level; 160 // that is, the level at which the cells are smallest. 161 func (ci CellID) IsLeaf() bool { return uint64(ci)&1 != 0 } 162 163 // ChildPosition returns the child position (0..3) of this cell's 164 // ancestor at the given level, relative to its parent. The argument 165 // should be in the range 1..kMaxLevel. For example, 166 // ChildPosition(1) returns the position of this cell's level-1 167 // ancestor within its top-level face cell. 168 func (ci CellID) ChildPosition(level int) int { 169 return int(uint64(ci)>>uint64(2*(MaxLevel-level)+1)) & 3 170 } 171 172 // lsbForLevel returns the lowest-numbered bit that is on for cells at the given level. 173 func lsbForLevel(level int) uint64 { return 1 << uint64(2*(MaxLevel-level)) } 174 175 // Parent returns the cell at the given level, which must be no greater than the current level. 176 func (ci CellID) Parent(level int) CellID { 177 lsb := lsbForLevel(level) 178 return CellID((uint64(ci) & -lsb) | lsb) 179 } 180 181 // immediateParent is cheaper than Parent, but assumes !ci.isFace(). 182 func (ci CellID) immediateParent() CellID { 183 nlsb := CellID(ci.lsb() << 2) 184 return (ci & -nlsb) | nlsb 185 } 186 187 // isFace returns whether this is a top-level (face) cell. 188 func (ci CellID) isFace() bool { return uint64(ci)&(lsbForLevel(0)-1) == 0 } 189 190 // lsb returns the least significant bit that is set. 191 func (ci CellID) lsb() uint64 { return uint64(ci) & -uint64(ci) } 192 193 // Children returns the four immediate children of this cell. 194 // If ci is a leaf cell, it returns four identical cells that are not the children. 195 func (ci CellID) Children() [4]CellID { 196 var ch [4]CellID 197 lsb := CellID(ci.lsb()) 198 ch[0] = ci - lsb + lsb>>2 199 lsb >>= 1 200 ch[1] = ch[0] + lsb 201 ch[2] = ch[1] + lsb 202 ch[3] = ch[2] + lsb 203 return ch 204 } 205 206 func sizeIJ(level int) int { 207 return 1 << uint(MaxLevel-level) 208 } 209 210 // EdgeNeighbors returns the four cells that are adjacent across the cell's four edges. 211 // Edges 0, 1, 2, 3 are in the down, right, up, left directions in the face space. 212 // All neighbors are guaranteed to be distinct. 213 func (ci CellID) EdgeNeighbors() [4]CellID { 214 level := ci.Level() 215 size := sizeIJ(level) 216 f, i, j, _ := ci.faceIJOrientation() 217 return [4]CellID{ 218 cellIDFromFaceIJWrap(f, i, j-size).Parent(level), 219 cellIDFromFaceIJWrap(f, i+size, j).Parent(level), 220 cellIDFromFaceIJWrap(f, i, j+size).Parent(level), 221 cellIDFromFaceIJWrap(f, i-size, j).Parent(level), 222 } 223 } 224 225 // VertexNeighbors returns the neighboring cellIDs with vertex closest to this cell at the given level. 226 // (Normally there are four neighbors, but the closest vertex may only have three neighbors if it is one of 227 // the 8 cube vertices.) 228 func (ci CellID) VertexNeighbors(level int) []CellID { 229 halfSize := sizeIJ(level + 1) 230 size := halfSize << 1 231 f, i, j, _ := ci.faceIJOrientation() 232 233 var isame, jsame bool 234 var ioffset, joffset int 235 if i&halfSize != 0 { 236 ioffset = size 237 isame = (i + size) < MaxSize 238 } else { 239 ioffset = -size 240 isame = (i - size) >= 0 241 } 242 if j&halfSize != 0 { 243 joffset = size 244 jsame = (j + size) < MaxSize 245 } else { 246 joffset = -size 247 jsame = (j - size) >= 0 248 } 249 250 results := []CellID{ 251 ci.Parent(level), 252 cellIDFromFaceIJSame(f, i+ioffset, j, isame).Parent(level), 253 cellIDFromFaceIJSame(f, i, j+joffset, jsame).Parent(level), 254 } 255 256 if isame || jsame { 257 results = append(results, cellIDFromFaceIJSame(f, i+ioffset, j+joffset, isame && jsame).Parent(level)) 258 } 259 260 return results 261 } 262 263 // AllNeighbors returns all neighbors of this cell at the given level. Two 264 // cells X and Y are neighbors if their boundaries intersect but their 265 // interiors do not. In particular, two cells that intersect at a single 266 // point are neighbors. Note that for cells adjacent to a face vertex, the 267 // same neighbor may be returned more than once. There could be up to eight 268 // neighbors including the diagonal ones that share the vertex. 269 // 270 // This requires level >= ci.Level(). 271 func (ci CellID) AllNeighbors(level int) []CellID { 272 var neighbors []CellID 273 274 face, i, j, _ := ci.faceIJOrientation() 275 276 // Find the coordinates of the lower left-hand leaf cell. We need to 277 // normalize (i,j) to a known position within the cell because level 278 // may be larger than this cell's level. 279 size := sizeIJ(ci.Level()) 280 i &= -size 281 j &= -size 282 283 nbrSize := sizeIJ(level) 284 285 // We compute the top-bottom, left-right, and diagonal neighbors in one 286 // pass. The loop test is at the end of the loop to avoid 32-bit overflow. 287 for k := -nbrSize; ; k += nbrSize { 288 var sameFace bool 289 if k < 0 { 290 sameFace = (j+k >= 0) 291 } else if k >= size { 292 sameFace = (j+k < MaxSize) 293 } else { 294 sameFace = true 295 // Top and bottom neighbors. 296 neighbors = append(neighbors, cellIDFromFaceIJSame(face, i+k, j-nbrSize, 297 j-size >= 0).Parent(level)) 298 neighbors = append(neighbors, cellIDFromFaceIJSame(face, i+k, j+size, 299 j+size < MaxSize).Parent(level)) 300 } 301 302 // Left, right, and diagonal neighbors. 303 neighbors = append(neighbors, cellIDFromFaceIJSame(face, i-nbrSize, j+k, 304 sameFace && i-size >= 0).Parent(level)) 305 neighbors = append(neighbors, cellIDFromFaceIJSame(face, i+size, j+k, 306 sameFace && i+size < MaxSize).Parent(level)) 307 308 if k >= size { 309 break 310 } 311 } 312 313 return neighbors 314 } 315 316 // RangeMin returns the minimum CellID that is contained within this cell. 317 func (ci CellID) RangeMin() CellID { return CellID(uint64(ci) - (ci.lsb() - 1)) } 318 319 // RangeMax returns the maximum CellID that is contained within this cell. 320 func (ci CellID) RangeMax() CellID { return CellID(uint64(ci) + (ci.lsb() - 1)) } 321 322 // Contains returns true iff the CellID contains oci. 323 func (ci CellID) Contains(oci CellID) bool { 324 return uint64(ci.RangeMin()) <= uint64(oci) && uint64(oci) <= uint64(ci.RangeMax()) 325 } 326 327 // Intersects returns true iff the CellID intersects oci. 328 func (ci CellID) Intersects(oci CellID) bool { 329 return uint64(oci.RangeMin()) <= uint64(ci.RangeMax()) && uint64(oci.RangeMax()) >= uint64(ci.RangeMin()) 330 } 331 332 // String returns the string representation of the cell ID in the form "1/3210". 333 func (ci CellID) String() string { 334 if !ci.IsValid() { 335 return "Invalid: " + strconv.FormatInt(int64(ci), 16) 336 } 337 var b bytes.Buffer 338 b.WriteByte("012345"[ci.Face()]) // values > 5 will have been picked off by !IsValid above 339 b.WriteByte('/') 340 for level := 1; level <= ci.Level(); level++ { 341 b.WriteByte("0123"[ci.ChildPosition(level)]) 342 } 343 return b.String() 344 } 345 346 // CellIDFromString returns a CellID from a string in the form "1/3210". 347 func CellIDFromString(s string) CellID { 348 level := len(s) - 2 349 if level < 0 || level > MaxLevel { 350 return CellID(0) 351 } 352 face := int(s[0] - '0') 353 if face < 0 || face > 5 || s[1] != '/' { 354 return CellID(0) 355 } 356 id := CellIDFromFace(face) 357 for i := 2; i < len(s); i++ { 358 childPos := s[i] - '0' 359 if childPos < 0 || childPos > 3 { 360 return CellID(0) 361 } 362 id = id.Children()[childPos] 363 } 364 return id 365 } 366 367 // Point returns the center of the s2 cell on the sphere as a Point. 368 // The maximum directional error in Point (compared to the exact 369 // mathematical result) is 1.5 * dblEpsilon radians, and the maximum length 370 // error is 2 * dblEpsilon (the same as Normalize). 371 func (ci CellID) Point() Point { return Point{ci.rawPoint().Normalize()} } 372 373 // LatLng returns the center of the s2 cell on the sphere as a LatLng. 374 func (ci CellID) LatLng() LatLng { return LatLngFromPoint(Point{ci.rawPoint()}) } 375 376 // ChildBegin returns the first child in a traversal of the children of this cell, in Hilbert curve order. 377 // 378 // for ci := c.ChildBegin(); ci != c.ChildEnd(); ci = ci.Next() { 379 // ... 380 // } 381 func (ci CellID) ChildBegin() CellID { 382 ol := ci.lsb() 383 return CellID(uint64(ci) - ol + ol>>2) 384 } 385 386 // ChildBeginAtLevel returns the first cell in a traversal of children a given level deeper than this cell, in 387 // Hilbert curve order. The given level must be no smaller than the cell's level. 388 // See ChildBegin for example use. 389 func (ci CellID) ChildBeginAtLevel(level int) CellID { 390 return CellID(uint64(ci) - ci.lsb() + lsbForLevel(level)) 391 } 392 393 // ChildEnd returns the first cell after a traversal of the children of this cell in Hilbert curve order. 394 // The returned cell may be invalid. 395 func (ci CellID) ChildEnd() CellID { 396 ol := ci.lsb() 397 return CellID(uint64(ci) + ol + ol>>2) 398 } 399 400 // ChildEndAtLevel returns the first cell after the last child in a traversal of children a given level deeper 401 // than this cell, in Hilbert curve order. 402 // The given level must be no smaller than the cell's level. 403 // The returned cell may be invalid. 404 func (ci CellID) ChildEndAtLevel(level int) CellID { 405 return CellID(uint64(ci) + ci.lsb() + lsbForLevel(level)) 406 } 407 408 // Next returns the next cell along the Hilbert curve. 409 // This is expected to be used with ChildBegin and ChildEnd, 410 // or ChildBeginAtLevel and ChildEndAtLevel. 411 func (ci CellID) Next() CellID { 412 return CellID(uint64(ci) + ci.lsb()<<1) 413 } 414 415 // Prev returns the previous cell along the Hilbert curve. 416 func (ci CellID) Prev() CellID { 417 return CellID(uint64(ci) - ci.lsb()<<1) 418 } 419 420 // NextWrap returns the next cell along the Hilbert curve, wrapping from last to 421 // first as necessary. This should not be used with ChildBegin and ChildEnd. 422 func (ci CellID) NextWrap() CellID { 423 n := ci.Next() 424 if uint64(n) < wrapOffset { 425 return n 426 } 427 return CellID(uint64(n) - wrapOffset) 428 } 429 430 // PrevWrap returns the previous cell along the Hilbert curve, wrapping around from 431 // first to last as necessary. This should not be used with ChildBegin and ChildEnd. 432 func (ci CellID) PrevWrap() CellID { 433 p := ci.Prev() 434 if uint64(p) < wrapOffset { 435 return p 436 } 437 return CellID(uint64(p) + wrapOffset) 438 } 439 440 // AdvanceWrap advances or retreats the indicated number of steps along the 441 // Hilbert curve at the current level and returns the new position. The 442 // position wraps between the first and last faces as necessary. 443 func (ci CellID) AdvanceWrap(steps int64) CellID { 444 if steps == 0 { 445 return ci 446 } 447 448 // We clamp the number of steps if necessary to ensure that we do not 449 // advance past the End() or before the Begin() of this level. 450 shift := uint(2*(MaxLevel-ci.Level()) + 1) 451 if steps < 0 { 452 if min := -int64(uint64(ci) >> shift); steps < min { 453 wrap := int64(wrapOffset >> shift) 454 steps %= wrap 455 if steps < min { 456 steps += wrap 457 } 458 } 459 } else { 460 // Unlike Advance(), we don't want to return End(level). 461 if max := int64((wrapOffset - uint64(ci)) >> shift); steps > max { 462 wrap := int64(wrapOffset >> shift) 463 steps %= wrap 464 if steps > max { 465 steps -= wrap 466 } 467 } 468 } 469 470 // If steps is negative, then shifting it left has undefined behavior. 471 // Cast to uint64 for a 2's complement answer. 472 return CellID(uint64(ci) + (uint64(steps) << shift)) 473 } 474 475 // Encode encodes the CellID. 476 func (ci CellID) Encode(w io.Writer) error { 477 e := &encoder{w: w} 478 ci.encode(e) 479 return e.err 480 } 481 482 func (ci CellID) encode(e *encoder) { 483 e.writeUint64(uint64(ci)) 484 } 485 486 // Decode decodes the CellID. 487 func (ci *CellID) Decode(r io.Reader) error { 488 d := &decoder{r: asByteReader(r)} 489 ci.decode(d) 490 return d.err 491 } 492 493 func (ci *CellID) decode(d *decoder) { 494 *ci = CellID(d.readUint64()) 495 } 496 497 // TODO: the methods below are not exported yet. Settle on the entire API design 498 // before doing this. Do we want to mirror the C++ one as closely as possible? 499 500 // distanceFromBegin returns the number of steps along the Hilbert curve that 501 // this cell is from the first node in the S2 hierarchy at our level. (i.e., 502 // FromFace(0).ChildBeginAtLevel(ci.Level())). This is analogous to Pos(), but 503 // for this cell's level. 504 // The return value is always non-negative. 505 func (ci CellID) distanceFromBegin() int64 { 506 return int64(ci >> uint64(2*(MaxLevel-ci.Level())+1)) 507 } 508 509 // rawPoint returns an unnormalized r3 vector from the origin through the center 510 // of the s2 cell on the sphere. 511 func (ci CellID) rawPoint() r3.Vector { 512 face, si, ti := ci.faceSiTi() 513 return faceUVToXYZ(face, stToUV((0.5/MaxSize)*float64(si)), stToUV((0.5/MaxSize)*float64(ti))) 514 } 515 516 // faceSiTi returns the Face/Si/Ti coordinates of the center of the cell. 517 func (ci CellID) faceSiTi() (face int, si, ti uint32) { 518 face, i, j, _ := ci.faceIJOrientation() 519 delta := 0 520 if ci.IsLeaf() { 521 delta = 1 522 } else { 523 if (i^(int(ci)>>2))&1 != 0 { 524 delta = 2 525 } 526 } 527 return face, uint32(2*i + delta), uint32(2*j + delta) 528 } 529 530 // faceIJOrientation uses the global lookupIJ table to unfiddle the bits of ci. 531 func (ci CellID) faceIJOrientation() (f, i, j, orientation int) { 532 f = ci.Face() 533 orientation = f & swapMask 534 nbits := MaxLevel - 7*lookupBits // first iteration 535 536 // Each iteration maps 8 bits of the Hilbert curve position into 537 // 4 bits of "i" and "j". The lookup table transforms a key of the 538 // form "ppppppppoo" to a value of the form "iiiijjjjoo", where the 539 // letters [ijpo] represents bits of "i", "j", the Hilbert curve 540 // position, and the Hilbert curve orientation respectively. 541 // 542 // On the first iteration we need to be careful to clear out the bits 543 // representing the cube face. 544 for k := 7; k >= 0; k-- { 545 orientation += (int(uint64(ci)>>uint64(k*2*lookupBits+1)) & ((1 << uint(2*nbits)) - 1)) << 2 546 orientation = lookupIJ[orientation] 547 i += (orientation >> (lookupBits + 2)) << uint(k*lookupBits) 548 j += ((orientation >> 2) & ((1 << lookupBits) - 1)) << uint(k*lookupBits) 549 orientation &= (swapMask | invertMask) 550 nbits = lookupBits // following iterations 551 } 552 553 // The position of a non-leaf cell at level "n" consists of a prefix of 554 // 2*n bits that identifies the cell, followed by a suffix of 555 // 2*(MaxLevel-n)+1 bits of the form 10*. If n==MaxLevel, the suffix is 556 // just "1" and has no effect. Otherwise, it consists of "10", followed 557 // by (MaxLevel-n-1) repetitions of "00", followed by "0". The "10" has 558 // no effect, while each occurrence of "00" has the effect of reversing 559 // the swapMask bit. 560 if ci.lsb()&0x1111111111111110 != 0 { 561 orientation ^= swapMask 562 } 563 564 return 565 } 566 567 // cellIDFromFaceIJ returns a leaf cell given its cube face (range 0..5) and IJ coordinates. 568 func cellIDFromFaceIJ(f, i, j int) CellID { 569 // Note that this value gets shifted one bit to the left at the end 570 // of the function. 571 n := uint64(f) << (PosBits - 1) 572 // Alternating faces have opposite Hilbert curve orientations; this 573 // is necessary in order for all faces to have a right-handed 574 // coordinate system. 575 bits := f & swapMask 576 // Each iteration maps 4 bits of "i" and "j" into 8 bits of the Hilbert 577 // curve position. The lookup table transforms a 10-bit key of the form 578 // "iiiijjjjoo" to a 10-bit value of the form "ppppppppoo", where the 579 // letters [ijpo] denote bits of "i", "j", Hilbert curve position, and 580 // Hilbert curve orientation respectively. 581 for k := 7; k >= 0; k-- { 582 mask := (1 << lookupBits) - 1 583 bits += ((i >> uint(k*lookupBits)) & mask) << (lookupBits + 2) 584 bits += ((j >> uint(k*lookupBits)) & mask) << 2 585 bits = lookupPos[bits] 586 n |= uint64(bits>>2) << (uint(k) * 2 * lookupBits) 587 bits &= (swapMask | invertMask) 588 } 589 return CellID(n*2 + 1) 590 } 591 592 func cellIDFromFaceIJWrap(f, i, j int) CellID { 593 // Convert i and j to the coordinates of a leaf cell just beyond the 594 // boundary of this face. This prevents 32-bit overflow in the case 595 // of finding the neighbors of a face cell. 596 i = clampInt(i, -1, MaxSize) 597 j = clampInt(j, -1, MaxSize) 598 599 // We want to wrap these coordinates onto the appropriate adjacent face. 600 // The easiest way to do this is to convert the (i,j) coordinates to (x,y,z) 601 // (which yields a point outside the normal face boundary), and then call 602 // xyzToFaceUV to project back onto the correct face. 603 // 604 // The code below converts (i,j) to (si,ti), and then (si,ti) to (u,v) using 605 // the linear projection (u=2*s-1 and v=2*t-1). (The code further below 606 // converts back using the inverse projection, s=0.5*(u+1) and t=0.5*(v+1). 607 // Any projection would work here, so we use the simplest.) We also clamp 608 // the (u,v) coordinates so that the point is barely outside the 609 // [-1,1]x[-1,1] face rectangle, since otherwise the reprojection step 610 // (which divides by the new z coordinate) might change the other 611 // coordinates enough so that we end up in the wrong leaf cell. 612 const scale = 1.0 / MaxSize 613 limit := math.Nextafter(1, 2) 614 u := math.Max(-limit, math.Min(limit, scale*float64((i<<1)+1-MaxSize))) 615 v := math.Max(-limit, math.Min(limit, scale*float64((j<<1)+1-MaxSize))) 616 617 // Find the leaf cell coordinates on the adjacent face, and convert 618 // them to a cell id at the appropriate level. 619 f, u, v = xyzToFaceUV(faceUVToXYZ(f, u, v)) 620 return cellIDFromFaceIJ(f, stToIJ(0.5*(u+1)), stToIJ(0.5*(v+1))) 621 } 622 623 func cellIDFromFaceIJSame(f, i, j int, sameFace bool) CellID { 624 if sameFace { 625 return cellIDFromFaceIJ(f, i, j) 626 } 627 return cellIDFromFaceIJWrap(f, i, j) 628 } 629 630 // ijToSTMin converts the i- or j-index of a leaf cell to the minimum corresponding 631 // s- or t-value contained by that cell. The argument must be in the range 632 // [0..2**30], i.e. up to one position beyond the normal range of valid leaf 633 // cell indices. 634 func ijToSTMin(i int) float64 { 635 return float64(i) / float64(MaxSize) 636 } 637 638 // stToIJ converts value in ST coordinates to a value in IJ coordinates. 639 func stToIJ(s float64) int { 640 return clampInt(int(math.Floor(MaxSize*s)), 0, MaxSize-1) 641 } 642 643 // cellIDFromPoint returns a leaf cell containing point p. Usually there is 644 // exactly one such cell, but for points along the edge of a cell, any 645 // adjacent cell may be (deterministically) chosen. This is because 646 // s2.CellIDs are considered to be closed sets. The returned cell will 647 // always contain the given point, i.e. 648 // 649 // CellFromPoint(p).ContainsPoint(p) 650 // 651 // is always true. 652 func cellIDFromPoint(p Point) CellID { 653 f, u, v := xyzToFaceUV(r3.Vector{p.X, p.Y, p.Z}) 654 i := stToIJ(uvToST(u)) 655 j := stToIJ(uvToST(v)) 656 return cellIDFromFaceIJ(f, i, j) 657 } 658 659 // ijLevelToBoundUV returns the bounds in (u,v)-space for the cell at the given 660 // level containing the leaf cell with the given (i,j)-coordinates. 661 func ijLevelToBoundUV(i, j, level int) r2.Rect { 662 cellSize := sizeIJ(level) 663 xLo := i & -cellSize 664 yLo := j & -cellSize 665 666 return r2.Rect{ 667 X: r1.Interval{ 668 Lo: stToUV(ijToSTMin(xLo)), 669 Hi: stToUV(ijToSTMin(xLo + cellSize)), 670 }, 671 Y: r1.Interval{ 672 Lo: stToUV(ijToSTMin(yLo)), 673 Hi: stToUV(ijToSTMin(yLo + cellSize)), 674 }, 675 } 676 } 677 678 // Constants related to the bit mangling in the Cell ID. 679 const ( 680 lookupBits = 4 681 swapMask = 0x01 682 invertMask = 0x02 683 ) 684 685 // The following lookup tables are used to convert efficiently between an 686 // (i,j) cell index and the corresponding position along the Hilbert curve. 687 // 688 // lookupPos maps 4 bits of "i", 4 bits of "j", and 2 bits representing the 689 // orientation of the current cell into 8 bits representing the order in which 690 // that subcell is visited by the Hilbert curve, plus 2 bits indicating the 691 // new orientation of the Hilbert curve within that subcell. (Cell 692 // orientations are represented as combination of swapMask and invertMask.) 693 // 694 // lookupIJ is an inverted table used for mapping in the opposite 695 // direction. 696 // 697 // We also experimented with looking up 16 bits at a time (14 bits of position 698 // plus 2 of orientation) but found that smaller lookup tables gave better 699 // performance. (2KB fits easily in the primary cache.) 700 var ( 701 ijToPos = [4][4]int{ 702 {0, 1, 3, 2}, // canonical order 703 {0, 3, 1, 2}, // axes swapped 704 {2, 3, 1, 0}, // bits inverted 705 {2, 1, 3, 0}, // swapped & inverted 706 } 707 posToIJ = [4][4]int{ 708 {0, 1, 3, 2}, // canonical order: (0,0), (0,1), (1,1), (1,0) 709 {0, 2, 3, 1}, // axes swapped: (0,0), (1,0), (1,1), (0,1) 710 {3, 2, 0, 1}, // bits inverted: (1,1), (1,0), (0,0), (0,1) 711 {3, 1, 0, 2}, // swapped & inverted: (1,1), (0,1), (0,0), (1,0) 712 } 713 posToOrientation = [4]int{swapMask, 0, 0, invertMask | swapMask} 714 lookupIJ [1 << (2*lookupBits + 2)]int 715 lookupPos [1 << (2*lookupBits + 2)]int 716 ) 717 718 func init() { 719 initLookupCell(0, 0, 0, 0, 0, 0) 720 initLookupCell(0, 0, 0, swapMask, 0, swapMask) 721 initLookupCell(0, 0, 0, invertMask, 0, invertMask) 722 initLookupCell(0, 0, 0, swapMask|invertMask, 0, swapMask|invertMask) 723 } 724 725 // initLookupCell initializes the lookupIJ table at init time. 726 func initLookupCell(level, i, j, origOrientation, pos, orientation int) { 727 if level == lookupBits { 728 ij := (i << lookupBits) + j 729 lookupPos[(ij<<2)+origOrientation] = (pos << 2) + orientation 730 lookupIJ[(pos<<2)+origOrientation] = (ij << 2) + orientation 731 return 732 } 733 734 level++ 735 i <<= 1 736 j <<= 1 737 pos <<= 2 738 r := posToIJ[orientation] 739 initLookupCell(level, i+(r[0]>>1), j+(r[0]&1), origOrientation, pos, orientation^posToOrientation[0]) 740 initLookupCell(level, i+(r[1]>>1), j+(r[1]&1), origOrientation, pos+1, orientation^posToOrientation[1]) 741 initLookupCell(level, i+(r[2]>>1), j+(r[2]&1), origOrientation, pos+2, orientation^posToOrientation[2]) 742 initLookupCell(level, i+(r[3]>>1), j+(r[3]&1), origOrientation, pos+3, orientation^posToOrientation[3]) 743 } 744 745 // CommonAncestorLevel returns the level of the common ancestor of the two S2 CellIDs. 746 func (ci CellID) CommonAncestorLevel(other CellID) (level int, ok bool) { 747 bits := uint64(ci ^ other) 748 if bits < ci.lsb() { 749 bits = ci.lsb() 750 } 751 if bits < other.lsb() { 752 bits = other.lsb() 753 } 754 755 msbPos := findMSBSetNonZero64(bits) 756 if msbPos > 60 { 757 return 0, false 758 } 759 return (60 - msbPos) >> 1, true 760 } 761 762 // Advance advances or retreats the indicated number of steps along the 763 // Hilbert curve at the current level, and returns the new position. The 764 // position is never advanced past End() or before Begin(). 765 func (ci CellID) Advance(steps int64) CellID { 766 if steps == 0 { 767 return ci 768 } 769 770 // We clamp the number of steps if necessary to ensure that we do not 771 // advance past the End() or before the Begin() of this level. Note that 772 // minSteps and maxSteps always fit in a signed 64-bit integer. 773 stepShift := uint(2*(MaxLevel-ci.Level()) + 1) 774 if steps < 0 { 775 minSteps := -int64(uint64(ci) >> stepShift) 776 if steps < minSteps { 777 steps = minSteps 778 } 779 } else { 780 maxSteps := int64((wrapOffset + ci.lsb() - uint64(ci)) >> stepShift) 781 if steps > maxSteps { 782 steps = maxSteps 783 } 784 } 785 return ci + CellID(steps)<<stepShift 786 } 787 788 // centerST return the center of the CellID in (s,t)-space. 789 func (ci CellID) centerST() r2.Point { 790 _, si, ti := ci.faceSiTi() 791 return r2.Point{siTiToST(si), siTiToST(ti)} 792 } 793 794 // sizeST returns the edge length of this CellID in (s,t)-space at the given level. 795 func (ci CellID) sizeST(level int) float64 { 796 return ijToSTMin(sizeIJ(level)) 797 } 798 799 // boundST returns the bound of this CellID in (s,t)-space. 800 func (ci CellID) boundST() r2.Rect { 801 s := ci.sizeST(ci.Level()) 802 return r2.RectFromCenterSize(ci.centerST(), r2.Point{s, s}) 803 } 804 805 // centerUV returns the center of this CellID in (u,v)-space. Note that 806 // the center of the cell is defined as the point at which it is recursively 807 // subdivided into four children; in general, it is not at the midpoint of 808 // the (u,v) rectangle covered by the cell. 809 func (ci CellID) centerUV() r2.Point { 810 _, si, ti := ci.faceSiTi() 811 return r2.Point{stToUV(siTiToST(si)), stToUV(siTiToST(ti))} 812 } 813 814 // boundUV returns the bound of this CellID in (u,v)-space. 815 func (ci CellID) boundUV() r2.Rect { 816 _, i, j, _ := ci.faceIJOrientation() 817 return ijLevelToBoundUV(i, j, ci.Level()) 818 } 819 820 // expandEndpoint returns a new u-coordinate u' such that the distance from the 821 // line u=u' to the given edge (u,v0)-(u,v1) is exactly the given distance 822 // (which is specified as the sine of the angle corresponding to the distance). 823 func expandEndpoint(u, maxV, sinDist float64) float64 { 824 // This is based on solving a spherical right triangle, similar to the 825 // calculation in Cap.RectBound. 826 // Given an edge of the form (u,v0)-(u,v1), let maxV = max(abs(v0), abs(v1)). 827 sinUShift := sinDist * math.Sqrt((1+u*u+maxV*maxV)/(1+u*u)) 828 cosUShift := math.Sqrt(1 - sinUShift*sinUShift) 829 // The following is an expansion of tan(atan(u) + asin(sinUShift)). 830 return (cosUShift*u + sinUShift) / (cosUShift - sinUShift*u) 831 } 832 833 // expandedByDistanceUV returns a rectangle expanded in (u,v)-space so that it 834 // contains all points within the given distance of the boundary, and return the 835 // smallest such rectangle. If the distance is negative, then instead shrink this 836 // rectangle so that it excludes all points within the given absolute distance 837 // of the boundary. 838 // 839 // Distances are measured *on the sphere*, not in (u,v)-space. For example, 840 // you can use this method to expand the (u,v)-bound of an CellID so that 841 // it contains all points within 5km of the original cell. You can then 842 // test whether a point lies within the expanded bounds like this: 843 // 844 // if u, v, ok := faceXYZtoUV(face, point); ok && bound.ContainsPoint(r2.Point{u,v}) { ... } 845 // 846 // Limitations: 847 // 848 // - Because the rectangle is drawn on one of the six cube-face planes 849 // (i.e., {x,y,z} = +/-1), it can cover at most one hemisphere. This 850 // limits the maximum amount that a rectangle can be expanded. For 851 // example, CellID bounds can be expanded safely by at most 45 degrees 852 // (about 5000 km on the Earth's surface). 853 // 854 // - The implementation is not exact for negative distances. The resulting 855 // rectangle will exclude all points within the given distance of the 856 // boundary but may be slightly smaller than necessary. 857 func expandedByDistanceUV(uv r2.Rect, distance s1.Angle) r2.Rect { 858 // Expand each of the four sides of the rectangle just enough to include all 859 // points within the given distance of that side. (The rectangle may be 860 // expanded by a different amount in (u,v)-space on each side.) 861 maxU := math.Max(math.Abs(uv.X.Lo), math.Abs(uv.X.Hi)) 862 maxV := math.Max(math.Abs(uv.Y.Lo), math.Abs(uv.Y.Hi)) 863 sinDist := math.Sin(float64(distance)) 864 return r2.Rect{ 865 X: r1.Interval{expandEndpoint(uv.X.Lo, maxV, -sinDist), 866 expandEndpoint(uv.X.Hi, maxV, sinDist)}, 867 Y: r1.Interval{expandEndpoint(uv.Y.Lo, maxU, -sinDist), 868 expandEndpoint(uv.Y.Hi, maxU, sinDist)}} 869 } 870 871 // MaxTile returns the largest cell with the same RangeMin such that 872 // RangeMax < limit.RangeMin. It returns limit if no such cell exists. 873 // This method can be used to generate a small set of CellIDs that covers 874 // a given range (a tiling). This example shows how to generate a tiling 875 // for a semi-open range of leaf cells [start, limit): 876 // 877 // for id := start.MaxTile(limit); id != limit; id = id.Next().MaxTile(limit)) { ... } 878 // 879 // Note that in general the cells in the tiling will be of different sizes; 880 // they gradually get larger (near the middle of the range) and then 881 // gradually get smaller as limit is approached. 882 func (ci CellID) MaxTile(limit CellID) CellID { 883 start := ci.RangeMin() 884 if start >= limit.RangeMin() { 885 return limit 886 } 887 888 if ci.RangeMax() >= limit { 889 // The cell is too large, shrink it. Note that when generating coverings 890 // of CellID ranges, this loop usually executes only once. Also because 891 // ci.RangeMin() < limit.RangeMin(), we will always exit the loop by the 892 // time we reach a leaf cell. 893 for { 894 ci = ci.Children()[0] 895 if ci.RangeMax() < limit { 896 break 897 } 898 } 899 return ci 900 } 901 902 // The cell may be too small. Grow it if necessary. Note that generally 903 // this loop only iterates once. 904 for !ci.isFace() { 905 parent := ci.immediateParent() 906 if parent.RangeMin() != start || parent.RangeMax() >= limit { 907 break 908 } 909 ci = parent 910 } 911 return ci 912 } 913 914 // centerFaceSiTi returns the (face, si, ti) coordinates of the center of the cell. 915 // Note that although (si,ti) coordinates span the range [0,2**31] in general, 916 // the cell center coordinates are always in the range [1,2**31-1] and 917 // therefore can be represented using a signed 32-bit integer. 918 func (ci CellID) centerFaceSiTi() (face, si, ti int) { 919 // First we compute the discrete (i,j) coordinates of a leaf cell contained 920 // within the given cell. Given that cells are represented by the Hilbert 921 // curve position corresponding at their center, it turns out that the cell 922 // returned by faceIJOrientation is always one of two leaf cells closest 923 // to the center of the cell (unless the given cell is a leaf cell itself, 924 // in which case there is only one possibility). 925 // 926 // Given a cell of size s >= 2 (i.e. not a leaf cell), and letting (imin, 927 // jmin) be the coordinates of its lower left-hand corner, the leaf cell 928 // returned by faceIJOrientation is either (imin + s/2, jmin + s/2) 929 // (imin + s/2 - 1, jmin + s/2 - 1). The first case is the one we want. 930 // We can distinguish these two cases by looking at the low bit of i or 931 // j. In the second case the low bit is one, unless s == 2 (i.e. the 932 // level just above leaf cells) in which case the low bit is zero. 933 // 934 // In the code below, the expression ((i ^ (int(id) >> 2)) & 1) is true 935 // if we are in the second case described above. 936 face, i, j, _ := ci.faceIJOrientation() 937 delta := 0 938 if ci.IsLeaf() { 939 delta = 1 940 } else if (int64(i)^(int64(ci)>>2))&1 == 1 { 941 delta = 2 942 } 943 944 // Note that (2 * {i,j} + delta) will never overflow a 32-bit integer. 945 return face, 2*i + delta, 2*j + delta 946 } 947