1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#include <algorithm>
4using std::min;
5using std::max;
6using std::swap;
7using std::reverse;
8
9#include <hash_map>
10using __gnu_cxx::hash_map;
11
12#include <set>
13using std::set;
14using std::multiset;
15
16#include <vector>
17using std::vector;
18
19
20// #include "base/commandlineflags.h"
21#include "s2polygon.h"
22
23#include "base/port.h" // for HASH_NAMESPACE_DECLARATION_START
24#include "util/coding/coder.h"
25#include "s2edgeindex.h"
26#include "s2cap.h"
27#include "s2cell.h"
28#include "s2cellunion.h"
29#include "s2latlngrect.h"
30#include "s2polygonbuilder.h"
31#include "s2polyline.h"
32
33// DECLARE_bool(s2debug); // defined in s2.cc
34
35static const unsigned char kCurrentEncodingVersionNumber = 1;
36
37typedef pair<S2Point, S2Point> S2Edge;
38
39S2Polygon::S2Polygon()
40 : loops_(),
41 bound_(S2LatLngRect::Empty()),
42 owns_loops_(true),
43 has_holes_(false),
44 num_vertices_(0) {
45}
46
47S2Polygon::S2Polygon(vector<S2Loop*>* loops)
48 : bound_(S2LatLngRect::Empty()),
49 owns_loops_(true) {
50 Init(loops);
51}
52
53S2Polygon::S2Polygon(S2Cell const& cell)
54 : bound_(S2LatLngRect::Empty()),
55 owns_loops_(true),
56 has_holes_(false),
57 num_vertices_(4) {
58 S2Loop* loop = new S2Loop(cell);
59 bound_ = loop->GetRectBound();
60 loops_.push_back(loop);
61}
62
63S2Polygon::S2Polygon(S2Loop* loop)
64 : bound_(loop->GetRectBound()),
65 owns_loops_(false),
66 has_holes_(false),
67 num_vertices_(loop->num_vertices()) {
68 loops_.push_back(loop);
69}
70
71void S2Polygon::Copy(S2Polygon const* src) {
72 DCHECK_EQ(0, num_loops());
73 for (int i = 0; i < src->num_loops(); ++i) {
74 loops_.push_back(src->loop(i)->Clone());
75 }
76 bound_ = src->bound_;
77 owns_loops_ = true;
78 has_holes_ = src->has_holes_;
79 num_vertices_ = src->num_vertices();
80}
81
82S2Polygon* S2Polygon::Clone() const {
83 S2Polygon* result = new S2Polygon;
84 result->Copy(this);
85 return result;
86}
87
88void S2Polygon::Release(vector<S2Loop*>* loops) {
89 if (loops != NULL) {
90 loops->insert(loops->end(), loops_.begin(), loops_.end());
91 }
92 loops_.clear();
93 bound_ = S2LatLngRect::Empty();
94 has_holes_ = false;
95}
96
97static void DeleteLoopsInVector(vector<S2Loop*>* loops) {
98 for (int i = 0; i < loops->size(); ++i) {
99 delete loops->at(i);
100 }
101 loops->clear();
102}
103
104S2Polygon::~S2Polygon() {
105 if (owns_loops_) DeleteLoopsInVector(&loops_);
106}
107
108typedef pair<S2Point, S2Point> S2PointPair;
109
110#include<hash_set>
111namespace __gnu_cxx {
112
113template<> struct hash<S2PointPair> {
114 size_t operator()(S2PointPair const& p) const {
115 hash<S2Point> h;
116 return h(p.first) + (h(p.second) << 1);
117 }
118};
119
120} // namespace __gnu_cxx
121
122
123bool S2Polygon::IsValid(const vector<S2Loop*>& loops) {
124 // If a loop contains an edge AB, then no other loop may contain AB or BA.
125 if (loops.size() > 1) {
126 hash_map<S2PointPair, pair<int, int> > edges;
127 for (int i = 0; i < loops.size(); ++i) {
128 S2Loop* lp = loops[i];
129 for (int j = 0; j < lp->num_vertices(); ++j) {
130 S2PointPair key = make_pair(lp->vertex(j), lp->vertex(j + 1));
131 if (edges.insert(make_pair(key, make_pair(i, j))).second) {
132 key = make_pair(lp->vertex(j + 1), lp->vertex(j));
133 if (edges.insert(make_pair(key, make_pair(i, j))).second)
134 continue;
135 }
136 pair<int, int> other = edges[key];
137 VLOG(2) << "Duplicate edge: loop " << i << ", edge " << j
138 << " and loop " << other.first << ", edge " << other.second;
139 return false;
140 }
141 }
142 }
143
144 // Verify that no loop covers more than half of the sphere, and that no
145 // two loops cross.
146 for (int i = 0; i < loops.size(); ++i) {
147 if (!loops[i]->IsNormalized()) {
148 VLOG(2) << "Loop " << i << " encloses more than half the sphere";
149 return false;
150 }
151 for (int j = i + 1; j < loops.size(); ++j) {
152 // This test not only checks for edge crossings, it also detects
153 // cases where the two boundaries cross at a shared vertex.
154 if (loops[i]->ContainsOrCrosses(loops[j]) < 0) {
155 VLOG(2) << "Loop " << i << " crosses loop " << j;
156 return false;
157 }
158 }
159 }
160
161 return true;
162}
163
164bool S2Polygon::IsValid() const {
165 for (int i = 0; i < num_loops(); ++i) {
166 if (!loop(i)->IsValid()) {
167 return false;
168 }
169 }
170 return IsValid(loops_);
171}
172
173bool S2Polygon::IsValid(bool check_loops, int max_adjacent) const {
174 return IsValid();
175}
176
177void S2Polygon::InsertLoop(S2Loop* new_loop, S2Loop* parent,
178 LoopMap* loop_map) {
179 vector<S2Loop*>* children = &(*loop_map)[parent];
180 for (int i = 0; i < children->size(); ++i) {
181 S2Loop* child = (*children)[i];
182 if (child->ContainsNested(new_loop)) {
183 InsertLoop(new_loop, child, loop_map);
184 return;
185 }
186 }
187 // No loop may contain the complement of another loop. (Handling this case
188 // is significantly more complicated.)
189 DCHECK(parent == NULL || !new_loop->ContainsNested(parent));
190
191 // Some of the children of the parent loop may now be children of
192 // the new loop.
193 vector<S2Loop*>* new_children = &(*loop_map)[new_loop];
194 for (int i = 0; i < children->size();) {
195 S2Loop* child = (*children)[i];
196 if (new_loop->ContainsNested(child)) {
197 new_children->push_back(child);
198 children->erase(children->begin() + i);
199 } else {
200 ++i;
201 }
202 }
203 children->push_back(new_loop);
204}
205
206void S2Polygon::InitLoop(S2Loop* loop, int depth, LoopMap* loop_map) {
207 if (loop) {
208 loop->set_depth(depth);
209 loops_.push_back(loop);
210 }
211 vector<S2Loop*> const& children = (*loop_map)[loop];
212 for (int i = 0; i < children.size(); ++i) {
213 InitLoop(children[i], depth + 1, loop_map);
214 }
215}
216
217bool S2Polygon::ContainsChild(S2Loop* a, S2Loop* b, LoopMap const& loop_map) {
218 // This function is just used to verify that the loop map was
219 // constructed correctly.
220
221 if (a == b) return true;
222 vector<S2Loop*> const& children = loop_map.find(a)->second;
223 for (int i = 0; i < children.size(); ++i) {
224 if (ContainsChild(children[i], b, loop_map)) return true;
225 }
226 return false;
227}
228
229void S2Polygon::Init(vector<S2Loop*>* loops) {
230 // if (FLAGS_s2debug) CHECK(IsValid(*loops));
231 if (S2::debug) CHECK(IsValid(*loops));
232 DCHECK(loops_.empty());
233 loops_.swap(*loops);
234
235 num_vertices_ = 0;
236 for (int i = 0; i < num_loops(); ++i) {
237 num_vertices_ += loop(i)->num_vertices();
238 }
239
240 LoopMap loop_map;
241 for (int i = 0; i < num_loops(); ++i) {
242 InsertLoop(loop(i), NULL, &loop_map);
243 }
244 // Reorder the loops in depth-first traversal order.
245 loops_.clear();
246 InitLoop(NULL, -1, &loop_map);
247
248 // if (FLAGS_s2debug) {
249 if (S2::debug) {
250 // Check that the LoopMap is correct (this is fairly cheap).
251 for (int i = 0; i < num_loops(); ++i) {
252 for (int j = 0; j < num_loops(); ++j) {
253 if (i == j) continue;
254 CHECK_EQ(ContainsChild(loop(i), loop(j), loop_map),
255 loop(i)->ContainsNested(loop(j)));
256 }
257 }
258 }
259
260 // Compute the bounding rectangle of the entire polygon.
261 has_holes_ = false;
262 bound_ = S2LatLngRect::Empty();
263 for (int i = 0; i < num_loops(); ++i) {
264 if (loop(i)->sign() < 0) {
265 has_holes_ = true;
266 } else {
267 bound_ = bound_.Union(loop(i)->GetRectBound());
268 }
269 }
270}
271
272int S2Polygon::GetParent(int k) const {
273 int depth = loop(k)->depth();
274 if (depth == 0) return -1; // Optimization.
275 while (--k >= 0 && loop(k)->depth() >= depth) continue;
276 return k;
277}
278
279int S2Polygon::GetLastDescendant(int k) const {
280 if (k < 0) return num_loops() - 1;
281 int depth = loop(k)->depth();
282 while (++k < num_loops() && loop(k)->depth() > depth) continue;
283 return k - 1;
284}
285
286double S2Polygon::GetArea() const {
287 double area = 0;
288 for (int i = 0; i < num_loops(); ++i) {
289 area += loop(i)->sign() * loop(i)->GetArea();
290 }
291 return area;
292}
293
294S2Point S2Polygon::GetCentroid() const {
295 S2Point centroid;
296 for (int i = 0; i < num_loops(); ++i) {
297 centroid += loop(i)->sign() * loop(i)->GetCentroid();
298 }
299 return centroid;
300}
301
302int S2Polygon::ContainsOrCrosses(S2Loop const* b) const {
303 bool inside = false;
304 for (int i = 0; i < num_loops(); ++i) {
305 int result = loop(i)->ContainsOrCrosses(b);
306 if (result < 0) return -1; // The loop boundaries intersect.
307 if (result > 0) inside ^= true;
308 }
309 return static_cast<int>(inside); // True if loop B is contained by the
310 // polygon.
311}
312
313bool S2Polygon::AnyLoopContains(S2Loop const* b) const {
314 // Return true if any loop contains the given loop.
315 for (int i = 0; i < num_loops(); ++i) {
316 if (loop(i)->Contains(b)) return true;
317 }
318 return false;
319}
320
321bool S2Polygon::ContainsAllShells(S2Polygon const* b) const {
322 // Return true if this polygon (A) contains all the shells of B.
323 for (int j = 0; j < b->num_loops(); ++j) {
324 if (b->loop(j)->sign() < 0) continue;
325 if (ContainsOrCrosses(b->loop(j)) <= 0) {
326 // Shell of B is not contained by A, or the boundaries intersect.
327 return false;
328 }
329 }
330 return true;
331}
332
333bool S2Polygon::ExcludesAllHoles(S2Polygon const* b) const {
334 // Return true if this polygon (A) excludes (i.e. does not intersect)
335 // all holes of B.
336 for (int j = 0; j < b->num_loops(); ++j) {
337 if (b->loop(j)->sign() > 0) continue;
338 if (ContainsOrCrosses(b->loop(j)) != 0) {
339 // Hole of B is contained by A, or the boundaries intersect.
340 return false;
341 }
342 }
343 return true;
344}
345
346bool S2Polygon::IntersectsAnyShell(S2Polygon const* b) const {
347 // Return true if this polygon (A) intersects any shell of B.
348 for (int j = 0; j < b->num_loops(); ++j) {
349 if (b->loop(j)->sign() < 0) continue;
350 if (IntersectsShell(b->loop(j)) != 0)
351 return true;
352 }
353 return false;
354}
355
356bool S2Polygon::IntersectsShell(S2Loop const* b) const {
357 bool inside = false;
358 for (int i = 0; i < num_loops(); ++i) {
359 if (loop(i)->Contains(b)) {
360 inside ^= true;
361 } else if (!b->Contains(loop(i)) && loop(i)->Intersects(b)) {
362 // We definitely have an intersection if the loops intersect AND one
363 // is not properly contained in the other. If A (this) is properly
364 // contained in a loop of B, we don't know yet if it may be actually
365 // inside a hole within B.
366 return true;
367 }
368 }
369 return inside;
370}
371
372bool S2Polygon::Contains(S2Polygon const* b) const {
373 // If both polygons have one loop, use the more efficient S2Loop method.
374 // Note that S2Loop::Contains does its own bounding rectangle check.
375 if (num_loops() == 1 && b->num_loops() == 1) {
376 return loop(0)->Contains(b->loop(0));
377 }
378
379 // Otherwise if neither polygon has holes, we can still use the more
380 // efficient S2Loop::Contains method (rather than ContainsOrCrosses),
381 // but it's worthwhile to do our own bounds check first.
382 if (!bound_.Contains(b->bound_)) {
383 // If the union of the bounding boxes spans the full longitude range,
384 // it is still possible that polygon A contains B. (This is only
385 // possible if at least one polygon has multiple shells.)
386 if (!bound_.lng().Union(b->bound_.lng()).is_full()) return false;
387 }
388 if (!has_holes_ && !b->has_holes_) {
389 for (int j = 0; j < b->num_loops(); ++j) {
390 if (!AnyLoopContains(b->loop(j))) return false;
391 }
392 return true;
393 }
394
395 // This could be implemented more efficiently for polygons with lots of
396 // holes by keeping a copy of the LoopMap computed during initialization.
397 // However, in practice most polygons are one loop, and multiloop polygons
398 // tend to consist of many shells rather than holes. In any case, the real
399 // way to get more efficiency is to implement a sub-quadratic algorithm
400 // such as building a trapezoidal map.
401
402 // Every shell of B must be contained by an odd number of loops of A,
403 // and every hole of A must be contained by an even number of loops of B.
404 return ContainsAllShells(b) && b->ExcludesAllHoles(this);
405}
406
407bool S2Polygon::Intersects(S2Polygon const* b) const {
408 // A.Intersects(B) if and only if !Complement(A).Contains(B). However,
409 // implementing a Complement() operation is trickier than it sounds,
410 // and in any case it's more efficient to test for intersection directly.
411
412 // If both polygons have one loop, use the more efficient S2Loop method.
413 // Note that S2Loop::Intersects does its own bounding rectangle check.
414 if (num_loops() == 1 && b->num_loops() == 1) {
415 return loop(0)->Intersects(b->loop(0));
416 }
417
418 // Otherwise if neither polygon has holes, we can still use the more
419 // efficient S2Loop::Intersects method. The polygons intersect if and
420 // only if some pair of loop regions intersect.
421 if (!bound_.Intersects(b->bound_)) return false;
422 if (!has_holes_ && !b->has_holes_) {
423 for (int i = 0; i < num_loops(); ++i) {
424 for (int j = 0; j < b->num_loops(); ++j) {
425 if (loop(i)->Intersects(b->loop(j))) return true;
426 }
427 }
428 return false;
429 }
430
431 // Otherwise if any shell of B is contained by an odd number of loops of A,
432 // or any shell of A is contained by an odd number of loops of B, or there is
433 // an intersection without containment, then there is an intersection.
434 return IntersectsAnyShell(b) || b->IntersectsAnyShell(this);
435}
436
437S2Cap S2Polygon::GetCapBound() const {
438 return bound_.GetCapBound();
439}
440
441bool S2Polygon::Contains(S2Cell const& cell) const {
442 if (num_loops() == 1) {
443 return loop(0)->Contains(cell);
444 }
445
446 // We can't check bound_.Contains(cell.GetRectBound()) because S2Cell's
447 // GetRectBound() calculation is not precise.
448 if (!bound_.Contains(cell.GetCenter())) return false;
449
450 S2Loop cell_loop(cell);
451 S2Polygon cell_poly(&cell_loop);
452 bool contains = Contains(&cell_poly);
453 if (contains) DCHECK(Contains(cell.GetCenter()));
454 return contains;
455}
456
457bool S2Polygon::ApproxContains(S2Polygon const* b,
458 S1Angle vertex_merge_radius) const {
459 S2Polygon difference;
460 difference.InitToDifferenceSloppy(b, this, vertex_merge_radius);
461 return difference.num_loops() == 0;
462}
463
464bool S2Polygon::MayIntersect(S2Cell const& cell) const {
465 if (num_loops() == 1) {
466 return loop(0)->MayIntersect(cell);
467 }
468 if (!bound_.Intersects(cell.GetRectBound())) return false;
469
470 S2Loop cell_loop(cell);
471 S2Polygon cell_poly(&cell_loop);
472 bool intersects = Intersects(&cell_poly);
473 if (!intersects) DCHECK(!Contains(cell.GetCenter()));
474 return intersects;
475}
476
477bool S2Polygon::VirtualContainsPoint(S2Point const& p) const {
478 return Contains(p); // The same as Contains() below, just virtual.
479}
480
481bool S2Polygon::Contains(S2Point const& p) const {
482 if (num_loops() == 1) {
483 return loop(0)->Contains(p); // Optimization.
484 }
485 if (!bound_.Contains(p)) return false;
486 bool inside = false;
487 for (int i = 0; i < num_loops(); ++i) {
488 inside ^= loop(i)->Contains(p);
489 if (inside && !has_holes_) break; // Shells are disjoint.
490 }
491 return inside;
492}
493
494void S2Polygon::Encode(Encoder* const encoder) const {
495 encoder->Ensure(10); // Sufficient
496 encoder->put8(kCurrentEncodingVersionNumber);
497 encoder->put8(owns_loops_);
498 encoder->put8(has_holes_);
499 encoder->put32(loops_.size());
500 DCHECK_GE(encoder->avail(), 0);
501
502 for (int i = 0; i < num_loops(); ++i) {
503 loop(i)->Encode(encoder);
504 }
505 bound_.Encode(encoder);
506}
507
508bool S2Polygon::Decode(Decoder* const decoder) {
509 return DecodeInternal(decoder, false);
510}
511
512bool S2Polygon::DecodeWithinScope(Decoder* const decoder) {
513 return DecodeInternal(decoder, true);
514}
515
516bool S2Polygon::DecodeInternal(Decoder* const decoder, bool within_scope) {
517 unsigned char version = decoder->get8();
518 if (version > kCurrentEncodingVersionNumber) return false;
519
520 if (owns_loops_) DeleteLoopsInVector(&loops_);
521
522 owns_loops_ = decoder->get8();
523 has_holes_ = decoder->get8();
524 int num_loops = decoder->get32();
525 loops_.clear();
526 loops_.reserve(num_loops);
527 num_vertices_ = 0;
528 for (int i = 0; i < num_loops; ++i) {
529 loops_.push_back(new S2Loop);
530 if (within_scope) {
531 if (!loops_.back()->DecodeWithinScope(decoder)) return false;
532 } else {
533 if (!loops_.back()->Decode(decoder)) return false;
534 }
535 num_vertices_ += loops_.back()->num_vertices();
536 }
537 if (!bound_.Decode(decoder)) return false;
538
539 DCHECK(IsValid(loops_));
540
541 return decoder->avail() >= 0;
542}
543
544// Indexing structure to efficiently ClipEdge() of a polygon. This is
545// an abstract class because we need to use if for both polygons (for
546// InitToIntersection() and friends) and for sets of vectors of points
547// (for InitToSimplified()).
548//
549// Usage -- in your subclass:
550// - Call AddLoop() for each of your loops -- and keep them accessible in
551// your subclass.
552// - Overwrite EdgeFromTo(), calling DecodeIndex() and accessing your
553// underlying data with the resulting two indices.
554class S2LoopSequenceIndex: public S2EdgeIndex {
555 public:
556 S2LoopSequenceIndex(): num_edges_(0), num_loops_(0) {}
557 ~S2LoopSequenceIndex() {}
558
559 void AddLoop(int num_vertices) {
560 int vertices_so_far = num_edges_;
561 loop_to_first_index_.push_back(vertices_so_far);
562 index_to_loop_.resize(vertices_so_far + num_vertices);
563 for (int i = 0; i < num_vertices; ++i) {
564 index_to_loop_[vertices_so_far] = num_loops_;
565 vertices_so_far++;
566 }
567 num_edges_ += num_vertices;
568 num_loops_++;
569 }
570
571 inline void DecodeIndex(int index,
572 int* loop_index, int* vertex_in_loop) const {
573 *loop_index = index_to_loop_[index];
574 *vertex_in_loop = index - loop_to_first_index_[*loop_index];
575 }
576
577 // It is faster to return both vertices at once. It makes a difference
578 // for small polygons.
579 virtual void EdgeFromTo(int index,
580 S2Point const* * from, S2Point const* * to) const = 0;
581
582 int num_edges() const { return num_edges_; }
583
584 virtual S2Point const* edge_from(int index) const {
585 S2Point const* from;
586 S2Point const* to;
587 EdgeFromTo(index, &from, &to);
588 return from;
589 }
590
591 virtual S2Point const* edge_to(int index) const {
592 S2Point const* from;
593 S2Point const* to;
594 EdgeFromTo(index, &from, &to);
595 return to;
596 }
597
598 protected:
599 // Map from the unidimensional edge index to the loop this edge
600 // belongs to.
601 vector<int> index_to_loop_;
602
603 // Reverse of index_to_loop_: maps a loop index to the
604 // unidimensional index of the first edge in the loop.
605 vector<int> loop_to_first_index_;
606
607 // Total number of edges.
608 int num_edges_;
609
610 // Total number of loops.
611 int num_loops_;
612};
613
614// Indexing structure for an S2Polygon.
615class S2PolygonIndex: public S2LoopSequenceIndex {
616 public:
617 S2PolygonIndex(S2Polygon const* poly, bool reverse):
618 poly_(poly),
619 reverse_(reverse) {
620 for (int iloop = 0; iloop < poly_->num_loops(); ++iloop) {
621 AddLoop(poly_->loop(iloop)->num_vertices());
622 }
623 }
624
625 virtual ~S2PolygonIndex() {}
626
627 virtual void EdgeFromTo(int index,
628 S2Point const* * from, S2Point const* * to) const {
629 int loop_index;
630 int vertex_in_loop;
631 DecodeIndex(index, &loop_index, &vertex_in_loop);
632 S2Loop const* loop(poly_->loop(loop_index));
633 int from_index, to_index;
634 if (loop->is_hole() ^ reverse_) {
635 from_index = loop->num_vertices() - 1 - vertex_in_loop;
636 to_index = 2 * loop->num_vertices() - 2 - vertex_in_loop;
637 } else {
638 from_index = vertex_in_loop;
639 to_index = vertex_in_loop + 1;
640 }
641 *from = &(loop->vertex(from_index));
642 *to = &(loop->vertex(to_index));
643 }
644
645 private:
646 S2Polygon const* poly_;
647 bool reverse_;
648};
649
650// Indexing structure for a sequence of loops (not in the sense of S2:
651// the loops can self-intersect). Each loop is given in a vector<S2Point>
652// where, as in S2Loop, the last vertex is implicitely joined to the first.
653// Add each loop individually with AddVector(). The caller owns
654// the vector<S2Point>'s.
655class S2LoopsAsVectorsIndex: public S2LoopSequenceIndex {
656 public:
657 S2LoopsAsVectorsIndex() {}
658 ~S2LoopsAsVectorsIndex() {}
659
660 void AddVector(vector<S2Point> const* v) {
661 loops_.push_back(v);
662 AddLoop(v->size());
663 }
664
665 virtual void EdgeFromTo(int index,
666 S2Point const* *from, S2Point const* *to) const {
667 int loop_index;
668 int vertex_in_loop;
669 DecodeIndex(index, &loop_index, &vertex_in_loop);
670 vector<S2Point> const* loop = loops_[loop_index];
671 *from = &loop->at(vertex_in_loop);
672 *to = &loop->at(vertex_in_loop == loop->size() - 1
673 ? 0
674 : vertex_in_loop + 1);
675 }
676
677 private:
678 vector< vector<S2Point> const* > loops_;
679};
680
681typedef vector<pair<double, S2Point> > IntersectionSet;
682
683static void AddIntersection(S2Point const& a0, S2Point const& a1,
684 S2Point const& b0, S2Point const& b1,
685 bool add_shared_edges, int crossing,
686 IntersectionSet* intersections) {
687 if (crossing > 0) {
688 // There is a proper edge crossing.
689 S2Point x = S2EdgeUtil::GetIntersection(a0, a1, b0, b1);
690 double t = S2EdgeUtil::GetDistanceFraction(x, a0, a1);
691 intersections->push_back(make_pair(t, x));
692
693 } else if (S2EdgeUtil::VertexCrossing(a0, a1, b0, b1)) {
694 // There is a crossing at one of the vertices. The basic rule is simple:
695 // if a0 equals one of the "b" vertices, the crossing occurs at t=0;
696 // otherwise, it occurs at t=1.
697 //
698 // This has the effect that when two symmetric edges are
699 // encountered (an edge an its reverse), neither one is included
700 // in the output. When two duplicate edges are encountered, both
701 // are included in the output. The "add_shared_edges" flag allows
702 // one of these two copies to be removed by changing its
703 // intersection parameter from 0 to 1.
704
705 double t = (a0 == b0 || a0 == b1) ? 0 : 1;
706 if (!add_shared_edges && a1 == b1) t = 1;
707 intersections->push_back(make_pair(t, t == 0 ? a0 : a1));
708 }
709}
710
711static void ClipEdge(S2Point const& a0, S2Point const& a1,
712 S2LoopSequenceIndex* b_index,
713 bool add_shared_edges, IntersectionSet* intersections) {
714 // Find all points where the polygon B intersects the edge (a0,a1),
715 // and add the corresponding parameter values (in the range [0,1]) to
716 // "intersections".
717 S2LoopSequenceIndex::Iterator it(b_index);
718 it.GetCandidates(a0, a1);
719 S2EdgeUtil::EdgeCrosser crosser(&a0, &a1, &a0);
720 S2Point const* from;
721 S2Point const* to = NULL;
722 for (; !it.Done(); it.Next()) {
723 S2Point const* const previous_to = to;
724 b_index->EdgeFromTo(it.Index(), &from, &to);
725 if (previous_to != from) crosser.RestartAt(from);
726 int crossing = crosser.RobustCrossing(to);
727 if (crossing < 0) continue;
728 AddIntersection(a0, a1, *from, *to,
729 add_shared_edges, crossing, intersections);
730 }
731}
732
733
734static void ClipBoundary(S2Polygon const* a, bool reverse_a,
735 S2Polygon const* b, bool reverse_b, bool invert_b,
736 bool add_shared_edges, S2PolygonBuilder* builder) {
737 // Clip the boundary of A to the interior of B, and add the resulting edges
738 // to "builder". Shells are directed CCW and holes are directed clockwise,
739 // unless "reverse_a" or "reverse_b" is true in which case these directions
740 // in the corresponding polygon are reversed. If "invert_b" is true, the
741 // boundary of A is clipped to the exterior rather than the interior of B.
742 // If "add_shared_edges" is true, then the output will include any edges
743 // that are shared between A and B (both edges must be in the same direction
744 // after any edge reversals are taken into account).
745
746 S2PolygonIndex b_index(b, reverse_b);
747 b_index.PredictAdditionalCalls(a->num_vertices());
748
749 IntersectionSet intersections;
750 for (int i = 0; i < a->num_loops(); ++i) {
751 S2Loop* a_loop = a->loop(i);
752 int n = a_loop->num_vertices();
753 int dir = (a_loop->is_hole() ^ reverse_a) ? -1 : 1;
754 bool inside = b->Contains(a_loop->vertex(0)) ^ invert_b;
755 for (int j = (dir > 0) ? 0 : n; n > 0; --n, j += dir) {
756 S2Point const& a0 = a_loop->vertex(j);
757 S2Point const& a1 = a_loop->vertex(j + dir);
758 intersections.clear();
759 ClipEdge(a0, a1, &b_index, add_shared_edges, &intersections);
760
761 if (inside) intersections.push_back(make_pair(0, a0));
762 inside = (intersections.size() & 1);
763 DCHECK_EQ((b->Contains(a1) ^ invert_b), inside);
764 if (inside) intersections.push_back(make_pair(1, a1));
765 sort(intersections.begin(), intersections.end());
766 for (int k = 0; k < intersections.size(); k += 2) {
767 if (intersections[k] == intersections[k+1]) continue;
768 builder->AddEdge(intersections[k].second, intersections[k+1].second);
769 }
770 }
771 }
772}
773
774void S2Polygon::InitToIntersection(S2Polygon const* a, S2Polygon const* b) {
775 InitToIntersectionSloppy(a, b, S2EdgeUtil::kIntersectionTolerance);
776}
777
778void S2Polygon::InitToIntersectionSloppy(S2Polygon const* a, S2Polygon const* b,
779 S1Angle vertex_merge_radius) {
780 DCHECK_EQ(0, num_loops());
781 if (!a->bound_.Intersects(b->bound_)) return;
782
783 // We want the boundary of A clipped to the interior of B,
784 // plus the boundary of B clipped to the interior of A,
785 // plus one copy of any directed edges that are in both boundaries.
786
787 S2PolygonBuilderOptions options(S2PolygonBuilderOptions::DIRECTED_XOR());
788 options.set_vertex_merge_radius(vertex_merge_radius);
789 S2PolygonBuilder builder(options);
790 ClipBoundary(a, false, b, false, false, true, &builder);
791 ClipBoundary(b, false, a, false, false, false, &builder);
792 if (!builder.AssemblePolygon(this, NULL)) {
793 LOG(DFATAL) << "Bad directed edges in InitToIntersection";
794 }
795}
796
797void S2Polygon::InitToUnion(S2Polygon const* a, S2Polygon const* b) {
798 InitToUnionSloppy(a, b, S2EdgeUtil::kIntersectionTolerance);
799}
800
801void S2Polygon::InitToUnionSloppy(S2Polygon const* a, S2Polygon const* b,
802 S1Angle vertex_merge_radius) {
803 DCHECK_EQ(0, num_loops());
804
805 // We want the boundary of A clipped to the exterior of B,
806 // plus the boundary of B clipped to the exterior of A,
807 // plus one copy of any directed edges that are in both boundaries.
808
809 S2PolygonBuilderOptions options(S2PolygonBuilderOptions::DIRECTED_XOR());
810 options.set_vertex_merge_radius(vertex_merge_radius);
811 S2PolygonBuilder builder(options);
812 ClipBoundary(a, false, b, false, true, true, &builder);
813 ClipBoundary(b, false, a, false, true, false, &builder);
814 if (!builder.AssemblePolygon(this, NULL)) {
815 LOG(DFATAL) << "Bad directed edges";
816 }
817}
818
819void S2Polygon::InitToDifference(S2Polygon const* a, S2Polygon const* b) {
820 InitToDifferenceSloppy(a, b, S2EdgeUtil::kIntersectionTolerance);
821}
822
823void S2Polygon::InitToDifferenceSloppy(S2Polygon const* a, S2Polygon const* b,
824 S1Angle vertex_merge_radius) {
825 DCHECK_EQ(0, num_loops());
826
827 // We want the boundary of A clipped to the exterior of B,
828 // plus the reversed boundary of B clipped to the interior of A,
829 // plus one copy of any edge in A that is also a reverse edge in B.
830
831 S2PolygonBuilderOptions options(S2PolygonBuilderOptions::DIRECTED_XOR());
832 options.set_vertex_merge_radius(vertex_merge_radius);
833 S2PolygonBuilder builder(options);
834 ClipBoundary(a, false, b, true, true, true, &builder);
835 ClipBoundary(b, true, a, false, false, false, &builder);
836 if (!builder.AssemblePolygon(this, NULL)) {
837 LOG(DFATAL) << "Bad directed edges in InitToDifference";
838 }
839}
840
841// Takes a loop and simplifies it. This may return a self-intersecting
842// polyline. Always keeps the first vertex from the loop.
843vector<S2Point>* SimplifyLoopAsPolyline(S2Loop const* loop, S1Angle tolerance) {
844 vector<S2Point> points(loop->num_vertices() + 1);
845 // Add the first vertex at the beginning and at the end.
846 for (int i = 0; i <= loop->num_vertices(); ++i) {
847 points[i] = loop->vertex(i);
848 }
849 S2Polyline line(points);
850 vector<int> indices;
851 line.SubsampleVertices(tolerance, &indices);
852 if (indices.size() <= 2) return NULL;
853 // Add them all except the last: it is the same as the first.
854 vector<S2Point>* simplified_line = new vector<S2Point>(indices.size() - 1);
855 VLOG(4) << "Now simplified to: ";
856 for (int i = 0; i + 1 < indices.size(); ++i) {
857 (*simplified_line)[i] = line.vertex(indices[i]);
858 VLOG(4) << S2LatLng(line.vertex(indices[i]));
859 }
860 return simplified_line;
861}
862
863// Takes a set of possibly intersecting edges, stored in an
864// S2EdgeIndex. Breaks the edges into small pieces so that there is
865// no intersection anymore, and adds all these edges to the builder.
866void BreakEdgesAndAddToBuilder(S2LoopsAsVectorsIndex* edge_index,
867 S2PolygonBuilder* builder) {
868 // If there are self intersections, we add the pieces separately.
869 for (int i = 0; i < edge_index->num_edges(); ++i) {
870 S2Point const* from;
871 S2Point const* to;
872 edge_index->EdgeFromTo(i, &from, &to);
873
874 IntersectionSet intersections;
875 intersections.push_back(make_pair(0, *from));
876 // add_shared_edges can be false or true: it makes no difference
877 // due to the way we call ClipEdge.
878 ClipEdge(*from, *to, edge_index, false, &intersections);
879 intersections.push_back(make_pair(1, *to));
880 sort(intersections.begin(), intersections.end());
881 for (int k = 0; k + 1 < intersections.size(); ++k) {
882 if (intersections[k] == intersections[k+1]) continue;
883 builder->AddEdge(intersections[k].second, intersections[k+1].second);
884 }
885 }
886}
887
888// Simplifies the polygon. The algorithm is straightforward and naive:
889// 1. Simplify each loop by removing points while staying in the
890// tolerance zone. This uses S2Polyline::SubsampleVertices(),
891// which is not guaranteed to be optimal in terms of number of
892// points.
893// 2. Break any edge in pieces such that no piece intersects any
894// other.
895// 3. Use the polygon builder to regenerate the full polygon.
896void S2Polygon::InitToSimplified(S2Polygon const* a, S1Angle tolerance) {
897 S2PolygonBuilderOptions builder_options =
898 S2PolygonBuilderOptions::UNDIRECTED_XOR();
899 builder_options.set_validate(false);
900 // Ideally, we would want to set the vertex_merge_radius of the
901 // builder roughly to tolerance (and in fact forego the edge
902 // splitting step). Alas, if we do that, we are liable to the
903 // 'chain effect', where vertices are merged with closeby vertices
904 // and so on, so that a vertex can move by an arbitrary distance.
905 // So we remain conservative:
906 builder_options.set_vertex_merge_radius(tolerance * 0.10);
907 S2PolygonBuilder builder(builder_options);
908
909 // Simplify each loop separately and add to the edge index
910 S2LoopsAsVectorsIndex index;
911 vector<vector<S2Point>*> simplified_loops;
912 for (int i = 0; i < a->num_loops(); ++i) {
913 vector<S2Point>* simpler = SimplifyLoopAsPolyline(a->loop(i), tolerance);
914 if (NULL == simpler) continue;
915 simplified_loops.push_back(simpler);
916 index.AddVector(simpler);
917 }
918 if (0 != index.num_edges()) {
919 BreakEdgesAndAddToBuilder(&index, &builder);
920
921 if (!builder.AssemblePolygon(this, NULL)) {
922 LOG(DFATAL) << "Bad edges in InitToSimplified.";
923 }
924 }
925
926 for (int i = 0; i < simplified_loops.size(); ++i) {
927 delete simplified_loops[i];
928 }
929 simplified_loops.clear();
930}
931
932void S2Polygon::InternalClipPolyline(bool invert,
933 S2Polyline const* a,
934 vector<S2Polyline*> *out,
935 S1Angle merge_radius) const {
936 // Clip the polyline A to the interior of this polygon.
937 // The resulting polyline(s) will be appended to 'out'.
938 // If invert is true, we clip A to the exterior of this polygon instead.
939 // Vertices will be dropped such that adjacent vertices will not
940 // be closer than 'merge_radius'.
941 //
942 // We do the intersection/subtraction by walking the polyline edges.
943 // For each edge, we compute all intersections with the polygon boundary
944 // and sort them in increasing order of distance along that edge.
945 // We then divide the intersection points into pairs, and output a
946 // clipped polyline segment for each one.
947 // We keep track of whether we're inside or outside of the polygon at
948 // all times to decide which segments to output.
949
950 CHECK(out->empty());
951
952 IntersectionSet intersections;
953 vector<S2Point> vertices;
954 S2PolygonIndex poly_index(this, false);
955 int n = a->num_vertices();
956 bool inside = Contains(a->vertex(0)) ^ invert;
957 for (int j = 0; j < n-1; j++) {
958 S2Point const& a0 = a->vertex(j);
959 S2Point const& a1 = a->vertex(j + 1);
960 ClipEdge(a0, a1, &poly_index, true, &intersections);
961 if (inside) intersections.push_back(make_pair(0, a0));
962 inside = (intersections.size() & 1);
963 DCHECK_EQ((Contains(a1) ^ invert), inside);
964 if (inside) intersections.push_back(make_pair(1, a1));
965 sort(intersections.begin(), intersections.end());
966 // At this point we have a sorted array of vertex pairs representing
967 // the edge(s) obtained after clipping (a0,a1) against the polygon.
968 for (int k = 0; k < intersections.size(); k += 2) {
969 if (intersections[k] == intersections[k+1]) continue;
970 S2Point const& v0 = intersections[k].second;
971 S2Point const& v1 = intersections[k+1].second;
972
973 // If the gap from the previous vertex to this one is large
974 // enough, start a new polyline.
975 if (!vertices.empty() &&
976 vertices.back().Angle(v0) > merge_radius.radians()) {
977 out->push_back(new S2Polyline(vertices));
978 vertices.clear();
979 }
980 // Append this segment to the current polyline, ignoring any
981 // vertices that are too close to the previous vertex.
982 if (vertices.empty()) vertices.push_back(v0);
983 if (vertices.back().Angle(v1) > merge_radius.radians()) {
984 vertices.push_back(v1);
985 }
986 }
987 intersections.clear();
988 }
989 if (!vertices.empty()) {
990 out->push_back(new S2Polyline(vertices));
991 }
992}
993
994void S2Polygon::IntersectWithPolyline(
995 S2Polyline const* a,
996 vector<S2Polyline*> *out) const {
997 IntersectWithPolylineSloppy(a, out, S2EdgeUtil::kIntersectionTolerance);
998}
999
1000void S2Polygon::IntersectWithPolylineSloppy(
1001 S2Polyline const* a,
1002 vector<S2Polyline*> *out,
1003 S1Angle vertex_merge_radius) const {
1004 InternalClipPolyline(false, a, out, vertex_merge_radius);
1005}
1006
1007void S2Polygon::SubtractFromPolyline(S2Polyline const* a,
1008 vector<S2Polyline*> *out) const {
1009 SubtractFromPolylineSloppy(a, out, S2EdgeUtil::kIntersectionTolerance);
1010}
1011
1012void S2Polygon::SubtractFromPolylineSloppy(
1013 S2Polyline const* a,
1014 vector<S2Polyline*> *out,
1015 S1Angle vertex_merge_radius) const {
1016 InternalClipPolyline(true, a, out, vertex_merge_radius);
1017}
1018
1019
1020S2Polygon* S2Polygon::DestructiveUnion(vector<S2Polygon*>* polygons) {
1021 return DestructiveUnionSloppy(polygons, S2EdgeUtil::kIntersectionTolerance);
1022}
1023
1024S2Polygon* S2Polygon::DestructiveUnionSloppy(vector<S2Polygon*>* polygons,
1025 S1Angle vertex_merge_radius) {
1026 // Effectively create a priority queue of polygons in order of number of
1027 // vertices. Repeatedly union the two smallest polygons and add the result
1028 // to the queue until we have a single polygon to return.
1029 typedef multimap<int, S2Polygon*> QueueType;
1030 QueueType queue; // Map from # of vertices to polygon.
1031 for (int i = 0; i < polygons->size(); ++i)
1032 queue.insert(make_pair((*polygons)[i]->num_vertices(), (*polygons)[i]));
1033 polygons->clear();
1034
1035 while (queue.size() > 1) {
1036 // Pop two simplest polygons from queue.
1037 QueueType::iterator smallest_it = queue.begin();
1038 int a_size = smallest_it->first;
1039 S2Polygon* a_polygon = smallest_it->second;
1040 queue.erase(smallest_it);
1041 smallest_it = queue.begin();
1042 int b_size = smallest_it->first;
1043 S2Polygon* b_polygon = smallest_it->second;
1044 queue.erase(smallest_it);
1045
1046 // Union and add result back to queue.
1047 S2Polygon* union_polygon = new S2Polygon();
1048 union_polygon->InitToUnionSloppy(a_polygon, b_polygon, vertex_merge_radius);
1049 delete a_polygon;
1050 delete b_polygon;
1051 queue.insert(make_pair(a_size + b_size, union_polygon));
1052 // We assume that the number of vertices in the union polygon is the
1053 // sum of the number of vertices in the original polygons, which is not
1054 // always true, but will almost always be a decent approximation, and
1055 // faster than recomputing.
1056 }
1057
1058 if (queue.empty())
1059 return new S2Polygon();
1060 else
1061 return queue.begin()->second;
1062}
1063
1064void S2Polygon::InitToCellUnionBorder(S2CellUnion const& cells) {
1065 // Use a polygon builder to union the cells in the union. Due to rounding
1066 // errors, we can't do an exact union - when a small cell is adjacent to a
1067 // larger cell, the shared edges can fail to line up exactly. Two cell edges
1068 // cannot come closer then kMinWidth, so if we have the polygon builder merge
1069 // edges within half that distance, we should always merge shared edges
1070 // without merging different edges.
1071 S2PolygonBuilderOptions options(S2PolygonBuilderOptions::DIRECTED_XOR());
1072 double min_cell_angle = S2::kMinWidth.GetValue(S2CellId::kMaxLevel);
1073 options.set_vertex_merge_radius(S1Angle::Radians(min_cell_angle / 2));
1074 S2PolygonBuilder builder(options);
1075 for (int i = 0; i < cells.num_cells(); ++i) {
1076 S2Loop cell_loop(S2Cell(cells.cell_id(i)));
1077 builder.AddLoop(&cell_loop);
1078 }
1079 if (!builder.AssemblePolygon(this, NULL)) {
1080 LOG(DFATAL) << "AssemblePolygon failed in InitToCellUnionBorder";
1081 }
1082}
1083
1084bool S2Polygon::IsNormalized() const {
1085 set<S2Point> vertices;
1086 S2Loop* last_parent = NULL;
1087 for (int i = 0; i < num_loops(); ++i) {
1088 S2Loop* child = loop(i);
1089 if (child->depth() == 0) continue;
1090 S2Loop* parent = loop(GetParent(i));
1091 if (parent != last_parent) {
1092 vertices.clear();
1093 for (int j = 0; j < parent->num_vertices(); ++j) {
1094 vertices.insert(parent->vertex(j));
1095 }
1096 last_parent = parent;
1097 }
1098 int count = 0;
1099 for (int j = 0; j < child->num_vertices(); ++j) {
1100 if (vertices.count(child->vertex(j)) > 0) ++count;
1101 }
1102 if (count > 1) return false;
1103 }
1104 return true;
1105}
1106
1107bool S2Polygon::BoundaryEquals(S2Polygon const* b) const {
1108 if (num_loops() != b->num_loops()) return false;
1109
1110 for (int i = 0; i < num_loops(); ++i) {
1111 S2Loop* a_loop = loop(i);
1112 bool success = false;
1113 for (int j = 0; j < num_loops(); ++j) {
1114 S2Loop* b_loop = b->loop(j);
1115 if ((b_loop->depth() == a_loop->depth()) &&
1116 b_loop->BoundaryEquals(a_loop)) {
1117 success = true;
1118 break;
1119 }
1120 }
1121 if (!success) return false;
1122 }
1123 return true;
1124}
1125
1126bool S2Polygon::BoundaryApproxEquals(S2Polygon const* b,
1127 double max_error) const {
1128 if (num_loops() != b->num_loops()) return false;
1129
1130 // For now, we assume that there is at most one candidate match for each
1131 // loop. (So far this method is just used for testing.)
1132
1133 for (int i = 0; i < num_loops(); ++i) {
1134 S2Loop* a_loop = loop(i);
1135 bool success = false;
1136 for (int j = 0; j < num_loops(); ++j) {
1137 S2Loop* b_loop = b->loop(j);
1138 if (b_loop->depth() == a_loop->depth() &&
1139 b_loop->BoundaryApproxEquals(a_loop, max_error)) {
1140 success = true;
1141 break;
1142 }
1143 }
1144 if (!success) return false;
1145 }
1146 return true;
1147}
1148
1149bool S2Polygon::BoundaryNear(S2Polygon const* b, double max_error) const {
1150 if (num_loops() != b->num_loops()) return false;
1151
1152 // For now, we assume that there is at most one candidate match for each
1153 // loop. (So far this method is just used for testing.)
1154
1155 for (int i = 0; i < num_loops(); ++i) {
1156 S2Loop* a_loop = loop(i);
1157 bool success = false;
1158 for (int j = 0; j < num_loops(); ++j) {
1159 S2Loop* b_loop = b->loop(j);
1160 if (b_loop->depth() == a_loop->depth() &&
1161 b_loop->BoundaryNear(a_loop, max_error)) {
1162 success = true;
1163 break;
1164 }
1165 }
1166 if (!success) return false;
1167 }
1168 return true;
1169}
1170
1171S2Point S2Polygon::Project(S2Point const& point) const {
1172 DCHECK(!loops_.empty());
1173
1174 if (Contains(point)) {
1175 return point;
1176 }
1177
1178 S1Angle min_distance = S1Angle::Radians(10);
1179 int min_loop_index = 0;
1180 int min_vertex_index = 0;
1181
1182 for (int l = 0; l < num_loops(); ++l) {
1183 S2Loop *a_loop = loop(l);
1184 for (int v = 0; v < a_loop->num_vertices(); ++v) {
1185 S1Angle distance_to_segment =
1186 S2EdgeUtil::GetDistance(point,
1187 a_loop->vertex(v),
1188 a_loop->vertex(v + 1));
1189 if (distance_to_segment < min_distance) {
1190 min_distance = distance_to_segment;
1191 min_loop_index = l;
1192 min_vertex_index = v;
1193 }
1194 }
1195 }
1196
1197 S2Point closest_point = S2EdgeUtil::GetClosestPoint(
1198 point,
1199 loop(min_loop_index)->vertex(min_vertex_index),
1200 loop(min_loop_index)->vertex(min_vertex_index + 1));
1201
1202 return closest_point;
1203}
1204