forked from golang/geo
/
polygon.go
1228 lines (1089 loc) · 38.3 KB
/
polygon.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// Copyright 2015 Google Inc. All rights reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
package s2
import (
"fmt"
"io"
"math"
)
// Polygon represents a sequence of zero or more loops; recall that the
// interior of a loop is defined to be its left-hand side (see Loop).
//
// When the polygon is initialized, the given loops are automatically converted
// into a canonical form consisting of "shells" and "holes". Shells and holes
// are both oriented CCW, and are nested hierarchically. The loops are
// reordered to correspond to a pre-order traversal of the nesting hierarchy.
//
// Polygons may represent any region of the sphere with a polygonal boundary,
// including the entire sphere (known as the "full" polygon). The full polygon
// consists of a single full loop (see Loop), whereas the empty polygon has no
// loops at all.
//
// Use FullPolygon() to construct a full polygon. The zero value of Polygon is
// treated as the empty polygon.
//
// Polygons have the following restrictions:
//
// - Loops may not cross, i.e. the boundary of a loop may not intersect
// both the interior and exterior of any other loop.
//
// - Loops may not share edges, i.e. if a loop contains an edge AB, then
// no other loop may contain AB or BA.
//
// - Loops may share vertices, however no vertex may appear twice in a
// single loop (see Loop).
//
// - No loop may be empty. The full loop may appear only in the full polygon.
type Polygon struct {
loops []*Loop
// index is a spatial index of all the polygon loops.
index *ShapeIndex
// hasHoles tracks if this polygon has at least one hole.
hasHoles bool
// numVertices keeps the running total of all of the vertices of the contained loops.
numVertices int
// numEdges tracks the total number of edges in all the loops in this polygon.
numEdges int
// bound is a conservative bound on all points contained by this loop.
// If l.ContainsPoint(P), then l.bound.ContainsPoint(P).
bound Rect
// Since bound is not exact, it is possible that a loop A contains
// another loop B whose bounds are slightly larger. subregionBound
// has been expanded sufficiently to account for this error, i.e.
// if A.Contains(B), then A.subregionBound.Contains(B.bound).
subregionBound Rect
// A slice where element i is the cumulative number of edges in the
// preceding loops in the polygon. This field is used for polygons that
// have a large number of loops, and may be empty for polygons with few loops.
cumulativeEdges []int
}
// PolygonFromLoops constructs a polygon from the given set of loops. The polygon
// interior consists of the points contained by an odd number of loops. (Recall
// that a loop contains the set of points on its left-hand side.)
//
// This method determines the loop nesting hierarchy and assigns every loop a
// depth. Shells have even depths, and holes have odd depths.
//
// Note: The given set of loops are reordered by this method so that the hierarchy
// can be traversed using Parent, LastDescendant and the loops depths.
func PolygonFromLoops(loops []*Loop) *Polygon {
p := &Polygon{}
// Empty polygons do not contain any loops, even the Empty loop.
if len(loops) == 1 && loops[0].IsEmpty() {
p.initLoopProperties()
return p
}
p.loops = loops
p.initNested()
return p
}
// PolygonFromOrientedLoops returns a Polygon from the given set of loops,
// like PolygonFromLoops. It expects loops to be oriented such that the polygon
// interior is on the left-hand side of all loops. This implies that shells
// and holes should have opposite orientations in the input to this method.
// (During initialization, loops representing holes will automatically be
// inverted.)
func PolygonFromOrientedLoops(loops []*Loop) *Polygon {
// Here is the algorithm:
//
// 1. Remember which of the given loops contain OriginPoint.
//
// 2. Invert loops as necessary to ensure that they are nestable (i.e., no
// loop contains the complement of any other loop). This may result in a
// set of loops corresponding to the complement of the given polygon, but
// we will fix that problem later.
//
// We make the loops nestable by first normalizing all the loops (i.e.,
// inverting any loops whose turning angle is negative). This handles
// all loops except those whose turning angle is very close to zero
// (within the maximum error tolerance). Any such loops are inverted if
// and only if they contain OriginPoint(). (In theory this step is only
// necessary if there are at least two such loops.) The resulting set of
// loops is guaranteed to be nestable.
//
// 3. Build the polygon. This yields either the desired polygon or its
// complement.
//
// 4. If there is at least one loop, we find a loop L that is adjacent to
// OriginPoint() (where "adjacent" means that there exists a path
// connecting OriginPoint() to some vertex of L such that the path does
// not cross any loop). There may be a single such adjacent loop, or
// there may be several (in which case they should all have the same
// contains_origin() value). We choose L to be the loop containing the
// origin whose depth is greatest, or loop(0) (a top-level shell) if no
// such loop exists.
//
// 5. If (L originally contained origin) != (polygon contains origin), we
// invert the polygon. This is done by inverting a top-level shell whose
// turning angle is minimal and then fixing the nesting hierarchy. Note
// that because we normalized all the loops initially, this step is only
// necessary if the polygon requires at least one non-normalized loop to
// represent it.
containedOrigin := make(map[*Loop]bool)
for _, l := range loops {
containedOrigin[l] = l.ContainsOrigin()
}
for _, l := range loops {
angle := l.TurningAngle()
if math.Abs(angle) > l.turningAngleMaxError() {
// Normalize the loop.
if angle < 0 {
l.Invert()
}
} else {
// Ensure that the loop does not contain the origin.
if l.ContainsOrigin() {
l.Invert()
}
}
}
p := PolygonFromLoops(loops)
if p.NumLoops() > 0 {
originLoop := p.Loop(0)
polygonContainsOrigin := false
for _, l := range p.Loops() {
if l.ContainsOrigin() {
polygonContainsOrigin = !polygonContainsOrigin
originLoop = l
}
}
if containedOrigin[originLoop] != polygonContainsOrigin {
p.Invert()
}
}
return p
}
// Invert inverts the polygon (replaces it by its complement).
func (p *Polygon) Invert() {
// Inverting any one loop will invert the polygon. The best loop to invert
// is the one whose area is largest, since this yields the smallest area
// after inversion. The loop with the largest area is always at depth 0.
// The descendents of this loop all have their depth reduced by 1, while the
// former siblings of this loop all have their depth increased by 1.
// The empty and full polygons are handled specially.
if p.IsEmpty() {
*p = *FullPolygon()
p.initLoopProperties()
return
}
if p.IsFull() {
*p = Polygon{}
p.initLoopProperties()
return
}
// Find the loop whose area is largest (i.e., whose turning angle is
// smallest), minimizing calls to TurningAngle(). In particular, for
// polygons with a single shell at level 0 there is no need to call
// TurningAngle() at all. (This method is relatively expensive.)
best := 0
const none = 10.0 // Flag that means "not computed yet"
bestAngle := none
for i := 1; i < p.NumLoops(); i++ {
if p.Loop(i).depth != 0 {
continue
}
// We defer computing the turning angle of loop 0 until we discover
// that the polygon has another top-level shell.
if bestAngle == none {
bestAngle = p.Loop(best).TurningAngle()
}
angle := p.Loop(i).TurningAngle()
// We break ties deterministically in order to avoid having the output
// depend on the input order of the loops.
if angle < bestAngle || (angle == bestAngle && compareLoops(p.Loop(i), p.Loop(best)) < 0) {
best = i
bestAngle = angle
}
}
// Build the new loops vector, starting with the inverted loop.
p.Loop(best).Invert()
newLoops := make([]*Loop, 0, p.NumLoops())
// Add the former siblings of this loop as descendants.
lastBest := p.LastDescendant(best)
newLoops = append(newLoops, p.Loop(best))
for i, l := range p.Loops() {
if i < best || i > lastBest {
l.depth++
newLoops = append(newLoops, l)
}
}
// Add the former children of this loop as siblings.
for i, l := range p.Loops() {
if i > best && i <= lastBest {
l.depth--
newLoops = append(newLoops, l)
}
}
p.loops = newLoops
p.initLoopProperties()
}
// Defines a total ordering on Loops that does not depend on the cyclic
// order of loop vertices. This function is used to choose which loop to
// invert in the case where several loops have exactly the same area.
func compareLoops(a, b *Loop) int {
if na, nb := a.NumVertices(), b.NumVertices(); na != nb {
return na - nb
}
ai, aDir := a.CanonicalFirstVertex()
bi, bDir := b.CanonicalFirstVertex()
if aDir != bDir {
return aDir - bDir
}
for n := a.NumVertices() - 1; n >= 0; n, ai, bi = n-1, ai+aDir, bi+bDir {
if cmp := a.Vertex(ai).Cmp(b.Vertex(bi).Vector); cmp != 0 {
return cmp
}
}
return 0
}
// PolygonFromCell returns a Polygon from a single loop created from the given Cell.
func PolygonFromCell(cell Cell) *Polygon {
return PolygonFromLoops([]*Loop{LoopFromCell(cell)})
}
// initNested takes the set of loops in this polygon and performs the nesting
// computations to set the proper nesting and parent/child relationships.
func (p *Polygon) initNested() {
if len(p.loops) == 1 {
p.initOneLoop()
return
}
lm := make(loopMap)
for _, l := range p.loops {
lm.insertLoop(l, nil)
}
// The loops have all been added to the loopMap for ordering. Clear the
// loops slice because we add all the loops in-order in initLoops.
p.loops = nil
// Reorder the loops in depth-first traversal order.
p.initLoops(lm)
p.initLoopProperties()
}
// loopMap is a map of a loop to its immediate children with respect to nesting.
// It is used to determine which loops are shells and which are holes.
type loopMap map[*Loop][]*Loop
// insertLoop adds the given loop to the loop map under the specified parent.
// All children of the new entry are checked to see if the need to move up to
// a different level.
func (lm loopMap) insertLoop(newLoop, parent *Loop) {
var children []*Loop
for done := false; !done; {
children = lm[parent]
done = true
for _, child := range children {
if child.ContainsNested(newLoop) {
parent = child
done = false
break
}
}
}
// Now, we have found a parent for this loop, it may be that some of the
// children of the parent of this loop may now be children of the new loop.
newChildren := lm[newLoop]
for i := 0; i < len(children); {
child := children[i]
if newLoop.ContainsNested(child) {
newChildren = append(newChildren, child)
children = append(children[0:i], children[i+1:]...)
} else {
i++
}
}
lm[newLoop] = newChildren
lm[parent] = append(children, newLoop)
}
// loopStack simplifies access to the loops while being initialized.
type loopStack []*Loop
func (s *loopStack) push(v *Loop) {
*s = append(*s, v)
}
func (s *loopStack) pop() *Loop {
l := len(*s)
r := (*s)[l-1]
*s = (*s)[:l-1]
return r
}
// initLoops walks the mapping of loops to all of their children, and adds them in
// order into to the polygons set of loops.
func (p *Polygon) initLoops(lm loopMap) {
var stack loopStack
stack.push(nil)
depth := -1
for len(stack) > 0 {
loop := stack.pop()
if loop != nil {
depth = loop.depth
p.loops = append(p.loops, loop)
}
children := lm[loop]
for i := len(children) - 1; i >= 0; i-- {
child := children[i]
child.depth = depth + 1
stack.push(child)
}
}
}
// initOneLoop set the properties for a polygon made of a single loop.
// TODO(roberts): Can this be merged with initLoopProperties
func (p *Polygon) initOneLoop() {
p.hasHoles = false
p.numVertices = len(p.loops[0].vertices)
p.bound = p.loops[0].RectBound()
p.subregionBound = ExpandForSubregions(p.bound)
// Ensure the loops depth is set correctly.
p.loops[0].depth = 0
p.initEdgesAndIndex()
}
// initLoopProperties sets the properties for polygons with multiple loops.
func (p *Polygon) initLoopProperties() {
p.numVertices = 0
// the loops depths are set by initNested/initOriented prior to this.
p.bound = EmptyRect()
p.hasHoles = false
for _, l := range p.loops {
if l.IsHole() {
p.hasHoles = true
} else {
p.bound = p.bound.Union(l.RectBound())
}
p.numVertices += l.NumVertices()
}
p.subregionBound = ExpandForSubregions(p.bound)
p.initEdgesAndIndex()
}
// initEdgesAndIndex performs the shape related initializations and adds the final
// polygon to the index.
func (p *Polygon) initEdgesAndIndex() {
p.numEdges = 0
p.cumulativeEdges = nil
if p.IsFull() {
return
}
const maxLinearSearchLoops = 12 // Based on benchmarks.
if len(p.loops) > maxLinearSearchLoops {
p.cumulativeEdges = make([]int, 0, len(p.loops))
}
for _, l := range p.loops {
if p.cumulativeEdges != nil {
p.cumulativeEdges = append(p.cumulativeEdges, p.numEdges)
}
p.numEdges += len(l.vertices)
}
p.index = NewShapeIndex()
p.index.Add(p)
}
// FullPolygon returns a special "full" polygon.
func FullPolygon() *Polygon {
ret := &Polygon{
loops: []*Loop{
FullLoop(),
},
numVertices: len(FullLoop().Vertices()),
bound: FullRect(),
subregionBound: FullRect(),
}
ret.initEdgesAndIndex()
return ret
}
// Validate checks whether this is a valid polygon,
// including checking whether all the loops are themselves valid.
func (p *Polygon) Validate() error {
for i, l := range p.loops {
// Check for loop errors that don't require building a ShapeIndex.
if err := l.findValidationErrorNoIndex(); err != nil {
return fmt.Errorf("loop %d: %v", i, err)
}
// Check that no loop is empty, and that the full loop only appears in the
// full polygon.
if l.IsEmpty() {
return fmt.Errorf("loop %d: empty loops are not allowed", i)
}
if l.IsFull() && len(p.loops) > 1 {
return fmt.Errorf("loop %d: full loop appears in non-full polygon", i)
}
}
// TODO(roberts): Uncomment the remaining checks when they are completed.
// Check for loop self-intersections and loop pairs that cross
// (including duplicate edges and vertices).
// if findSelfIntersection(p.index) {
// return fmt.Errorf("polygon has loop pairs that cross")
// }
// Check whether initOriented detected inconsistent loop orientations.
// if p.hasInconsistentLoopOrientations {
// return fmt.Errorf("inconsistent loop orientations detected")
// }
// Finally, verify the loop nesting hierarchy.
return p.findLoopNestingError()
}
// findLoopNestingError reports if there is an error in the loop nesting hierarchy.
func (p *Polygon) findLoopNestingError() error {
// First check that the loop depths make sense.
lastDepth := -1
for i, l := range p.loops {
depth := l.depth
if depth < 0 || depth > lastDepth+1 {
return fmt.Errorf("loop %d: invalid loop depth (%d)", i, depth)
}
lastDepth = depth
}
// Then check that they correspond to the actual loop nesting. This test
// is quadratic in the number of loops but the cost per iteration is small.
for i, l := range p.loops {
last := p.LastDescendant(i)
for j, l2 := range p.loops {
if i == j {
continue
}
nested := (j >= i+1) && (j <= last)
const reverseB = false
if l.containsNonCrossingBoundary(l2, reverseB) != nested {
nestedStr := ""
if !nested {
nestedStr = "not "
}
return fmt.Errorf("invalid nesting: loop %d should %scontain loop %d", i, nestedStr, j)
}
}
}
return nil
}
// IsEmpty reports whether this is the special "empty" polygon (consisting of no loops).
func (p *Polygon) IsEmpty() bool {
return len(p.loops) == 0
}
// IsFull reports whether this is the special "full" polygon (consisting of a
// single loop that encompasses the entire sphere).
func (p *Polygon) IsFull() bool {
return len(p.loops) == 1 && p.loops[0].IsFull()
}
// NumLoops returns the number of loops in this polygon.
func (p *Polygon) NumLoops() int {
return len(p.loops)
}
// Loops returns the loops in this polygon.
func (p *Polygon) Loops() []*Loop {
return p.loops
}
// Loop returns the loop at the given index. Note that during initialization,
// the given loops are reordered according to a pre-order traversal of the loop
// nesting hierarchy. This implies that every loop is immediately followed by
// its descendants. This hierarchy can be traversed using the methods Parent,
// LastDescendant, and Loop.depth.
func (p *Polygon) Loop(k int) *Loop {
return p.loops[k]
}
// Parent returns the index of the parent of loop k.
// If the loop does not have a parent, ok=false is returned.
func (p *Polygon) Parent(k int) (index int, ok bool) {
// See where we are on the depth hierarchy.
depth := p.loops[k].depth
if depth == 0 {
return -1, false
}
// There may be several loops at the same nesting level as us that share a
// parent loop with us. (Imagine a slice of swiss cheese, of which we are one loop.
// we don't know how many may be next to us before we get back to our parent loop.)
// Move up one position from us, and then begin traversing back through the set of loops
// until we find the one that is our parent or we get to the top of the polygon.
for k--; k >= 0 && p.loops[k].depth <= depth; k-- {
}
return k, true
}
// LastDescendant returns the index of the last loop that is contained within loop k.
// If k is negative, it returns the last loop in the polygon.
// Note that loops are indexed according to a pre-order traversal of the nesting
// hierarchy, so the immediate children of loop k can be found by iterating over
// the loops (k+1)..LastDescendant(k) and selecting those whose depth is equal
// to Loop(k).depth+1.
func (p *Polygon) LastDescendant(k int) int {
if k < 0 {
return len(p.loops) - 1
}
depth := p.loops[k].depth
// Find the next loop immediately past us in the set of loops, and then start
// moving down the list until we either get to the end or find the next loop
// that is higher up the hierarchy than we are.
for k++; k < len(p.loops) && p.loops[k].depth > depth; k++ {
}
return k - 1
}
// CapBound returns a bounding spherical cap.
func (p *Polygon) CapBound() Cap { return p.bound.CapBound() }
// RectBound returns a bounding latitude-longitude rectangle.
func (p *Polygon) RectBound() Rect { return p.bound }
// ContainsPoint reports whether the polygon contains the point.
func (p *Polygon) ContainsPoint(point Point) bool {
// NOTE: A bounds check slows down this function by about 50%. It is
// worthwhile only when it might allow us to delay building the index.
if !p.index.IsFresh() && !p.bound.ContainsPoint(point) {
return false
}
// For small polygons, and during initial construction, it is faster to just
// check all the crossing.
const maxBruteForceVertices = 32
if p.numVertices < maxBruteForceVertices || p.index == nil {
inside := false
for _, l := range p.loops {
// use loops bruteforce to avoid building the index on each loop.
inside = inside != l.bruteForceContainsPoint(point)
}
return inside
}
// Otherwise we look up the ShapeIndex cell containing this point.
return NewContainsPointQuery(p.index, VertexModelSemiOpen).Contains(point)
}
// ContainsCell reports whether the polygon contains the given cell.
func (p *Polygon) ContainsCell(cell Cell) bool {
it := p.index.Iterator()
relation := it.LocateCellID(cell.ID())
// If "cell" is disjoint from all index cells, it is not contained.
// Similarly, if "cell" is subdivided into one or more index cells then it
// is not contained, since index cells are subdivided only if they (nearly)
// intersect a sufficient number of edges. (But note that if "cell" itself
// is an index cell then it may be contained, since it could be a cell with
// no edges in the loop interior.)
if relation != Indexed {
return false
}
// Otherwise check if any edges intersect "cell".
if p.boundaryApproxIntersects(it, cell) {
return false
}
// Otherwise check if the loop contains the center of "cell".
return p.iteratorContainsPoint(it, cell.Center())
}
// IntersectsCell reports whether the polygon intersects the given cell.
func (p *Polygon) IntersectsCell(cell Cell) bool {
it := p.index.Iterator()
relation := it.LocateCellID(cell.ID())
// If cell does not overlap any index cell, there is no intersection.
if relation == Disjoint {
return false
}
// If cell is subdivided into one or more index cells, there is an
// intersection to within the S2ShapeIndex error bound (see Contains).
if relation == Subdivided {
return true
}
// If cell is an index cell, there is an intersection because index cells
// are created only if they have at least one edge or they are entirely
// contained by the loop.
if it.CellID() == cell.id {
return true
}
// Otherwise check if any edges intersect cell.
if p.boundaryApproxIntersects(it, cell) {
return true
}
// Otherwise check if the loop contains the center of cell.
return p.iteratorContainsPoint(it, cell.Center())
}
// CellUnionBound computes a covering of the Polygon.
func (p *Polygon) CellUnionBound() []CellID {
// TODO(roberts): Use ShapeIndexRegion when it's available.
return p.CapBound().CellUnionBound()
}
// boundaryApproxIntersects reports whether the loop's boundary intersects cell.
// It may also return true when the loop boundary does not intersect cell but
// some edge comes within the worst-case error tolerance.
//
// This requires that it.Locate(cell) returned Indexed.
func (p *Polygon) boundaryApproxIntersects(it *ShapeIndexIterator, cell Cell) bool {
aClipped := it.IndexCell().findByShapeID(0)
// If there are no edges, there is no intersection.
if len(aClipped.edges) == 0 {
return false
}
// We can save some work if cell is the index cell itself.
if it.CellID() == cell.ID() {
return true
}
// Otherwise check whether any of the edges intersect cell.
maxError := (faceClipErrorUVCoord + intersectsRectErrorUVDist)
bound := cell.BoundUV().ExpandedByMargin(maxError)
for _, e := range aClipped.edges {
edge := p.index.Shape(0).Edge(e)
v0, v1, ok := ClipToPaddedFace(edge.V0, edge.V1, cell.Face(), maxError)
if ok && edgeIntersectsRect(v0, v1, bound) {
return true
}
}
return false
}
// iteratorContainsPoint reports whether the iterator that is positioned at the
// ShapeIndexCell that may contain p, contains the point p.
func (p *Polygon) iteratorContainsPoint(it *ShapeIndexIterator, point Point) bool {
// Test containment by drawing a line segment from the cell center to the
// given point and counting edge crossings.
aClipped := it.IndexCell().findByShapeID(0)
inside := aClipped.containsCenter
if len(aClipped.edges) == 0 {
return inside
}
// This block requires ShapeIndex.
crosser := NewEdgeCrosser(it.Center(), point)
shape := p.index.Shape(0)
for _, e := range aClipped.edges {
edge := shape.Edge(e)
inside = inside != crosser.EdgeOrVertexCrossing(edge.V0, edge.V1)
}
return inside
}
// Shape Interface
// NumEdges returns the number of edges in this shape.
func (p *Polygon) NumEdges() int {
return p.numEdges
}
// Edge returns endpoints for the given edge index.
func (p *Polygon) Edge(e int) Edge {
var i int
if len(p.cumulativeEdges) > 0 {
for i = range p.cumulativeEdges {
if i+1 >= len(p.cumulativeEdges) || e < p.cumulativeEdges[i+1] {
e -= p.cumulativeEdges[i]
break
}
}
} else {
// When the number of loops is small, use linear search. Most often
// there is exactly one loop and the code below executes zero times.
for i = 0; e >= len(p.Loop(i).vertices); i++ {
e -= len(p.Loop(i).vertices)
}
}
return Edge{p.Loop(i).OrientedVertex(e), p.Loop(i).OrientedVertex(e + 1)}
}
// ReferencePoint returns the reference point for this polygon.
func (p *Polygon) ReferencePoint() ReferencePoint {
containsOrigin := false
for _, l := range p.loops {
containsOrigin = containsOrigin != l.ContainsOrigin()
}
return OriginReferencePoint(containsOrigin)
}
// NumChains reports the number of contiguous edge chains in the Polygon.
func (p *Polygon) NumChains() int {
return p.NumLoops()
}
// Chain returns the i-th edge Chain (loop) in the Shape.
func (p *Polygon) Chain(chainID int) Chain {
if p.cumulativeEdges != nil {
return Chain{p.cumulativeEdges[chainID], len(p.Loop(chainID).vertices)}
}
e := 0
for j := 0; j < chainID; j++ {
e += len(p.Loop(j).vertices)
}
// Polygon represents a full loop as a loop with one vertex, while
// Shape represents a full loop as a chain with no vertices.
if numVertices := p.Loop(chainID).NumVertices(); numVertices != 1 {
return Chain{e, numVertices}
}
return Chain{e, 0}
}
// ChainEdge returns the j-th edge of the i-th edge Chain (loop).
func (p *Polygon) ChainEdge(i, j int) Edge {
return Edge{p.Loop(i).OrientedVertex(j), p.Loop(i).OrientedVertex(j + 1)}
}
// ChainPosition returns a pair (i, j) such that edgeID is the j-th edge
// of the i-th edge Chain.
func (p *Polygon) ChainPosition(edgeID int) ChainPosition {
var i int
if len(p.cumulativeEdges) > 0 {
for i = range p.cumulativeEdges {
if i+1 >= len(p.cumulativeEdges) || edgeID < p.cumulativeEdges[i+1] {
edgeID -= p.cumulativeEdges[i]
break
}
}
} else {
// When the number of loops is small, use linear search. Most often
// there is exactly one loop and the code below executes zero times.
for i = 0; edgeID >= len(p.Loop(i).vertices); i++ {
edgeID -= len(p.Loop(i).vertices)
}
}
// TODO(roberts): unify this and Edge since they are mostly identical.
return ChainPosition{i, edgeID}
}
// Dimension returns the dimension of the geometry represented by this Polygon.
func (p *Polygon) Dimension() int { return 2 }
func (p *Polygon) typeTag() typeTag { return typeTagPolygon }
func (p *Polygon) privateInterface() {}
// Contains reports whether this polygon contains the other polygon.
// Specifically, it reports whether all the points in the other polygon
// are also in this polygon.
func (p *Polygon) Contains(o *Polygon) bool {
// If both polygons have one loop, use the more efficient Loop method.
// Note that Loop's Contains does its own bounding rectangle check.
if len(p.loops) == 1 && len(o.loops) == 1 {
return p.loops[0].Contains(o.loops[0])
}
// Otherwise if neither polygon has holes, we can still use the more
// efficient Loop's Contains method (rather than compareBoundary),
// but it's worthwhile to do our own bounds check first.
if !p.subregionBound.Contains(o.bound) {
// Even though Bound(A) does not contain Bound(B), it is still possible
// that A contains B. This can only happen when union of the two bounds
// spans all longitudes. For example, suppose that B consists of two
// shells with a longitude gap between them, while A consists of one shell
// that surrounds both shells of B but goes the other way around the
// sphere (so that it does not intersect the longitude gap).
if !p.bound.Lng.Union(o.bound.Lng).IsFull() {
return false
}
}
if !p.hasHoles && !o.hasHoles {
for _, l := range o.loops {
if !p.anyLoopContains(l) {
return false
}
}
return true
}
// Polygon A contains B iff B does not intersect the complement of A. From
// the intersection algorithm below, this means that the complement of A
// must exclude the entire boundary of B, and B must exclude all shell
// boundaries of the complement of A. (It can be shown that B must then
// exclude the entire boundary of the complement of A.) The first call
// below returns false if the boundaries cross, therefore the second call
// does not need to check for any crossing edges (which makes it cheaper).
return p.containsBoundary(o) && o.excludesNonCrossingComplementShells(p)
}
// Intersects reports whether this polygon intersects the other polygon, i.e.
// if there is a point that is contained by both polygons.
func (p *Polygon) Intersects(o *Polygon) bool {
// If both polygons have one loop, use the more efficient Loop method.
// Note that Loop Intersects does its own bounding rectangle check.
if len(p.loops) == 1 && len(o.loops) == 1 {
return p.loops[0].Intersects(o.loops[0])
}
// Otherwise if neither polygon has holes, we can still use the more
// efficient Loop.Intersects method. The polygons intersect if and
// only if some pair of loop regions intersect.
if !p.bound.Intersects(o.bound) {
return false
}
if !p.hasHoles && !o.hasHoles {
for _, l := range o.loops {
if p.anyLoopIntersects(l) {
return true
}
}
return false
}
// Polygon A is disjoint from B if A excludes the entire boundary of B and B
// excludes all shell boundaries of A. (It can be shown that B must then
// exclude the entire boundary of A.) The first call below returns false if
// the boundaries cross, therefore the second call does not need to check
// for crossing edges.
return !p.excludesBoundary(o) || !o.excludesNonCrossingShells(p)
}
// compareBoundary returns +1 if this polygon contains the boundary of B, -1 if A
// excludes the boundary of B, and 0 if the boundaries of A and B cross.
func (p *Polygon) compareBoundary(o *Loop) int {
result := -1
for i := 0; i < len(p.loops) && result != 0; i++ {
// If B crosses any loop of A, the result is 0. Otherwise the result
// changes sign each time B is contained by a loop of A.
result *= -p.loops[i].compareBoundary(o)
}
return result
}
// containsBoundary reports whether this polygon contains the entire boundary of B.
func (p *Polygon) containsBoundary(o *Polygon) bool {
for _, l := range o.loops {
if p.compareBoundary(l) <= 0 {
return false
}
}
return true
}
// excludesBoundary reports whether this polygon excludes the entire boundary of B.
func (p *Polygon) excludesBoundary(o *Polygon) bool {
for _, l := range o.loops {
if p.compareBoundary(l) >= 0 {
return false
}
}
return true
}
// containsNonCrossingBoundary reports whether polygon A contains the boundary of
// loop B. Shared edges are handled according to the rule described in loops
// containsNonCrossingBoundary.
func (p *Polygon) containsNonCrossingBoundary(o *Loop, reverse bool) bool {
var inside bool
for _, l := range p.loops {
x := l.containsNonCrossingBoundary(o, reverse)
inside = (inside != x)
}
return inside
}
// excludesNonCrossingShells reports wheterh given two polygons A and B such that the
// boundary of A does not cross any loop of B, if A excludes all shell boundaries of B.
func (p *Polygon) excludesNonCrossingShells(o *Polygon) bool {
for _, l := range o.loops {
if l.IsHole() {
continue
}
if p.containsNonCrossingBoundary(l, false) {
return false
}
}
return true
}
// excludesNonCrossingComplementShells reports whether given two polygons A and B
// such that the boundary of A does not cross any loop of B, if A excludes all
// shell boundaries of the complement of B.
func (p *Polygon) excludesNonCrossingComplementShells(o *Polygon) bool {
// Special case to handle the complement of the empty or full polygons.
if o.IsEmpty() {
return !p.IsFull()
}
if o.IsFull() {
return true
}
// Otherwise the complement of B may be obtained by inverting loop(0) and
// then swapping the shell/hole status of all other loops. This implies
// that the shells of the complement consist of loop 0 plus all the holes of
// the original polygon.
for j, l := range o.loops {
if j > 0 && !l.IsHole() {
continue
}
// The interior of the complement is to the right of loop 0, and to the
// left of the loops that were originally holes.
if p.containsNonCrossingBoundary(l, j == 0) {
return false
}
}
return true
}
// anyLoopContains reports whether any loop in this polygon contains the given loop.
func (p *Polygon) anyLoopContains(o *Loop) bool {
for _, l := range p.loops {
if l.Contains(o) {
return true
}
}
return false
}
// anyLoopIntersects reports whether any loop in this polygon intersects the given loop.
func (p *Polygon) anyLoopIntersects(o *Loop) bool {
for _, l := range p.loops {
if l.Intersects(o) {