1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#include "s2cellid.h"
4
5#include <pthread.h>
6
7#include <algorithm>
8using std::min;
9using std::max;
10using std::swap;
11using std::reverse;
12
13#include <iomanip>
14using std::setprecision;
15
16#include <vector>
17using std::vector;
18
19
20#include "base/integral_types.h"
21#include "base/logging.h"
22#include "strings/strutil.h"
23#include "s2.h"
24#include "s2latlng.h"
25#include "util/math/mathutil.h"
26#include "util/math/vector2-inl.h"
27
28// The following lookup tables are used to convert efficiently between an
29// (i,j) cell index and the corresponding position along the Hilbert curve.
30// "lookup_pos" maps 4 bits of "i", 4 bits of "j", and 2 bits representing the
31// orientation of the current cell into 8 bits representing the order in which
32// that subcell is visited by the Hilbert curve, plus 2 bits indicating the
33// new orientation of the Hilbert curve within that subcell. (Cell
34// orientations are represented as combination of kSwapMask and kInvertMask.)
35//
36// "lookup_ij" is an inverted table used for mapping in the opposite
37// direction.
38//
39// We also experimented with looking up 16 bits at a time (14 bits of position
40// plus 2 of orientation) but found that smaller lookup tables gave better
41// performance. (2KB fits easily in the primary cache.)
42
43
44// Values for these constants are *declared* in the *.h file. Even though
45// the declaration specifies a value for the constant, that declaration
46// is not a *definition* of storage for the value. Because the values are
47// supplied in the declaration, we don't need the values here. Failing to
48// define storage causes link errors for any code that tries to take the
49// address of one of these values.
50int const S2CellId::kFaceBits;
51int const S2CellId::kNumFaces;
52int const S2CellId::kMaxLevel;
53int const S2CellId::kPosBits;
54int const S2CellId::kMaxSize;
55
56static int const kLookupBits = 4;
57static int const kSwapMask = 0x01;
58static int const kInvertMask = 0x02;
59
60static uint16 lookup_pos[1 << (2 * kLookupBits + 2)];
61static uint16 lookup_ij[1 << (2 * kLookupBits + 2)];
62
63static void InitLookupCell(int level, int i, int j, int orig_orientation,
64 int pos, int orientation) {
65 if (level == kLookupBits) {
66 int ij = (i << kLookupBits) + j;
67 lookup_pos[(ij << 2) + orig_orientation] = (pos << 2) + orientation;
68 lookup_ij[(pos << 2) + orig_orientation] = (ij << 2) + orientation;
69 } else {
70 level++;
71 i <<= 1;
72 j <<= 1;
73 pos <<= 2;
74 int const* r = S2::kPosToIJ[orientation];
75 InitLookupCell(level, i + (r[0] >> 1), j + (r[0] & 1), orig_orientation,
76 pos, orientation ^ S2::kPosToOrientation[0]);
77 InitLookupCell(level, i + (r[1] >> 1), j + (r[1] & 1), orig_orientation,
78 pos + 1, orientation ^ S2::kPosToOrientation[1]);
79 InitLookupCell(level, i + (r[2] >> 1), j + (r[2] & 1), orig_orientation,
80 pos + 2, orientation ^ S2::kPosToOrientation[2]);
81 InitLookupCell(level, i + (r[3] >> 1), j + (r[3] & 1), orig_orientation,
82 pos + 3, orientation ^ S2::kPosToOrientation[3]);
83 }
84}
85
86static void Init() {
87 InitLookupCell(0, 0, 0, 0, 0, 0);
88 InitLookupCell(0, 0, 0, kSwapMask, 0, kSwapMask);
89 InitLookupCell(0, 0, 0, kInvertMask, 0, kInvertMask);
90 InitLookupCell(0, 0, 0, kSwapMask|kInvertMask, 0, kSwapMask|kInvertMask);
91}
92
93static pthread_once_t init_once = PTHREAD_ONCE_INIT;
94inline static void MaybeInit() {
95 pthread_once(&init_once, Init);
96}
97
98int S2CellId::level() const {
99 // Fast path for leaf cells.
100 if (is_leaf()) return kMaxLevel;
101
102 uint32 x = static_cast<uint32>(id_);
103 int level = -1;
104 if (x != 0) {
105 level += 16;
106 } else {
107 x = static_cast<uint32>(id_ >> 32);
108 }
109 // We only need to look at even-numbered bits to determine the
110 // level of a valid cell id.
111 x &= -x; // Get lowest bit.
112 if (x & 0x00005555) level += 8;
113 if (x & 0x00550055) level += 4;
114 if (x & 0x05050505) level += 2;
115 if (x & 0x11111111) level += 1;
116 DCHECK_GE(level, 0);
117 DCHECK_LE(level, kMaxLevel);
118 return level;
119}
120
121S2CellId S2CellId::advance(int64 steps) const {
122 if (steps == 0) return *this;
123
124 // We clamp the number of steps if necessary to ensure that we do not
125 // advance past the End() or before the Begin() of this level. Note that
126 // min_steps and max_steps always fit in a signed 64-bit integer.
127
128 int step_shift = 2 * (kMaxLevel - level()) + 1;
129 if (steps < 0) {
130 int64 min_steps = -static_cast<int64>(id_ >> step_shift);
131 if (steps < min_steps) steps = min_steps;
132 } else {
133 int64 max_steps = (kWrapOffset + lsb() - id_) >> step_shift;
134 if (steps > max_steps) steps = max_steps;
135 }
136 return S2CellId(id_ + (steps << step_shift));
137}
138
139S2CellId S2CellId::advance_wrap(int64 steps) const {
140 DCHECK(is_valid());
141 if (steps == 0) return *this;
142
143 int step_shift = 2 * (kMaxLevel - level()) + 1;
144 if (steps < 0) {
145 int64 min_steps = -static_cast<int64>(id_ >> step_shift);
146 if (steps < min_steps) {
147 int64 step_wrap = kWrapOffset >> step_shift;
148 steps %= step_wrap;
149 if (steps < min_steps) steps += step_wrap;
150 }
151 } else {
152 // Unlike advance(), we don't want to return End(level).
153 int64 max_steps = (kWrapOffset - id_) >> step_shift;
154 if (steps > max_steps) {
155 int64 step_wrap = kWrapOffset >> step_shift;
156 steps %= step_wrap;
157 if (steps > max_steps) steps -= step_wrap;
158 }
159 }
160 return S2CellId(id_ + (steps << step_shift));
161}
162
163S2CellId S2CellId::FromFacePosLevel(int face, uint64 pos, int level) {
164 S2CellId cell((static_cast<uint64>(face) << kPosBits) + (pos | 1));
165 return cell.parent(level);
166}
167
168string S2CellId::ToToken() const {
169 // Simple implementation: convert the id to hex and strip trailing zeros.
170 // Using hex has the advantage that the tokens are case-insensitive, all
171 // characters are alphanumeric, no characters require any special escaping
172 // in Mustang queries, and it's easy to compare cell tokens against the
173 // feature ids of the corresponding features.
174 //
175 // Using base 64 would produce slightly shorter tokens, but for typical cell
176 // sizes used during indexing (up to level 15 or so) the average savings
177 // would be less than 2 bytes per cell which doesn't seem worth it.
178
179 char digits[17];
180 FastHex64ToBuffer(id_, digits);
181 for (int len = 16; len > 0; --len) {
182 if (digits[len-1] != '0') {
183 return string(digits, len);
184 }
185 }
186 return "X"; // Invalid hex string.
187}
188
189S2CellId S2CellId::FromToken(string const& token) {
190 if (token.size() > 16) return S2CellId::None();
191 char digits[17] = "0000000000000000";
192 memcpy(digits, token.data(), token.size());
193 return S2CellId(ParseLeadingHex64Value(digits, 0));
194}
195
196inline int S2CellId::STtoIJ(double s) {
197 // Converting from floating-point to integers via static_cast is very slow
198 // on Intel processors because it requires changing the rounding mode.
199 // Rounding to the nearest integer using FastIntRound() is much faster.
200
201 return max(0, min(kMaxSize - 1, MathUtil::FastIntRound(kMaxSize * s - 0.5)));
202}
203
204
205S2CellId S2CellId::FromFaceIJ(int face, int i, int j) {
206 // Initialization if not done yet
207 MaybeInit();
208
209 // Optimization notes:
210 // - Non-overlapping bit fields can be combined with either "+" or "|".
211 // Generally "+" seems to produce better code, but not always.
212
213 // gcc doesn't have very good code generation for 64-bit operations.
214 // We optimize this by computing the result as two 32-bit integers
215 // and combining them at the end. Declaring the result as an array
216 // rather than local variables helps the compiler to do a better job
217 // of register allocation as well. Note that the two 32-bits halves
218 // get shifted one bit to the left when they are combined.
219 uint32 n[2] = { 0, static_cast<uint32>(face << (kPosBits - 33)) };
220
221 // Alternating faces have opposite Hilbert curve orientations; this
222 // is necessary in order for all faces to have a right-handed
223 // coordinate system.
224 int bits = (face & kSwapMask);
225
226 // Each iteration maps 4 bits of "i" and "j" into 8 bits of the Hilbert
227 // curve position. The lookup table transforms a 10-bit key of the form
228 // "iiiijjjjoo" to a 10-bit value of the form "ppppppppoo", where the
229 // letters [ijpo] denote bits of "i", "j", Hilbert curve position, and
230 // Hilbert curve orientation respectively.
231#define GET_BITS(k) do { \
232 int const mask = (1 << kLookupBits) - 1; \
233 bits += ((i >> (k * kLookupBits)) & mask) << (kLookupBits + 2); \
234 bits += ((j >> (k * kLookupBits)) & mask) << 2; \
235 bits = lookup_pos[bits]; \
236 n[k >> 2] |= (bits >> 2) << ((k & 3) * 2 * kLookupBits); \
237 bits &= (kSwapMask | kInvertMask); \
238 } while (0)
239
240 GET_BITS(7);
241 GET_BITS(6);
242 GET_BITS(5);
243 GET_BITS(4);
244 GET_BITS(3);
245 GET_BITS(2);
246 GET_BITS(1);
247 GET_BITS(0);
248#undef GET_BITS
249
250 return S2CellId(((static_cast<uint64>(n[1]) << 32) + n[0]) * 2 + 1);
251}
252
253S2CellId S2CellId::FromPoint(S2Point const& p) {
254 double u, v;
255 int face = S2::XYZtoFaceUV(p, &u, &v);
256 int i = STtoIJ(S2::UVtoST(u));
257 int j = STtoIJ(S2::UVtoST(v));
258 return FromFaceIJ(face, i, j);
259}
260
261S2CellId S2CellId::FromLatLng(S2LatLng const& ll) {
262 return FromPoint(ll.ToPoint());
263}
264
265int S2CellId::ToFaceIJOrientation(int* pi, int* pj, int* orientation) const {
266 // Initialization if not done yet
267 MaybeInit();
268
269 int i = 0, j = 0;
270 int face = this->face();
271 int bits = (face & kSwapMask);
272
273 // Each iteration maps 8 bits of the Hilbert curve position into
274 // 4 bits of "i" and "j". The lookup table transforms a key of the
275 // form "ppppppppoo" to a value of the form "iiiijjjjoo", where the
276 // letters [ijpo] represents bits of "i", "j", the Hilbert curve
277 // position, and the Hilbert curve orientation respectively.
278 //
279 // On the first iteration we need to be careful to clear out the bits
280 // representing the cube face.
281#define GET_BITS(k) do { \
282 int const nbits = (k == 7) ? (kMaxLevel - 7 * kLookupBits) : kLookupBits; \
283 bits += (static_cast<int>(id_ >> (k * 2 * kLookupBits + 1)) \
284 & ((1 << (2 * nbits)) - 1)) << 2; \
285 bits = lookup_ij[bits]; \
286 i += (bits >> (kLookupBits + 2)) << (k * kLookupBits); \
287 j += ((bits >> 2) & ((1 << kLookupBits) - 1)) << (k * kLookupBits); \
288 bits &= (kSwapMask | kInvertMask); \
289 } while (0)
290
291 GET_BITS(7);
292 GET_BITS(6);
293 GET_BITS(5);
294 GET_BITS(4);
295 GET_BITS(3);
296 GET_BITS(2);
297 GET_BITS(1);
298 GET_BITS(0);
299#undef GET_BITS
300
301 *pi = i;
302 *pj = j;
303
304 if (orientation != NULL) {
305 // The position of a non-leaf cell at level "n" consists of a prefix of
306 // 2*n bits that identifies the cell, followed by a suffix of
307 // 2*(kMaxLevel-n)+1 bits of the form 10*. If n==kMaxLevel, the suffix is
308 // just "1" and has no effect. Otherwise, it consists of "10", followed
309 // by (kMaxLevel-n-1) repetitions of "00", followed by "0". The "10" has
310 // no effect, while each occurrence of "00" has the effect of reversing
311 // the kSwapMask bit.
312 DCHECK_EQ(0, S2::kPosToOrientation[2]);
313 DCHECK_EQ(S2::kSwapMask, S2::kPosToOrientation[0]);
314 if (lsb() & GG_ULONGLONG(0x1111111111111110)) {
315 bits ^= S2::kSwapMask;
316 }
317 *orientation = bits;
318 }
319 return face;
320}
321
322inline int S2CellId::GetCenterSiTi(int* psi, int* pti) const {
323 // First we compute the discrete (i,j) coordinates of a leaf cell contained
324 // within the given cell. Given that cells are represented by the Hilbert
325 // curve position corresponding at their center, it turns out that the cell
326 // returned by ToFaceIJOrientation is always one of two leaf cells closest
327 // to the center of the cell (unless the given cell is a leaf cell itself,
328 // in which case there is only one possibility).
329 //
330 // Given a cell of size s >= 2 (i.e. not a leaf cell), and letting (imin,
331 // jmin) be the coordinates of its lower left-hand corner, the leaf cell
332 // returned by ToFaceIJOrientation() is either (imin + s/2, jmin + s/2)
333 // (imin + s/2 - 1, jmin + s/2 - 1). The first case is the one we want.
334 // We can distinguish these two cases by looking at the low bit of "i" or
335 // "j". In the second case the low bit is one, unless s == 2 (i.e. the
336 // level just above leaf cells) in which case the low bit is zero.
337 //
338 // In the code below, the expression ((i ^ (int(id_) >> 2)) & 1) is true
339 // if we are in the second case described above.
340 int i, j;
341 int face = ToFaceIJOrientation(&i, &j, NULL);
342 int delta = is_leaf() ? 1 : ((i ^ (static_cast<int>(id_) >> 2)) & 1) ? 2 : 0;
343
344 // Note that (2 * {i,j} + delta) will never overflow a 32-bit integer.
345 *psi = 2 * i + delta;
346 *pti = 2 * j + delta;
347 return face;
348}
349
350S2Point S2CellId::ToPointRaw() const {
351 // This code would be slightly shorter if we called GetCenterUV(),
352 // but this method is heavily used and it's 25% faster to include
353 // the method inline.
354 int si, ti;
355 int face = GetCenterSiTi(&si, &ti);
356 return S2::FaceUVtoXYZ(face,
357 S2::STtoUV((0.5 / kMaxSize) * si),
358 S2::STtoUV((0.5 / kMaxSize) * ti));
359}
360
361S2LatLng S2CellId::ToLatLng() const {
362 return S2LatLng(ToPointRaw());
363}
364
365Vector2_d S2CellId::GetCenterST() const {
366 int si, ti;
367 GetCenterSiTi(&si, &ti);
368 return Vector2_d((0.5 / kMaxSize) * si, (0.5 / kMaxSize) * ti);
369}
370
371Vector2_d S2CellId::GetCenterUV() const {
372 int si, ti;
373 GetCenterSiTi(&si, &ti);
374 return Vector2_d(S2::STtoUV((0.5 / kMaxSize) * si),
375 S2::STtoUV((0.5 / kMaxSize) * ti));
376}
377
378S2CellId S2CellId::FromFaceIJWrap(int face, int i, int j) {
379 // Convert i and j to the coordinates of a leaf cell just beyond the
380 // boundary of this face. This prevents 32-bit overflow in the case
381 // of finding the neighbors of a face cell.
382 i = max(-1, min(kMaxSize, i));
383 j = max(-1, min(kMaxSize, j));
384
385 // Find the (u,v) coordinates corresponding to the center of cell (i,j).
386 // For our purposes it's sufficient to always use the linear projection
387 // from (s,t) to (u,v): u=2*s-1 and v=2*t-1.
388 static const double kScale = 1.0 / kMaxSize;
389 double u = kScale * ((i << 1) + 1 - kMaxSize);
390 double v = kScale * ((j << 1) + 1 - kMaxSize);
391
392 // Find the leaf cell coordinates on the adjacent face, and convert
393 // them to a cell id at the appropriate level. We convert from (u,v)
394 // back to (s,t) using s=0.5*(u+1), t=0.5*(v+1).
395 face = S2::XYZtoFaceUV(S2::FaceUVtoXYZ(face, u, v), &u, &v);
396 return FromFaceIJ(face, STtoIJ(0.5*(u+1)), STtoIJ(0.5*(v+1)));
397}
398
399inline S2CellId S2CellId::FromFaceIJSame(int face, int i, int j,
400 bool same_face) {
401 if (same_face)
402 return S2CellId::FromFaceIJ(face, i, j);
403 else
404 return S2CellId::FromFaceIJWrap(face, i, j);
405}
406
407void S2CellId::GetEdgeNeighbors(S2CellId neighbors[4]) const {
408 int i, j;
409 int level = this->level();
410 int size = GetSizeIJ(level);
411 int face = ToFaceIJOrientation(&i, &j, NULL);
412
413 // Edges 0, 1, 2, 3 are in the S, E, N, W directions.
414 neighbors[0] = FromFaceIJSame(face, i, j - size, j - size >= 0)
415 .parent(level);
416 neighbors[1] = FromFaceIJSame(face, i + size, j, i + size < kMaxSize)
417 .parent(level);
418 neighbors[2] = FromFaceIJSame(face, i, j + size, j + size < kMaxSize)
419 .parent(level);
420 neighbors[3] = FromFaceIJSame(face, i - size, j, i - size >= 0)
421 .parent(level);
422}
423
424void S2CellId::AppendVertexNeighbors(int level,
425 vector<S2CellId>* output) const {
426 // "level" must be strictly less than this cell's level so that we can
427 // determine which vertex this cell is closest to.
428 DCHECK_LT(level, this->level());
429 int i, j;
430 int face = ToFaceIJOrientation(&i, &j, NULL);
431
432 // Determine the i- and j-offsets to the closest neighboring cell in each
433 // direction. This involves looking at the next bit of "i" and "j" to
434 // determine which quadrant of this->parent(level) this cell lies in.
435 int halfsize = GetSizeIJ(level + 1);
436 int size = halfsize << 1;
437 bool isame, jsame;
438 int ioffset, joffset;
439 if (i & halfsize) {
440 ioffset = size;
441 isame = (i + size) < kMaxSize;
442 } else {
443 ioffset = -size;
444 isame = (i - size) >= 0;
445 }
446 if (j & halfsize) {
447 joffset = size;
448 jsame = (j + size) < kMaxSize;
449 } else {
450 joffset = -size;
451 jsame = (j - size) >= 0;
452 }
453
454 output->push_back(parent(level));
455 output->push_back(FromFaceIJSame(face, i + ioffset, j, isame).parent(level));
456 output->push_back(FromFaceIJSame(face, i, j + joffset, jsame).parent(level));
457 // If i- and j- edge neighbors are *both* on a different face, then this
458 // vertex only has three neighbors (it is one of the 8 cube vertices).
459 if (isame || jsame) {
460 output->push_back(FromFaceIJSame(face, i + ioffset, j + joffset,
461 isame && jsame).parent(level));
462 }
463}
464
465void S2CellId::AppendAllNeighbors(int nbr_level,
466 vector<S2CellId>* output) const {
467 int i, j;
468 int face = ToFaceIJOrientation(&i, &j, NULL);
469
470 // Find the coordinates of the lower left-hand leaf cell. We need to
471 // normalize (i,j) to a known position within the cell because nbr_level
472 // may be larger than this cell's level.
473 int size = GetSizeIJ();
474 i &= -size;
475 j &= -size;
476
477 int nbr_size = GetSizeIJ(nbr_level);
478 DCHECK_LE(nbr_size, size);
479
480 // We compute the N-S, E-W, and diagonal neighbors in one pass.
481 // The loop test is at the end of the loop to avoid 32-bit overflow.
482 for (int k = -nbr_size; ; k += nbr_size) {
483 bool same_face;
484 if (k < 0) {
485 same_face = (j + k >= 0);
486 } else if (k >= size) {
487 same_face = (j + k < kMaxSize);
488 } else {
489 same_face = true;
490 // North and South neighbors.
491 output->push_back(FromFaceIJSame(face, i + k, j - nbr_size,
492 j - size >= 0).parent(nbr_level));
493 output->push_back(FromFaceIJSame(face, i + k, j + size,
494 j + size < kMaxSize).parent(nbr_level));
495 }
496 // East, West, and Diagonal neighbors.
497 output->push_back(FromFaceIJSame(face, i - nbr_size, j + k,
498 same_face && i - size >= 0)
499 .parent(nbr_level));
500 output->push_back(FromFaceIJSame(face, i + size, j + k,
501 same_face && i + size < kMaxSize)
502 .parent(nbr_level));
503 if (k >= size) break;
504 }
505}
506
507string S2CellId::ToString() const {
508 if (!is_valid()) {
509 return StringPrintf("Invalid: %016llx", id());
510 }
511 string out = IntToString(face(), "%d/");
512 for (int current_level = 1; current_level <= level(); ++current_level) {
513 out += IntToString(child_position(current_level), "%d");
514 }
515 return out;
516}
517
518ostream& operator<<(ostream& os, S2CellId const& id) {
519 return os << id.ToString();
520}
521