1
// Copyright 2005 Google Inc. All Rights Reserved.
9
#include "base/definer.h"
20
#include "s2polygon.h"
22
#include "base/port.h" // for HASH_NAMESPACE_DECLARATION_START
23
#include "util/coding/coder.h"
24
#include "s2edgeindex.h"
27
#include "s2cellunion.h"
28
#include "s2latlngrect.h"
29
#include "s2polygonbuilder.h"
30
#include "s2polyline.h"
32
static const unsigned char kCurrentEncodingVersionNumber = 1;
34
typedef pair<S2Point, S2Point> S2Edge;
36
S2Polygon::S2Polygon()
38
bound_(S2LatLngRect::Empty()),
44
S2Polygon::S2Polygon(vector<S2Loop*>* loops)
45
: bound_(S2LatLngRect::Empty()),
50
S2Polygon::S2Polygon(S2Cell const& cell)
51
: bound_(S2LatLngRect::Empty()),
55
S2Loop* loop = new S2Loop(cell);
56
bound_ = loop->GetRectBound();
57
loops_.push_back(loop);
60
S2Polygon::S2Polygon(S2Loop* loop)
61
: bound_(loop->GetRectBound()),
64
num_vertices_(loop->num_vertices()) {
65
loops_.push_back(loop);
68
void S2Polygon::Copy(S2Polygon const* src) {
69
DCHECK_EQ(0, num_loops());
70
for (int i = 0; i < src->num_loops(); ++i) {
71
loops_.push_back(src->loop(i)->Clone());
75
has_holes_ = src->has_holes_;
76
num_vertices_ = src->num_vertices();
79
S2Polygon* S2Polygon::Clone() const {
80
S2Polygon* result = new S2Polygon;
85
void S2Polygon::Release(vector<S2Loop*>* loops) {
87
loops->insert(loops->end(), loops_.begin(), loops_.end());
90
bound_ = S2LatLngRect::Empty();
94
static void DeleteLoopsInVector(vector<S2Loop*>* loops) {
95
for (size_t i = 0; i < loops->size(); ++i) {
101
S2Polygon::~S2Polygon() {
102
if (owns_loops_) DeleteLoopsInVector(&loops_);
105
typedef pair<S2Point, S2Point> S2PointPair;
108
template<> class hash<S2PointPair> {
110
size_t operator()(S2PointPair const& p) const {
112
return h(p.first) + (h(p.second) << 1);
117
bool S2Polygon::IsValid(const vector<S2Loop*>& loops) {
118
// If a loop contains an edge AB, then no other loop may contain AB or BA.
119
if (loops.size() > 1) {
120
hash_map<S2PointPair, pair<int, int> > edges;
121
for (size_t i = 0; i < loops.size(); ++i) {
122
S2Loop* lp = loops[i];
123
for (int j = 0; j < lp->num_vertices(); ++j) {
124
S2PointPair key = make_pair(lp->vertex(j), lp->vertex(j + 1));
125
if (edges.insert(make_pair(key, make_pair(i, j))).second) {
126
key = make_pair(lp->vertex(j + 1), lp->vertex(j));
127
if (edges.insert(make_pair(key, make_pair(i, j))).second)
130
pair<int, int> other = edges[key];
131
VLOG(2) << "Duplicate edge: loop " << i << ", edge " << j
132
<< " and loop " << other.first << ", edge " << other.second;
138
// Verify that no loop covers more than half of the sphere, and that no
140
for (size_t i = 0; i < loops.size(); ++i) {
141
if (!loops[i]->IsNormalized()) {
142
VLOG(2) << "Loop " << i << " encloses more than half the sphere";
145
for (size_t j = i + 1; j < loops.size(); ++j) {
146
// This test not only checks for edge crossings, it also detects
147
// cases where the two boundaries cross at a shared vertex.
148
if (loops[i]->ContainsOrCrosses(loops[j]) < 0) {
149
VLOG(2) << "Loop " << i << " crosses loop " << j;
158
bool S2Polygon::IsValid() const {
159
for (int i = 0; i < num_loops(); ++i) {
160
if (!loop(i)->IsValid()) {
164
return IsValid(loops_);
167
bool S2Polygon::IsValid(bool check_loops, int max_adjacent) const {
171
void S2Polygon::InsertLoop(S2Loop* new_loop, S2Loop* parent,
173
vector<S2Loop*>* children = &(*loop_map)[parent];
174
for (size_t i = 0; i < children->size(); ++i) {
175
S2Loop* child = (*children)[i];
176
if (child->ContainsNested(new_loop)) {
177
InsertLoop(new_loop, child, loop_map);
181
// No loop may contain the complement of another loop. (Handling this case
182
// is significantly more complicated.)
183
DCHECK(parent == NULL || !new_loop->ContainsNested(parent));
185
// Some of the children of the parent loop may now be children of
187
vector<S2Loop*>* new_children = &(*loop_map)[new_loop];
188
for (size_t i = 0; i < children->size();) {
189
S2Loop* child = (*children)[i];
190
if (new_loop->ContainsNested(child)) {
191
new_children->push_back(child);
192
children->erase(children->begin() + i);
197
children->push_back(new_loop);
200
void S2Polygon::InitLoop(S2Loop* loop, int depth, LoopMap* loop_map) {
202
loop->set_depth(depth);
203
loops_.push_back(loop);
205
vector<S2Loop*> const& children = (*loop_map)[loop];
206
for (size_t i = 0; i < children.size(); ++i) {
207
InitLoop(children[i], depth + 1, loop_map);
211
bool S2Polygon::ContainsChild(S2Loop* a, S2Loop* b, LoopMap const& loop_map) {
212
// This function is just used to verify that the loop map was
213
// constructed correctly.
215
if (a == b) return true;
216
vector<S2Loop*> const& children = loop_map.find(a)->second;
217
for (size_t i = 0; i < children.size(); ++i) {
218
if (ContainsChild(children[i], b, loop_map)) return true;
223
void S2Polygon::Init(vector<S2Loop*>* loops) {
225
CHECK(IsValid(*loops));
227
DCHECK(loops_.empty());
231
for (int i = 0; i < num_loops(); ++i) {
232
num_vertices_ += loop(i)->num_vertices();
236
for (int i = 0; i < num_loops(); ++i) {
237
InsertLoop(loop(i), NULL, &loop_map);
239
// Reorder the loops in depth-first traversal order.
241
InitLoop(NULL, -1, &loop_map);
244
// Check that the LoopMap is correct (this is fairly cheap).
245
for (int i = 0; i < num_loops(); ++i) {
246
for (int j = 0; j < num_loops(); ++j) {
247
if (i == j) continue;
248
CHECK_EQ(ContainsChild(loop(i), loop(j), loop_map),
249
loop(i)->ContainsNested(loop(j)));
254
// Compute the bounding rectangle of the entire polygon.
256
bound_ = S2LatLngRect::Empty();
257
for (int i = 0; i < num_loops(); ++i) {
258
if (loop(i)->sign() < 0) {
261
bound_ = bound_.Union(loop(i)->GetRectBound());
266
int S2Polygon::GetParent(int k) const {
267
int depth = loop(k)->depth();
268
if (depth == 0) return -1; // Optimization.
269
while (--k >= 0 && loop(k)->depth() >= depth) continue;
273
int S2Polygon::GetLastDescendant(int k) const {
274
if (k < 0) return num_loops() - 1;
275
int depth = loop(k)->depth();
276
while (++k < num_loops() && loop(k)->depth() > depth) continue;
280
double S2Polygon::GetArea() const {
282
for (int i = 0; i < num_loops(); ++i) {
283
area += loop(i)->sign() * loop(i)->GetArea();
288
S2Point S2Polygon::GetCentroid() const {
290
for (int i = 0; i < num_loops(); ++i) {
291
centroid += loop(i)->sign() * loop(i)->GetCentroid();
296
int S2Polygon::ContainsOrCrosses(S2Loop const* b) const {
298
for (int i = 0; i < num_loops(); ++i) {
299
int result = loop(i)->ContainsOrCrosses(b);
300
if (result < 0) return -1; // The loop boundaries intersect.
301
if (result > 0) inside ^= true;
303
return static_cast<int>(inside); // True if loop B is contained by the
307
bool S2Polygon::AnyLoopContains(S2Loop const* b) const {
308
// Return true if any loop contains the given loop.
309
for (int i = 0; i < num_loops(); ++i) {
310
if (loop(i)->Contains(b)) return true;
315
bool S2Polygon::ContainsAllShells(S2Polygon const* b) const {
316
// Return true if this polygon (A) contains all the shells of B.
317
for (int j = 0; j < b->num_loops(); ++j) {
318
if (b->loop(j)->sign() < 0) continue;
319
if (ContainsOrCrosses(b->loop(j)) <= 0) {
320
// Shell of B is not contained by A, or the boundaries intersect.
327
bool S2Polygon::ExcludesAllHoles(S2Polygon const* b) const {
328
// Return true if this polygon (A) excludes (i.e. does not intersect)
330
for (int j = 0; j < b->num_loops(); ++j) {
331
if (b->loop(j)->sign() > 0) continue;
332
if (ContainsOrCrosses(b->loop(j)) != 0) {
333
// Hole of B is contained by A, or the boundaries intersect.
340
bool S2Polygon::IntersectsAnyShell(S2Polygon const* b) const {
341
// Return true if this polygon (A) intersects any shell of B.
342
for (int j = 0; j < b->num_loops(); ++j) {
343
if (b->loop(j)->sign() < 0) continue;
344
if (IntersectsShell(b->loop(j)) != 0)
350
bool S2Polygon::IntersectsShell(S2Loop const* b) const {
352
for (int i = 0; i < num_loops(); ++i) {
353
if (loop(i)->Contains(b)) {
355
} else if (!b->Contains(loop(i)) && loop(i)->Intersects(b)) {
356
// We definitely have an intersection if the loops intersect AND one
357
// is not properly contained in the other. If A (this) is properly
358
// contained in a loop of B, we don't know yet if it may be actually
359
// inside a hole within B.
366
bool S2Polygon::Contains(S2Polygon const* b) const {
367
// If both polygons have one loop, use the more efficient S2Loop method.
368
// Note that S2Loop::Contains does its own bounding rectangle check.
369
if (num_loops() == 1 && b->num_loops() == 1) {
370
return loop(0)->Contains(b->loop(0));
373
// Otherwise if neither polygon has holes, we can still use the more
374
// efficient S2Loop::Contains method (rather than ContainsOrCrosses),
375
// but it's worthwhile to do our own bounds check first.
376
if (!bound_.Contains(b->bound_)) {
377
// If the union of the bounding boxes spans the full longitude range,
378
// it is still possible that polygon A contains B. (This is only
379
// possible if at least one polygon has multiple shells.)
380
if (!bound_.lng().Union(b->bound_.lng()).is_full()) return false;
382
if (!has_holes_ && !b->has_holes_) {
383
for (int j = 0; j < b->num_loops(); ++j) {
384
if (!AnyLoopContains(b->loop(j))) return false;
389
// This could be implemented more efficiently for polygons with lots of
390
// holes by keeping a copy of the LoopMap computed during initialization.
391
// However, in practice most polygons are one loop, and multiloop polygons
392
// tend to consist of many shells rather than holes. In any case, the real
393
// way to get more efficiency is to implement a sub-quadratic algorithm
394
// such as building a trapezoidal map.
396
// Every shell of B must be contained by an odd number of loops of A,
397
// and every hole of A must be contained by an even number of loops of B.
398
return ContainsAllShells(b) && b->ExcludesAllHoles(this);
401
bool S2Polygon::Intersects(S2Polygon const* b) const {
402
// A.Intersects(B) if and only if !Complement(A).Contains(B). However,
403
// implementing a Complement() operation is trickier than it sounds,
404
// and in any case it's more efficient to test for intersection directly.
406
// If both polygons have one loop, use the more efficient S2Loop method.
407
// Note that S2Loop::Intersects does its own bounding rectangle check.
408
if (num_loops() == 1 && b->num_loops() == 1) {
409
return loop(0)->Intersects(b->loop(0));
412
// Otherwise if neither polygon has holes, we can still use the more
413
// efficient S2Loop::Intersects method. The polygons intersect if and
414
// only if some pair of loop regions intersect.
415
if (!bound_.Intersects(b->bound_)) return false;
416
if (!has_holes_ && !b->has_holes_) {
417
for (int i = 0; i < num_loops(); ++i) {
418
for (int j = 0; j < b->num_loops(); ++j) {
419
if (loop(i)->Intersects(b->loop(j))) return true;
425
// Otherwise if any shell of B is contained by an odd number of loops of A,
426
// or any shell of A is contained by an odd number of loops of B, or there is
427
// an intersection without containment, then there is an intersection.
428
return IntersectsAnyShell(b) || b->IntersectsAnyShell(this);
431
S2Cap S2Polygon::GetCapBound() const {
432
return bound_.GetCapBound();
435
bool S2Polygon::Contains(S2Cell const& cell) const {
436
if (num_loops() == 1) {
437
return loop(0)->Contains(cell);
440
// We can't check bound_.Contains(cell.GetRectBound()) because S2Cell's
441
// GetRectBound() calculation is not precise.
442
if (!bound_.Contains(cell.GetCenter())) return false;
444
S2Loop cell_loop(cell);
445
S2Polygon cell_poly(&cell_loop);
446
bool contains = Contains(&cell_poly);
448
DCHECK(Contains(cell.GetCenter()));
453
bool S2Polygon::ApproxContains(S2Polygon const* b,
454
S1Angle vertex_merge_radius) const {
455
S2Polygon difference;
456
difference.InitToDifferenceSloppy(b, this, vertex_merge_radius);
457
return difference.num_loops() == 0;
460
bool S2Polygon::MayIntersect(S2Cell const& cell) const {
461
if (num_loops() == 1) {
462
return loop(0)->MayIntersect(cell);
464
if (!bound_.Intersects(cell.GetRectBound())) return false;
466
S2Loop cell_loop(cell);
467
S2Polygon cell_poly(&cell_loop);
468
bool intersects = Intersects(&cell_poly);
470
DCHECK(!Contains(cell.GetCenter()));
475
bool S2Polygon::VirtualContainsPoint(S2Point const& p) const {
476
return Contains(p); // The same as Contains() below, just virtual.
479
bool S2Polygon::Contains(S2Point const& p) const {
480
if (num_loops() == 1) {
481
return loop(0)->Contains(p); // Optimization.
483
if (!bound_.Contains(p)) return false;
485
for (int i = 0; i < num_loops(); ++i) {
486
inside ^= loop(i)->Contains(p);
487
if (inside && !has_holes_) break; // Shells are disjoint.
492
void S2Polygon::Encode(Encoder* const encoder) const {
493
encoder->Ensure(10); // Sufficient
494
encoder->put8(kCurrentEncodingVersionNumber);
495
encoder->put8(owns_loops_);
496
encoder->put8(has_holes_);
497
encoder->put32(loops_.size());
498
DCHECK_GE(encoder->avail(), 0);
500
for (int i = 0; i < num_loops(); ++i) {
501
loop(i)->Encode(encoder);
503
bound_.Encode(encoder);
506
bool S2Polygon::Decode(Decoder* const decoder) {
507
return DecodeInternal(decoder, false);
510
bool S2Polygon::DecodeWithinScope(Decoder* const decoder) {
511
return DecodeInternal(decoder, true);
514
bool S2Polygon::DecodeInternal(Decoder* const decoder, bool within_scope) {
515
unsigned char version = decoder->get8();
516
if (version > kCurrentEncodingVersionNumber) return false;
518
if (owns_loops_) DeleteLoopsInVector(&loops_);
520
owns_loops_ = decoder->get8();
521
has_holes_ = decoder->get8();
522
int num_loops = decoder->get32();
524
loops_.reserve(num_loops);
526
for (int i = 0; i < num_loops; ++i) {
527
loops_.push_back(new S2Loop);
529
if (!loops_.back()->DecodeWithinScope(decoder)) return false;
531
if (!loops_.back()->Decode(decoder)) return false;
533
num_vertices_ += loops_.back()->num_vertices();
535
if (!bound_.Decode(decoder)) return false;
537
DCHECK(IsValid(loops_));
539
return decoder->avail() >= 0;
542
// Indexing structure to efficiently ClipEdge() of a polygon. This is
543
// an abstract class because we need to use if for both polygons (for
544
// InitToIntersection() and friends) and for sets of vectors of points
545
// (for InitToSimplified()).
547
// Usage -- in your subclass:
548
// - Call AddLoop() for each of your loops -- and keep them accessible in
550
// - Overwrite EdgeFromTo(), calling DecodeIndex() and accessing your
551
// underlying data with the resulting two indices.
552
class S2LoopSequenceIndex: public S2EdgeIndex {
554
S2LoopSequenceIndex(): num_edges_(0), num_loops_(0) {}
555
~S2LoopSequenceIndex() {}
557
void AddLoop(int num_vertices) {
558
int vertices_so_far = num_edges_;
559
loop_to_first_index_.push_back(vertices_so_far);
560
index_to_loop_.resize(vertices_so_far + num_vertices);
561
for (int i = 0; i < num_vertices; ++i) {
562
index_to_loop_[vertices_so_far] = num_loops_;
565
num_edges_ += num_vertices;
569
inline void DecodeIndex(int index,
570
int* loop_index, int* vertex_in_loop) const {
571
*loop_index = index_to_loop_[index];
572
*vertex_in_loop = index - loop_to_first_index_[*loop_index];
575
// It is faster to return both vertices at once. It makes a difference
576
// for small polygons.
577
virtual void EdgeFromTo(int index,
578
S2Point const* * from, S2Point const* * to) const = 0;
580
int num_edges() const { return num_edges_; }
582
virtual S2Point const* edge_from(int index) const {
585
EdgeFromTo(index, &from, &to);
589
virtual S2Point const* edge_to(int index) const {
592
EdgeFromTo(index, &from, &to);
597
// Map from the unidimensional edge index to the loop this edge
599
vector<int> index_to_loop_;
601
// Reverse of index_to_loop_: maps a loop index to the
602
// unidimensional index of the first edge in the loop.
603
vector<int> loop_to_first_index_;
605
// Total number of edges.
608
// Total number of loops.
612
// Indexing structure for an S2Polygon.
613
class S2PolygonIndex: public S2LoopSequenceIndex {
615
S2PolygonIndex(S2Polygon const* poly, bool reverse):
618
for (int iloop = 0; iloop < poly_->num_loops(); ++iloop) {
619
AddLoop(poly_->loop(iloop)->num_vertices());
623
virtual ~S2PolygonIndex() {}
625
virtual void EdgeFromTo(int index,
626
S2Point const* * from, S2Point const* * to) const {
629
DecodeIndex(index, &loop_index, &vertex_in_loop);
630
S2Loop const* loop(poly_->loop(loop_index));
631
int from_index, to_index;
632
if (loop->is_hole() ^ reverse_) {
633
from_index = loop->num_vertices() - 1 - vertex_in_loop;
634
to_index = 2 * loop->num_vertices() - 2 - vertex_in_loop;
636
from_index = vertex_in_loop;
637
to_index = vertex_in_loop + 1;
639
*from = &(loop->vertex(from_index));
640
*to = &(loop->vertex(to_index));
644
S2Polygon const* poly_;
648
// Indexing structure for a sequence of loops (not in the sense of S2:
649
// the loops can self-intersect). Each loop is given in a vector<S2Point>
650
// where, as in S2Loop, the last vertex is implicitely joined to the first.
651
// Add each loop individually with AddVector(). The caller owns
652
// the vector<S2Point>'s.
653
class S2LoopsAsVectorsIndex: public S2LoopSequenceIndex {
655
S2LoopsAsVectorsIndex() {}
656
~S2LoopsAsVectorsIndex() {}
658
void AddVector(vector<S2Point> const* v) {
663
virtual void EdgeFromTo(int index,
664
S2Point const* *from, S2Point const* *to) const {
667
DecodeIndex(index, &loop_index, &vertex_in_loop);
668
vector<S2Point> const* loop = loops_[loop_index];
669
*from = &loop->at(vertex_in_loop);
670
*to = &loop->at((size_t)vertex_in_loop == loop->size() - 1
672
: vertex_in_loop + 1);
676
vector< vector<S2Point> const* > loops_;
679
typedef vector<pair<double, S2Point> > IntersectionSet;
681
static void AddIntersection(S2Point const& a0, S2Point const& a1,
682
S2Point const& b0, S2Point const& b1,
683
bool add_shared_edges, int crossing,
684
IntersectionSet* intersections) {
686
// There is a proper edge crossing.
687
S2Point x = S2EdgeUtil::GetIntersection(a0, a1, b0, b1);
688
double t = S2EdgeUtil::GetDistanceFraction(x, a0, a1);
689
intersections->push_back(make_pair(t, x));
691
} else if (S2EdgeUtil::VertexCrossing(a0, a1, b0, b1)) {
692
// There is a crossing at one of the vertices. The basic rule is simple:
693
// if a0 equals one of the "b" vertices, the crossing occurs at t=0;
694
// otherwise, it occurs at t=1.
696
// This has the effect that when two symmetric edges are
697
// encountered (an edge an its reverse), neither one is included
698
// in the output. When two duplicate edges are encountered, both
699
// are included in the output. The "add_shared_edges" flag allows
700
// one of these two copies to be removed by changing its
701
// intersection parameter from 0 to 1.
703
double t = (a0 == b0 || a0 == b1) ? 0 : 1;
704
if (!add_shared_edges && a1 == b1) t = 1;
705
intersections->push_back(make_pair(t, t == 0 ? a0 : a1));
709
static void ClipEdge(S2Point const& a0, S2Point const& a1,
710
S2LoopSequenceIndex* b_index,
711
bool add_shared_edges, IntersectionSet* intersections) {
712
// Find all points where the polygon B intersects the edge (a0,a1),
713
// and add the corresponding parameter values (in the range [0,1]) to
715
S2LoopSequenceIndex::Iterator it(b_index);
716
it.GetCandidates(a0, a1);
717
S2EdgeUtil::EdgeCrosser crosser(&a0, &a1, &a0);
719
S2Point const* to = NULL;
720
for (; !it.Done(); it.Next()) {
721
S2Point const* const previous_to = to;
722
b_index->EdgeFromTo(it.Index(), &from, &to);
723
if (previous_to != from) crosser.RestartAt(from);
724
int crossing = crosser.RobustCrossing(to);
725
if (crossing < 0) continue;
726
AddIntersection(a0, a1, *from, *to,
727
add_shared_edges, crossing, intersections);
732
static void ClipBoundary(S2Polygon const* a, bool reverse_a,
733
S2Polygon const* b, bool reverse_b, bool invert_b,
734
bool add_shared_edges, S2PolygonBuilder* builder) {
735
// Clip the boundary of A to the interior of B, and add the resulting edges
736
// to "builder". Shells are directed CCW and holes are directed clockwise,
737
// unless "reverse_a" or "reverse_b" is true in which case these directions
738
// in the corresponding polygon are reversed. If "invert_b" is true, the
739
// boundary of A is clipped to the exterior rather than the interior of B.
740
// If "add_shared_edges" is true, then the output will include any edges
741
// that are shared between A and B (both edges must be in the same direction
742
// after any edge reversals are taken into account).
744
S2PolygonIndex b_index(b, reverse_b);
745
b_index.PredictAdditionalCalls(a->num_vertices());
747
IntersectionSet intersections;
748
for (int i = 0; i < a->num_loops(); ++i) {
749
S2Loop* a_loop = a->loop(i);
750
int n = a_loop->num_vertices();
751
int dir = (a_loop->is_hole() ^ reverse_a) ? -1 : 1;
752
bool inside = b->Contains(a_loop->vertex(0)) ^ invert_b;
753
for (int j = (dir > 0) ? 0 : n; n > 0; --n, j += dir) {
754
S2Point const& a0 = a_loop->vertex(j);
755
S2Point const& a1 = a_loop->vertex(j + dir);
756
intersections.clear();
757
ClipEdge(a0, a1, &b_index, add_shared_edges, &intersections);
759
if (inside) intersections.push_back(make_pair(0, a0));
760
inside = (intersections.size() & 1);
761
DCHECK_EQ((b->Contains(a1) ^ invert_b), inside);
762
if (inside) intersections.push_back(make_pair(1, a1));
763
sort(intersections.begin(), intersections.end());
764
for (size_t k = 0; k < intersections.size(); k += 2) {
765
if (intersections[k] == intersections[k+1]) continue;
766
builder->AddEdge(intersections[k].second, intersections[k+1].second);
772
void S2Polygon::InitToIntersection(S2Polygon const* a, S2Polygon const* b) {
773
InitToIntersectionSloppy(a, b, S2EdgeUtil::kIntersectionTolerance);
776
void S2Polygon::InitToIntersectionSloppy(S2Polygon const* a, S2Polygon const* b,
777
S1Angle vertex_merge_radius) {
778
DCHECK_EQ(0, num_loops());
779
if (!a->bound_.Intersects(b->bound_)) return;
781
// We want the boundary of A clipped to the interior of B,
782
// plus the boundary of B clipped to the interior of A,
783
// plus one copy of any directed edges that are in both boundaries.
785
S2PolygonBuilderOptions options(S2PolygonBuilderOptions::DIRECTED_XOR());
786
options.set_vertex_merge_radius(vertex_merge_radius);
787
S2PolygonBuilder builder(options);
788
ClipBoundary(a, false, b, false, false, true, &builder);
789
ClipBoundary(b, false, a, false, false, false, &builder);
790
if (!builder.AssemblePolygon(this, NULL)) {
791
S2LOG(DFATAL) << "Bad directed edges in InitToIntersection";
795
void S2Polygon::InitToUnion(S2Polygon const* a, S2Polygon const* b) {
796
InitToUnionSloppy(a, b, S2EdgeUtil::kIntersectionTolerance);
799
void S2Polygon::InitToUnionSloppy(S2Polygon const* a, S2Polygon const* b,
800
S1Angle vertex_merge_radius) {
801
DCHECK_EQ(0, num_loops());
803
// We want the boundary of A clipped to the exterior of B,
804
// plus the boundary of B clipped to the exterior of A,
805
// plus one copy of any directed edges that are in both boundaries.
807
S2PolygonBuilderOptions options(S2PolygonBuilderOptions::DIRECTED_XOR());
808
options.set_vertex_merge_radius(vertex_merge_radius);
809
S2PolygonBuilder builder(options);
810
ClipBoundary(a, false, b, false, true, true, &builder);
811
ClipBoundary(b, false, a, false, true, false, &builder);
812
if (!builder.AssemblePolygon(this, NULL)) {
813
S2LOG(DFATAL) << "Bad directed edges";
817
void S2Polygon::InitToDifference(S2Polygon const* a, S2Polygon const* b) {
818
InitToDifferenceSloppy(a, b, S2EdgeUtil::kIntersectionTolerance);
821
void S2Polygon::InitToDifferenceSloppy(S2Polygon const* a, S2Polygon const* b,
822
S1Angle vertex_merge_radius) {
823
DCHECK_EQ(0, num_loops());
825
// We want the boundary of A clipped to the exterior of B,
826
// plus the reversed boundary of B clipped to the interior of A,
827
// plus one copy of any edge in A that is also a reverse edge in B.
829
S2PolygonBuilderOptions options(S2PolygonBuilderOptions::DIRECTED_XOR());
830
options.set_vertex_merge_radius(vertex_merge_radius);
831
S2PolygonBuilder builder(options);
832
ClipBoundary(a, false, b, true, true, true, &builder);
833
ClipBoundary(b, true, a, false, false, false, &builder);
834
if (!builder.AssemblePolygon(this, NULL)) {
835
S2LOG(DFATAL) << "Bad directed edges in InitToDifference";
839
// Takes a loop and simplifies it. This may return a self-intersecting
840
// polyline. Always keeps the first vertex from the loop.
841
vector<S2Point>* SimplifyLoopAsPolyline(S2Loop const* loop, S1Angle tolerance) {
842
vector<S2Point> points(loop->num_vertices() + 1);
843
// Add the first vertex at the beginning and at the end.
844
for (int i = 0; i <= loop->num_vertices(); ++i) {
845
points[i] = loop->vertex(i);
847
S2Polyline line(points);
849
line.SubsampleVertices(tolerance, &indices);
850
if (indices.size() <= 2) return NULL;
851
// Add them all except the last: it is the same as the first.
852
vector<S2Point>* simplified_line = new vector<S2Point>(indices.size() - 1);
853
VLOG(4) << "Now simplified to: ";
854
for (size_t i = 0; i + 1 < indices.size(); ++i) {
855
(*simplified_line)[i] = line.vertex(indices[i]);
856
VLOG(4) << S2LatLng(line.vertex(indices[i]));
858
return simplified_line;
861
// Takes a set of possibly intersecting edges, stored in an
862
// S2EdgeIndex. Breaks the edges into small pieces so that there is
863
// no intersection anymore, and adds all these edges to the builder.
864
void BreakEdgesAndAddToBuilder(S2LoopsAsVectorsIndex* edge_index,
865
S2PolygonBuilder* builder) {
866
// If there are self intersections, we add the pieces separately.
867
for (int i = 0; i < edge_index->num_edges(); ++i) {
870
edge_index->EdgeFromTo(i, &from, &to);
872
IntersectionSet intersections;
873
intersections.push_back(make_pair(0, *from));
874
// add_shared_edges can be false or true: it makes no difference
875
// due to the way we call ClipEdge.
876
ClipEdge(*from, *to, edge_index, false, &intersections);
877
intersections.push_back(make_pair(1, *to));
878
sort(intersections.begin(), intersections.end());
879
for (size_t k = 0; k + 1 < intersections.size(); ++k) {
880
if (intersections[k] == intersections[k+1]) continue;
881
builder->AddEdge(intersections[k].second, intersections[k+1].second);
886
// Simplifies the polygon. The algorithm is straightforward and naive:
887
// 1. Simplify each loop by removing points while staying in the
888
// tolerance zone. This uses S2Polyline::SubsampleVertices(),
889
// which is not guaranteed to be optimal in terms of number of
891
// 2. Break any edge in pieces such that no piece intersects any
893
// 3. Use the polygon builder to regenerate the full polygon.
894
void S2Polygon::InitToSimplified(S2Polygon const* a, S1Angle tolerance) {
895
S2PolygonBuilderOptions builder_options =
896
S2PolygonBuilderOptions::UNDIRECTED_XOR();
897
builder_options.set_validate(false);
898
// Ideally, we would want to set the vertex_merge_radius of the
899
// builder roughly to tolerance (and in fact forego the edge
900
// splitting step). Alas, if we do that, we are liable to the
901
// 'chain effect', where vertices are merged with closeby vertices
902
// and so on, so that a vertex can move by an arbitrary distance.
903
// So we remain conservative:
904
builder_options.set_vertex_merge_radius(tolerance * 0.10);
905
S2PolygonBuilder builder(builder_options);
907
// Simplify each loop separately and add to the edge index
908
S2LoopsAsVectorsIndex index;
909
vector<vector<S2Point>*> simplified_loops;
910
for (int i = 0; i < a->num_loops(); ++i) {
911
vector<S2Point>* simpler = SimplifyLoopAsPolyline(a->loop(i), tolerance);
912
if (NULL == simpler) continue;
913
simplified_loops.push_back(simpler);
914
index.AddVector(simpler);
916
if (0 != index.num_edges()) {
917
BreakEdgesAndAddToBuilder(&index, &builder);
919
if (!builder.AssemblePolygon(this, NULL)) {
920
S2LOG(DFATAL) << "Bad edges in InitToSimplified.";
924
for (size_t i = 0; i < simplified_loops.size(); ++i) {
925
delete simplified_loops[i];
927
simplified_loops.clear();
930
void S2Polygon::InternalClipPolyline(bool invert,
932
vector<S2Polyline*> *out,
933
S1Angle merge_radius) const {
934
// Clip the polyline A to the interior of this polygon.
935
// The resulting polyline(s) will be appended to 'out'.
936
// If invert is true, we clip A to the exterior of this polygon instead.
937
// Vertices will be dropped such that adjacent vertices will not
938
// be closer than 'merge_radius'.
940
// We do the intersection/subtraction by walking the polyline edges.
941
// For each edge, we compute all intersections with the polygon boundary
942
// and sort them in increasing order of distance along that edge.
943
// We then divide the intersection points into pairs, and output a
944
// clipped polyline segment for each one.
945
// We keep track of whether we're inside or outside of the polygon at
946
// all times to decide which segments to output.
950
IntersectionSet intersections;
951
vector<S2Point> vertices;
952
S2PolygonIndex poly_index(this, false);
953
int n = a->num_vertices();
954
bool inside = Contains(a->vertex(0)) ^ invert;
955
for (int j = 0; j < n-1; j++) {
956
S2Point const& a0 = a->vertex(j);
957
S2Point const& a1 = a->vertex(j + 1);
958
ClipEdge(a0, a1, &poly_index, true, &intersections);
959
if (inside) intersections.push_back(make_pair(0, a0));
960
inside = (intersections.size() & 1);
961
DCHECK_EQ((Contains(a1) ^ invert), inside);
962
if (inside) intersections.push_back(make_pair(1, a1));
963
sort(intersections.begin(), intersections.end());
964
// At this point we have a sorted array of vertex pairs representing
965
// the edge(s) obtained after clipping (a0,a1) against the polygon.
966
for (size_t k = 0; k < intersections.size(); k += 2) {
967
if (intersections[k] == intersections[k+1]) continue;
968
S2Point const& v0 = intersections[k].second;
969
S2Point const& v1 = intersections[k+1].second;
971
// If the gap from the previous vertex to this one is large
972
// enough, start a new polyline.
973
if (!vertices.empty() &&
974
vertices.back().Angle(v0) > merge_radius.radians()) {
975
out->push_back(new S2Polyline(vertices));
978
// Append this segment to the current polyline, ignoring any
979
// vertices that are too close to the previous vertex.
980
if (vertices.empty()) vertices.push_back(v0);
981
if (vertices.back().Angle(v1) > merge_radius.radians()) {
982
vertices.push_back(v1);
985
intersections.clear();
987
if (!vertices.empty()) {
988
out->push_back(new S2Polyline(vertices));
992
void S2Polygon::IntersectWithPolyline(
994
vector<S2Polyline*> *out) const {
995
IntersectWithPolylineSloppy(a, out, S2EdgeUtil::kIntersectionTolerance);
998
void S2Polygon::IntersectWithPolylineSloppy(
1000
vector<S2Polyline*> *out,
1001
S1Angle vertex_merge_radius) const {
1002
InternalClipPolyline(false, a, out, vertex_merge_radius);
1005
void S2Polygon::SubtractFromPolyline(S2Polyline const* a,
1006
vector<S2Polyline*> *out) const {
1007
SubtractFromPolylineSloppy(a, out, S2EdgeUtil::kIntersectionTolerance);
1010
void S2Polygon::SubtractFromPolylineSloppy(
1011
S2Polyline const* a,
1012
vector<S2Polyline*> *out,
1013
S1Angle vertex_merge_radius) const {
1014
InternalClipPolyline(true, a, out, vertex_merge_radius);
1018
S2Polygon* S2Polygon::DestructiveUnion(vector<S2Polygon*>* polygons) {
1019
return DestructiveUnionSloppy(polygons, S2EdgeUtil::kIntersectionTolerance);
1022
S2Polygon* S2Polygon::DestructiveUnionSloppy(vector<S2Polygon*>* polygons,
1023
S1Angle vertex_merge_radius) {
1024
// Effectively create a priority queue of polygons in order of number of
1025
// vertices. Repeatedly union the two smallest polygons and add the result
1026
// to the queue until we have a single polygon to return.
1027
typedef multimap<int, S2Polygon*> QueueType;
1028
QueueType queue; // Map from # of vertices to polygon.
1029
for (size_t i = 0; i < polygons->size(); ++i)
1030
queue.insert(make_pair((*polygons)[i]->num_vertices(), (*polygons)[i]));
1033
while (queue.size() > 1) {
1034
// Pop two simplest polygons from queue.
1035
QueueType::iterator smallest_it = queue.begin();
1036
int a_size = smallest_it->first;
1037
S2Polygon* a_polygon = smallest_it->second;
1038
queue.erase(smallest_it);
1039
smallest_it = queue.begin();
1040
int b_size = smallest_it->first;
1041
S2Polygon* b_polygon = smallest_it->second;
1042
queue.erase(smallest_it);
1044
// Union and add result back to queue.
1045
S2Polygon* union_polygon = new S2Polygon();
1046
union_polygon->InitToUnionSloppy(a_polygon, b_polygon, vertex_merge_radius);
1049
queue.insert(make_pair(a_size + b_size, union_polygon));
1050
// We assume that the number of vertices in the union polygon is the
1051
// sum of the number of vertices in the original polygons, which is not
1052
// always true, but will almost always be a decent approximation, and
1053
// faster than recomputing.
1057
return new S2Polygon();
1059
return queue.begin()->second;
1062
void S2Polygon::InitToCellUnionBorder(S2CellUnion const& cells) {
1063
// Use a polygon builder to union the cells in the union. Due to rounding
1064
// errors, we can't do an exact union - when a small cell is adjacent to a
1065
// larger cell, the shared edges can fail to line up exactly. Two cell edges
1066
// cannot come closer then kMinWidth, so if we have the polygon builder merge
1067
// edges within half that distance, we should always merge shared edges
1068
// without merging different edges.
1069
S2PolygonBuilderOptions options(S2PolygonBuilderOptions::DIRECTED_XOR());
1070
double min_cell_angle = S2::kMinWidth.GetValue(S2CellId::kMaxLevel);
1071
options.set_vertex_merge_radius(S1Angle::Radians(min_cell_angle / 2));
1072
S2PolygonBuilder builder(options);
1073
for (int i = 0; i < cells.num_cells(); ++i) {
1074
S2Loop cell_loop(S2Cell(cells.cell_id(i)));
1075
builder.AddLoop(&cell_loop);
1077
if (!builder.AssemblePolygon(this, NULL)) {
1078
S2LOG(DFATAL) << "AssemblePolygon failed in InitToCellUnionBorder";
1082
bool S2Polygon::IsNormalized() const {
1083
set<S2Point> vertices;
1084
S2Loop* last_parent = NULL;
1085
for (int i = 0; i < num_loops(); ++i) {
1086
S2Loop* child = loop(i);
1087
if (child->depth() == 0) continue;
1088
S2Loop* parent = loop(GetParent(i));
1089
if (parent != last_parent) {
1091
for (int j = 0; j < parent->num_vertices(); ++j) {
1092
vertices.insert(parent->vertex(j));
1094
last_parent = parent;
1097
for (int j = 0; j < child->num_vertices(); ++j) {
1098
if (vertices.count(child->vertex(j)) > 0) ++count;
1100
if (count > 1) return false;
1105
bool S2Polygon::BoundaryEquals(S2Polygon const* b) const {
1106
if (num_loops() != b->num_loops()) return false;
1108
for (int i = 0; i < num_loops(); ++i) {
1109
S2Loop* a_loop = loop(i);
1110
bool success = false;
1111
for (int j = 0; j < num_loops(); ++j) {
1112
S2Loop* b_loop = b->loop(j);
1113
if ((b_loop->depth() == a_loop->depth()) &&
1114
b_loop->BoundaryEquals(a_loop)) {
1119
if (!success) return false;
1124
bool S2Polygon::BoundaryApproxEquals(S2Polygon const* b,
1125
double max_error) const {
1126
if (num_loops() != b->num_loops()) return false;
1128
// For now, we assume that there is at most one candidate match for each
1129
// loop. (So far this method is just used for testing.)
1131
for (int i = 0; i < num_loops(); ++i) {
1132
S2Loop* a_loop = loop(i);
1133
bool success = false;
1134
for (int j = 0; j < num_loops(); ++j) {
1135
S2Loop* b_loop = b->loop(j);
1136
if (b_loop->depth() == a_loop->depth() &&
1137
b_loop->BoundaryApproxEquals(a_loop, max_error)) {
1142
if (!success) return false;
1147
bool S2Polygon::BoundaryNear(S2Polygon const* b, double max_error) const {
1148
if (num_loops() != b->num_loops()) return false;
1150
// For now, we assume that there is at most one candidate match for each
1151
// loop. (So far this method is just used for testing.)
1153
for (int i = 0; i < num_loops(); ++i) {
1154
S2Loop* a_loop = loop(i);
1155
bool success = false;
1156
for (int j = 0; j < num_loops(); ++j) {
1157
S2Loop* b_loop = b->loop(j);
1158
if (b_loop->depth() == a_loop->depth() &&
1159
b_loop->BoundaryNear(a_loop, max_error)) {
1164
if (!success) return false;
1169
S2Point S2Polygon::Project(S2Point const& point) const {
1170
DCHECK(!loops_.empty());
1172
if (Contains(point)) {
1176
S1Angle min_distance = S1Angle::Radians(10);
1177
int min_loop_index = 0;
1178
int min_vertex_index = 0;
1180
for (int l = 0; l < num_loops(); ++l) {
1181
S2Loop *a_loop = loop(l);
1182
for (int v = 0; v < a_loop->num_vertices(); ++v) {
1183
S1Angle distance_to_segment =
1184
S2EdgeUtil::GetDistance(point,
1186
a_loop->vertex(v + 1));
1187
if (distance_to_segment < min_distance) {
1188
min_distance = distance_to_segment;
1190
min_vertex_index = v;
1195
S2Point closest_point = S2EdgeUtil::GetClosestPoint(
1197
loop(min_loop_index)->vertex(min_vertex_index),
1198
loop(min_loop_index)->vertex(min_vertex_index + 1));
1200
return closest_point;
added
removed
src/third_party/s2/s2polygon.cc
