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