1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#include "s2cellunion.h"
4
5#include <algorithm>
6using std::min;
7using std::max;
8using std::swap;
9using std::reverse;
10
11#include <vector>
12using std::vector;
13
14
15#include "base/integral_types.h"
16#include "base/logging.h"
17#include "s2.h"
18#include "s2cap.h"
19#include "s2cell.h"
20#include "s2cellid.h"
21#include "s2latlngrect.h"
22
23// Returns true if the vector of cell_ids is sorted. Used only in
24// DCHECKs.
25extern bool IsSorted(vector<S2CellId> const& cell_ids) {
26 for (int i = 0; i + 1 < cell_ids.size(); ++i) {
27 if (cell_ids[i + 1] < cell_ids[i]) {
28 return false;
29 }
30 }
31 return true;
32}
33
34void S2CellUnion::Init(vector<S2CellId> const& cell_ids) {
35 InitRaw(cell_ids);
36 Normalize();
37}
38
39void S2CellUnion::Init(vector<uint64> const& cell_ids) {
40 InitRaw(cell_ids);
41 Normalize();
42}
43
44void S2CellUnion::InitSwap(vector<S2CellId>* cell_ids) {
45 InitRawSwap(cell_ids);
46 Normalize();
47}
48
49void S2CellUnion::InitRaw(vector<S2CellId> const& cell_ids) {
50 cell_ids_ = cell_ids;
51}
52
53void S2CellUnion::InitRaw(vector<uint64> const& cell_ids) {
54 cell_ids_.resize(cell_ids.size());
55 for (int i = 0; i < num_cells(); ++i) {
56 cell_ids_[i] = S2CellId(cell_ids[i]);
57 }
58}
59
60void S2CellUnion::InitRawSwap(vector<S2CellId>* cell_ids) {
61 cell_ids_.swap(*cell_ids);
62 cell_ids->clear();
63}
64
65void S2CellUnion::Detach(vector<S2CellId>* cell_ids) {
66 cell_ids_.swap(*cell_ids);
67 cell_ids_.clear();
68}
69
70void S2CellUnion::Pack(int excess) {
71 if (cell_ids_.capacity() - cell_ids_.size() > excess) {
72 vector<S2CellId> packed = cell_ids_;
73 cell_ids_.swap(packed);
74 }
75}
76
77S2CellUnion* S2CellUnion::Clone() const {
78 S2CellUnion* copy = new S2CellUnion;
79 copy->InitRaw(cell_ids_);
80 return copy;
81}
82
83bool S2CellUnion::Normalize() {
84 // Optimize the representation by looking for cases where all subcells
85 // of a parent cell are present.
86
87 vector<S2CellId> output;
88 output.reserve(cell_ids_.size());
89 sort(cell_ids_.begin(), cell_ids_.end());
90
91 for (int i = 0; i < num_cells(); ++i) {
92 S2CellId id = cell_id(i);
93
94 // Check whether this cell is contained by the previous cell.
95 if (!output.empty() && output.back().contains(id)) continue;
96
97 // Discard any previous cells contained by this cell.
98 while (!output.empty() && id.contains(output.back())) {
99 output.pop_back();
100 }
101
102 // Check whether the last 3 elements of "output" plus "id" can be
103 // collapsed into a single parent cell.
104 while (output.size() >= 3) {
105 // A necessary (but not sufficient) condition is that the XOR of the
106 // four cells must be zero. This is also very fast to test.
107 if ((output.end()[-3].id() ^ output.end()[-2].id() ^ output.back().id())
108 != id.id())
109 break;
110
111 // Now we do a slightly more expensive but exact test. First, compute a
112 // mask that blocks out the two bits that encode the child position of
113 // "id" with respect to its parent, then check that the other three
114 // children all agree with "mask.
115 uint64 mask = id.lsb() << 1;
116 mask = ~(mask + (mask << 1));
117 uint64 id_masked = (id.id() & mask);
118 if ((output.end()[-3].id() & mask) != id_masked ||
119 (output.end()[-2].id() & mask) != id_masked ||
120 (output.end()[-1].id() & mask) != id_masked ||
121 id.is_face())
122 break;
123
124 // Replace four children by their parent cell.
125 output.erase(output.end() - 3, output.end());
126 id = id.parent();
127 }
128 output.push_back(id);
129 }
130 if (output.size() < num_cells()) {
131 InitRawSwap(&output);
132 return true;
133 }
134 return false;
135}
136
137void S2CellUnion::Denormalize(int min_level, int level_mod,
138 vector<S2CellId>* output) const {
139 DCHECK_GE(min_level, 0);
140 DCHECK_LE(min_level, S2CellId::kMaxLevel);
141 DCHECK_GE(level_mod, 1);
142 DCHECK_LE(level_mod, 3);
143
144 output->clear();
145 output->reserve(num_cells());
146 for (int i = 0; i < num_cells(); ++i) {
147 S2CellId id = cell_id(i);
148 int level = id.level();
149 int new_level = max(min_level, level);
150 if (level_mod > 1) {
151 // Round up so that (new_level - min_level) is a multiple of level_mod.
152 // (Note that S2CellId::kMaxLevel is a multiple of 1, 2, and 3.)
153 new_level += (S2CellId::kMaxLevel - (new_level - min_level)) % level_mod;
154 new_level = min(S2CellId::kMaxLevel, new_level);
155 }
156 if (new_level == level) {
157 output->push_back(id);
158 } else {
159 S2CellId end = id.child_end(new_level);
160 for (id = id.child_begin(new_level); id != end; id = id.next()) {
161 output->push_back(id);
162 }
163 }
164 }
165}
166
167S2Cap S2CellUnion::GetCapBound() const {
168 // Compute the approximate centroid of the region. This won't produce the
169 // bounding cap of minimal area, but it should be close enough.
170 if (cell_ids_.empty()) return S2Cap::Empty();
171 S2Point centroid(0, 0, 0);
172 for (int i = 0; i < num_cells(); ++i) {
173 double area = S2Cell::AverageArea(cell_id(i).level());
174 centroid += area * cell_id(i).ToPoint();
175 }
176 if (centroid == S2Point(0, 0, 0)) {
177 centroid = S2Point(1, 0, 0);
178 } else {
179 centroid = centroid.Normalize();
180 }
181
182 // Use the centroid as the cap axis, and expand the cap angle so that it
183 // contains the bounding caps of all the individual cells. Note that it is
184 // *not* sufficient to just bound all the cell vertices because the bounding
185 // cap may be concave (i.e. cover more than one hemisphere).
186 S2Cap cap = S2Cap::FromAxisHeight(centroid, 0);
187 for (int i = 0; i < num_cells(); ++i) {
188 cap.AddCap(S2Cell(cell_id(i)).GetCapBound());
189 }
190 return cap;
191}
192
193S2LatLngRect S2CellUnion::GetRectBound() const {
194 S2LatLngRect bound = S2LatLngRect::Empty();
195 for (int i = 0; i < num_cells(); ++i) {
196 bound = bound.Union(S2Cell(cell_id(i)).GetRectBound());
197 }
198 return bound;
199}
200
201bool S2CellUnion::Contains(S2CellId const& id) const {
202 // This function requires that Normalize has been called first.
203 //
204 // This is an exact test. Each cell occupies a linear span of the S2
205 // space-filling curve, and the cell id is simply the position at the center
206 // of this span. The cell union ids are sorted in increasing order along
207 // the space-filling curve. So we simply find the pair of cell ids that
208 // surround the given cell id (using binary search). There is containment
209 // if and only if one of these two cell ids contains this cell.
210
211 vector<S2CellId>::const_iterator i =
212 lower_bound(cell_ids_.begin(), cell_ids_.end(), id);
213 if (i != cell_ids_.end() && i->range_min() <= id) return true;
214 return i != cell_ids_.begin() && (--i)->range_max() >= id;
215}
216
217bool S2CellUnion::Intersects(S2CellId const& id) const {
218 // This function requires that Normalize has been called first.
219 // This is an exact test; see the comments for Contains() above.
220
221 vector<S2CellId>::const_iterator i =
222 lower_bound(cell_ids_.begin(), cell_ids_.end(), id);
223 if (i != cell_ids_.end() && i->range_min() <= id.range_max()) return true;
224 return i != cell_ids_.begin() && (--i)->range_max() >= id.range_min();
225}
226
227bool S2CellUnion::Contains(S2CellUnion const* y) const {
228 // TODO: A divide-and-conquer or alternating-skip-search approach may be
229 // sigificantly faster in both the average and worst case.
230
231 for (int i = 0; i < y->num_cells(); ++i) {
232 if (!Contains(y->cell_id(i))) return false;
233 }
234 return true;
235}
236
237bool S2CellUnion::Intersects(S2CellUnion const* y) const {
238 // TODO: A divide-and-conquer or alternating-skip-search approach may be
239 // sigificantly faster in both the average and worst case.
240
241 for (int i = 0; i < y->num_cells(); ++i) {
242 if (Intersects(y->cell_id(i))) return true;
243 }
244 return false;
245}
246
247void S2CellUnion::GetUnion(S2CellUnion const* x, S2CellUnion const* y) {
248 DCHECK_NE(this, x);
249 DCHECK_NE(this, y);
250 cell_ids_.reserve(x->num_cells() + y->num_cells());
251 cell_ids_ = x->cell_ids_;
252 cell_ids_.insert(cell_ids_.end(), y->cell_ids_.begin(), y->cell_ids_.end());
253 Normalize();
254}
255
256void S2CellUnion::GetIntersection(S2CellUnion const* x, S2CellId const& id) {
257 DCHECK_NE(this, x);
258 cell_ids_.clear();
259 if (x->Contains(id)) {
260 cell_ids_.push_back(id);
261 } else {
262 vector<S2CellId>::const_iterator i =
263 lower_bound(x->cell_ids_.begin(), x->cell_ids_.end(), id.range_min());
264 S2CellId idmax = id.range_max();
265 while (i != x->cell_ids_.end() && *i <= idmax)
266 cell_ids_.push_back(*i++);
267 }
268}
269
270void S2CellUnion::GetIntersection(S2CellUnion const* x, S2CellUnion const* y) {
271 DCHECK_NE(this, x);
272 DCHECK_NE(this, y);
273
274 // This is a fairly efficient calculation that uses binary search to skip
275 // over sections of both input vectors. It takes constant time if all the
276 // cells of "x" come before or after all the cells of "y" in S2CellId order.
277
278 cell_ids_.clear();
279 vector<S2CellId>::const_iterator i = x->cell_ids_.begin();
280 vector<S2CellId>::const_iterator j = y->cell_ids_.begin();
281 while (i != x->cell_ids_.end() && j != y->cell_ids_.end()) {
282 S2CellId imin = i->range_min();
283 S2CellId jmin = j->range_min();
284 if (imin > jmin) {
285 // Either j->contains(*i) or the two cells are disjoint.
286 if (*i <= j->range_max()) {
287 cell_ids_.push_back(*i++);
288 } else {
289 // Advance "j" to the first cell possibly contained by *i.
290 j = lower_bound(j + 1, y->cell_ids_.end(), imin);
291 // The previous cell *(j-1) may now contain *i.
292 if (*i <= (j - 1)->range_max()) --j;
293 }
294 } else if (jmin > imin) {
295 // Identical to the code above with "i" and "j" reversed.
296 if (*j <= i->range_max()) {
297 cell_ids_.push_back(*j++);
298 } else {
299 i = lower_bound(i + 1, x->cell_ids_.end(), jmin);
300 if (*j <= (i - 1)->range_max()) --i;
301 }
302 } else {
303 // "i" and "j" have the same range_min(), so one contains the other.
304 if (*i < *j)
305 cell_ids_.push_back(*i++);
306 else
307 cell_ids_.push_back(*j++);
308 }
309 }
310 // The output is generated in sorted order, and there should not be any
311 // cells that can be merged (provided that both inputs were normalized).
312 DCHECK(IsSorted(cell_ids_));
313 DCHECK(!Normalize());
314}
315
316static void GetDifferenceInternal(S2CellId cell,
317 S2CellUnion const* y,
318 vector<S2CellId>* cell_ids) {
319 // Add the difference between cell and y to cell_ids.
320 // If they intersect but the difference is non-empty, divides and conquers.
321
322 if (!y->Intersects(cell)) {
323 cell_ids->push_back(cell);
324 } else if (!y->Contains(cell)) {
325 S2CellId child = cell.child_begin();
326 for (int i = 0; ; ++i) {
327 GetDifferenceInternal(child, y, cell_ids);
328 if (i == 3) break; // Avoid unnecessary next() computation.
329 child = child.next();
330 }
331 }
332}
333
334void S2CellUnion::GetDifference(S2CellUnion const* x, S2CellUnion const* y) {
335 DCHECK_NE(this, x);
336 DCHECK_NE(this, y);
337 // TODO: this is approximately O(N*log(N)), but could probably use similar
338 // techniques as GetIntersection() to be more efficient.
339
340 cell_ids_.clear();
341 for (int i = 0; i < x->num_cells(); ++i) {
342 GetDifferenceInternal(x->cell_id(i), y, &cell_ids_);
343 }
344 // The output is generated in sorted order, and there should not be any
345 // cells that can be merged (provided that both inputs were normalized).
346 DCHECK(IsSorted(cell_ids_));
347 DCHECK(!Normalize());
348}
349
350void S2CellUnion::Expand(int level) {
351 vector<S2CellId> output;
352 uint64 level_lsb = S2CellId::lsb_for_level(level);
353 for (int i = num_cells() - 1; i >= 0; --i) {
354 S2CellId id = cell_id(i);
355 if (id.lsb() < level_lsb) {
356 id = id.parent(level);
357 // Optimization: skip over any cells contained by this one. This is
358 // especially important when very small regions are being expanded.
359 while (i > 0 && id.contains(cell_id(i-1))) --i;
360 }
361 output.push_back(id);
362 id.AppendAllNeighbors(level, &output);
363 }
364 InitSwap(&output);
365}
366
367void S2CellUnion::Expand(S1Angle const& min_radius, int max_level_diff) {
368 int min_level = S2CellId::kMaxLevel;
369 for (int i = 0; i < num_cells(); ++i) {
370 min_level = min(min_level, cell_id(i).level());
371 }
372 // Find the maximum level such that all cells are at least "min_radius" wide.
373 int radius_level = S2::kMinWidth.GetMaxLevel(min_radius.radians());
374 if (radius_level == 0 && min_radius.radians() > S2::kMinWidth.GetValue(0)) {
375 // The requested expansion is greater than the width of a face cell.
376 // The easiest way to handle this is to expand twice.
377 Expand(0);
378 }
379 Expand(min(min_level + max_level_diff, radius_level));
380}
381
382void S2CellUnion::InitFromRange(S2CellId const& min_id,
383 S2CellId const& max_id) {
384 DCHECK(min_id.is_leaf());
385 DCHECK(max_id.is_leaf());
386 DCHECK_LE(min_id, max_id);
387
388 // We repeatedly add the largest cell we can.
389 cell_ids_.clear();
390 for (S2CellId next_min_id = min_id; next_min_id <= max_id; ) {
391 DCHECK(next_min_id.is_leaf());
392
393 // Find the largest cell that starts at "next_min_id" and doesn't extend
394 // beyond "max_id".
395 S2CellId next_id = next_min_id;
396 while (!next_id.is_face() &&
397 next_id.parent().range_min() == next_min_id &&
398 next_id.parent().range_max() <= max_id) {
399 next_id = next_id.parent();
400 }
401 cell_ids_.push_back(next_id);
402 next_min_id = next_id.range_max().next();
403 }
404
405 // The output is already normalized.
406 DCHECK(IsSorted(cell_ids_));
407 DCHECK(!Normalize());
408}
409
410uint64 S2CellUnion::LeafCellsCovered() const {
411 uint64 num_leaves = 0;
412 for (int i = 0; i < num_cells(); ++i) {
413 const int inverted_level =
414 S2CellId::kMaxLevel - cell_id(i).level();
415 num_leaves += (1ULL << (inverted_level << 1));
416 }
417 return num_leaves;
418}
419
420double S2CellUnion::AverageBasedArea() const {
421 return S2Cell::AverageArea(S2CellId::kMaxLevel) * LeafCellsCovered();
422}
423
424double S2CellUnion::ApproxArea() const {
425 double area = 0;
426 for (int i = 0; i < num_cells(); ++i) {
427 area += S2Cell(cell_id(i)).ApproxArea();
428 }
429 return area;
430}
431
432double S2CellUnion::ExactArea() const {
433 double area = 0;
434 for (int i = 0; i < num_cells(); ++i) {
435 area += S2Cell(cell_id(i)).ExactArea();
436 }
437 return area;
438}
439
440bool operator==(S2CellUnion const& x, S2CellUnion const& y) {
441 return x.cell_ids() == y.cell_ids();
442}
443
444bool S2CellUnion::Contains(S2Cell const& cell) const {
445 return Contains(cell.id());
446}
447
448bool S2CellUnion::MayIntersect(S2Cell const& cell) const {
449 return Intersects(cell.id());
450}
451
452bool S2CellUnion::Contains(S2Point const& p) const {
453 return Contains(S2CellId::FromPoint(p));
454}
455