1 // Copyright 2005 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
16 // Author: ericv@google.com (Eric Veach)
17
18 #include "s2/s2cell_union.h"
19
20 #include <algorithm>
21 #include <vector>
22
23 #include "s2/base/integral_types.h"
24 #include "s2/base/logging.h"
25 #include "s2/s1angle.h"
26 #include "s2/s2cap.h"
27 #include "s2/s2cell.h"
28 #include "s2/s2cell_id.h"
29 #include "s2/s2latlng_rect.h"
30 #include "s2/s2metrics.h"
31 #include "s2/util/coding/coder.h"
32
33 using std::is_sorted;
34 using std::max;
35 using std::min;
36 using std::vector;
37
38 DEFINE_int32(s2cell_union_decode_max_num_cells, 1000000,
39 "The maximum number of cells allowed by S2CellUnion::Decode");
40
41 static const unsigned char kCurrentLosslessEncodingVersionNumber = 1;
42
ToS2CellIds(const vector<uint64> & ids)43 vector<S2CellId> S2CellUnion::ToS2CellIds(const vector<uint64>& ids) {
44 vector<S2CellId> cell_ids;
45 cell_ids.reserve(ids.size());
46 for (auto id : ids) cell_ids.push_back(S2CellId(id));
47 return cell_ids;
48 }
49
S2CellUnion(const vector<uint64> & cell_ids)50 S2CellUnion::S2CellUnion(const vector<uint64>& cell_ids)
51 : cell_ids_(ToS2CellIds(cell_ids)) {
52 Normalize();
53 }
54
WholeSphere()55 S2CellUnion S2CellUnion::WholeSphere() {
56 return S2CellUnion({S2CellId::FromFace(0), S2CellId::FromFace(1),
57 S2CellId::FromFace(2), S2CellId::FromFace(3),
58 S2CellId::FromFace(4), S2CellId::FromFace(5)});
59 }
60
FromMinMax(S2CellId min_id,S2CellId max_id)61 S2CellUnion S2CellUnion::FromMinMax(S2CellId min_id, S2CellId max_id) {
62 S2CellUnion result;
63 result.InitFromMinMax(min_id, max_id);
64 return result;
65 }
66
FromBeginEnd(S2CellId begin,S2CellId end)67 S2CellUnion S2CellUnion::FromBeginEnd(S2CellId begin, S2CellId end) {
68 S2CellUnion result;
69 result.InitFromBeginEnd(begin, end);
70 return result;
71 }
72
Init(const vector<uint64> & cell_ids)73 void S2CellUnion::Init(const vector<uint64>& cell_ids) {
74 cell_ids_ = ToS2CellIds(cell_ids);
75 Normalize();
76 }
77
InitFromMinMax(S2CellId min_id,S2CellId max_id)78 void S2CellUnion::InitFromMinMax(S2CellId min_id, S2CellId max_id) {
79 S2_DCHECK(max_id.is_valid());
80 InitFromBeginEnd(min_id, max_id.next());
81 }
82
InitFromBeginEnd(S2CellId begin,S2CellId end)83 void S2CellUnion::InitFromBeginEnd(S2CellId begin, S2CellId end) {
84 S2_DCHECK(begin.is_leaf());
85 S2_DCHECK(end.is_leaf());
86 S2_DCHECK_LE(begin, end);
87
88 // We repeatedly add the largest cell we can.
89 cell_ids_.clear();
90 for (S2CellId id = begin.maximum_tile(end);
91 id != end; id = id.next().maximum_tile(end)) {
92 cell_ids_.push_back(id);
93 }
94 // The output is already normalized.
95 S2_DCHECK(IsNormalized());
96 }
97
Pack(int excess)98 void S2CellUnion::Pack(int excess) {
99 if (cell_ids_.capacity() - cell_ids_.size() > excess) {
100 cell_ids_.shrink_to_fit();
101 }
102 }
103
Clone() const104 S2CellUnion* S2CellUnion::Clone() const {
105 return new S2CellUnion(cell_ids_, VERBATIM);
106 }
107
108 // Returns true if the given four cells have a common parent.
109 // REQUIRES: The four cells are distinct.
AreSiblings(S2CellId a,S2CellId b,S2CellId c,S2CellId d)110 inline static bool AreSiblings(S2CellId a, S2CellId b, S2CellId c, S2CellId d) {
111 // A necessary (but not sufficient) condition is that the XOR of the
112 // four cells must be zero. This is also very fast to test.
113 if ((a.id() ^ b.id() ^ c.id()) != d.id()) return false;
114
115 // Now we do a slightly more expensive but exact test. First, compute a
116 // mask that blocks out the two bits that encode the child position of
117 // "id" with respect to its parent, then check that the other three
118 // children all agree with "mask".
119 uint64 mask = d.lsb() << 1;
120 mask = ~(mask + (mask << 1));
121 uint64 id_masked = (d.id() & mask);
122 return ((a.id() & mask) == id_masked &&
123 (b.id() & mask) == id_masked &&
124 (c.id() & mask) == id_masked &&
125 !d.is_face());
126 }
127
IsValid() const128 bool S2CellUnion::IsValid() const {
129 if (num_cells() > 0 && !cell_id(0).is_valid()) return false;
130 for (int i = 1; i < num_cells(); ++i) {
131 if (!cell_id(i).is_valid()) return false;
132 if (cell_id(i - 1).range_max() >= cell_id(i).range_min()) return false;
133 }
134 return true;
135 }
136
IsNormalized() const137 bool S2CellUnion::IsNormalized() const {
138 if (num_cells() > 0 && !cell_id(0).is_valid()) return false;
139 for (int i = 1; i < num_cells(); ++i) {
140 if (!cell_id(i).is_valid()) return false;
141 if (cell_id(i - 1).range_max() >= cell_id(i).range_min()) return false;
142 if (i >= 3 && AreSiblings(cell_id(i - 3), cell_id(i - 2),
143 cell_id(i - 1), cell_id(i))) {
144 return false;
145 }
146 }
147 return true;
148 }
149
Normalize()150 bool S2CellUnion::Normalize() {
151 return Normalize(&cell_ids_);
152 }
153
Normalize(vector<S2CellId> * ids)154 /*static*/ bool S2CellUnion::Normalize(vector<S2CellId>* ids) {
155 // Optimize the representation by discarding cells contained by other cells,
156 // and looking for cases where all subcells of a parent cell are present.
157
158 std::sort(ids->begin(), ids->end());
159 int out = 0;
160 for (S2CellId id : *ids) {
161 // Check whether this cell is contained by the previous cell.
162 if (out > 0 && (*ids)[out-1].contains(id)) continue;
163
164 // Discard any previous cells contained by this cell.
165 while (out > 0 && id.contains((*ids)[out-1])) --out;
166
167 // Check whether the last 3 elements plus "id" can be collapsed into a
168 // single parent cell.
169 while (out >= 3 && AreSiblings((*ids)[out - 3], (*ids)[out - 2],
170 (*ids)[out - 1], id)) {
171 // Replace four children by their parent cell.
172 id = id.parent();
173 out -= 3;
174 }
175 (*ids)[out++] = id;
176 }
177 if (ids->size() == out) return false;
178 ids->resize(out);
179 return true;
180 }
181
Denormalize(int min_level,int level_mod,vector<S2CellId> * out) const182 void S2CellUnion::Denormalize(int min_level, int level_mod,
183 vector<S2CellId>* out) const {
184 Denormalize(cell_ids_, min_level, level_mod, out);
185 }
186
Denormalize(const vector<S2CellId> & in,int min_level,int level_mod,vector<S2CellId> * out)187 void S2CellUnion::Denormalize(const vector<S2CellId>& in,
188 int min_level, int level_mod,
189 vector<S2CellId>* out) {
190 S2_DCHECK_GE(min_level, 0);
191 S2_DCHECK_LE(min_level, S2CellId::kMaxLevel);
192 S2_DCHECK_GE(level_mod, 1);
193 S2_DCHECK_LE(level_mod, 3);
194 S2_DCHECK_NE(out, &in);
195
196 out->clear();
197 out->reserve(in.size());
198 for (S2CellId id : in) {
199 int level = id.level();
200 int new_level = max(min_level, level);
201 if (level_mod > 1) {
202 // Round up so that (new_level - min_level) is a multiple of level_mod.
203 // (Note that S2CellId::kMaxLevel is a multiple of 1, 2, and 3.)
204 new_level += (S2CellId::kMaxLevel - (new_level - min_level)) % level_mod;
205 new_level = min(S2CellId::kMaxLevel, new_level);
206 }
207 if (new_level == level) {
208 out->push_back(id);
209 } else {
210 S2CellId end = id.child_end(new_level);
211 for (id = id.child_begin(new_level); id != end; id = id.next()) {
212 out->push_back(id);
213 }
214 }
215 }
216 }
217
GetCapBound() const218 S2Cap S2CellUnion::GetCapBound() const {
219 // Compute the approximate centroid of the region. This won't produce the
220 // bounding cap of minimal area, but it should be close enough.
221 if (cell_ids_.empty()) return S2Cap::Empty();
222 S2Point centroid(0, 0, 0);
223 for (S2CellId id : *this) {
224 double area = S2Cell::AverageArea(id.level());
225 centroid += area * id.ToPoint();
226 }
227 if (centroid == S2Point(0, 0, 0)) {
228 centroid = S2Point(1, 0, 0);
229 } else {
230 centroid = centroid.Normalize();
231 }
232
233 // Use the centroid as the cap axis, and expand the cap angle so that it
234 // contains the bounding caps of all the individual cells. Note that it is
235 // *not* sufficient to just bound all the cell vertices because the bounding
236 // cap may be concave (i.e. cover more than one hemisphere).
237 S2Cap cap = S2Cap::FromPoint(centroid);
238 for (S2CellId id : *this) {
239 cap.AddCap(S2Cell(id).GetCapBound());
240 }
241 return cap;
242 }
243
GetRectBound() const244 S2LatLngRect S2CellUnion::GetRectBound() const {
245 S2LatLngRect bound = S2LatLngRect::Empty();
246 for (S2CellId id : *this) {
247 bound = bound.Union(S2Cell(id).GetRectBound());
248 }
249 return bound;
250 }
251
Contains(S2CellId id) const252 bool S2CellUnion::Contains(S2CellId id) const {
253 // This is an exact test. Each cell occupies a linear span of the S2
254 // space-filling curve, and the cell id is simply the position at the center
255 // of this span. The cell union ids are sorted in increasing order along
256 // the space-filling curve. So we simply find the pair of cell ids that
257 // surround the given cell id (using binary search). There is containment
258 // if and only if one of these two cell ids contains this cell.
259
260 vector<S2CellId>::const_iterator i =
261 std::lower_bound(cell_ids_.begin(), cell_ids_.end(), id);
262 if (i != cell_ids_.end() && i->range_min() <= id) return true;
263 return i != cell_ids_.begin() && (--i)->range_max() >= id;
264 }
265
Intersects(S2CellId id) const266 bool S2CellUnion::Intersects(S2CellId id) const {
267 // This is an exact test; see the comments for Contains() above.
268
269 vector<S2CellId>::const_iterator i =
270 std::lower_bound(cell_ids_.begin(), cell_ids_.end(), id);
271 if (i != cell_ids_.end() && i->range_min() <= id.range_max()) return true;
272 return i != cell_ids_.begin() && (--i)->range_max() >= id.range_min();
273 }
274
Contains(const S2CellUnion & y) const275 bool S2CellUnion::Contains(const S2CellUnion& y) const {
276 // TODO(ericv): A divide-and-conquer or alternating-skip-search
277 // approach may be sigificantly faster in both the average and worst case.
278
279 for (S2CellId y_id : y) {
280 if (!Contains(y_id)) return false;
281 }
282 return true;
283 }
284
Intersects(const S2CellUnion & y) const285 bool S2CellUnion::Intersects(const S2CellUnion& y) const {
286 // TODO(ericv): A divide-and-conquer or alternating-skip-search
287 // approach may be sigificantly faster in both the average and worst case.
288
289 for (S2CellId y_id : y) {
290 if (Intersects(y_id)) return true;
291 }
292 return false;
293 }
294
Union(const S2CellUnion & y) const295 S2CellUnion S2CellUnion::Union(const S2CellUnion& y) const {
296 vector<S2CellId> cell_ids;
297 cell_ids.reserve(num_cells() + y.num_cells());
298 cell_ids = cell_ids_;
299 cell_ids.insert(cell_ids.end(), y.cell_ids_.begin(), y.cell_ids_.end());
300 return S2CellUnion(std::move(cell_ids));
301 }
302
Intersection(S2CellId id) const303 S2CellUnion S2CellUnion::Intersection(S2CellId id) const {
304 S2CellUnion result;
305 if (Contains(id)) {
306 result.cell_ids_.push_back(id);
307 } else {
308 vector<S2CellId>::const_iterator i = std::lower_bound(
309 cell_ids_.begin(), cell_ids_.end(), id.range_min());
310 S2CellId id_max = id.range_max();
311 while (i != cell_ids_.end() && *i <= id_max)
312 result.cell_ids_.push_back(*i++);
313 }
314 S2_DCHECK(result.IsNormalized() || !IsNormalized());
315 return result;
316 }
317
Intersection(const S2CellUnion & y) const318 S2CellUnion S2CellUnion::Intersection(const S2CellUnion& y) const {
319 S2CellUnion result;
320 GetIntersection(cell_ids_, y.cell_ids_, &result.cell_ids_);
321 // The output is normalized as long as at least one input is normalized.
322 S2_DCHECK(result.IsNormalized() || (!IsNormalized() && !y.IsNormalized()));
323 return result;
324 }
325
GetIntersection(const vector<S2CellId> & x,const vector<S2CellId> & y,vector<S2CellId> * out)326 /*static*/ void S2CellUnion::GetIntersection(const vector<S2CellId>& x,
327 const vector<S2CellId>& y,
328 vector<S2CellId>* out) {
329 S2_DCHECK_NE(out, &x);
330 S2_DCHECK_NE(out, &y);
331 S2_DCHECK(is_sorted(x.begin(), x.end()));
332 S2_DCHECK(is_sorted(y.begin(), y.end()));
333
334 // This is a fairly efficient calculation that uses binary search to skip
335 // over sections of both input vectors. It takes logarithmic time if all the
336 // cells of "x" come before or after all the cells of "y" in S2CellId order.
337
338 out->clear();
339 vector<S2CellId>::const_iterator i = x.begin();
340 vector<S2CellId>::const_iterator j = y.begin();
341 while (i != x.end() && j != y.end()) {
342 S2CellId imin = i->range_min();
343 S2CellId jmin = j->range_min();
344 if (imin > jmin) {
345 // Either j->contains(*i) or the two cells are disjoint.
346 if (*i <= j->range_max()) {
347 out->push_back(*i++);
348 } else {
349 // Advance "j" to the first cell possibly contained by *i.
350 j = std::lower_bound(j + 1, y.end(), imin);
351 // The previous cell *(j-1) may now contain *i.
352 if (*i <= (j - 1)->range_max()) --j;
353 }
354 } else if (jmin > imin) {
355 // Identical to the code above with "i" and "j" reversed.
356 if (*j <= i->range_max()) {
357 out->push_back(*j++);
358 } else {
359 i = std::lower_bound(i + 1, x.end(), jmin);
360 if (*j <= (i - 1)->range_max()) --i;
361 }
362 } else {
363 // "i" and "j" have the same range_min(), so one contains the other.
364 if (*i < *j)
365 out->push_back(*i++);
366 else
367 out->push_back(*j++);
368 }
369 }
370 // The output is generated in sorted order.
371 S2_DCHECK(is_sorted(out->begin(), out->end()));
372 }
373
GetDifferenceInternal(S2CellId cell,const S2CellUnion & y,vector<S2CellId> * cell_ids)374 static void GetDifferenceInternal(S2CellId cell,
375 const S2CellUnion& y,
376 vector<S2CellId>* cell_ids) {
377 // Add the difference between cell and y to cell_ids.
378 // If they intersect but the difference is non-empty, divide and conquer.
379 if (!y.Intersects(cell)) {
380 cell_ids->push_back(cell);
381 } else if (!y.Contains(cell)) {
382 S2CellId child = cell.child_begin();
383 for (int i = 0; ; ++i) {
384 GetDifferenceInternal(child, y, cell_ids);
385 if (i == 3) break; // Avoid unnecessary next() computation.
386 child = child.next();
387 }
388 }
389 }
390
Difference(const S2CellUnion & y) const391 S2CellUnion S2CellUnion::Difference(const S2CellUnion& y) const {
392 // TODO(ericv): this is approximately O(N*log(N)), but could probably
393 // use similar techniques as GetIntersection() to be more efficient.
394
395 S2CellUnion result;
396 for (S2CellId id : *this) {
397 GetDifferenceInternal(id, y, &result.cell_ids_);
398 }
399 // The output is normalized as long as the first argument is normalized.
400 S2_DCHECK(result.IsNormalized() || !IsNormalized());
401 return result;
402 }
403
Expand(int expand_level)404 void S2CellUnion::Expand(int expand_level) {
405 vector<S2CellId> output;
406 uint64 level_lsb = S2CellId::lsb_for_level(expand_level);
407 for (int i = num_cells(); --i >= 0; ) {
408 S2CellId id = cell_id(i);
409 if (id.lsb() < level_lsb) {
410 id = id.parent(expand_level);
411 // Optimization: skip over any cells contained by this one. This is
412 // especially important when very small regions are being expanded.
413 while (i > 0 && id.contains(cell_id(i - 1))) --i;
414 }
415 output.push_back(id);
416 id.AppendAllNeighbors(expand_level, &output);
417 }
418 Init(std::move(output));
419 }
420
Expand(S1Angle min_radius,int max_level_diff)421 void S2CellUnion::Expand(S1Angle min_radius, int max_level_diff) {
422 int min_level = S2CellId::kMaxLevel;
423 for (S2CellId id : *this) {
424 min_level = min(min_level, id.level());
425 }
426 // Find the maximum level such that all cells are at least "min_radius" wide.
427 int radius_level = S2::kMinWidth.GetLevelForMinValue(min_radius.radians());
428 if (radius_level == 0 && min_radius.radians() > S2::kMinWidth.GetValue(0)) {
429 // The requested expansion is greater than the width of a face cell.
430 // The easiest way to handle this is to expand twice.
431 Expand(0);
432 }
433 Expand(min(min_level + max_level_diff, radius_level));
434 }
435
LeafCellsCovered() const436 uint64 S2CellUnion::LeafCellsCovered() const {
437 uint64 num_leaves = 0;
438 for (S2CellId id : *this) {
439 const int inverted_level = S2CellId::kMaxLevel - id.level();
440 num_leaves += (1ULL << (inverted_level << 1));
441 }
442 return num_leaves;
443 }
444
AverageBasedArea() const445 double S2CellUnion::AverageBasedArea() const {
446 return S2Cell::AverageArea(S2CellId::kMaxLevel) * LeafCellsCovered();
447 }
448
ApproxArea() const449 double S2CellUnion::ApproxArea() const {
450 double area = 0;
451 for (S2CellId id : *this) {
452 area += S2Cell(id).ApproxArea();
453 }
454 return area;
455 }
456
ExactArea() const457 double S2CellUnion::ExactArea() const {
458 double area = 0;
459 for (S2CellId id : *this) {
460 area += S2Cell(id).ExactArea();
461 }
462 return area;
463 }
464
operator ==(const S2CellUnion & x,const S2CellUnion & y)465 bool operator==(const S2CellUnion& x, const S2CellUnion& y) {
466 return x.cell_ids() == y.cell_ids();
467 }
468
operator !=(const S2CellUnion & x,const S2CellUnion & y)469 bool operator!=(const S2CellUnion& x, const S2CellUnion& y) {
470 return x.cell_ids() != y.cell_ids();
471 }
472
Contains(const S2Cell & cell) const473 bool S2CellUnion::Contains(const S2Cell& cell) const {
474 return Contains(cell.id());
475 }
476
MayIntersect(const S2Cell & cell) const477 bool S2CellUnion::MayIntersect(const S2Cell& cell) const {
478 return Intersects(cell.id());
479 }
480
Encode(Encoder * const encoder) const481 void S2CellUnion::Encode(Encoder* const encoder) const {
482 // Unsigned char for version number, and N+1 uint64's for N cell_ids
483 // (1 for vector length, N for the ids).
484 encoder->Ensure(sizeof(unsigned char) +
485 sizeof(uint64) * (1 + cell_ids_.size()));
486
487 encoder->put8(kCurrentLosslessEncodingVersionNumber);
488 encoder->put64(uint64{cell_ids_.size()});
489 for (const S2CellId& cell_id : cell_ids_) {
490 cell_id.Encode(encoder);
491 }
492 }
493
Decode(Decoder * const decoder)494 bool S2CellUnion::Decode(Decoder* const decoder) {
495 // Should contain at least version and vector length.
496 if (decoder->avail() < sizeof(unsigned char) + sizeof(uint64)) return false;
497 unsigned char version = decoder->get8();
498 if (version > kCurrentLosslessEncodingVersionNumber) return false;
499
500 uint64 num_cells = decoder->get64();
501 if (num_cells > FLAGS_s2cell_union_decode_max_num_cells) {
502 return false;
503 }
504
505 vector<S2CellId> temp_cell_ids(num_cells);
506 for (int i = 0; i < num_cells; ++i) {
507 if (!temp_cell_ids[i].Decode(decoder)) return false;
508 }
509 cell_ids_.swap(temp_cell_ids);
510 return true;
511 }
512
Contains(const S2Point & p) const513 bool S2CellUnion::Contains(const S2Point& p) const {
514 return Contains(S2CellId(p));
515 }
516