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 <set>
10using std::set;
11using std::multiset;
12
13#include <vector>
14using std::vector;
15
16#include <hash_map>
17using __gnu_cxx::hash_map;
18
19#include <utility>
20using std::pair;
21using std::make_pair;
22
23
24#include "s2loop.h"
25
26// #include "base/commandlineflags.h"
27#include "base/logging.h"
28#include "base/scoped_ptr.h"
29#include "util/coding/coder.h"
30#include "s2cap.h"
31#include "s2cell.h"
32#include "s2edgeindex.h"
33
34static const unsigned char kCurrentEncodingVersionNumber = 1;
35
36// DECLARE_bool(s2debug); // defined in s2.cc
37
38S2Point const* S2LoopIndex::edge_from(int index) const {
39 return &loop_->vertex(index);
40}
41
42S2Point const* S2LoopIndex::edge_to(int index) const {
43 return &loop_->vertex(index+1);
44}
45
46int S2LoopIndex::num_edges() const {
47 return loop_->num_vertices();
48}
49
50S2Loop::S2Loop()
51 : num_vertices_(0),
52 vertices_(NULL),
53 owns_vertices_(false),
54 bound_(S2LatLngRect::Empty()),
55 depth_(0),
56 index_(this),
57 num_find_vertex_calls_(0) {
58}
59
60S2Loop::S2Loop(vector<S2Point> const& vertices)
61 : num_vertices_(0),
62 vertices_(NULL),
63 owns_vertices_(false),
64 bound_(S2LatLngRect::Full()),
65 depth_(0),
66 index_(this),
67 num_find_vertex_calls_(0) {
68 Init(vertices);
69}
70
71void S2Loop::ResetMutableFields() {
72 index_.Reset();
73 num_find_vertex_calls_ = 0;
74 vertex_to_index_.clear();
75}
76
77void S2Loop::Init(vector<S2Point> const& vertices) {
78 ResetMutableFields();
79
80 if (owns_vertices_) delete[] vertices_;
81 num_vertices_ = vertices.size();
82 if (vertices.empty()) {
83 vertices_ = NULL;
84 } else {
85 vertices_ = new S2Point[num_vertices_];
86 memcpy(vertices_, &vertices[0], num_vertices_ * sizeof(vertices_[0]));
87 }
88 owns_vertices_ = true;
89 bound_ = S2LatLngRect::Full();
90
91 // InitOrigin() must be called before InitBound() because the latter
92 // function expects Contains() to work properly.
93 InitOrigin();
94 InitBound();
95}
96
97bool S2Loop::IsValid() const {
98 // Loops must have at least 3 vertices.
99 if (num_vertices() < 3) {
100 VLOG(2) << "Degenerate loop";
101 return false;
102 }
103 // All vertices must be unit length.
104 for (int i = 0; i < num_vertices(); ++i) {
105 if (!S2::IsUnitLength(vertex(i))) {
106 VLOG(2) << "Vertex " << i << " is not unit length";
107 return false;
108 }
109 }
110 // Loops are not allowed to have any duplicate vertices.
111 hash_map<S2Point, int> vmap;
112 for (int i = 0; i < num_vertices(); ++i) {
113 if (!vmap.insert(make_pair(vertex(i), i)).second) {
114 VLOG(2) << "Duplicate vertices: " << vmap[vertex(i)] << " and " << i;
115 return false;
116 }
117 }
118 // Non-adjacent edges are not allowed to intersect.
119 bool crosses = false;
120 index_.PredictAdditionalCalls(num_vertices());
121 S2EdgeIndex::Iterator it(&index_);
122 for (int i = 0; i < num_vertices(); ++i) {
123 S2EdgeUtil::EdgeCrosser crosser(&vertex(i), &vertex(i+1), &vertex(0));
124 int previous_index = -2;
125 for (it.GetCandidates(vertex(i), vertex(i+1)); !it.Done(); it.Next()) {
126 int ai = it.Index();
127 // There is no need to test the same thing twice. Moreover, two edges
128 // that abut at ai+1 will have been tested for equality above.
129 if (ai > i+1) {
130 if (previous_index != ai) crosser.RestartAt(&vertex(ai));
131 // Beware, this may return the loop is valid if there is a
132 // "vertex crossing".
133 // TODO(user): Fix that.
134 crosses = crosser.RobustCrossing(&vertex(ai+1)) > 0;
135 previous_index = ai + 1;
136 if (crosses) {
137 VLOG(2) << "Edges " << i << " and " << ai << " cross";
138 // additional debugging information:
139 VLOG(2) << "Edge locations in degrees: "
140 << S2LatLng(vertex(i)) << "-" << S2LatLng(vertex(i+1))
141 << " and "
142 << S2LatLng(vertex(ai)) << "-" << S2LatLng(vertex(ai+1));
143 break;
144 }
145 }
146 }
147 if (crosses) break;
148 }
149
150 return !crosses;
151}
152
153bool S2Loop::IsValid(vector<S2Point> const& vertices, int max_adjacent) {
154 if (vertices.size() < 3) return false;
155
156 S2Loop loop(vertices);
157 return loop.IsValid();
158}
159
160bool S2Loop::IsValid(int max_adjacent) const {
161 return IsValid();
162}
163
164void S2Loop::InitOrigin() {
165 // The bounding box does not need to be correct before calling this
166 // function, but it must at least contain vertex(1) since we need to
167 // do a Contains() test on this point below.
168 DCHECK(bound_.Contains(vertex(1)));
169
170 // To ensure that every point is contained in exactly one face of a
171 // subdivision of the sphere, all containment tests are done by counting the
172 // edge crossings starting at a fixed point on the sphere (S2::Origin()).
173 // We need to know whether this point is inside or outside of the loop.
174 // We do this by first guessing that it is outside, and then seeing whether
175 // we get the correct containment result for vertex 1. If the result is
176 // incorrect, the origin must be inside the loop.
177 //
178 // A loop with consecutive vertices A,B,C contains vertex B if and only if
179 // the fixed vector R = S2::Ortho(B) is on the left side of the wedge ABC.
180 // The test below is written so that B is inside if C=R but not if A=R.
181
182 origin_inside_ = false; // Initialize before calling Contains().
183 bool v1_inside = S2::OrderedCCW(S2::Ortho(vertex(1)), vertex(0), vertex(2),
184 vertex(1));
185 if (v1_inside != Contains(vertex(1)))
186 origin_inside_ = true;
187}
188
189void S2Loop::InitBound() {
190 // The bounding rectangle of a loop is not necessarily the same as the
191 // bounding rectangle of its vertices. First, the loop may wrap entirely
192 // around the sphere (e.g. a loop that defines two revolutions of a
193 // candy-cane stripe). Second, the loop may include one or both poles.
194 // Note that a small clockwise loop near the equator contains both poles.
195
196 S2EdgeUtil::RectBounder bounder;
197 for (int i = 0; i <= num_vertices(); ++i) {
198 bounder.AddPoint(&vertex(i));
199 }
200 S2LatLngRect b = bounder.GetBound();
201 // Note that we need to initialize bound_ with a temporary value since
202 // Contains() does a bounding rectangle check before doing anything else.
203 bound_ = S2LatLngRect::Full();
204 if (Contains(S2Point(0, 0, 1))) {
205 b = S2LatLngRect(R1Interval(b.lat().lo(), M_PI_2), S1Interval::Full());
206 }
207 // If a loop contains the south pole, then either it wraps entirely
208 // around the sphere (full longitude range), or it also contains the
209 // north pole in which case b.lng().is_full() due to the test above.
210 // Either way, we only need to do the south pole containment test if
211 // b.lng().is_full().
212 if (b.lng().is_full() && Contains(S2Point(0, 0, -1))) {
213 b.mutable_lat()->set_lo(-M_PI_2);
214 }
215 bound_ = b;
216}
217
218S2Loop::S2Loop(S2Cell const& cell)
219 : bound_(cell.GetRectBound()),
220 index_(this),
221 num_find_vertex_calls_(0) {
222 num_vertices_ = 4;
223 vertices_ = new S2Point[num_vertices_];
224 depth_ = 0;
225 for (int i = 0; i < 4; ++i) {
226 vertices_[i] = cell.GetVertex(i);
227 }
228 owns_vertices_ = true;
229 InitOrigin();
230 InitBound();
231}
232
233S2Loop::~S2Loop() {
234 if (owns_vertices_) {
235 delete[] vertices_;
236 }
237}
238
239S2Loop::S2Loop(S2Loop const* src)
240 : num_vertices_(src->num_vertices_),
241 vertices_(new S2Point[num_vertices_]),
242 owns_vertices_(true),
243 bound_(src->bound_),
244 origin_inside_(src->origin_inside_),
245 depth_(src->depth_),
246 index_(this),
247 num_find_vertex_calls_(0) {
248 memcpy(vertices_, src->vertices_, num_vertices_ * sizeof(vertices_[0]));
249}
250
251S2Loop* S2Loop::Clone() const {
252 return new S2Loop(this);
253}
254
255int S2Loop::FindVertex(S2Point const& p) const {
256 num_find_vertex_calls_++;
257 if (num_vertices() < 10 || num_find_vertex_calls_ < 20) {
258 // Exhaustive search
259 for (int i = 1; i <= num_vertices(); ++i) {
260 if (vertex(i) == p) return i;
261 }
262 return -1;
263 }
264
265 if (vertex_to_index_.empty()) { // We haven't computed it yet.
266 for (int i = num_vertices(); i > 0; --i) {
267 vertex_to_index_[vertex(i)] = i;
268 }
269 }
270
271 map<S2Point, int>::const_iterator it;
272 it = vertex_to_index_.find(p);
273 if (it == vertex_to_index_.end()) return -1;
274 return it->second;
275}
276
277
278bool S2Loop::IsNormalized() const {
279 // Optimization: if the longitude span is less than 180 degrees, then the
280 // loop covers less than half the sphere and is therefore normalized.
281 if (bound_.lng().GetLength() < M_PI) return true;
282
283 // We allow some error so that hemispheres are always considered normalized.
284 // TODO(user): This might not be necessary once S2Polygon is enhanced so
285 // that it does not require its input loops to be normalized.
286 return GetTurningAngle() >= -1e-14;
287}
288
289void S2Loop::Normalize() {
290 CHECK(owns_vertices_);
291 if (!IsNormalized()) Invert();
292 DCHECK(IsNormalized());
293}
294
295void S2Loop::Invert() {
296 CHECK(owns_vertices_);
297
298 ResetMutableFields();
299 reverse(vertices_, vertices_ + num_vertices());
300 origin_inside_ ^= true;
301 if (bound_.lat().lo() > -M_PI_2 && bound_.lat().hi() < M_PI_2) {
302 // The complement of this loop contains both poles.
303 bound_ = S2LatLngRect::Full();
304 } else {
305 InitBound();
306 }
307}
308
309double S2Loop::GetArea() const {
310 double area = GetSurfaceIntegral(S2::SignedArea);
311 // The signed area should be between approximately -4*Pi and 4*Pi.
312 DCHECK_LE(fabs(area), 4 * M_PI + 1e-12);
313 if (area < 0) {
314 // We have computed the negative of the area of the loop exterior.
315 area += 4 * M_PI;
316 }
317 return max(0.0, min(4 * M_PI, area));
318}
319
320S2Point S2Loop::GetCentroid() const {
321 // GetSurfaceIntegral() returns either the integral of position over loop
322 // interior, or the negative of the integral of position over the loop
323 // exterior. But these two values are the same (!), because the integral of
324 // position over the entire sphere is (0, 0, 0).
325 return GetSurfaceIntegral(S2::TrueCentroid);
326}
327
328// Return (first, dir) such that first..first+n*dir are valid indices.
329int S2Loop::GetCanonicalFirstVertex(int* dir) const {
330 int first = 0;
331 int n = num_vertices();
332 for (int i = 1; i < n; ++i) {
333 if (vertex(i) < vertex(first)) first = i;
334 }
335 if (vertex(first + 1) < vertex(first + n - 1)) {
336 *dir = 1;
337 // 0 <= first <= n-1, so (first+n*dir) <= 2*n-1.
338 } else {
339 *dir = -1;
340 first += n;
341 // n <= first <= 2*n-1, so (first+n*dir) >= 0.
342 }
343 return first;
344}
345
346double S2Loop::GetTurningAngle() const {
347 // Don't crash even if the loop is not well-defined.
348 if (num_vertices() < 3) return 0;
349
350 // To ensure that we get the same result when the loop vertex order is
351 // rotated, and that we get the same result with the opposite sign when the
352 // vertices are reversed, we need to be careful to add up the individual
353 // turn angles in a consistent order. In general, adding up a set of
354 // numbers in a different order can change the sum due to rounding errors.
355 int n = num_vertices();
356 int dir, i = GetCanonicalFirstVertex(&dir);
357 double angle = S2::TurnAngle(vertex((i + n - dir) % n), vertex(i),
358 vertex((i + dir) % n));
359 while (--n > 0) {
360 i += dir;
361 angle += S2::TurnAngle(vertex(i - dir), vertex(i), vertex(i + dir));
362 }
363 return dir * angle;
364}
365
366S2Cap S2Loop::GetCapBound() const {
367 return bound_.GetCapBound();
368}
369
370bool S2Loop::Contains(S2Cell const& cell) const {
371 // A future optimization could also take advantage of the fact than an S2Cell
372 // is convex.
373
374 // It's not necessarily true that bound_.Contains(cell.GetRectBound())
375 // because S2Cell bounds are slightly conservative.
376 if (!bound_.Contains(cell.GetCenter())) return false;
377 S2Loop cell_loop(cell);
378 return Contains(&cell_loop);
379}
380
381bool S2Loop::MayIntersect(S2Cell const& cell) const {
382 // It is faster to construct a bounding rectangle for an S2Cell than for
383 // a general polygon. A future optimization could also take advantage of
384 // the fact than an S2Cell is convex.
385
386 if (!bound_.Intersects(cell.GetRectBound())) return false;
387 return S2Loop(cell).Intersects(this);
388}
389
390bool S2Loop::Contains(S2Point const& p) const {
391 if (!bound_.Contains(p)) return false;
392
393 bool inside = origin_inside_;
394 S2Point origin = S2::Origin();
395 S2EdgeUtil::EdgeCrosser crosser(&origin, &p, &vertex(0));
396
397 // The s2edgeindex library is not optimized yet for long edges,
398 // so the tradeoff to using it comes later.
399 if (num_vertices() < 2000) {
400 for (int i = 1; i <= num_vertices(); ++i) {
401 inside ^= crosser.EdgeOrVertexCrossing(&vertex(i));
402 }
403 return inside;
404 }
405
406 S2EdgeIndex::Iterator it(&index_);
407 int previous_index = -2;
408 for (it.GetCandidates(origin, p); !it.Done(); it.Next()) {
409 int ai = it.Index();
410 if (previous_index != ai - 1) crosser.RestartAt(&vertex(ai));
411 previous_index = ai;
412 inside ^= crosser.EdgeOrVertexCrossing(&vertex(ai+1));
413 }
414 return inside;
415}
416
417void S2Loop::Encode(Encoder* const encoder) const {
418 encoder->Ensure(num_vertices_ * sizeof(*vertices_) + 20); // sufficient
419
420 encoder->put8(kCurrentEncodingVersionNumber);
421 encoder->put32(num_vertices_);
422 encoder->putn(vertices_, sizeof(*vertices_) * num_vertices_);
423 encoder->put8(origin_inside_);
424 encoder->put32(depth_);
425 DCHECK_GE(encoder->avail(), 0);
426
427 bound_.Encode(encoder);
428}
429
430bool S2Loop::Decode(Decoder* const decoder) {
431 return DecodeInternal(decoder, false);
432}
433
434bool S2Loop::DecodeWithinScope(Decoder* const decoder) {
435 return DecodeInternal(decoder, true);
436}
437
438bool S2Loop::DecodeInternal(Decoder* const decoder,
439 bool within_scope) {
440 unsigned char version = decoder->get8();
441 if (version > kCurrentEncodingVersionNumber) return false;
442
443 num_vertices_ = decoder->get32();
444 if (owns_vertices_) delete[] vertices_;
445 if (within_scope) {
446 vertices_ = const_cast<S2Point *>(reinterpret_cast<S2Point const*>(
447 decoder->ptr()));
448 decoder->skip(num_vertices_ * sizeof(*vertices_));
449 owns_vertices_ = false;
450 } else {
451 vertices_ = new S2Point[num_vertices_];
452 decoder->getn(vertices_, num_vertices_ * sizeof(*vertices_));
453 owns_vertices_ = true;
454 }
455 origin_inside_ = decoder->get8();
456 depth_ = decoder->get32();
457 if (!bound_.Decode(decoder)) return false;
458
459 DCHECK(IsValid());
460
461 return decoder->avail() >= 0;
462}
463
464// This is a helper class for the AreBoundariesCrossing function.
465// Each time there is a point in common between the two loops passed
466// as parameters, the two associated wedges centered at this point are
467// passed to the ProcessWedge function of this processor. The function
468// updates an internal state based on the wedges, and returns true to
469// signal that no further processing is needed.
470//
471// To use AreBoundariesCrossing, subclass this class and keep an
472// internal state that you update each time ProcessWedge is called,
473// then query this internal state in the function that called
474// AreBoundariesCrossing.
475class WedgeProcessor {
476 public:
477 virtual ~WedgeProcessor() { }
478
479 virtual bool ProcessWedge(S2Point const& a0, S2Point const& ab1,
480 S2Point const& a2, S2Point const& b0,
481 S2Point const& b2) = 0;
482};
483
484bool S2Loop::AreBoundariesCrossing(
485 S2Loop const* b, WedgeProcessor* wedge_processor) const {
486 // See the header file for a description of what this method does.
487 index_.PredictAdditionalCalls(b->num_vertices());
488 S2EdgeIndex::Iterator it(&index_);
489 for (int j = 0; j < b->num_vertices(); ++j) {
490 S2EdgeUtil::EdgeCrosser crosser(&b->vertex(j), &b->vertex(j+1),
491 &b->vertex(0));
492 int previous_index = -2;
493 for (it.GetCandidates(b->vertex(j), b->vertex(j+1));
494 !it.Done(); it.Next()) {
495 int ai = it.Index();
496 if (previous_index != ai - 1) crosser.RestartAt(&vertex(ai));
497 previous_index = ai;
498 int crossing = crosser.RobustCrossing(&vertex(ai + 1));
499 if (crossing < 0) continue;
500 if (crossing > 0) return true;
501 // We only need to check each shared vertex once, so we only
502 // consider the case where vertex(i+1) == b->vertex(j+1).
503 if (vertex(ai+1) == b->vertex(j+1) &&
504 wedge_processor->ProcessWedge(vertex(ai), vertex(ai+1), vertex(ai+2),
505 b->vertex(j), b->vertex(j+2))) {
506 return false;
507 }
508 }
509 }
510 return false;
511}
512
513// WedgeProcessor to be used to check if loop A contains loop B.
514// DoesntContain() then returns true if there is a wedge of B not
515// contained in the associated wedge of A (and hence loop B is not
516// contained in loop A).
517class ContainsWedgeProcessor: public WedgeProcessor {
518 public:
519 ContainsWedgeProcessor(): doesnt_contain_(false) {}
520 bool DoesntContain() { return doesnt_contain_; }
521
522 protected:
523 virtual bool ProcessWedge(S2Point const& a0, S2Point const& ab1,
524 S2Point const& a2, S2Point const& b0,
525 S2Point const& b2) {
526 doesnt_contain_ = !S2EdgeUtil::WedgeContains(a0, ab1, a2, b0, b2);
527 return doesnt_contain_;
528 }
529
530 private:
531 bool doesnt_contain_;
532};
533
534bool S2Loop::Contains(S2Loop const* b) const {
535 // For this loop A to contains the given loop B, all of the following must
536 // be true:
537 //
538 // (1) There are no edge crossings between A and B except at vertices.
539 //
540 // (2) At every vertex that is shared between A and B, the local edge
541 // ordering implies that A contains B.
542 //
543 // (3) If there are no shared vertices, then A must contain a vertex of B
544 // and B must not contain a vertex of A. (An arbitrary vertex may be
545 // chosen in each case.)
546 //
547 // The second part of (3) is necessary to detect the case of two loops whose
548 // union is the entire sphere, i.e. two loops that contains each other's
549 // boundaries but not each other's interiors.
550
551 if (!bound_.Contains(b->bound_)) return false;
552
553 // Unless there are shared vertices, we need to check whether A contains a
554 // vertex of B. Since shared vertices are rare, it is more efficient to do
555 // this test up front as a quick rejection test.
556 if (!Contains(b->vertex(0)) && FindVertex(b->vertex(0)) < 0)
557 return false;
558
559 // Now check whether there are any edge crossings, and also check the loop
560 // relationship at any shared vertices.
561 ContainsWedgeProcessor p_wedge;
562 if (AreBoundariesCrossing(b, &p_wedge) || p_wedge.DoesntContain()) {
563 return false;
564 }
565
566 // At this point we know that the boundaries of A and B do not intersect,
567 // and that A contains a vertex of B. However we still need to check for
568 // the case mentioned above, where (A union B) is the entire sphere.
569 // Normally this check is very cheap due to the bounding box precondition.
570 if (bound_.Union(b->bound_).is_full()) {
571 if (b->Contains(vertex(0)) && b->FindVertex(vertex(0)) < 0) return false;
572 }
573 return true;
574}
575
576// WedgeProcessor to be used to check if loop A intersects loop B.
577// Intersects() then returns true when A and B have at least one pair
578// of associated wedges that intersect.
579class IntersectsWedgeProcessor: public WedgeProcessor {
580 public:
581 IntersectsWedgeProcessor(): intersects_(false) {}
582 bool Intersects() { return intersects_; }
583
584 protected:
585 virtual bool ProcessWedge(S2Point const& a0, S2Point const& ab1,
586 S2Point const& a2, S2Point const& b0,
587 S2Point const& b2) {
588 intersects_ = S2EdgeUtil::WedgeIntersects(a0, ab1, a2, b0, b2);
589 return intersects_;
590 }
591
592 private:
593 bool intersects_;
594};
595
596bool S2Loop::Intersects(S2Loop const* b) const {
597 // a->Intersects(b) if and only if !a->Complement()->Contains(b).
598 // This code is similar to Contains(), but is optimized for the case
599 // where both loops enclose less than half of the sphere.
600
601 // The largest of the two loops should be edgeindex'd.
602 if (b->num_vertices() > num_vertices()) return b->Intersects(this);
603
604 if (!bound_.Intersects(b->bound_)) return false;
605
606 // Unless there are shared vertices, we need to check whether A contains a
607 // vertex of B. Since shared vertices are rare, it is more efficient to do
608 // this test up front as a quick acceptance test.
609 if (Contains(b->vertex(0)) && FindVertex(b->vertex(0)) < 0)
610 return true;
611
612 // Now check whether there are any edge crossings, and also check the loop
613 // relationship at any shared vertices.
614 IntersectsWedgeProcessor p_wedge;
615 if (AreBoundariesCrossing(b, &p_wedge) || p_wedge.Intersects()) {
616 return true;
617 }
618
619 // We know that A does not contain a vertex of B, and that there are no edge
620 // crossings. Therefore the only way that A can intersect B is if B
621 // entirely contains A. We can check this by testing whether B contains an
622 // arbitrary non-shared vertex of A. Note that this check is usually cheap
623 // because of the bounding box precondition.
624 if (b->bound_.Contains(bound_)) {
625 if (b->Contains(vertex(0)) && b->FindVertex(vertex(0)) < 0) return true;
626 }
627 return false;
628}
629
630// WedgeProcessor to be used to check if the interior of loop A
631// contains the interior of loop B, or their boundaries cross each
632// other (therefore they have a proper intersection).
633// CrossesOrMayContain() then returns -1 if A crossed B, 0 if it is
634// not possible for A to contain B, and 1 otherwise.
635class ContainsOrCrossesProcessor: public WedgeProcessor {
636 public:
637 ContainsOrCrossesProcessor():
638 has_boundary_crossing_(false),
639 a_has_strictly_super_wedge_(false), b_has_strictly_super_wedge_(false),
640 has_disjoint_wedge_(false) {}
641
642 int CrossesOrMayContain() {
643 if (has_boundary_crossing_) return -1;
644 if (has_disjoint_wedge_ || b_has_strictly_super_wedge_) return 0;
645 return 1;
646 }
647
648 protected:
649 virtual bool ProcessWedge(S2Point const& a0, S2Point const& ab1,
650 S2Point const& a2, S2Point const& b0,
651 S2Point const& b2) {
652 const S2EdgeUtil::WedgeRelation wedge_relation =
653 S2EdgeUtil::GetWedgeRelation(a0, ab1, a2, b0, b2);
654 if (wedge_relation == S2EdgeUtil::WEDGE_PROPERLY_OVERLAPS) {
655 has_boundary_crossing_ = true;
656 return true;
657 }
658
659 a_has_strictly_super_wedge_ |=
660 (wedge_relation == S2EdgeUtil::WEDGE_PROPERLY_CONTAINS);
661 b_has_strictly_super_wedge_ |=
662 (wedge_relation == S2EdgeUtil::WEDGE_IS_PROPERLY_CONTAINED);
663 if (a_has_strictly_super_wedge_ && b_has_strictly_super_wedge_) {
664 has_boundary_crossing_ = true;
665 return true;
666 }
667
668 has_disjoint_wedge_ |= (wedge_relation == S2EdgeUtil::WEDGE_IS_DISJOINT);
669 return false;
670 }
671
672 private:
673 // True if any crossing on the boundary is discovered.
674 bool has_boundary_crossing_;
675 // True if A (B) has a strictly superwedge on a pair of wedges that
676 // share a common center point.
677 bool a_has_strictly_super_wedge_;
678 bool b_has_strictly_super_wedge_;
679 // True if there is a pair of disjoint wedges with common center
680 // point.
681 bool has_disjoint_wedge_;
682};
683
684int S2Loop::ContainsOrCrosses(S2Loop const* b) const {
685 // There can be containment or crossing only if the bounds intersect.
686 if (!bound_.Intersects(b->bound_)) return 0;
687
688 // Now check whether there are any edge crossings, and also check the loop
689 // relationship at any shared vertices. Note that unlike Contains() or
690 // Intersects(), we can't do a point containment test as a shortcut because
691 // we need to detect whether there are any edge crossings.
692 ContainsOrCrossesProcessor p_wedge;
693 if (AreBoundariesCrossing(b, &p_wedge)) {
694 return -1;
695 }
696 const int result = p_wedge.CrossesOrMayContain();
697 if (result <= 0) return result;
698
699 // At this point we know that the boundaries do not intersect, and we are
700 // given that (A union B) is a proper subset of the sphere. Furthermore
701 // either A contains B, or there are no shared vertices (due to the check
702 // above). So now we just need to distinguish the case where A contains B
703 // from the case where B contains A or the two loops are disjoint.
704 if (!bound_.Contains(b->bound_)) return 0;
705 if (!Contains(b->vertex(0)) && FindVertex(b->vertex(0)) < 0) return 0;
706 return 1;
707}
708
709bool S2Loop::ContainsNested(S2Loop const* b) const {
710 if (!bound_.Contains(b->bound_)) return false;
711
712 // We are given that A and B do not share any edges, and that either one
713 // loop contains the other or they do not intersect.
714 int m = FindVertex(b->vertex(1));
715 if (m < 0) {
716 // Since b->vertex(1) is not shared, we can check whether A contains it.
717 return Contains(b->vertex(1));
718 }
719 // Check whether the edge order around b->vertex(1) is compatible with
720 // A containing B.
721 return S2EdgeUtil::WedgeContains(vertex(m-1), vertex(m), vertex(m+1),
722 b->vertex(0), b->vertex(2));
723}
724
725bool S2Loop::BoundaryEquals(S2Loop const* b) const {
726 if (num_vertices() != b->num_vertices()) return false;
727 for (int offset = 0; offset < num_vertices(); ++offset) {
728 if (vertex(offset) == b->vertex(0)) {
729 // There is at most one starting offset since loop vertices are unique.
730 for (int i = 0; i < num_vertices(); ++i) {
731 if (vertex(i + offset) != b->vertex(i)) return false;
732 }
733 return true;
734 }
735 }
736 return false;
737}
738
739bool S2Loop::BoundaryApproxEquals(S2Loop const* b, double max_error) const {
740 if (num_vertices() != b->num_vertices()) return false;
741 for (int offset = 0; offset < num_vertices(); ++offset) {
742 if (S2::ApproxEquals(vertex(offset), b->vertex(0), max_error)) {
743 bool success = true;
744 for (int i = 0; i < num_vertices(); ++i) {
745 if (!S2::ApproxEquals(vertex(i + offset), b->vertex(i), max_error)) {
746 success = false;
747 break;
748 }
749 }
750 if (success) return true;
751 // Otherwise continue looping. There may be more than one candidate
752 // starting offset since vertices are only matched approximately.
753 }
754 }
755 return false;
756}
757
758static bool MatchBoundaries(S2Loop const* a, S2Loop const* b, int a_offset,
759 double max_error) {
760 // The state consists of a pair (i,j). A state transition consists of
761 // incrementing either "i" or "j". "i" can be incremented only if
762 // a(i+1+a_offset) is near the edge from b(j) to b(j+1), and a similar rule
763 // applies to "j". The function returns true iff we can proceed all the way
764 // around both loops in this way.
765 //
766 // Note that when "i" and "j" can both be incremented, sometimes only one
767 // choice leads to a solution. We handle this using a stack and
768 // backtracking. We also keep track of which states have already been
769 // explored to avoid duplicating work.
770
771 vector<pair<int, int> > pending;
772 set<pair<int, int> > done;
773 pending.push_back(make_pair(0, 0));
774 while (!pending.empty()) {
775 int i = pending.back().first;
776 int j = pending.back().second;
777 pending.pop_back();
778 if (i == a->num_vertices() && j == b->num_vertices()) {
779 return true;
780 }
781 done.insert(make_pair(i, j));
782
783 // If (i == na && offset == na-1) where na == a->num_vertices(), then
784 // then (i+1+offset) overflows the [0, 2*na-1] range allowed by vertex().
785 // So we reduce the range if necessary.
786 int io = i + a_offset;
787 if (io >= a->num_vertices()) io -= a->num_vertices();
788
789 if (i < a->num_vertices() && done.count(make_pair(i+1, j)) == 0 &&
790 S2EdgeUtil::GetDistance(a->vertex(io+1),
791 b->vertex(j),
792 b->vertex(j+1)).radians() <= max_error) {
793 pending.push_back(make_pair(i+1, j));
794 }
795 if (j < b->num_vertices() && done.count(make_pair(i, j+1)) == 0 &&
796 S2EdgeUtil::GetDistance(b->vertex(j+1),
797 a->vertex(io),
798 a->vertex(io+1)).radians() <= max_error) {
799 pending.push_back(make_pair(i, j+1));
800 }
801 }
802 return false;
803}
804
805bool S2Loop::BoundaryNear(S2Loop const* b, double max_error) const {
806 for (int a_offset = 0; a_offset < num_vertices(); ++a_offset) {
807 if (MatchBoundaries(this, b, a_offset, max_error)) return true;
808 }
809 return false;
810}
811