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 "github.com/golang/geo/r1" 23 "github.com/golang/geo/r3" 24 "github.com/golang/geo/s1" 25 ) 26 27 // Loop represents a simple spherical polygon. It consists of a sequence 28 // of vertices where the first vertex is implicitly connected to the 29 // last. All loops are defined to have a CCW orientation, i.e. the interior of 30 // the loop is on the left side of the edges. This implies that a clockwise 31 // loop enclosing a small area is interpreted to be a CCW loop enclosing a 32 // very large area. 33 // 34 // Loops are not allowed to have any duplicate vertices (whether adjacent or 35 // not). Non-adjacent edges are not allowed to intersect, and furthermore edges 36 // of length 180 degrees are not allowed (i.e., adjacent vertices cannot be 37 // antipodal). Loops must have at least 3 vertices (except for the "empty" and 38 // "full" loops discussed below). 39 // 40 // There are two special loops: the "empty" loop contains no points and the 41 // "full" loop contains all points. These loops do not have any edges, but to 42 // preserve the invariant that every loop can be represented as a vertex 43 // chain, they are defined as having exactly one vertex each (see EmptyLoop 44 // and FullLoop). 45 type Loop struct { 46 vertices []Point 47 48 // originInside keeps a precomputed value whether this loop contains the origin 49 // versus computing from the set of vertices every time. 50 originInside bool 51 52 // depth is the nesting depth of this Loop if it is contained by a Polygon 53 // or other shape and is used to determine if this loop represents a hole 54 // or a filled in portion. 55 depth int 56 57 // bound is a conservative bound on all points contained by this loop. 58 // If l.ContainsPoint(P), then l.bound.ContainsPoint(P). 59 bound Rect 60 61 // Since bound is not exact, it is possible that a loop A contains 62 // another loop B whose bounds are slightly larger. subregionBound 63 // has been expanded sufficiently to account for this error, i.e. 64 // if A.Contains(B), then A.subregionBound.Contains(B.bound). 65 subregionBound Rect 66 67 // index is the spatial index for this Loop. 68 index *ShapeIndex 69 } 70 71 // LoopFromPoints constructs a loop from the given points. 72 func LoopFromPoints(pts []Point) *Loop { 73 l := &Loop{ 74 vertices: pts, 75 index: NewShapeIndex(), 76 } 77 78 l.initOriginAndBound() 79 return l 80 } 81 82 // LoopFromCell constructs a loop corresponding to the given cell. 83 // 84 // Note that the loop and cell *do not* contain exactly the same set of 85 // points, because Loop and Cell have slightly different definitions of 86 // point containment. For example, a Cell vertex is contained by all 87 // four neighboring Cells, but it is contained by exactly one of four 88 // Loops constructed from those cells. As another example, the cell 89 // coverings of cell and LoopFromCell(cell) will be different, because the 90 // loop contains points on its boundary that actually belong to other cells 91 // (i.e., the covering will include a layer of neighboring cells). 92 func LoopFromCell(c Cell) *Loop { 93 l := &Loop{ 94 vertices: []Point{ 95 c.Vertex(0), 96 c.Vertex(1), 97 c.Vertex(2), 98 c.Vertex(3), 99 }, 100 index: NewShapeIndex(), 101 } 102 103 l.initOriginAndBound() 104 return l 105 } 106 107 // These two points are used for the special Empty and Full loops. 108 var ( 109 emptyLoopPoint = Point{r3.Vector{X: 0, Y: 0, Z: 1}} 110 fullLoopPoint = Point{r3.Vector{X: 0, Y: 0, Z: -1}} 111 ) 112 113 // EmptyLoop returns a special "empty" loop. 114 func EmptyLoop() *Loop { 115 return LoopFromPoints([]Point{emptyLoopPoint}) 116 } 117 118 // FullLoop returns a special "full" loop. 119 func FullLoop() *Loop { 120 return LoopFromPoints([]Point{fullLoopPoint}) 121 } 122 123 // initOriginAndBound sets the origin containment for the given point and then calls 124 // the initialization for the bounds objects and the internal index. 125 func (l *Loop) initOriginAndBound() { 126 if len(l.vertices) < 3 { 127 // Check for the special "empty" and "full" loops (which have one vertex). 128 if !l.isEmptyOrFull() { 129 l.originInside = false 130 return 131 } 132 133 // This is the special empty or full loop, so the origin depends on if 134 // the vertex is in the southern hemisphere or not. 135 l.originInside = l.vertices[0].Z < 0 136 } else { 137 // The brute force point containment algorithm works by counting edge 138 // crossings starting at a fixed reference point (chosen as OriginPoint() 139 // for historical reasons). Loop initialization would be more efficient 140 // if we used a loop vertex such as vertex(0) as the reference point 141 // instead, however making this change would be a lot of work because 142 // originInside is currently part of the Encode() format. 143 // 144 // In any case, we initialize originInside by first guessing that it is 145 // outside, and then seeing whether we get the correct containment result 146 // for vertex 1. If the result is incorrect, the origin must be inside 147 // the loop instead. Note that the Loop is not necessarily valid and so 148 // we need to check the requirements of AngleContainsVertex first. 149 v1Inside := l.vertices[0] != l.vertices[1] && 150 l.vertices[2] != l.vertices[1] && 151 AngleContainsVertex(l.vertices[0], l.vertices[1], l.vertices[2]) 152 153 // initialize before calling ContainsPoint 154 l.originInside = false 155 156 // Note that ContainsPoint only does a bounds check once initIndex 157 // has been called, so it doesn't matter that bound is undefined here. 158 if v1Inside != l.ContainsPoint(l.vertices[1]) { 159 l.originInside = true 160 } 161 162 } 163 164 // We *must* call initBound before initializing the index, because 165 // initBound calls ContainsPoint which does a bounds check before using 166 // the index. 167 l.initBound() 168 169 // Create a new index and add us to it. 170 l.index = NewShapeIndex() 171 l.index.Add(l) 172 } 173 174 // initBound sets up the approximate bounding Rects for this loop. 175 func (l *Loop) initBound() { 176 if len(l.vertices) == 0 { 177 *l = *EmptyLoop() 178 return 179 } 180 // Check for the special "empty" and "full" loops. 181 if l.isEmptyOrFull() { 182 if l.IsEmpty() { 183 l.bound = EmptyRect() 184 } else { 185 l.bound = FullRect() 186 } 187 l.subregionBound = l.bound 188 return 189 } 190 191 // The bounding rectangle of a loop is not necessarily the same as the 192 // bounding rectangle of its vertices. First, the maximal latitude may be 193 // attained along the interior of an edge. Second, the loop may wrap 194 // entirely around the sphere (e.g. a loop that defines two revolutions of a 195 // candy-cane stripe). Third, the loop may include one or both poles. 196 // Note that a small clockwise loop near the equator contains both poles. 197 bounder := NewRectBounder() 198 for i := 0; i <= len(l.vertices); i++ { // add vertex 0 twice 199 bounder.AddPoint(l.Vertex(i)) 200 } 201 b := bounder.RectBound() 202 203 if l.ContainsPoint(Point{r3.Vector{0, 0, 1}}) { 204 b = Rect{r1.Interval{b.Lat.Lo, math.Pi / 2}, s1.FullInterval()} 205 } 206 // If a loop contains the south pole, then either it wraps entirely 207 // around the sphere (full longitude range), or it also contains the 208 // north pole in which case b.Lng.IsFull() due to the test above. 209 // Either way, we only need to do the south pole containment test if 210 // b.Lng.IsFull(). 211 if b.Lng.IsFull() && l.ContainsPoint(Point{r3.Vector{0, 0, -1}}) { 212 b.Lat.Lo = -math.Pi / 2 213 } 214 l.bound = b 215 l.subregionBound = ExpandForSubregions(l.bound) 216 } 217 218 // Validate checks whether this is a valid loop. 219 func (l *Loop) Validate() error { 220 if err := l.findValidationErrorNoIndex(); err != nil { 221 return err 222 } 223 224 // Check for intersections between non-adjacent edges (including at vertices) 225 // TODO(roberts): Once shapeutil gets findAnyCrossing uncomment this. 226 // return findAnyCrossing(l.index) 227 228 return nil 229 } 230 231 // findValidationErrorNoIndex reports whether this is not a valid loop, but 232 // skips checks that would require a ShapeIndex to be built for the loop. This 233 // is primarily used by Polygon to do validation so it doesn't trigger the 234 // creation of unneeded ShapeIndices. 235 func (l *Loop) findValidationErrorNoIndex() error { 236 // All vertices must be unit length. 237 for i, v := range l.vertices { 238 if !v.IsUnit() { 239 return fmt.Errorf("vertex %d is not unit length", i) 240 } 241 } 242 243 // Loops must have at least 3 vertices (except for empty and full). 244 if len(l.vertices) < 3 { 245 if l.isEmptyOrFull() { 246 return nil // Skip remaining tests. 247 } 248 return fmt.Errorf("non-empty, non-full loops must have at least 3 vertices") 249 } 250 251 // Loops are not allowed to have any duplicate vertices or edge crossings. 252 // We split this check into two parts. First we check that no edge is 253 // degenerate (identical endpoints). Then we check that there are no 254 // intersections between non-adjacent edges (including at vertices). The 255 // second check needs the ShapeIndex, so it does not fall within the scope 256 // of this method. 257 for i, v := range l.vertices { 258 if v == l.Vertex(i+1) { 259 return fmt.Errorf("edge %d is degenerate (duplicate vertex)", i) 260 } 261 262 // Antipodal vertices are not allowed. 263 if other := (Point{l.Vertex(i + 1).Mul(-1)}); v == other { 264 return fmt.Errorf("vertices %d and %d are antipodal", i, 265 (i+1)%len(l.vertices)) 266 } 267 } 268 269 return nil 270 } 271 272 // Contains reports whether the region contained by this loop is a superset of the 273 // region contained by the given other loop. 274 func (l *Loop) Contains(o *Loop) bool { 275 // For a loop A to contain the loop B, all of the following must 276 // be true: 277 // 278 // (1) There are no edge crossings between A and B except at vertices. 279 // 280 // (2) At every vertex that is shared between A and B, the local edge 281 // ordering implies that A contains B. 282 // 283 // (3) If there are no shared vertices, then A must contain a vertex of B 284 // and B must not contain a vertex of A. (An arbitrary vertex may be 285 // chosen in each case.) 286 // 287 // The second part of (3) is necessary to detect the case of two loops whose 288 // union is the entire sphere, i.e. two loops that contains each other's 289 // boundaries but not each other's interiors. 290 if !l.subregionBound.Contains(o.bound) { 291 return false 292 } 293 294 // Special cases to handle either loop being empty or full. 295 if l.isEmptyOrFull() || o.isEmptyOrFull() { 296 return l.IsFull() || o.IsEmpty() 297 } 298 299 // Check whether there are any edge crossings, and also check the loop 300 // relationship at any shared vertices. 301 relation := &containsRelation{} 302 if hasCrossingRelation(l, o, relation) { 303 return false 304 } 305 306 // There are no crossings, and if there are any shared vertices then A 307 // contains B locally at each shared vertex. 308 if relation.foundSharedVertex { 309 return true 310 } 311 312 // Since there are no edge intersections or shared vertices, we just need to 313 // test condition (3) above. We can skip this test if we discovered that A 314 // contains at least one point of B while checking for edge crossings. 315 if !l.ContainsPoint(o.Vertex(0)) { 316 return false 317 } 318 319 // We still need to check whether (A union B) is the entire sphere. 320 // Normally this check is very cheap due to the bounding box precondition. 321 if (o.subregionBound.Contains(l.bound) || o.bound.Union(l.bound).IsFull()) && 322 o.ContainsPoint(l.Vertex(0)) { 323 return false 324 } 325 return true 326 } 327 328 // Intersects reports whether the region contained by this loop intersects the region 329 // contained by the other loop. 330 func (l *Loop) Intersects(o *Loop) bool { 331 // Given two loops, A and B, A.Intersects(B) if and only if !A.Complement().Contains(B). 332 // 333 // This code is similar to Contains, but is optimized for the case 334 // where both loops enclose less than half of the sphere. 335 if !l.bound.Intersects(o.bound) { 336 return false 337 } 338 339 // Check whether there are any edge crossings, and also check the loop 340 // relationship at any shared vertices. 341 relation := &intersectsRelation{} 342 if hasCrossingRelation(l, o, relation) { 343 return true 344 } 345 if relation.foundSharedVertex { 346 return false 347 } 348 349 // Since there are no edge intersections or shared vertices, the loops 350 // intersect only if A contains B, B contains A, or the two loops contain 351 // each other's boundaries. These checks are usually cheap because of the 352 // bounding box preconditions. Note that neither loop is empty (because of 353 // the bounding box check above), so it is safe to access vertex(0). 354 355 // Check whether A contains B, or A and B contain each other's boundaries. 356 // (Note that A contains all the vertices of B in either case.) 357 if l.subregionBound.Contains(o.bound) || l.bound.Union(o.bound).IsFull() { 358 if l.ContainsPoint(o.Vertex(0)) { 359 return true 360 } 361 } 362 // Check whether B contains A. 363 if o.subregionBound.Contains(l.bound) { 364 if o.ContainsPoint(l.Vertex(0)) { 365 return true 366 } 367 } 368 return false 369 } 370 371 // Equal reports whether two loops have the same vertices in the same linear order 372 // (i.e., cyclic rotations are not allowed). 373 func (l *Loop) Equal(other *Loop) bool { 374 if len(l.vertices) != len(other.vertices) { 375 return false 376 } 377 378 for i, v := range l.vertices { 379 if v != other.Vertex(i) { 380 return false 381 } 382 } 383 return true 384 } 385 386 // BoundaryEqual reports whether the two loops have the same boundary. This is 387 // true if and only if the loops have the same vertices in the same cyclic order 388 // (i.e., the vertices may be cyclically rotated). The empty and full loops are 389 // considered to have different boundaries. 390 func (l *Loop) BoundaryEqual(o *Loop) bool { 391 if len(l.vertices) != len(o.vertices) { 392 return false 393 } 394 395 // Special case to handle empty or full loops. Since they have the same 396 // number of vertices, if one loop is empty/full then so is the other. 397 if l.isEmptyOrFull() { 398 return l.IsEmpty() == o.IsEmpty() 399 } 400 401 // Loop through the vertices to find the first of ours that matches the 402 // starting vertex of the other loop. Use that offset to then 'align' the 403 // vertices for comparison. 404 for offset, vertex := range l.vertices { 405 if vertex == o.Vertex(0) { 406 // There is at most one starting offset since loop vertices are unique. 407 for i := 0; i < len(l.vertices); i++ { 408 if l.Vertex(i+offset) != o.Vertex(i) { 409 return false 410 } 411 } 412 return true 413 } 414 } 415 return false 416 } 417 418 // compareBoundary returns +1 if this loop contains the boundary of the other loop, 419 // -1 if it excludes the boundary of the other, and 0 if the boundaries of the two 420 // loops cross. Shared edges are handled as follows: 421 // 422 // If XY is a shared edge, define Reversed(XY) to be true if XY 423 // appears in opposite directions in both loops. 424 // Then this loop contains XY if and only if Reversed(XY) == the other loop is a hole. 425 // (Intuitively, this checks whether this loop contains a vanishingly small region 426 // extending from the boundary of the other toward the interior of the polygon to 427 // which the other belongs.) 428 // 429 // This function is used for testing containment and intersection of 430 // multi-loop polygons. Note that this method is not symmetric, since the 431 // result depends on the direction of this loop but not on the direction of 432 // the other loop (in the absence of shared edges). 433 // 434 // This requires that neither loop is empty, and if other loop IsFull, then it must not 435 // be a hole. 436 func (l *Loop) compareBoundary(o *Loop) int { 437 // The bounds must intersect for containment or crossing. 438 if !l.bound.Intersects(o.bound) { 439 return -1 440 } 441 442 // Full loops are handled as though the loop surrounded the entire sphere. 443 if l.IsFull() { 444 return 1 445 } 446 if o.IsFull() { 447 return -1 448 } 449 450 // Check whether there are any edge crossings, and also check the loop 451 // relationship at any shared vertices. 452 relation := newCompareBoundaryRelation(o.IsHole()) 453 if hasCrossingRelation(l, o, relation) { 454 return 0 455 } 456 if relation.foundSharedVertex { 457 if relation.containsEdge { 458 return 1 459 } 460 return -1 461 } 462 463 // There are no edge intersections or shared vertices, so we can check 464 // whether A contains an arbitrary vertex of B. 465 if l.ContainsPoint(o.Vertex(0)) { 466 return 1 467 } 468 return -1 469 } 470 471 // ContainsOrigin reports true if this loop contains s2.OriginPoint(). 472 func (l *Loop) ContainsOrigin() bool { 473 return l.originInside 474 } 475 476 // ReferencePoint returns the reference point for this loop. 477 func (l *Loop) ReferencePoint() ReferencePoint { 478 return OriginReferencePoint(l.originInside) 479 } 480 481 // NumEdges returns the number of edges in this shape. 482 func (l *Loop) NumEdges() int { 483 if l.isEmptyOrFull() { 484 return 0 485 } 486 return len(l.vertices) 487 } 488 489 // Edge returns the endpoints for the given edge index. 490 func (l *Loop) Edge(i int) Edge { 491 return Edge{l.Vertex(i), l.Vertex(i + 1)} 492 } 493 494 // NumChains reports the number of contiguous edge chains in the Loop. 495 func (l *Loop) NumChains() int { 496 if l.IsEmpty() { 497 return 0 498 } 499 return 1 500 } 501 502 // Chain returns the i-th edge chain in the Shape. 503 func (l *Loop) Chain(chainID int) Chain { 504 return Chain{0, l.NumEdges()} 505 } 506 507 // ChainEdge returns the j-th edge of the i-th edge chain. 508 func (l *Loop) ChainEdge(chainID, offset int) Edge { 509 return Edge{l.Vertex(offset), l.Vertex(offset + 1)} 510 } 511 512 // ChainPosition returns a ChainPosition pair (i, j) such that edgeID is the 513 // j-th edge of the Loop. 514 func (l *Loop) ChainPosition(edgeID int) ChainPosition { 515 return ChainPosition{0, edgeID} 516 } 517 518 // Dimension returns the dimension of the geometry represented by this Loop. 519 func (l *Loop) Dimension() int { return 2 } 520 521 func (l *Loop) typeTag() typeTag { return typeTagNone } 522 523 func (l *Loop) privateInterface() {} 524 525 // IsEmpty reports true if this is the special empty loop that contains no points. 526 func (l *Loop) IsEmpty() bool { 527 return l.isEmptyOrFull() && !l.ContainsOrigin() 528 } 529 530 // IsFull reports true if this is the special full loop that contains all points. 531 func (l *Loop) IsFull() bool { 532 return l.isEmptyOrFull() && l.ContainsOrigin() 533 } 534 535 // isEmptyOrFull reports true if this loop is either the "empty" or "full" special loops. 536 func (l *Loop) isEmptyOrFull() bool { 537 return len(l.vertices) == 1 538 } 539 540 // Vertices returns the vertices in the loop. 541 func (l *Loop) Vertices() []Point { 542 return l.vertices 543 } 544 545 // RectBound returns a tight bounding rectangle. If the loop contains the point, 546 // the bound also contains it. 547 func (l *Loop) RectBound() Rect { 548 return l.bound 549 } 550 551 // CapBound returns a bounding cap that may have more padding than the corresponding 552 // RectBound. The bound is conservative such that if the loop contains a point P, 553 // the bound also contains it. 554 func (l *Loop) CapBound() Cap { 555 return l.bound.CapBound() 556 } 557 558 // Vertex returns the vertex for the given index. For convenience, the vertex indices 559 // wrap automatically for methods that do index math such as Edge. 560 // i.e., Vertex(NumEdges() + n) is the same as Vertex(n). 561 func (l *Loop) Vertex(i int) Point { 562 return l.vertices[i%len(l.vertices)] 563 } 564 565 // OrientedVertex returns the vertex in reverse order if the loop represents a polygon 566 // hole. For example, arguments 0, 1, 2 are mapped to vertices n-1, n-2, n-3, where 567 // n == len(vertices). This ensures that the interior of the polygon is always to 568 // the left of the vertex chain. 569 // 570 // This requires: 0 <= i < 2 * len(vertices) 571 func (l *Loop) OrientedVertex(i int) Point { 572 j := i - len(l.vertices) 573 if j < 0 { 574 j = i 575 } 576 if l.IsHole() { 577 j = len(l.vertices) - 1 - j 578 } 579 return l.Vertex(j) 580 } 581 582 // NumVertices returns the number of vertices in this loop. 583 func (l *Loop) NumVertices() int { 584 return len(l.vertices) 585 } 586 587 // bruteForceContainsPoint reports if the given point is contained by this loop. 588 // This method does not use the ShapeIndex, so it is only preferable below a certain 589 // size of loop. 590 func (l *Loop) bruteForceContainsPoint(p Point) bool { 591 origin := OriginPoint() 592 inside := l.originInside 593 crosser := NewChainEdgeCrosser(origin, p, l.Vertex(0)) 594 for i := 1; i <= len(l.vertices); i++ { // add vertex 0 twice 595 inside = inside != crosser.EdgeOrVertexChainCrossing(l.Vertex(i)) 596 } 597 return inside 598 } 599 600 // ContainsPoint returns true if the loop contains the point. 601 func (l *Loop) ContainsPoint(p Point) bool { 602 if !l.index.IsFresh() && !l.bound.ContainsPoint(p) { 603 return false 604 } 605 606 // For small loops it is faster to just check all the crossings. We also 607 // use this method during loop initialization because InitOriginAndBound() 608 // calls Contains() before InitIndex(). Otherwise, we keep track of the 609 // number of calls to Contains() and only build the index when enough calls 610 // have been made so that we think it is worth the effort. Note that the 611 // code below is structured so that if many calls are made in parallel only 612 // one thread builds the index, while the rest continue using brute force 613 // until the index is actually available. 614 615 const maxBruteForceVertices = 32 616 // TODO(roberts): add unindexed contains calls tracking 617 618 if len(l.index.shapes) == 0 || // Index has not been initialized yet. 619 len(l.vertices) <= maxBruteForceVertices { 620 return l.bruteForceContainsPoint(p) 621 } 622 623 // Otherwise, look up the point in the index. 624 it := l.index.Iterator() 625 if !it.LocatePoint(p) { 626 return false 627 } 628 return l.iteratorContainsPoint(it, p) 629 } 630 631 // ContainsCell reports whether the given Cell is contained by this Loop. 632 func (l *Loop) ContainsCell(target Cell) bool { 633 it := l.index.Iterator() 634 relation := it.LocateCellID(target.ID()) 635 636 // If "target" is disjoint from all index cells, it is not contained. 637 // Similarly, if "target" is subdivided into one or more index cells then it 638 // is not contained, since index cells are subdivided only if they (nearly) 639 // intersect a sufficient number of edges. (But note that if "target" itself 640 // is an index cell then it may be contained, since it could be a cell with 641 // no edges in the loop interior.) 642 if relation != Indexed { 643 return false 644 } 645 646 // Otherwise check if any edges intersect "target". 647 if l.boundaryApproxIntersects(it, target) { 648 return false 649 } 650 651 // Otherwise check if the loop contains the center of "target". 652 return l.iteratorContainsPoint(it, target.Center()) 653 } 654 655 // IntersectsCell reports whether this Loop intersects the given cell. 656 func (l *Loop) IntersectsCell(target Cell) bool { 657 it := l.index.Iterator() 658 relation := it.LocateCellID(target.ID()) 659 660 // If target does not overlap any index cell, there is no intersection. 661 if relation == Disjoint { 662 return false 663 } 664 // If target is subdivided into one or more index cells, there is an 665 // intersection to within the ShapeIndex error bound (see Contains). 666 if relation == Subdivided { 667 return true 668 } 669 // If target is an index cell, there is an intersection because index cells 670 // are created only if they have at least one edge or they are entirely 671 // contained by the loop. 672 if it.CellID() == target.id { 673 return true 674 } 675 // Otherwise check if any edges intersect target. 676 if l.boundaryApproxIntersects(it, target) { 677 return true 678 } 679 // Otherwise check if the loop contains the center of target. 680 return l.iteratorContainsPoint(it, target.Center()) 681 } 682 683 // CellUnionBound computes a covering of the Loop. 684 func (l *Loop) CellUnionBound() []CellID { 685 return l.CapBound().CellUnionBound() 686 } 687 688 // boundaryApproxIntersects reports if the loop's boundary intersects target. 689 // It may also return true when the loop boundary does not intersect target but 690 // some edge comes within the worst-case error tolerance. 691 // 692 // This requires that it.Locate(target) returned Indexed. 693 func (l *Loop) boundaryApproxIntersects(it *ShapeIndexIterator, target Cell) bool { 694 aClipped := it.IndexCell().findByShapeID(0) 695 696 // If there are no edges, there is no intersection. 697 if len(aClipped.edges) == 0 { 698 return false 699 } 700 701 // We can save some work if target is the index cell itself. 702 if it.CellID() == target.ID() { 703 return true 704 } 705 706 // Otherwise check whether any of the edges intersect target. 707 maxError := (faceClipErrorUVCoord + intersectsRectErrorUVDist) 708 bound := target.BoundUV().ExpandedByMargin(maxError) 709 for _, ai := range aClipped.edges { 710 v0, v1, ok := ClipToPaddedFace(l.Vertex(ai), l.Vertex(ai+1), target.Face(), maxError) 711 if ok && edgeIntersectsRect(v0, v1, bound) { 712 return true 713 } 714 } 715 return false 716 } 717 718 // iteratorContainsPoint reports if the iterator that is positioned at the ShapeIndexCell 719 // that may contain p, contains the point p. 720 func (l *Loop) iteratorContainsPoint(it *ShapeIndexIterator, p Point) bool { 721 // Test containment by drawing a line segment from the cell center to the 722 // given point and counting edge crossings. 723 aClipped := it.IndexCell().findByShapeID(0) 724 inside := aClipped.containsCenter 725 if len(aClipped.edges) > 0 { 726 center := it.Center() 727 crosser := NewEdgeCrosser(center, p) 728 aiPrev := -2 729 for _, ai := range aClipped.edges { 730 if ai != aiPrev+1 { 731 crosser.RestartAt(l.Vertex(ai)) 732 } 733 aiPrev = ai 734 inside = inside != crosser.EdgeOrVertexChainCrossing(l.Vertex(ai+1)) 735 } 736 } 737 return inside 738 } 739 740 // RegularLoop creates a loop with the given number of vertices, all 741 // located on a circle of the specified radius around the given center. 742 func RegularLoop(center Point, radius s1.Angle, numVertices int) *Loop { 743 return RegularLoopForFrame(getFrame(center), radius, numVertices) 744 } 745 746 // RegularLoopForFrame creates a loop centered around the z-axis of the given 747 // coordinate frame, with the first vertex in the direction of the positive x-axis. 748 func RegularLoopForFrame(frame matrix3x3, radius s1.Angle, numVertices int) *Loop { 749 return LoopFromPoints(regularPointsForFrame(frame, radius, numVertices)) 750 } 751 752 // CanonicalFirstVertex returns a first index and a direction (either +1 or -1) 753 // such that the vertex sequence (first, first+dir, ..., first+(n-1)*dir) does 754 // not change when the loop vertex order is rotated or inverted. This allows the 755 // loop vertices to be traversed in a canonical order. The return values are 756 // chosen such that (first, ..., first+n*dir) are in the range [0, 2*n-1] as 757 // expected by the Vertex method. 758 func (l *Loop) CanonicalFirstVertex() (firstIdx, direction int) { 759 firstIdx = 0 760 n := len(l.vertices) 761 for i := 1; i < n; i++ { 762 if l.Vertex(i).Cmp(l.Vertex(firstIdx).Vector) == -1 { 763 firstIdx = i 764 } 765 } 766 767 // 0 <= firstIdx <= n-1, so (firstIdx+n*dir) <= 2*n-1. 768 if l.Vertex(firstIdx+1).Cmp(l.Vertex(firstIdx+n-1).Vector) == -1 { 769 return firstIdx, 1 770 } 771 772 // n <= firstIdx <= 2*n-1, so (firstIdx+n*dir) >= 0. 773 firstIdx += n 774 return firstIdx, -1 775 } 776 777 // TurningAngle returns the sum of the turning angles at each vertex. The return 778 // value is positive if the loop is counter-clockwise, negative if the loop is 779 // clockwise, and zero if the loop is a great circle. Degenerate and 780 // nearly-degenerate loops are handled consistently with Sign. So for example, 781 // if a loop has zero area (i.e., it is a very small CCW loop) then the turning 782 // angle will always be negative. 783 // 784 // This quantity is also called the "geodesic curvature" of the loop. 785 func (l *Loop) TurningAngle() float64 { 786 // For empty and full loops, we return the limit value as the loop area 787 // approaches 0 or 4*Pi respectively. 788 if l.isEmptyOrFull() { 789 if l.ContainsOrigin() { 790 return -2 * math.Pi 791 } 792 return 2 * math.Pi 793 } 794 795 // Don't crash even if the loop is not well-defined. 796 if len(l.vertices) < 3 { 797 return 0 798 } 799 800 // To ensure that we get the same result when the vertex order is rotated, 801 // and that the result is negated when the vertex order is reversed, we need 802 // to add up the individual turn angles in a consistent order. (In general, 803 // adding up a set of numbers in a different order can change the sum due to 804 // rounding errors.) 805 // 806 // Furthermore, if we just accumulate an ordinary sum then the worst-case 807 // error is quadratic in the number of vertices. (This can happen with 808 // spiral shapes, where the partial sum of the turning angles can be linear 809 // in the number of vertices.) To avoid this we use the Kahan summation 810 // algorithm (http://en.wikipedia.org/wiki/Kahan_summation_algorithm). 811 n := len(l.vertices) 812 i, dir := l.CanonicalFirstVertex() 813 sum := TurnAngle(l.Vertex((i+n-dir)%n), l.Vertex(i), l.Vertex((i+dir)%n)) 814 815 compensation := s1.Angle(0) 816 for n-1 > 0 { 817 i += dir 818 angle := TurnAngle(l.Vertex(i-dir), l.Vertex(i), l.Vertex(i+dir)) 819 oldSum := sum 820 angle += compensation 821 sum += angle 822 compensation = (oldSum - sum) + angle 823 n-- 824 } 825 826 const maxCurvature = 2*math.Pi - 4*dblEpsilon 827 828 return math.Max(-maxCurvature, math.Min(maxCurvature, float64(dir)*float64(sum+compensation))) 829 } 830 831 // turningAngleMaxError return the maximum error in TurningAngle. The value is not 832 // constant; it depends on the loop. 833 func (l *Loop) turningAngleMaxError() float64 { 834 // The maximum error can be bounded as follows: 835 // 3.00 * dblEpsilon for RobustCrossProd(b, a) 836 // 3.00 * dblEpsilon for RobustCrossProd(c, b) 837 // 3.25 * dblEpsilon for Angle() 838 // 2.00 * dblEpsilon for each addition in the Kahan summation 839 // ------------------ 840 // 11.25 * dblEpsilon 841 maxErrorPerVertex := 11.25 * dblEpsilon 842 return maxErrorPerVertex * float64(len(l.vertices)) 843 } 844 845 // IsHole reports whether this loop represents a hole in its containing polygon. 846 func (l *Loop) IsHole() bool { return l.depth&1 != 0 } 847 848 // Sign returns -1 if this Loop represents a hole in its containing polygon, and +1 otherwise. 849 func (l *Loop) Sign() int { 850 if l.IsHole() { 851 return -1 852 } 853 return 1 854 } 855 856 // IsNormalized reports whether the loop area is at most 2*pi. Degenerate loops are 857 // handled consistently with Sign, i.e., if a loop can be 858 // expressed as the union of degenerate or nearly-degenerate CCW triangles, 859 // then it will always be considered normalized. 860 func (l *Loop) IsNormalized() bool { 861 // Optimization: if the longitude span is less than 180 degrees, then the 862 // loop covers less than half the sphere and is therefore normalized. 863 if l.bound.Lng.Length() < math.Pi { 864 return true 865 } 866 867 // We allow some error so that hemispheres are always considered normalized. 868 // TODO(roberts): This is no longer required by the Polygon implementation, 869 // so alternatively we could create the invariant that a loop is normalized 870 // if and only if its complement is not normalized. 871 return l.TurningAngle() >= -l.turningAngleMaxError() 872 } 873 874 // Normalize inverts the loop if necessary so that the area enclosed by the loop 875 // is at most 2*pi. 876 func (l *Loop) Normalize() { 877 if !l.IsNormalized() { 878 l.Invert() 879 } 880 } 881 882 // Invert reverses the order of the loop vertices, effectively complementing the 883 // region represented by the loop. For example, the loop ABCD (with edges 884 // AB, BC, CD, DA) becomes the loop DCBA (with edges DC, CB, BA, AD). 885 // Notice that the last edge is the same in both cases except that its 886 // direction has been reversed. 887 func (l *Loop) Invert() { 888 l.index.Reset() 889 if l.isEmptyOrFull() { 890 if l.IsFull() { 891 l.vertices[0] = emptyLoopPoint 892 } else { 893 l.vertices[0] = fullLoopPoint 894 } 895 } else { 896 // For non-special loops, reverse the slice of vertices. 897 for i := len(l.vertices)/2 - 1; i >= 0; i-- { 898 opp := len(l.vertices) - 1 - i 899 l.vertices[i], l.vertices[opp] = l.vertices[opp], l.vertices[i] 900 } 901 } 902 903 // originInside must be set correctly before building the ShapeIndex. 904 l.originInside = !l.originInside 905 if l.bound.Lat.Lo > -math.Pi/2 && l.bound.Lat.Hi < math.Pi/2 { 906 // The complement of this loop contains both poles. 907 l.bound = FullRect() 908 l.subregionBound = l.bound 909 } else { 910 l.initBound() 911 } 912 l.index.Add(l) 913 } 914 915 // findVertex returns the index of the vertex at the given Point in the range 916 // 1..numVertices, and a boolean indicating if a vertex was found. 917 func (l *Loop) findVertex(p Point) (index int, ok bool) { 918 const notFound = 0 919 if len(l.vertices) < 10 { 920 // Exhaustive search for loops below a small threshold. 921 for i := 1; i <= len(l.vertices); i++ { 922 if l.Vertex(i) == p { 923 return i, true 924 } 925 } 926 return notFound, false 927 } 928 929 it := l.index.Iterator() 930 if !it.LocatePoint(p) { 931 return notFound, false 932 } 933 934 aClipped := it.IndexCell().findByShapeID(0) 935 for i := aClipped.numEdges() - 1; i >= 0; i-- { 936 ai := aClipped.edges[i] 937 if l.Vertex(ai) == p { 938 if ai == 0 { 939 return len(l.vertices), true 940 } 941 return ai, true 942 } 943 944 if l.Vertex(ai+1) == p { 945 return ai + 1, true 946 } 947 } 948 return notFound, false 949 } 950 951 // ContainsNested reports whether the given loops is contained within this loop. 952 // This function does not test for edge intersections. The two loops must meet 953 // all of the Polygon requirements; for example this implies that their 954 // boundaries may not cross or have any shared edges (although they may have 955 // shared vertices). 956 func (l *Loop) ContainsNested(other *Loop) bool { 957 if !l.subregionBound.Contains(other.bound) { 958 return false 959 } 960 961 // Special cases to handle either loop being empty or full. Also bail out 962 // when B has no vertices to avoid heap overflow on the vertex(1) call 963 // below. (This method is called during polygon initialization before the 964 // client has an opportunity to call IsValid().) 965 if l.isEmptyOrFull() || other.NumVertices() < 2 { 966 return l.IsFull() || other.IsEmpty() 967 } 968 969 // We are given that A and B do not share any edges, and that either one 970 // loop contains the other or they do not intersect. 971 m, ok := l.findVertex(other.Vertex(1)) 972 if !ok { 973 // Since other.vertex(1) is not shared, we can check whether A contains it. 974 return l.ContainsPoint(other.Vertex(1)) 975 } 976 977 // Check whether the edge order around other.Vertex(1) is compatible with 978 // A containing B. 979 return WedgeContains(l.Vertex(m-1), l.Vertex(m), l.Vertex(m+1), other.Vertex(0), other.Vertex(2)) 980 } 981 982 // surfaceIntegralFloat64 computes the oriented surface integral of some quantity f(x) 983 // over the loop interior, given a function f(A,B,C) that returns the 984 // corresponding integral over the spherical triangle ABC. Here "oriented 985 // surface integral" means: 986 // 987 // (1) f(A,B,C) must be the integral of f if ABC is counterclockwise, 988 // 989 // and the integral of -f if ABC is clockwise. 990 // 991 // (2) The result of this function is *either* the integral of f over the 992 // 993 // loop interior, or the integral of (-f) over the loop exterior. 994 // 995 // Note that there are at least two common situations where it easy to work 996 // around property (2) above: 997 // 998 // - If the integral of f over the entire sphere is zero, then it doesn't 999 // matter which case is returned because they are always equal. 1000 // 1001 // - If f is non-negative, then it is easy to detect when the integral over 1002 // the loop exterior has been returned, and the integral over the loop 1003 // interior can be obtained by adding the integral of f over the entire 1004 // unit sphere (a constant) to the result. 1005 // 1006 // Any changes to this method may need corresponding changes to surfaceIntegralPoint as well. 1007 func (l *Loop) surfaceIntegralFloat64(f func(a, b, c Point) float64) float64 { 1008 // We sum f over a collection T of oriented triangles, possibly 1009 // overlapping. Let the sign of a triangle be +1 if it is CCW and -1 1010 // otherwise, and let the sign of a point x be the sum of the signs of the 1011 // triangles containing x. Then the collection of triangles T is chosen 1012 // such that either: 1013 // 1014 // (1) Each point in the loop interior has sign +1, and sign 0 otherwise; or 1015 // (2) Each point in the loop exterior has sign -1, and sign 0 otherwise. 1016 // 1017 // The triangles basically consist of a fan from vertex 0 to every loop 1018 // edge that does not include vertex 0. These triangles will always satisfy 1019 // either (1) or (2). However, what makes this a bit tricky is that 1020 // spherical edges become numerically unstable as their length approaches 1021 // 180 degrees. Of course there is not much we can do if the loop itself 1022 // contains such edges, but we would like to make sure that all the triangle 1023 // edges under our control (i.e., the non-loop edges) are stable. For 1024 // example, consider a loop around the equator consisting of four equally 1025 // spaced points. This is a well-defined loop, but we cannot just split it 1026 // into two triangles by connecting vertex 0 to vertex 2. 1027 // 1028 // We handle this type of situation by moving the origin of the triangle fan 1029 // whenever we are about to create an unstable edge. We choose a new 1030 // location for the origin such that all relevant edges are stable. We also 1031 // create extra triangles with the appropriate orientation so that the sum 1032 // of the triangle signs is still correct at every point. 1033 1034 // The maximum length of an edge for it to be considered numerically stable. 1035 // The exact value is fairly arbitrary since it depends on the stability of 1036 // the function f. The value below is quite conservative but could be 1037 // reduced further if desired. 1038 const maxLength = math.Pi - 1e-5 1039 1040 var sum float64 1041 origin := l.Vertex(0) 1042 for i := 1; i+1 < len(l.vertices); i++ { 1043 // Let V_i be vertex(i), let O be the current origin, and let length(A,B) 1044 // be the length of edge (A,B). At the start of each loop iteration, the 1045 // "leading edge" of the triangle fan is (O,V_i), and we want to extend 1046 // the triangle fan so that the leading edge is (O,V_i+1). 1047 // 1048 // Invariants: 1049 // 1. length(O,V_i) < maxLength for all (i > 1). 1050 // 2. Either O == V_0, or O is approximately perpendicular to V_0. 1051 // 3. "sum" is the oriented integral of f over the area defined by 1052 // (O, V_0, V_1, ..., V_i). 1053 if l.Vertex(i+1).Angle(origin.Vector) > maxLength { 1054 // We are about to create an unstable edge, so choose a new origin O' 1055 // for the triangle fan. 1056 oldOrigin := origin 1057 if origin == l.Vertex(0) { 1058 // The following point is well-separated from V_i and V_0 (and 1059 // therefore V_i+1 as well). 1060 origin = Point{l.Vertex(0).PointCross(l.Vertex(i)).Normalize()} 1061 } else if l.Vertex(i).Angle(l.Vertex(0).Vector) < maxLength { 1062 // All edges of the triangle (O, V_0, V_i) are stable, so we can 1063 // revert to using V_0 as the origin. 1064 origin = l.Vertex(0) 1065 } else { 1066 // (O, V_i+1) and (V_0, V_i) are antipodal pairs, and O and V_0 are 1067 // perpendicular. Therefore V_0.CrossProd(O) is approximately 1068 // perpendicular to all of {O, V_0, V_i, V_i+1}, and we can choose 1069 // this point O' as the new origin. 1070 origin = Point{l.Vertex(0).Cross(oldOrigin.Vector)} 1071 1072 // Advance the edge (V_0,O) to (V_0,O'). 1073 sum += f(l.Vertex(0), oldOrigin, origin) 1074 } 1075 // Advance the edge (O,V_i) to (O',V_i). 1076 sum += f(oldOrigin, l.Vertex(i), origin) 1077 } 1078 // Advance the edge (O,V_i) to (O,V_i+1). 1079 sum += f(origin, l.Vertex(i), l.Vertex(i+1)) 1080 } 1081 // If the origin is not V_0, we need to sum one more triangle. 1082 if origin != l.Vertex(0) { 1083 // Advance the edge (O,V_n-1) to (O,V_0). 1084 sum += f(origin, l.Vertex(len(l.vertices)-1), l.Vertex(0)) 1085 } 1086 return sum 1087 } 1088 1089 // surfaceIntegralPoint mirrors the surfaceIntegralFloat64 method but over Points; 1090 // see that method for commentary. The C++ version uses a templated method. 1091 // Any changes to this method may need corresponding changes to surfaceIntegralFloat64 as well. 1092 func (l *Loop) surfaceIntegralPoint(f func(a, b, c Point) Point) Point { 1093 const maxLength = math.Pi - 1e-5 1094 var sum r3.Vector 1095 1096 origin := l.Vertex(0) 1097 for i := 1; i+1 < len(l.vertices); i++ { 1098 if l.Vertex(i+1).Angle(origin.Vector) > maxLength { 1099 oldOrigin := origin 1100 if origin == l.Vertex(0) { 1101 origin = Point{l.Vertex(0).PointCross(l.Vertex(i)).Normalize()} 1102 } else if l.Vertex(i).Angle(l.Vertex(0).Vector) < maxLength { 1103 origin = l.Vertex(0) 1104 } else { 1105 origin = Point{l.Vertex(0).Cross(oldOrigin.Vector)} 1106 sum = sum.Add(f(l.Vertex(0), oldOrigin, origin).Vector) 1107 } 1108 sum = sum.Add(f(oldOrigin, l.Vertex(i), origin).Vector) 1109 } 1110 sum = sum.Add(f(origin, l.Vertex(i), l.Vertex(i+1)).Vector) 1111 } 1112 if origin != l.Vertex(0) { 1113 sum = sum.Add(f(origin, l.Vertex(len(l.vertices)-1), l.Vertex(0)).Vector) 1114 } 1115 return Point{sum} 1116 } 1117 1118 // Area returns the area of the loop interior, i.e. the region on the left side of 1119 // the loop. The return value is between 0 and 4*pi. (Note that the return 1120 // value is not affected by whether this loop is a "hole" or a "shell".) 1121 func (l *Loop) Area() float64 { 1122 // It is surprisingly difficult to compute the area of a loop robustly. The 1123 // main issues are (1) whether degenerate loops are considered to be CCW or 1124 // not (i.e., whether their area is close to 0 or 4*pi), and (2) computing 1125 // the areas of small loops with good relative accuracy. 1126 // 1127 // With respect to degeneracies, we would like Area to be consistent 1128 // with ContainsPoint in that loops that contain many points 1129 // should have large areas, and loops that contain few points should have 1130 // small areas. For example, if a degenerate triangle is considered CCW 1131 // according to s2predicates Sign, then it will contain very few points and 1132 // its area should be approximately zero. On the other hand if it is 1133 // considered clockwise, then it will contain virtually all points and so 1134 // its area should be approximately 4*pi. 1135 // 1136 // More precisely, let U be the set of Points for which IsUnitLength 1137 // is true, let P(U) be the projection of those points onto the mathematical 1138 // unit sphere, and let V(P(U)) be the Voronoi diagram of the projected 1139 // points. Then for every loop x, we would like Area to approximately 1140 // equal the sum of the areas of the Voronoi regions of the points p for 1141 // which x.ContainsPoint(p) is true. 1142 // 1143 // The second issue is that we want to compute the area of small loops 1144 // accurately. This requires having good relative precision rather than 1145 // good absolute precision. For example, if the area of a loop is 1e-12 and 1146 // the error is 1e-15, then the area only has 3 digits of accuracy. (For 1147 // reference, 1e-12 is about 40 square meters on the surface of the earth.) 1148 // We would like to have good relative accuracy even for small loops. 1149 // 1150 // To achieve these goals, we combine two different methods of computing the 1151 // area. This first method is based on the Gauss-Bonnet theorem, which says 1152 // that the area enclosed by the loop equals 2*pi minus the total geodesic 1153 // curvature of the loop (i.e., the sum of the "turning angles" at all the 1154 // loop vertices). The big advantage of this method is that as long as we 1155 // use Sign to compute the turning angle at each vertex, then 1156 // degeneracies are always handled correctly. In other words, if a 1157 // degenerate loop is CCW according to the symbolic perturbations used by 1158 // Sign, then its turning angle will be approximately 2*pi. 1159 // 1160 // The disadvantage of the Gauss-Bonnet method is that its absolute error is 1161 // about 2e-15 times the number of vertices (see turningAngleMaxError). 1162 // So, it cannot compute the area of small loops accurately. 1163 // 1164 // The second method is based on splitting the loop into triangles and 1165 // summing the area of each triangle. To avoid the difficulty and expense 1166 // of decomposing the loop into a union of non-overlapping triangles, 1167 // instead we compute a signed sum over triangles that may overlap (see the 1168 // comments for surfaceIntegral). The advantage of this method 1169 // is that the area of each triangle can be computed with much better 1170 // relative accuracy (using l'Huilier's theorem). The disadvantage is that 1171 // the result is a signed area: CCW loops may yield a small positive value, 1172 // while CW loops may yield a small negative value (which is converted to a 1173 // positive area by adding 4*pi). This means that small errors in computing 1174 // the signed area may translate into a very large error in the result (if 1175 // the sign of the sum is incorrect). 1176 // 1177 // So, our strategy is to combine these two methods as follows. First we 1178 // compute the area using the "signed sum over triangles" approach (since it 1179 // is generally more accurate). We also estimate the maximum error in this 1180 // result. If the signed area is too close to zero (i.e., zero is within 1181 // the error bounds), then we double-check the sign of the result using the 1182 // Gauss-Bonnet method. (In fact we just call IsNormalized, which is 1183 // based on this method.) If the two methods disagree, we return either 0 1184 // or 4*pi based on the result of IsNormalized. Otherwise we return the 1185 // area that we computed originally. 1186 if l.isEmptyOrFull() { 1187 if l.ContainsOrigin() { 1188 return 4 * math.Pi 1189 } 1190 return 0 1191 } 1192 area := l.surfaceIntegralFloat64(SignedArea) 1193 1194 // TODO(roberts): This error estimate is very approximate. There are two 1195 // issues: (1) SignedArea needs some improvements to ensure that its error 1196 // is actually never higher than GirardArea, and (2) although the number of 1197 // triangles in the sum is typically N-2, in theory it could be as high as 1198 // 2*N for pathological inputs. But in other respects this error bound is 1199 // very conservative since it assumes that the maximum error is achieved on 1200 // every triangle. 1201 maxError := l.turningAngleMaxError() 1202 1203 // The signed area should be between approximately -4*pi and 4*pi. 1204 if area < 0 { 1205 // We have computed the negative of the area of the loop exterior. 1206 area += 4 * math.Pi 1207 } 1208 1209 if area > 4*math.Pi { 1210 area = 4 * math.Pi 1211 } 1212 if area < 0 { 1213 area = 0 1214 } 1215 1216 // If the area is close enough to zero or 4*pi so that the loop orientation 1217 // is ambiguous, then we compute the loop orientation explicitly. 1218 if area < maxError && !l.IsNormalized() { 1219 return 4 * math.Pi 1220 } else if area > (4*math.Pi-maxError) && l.IsNormalized() { 1221 return 0 1222 } 1223 1224 return area 1225 } 1226 1227 // Centroid returns the true centroid of the loop multiplied by the area of the 1228 // loop. The result is not unit length, so you may want to normalize it. Also 1229 // note that in general, the centroid may not be contained by the loop. 1230 // 1231 // We prescale by the loop area for two reasons: (1) it is cheaper to 1232 // compute this way, and (2) it makes it easier to compute the centroid of 1233 // more complicated shapes (by splitting them into disjoint regions and 1234 // adding their centroids). 1235 // 1236 // Note that the return value is not affected by whether this loop is a 1237 // "hole" or a "shell". 1238 func (l *Loop) Centroid() Point { 1239 // surfaceIntegralPoint() returns either the integral of position over loop 1240 // interior, or the negative of the integral of position over the loop 1241 // exterior. But these two values are the same (!), because the integral of 1242 // position over the entire sphere is (0, 0, 0). 1243 return l.surfaceIntegralPoint(TrueCentroid) 1244 } 1245 1246 // Encode encodes the Loop. 1247 func (l Loop) Encode(w io.Writer) error { 1248 e := &encoder{w: w} 1249 l.encode(e) 1250 return e.err 1251 } 1252 1253 func (l Loop) encode(e *encoder) { 1254 e.writeInt8(encodingVersion) 1255 e.writeUint32(uint32(len(l.vertices))) 1256 for _, v := range l.vertices { 1257 e.writeFloat64(v.X) 1258 e.writeFloat64(v.Y) 1259 e.writeFloat64(v.Z) 1260 } 1261 1262 e.writeBool(l.originInside) 1263 e.writeInt32(int32(l.depth)) 1264 1265 // Encode the bound. 1266 l.bound.encode(e) 1267 } 1268 1269 // Decode decodes a loop. 1270 func (l *Loop) Decode(r io.Reader) error { 1271 *l = Loop{} 1272 d := &decoder{r: asByteReader(r)} 1273 l.decode(d) 1274 return d.err 1275 } 1276 1277 func (l *Loop) decode(d *decoder) { 1278 version := int8(d.readUint8()) 1279 if d.err != nil { 1280 return 1281 } 1282 if version != encodingVersion { 1283 d.err = fmt.Errorf("cannot decode version %d", version) 1284 return 1285 } 1286 1287 // Empty loops are explicitly allowed here: a newly created loop has zero vertices 1288 // and such loops encode and decode properly. 1289 nvertices := d.readUint32() 1290 if nvertices > maxEncodedVertices { 1291 if d.err == nil { 1292 d.err = fmt.Errorf("too many vertices (%d; max is %d)", nvertices, maxEncodedVertices) 1293 1294 } 1295 return 1296 } 1297 l.vertices = make([]Point, nvertices) 1298 for i := range l.vertices { 1299 l.vertices[i].X = d.readFloat64() 1300 l.vertices[i].Y = d.readFloat64() 1301 l.vertices[i].Z = d.readFloat64() 1302 } 1303 l.index = NewShapeIndex() 1304 l.originInside = d.readBool() 1305 l.depth = int(d.readUint32()) 1306 l.bound.decode(d) 1307 l.subregionBound = ExpandForSubregions(l.bound) 1308 1309 l.index.Add(l) 1310 } 1311 1312 // Bitmasks to read from properties. 1313 const ( 1314 originInside = 1 << iota 1315 boundEncoded 1316 ) 1317 1318 func (l *Loop) xyzFaceSiTiVertices() []xyzFaceSiTi { 1319 ret := make([]xyzFaceSiTi, len(l.vertices)) 1320 for i, v := range l.vertices { 1321 ret[i].xyz = v 1322 ret[i].face, ret[i].si, ret[i].ti, ret[i].level = xyzToFaceSiTi(v) 1323 } 1324 return ret 1325 } 1326 1327 func (l *Loop) encodeCompressed(e *encoder, snapLevel int, vertices []xyzFaceSiTi) { 1328 if len(l.vertices) != len(vertices) { 1329 panic("encodeCompressed: vertices must be the same length as l.vertices") 1330 } 1331 if len(vertices) > maxEncodedVertices { 1332 if e.err == nil { 1333 e.err = fmt.Errorf("too many vertices (%d; max is %d)", len(vertices), maxEncodedVertices) 1334 } 1335 return 1336 } 1337 e.writeUvarint(uint64(len(vertices))) 1338 encodePointsCompressed(e, vertices, snapLevel) 1339 1340 props := l.compressedEncodingProperties() 1341 e.writeUvarint(props) 1342 e.writeUvarint(uint64(l.depth)) 1343 if props&boundEncoded != 0 { 1344 l.bound.encode(e) 1345 } 1346 } 1347 1348 func (l *Loop) compressedEncodingProperties() uint64 { 1349 var properties uint64 1350 if l.originInside { 1351 properties |= originInside 1352 } 1353 1354 // Write whether there is a bound so we can change the threshold later. 1355 // Recomputing the bound multiplies the decode time taken per vertex 1356 // by a factor of about 3.5. Without recomputing the bound, decode 1357 // takes approximately 125 ns / vertex. A loop with 63 vertices 1358 // encoded without the bound will take ~30us to decode, which is 1359 // acceptable. At ~3.5 bytes / vertex without the bound, adding 1360 // the bound will increase the size by <15%, which is also acceptable. 1361 const minVerticesForBound = 64 1362 if len(l.vertices) >= minVerticesForBound { 1363 properties |= boundEncoded 1364 } 1365 1366 return properties 1367 } 1368 1369 func (l *Loop) decodeCompressed(d *decoder, snapLevel int) { 1370 nvertices := d.readUvarint() 1371 if d.err != nil { 1372 return 1373 } 1374 if nvertices > maxEncodedVertices { 1375 d.err = fmt.Errorf("too many vertices (%d; max is %d)", nvertices, maxEncodedVertices) 1376 return 1377 } 1378 l.vertices = make([]Point, nvertices) 1379 decodePointsCompressed(d, snapLevel, l.vertices) 1380 properties := d.readUvarint() 1381 1382 // Make sure values are valid before using. 1383 if d.err != nil { 1384 return 1385 } 1386 1387 l.index = NewShapeIndex() 1388 l.originInside = (properties & originInside) != 0 1389 1390 l.depth = int(d.readUvarint()) 1391 1392 if (properties & boundEncoded) != 0 { 1393 l.bound.decode(d) 1394 if d.err != nil { 1395 return 1396 } 1397 l.subregionBound = ExpandForSubregions(l.bound) 1398 } else { 1399 l.initBound() 1400 } 1401 1402 l.index.Add(l) 1403 } 1404 1405 // crossingTarget is an enum representing the possible crossing target cases for relations. 1406 type crossingTarget int 1407 1408 const ( 1409 crossingTargetDontCare crossingTarget = iota 1410 crossingTargetDontCross 1411 crossingTargetCross 1412 ) 1413 1414 // loopRelation defines the interface for checking a type of relationship between two loops. 1415 // Some examples of relations are Contains, Intersects, or CompareBoundary. 1416 type loopRelation interface { 1417 // Optionally, aCrossingTarget and bCrossingTarget can specify an early-exit 1418 // condition for the loop relation. If any point P is found such that 1419 // 1420 // A.ContainsPoint(P) == aCrossingTarget() && 1421 // B.ContainsPoint(P) == bCrossingTarget() 1422 // 1423 // then the loop relation is assumed to be the same as if a pair of crossing 1424 // edges were found. For example, the ContainsPoint relation has 1425 // 1426 // aCrossingTarget() == crossingTargetDontCross 1427 // bCrossingTarget() == crossingTargetCross 1428 // 1429 // because if A.ContainsPoint(P) == false and B.ContainsPoint(P) == true 1430 // for any point P, then it is equivalent to finding an edge crossing (i.e., 1431 // since Contains returns false in both cases). 1432 // 1433 // Loop relations that do not have an early-exit condition of this form 1434 // should return crossingTargetDontCare for both crossing targets. 1435 1436 // aCrossingTarget reports whether loop A crosses the target point with 1437 // the given relation type. 1438 aCrossingTarget() crossingTarget 1439 // bCrossingTarget reports whether loop B crosses the target point with 1440 // the given relation type. 1441 bCrossingTarget() crossingTarget 1442 1443 // wedgesCross reports if a shared vertex ab1 and the two associated wedges 1444 // (a0, ab1, b2) and (b0, ab1, b2) are equivalent to an edge crossing. 1445 // The loop relation is also allowed to maintain its own internal state, and 1446 // can return true if it observes any sequence of wedges that are equivalent 1447 // to an edge crossing. 1448 wedgesCross(a0, ab1, a2, b0, b2 Point) bool 1449 } 1450 1451 // loopCrosser is a helper type for determining whether two loops cross. 1452 // It is instantiated twice for each pair of loops to be tested, once for the 1453 // pair (A,B) and once for the pair (B,A), in order to be able to process 1454 // edges in either loop nesting order. 1455 type loopCrosser struct { 1456 a, b *Loop 1457 relation loopRelation 1458 swapped bool 1459 aCrossingTarget crossingTarget 1460 bCrossingTarget crossingTarget 1461 1462 // state maintained by startEdge and edgeCrossesCell. 1463 crosser *EdgeCrosser 1464 aj, bjPrev int 1465 1466 // temporary data declared here to avoid repeated memory allocations. 1467 bQuery *CrossingEdgeQuery 1468 bCells []*ShapeIndexCell 1469 } 1470 1471 // newLoopCrosser creates a loopCrosser from the given values. If swapped is true, 1472 // the loops A and B have been swapped. This affects how arguments are passed to 1473 // the given loop relation, since for example A.Contains(B) is not the same as 1474 // B.Contains(A). 1475 func newLoopCrosser(a, b *Loop, relation loopRelation, swapped bool) *loopCrosser { 1476 l := &loopCrosser{ 1477 a: a, 1478 b: b, 1479 relation: relation, 1480 swapped: swapped, 1481 aCrossingTarget: relation.aCrossingTarget(), 1482 bCrossingTarget: relation.bCrossingTarget(), 1483 bQuery: NewCrossingEdgeQuery(b.index), 1484 } 1485 if swapped { 1486 l.aCrossingTarget, l.bCrossingTarget = l.bCrossingTarget, l.aCrossingTarget 1487 } 1488 1489 return l 1490 } 1491 1492 // startEdge sets the crossers state for checking the given edge of loop A. 1493 func (l *loopCrosser) startEdge(aj int) { 1494 l.crosser = NewEdgeCrosser(l.a.Vertex(aj), l.a.Vertex(aj+1)) 1495 l.aj = aj 1496 l.bjPrev = -2 1497 } 1498 1499 // edgeCrossesCell reports whether the current edge of loop A has any crossings with 1500 // edges of the index cell of loop B. 1501 func (l *loopCrosser) edgeCrossesCell(bClipped *clippedShape) bool { 1502 // Test the current edge of A against all edges of bClipped 1503 bNumEdges := bClipped.numEdges() 1504 for j := 0; j < bNumEdges; j++ { 1505 bj := bClipped.edges[j] 1506 if bj != l.bjPrev+1 { 1507 l.crosser.RestartAt(l.b.Vertex(bj)) 1508 } 1509 l.bjPrev = bj 1510 if crossing := l.crosser.ChainCrossingSign(l.b.Vertex(bj + 1)); crossing == DoNotCross { 1511 continue 1512 } else if crossing == Cross { 1513 return true 1514 } 1515 1516 // We only need to check each shared vertex once, so we only 1517 // consider the case where l.aVertex(l.aj+1) == l.b.Vertex(bj+1). 1518 if l.a.Vertex(l.aj+1) == l.b.Vertex(bj+1) { 1519 if l.swapped { 1520 if l.relation.wedgesCross(l.b.Vertex(bj), l.b.Vertex(bj+1), l.b.Vertex(bj+2), l.a.Vertex(l.aj), l.a.Vertex(l.aj+2)) { 1521 return true 1522 } 1523 } else { 1524 if l.relation.wedgesCross(l.a.Vertex(l.aj), l.a.Vertex(l.aj+1), l.a.Vertex(l.aj+2), l.b.Vertex(bj), l.b.Vertex(bj+2)) { 1525 return true 1526 } 1527 } 1528 } 1529 } 1530 1531 return false 1532 } 1533 1534 // cellCrossesCell reports whether there are any edge crossings or wedge crossings 1535 // within the two given cells. 1536 func (l *loopCrosser) cellCrossesCell(aClipped, bClipped *clippedShape) bool { 1537 // Test all edges of aClipped against all edges of bClipped. 1538 for _, edge := range aClipped.edges { 1539 l.startEdge(edge) 1540 if l.edgeCrossesCell(bClipped) { 1541 return true 1542 } 1543 } 1544 1545 return false 1546 } 1547 1548 // cellCrossesAnySubcell reports whether given an index cell of A, if there are any 1549 // edge or wedge crossings with any index cell of B contained within bID. 1550 func (l *loopCrosser) cellCrossesAnySubcell(aClipped *clippedShape, bID CellID) bool { 1551 // Test all edges of aClipped against all edges of B. The relevant B 1552 // edges are guaranteed to be children of bID, which lets us find the 1553 // correct index cells more efficiently. 1554 bRoot := PaddedCellFromCellID(bID, 0) 1555 for _, aj := range aClipped.edges { 1556 // Use an CrossingEdgeQuery starting at bRoot to find the index cells 1557 // of B that might contain crossing edges. 1558 l.bCells = l.bQuery.getCells(l.a.Vertex(aj), l.a.Vertex(aj+1), bRoot) 1559 if len(l.bCells) == 0 { 1560 continue 1561 } 1562 l.startEdge(aj) 1563 for c := 0; c < len(l.bCells); c++ { 1564 if l.edgeCrossesCell(l.bCells[c].shapes[0]) { 1565 return true 1566 } 1567 } 1568 } 1569 1570 return false 1571 } 1572 1573 // hasCrossing reports whether given two iterators positioned such that 1574 // ai.cellID().ContainsCellID(bi.cellID()), there is an edge or wedge crossing 1575 // anywhere within ai.cellID(). This function advances bi only past ai.cellID(). 1576 func (l *loopCrosser) hasCrossing(ai, bi *rangeIterator) bool { 1577 // If ai.CellID() intersects many edges of B, then it is faster to use 1578 // CrossingEdgeQuery to narrow down the candidates. But if it intersects 1579 // only a few edges, it is faster to check all the crossings directly. 1580 // We handle this by advancing bi and keeping track of how many edges we 1581 // would need to test. 1582 const edgeQueryMinEdges = 20 // Tuned from benchmarks. 1583 var totalEdges int 1584 l.bCells = nil 1585 1586 for { 1587 if n := bi.it.IndexCell().shapes[0].numEdges(); n > 0 { 1588 totalEdges += n 1589 if totalEdges >= edgeQueryMinEdges { 1590 // There are too many edges to test them directly, so use CrossingEdgeQuery. 1591 if l.cellCrossesAnySubcell(ai.it.IndexCell().shapes[0], ai.cellID()) { 1592 return true 1593 } 1594 bi.seekBeyond(ai) 1595 return false 1596 } 1597 l.bCells = append(l.bCells, bi.indexCell()) 1598 } 1599 bi.next() 1600 if bi.cellID() > ai.rangeMax { 1601 break 1602 } 1603 } 1604 1605 // Test all the edge crossings directly. 1606 for _, c := range l.bCells { 1607 if l.cellCrossesCell(ai.it.IndexCell().shapes[0], c.shapes[0]) { 1608 return true 1609 } 1610 } 1611 1612 return false 1613 } 1614 1615 // containsCenterMatches reports if the clippedShapes containsCenter boolean corresponds 1616 // to the crossing target type given. (This is to work around C++ allowing false == 0, 1617 // true == 1 type implicit conversions and comparisons) 1618 func containsCenterMatches(a *clippedShape, target crossingTarget) bool { 1619 return (!a.containsCenter && target == crossingTargetDontCross) || 1620 (a.containsCenter && target == crossingTargetCross) 1621 } 1622 1623 // hasCrossingRelation reports whether given two iterators positioned such that 1624 // ai.cellID().ContainsCellID(bi.cellID()), there is a crossing relationship 1625 // anywhere within ai.cellID(). Specifically, this method returns true if there 1626 // is an edge crossing, a wedge crossing, or a point P that matches both relations 1627 // crossing targets. This function advances both iterators past ai.cellID. 1628 func (l *loopCrosser) hasCrossingRelation(ai, bi *rangeIterator) bool { 1629 aClipped := ai.it.IndexCell().shapes[0] 1630 if aClipped.numEdges() != 0 { 1631 // The current cell of A has at least one edge, so check for crossings. 1632 if l.hasCrossing(ai, bi) { 1633 return true 1634 } 1635 ai.next() 1636 return false 1637 } 1638 1639 if containsCenterMatches(aClipped, l.aCrossingTarget) { 1640 // The crossing target for A is not satisfied, so we skip over these cells of B. 1641 bi.seekBeyond(ai) 1642 ai.next() 1643 return false 1644 } 1645 1646 // All points within ai.cellID() satisfy the crossing target for A, so it's 1647 // worth iterating through the cells of B to see whether any cell 1648 // centers also satisfy the crossing target for B. 1649 for bi.cellID() <= ai.rangeMax { 1650 bClipped := bi.it.IndexCell().shapes[0] 1651 if containsCenterMatches(bClipped, l.bCrossingTarget) { 1652 return true 1653 } 1654 bi.next() 1655 } 1656 ai.next() 1657 return false 1658 } 1659 1660 // hasCrossingRelation checks all edges of loop A for intersection against all edges 1661 // of loop B and reports if there are any that satisfy the given relation. If there 1662 // is any shared vertex, the wedges centered at this vertex are sent to the given 1663 // relation to be tested. 1664 // 1665 // If the two loop boundaries cross, this method is guaranteed to return 1666 // true. It also returns true in certain cases if the loop relationship is 1667 // equivalent to crossing. For example, if the relation is Contains and a 1668 // point P is found such that B contains P but A does not contain P, this 1669 // method will return true to indicate that the result is the same as though 1670 // a pair of crossing edges were found (since Contains returns false in 1671 // both cases). 1672 // 1673 // See Contains, Intersects and CompareBoundary for the three uses of this function. 1674 func hasCrossingRelation(a, b *Loop, relation loopRelation) bool { 1675 // We look for CellID ranges where the indexes of A and B overlap, and 1676 // then test those edges for crossings. 1677 ai := newRangeIterator(a.index) 1678 bi := newRangeIterator(b.index) 1679 1680 ab := newLoopCrosser(a, b, relation, false) // Tests edges of A against B 1681 ba := newLoopCrosser(b, a, relation, true) // Tests edges of B against A 1682 1683 for !ai.done() || !bi.done() { 1684 if ai.rangeMax < bi.rangeMin { 1685 // The A and B cells don't overlap, and A precedes B. 1686 ai.seekTo(bi) 1687 } else if bi.rangeMax < ai.rangeMin { 1688 // The A and B cells don't overlap, and B precedes A. 1689 bi.seekTo(ai) 1690 } else { 1691 // One cell contains the other. Determine which cell is larger. 1692 abRelation := int64(ai.it.CellID().lsb() - bi.it.CellID().lsb()) 1693 if abRelation > 0 { 1694 // A's index cell is larger. 1695 if ab.hasCrossingRelation(ai, bi) { 1696 return true 1697 } 1698 } else if abRelation < 0 { 1699 // B's index cell is larger. 1700 if ba.hasCrossingRelation(bi, ai) { 1701 return true 1702 } 1703 } else { 1704 // The A and B cells are the same. Since the two cells 1705 // have the same center point P, check whether P satisfies 1706 // the crossing targets. 1707 aClipped := ai.it.IndexCell().shapes[0] 1708 bClipped := bi.it.IndexCell().shapes[0] 1709 if containsCenterMatches(aClipped, ab.aCrossingTarget) && 1710 containsCenterMatches(bClipped, ab.bCrossingTarget) { 1711 return true 1712 } 1713 // Otherwise test all the edge crossings directly. 1714 if aClipped.numEdges() > 0 && bClipped.numEdges() > 0 && ab.cellCrossesCell(aClipped, bClipped) { 1715 return true 1716 } 1717 ai.next() 1718 bi.next() 1719 } 1720 } 1721 } 1722 return false 1723 } 1724 1725 // containsRelation implements loopRelation for a contains operation. If 1726 // A.ContainsPoint(P) == false && B.ContainsPoint(P) == true, it is equivalent 1727 // to having an edge crossing (i.e., Contains returns false). 1728 type containsRelation struct { 1729 foundSharedVertex bool 1730 } 1731 1732 func (c *containsRelation) aCrossingTarget() crossingTarget { return crossingTargetDontCross } 1733 func (c *containsRelation) bCrossingTarget() crossingTarget { return crossingTargetCross } 1734 func (c *containsRelation) wedgesCross(a0, ab1, a2, b0, b2 Point) bool { 1735 c.foundSharedVertex = true 1736 return !WedgeContains(a0, ab1, a2, b0, b2) 1737 } 1738 1739 // intersectsRelation implements loopRelation for an intersects operation. Given 1740 // two loops, A and B, if A.ContainsPoint(P) == true && B.ContainsPoint(P) == true, 1741 // it is equivalent to having an edge crossing (i.e., Intersects returns true). 1742 type intersectsRelation struct { 1743 foundSharedVertex bool 1744 } 1745 1746 func (i *intersectsRelation) aCrossingTarget() crossingTarget { return crossingTargetCross } 1747 func (i *intersectsRelation) bCrossingTarget() crossingTarget { return crossingTargetCross } 1748 func (i *intersectsRelation) wedgesCross(a0, ab1, a2, b0, b2 Point) bool { 1749 i.foundSharedVertex = true 1750 return WedgeIntersects(a0, ab1, a2, b0, b2) 1751 } 1752 1753 // compareBoundaryRelation implements loopRelation for comparing boundaries. 1754 // 1755 // The compare boundary relation does not have a useful early-exit condition, 1756 // so we return crossingTargetDontCare for both crossing targets. 1757 // 1758 // Aside: A possible early exit condition could be based on the following. 1759 // 1760 // If A contains a point of both B and ~B, then A intersects Boundary(B). 1761 // If ~A contains a point of both B and ~B, then ~A intersects Boundary(B). 1762 // So if the intersections of {A, ~A} with {B, ~B} are all non-empty, 1763 // the return value is 0, i.e., Boundary(A) intersects Boundary(B). 1764 // 1765 // Unfortunately it isn't worth detecting this situation because by the 1766 // time we have seen a point in all four intersection regions, we are also 1767 // guaranteed to have seen at least one pair of crossing edges. 1768 type compareBoundaryRelation struct { 1769 reverse bool // True if the other loop should be reversed. 1770 foundSharedVertex bool // True if any wedge was processed. 1771 containsEdge bool // True if any edge of the other loop is contained by this loop. 1772 excludesEdge bool // True if any edge of the other loop is excluded by this loop. 1773 } 1774 1775 func newCompareBoundaryRelation(reverse bool) *compareBoundaryRelation { 1776 return &compareBoundaryRelation{reverse: reverse} 1777 } 1778 1779 func (c *compareBoundaryRelation) aCrossingTarget() crossingTarget { return crossingTargetDontCare } 1780 func (c *compareBoundaryRelation) bCrossingTarget() crossingTarget { return crossingTargetDontCare } 1781 func (c *compareBoundaryRelation) wedgesCross(a0, ab1, a2, b0, b2 Point) bool { 1782 // Because we don't care about the interior of the other, only its boundary, 1783 // it is sufficient to check whether this one contains the semiwedge (ab1, b2). 1784 c.foundSharedVertex = true 1785 if wedgeContainsSemiwedge(a0, ab1, a2, b2, c.reverse) { 1786 c.containsEdge = true 1787 } else { 1788 c.excludesEdge = true 1789 } 1790 return c.containsEdge && c.excludesEdge 1791 } 1792 1793 // wedgeContainsSemiwedge reports whether the wedge (a0, ab1, a2) contains the 1794 // "semiwedge" defined as any non-empty open set of rays immediately CCW from 1795 // the edge (ab1, b2). If reverse is true, then substitute clockwise for CCW; 1796 // this simulates what would happen if the direction of the other loop was reversed. 1797 func wedgeContainsSemiwedge(a0, ab1, a2, b2 Point, reverse bool) bool { 1798 if b2 == a0 || b2 == a2 { 1799 // We have a shared or reversed edge. 1800 return (b2 == a0) == reverse 1801 } 1802 return OrderedCCW(a0, a2, b2, ab1) 1803 } 1804 1805 // containsNonCrossingBoundary reports whether given two loops whose boundaries 1806 // do not cross (see compareBoundary), if this loop contains the boundary of the 1807 // other loop. If reverse is true, the boundary of the other loop is reversed 1808 // first (which only affects the result when there are shared edges). This method 1809 // is cheaper than compareBoundary because it does not test for edge intersections. 1810 // 1811 // This function requires that neither loop is empty, and that if the other is full, 1812 // then reverse == false. 1813 func (l *Loop) containsNonCrossingBoundary(other *Loop, reverseOther bool) bool { 1814 // The bounds must intersect for containment. 1815 if !l.bound.Intersects(other.bound) { 1816 return false 1817 } 1818 1819 // Full loops are handled as though the loop surrounded the entire sphere. 1820 if l.IsFull() { 1821 return true 1822 } 1823 if other.IsFull() { 1824 return false 1825 } 1826 1827 m, ok := l.findVertex(other.Vertex(0)) 1828 if !ok { 1829 // Since the other loops vertex 0 is not shared, we can check if this contains it. 1830 return l.ContainsPoint(other.Vertex(0)) 1831 } 1832 // Otherwise check whether the edge (b0, b1) is contained by this loop. 1833 return wedgeContainsSemiwedge(l.Vertex(m-1), l.Vertex(m), l.Vertex(m+1), 1834 other.Vertex(1), reverseOther) 1835 } 1836 1837 // TODO(roberts): Differences from the C++ version: 1838 // DistanceToPoint 1839 // DistanceToBoundary 1840 // Project 1841 // ProjectToBoundary 1842 // BoundaryApproxEqual 1843 // BoundaryNear 1844