1 /* 2 * Copyright 2005 Google Inc. 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 */ 16 package com.google.common.geometry; 17 18 import java.util.List; 19 import java.util.Locale; 20 21 /** 22 * An S2CellId is a 64-bit unsigned integer that uniquely identifies a cell in 23 * the S2 cell decomposition. It has the following format: 24 * 25 * <pre> 26 * id = [face][face_pos] 27 * </pre> 28 * 29 * face: a 3-bit number (range 0..5) encoding the cube face. 30 * 31 * face_pos: a 61-bit number encoding the position of the center of this cell 32 * along the Hilbert curve over this face (see the Wiki pages for details). 33 * 34 * Sequentially increasing cell ids follow a continuous space-filling curve over 35 * the entire sphere. They have the following properties: 36 * - The id of a cell at level k consists of a 3-bit face number followed by k 37 * bit pairs that recursively select one of the four children of each cell. The 38 * next bit is always 1, and all other bits are 0. Therefore, the level of a 39 * cell is determined by the position of its lowest-numbered bit that is turned 40 * on (for a cell at level k, this position is 2 * (MAX_LEVEL - k).) 41 * - The id of a parent cell is at the midpoint of the range of ids spanned by 42 * its children (or by its descendants at any level). 43 * 44 * Leaf cells are often used to represent points on the unit sphere, and this 45 * class provides methods for converting directly between these two 46 * representations. For cells that represent 2D regions rather than discrete 47 * point, it is better to use the S2Cell class. 48 * 49 * 50 */ 51 public final strictfp class S2CellId implements Comparable<S2CellId> { 52 53 // Although only 60 bits are needed to represent the index of a leaf 54 // cell, we need an extra bit in order to represent the position of 55 // the center of the leaf cell along the Hilbert curve. 56 public static final int FACE_BITS = 3; 57 public static final int NUM_FACES = 6; 58 public static final int MAX_LEVEL = 30; // Valid levels: 0..MAX_LEVEL 59 public static final int POS_BITS = 2 * MAX_LEVEL + 1; 60 public static final int MAX_SIZE = 1 << MAX_LEVEL; 61 62 // Constant related to unsigned long's 63 public static final long MAX_UNSIGNED = -1L; // Equivalent to 0xffffffffffffffffL 64 65 // The following lookup tables are used to convert efficiently between an 66 // (i,j) cell index and the corresponding position along the Hilbert curve. 67 // "lookup_pos" maps 4 bits of "i", 4 bits of "j", and 2 bits representing the 68 // orientation of the current cell into 8 bits representing the order in which 69 // that subcell is visited by the Hilbert curve, plus 2 bits indicating the 70 // new orientation of the Hilbert curve within that subcell. (Cell 71 // orientations are represented as combination of kSwapMask and kInvertMask.) 72 // 73 // "lookup_ij" is an inverted table used for mapping in the opposite 74 // direction. 75 // 76 // We also experimented with looking up 16 bits at a time (14 bits of position 77 // plus 2 of orientation) but found that smaller lookup tables gave better 78 // performance. (2KB fits easily in the primary cache.) 79 80 81 // Values for these constants are *declared* in the *.h file. Even though 82 // the declaration specifies a value for the constant, that declaration 83 // is not a *definition* of storage for the value. Because the values are 84 // supplied in the declaration, we don't need the values here. Failing to 85 // define storage causes link errors for any code that tries to take the 86 // address of one of these values. 87 private static final int LOOKUP_BITS = 4; 88 private static final int SWAP_MASK = 0x01; 89 private static final int INVERT_MASK = 0x02; 90 91 private static final int[] LOOKUP_POS = new int[1 << (2 * LOOKUP_BITS + 2)]; 92 private static final int[] LOOKUP_IJ = new int[1 << (2 * LOOKUP_BITS + 2)]; 93 94 /** 95 * This is the offset required to wrap around from the beginning of the 96 * Hilbert curve to the end or vice versa; see next_wrap() and prev_wrap(). 97 */ 98 private static final long WRAP_OFFSET = (long) (NUM_FACES) << POS_BITS; 99 100 static { 101 initLookupCell(0, 0, 0, 0, 0, 0); 102 initLookupCell(0, 0, 0, SWAP_MASK, 0, SWAP_MASK); 103 initLookupCell(0, 0, 0, INVERT_MASK, 0, INVERT_MASK); 104 initLookupCell(0, 0, 0, SWAP_MASK | INVERT_MASK, 0, SWAP_MASK | INVERT_MASK); 105 } 106 107 /** 108 * The id of the cell. 109 */ 110 private final long id; 111 S2CellId(long id)112 public S2CellId(long id) { 113 this.id = id; 114 } 115 S2CellId()116 public S2CellId() { 117 this.id = 0; 118 } 119 120 /** The default constructor returns an invalid cell id. */ none()121 public static S2CellId none() { 122 return new S2CellId(); 123 } 124 125 /** 126 * Returns an invalid cell id guaranteed to be larger than any valid cell id. 127 * Useful for creating indexes. 128 */ sentinel()129 public static S2CellId sentinel() { 130 return new S2CellId(MAX_UNSIGNED); // -1 131 } 132 133 /** 134 * Return a cell given its face (range 0..5), 61-bit Hilbert curve position 135 * within that face, and level (range 0..MAX_LEVEL). The given position will 136 * be modified to correspond to the Hilbert curve position at the center of 137 * the returned cell. This is a static function rather than a constructor in 138 * order to give names to the arguments. 139 */ fromFacePosLevel(int face, long pos, int level)140 public static S2CellId fromFacePosLevel(int face, long pos, int level) { 141 return new S2CellId((((long) face) << POS_BITS) + (pos | 1)).parent(level); 142 } 143 144 /** 145 * Return the leaf cell containing the given point (a direction vector, not 146 * necessarily unit length). 147 */ fromPoint(S2Point p)148 public static S2CellId fromPoint(S2Point p) { 149 int face = S2Projections.xyzToFace(p); 150 R2Vector uv = S2Projections.validFaceXyzToUv(face, p); 151 int i = stToIJ(S2Projections.uvToST(uv.x())); 152 int j = stToIJ(S2Projections.uvToST(uv.y())); 153 return fromFaceIJ(face, i, j); 154 } 155 156 157 /** Return the leaf cell containing the given S2LatLng. */ fromLatLng(S2LatLng ll)158 public static S2CellId fromLatLng(S2LatLng ll) { 159 return fromPoint(ll.toPoint()); 160 } 161 toPoint()162 public S2Point toPoint() { 163 return S2Point.normalize(toPointRaw()); 164 } 165 166 /** 167 * Return the direction vector corresponding to the center of the given cell. 168 * The vector returned by ToPointRaw is not necessarily unit length. 169 */ toPointRaw()170 public S2Point toPointRaw() { 171 // First we compute the discrete (i,j) coordinates of a leaf cell contained 172 // within the given cell. Given that cells are represented by the Hilbert 173 // curve position corresponding at their center, it turns out that the cell 174 // returned by ToFaceIJOrientation is always one of two leaf cells closest 175 // to the center of the cell (unless the given cell is a leaf cell itself, 176 // in which case there is only one possibility). 177 // 178 // Given a cell of size s >= 2 (i.e. not a leaf cell), and letting (imin, 179 // jmin) be the coordinates of its lower left-hand corner, the leaf cell 180 // returned by ToFaceIJOrientation() is either (imin + s/2, jmin + s/2) 181 // (imin + s/2 - 1, jmin + s/2 - 1). We can distinguish these two cases by 182 // looking at the low bit of "i" or "j". In the first case the low bit is 183 // zero, unless s == 2 (i.e. the level just above leaf cells) in which case 184 // the low bit is one. 185 // 186 // The following calculation converts (i,j) to the (si,ti) coordinates of 187 // the cell center. (We need to multiply the coordinates by a factor of 2 188 // so that the center of leaf cells can be represented exactly.) 189 190 MutableInteger i = new MutableInteger(0); 191 MutableInteger j = new MutableInteger(0); 192 int face = toFaceIJOrientation(i, j, null); 193 // System.out.println("i= " + i.intValue() + " j = " + j.intValue()); 194 int delta = isLeaf() ? 1 : (((i.intValue() ^ (((int) id) >>> 2)) & 1) != 0) 195 ? 2 : 0; 196 int si = (i.intValue() << 1) + delta - MAX_SIZE; 197 int ti = (j.intValue() << 1) + delta - MAX_SIZE; 198 return faceSiTiToXYZ(face, si, ti); 199 } 200 201 /** Return the S2LatLng corresponding to the center of the given cell. */ toLatLng()202 public S2LatLng toLatLng() { 203 return new S2LatLng(toPointRaw()); 204 } 205 206 207 /** The 64-bit unique identifier for this cell. */ id()208 public long id() { 209 return id; 210 } 211 212 /** Return true if id() represents a valid cell. */ isValid()213 public boolean isValid() { 214 return face() < NUM_FACES && ((lowestOnBit() & (0x1555555555555555L)) != 0); 215 } 216 217 /** Which cube face this cell belongs to, in the range 0..5. */ face()218 public int face() { 219 return (int) (id >>> POS_BITS); 220 } 221 222 /** 223 * The position of the cell center along the Hilbert curve over this face, in 224 * the range 0..(2**kPosBits-1). 225 */ pos()226 public long pos() { 227 return (id & (-1L >>> FACE_BITS)); 228 } 229 230 /** Return the subdivision level of the cell (range 0..MAX_LEVEL). */ level()231 public int level() { 232 // Fast path for leaf cells. 233 if (isLeaf()) { 234 return MAX_LEVEL; 235 } 236 int x = ((int) id); 237 int level = -1; 238 if (x != 0) { 239 level += 16; 240 } else { 241 x = (int) (id >>> 32); 242 } 243 // We only need to look at even-numbered bits to determine the 244 // level of a valid cell id. 245 x &= -x; // Get lowest bit. 246 if ((x & 0x00005555) != 0) { 247 level += 8; 248 } 249 if ((x & 0x00550055) != 0) { 250 level += 4; 251 } 252 if ((x & 0x05050505) != 0) { 253 level += 2; 254 } 255 if ((x & 0x11111111) != 0) { 256 level += 1; 257 } 258 // assert (level >= 0 && level <= MAX_LEVEL); 259 return level; 260 } 261 262 263 264 /** 265 * Return true if this is a leaf cell (more efficient than checking whether 266 * level() == MAX_LEVEL). 267 */ isLeaf()268 public boolean isLeaf() { 269 return ((int) id & 1) != 0; 270 } 271 272 /** 273 * Return true if this is a top-level face cell (more efficient than checking 274 * whether level() == 0). 275 */ isFace()276 public boolean isFace() { 277 return (id & (lowestOnBitForLevel(0) - 1)) == 0; 278 } 279 280 /** 281 * Return the child position (0..3) of this cell's ancestor at the given 282 * level, relative to its parent. The argument should be in the range 283 * 1..MAX_LEVEL. For example, child_position(1) returns the position of this 284 * cell's level-1 ancestor within its top-level face cell. 285 */ childPosition(int level)286 public int childPosition(int level) { 287 return (int) (id >>> (2 * (MAX_LEVEL - level) + 1)) & 3; 288 } 289 290 // Methods that return the range of cell ids that are contained 291 // within this cell (including itself). The range is *inclusive* 292 // (i.e. test using >= and <=) and the return values of both 293 // methods are valid leaf cell ids. 294 // 295 // These methods should not be used for iteration. If you want to 296 // iterate through all the leaf cells, call child_begin(MAX_LEVEL) and 297 // child_end(MAX_LEVEL) instead. 298 // 299 // It would in fact be error-prone to define a range_end() method, 300 // because (range_max().id() + 1) is not always a valid cell id, and the 301 // iterator would need to be tested using "<" rather that the usual "!=". rangeMin()302 public S2CellId rangeMin() { 303 return new S2CellId(id - (lowestOnBit() - 1)); 304 } 305 rangeMax()306 public S2CellId rangeMax() { 307 return new S2CellId(id + (lowestOnBit() - 1)); 308 } 309 310 311 /** Return true if the given cell is contained within this one. */ contains(S2CellId other)312 public boolean contains(S2CellId other) { 313 // assert (isValid() && other.isValid()); 314 return other.greaterOrEquals(rangeMin()) && other.lessOrEquals(rangeMax()); 315 } 316 317 /** Return true if the given cell intersects this one. */ intersects(S2CellId other)318 public boolean intersects(S2CellId other) { 319 // assert (isValid() && other.isValid()); 320 return other.rangeMin().lessOrEquals(rangeMax()) 321 && other.rangeMax().greaterOrEquals(rangeMin()); 322 } 323 parent()324 public S2CellId parent() { 325 // assert (isValid() && level() > 0); 326 long newLsb = lowestOnBit() << 2; 327 return new S2CellId((id & -newLsb) | newLsb); 328 } 329 330 /** 331 * Return the cell at the previous level or at the given level (which must be 332 * less than or equal to the current level). 333 */ parent(int level)334 public S2CellId parent(int level) { 335 // assert (isValid() && level >= 0 && level <= this.level()); 336 long newLsb = lowestOnBitForLevel(level); 337 return new S2CellId((id & -newLsb) | newLsb); 338 } 339 childBegin()340 public S2CellId childBegin() { 341 // assert (isValid() && level() < MAX_LEVEL); 342 long oldLsb = lowestOnBit(); 343 return new S2CellId(id - oldLsb + (oldLsb >>> 2)); 344 } 345 childBegin(int level)346 public S2CellId childBegin(int level) { 347 // assert (isValid() && level >= this.level() && level <= MAX_LEVEL); 348 return new S2CellId(id - lowestOnBit() + lowestOnBitForLevel(level)); 349 } 350 childEnd()351 public S2CellId childEnd() { 352 // assert (isValid() && level() < MAX_LEVEL); 353 long oldLsb = lowestOnBit(); 354 return new S2CellId(id + oldLsb + (oldLsb >>> 2)); 355 } 356 childEnd(int level)357 public S2CellId childEnd(int level) { 358 // assert (isValid() && level >= this.level() && level <= MAX_LEVEL); 359 return new S2CellId(id + lowestOnBit() + lowestOnBitForLevel(level)); 360 } 361 362 // Iterator-style methods for traversing the immediate children of a cell or 363 // all of the children at a given level (greater than or equal to the current 364 // level). Note that the end value is exclusive, just like standard STL 365 // iterators, and may not even be a valid cell id. You should iterate using 366 // code like this: 367 // 368 // for(S2CellId c = id.childBegin(); !c.equals(id.childEnd()); c = c.next()) 369 // ... 370 // 371 // The convention for advancing the iterator is "c = c.next()", so be sure 372 // to use 'equals()' in the loop guard, or compare 64-bit cell id's, 373 // rather than "c != id.childEnd()". 374 375 /** 376 * Return the next cell at the same level along the Hilbert curve. Works 377 * correctly when advancing from one face to the next, but does *not* wrap 378 * around from the last face to the first or vice versa. 379 */ next()380 public S2CellId next() { 381 return new S2CellId(id + (lowestOnBit() << 1)); 382 } 383 384 /** 385 * Return the previous cell at the same level along the Hilbert curve. Works 386 * correctly when advancing from one face to the next, but does *not* wrap 387 * around from the last face to the first or vice versa. 388 */ prev()389 public S2CellId prev() { 390 return new S2CellId(id - (lowestOnBit() << 1)); 391 } 392 393 394 /** 395 * Like next(), but wraps around from the last face to the first and vice 396 * versa. Should *not* be used for iteration in conjunction with 397 * child_begin(), child_end(), Begin(), or End(). 398 */ nextWrap()399 public S2CellId nextWrap() { 400 S2CellId n = next(); 401 if (unsignedLongLessThan(n.id, WRAP_OFFSET)) { 402 return n; 403 } 404 return new S2CellId(n.id - WRAP_OFFSET); 405 } 406 407 /** 408 * Like prev(), but wraps around from the last face to the first and vice 409 * versa. Should *not* be used for iteration in conjunction with 410 * child_begin(), child_end(), Begin(), or End(). 411 */ prevWrap()412 public S2CellId prevWrap() { 413 S2CellId p = prev(); 414 if (p.id < WRAP_OFFSET) { 415 return p; 416 } 417 return new S2CellId(p.id + WRAP_OFFSET); 418 } 419 420 begin(int level)421 public static S2CellId begin(int level) { 422 return fromFacePosLevel(0, 0, 0).childBegin(level); 423 } 424 end(int level)425 public static S2CellId end(int level) { 426 return fromFacePosLevel(5, 0, 0).childEnd(level); 427 } 428 429 430 /** 431 * Decodes the cell id from a compact text string suitable for display or 432 * indexing. Cells at lower levels (i.e. larger cells) are encoded into 433 * fewer characters. The maximum token length is 16. 434 * 435 * @param token the token to decode 436 * @return the S2CellId for that token 437 * @throws NumberFormatException if the token is not formatted correctly 438 */ fromToken(String token)439 public static S2CellId fromToken(String token) { 440 if (token == null) { 441 throw new NumberFormatException("Null string in S2CellId.fromToken"); 442 } 443 if (token.length() == 0) { 444 throw new NumberFormatException("Empty string in S2CellId.fromToken"); 445 } 446 if (token.length() > 16 || "X".equals(token)) { 447 return none(); 448 } 449 450 long value = 0; 451 for (int pos = 0; pos < 16; pos++) { 452 int digit = 0; 453 if (pos < token.length()) { 454 digit = Character.digit(token.charAt(pos), 16); 455 if (digit == -1) { 456 throw new NumberFormatException(token); 457 } 458 if (overflowInParse(value, digit)) { 459 throw new NumberFormatException("Too large for unsigned long: " + token); 460 } 461 } 462 value = (value * 16) + digit; 463 } 464 465 return new S2CellId(value); 466 } 467 468 /** 469 * Encodes the cell id to compact text strings suitable for display or indexing. 470 * Cells at lower levels (i.e. larger cells) are encoded into fewer characters. 471 * The maximum token length is 16. 472 * 473 * Simple implementation: convert the id to hex and strip trailing zeros. We 474 * could use base-32 or base-64, but assuming the cells used for indexing 475 * regions are at least 100 meters across (level 16 or less), the savings 476 * would be at most 3 bytes (9 bytes hex vs. 6 bytes base-64). 477 * 478 * @return the encoded cell id 479 */ toToken()480 public String toToken() { 481 if (id == 0) { 482 return "X"; 483 } 484 485 String hex = Long.toHexString(id).toLowerCase(Locale.ENGLISH); 486 StringBuilder sb = new StringBuilder(16); 487 for (int i = hex.length(); i < 16; i++) { 488 sb.append('0'); 489 } 490 sb.append(hex); 491 for (int len = 16; len > 0; len--) { 492 if (sb.charAt(len - 1) != '0') { 493 return sb.substring(0, len); 494 } 495 } 496 497 throw new RuntimeException("Shouldn't make it here"); 498 } 499 500 /** 501 * Returns true if (current * 10) + digit is a number too large to be 502 * represented by an unsigned long. This is useful for detecting overflow 503 * while parsing a string representation of a number. 504 */ overflowInParse(long current, int digit)505 private static boolean overflowInParse(long current, int digit) { 506 return overflowInParse(current, digit, 10); 507 } 508 509 /** 510 * Returns true if (current * radix) + digit is a number too large to be 511 * represented by an unsigned long. This is useful for detecting overflow 512 * while parsing a string representation of a number. 513 * Does not verify whether supplied radix is valid, passing an invalid radix 514 * will give undefined results or an ArrayIndexOutOfBoundsException. 515 */ overflowInParse(long current, int digit, int radix)516 private static boolean overflowInParse(long current, int digit, int radix) { 517 if (current >= 0) { 518 if (current < maxValueDivs[radix]) { 519 return false; 520 } 521 if (current > maxValueDivs[radix]) { 522 return true; 523 } 524 // current == maxValueDivs[radix] 525 return (digit > maxValueMods[radix]); 526 } 527 528 // current < 0: high bit is set 529 return true; 530 } 531 532 // calculated as 0xffffffffffffffff / radix 533 private static final long maxValueDivs[] = {0, 0, // 0 and 1 are invalid 534 9223372036854775807L, 6148914691236517205L, 4611686018427387903L, // 2-4 535 3689348814741910323L, 3074457345618258602L, 2635249153387078802L, // 5-7 536 2305843009213693951L, 2049638230412172401L, 1844674407370955161L, // 8-10 537 1676976733973595601L, 1537228672809129301L, 1418980313362273201L, // 11-13 538 1317624576693539401L, 1229782938247303441L, 1152921504606846975L, // 14-16 539 1085102592571150095L, 1024819115206086200L, 970881267037344821L, // 17-19 540 922337203685477580L, 878416384462359600L, 838488366986797800L, // 20-22 541 802032351030850070L, 768614336404564650L, 737869762948382064L, // 23-25 542 709490156681136600L, 683212743470724133L, 658812288346769700L, // 26-28 543 636094623231363848L, 614891469123651720L, 595056260442243600L, // 29-31 544 576460752303423487L, 558992244657865200L, 542551296285575047L, // 32-34 545 527049830677415760L, 512409557603043100L }; // 35-36 546 547 // calculated as 0xffffffffffffffff % radix 548 private static final int maxValueMods[] = {0, 0, // 0 and 1 are invalid 549 1, 0, 3, 0, 3, 1, 7, 6, 5, 4, 3, 2, 1, 0, 15, 0, 15, 16, 15, 15, // 2-21 550 15, 5, 15, 15, 15, 24, 15, 23, 15, 15, 31, 15, 17, 15, 15 }; // 22-36 551 552 /** 553 * Return the four cells that are adjacent across the cell's four edges. 554 * Neighbors are returned in the order defined by S2Cell::GetEdge. All 555 * neighbors are guaranteed to be distinct. 556 */ getEdgeNeighbors(S2CellId neighbors[])557 public void getEdgeNeighbors(S2CellId neighbors[]) { 558 559 MutableInteger i = new MutableInteger(0); 560 MutableInteger j = new MutableInteger(0); 561 562 int level = this.level(); 563 int size = 1 << (MAX_LEVEL - level); 564 int face = toFaceIJOrientation(i, j, null); 565 566 // Edges 0, 1, 2, 3 are in the S, E, N, W directions. 567 neighbors[0] = fromFaceIJSame(face, i.intValue(), j.intValue() - size, 568 j.intValue() - size >= 0).parent(level); 569 neighbors[1] = fromFaceIJSame(face, i.intValue() + size, j.intValue(), 570 i.intValue() + size < MAX_SIZE).parent(level); 571 neighbors[2] = fromFaceIJSame(face, i.intValue(), j.intValue() + size, 572 j.intValue() + size < MAX_SIZE).parent(level); 573 neighbors[3] = fromFaceIJSame(face, i.intValue() - size, j.intValue(), 574 i.intValue() - size >= 0).parent(level); 575 } 576 577 /** 578 * Return the neighbors of closest vertex to this cell at the given level, by 579 * appending them to "output". Normally there are four neighbors, but the 580 * closest vertex may only have three neighbors if it is one of the 8 cube 581 * vertices. 582 * 583 * Requires: level < this.evel(), so that we can determine which vertex is 584 * closest (in particular, level == MAX_LEVEL is not allowed). 585 */ getVertexNeighbors(int level, List<S2CellId> output)586 public void getVertexNeighbors(int level, List<S2CellId> output) { 587 // "level" must be strictly less than this cell's level so that we can 588 // determine which vertex this cell is closest to. 589 // assert (level < this.level()); 590 MutableInteger i = new MutableInteger(0); 591 MutableInteger j = new MutableInteger(0); 592 int face = toFaceIJOrientation(i, j, null); 593 594 // Determine the i- and j-offsets to the closest neighboring cell in each 595 // direction. This involves looking at the next bit of "i" and "j" to 596 // determine which quadrant of this->parent(level) this cell lies in. 597 int halfsize = 1 << (MAX_LEVEL - (level + 1)); 598 int size = halfsize << 1; 599 boolean isame, jsame; 600 int ioffset, joffset; 601 if ((i.intValue() & halfsize) != 0) { 602 ioffset = size; 603 isame = (i.intValue() + size) < MAX_SIZE; 604 } else { 605 ioffset = -size; 606 isame = (i.intValue() - size) >= 0; 607 } 608 if ((j.intValue() & halfsize) != 0) { 609 joffset = size; 610 jsame = (j.intValue() + size) < MAX_SIZE; 611 } else { 612 joffset = -size; 613 jsame = (j.intValue() - size) >= 0; 614 } 615 616 output.add(parent(level)); 617 output 618 .add(fromFaceIJSame(face, i.intValue() + ioffset, j.intValue(), isame) 619 .parent(level)); 620 output 621 .add(fromFaceIJSame(face, i.intValue(), j.intValue() + joffset, jsame) 622 .parent(level)); 623 // If i- and j- edge neighbors are *both* on a different face, then this 624 // vertex only has three neighbors (it is one of the 8 cube vertices). 625 if (isame || jsame) { 626 output.add(fromFaceIJSame(face, i.intValue() + ioffset, 627 j.intValue() + joffset, isame && jsame).parent(level)); 628 } 629 } 630 631 /** 632 * Append all neighbors of this cell at the given level to "output". Two cells 633 * X and Y are neighbors if their boundaries intersect but their interiors do 634 * not. In particular, two cells that intersect at a single point are 635 * neighbors. 636 * 637 * Requires: nbr_level >= this->level(). Note that for cells adjacent to a 638 * face vertex, the same neighbor may be appended more than once. 639 */ getAllNeighbors(int nbrLevel, List<S2CellId> output)640 public void getAllNeighbors(int nbrLevel, List<S2CellId> output) { 641 MutableInteger i = new MutableInteger(0); 642 MutableInteger j = new MutableInteger(0); 643 644 int face = toFaceIJOrientation(i, j, null); 645 646 // Find the coordinates of the lower left-hand leaf cell. We need to 647 // normalize (i,j) to a known position within the cell because nbr_level 648 // may be larger than this cell's level. 649 int size = 1 << (MAX_LEVEL - level()); 650 i.setValue(i.intValue() & -size); 651 j.setValue(j.intValue() & -size); 652 653 int nbrSize = 1 << (MAX_LEVEL - nbrLevel); 654 // assert (nbrSize <= size); 655 656 // We compute the N-S, E-W, and diagonal neighbors in one pass. 657 // The loop test is at the end of the loop to avoid 32-bit overflow. 658 for (int k = -nbrSize;; k += nbrSize) { 659 boolean sameFace; 660 if (k < 0) { 661 sameFace = (j.intValue() + k >= 0); 662 } else if (k >= size) { 663 sameFace = (j.intValue() + k < MAX_SIZE); 664 } else { 665 sameFace = true; 666 // North and South neighbors. 667 output.add(fromFaceIJSame(face, i.intValue() + k, 668 j.intValue() - nbrSize, j.intValue() - size >= 0).parent(nbrLevel)); 669 output.add(fromFaceIJSame(face, i.intValue() + k, j.intValue() + size, 670 j.intValue() + size < MAX_SIZE).parent(nbrLevel)); 671 } 672 // East, West, and Diagonal neighbors. 673 output.add(fromFaceIJSame(face, i.intValue() - nbrSize, 674 j.intValue() + k, sameFace && i.intValue() - size >= 0).parent( 675 nbrLevel)); 676 output.add(fromFaceIJSame(face, i.intValue() + size, j.intValue() + k, 677 sameFace && i.intValue() + size < MAX_SIZE).parent(nbrLevel)); 678 if (k >= size) { 679 break; 680 } 681 } 682 } 683 684 // /////////////////////////////////////////////////////////////////// 685 // Low-level methods. 686 687 /** 688 * Return a leaf cell given its cube face (range 0..5) and i- and 689 * j-coordinates (see s2.h). 690 */ fromFaceIJ(int face, int i, int j)691 public static S2CellId fromFaceIJ(int face, int i, int j) { 692 // Optimization notes: 693 // - Non-overlapping bit fields can be combined with either "+" or "|". 694 // Generally "+" seems to produce better code, but not always. 695 696 // gcc doesn't have very good code generation for 64-bit operations. 697 // We optimize this by computing the result as two 32-bit integers 698 // and combining them at the end. Declaring the result as an array 699 // rather than local variables helps the compiler to do a better job 700 // of register allocation as well. Note that the two 32-bits halves 701 // get shifted one bit to the left when they are combined. 702 long n[] = {0, face << (POS_BITS - 33)}; 703 704 // Alternating faces have opposite Hilbert curve orientations; this 705 // is necessary in order for all faces to have a right-handed 706 // coordinate system. 707 int bits = (face & SWAP_MASK); 708 709 // Each iteration maps 4 bits of "i" and "j" into 8 bits of the Hilbert 710 // curve position. The lookup table transforms a 10-bit key of the form 711 // "iiiijjjjoo" to a 10-bit value of the form "ppppppppoo", where the 712 // letters [ijpo] denote bits of "i", "j", Hilbert curve position, and 713 // Hilbert curve orientation respectively. 714 715 for (int k = 7; k >= 0; --k) { 716 bits = getBits(n, i, j, k, bits); 717 } 718 719 S2CellId s = new S2CellId((((n[1] << 32) + n[0]) << 1) + 1); 720 return s; 721 } 722 getBits(long[] n, int i, int j, int k, int bits)723 private static int getBits(long[] n, int i, int j, int k, int bits) { 724 final int mask = (1 << LOOKUP_BITS) - 1; 725 bits += (((i >> (k * LOOKUP_BITS)) & mask) << (LOOKUP_BITS + 2)); 726 bits += (((j >> (k * LOOKUP_BITS)) & mask) << 2); 727 bits = LOOKUP_POS[bits]; 728 n[k >> 2] |= ((((long) bits) >> 2) << ((k & 3) * 2 * LOOKUP_BITS)); 729 bits &= (SWAP_MASK | INVERT_MASK); 730 return bits; 731 } 732 733 734 /** 735 * Return the (face, i, j) coordinates for the leaf cell corresponding to this 736 * cell id. Since cells are represented by the Hilbert curve position at the 737 * center of the cell, the returned (i,j) for non-leaf cells will be a leaf 738 * cell adjacent to the cell center. If "orientation" is non-NULL, also return 739 * the Hilbert curve orientation for the current cell. 740 */ toFaceIJOrientation(MutableInteger pi, MutableInteger pj, MutableInteger orientation)741 public int toFaceIJOrientation(MutableInteger pi, MutableInteger pj, 742 MutableInteger orientation) { 743 // System.out.println("Entering toFaceIjorientation"); 744 int face = this.face(); 745 int bits = (face & SWAP_MASK); 746 747 // System.out.println("face = " + face + " bits = " + bits); 748 749 // Each iteration maps 8 bits of the Hilbert curve position into 750 // 4 bits of "i" and "j". The lookup table transforms a key of the 751 // form "ppppppppoo" to a value of the form "iiiijjjjoo", where the 752 // letters [ijpo] represents bits of "i", "j", the Hilbert curve 753 // position, and the Hilbert curve orientation respectively. 754 // 755 // On the first iteration we need to be careful to clear out the bits 756 // representing the cube face. 757 for (int k = 7; k >= 0; --k) { 758 bits = getBits1(pi, pj, k, bits); 759 // System.out.println("pi = " + pi + " pj= " + pj + " bits = " + bits); 760 } 761 762 if (orientation != null) { 763 // The position of a non-leaf cell at level "n" consists of a prefix of 764 // 2*n bits that identifies the cell, followed by a suffix of 765 // 2*(MAX_LEVEL-n)+1 bits of the form 10*. If n==MAX_LEVEL, the suffix is 766 // just "1" and has no effect. Otherwise, it consists of "10", followed 767 // by (MAX_LEVEL-n-1) repetitions of "00", followed by "0". The "10" has 768 // no effect, while each occurrence of "00" has the effect of reversing 769 // the kSwapMask bit. 770 // assert (S2.POS_TO_ORIENTATION[2] == 0); 771 // assert (S2.POS_TO_ORIENTATION[0] == S2.SWAP_MASK); 772 if ((lowestOnBit() & 0x1111111111111110L) != 0) { 773 bits ^= S2.SWAP_MASK; 774 } 775 orientation.setValue(bits); 776 } 777 return face; 778 } 779 getBits1(MutableInteger i, MutableInteger j, int k, int bits)780 private int getBits1(MutableInteger i, MutableInteger j, int k, int bits) { 781 final int nbits = (k == 7) ? (MAX_LEVEL - 7 * LOOKUP_BITS) : LOOKUP_BITS; 782 783 bits += (((int) (id >>> (k * 2 * LOOKUP_BITS + 1)) & 784 ((1 << (2 * nbits)) - 1))) << 2; 785 /* 786 * System.out.println("id is: " + id_); System.out.println("bits is " + 787 * bits); System.out.println("lookup_ij[bits] is " + lookup_ij[bits]); 788 */ 789 bits = LOOKUP_IJ[bits]; 790 i.setValue(i.intValue() 791 + ((bits >> (LOOKUP_BITS + 2)) << (k * LOOKUP_BITS))); 792 /* 793 * System.out.println("left is " + ((bits >> 2) & ((1 << kLookupBits) - 794 * 1))); System.out.println("right is " + (k * kLookupBits)); 795 * System.out.println("j is: " + j.intValue()); System.out.println("addition 796 * is: " + ((((bits >> 2) & ((1 << kLookupBits) - 1))) << (k * 797 * kLookupBits))); 798 */ 799 j.setValue(j.intValue() 800 + ((((bits >> 2) & ((1 << LOOKUP_BITS) - 1))) << (k * LOOKUP_BITS))); 801 bits &= (SWAP_MASK | INVERT_MASK); 802 return bits; 803 } 804 805 /** Return the lowest-numbered bit that is on for cells at the given level. */ lowestOnBit()806 public long lowestOnBit() { 807 return id & -id; 808 } 809 810 /** 811 * Return the lowest-numbered bit that is on for this cell id, which is equal 812 * to (uint64(1) << (2 * (MAX_LEVEL - level))). So for example, a.lsb() <= 813 * b.lsb() if and only if a.level() >= b.level(), but the first test is more 814 * efficient. 815 */ lowestOnBitForLevel(int level)816 public static long lowestOnBitForLevel(int level) { 817 return 1L << (2 * (MAX_LEVEL - level)); 818 } 819 820 821 /** 822 * Return the i- or j-index of the leaf cell containing the given s- or 823 * t-value. 824 */ stToIJ(double s)825 private static int stToIJ(double s) { 826 // Converting from floating-point to integers via static_cast is very slow 827 // on Intel processors because it requires changing the rounding mode. 828 // Rounding to the nearest integer using FastIntRound() is much faster. 829 830 final int m = MAX_SIZE / 2; // scaling multiplier 831 return (int) Math 832 .max(0, Math.min(2 * m - 1, Math.round(m * s + (m - 0.5)))); 833 } 834 835 /** 836 * Convert (face, si, ti) coordinates (see s2.h) to a direction vector (not 837 * necessarily unit length). 838 */ faceSiTiToXYZ(int face, int si, int ti)839 private static S2Point faceSiTiToXYZ(int face, int si, int ti) { 840 final double kScale = 1.0 / MAX_SIZE; 841 double u = S2Projections.stToUV(kScale * si); 842 double v = S2Projections.stToUV(kScale * ti); 843 return S2Projections.faceUvToXyz(face, u, v); 844 } 845 846 /** 847 * Given (i, j) coordinates that may be out of bounds, normalize them by 848 * returning the corresponding neighbor cell on an adjacent face. 849 */ fromFaceIJWrap(int face, int i, int j)850 private static S2CellId fromFaceIJWrap(int face, int i, int j) { 851 // Convert i and j to the coordinates of a leaf cell just beyond the 852 // boundary of this face. This prevents 32-bit overflow in the case 853 // of finding the neighbors of a face cell, and also means that we 854 // don't need to worry about the distinction between (s,t) and (u,v). 855 i = Math.max(-1, Math.min(MAX_SIZE, i)); 856 j = Math.max(-1, Math.min(MAX_SIZE, j)); 857 858 // Find the (s,t) coordinates corresponding to (i,j). At least one 859 // of these coordinates will be just outside the range [0, 1]. 860 final double kScale = 1.0 / MAX_SIZE; 861 double s = kScale * ((i << 1) + 1 - MAX_SIZE); 862 double t = kScale * ((j << 1) + 1 - MAX_SIZE); 863 864 // Find the leaf cell coordinates on the adjacent face, and convert 865 // them to a cell id at the appropriate level. 866 S2Point p = S2Projections.faceUvToXyz(face, s, t); 867 face = S2Projections.xyzToFace(p); 868 R2Vector st = S2Projections.validFaceXyzToUv(face, p); 869 return fromFaceIJ(face, stToIJ(st.x()), stToIJ(st.y())); 870 } 871 872 /** 873 * Public helper function that calls FromFaceIJ if sameFace is true, or 874 * FromFaceIJWrap if sameFace is false. 875 */ fromFaceIJSame(int face, int i, int j, boolean sameFace)876 public static S2CellId fromFaceIJSame(int face, int i, int j, 877 boolean sameFace) { 878 if (sameFace) { 879 return S2CellId.fromFaceIJ(face, i, j); 880 } else { 881 return S2CellId.fromFaceIJWrap(face, i, j); 882 } 883 } 884 885 @Override equals(Object that)886 public boolean equals(Object that) { 887 if (!(that instanceof S2CellId)) { 888 return false; 889 } 890 S2CellId x = (S2CellId) that; 891 return id() == x.id(); 892 } 893 894 /** 895 * Returns true if x1 < x2, when both values are treated as unsigned. 896 */ unsignedLongLessThan(long x1, long x2)897 public static boolean unsignedLongLessThan(long x1, long x2) { 898 return (x1 + Long.MIN_VALUE) < (x2 + Long.MIN_VALUE); 899 } 900 901 /** 902 * Returns true if x1 > x2, when both values are treated as unsigned. 903 */ unsignedLongGreaterThan(long x1, long x2)904 public static boolean unsignedLongGreaterThan(long x1, long x2) { 905 return (x1 + Long.MIN_VALUE) > (x2 + Long.MIN_VALUE); 906 } 907 lessThan(S2CellId x)908 public boolean lessThan(S2CellId x) { 909 return unsignedLongLessThan(id, x.id); 910 } 911 greaterThan(S2CellId x)912 public boolean greaterThan(S2CellId x) { 913 return unsignedLongGreaterThan(id, x.id); 914 } 915 lessOrEquals(S2CellId x)916 public boolean lessOrEquals(S2CellId x) { 917 return unsignedLongLessThan(id, x.id) || id == x.id; 918 } 919 greaterOrEquals(S2CellId x)920 public boolean greaterOrEquals(S2CellId x) { 921 return unsignedLongGreaterThan(id, x.id) || id == x.id; 922 } 923 924 @Override hashCode()925 public int hashCode() { 926 return (int) ((id >>> 32) + id); 927 } 928 929 930 @Override toString()931 public String toString() { 932 return "(face=" + face() + ", pos=" + Long.toHexString(pos()) + ", level=" 933 + level() + ")"; 934 } 935 initLookupCell(int level, int i, int j, int origOrientation, int pos, int orientation)936 private static void initLookupCell(int level, int i, int j, 937 int origOrientation, int pos, int orientation) { 938 if (level == LOOKUP_BITS) { 939 int ij = (i << LOOKUP_BITS) + j; 940 LOOKUP_POS[(ij << 2) + origOrientation] = (pos << 2) + orientation; 941 LOOKUP_IJ[(pos << 2) + origOrientation] = (ij << 2) + orientation; 942 } else { 943 level++; 944 i <<= 1; 945 j <<= 1; 946 pos <<= 2; 947 // Initialize each sub-cell recursively. 948 for (int subPos = 0; subPos < 4; subPos++) { 949 int ij = S2.posToIJ(orientation, subPos); 950 int orientationMask = S2.posToOrientation(subPos); 951 initLookupCell(level, i + (ij >>> 1), j + (ij & 1), origOrientation, 952 pos + subPos, orientation ^ orientationMask); 953 } 954 } 955 } 956 957 @Override compareTo(S2CellId that)958 public int compareTo(S2CellId that) { 959 return unsignedLongLessThan(this.id, that.id) ? -1 : 960 unsignedLongGreaterThan(this.id, that.id) ? 1 : 0; 961 } 962 } 963