1
// Copyright 2005 Google Inc. All Rights Reserved.
16
using std::setprecision;
22
#include "base/integral_types.h"
23
#include "base/logging.h"
24
#include "strings/strutil.h"
27
#include "util/math/mathutil.h"
28
#include "util/math/vector2-inl.h"
30
#include "mongo/base/init.h"
32
// The following lookup tables are used to convert efficiently between an
33
// (i,j) cell index and the corresponding position along the Hilbert curve.
34
// "lookup_pos" maps 4 bits of "i", 4 bits of "j", and 2 bits representing the
35
// orientation of the current cell into 8 bits representing the order in which
36
// that subcell is visited by the Hilbert curve, plus 2 bits indicating the
37
// new orientation of the Hilbert curve within that subcell. (Cell
38
// orientations are represented as combination of kSwapMask and kInvertMask.)
40
// "lookup_ij" is an inverted table used for mapping in the opposite
43
// We also experimented with looking up 16 bits at a time (14 bits of position
44
// plus 2 of orientation) but found that smaller lookup tables gave better
45
// performance. (2KB fits easily in the primary cache.)
48
// Values for these constants are *declared* in the *.h file. Even though
49
// the declaration specifies a value for the constant, that declaration
50
// is not a *definition* of storage for the value. Because the values are
51
// supplied in the declaration, we don't need the values here. Failing to
52
// define storage causes link errors for any code that tries to take the
53
// address of one of these values.
54
int const S2CellId::kFaceBits = 3;
55
int const S2CellId::kNumFaces = 6;
56
int const S2CellId::kMaxLevel = S2::kMaxCellLevel;
57
int const S2CellId::kPosBits = 2 * kMaxLevel + 1;
58
int const S2CellId::kMaxSize = 1 << kMaxLevel;
59
uint64 const S2CellId::kWrapOffset = uint64(kNumFaces) << kPosBits;
61
static int const kLookupBits = 4;
62
static int const kSwapMask = 0x01;
63
static int const kInvertMask = 0x02;
65
static uint16 lookup_pos[1 << (2 * kLookupBits + 2)];
66
static uint16 lookup_ij[1 << (2 * kLookupBits + 2)];
68
static void InitLookupCell(int level, int i, int j, int orig_orientation,
69
int pos, int orientation) {
70
if (level == kLookupBits) {
71
int ij = (i << kLookupBits) + j;
72
lookup_pos[(ij << 2) + orig_orientation] = (pos << 2) + orientation;
73
lookup_ij[(pos << 2) + orig_orientation] = (ij << 2) + orientation;
79
int const* r = S2::kPosToIJ[orientation];
80
InitLookupCell(level, i + (r[0] >> 1), j + (r[0] & 1), orig_orientation,
81
pos, orientation ^ S2::kPosToOrientation[0]);
82
InitLookupCell(level, i + (r[1] >> 1), j + (r[1] & 1), orig_orientation,
83
pos + 1, orientation ^ S2::kPosToOrientation[1]);
84
InitLookupCell(level, i + (r[2] >> 1), j + (r[2] & 1), orig_orientation,
85
pos + 2, orientation ^ S2::kPosToOrientation[2]);
86
InitLookupCell(level, i + (r[3] >> 1), j + (r[3] & 1), orig_orientation,
87
pos + 3, orientation ^ S2::kPosToOrientation[3]);
92
InitLookupCell(0, 0, 0, 0, 0, 0);
93
InitLookupCell(0, 0, 0, kSwapMask, 0, kSwapMask);
94
InitLookupCell(0, 0, 0, kInvertMask, 0, kInvertMask);
95
InitLookupCell(0, 0, 0, kSwapMask|kInvertMask, 0, kSwapMask|kInvertMask);
98
MONGO_INITIALIZER(S2CellIdInit)(mongo::InitializerContext *context) {
100
return mongo::Status::OK();
103
int S2CellId::level() const {
104
// Fast path for leaf cells.
105
if (is_leaf()) return kMaxLevel;
107
uint32 x = static_cast<uint32>(id_);
112
x = static_cast<uint32>(id_ >> 32);
114
// We only need to look at even-numbered bits to determine the
115
// level of a valid cell id.
116
#pragma warning(push)
117
#pragma warning( disable: 4146 )
118
x &= -x; // Get lowest bit.
120
if (x & 0x00005555) level += 8;
121
if (x & 0x00550055) level += 4;
122
if (x & 0x05050505) level += 2;
123
if (x & 0x11111111) level += 1;
125
DCHECK_LE(level, kMaxLevel);
129
S2CellId S2CellId::advance(int64 steps) const {
130
if (steps == 0) return *this;
132
// We clamp the number of steps if necessary to ensure that we do not
133
// advance past the End() or before the Begin() of this level. Note that
134
// min_steps and max_steps always fit in a signed 64-bit integer.
136
int step_shift = 2 * (kMaxLevel - level()) + 1;
138
int64 min_steps = -static_cast<int64>(id_ >> step_shift);
139
if (steps < min_steps) steps = min_steps;
141
int64 max_steps = (kWrapOffset + lsb() - id_) >> step_shift;
142
if (steps > max_steps) steps = max_steps;
144
return S2CellId(id_ + (steps << step_shift));
147
S2CellId S2CellId::advance_wrap(int64 steps) const {
149
if (steps == 0) return *this;
151
int step_shift = 2 * (kMaxLevel - level()) + 1;
153
int64 min_steps = -static_cast<int64>(id_ >> step_shift);
154
if (steps < min_steps) {
155
int64 step_wrap = kWrapOffset >> step_shift;
157
if (steps < min_steps) steps += step_wrap;
160
// Unlike advance(), we don't want to return End(level).
161
int64 max_steps = (kWrapOffset - id_) >> step_shift;
162
if (steps > max_steps) {
163
int64 step_wrap = kWrapOffset >> step_shift;
165
if (steps > max_steps) steps -= step_wrap;
168
return S2CellId(id_ + (steps << step_shift));
171
S2CellId S2CellId::FromFacePosLevel(int face, uint64 pos, int level) {
172
S2CellId cell((static_cast<uint64>(face) << kPosBits) + (pos | 1));
173
return cell.parent(level);
176
string S2CellId::ToToken() const {
177
// Simple implementation: convert the id to hex and strip trailing zeros.
178
// Using hex has the advantage that the tokens are case-insensitive, all
179
// characters are alphanumeric, no characters require any special escaping
180
// in Mustang queries, and it's easy to compare cell tokens against the
181
// feature ids of the corresponding features.
183
// Using base 64 would produce slightly shorter tokens, but for typical cell
184
// sizes used during indexing (up to level 15 or so) the average savings
185
// would be less than 2 bytes per cell which doesn't seem worth it.
188
FastHex64ToBuffer(id_, digits);
189
for (int len = 16; len > 0; --len) {
190
if (digits[len-1] != '0') {
191
return string(digits, len);
194
return "X"; // Invalid hex string.
197
S2CellId S2CellId::FromToken(string const& token) {
198
if (token.size() > 16) return S2CellId::None();
199
char digits[17] = "0000000000000000";
200
memcpy(digits, token.data(), token.size());
201
return S2CellId(ParseLeadingHex64Value(digits, 0));
204
inline int S2CellId::STtoIJ(double s) {
205
// Converting from floating-point to integers via static_cast is very slow
206
// on Intel processors because it requires changing the rounding mode.
207
// Rounding to the nearest integer using FastIntRound() is much faster.
209
return max(0, min(kMaxSize - 1, MathUtil::FastIntRound(kMaxSize * s - 0.5)));
213
S2CellId S2CellId::FromFaceIJ(int face, int i, int j) {
214
// Optimization notes:
215
// - Non-overlapping bit fields can be combined with either "+" or "|".
216
// Generally "+" seems to produce better code, but not always.
218
// gcc doesn't have very good code generation for 64-bit operations.
219
// We optimize this by computing the result as two 32-bit integers
220
// and combining them at the end. Declaring the result as an array
221
// rather than local variables helps the compiler to do a better job
222
// of register allocation as well. Note that the two 32-bits halves
223
// get shifted one bit to the left when they are combined.
226
n[1] = face << (kPosBits - 33);
228
// Alternating faces have opposite Hilbert curve orientations; this
229
// is necessary in order for all faces to have a right-handed
230
// coordinate system.
231
int bits = (face & kSwapMask);
233
// Each iteration maps 4 bits of "i" and "j" into 8 bits of the Hilbert
234
// curve position. The lookup table transforms a 10-bit key of the form
235
// "iiiijjjjoo" to a 10-bit value of the form "ppppppppoo", where the
236
// letters [ijpo] denote bits of "i", "j", Hilbert curve position, and
237
// Hilbert curve orientation respectively.
238
#define GET_BITS(k) do { \
239
int const mask = (1 << kLookupBits) - 1; \
240
bits += ((i >> (k * kLookupBits)) & mask) << (kLookupBits + 2); \
241
bits += ((j >> (k * kLookupBits)) & mask) << 2; \
242
bits = lookup_pos[bits]; \
243
n[k >> 2] |= (bits >> 2) << ((k & 3) * 2 * kLookupBits); \
244
bits &= (kSwapMask | kInvertMask); \
257
return S2CellId(((static_cast<uint64>(n[1]) << 32) + n[0]) * 2 + 1);
260
S2CellId S2CellId::FromPoint(S2Point const& p) {
262
int face = S2::XYZtoFaceUV(p, &u, &v);
263
int i = STtoIJ(S2::UVtoST(u));
264
int j = STtoIJ(S2::UVtoST(v));
265
return FromFaceIJ(face, i, j);
268
S2CellId S2CellId::FromLatLng(S2LatLng const& ll) {
269
return FromPoint(ll.ToPoint());
272
int S2CellId::ToFaceIJOrientation(int* pi, int* pj, int* orientation) const {
274
int face = this->face();
275
int bits = (face & kSwapMask);
277
// Each iteration maps 8 bits of the Hilbert curve position into
278
// 4 bits of "i" and "j". The lookup table transforms a key of the
279
// form "ppppppppoo" to a value of the form "iiiijjjjoo", where the
280
// letters [ijpo] represents bits of "i", "j", the Hilbert curve
281
// position, and the Hilbert curve orientation respectively.
283
// On the first iteration we need to be careful to clear out the bits
284
// representing the cube face.
285
#define GET_BITS(k) do { \
286
int const nbits = (k == 7) ? (kMaxLevel - 7 * kLookupBits) : kLookupBits; \
287
bits += (static_cast<int>(id_ >> (k * 2 * kLookupBits + 1)) \
288
& ((1 << (2 * nbits)) - 1)) << 2; \
289
bits = lookup_ij[bits]; \
290
i += (bits >> (kLookupBits + 2)) << (k * kLookupBits); \
291
j += ((bits >> 2) & ((1 << kLookupBits) - 1)) << (k * kLookupBits); \
292
bits &= (kSwapMask | kInvertMask); \
308
if (orientation != NULL) {
309
// The position of a non-leaf cell at level "n" consists of a prefix of
310
// 2*n bits that identifies the cell, followed by a suffix of
311
// 2*(kMaxLevel-n)+1 bits of the form 10*. If n==kMaxLevel, the suffix is
312
// just "1" and has no effect. Otherwise, it consists of "10", followed
313
// by (kMaxLevel-n-1) repetitions of "00", followed by "0". The "10" has
314
// no effect, while each occurrence of "00" has the effect of reversing
315
// the kSwapMask bit.
316
DCHECK_EQ(0, S2::kPosToOrientation[2]);
317
DCHECK_EQ(S2::kSwapMask, S2::kPosToOrientation[0]);
318
if (lsb() & GG_ULONGLONG(0x1111111111111110)) {
319
bits ^= S2::kSwapMask;
326
inline int S2CellId::GetCenterSiTi(int* psi, int* pti) const {
327
// First we compute the discrete (i,j) coordinates of a leaf cell contained
328
// within the given cell. Given that cells are represented by the Hilbert
329
// curve position corresponding at their center, it turns out that the cell
330
// returned by ToFaceIJOrientation is always one of two leaf cells closest
331
// to the center of the cell (unless the given cell is a leaf cell itself,
332
// in which case there is only one possibility).
334
// Given a cell of size s >= 2 (i.e. not a leaf cell), and letting (imin,
335
// jmin) be the coordinates of its lower left-hand corner, the leaf cell
336
// returned by ToFaceIJOrientation() is either (imin + s/2, jmin + s/2)
337
// (imin + s/2 - 1, jmin + s/2 - 1). The first case is the one we want.
338
// We can distinguish these two cases by looking at the low bit of "i" or
339
// "j". In the second case the low bit is one, unless s == 2 (i.e. the
340
// level just above leaf cells) in which case the low bit is zero.
342
// In the code below, the expression ((i ^ (int(id_) >> 2)) & 1) is true
343
// if we are in the second case described above.
345
int face = ToFaceIJOrientation(&i, &j, NULL);
346
int delta = is_leaf() ? 1 : ((i ^ (static_cast<int>(id_) >> 2)) & 1) ? 2 : 0;
348
// Note that (2 * {i,j} + delta) will never overflow a 32-bit integer.
349
*psi = 2 * i + delta;
350
*pti = 2 * j + delta;
354
S2Point S2CellId::ToPointRaw() const {
355
// This code would be slightly shorter if we called GetCenterUV(),
356
// but this method is heavily used and it's 25% faster to include
357
// the method inline.
359
int face = GetCenterSiTi(&si, &ti);
360
return S2::FaceUVtoXYZ(face,
361
S2::STtoUV((0.5 / kMaxSize) * si),
362
S2::STtoUV((0.5 / kMaxSize) * ti));
365
S2LatLng S2CellId::ToLatLng() const {
366
return S2LatLng(ToPointRaw());
369
Vector2_d S2CellId::GetCenterST() const {
371
GetCenterSiTi(&si, &ti);
372
return Vector2_d((0.5 / kMaxSize) * si, (0.5 / kMaxSize) * ti);
375
Vector2_d S2CellId::GetCenterUV() const {
377
GetCenterSiTi(&si, &ti);
378
return Vector2_d(S2::STtoUV((0.5 / kMaxSize) * si),
379
S2::STtoUV((0.5 / kMaxSize) * ti));
382
S2CellId S2CellId::FromFaceIJWrap(int face, int i, int j) {
383
// Convert i and j to the coordinates of a leaf cell just beyond the
384
// boundary of this face. This prevents 32-bit overflow in the case
385
// of finding the neighbors of a face cell.
386
i = max(-1, min(kMaxSize, i));
387
j = max(-1, min(kMaxSize, j));
389
// Find the (u,v) coordinates corresponding to the center of cell (i,j).
390
// For our purposes it's sufficient to always use the linear projection
391
// from (s,t) to (u,v): u=2*s-1 and v=2*t-1.
392
static const double kScale = 1.0 / kMaxSize;
393
double u = kScale * ((i << 1) + 1 - kMaxSize);
394
double v = kScale * ((j << 1) + 1 - kMaxSize);
396
// Find the leaf cell coordinates on the adjacent face, and convert
397
// them to a cell id at the appropriate level. We convert from (u,v)
398
// back to (s,t) using s=0.5*(u+1), t=0.5*(v+1).
399
face = S2::XYZtoFaceUV(S2::FaceUVtoXYZ(face, u, v), &u, &v);
400
return FromFaceIJ(face, STtoIJ(0.5*(u+1)), STtoIJ(0.5*(v+1)));
403
inline S2CellId S2CellId::FromFaceIJSame(int face, int i, int j,
406
return S2CellId::FromFaceIJ(face, i, j);
408
return S2CellId::FromFaceIJWrap(face, i, j);
411
void S2CellId::GetEdgeNeighbors(S2CellId neighbors[4]) const {
413
int level = this->level();
414
int size = GetSizeIJ(level);
415
int face = ToFaceIJOrientation(&i, &j, NULL);
417
// Edges 0, 1, 2, 3 are in the S, E, N, W directions.
418
neighbors[0] = FromFaceIJSame(face, i, j - size, j - size >= 0)
420
neighbors[1] = FromFaceIJSame(face, i + size, j, i + size < kMaxSize)
422
neighbors[2] = FromFaceIJSame(face, i, j + size, j + size < kMaxSize)
424
neighbors[3] = FromFaceIJSame(face, i - size, j, i - size >= 0)
428
void S2CellId::AppendVertexNeighbors(int level,
429
vector<S2CellId>* output) const {
430
// "level" must be strictly less than this cell's level so that we can
431
// determine which vertex this cell is closest to.
432
DCHECK_LT(level, this->level());
434
int face = ToFaceIJOrientation(&i, &j, NULL);
436
// Determine the i- and j-offsets to the closest neighboring cell in each
437
// direction. This involves looking at the next bit of "i" and "j" to
438
// determine which quadrant of this->parent(level) this cell lies in.
439
int halfsize = GetSizeIJ(level + 1);
440
int size = halfsize << 1;
442
int ioffset, joffset;
445
isame = (i + size) < kMaxSize;
448
isame = (i - size) >= 0;
452
jsame = (j + size) < kMaxSize;
455
jsame = (j - size) >= 0;
458
output->push_back(parent(level));
459
output->push_back(FromFaceIJSame(face, i + ioffset, j, isame).parent(level));
460
output->push_back(FromFaceIJSame(face, i, j + joffset, jsame).parent(level));
461
// If i- and j- edge neighbors are *both* on a different face, then this
462
// vertex only has three neighbors (it is one of the 8 cube vertices).
463
if (isame || jsame) {
464
output->push_back(FromFaceIJSame(face, i + ioffset, j + joffset,
465
isame && jsame).parent(level));
469
void S2CellId::AppendAllNeighbors(int nbr_level,
470
vector<S2CellId>* output) const {
472
int face = ToFaceIJOrientation(&i, &j, NULL);
474
// Find the coordinates of the lower left-hand leaf cell. We need to
475
// normalize (i,j) to a known position within the cell because nbr_level
476
// may be larger than this cell's level.
477
int size = GetSizeIJ();
481
int nbr_size = GetSizeIJ(nbr_level);
482
DCHECK_LE(nbr_size, size);
484
// We compute the N-S, E-W, and diagonal neighbors in one pass.
485
// The loop test is at the end of the loop to avoid 32-bit overflow.
486
for (int k = -nbr_size; ; k += nbr_size) {
489
same_face = (j + k >= 0);
490
} else if (k >= size) {
491
same_face = (j + k < kMaxSize);
494
// North and South neighbors.
495
output->push_back(FromFaceIJSame(face, i + k, j - nbr_size,
496
j - size >= 0).parent(nbr_level));
497
output->push_back(FromFaceIJSame(face, i + k, j + size,
498
j + size < kMaxSize).parent(nbr_level));
500
// East, West, and Diagonal neighbors.
501
output->push_back(FromFaceIJSame(face, i - nbr_size, j + k,
502
same_face && i - size >= 0)
504
output->push_back(FromFaceIJSame(face, i + size, j + k,
505
same_face && i + size < kMaxSize)
507
if (k >= size) break;
511
inline char intToChar(int i) {
515
inline uint64 charToInt(char c) {
516
return static_cast<uint64>(c - '0');
520
S2CellId S2CellId::FromString(const string& str) {
521
int face = charToInt(str[0]);
522
int level = str.length() - 2;
524
for (size_t i = 2; i < str.length(); ++i) {
525
pos |= charToInt(str[i]) << (2 * (kMaxLevel - (i - 1)) + 1);
527
pos |= (uint64)1 << (2 * (kMaxLevel - level));
528
return FromFacePosLevel(face, pos, level);
531
string S2CellId::ToString() const {
533
return StringPrintf("Invalid: %016llx", id());
536
out.reserve(2 + level());
537
out.push_back(intToChar(face()));
539
for (int current_level = 1; current_level <= level(); ++current_level) {
540
out.push_back(intToChar(child_position(current_level)));
545
string S2CellId::slowToString() const {
547
return StringPrintf("Invalid: %016llx", id());
549
string out = IntToString(face(), "%df");
550
for (int current_level = 1; current_level <= level(); ++current_level) {
551
out += IntToString(child_position(current_level), "%d");
556
ostream& operator<<(ostream& os, S2CellId const& id) {
557
return os << id.ToString();
561
template<> size_t stdext::hash_value<S2CellId>(const S2CellId &id) {
562
return static_cast<size_t>(id.id() >> 32) + static_cast<size_t>(id.id());
added
removed
src/third_party/s2/s2cellid.cc
