mirror of
https://github.com/superseriousbusiness/gotosocial.git
synced 2024-12-05 22:50:42 +01:00
94e87610c4
* add back exif-terminator and use only for jpeg,png,webp * fix arguments passed to terminateExif() * pull in latest exif-terminator * fix test * update processed img --------- Co-authored-by: tobi <tobi.smethurst@protonmail.com>
1817 lines
64 KiB
Go
1817 lines
64 KiB
Go
// Copyright 2015 Google Inc. All rights reserved.
|
|
//
|
|
// Licensed under the Apache License, Version 2.0 (the "License");
|
|
// you may not use this file except in compliance with the License.
|
|
// You may obtain a copy of the License at
|
|
//
|
|
// http://www.apache.org/licenses/LICENSE-2.0
|
|
//
|
|
// Unless required by applicable law or agreed to in writing, software
|
|
// distributed under the License is distributed on an "AS IS" BASIS,
|
|
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|
// See the License for the specific language governing permissions and
|
|
// limitations under the License.
|
|
|
|
package s2
|
|
|
|
import (
|
|
"fmt"
|
|
"io"
|
|
"math"
|
|
|
|
"github.com/golang/geo/r1"
|
|
"github.com/golang/geo/r3"
|
|
"github.com/golang/geo/s1"
|
|
)
|
|
|
|
// Loop represents a simple spherical polygon. It consists of a sequence
|
|
// of vertices where the first vertex is implicitly connected to the
|
|
// last. All loops are defined to have a CCW orientation, i.e. the interior of
|
|
// the loop is on the left side of the edges. This implies that a clockwise
|
|
// loop enclosing a small area is interpreted to be a CCW loop enclosing a
|
|
// very large area.
|
|
//
|
|
// Loops are not allowed to have any duplicate vertices (whether adjacent or
|
|
// not). Non-adjacent edges are not allowed to intersect, and furthermore edges
|
|
// of length 180 degrees are not allowed (i.e., adjacent vertices cannot be
|
|
// antipodal). Loops must have at least 3 vertices (except for the "empty" and
|
|
// "full" loops discussed below).
|
|
//
|
|
// There are two special loops: the "empty" loop contains no points and the
|
|
// "full" loop contains all points. These loops do not have any edges, but to
|
|
// preserve the invariant that every loop can be represented as a vertex
|
|
// chain, they are defined as having exactly one vertex each (see EmptyLoop
|
|
// and FullLoop).
|
|
type Loop struct {
|
|
vertices []Point
|
|
|
|
// originInside keeps a precomputed value whether this loop contains the origin
|
|
// versus computing from the set of vertices every time.
|
|
originInside bool
|
|
|
|
// depth is the nesting depth of this Loop if it is contained by a Polygon
|
|
// or other shape and is used to determine if this loop represents a hole
|
|
// or a filled in portion.
|
|
depth int
|
|
|
|
// bound is a conservative bound on all points contained by this loop.
|
|
// If l.ContainsPoint(P), then l.bound.ContainsPoint(P).
|
|
bound Rect
|
|
|
|
// Since bound is not exact, it is possible that a loop A contains
|
|
// another loop B whose bounds are slightly larger. subregionBound
|
|
// has been expanded sufficiently to account for this error, i.e.
|
|
// if A.Contains(B), then A.subregionBound.Contains(B.bound).
|
|
subregionBound Rect
|
|
|
|
// index is the spatial index for this Loop.
|
|
index *ShapeIndex
|
|
}
|
|
|
|
// LoopFromPoints constructs a loop from the given points.
|
|
func LoopFromPoints(pts []Point) *Loop {
|
|
l := &Loop{
|
|
vertices: pts,
|
|
}
|
|
|
|
l.initOriginAndBound()
|
|
return l
|
|
}
|
|
|
|
// LoopFromCell constructs a loop corresponding to the given cell.
|
|
//
|
|
// Note that the loop and cell *do not* contain exactly the same set of
|
|
// points, because Loop and Cell have slightly different definitions of
|
|
// point containment. For example, a Cell vertex is contained by all
|
|
// four neighboring Cells, but it is contained by exactly one of four
|
|
// Loops constructed from those cells. As another example, the cell
|
|
// coverings of cell and LoopFromCell(cell) will be different, because the
|
|
// loop contains points on its boundary that actually belong to other cells
|
|
// (i.e., the covering will include a layer of neighboring cells).
|
|
func LoopFromCell(c Cell) *Loop {
|
|
l := &Loop{
|
|
vertices: []Point{
|
|
c.Vertex(0),
|
|
c.Vertex(1),
|
|
c.Vertex(2),
|
|
c.Vertex(3),
|
|
},
|
|
}
|
|
|
|
l.initOriginAndBound()
|
|
return l
|
|
}
|
|
|
|
// These two points are used for the special Empty and Full loops.
|
|
var (
|
|
emptyLoopPoint = Point{r3.Vector{X: 0, Y: 0, Z: 1}}
|
|
fullLoopPoint = Point{r3.Vector{X: 0, Y: 0, Z: -1}}
|
|
)
|
|
|
|
// EmptyLoop returns a special "empty" loop.
|
|
func EmptyLoop() *Loop {
|
|
return LoopFromPoints([]Point{emptyLoopPoint})
|
|
}
|
|
|
|
// FullLoop returns a special "full" loop.
|
|
func FullLoop() *Loop {
|
|
return LoopFromPoints([]Point{fullLoopPoint})
|
|
}
|
|
|
|
// initOriginAndBound sets the origin containment for the given point and then calls
|
|
// the initialization for the bounds objects and the internal index.
|
|
func (l *Loop) initOriginAndBound() {
|
|
if len(l.vertices) < 3 {
|
|
// Check for the special "empty" and "full" loops (which have one vertex).
|
|
if !l.isEmptyOrFull() {
|
|
l.originInside = false
|
|
return
|
|
}
|
|
|
|
// This is the special empty or full loop, so the origin depends on if
|
|
// the vertex is in the southern hemisphere or not.
|
|
l.originInside = l.vertices[0].Z < 0
|
|
} else {
|
|
// Point containment testing is done by counting edge crossings starting
|
|
// at a fixed point on the sphere (OriginPoint). We need to know whether
|
|
// the reference point (OriginPoint) is inside or outside the loop before
|
|
// we can construct the ShapeIndex. We do this by first guessing that
|
|
// it is outside, and then seeing whether we get the correct containment
|
|
// result for vertex 1. If the result is incorrect, the origin must be
|
|
// inside the loop.
|
|
//
|
|
// A loop with consecutive vertices A,B,C contains vertex B if and only if
|
|
// the fixed vector R = B.Ortho is contained by the wedge ABC. The
|
|
// wedge is closed at A and open at C, i.e. the point B is inside the loop
|
|
// if A = R but not if C = R. This convention is required for compatibility
|
|
// with VertexCrossing. (Note that we can't use OriginPoint
|
|
// as the fixed vector because of the possibility that B == OriginPoint.)
|
|
l.originInside = false
|
|
v1Inside := OrderedCCW(Point{l.vertices[1].Ortho()}, l.vertices[0], l.vertices[2], l.vertices[1])
|
|
if v1Inside != l.ContainsPoint(l.vertices[1]) {
|
|
l.originInside = true
|
|
}
|
|
}
|
|
|
|
// We *must* call initBound before initializing the index, because
|
|
// initBound calls ContainsPoint which does a bounds check before using
|
|
// the index.
|
|
l.initBound()
|
|
|
|
// Create a new index and add us to it.
|
|
l.index = NewShapeIndex()
|
|
l.index.Add(l)
|
|
}
|
|
|
|
// initBound sets up the approximate bounding Rects for this loop.
|
|
func (l *Loop) initBound() {
|
|
// Check for the special "empty" and "full" loops.
|
|
if l.isEmptyOrFull() {
|
|
if l.IsEmpty() {
|
|
l.bound = EmptyRect()
|
|
} else {
|
|
l.bound = FullRect()
|
|
}
|
|
l.subregionBound = l.bound
|
|
return
|
|
}
|
|
|
|
// The bounding rectangle of a loop is not necessarily the same as the
|
|
// bounding rectangle of its vertices. First, the maximal latitude may be
|
|
// attained along the interior of an edge. Second, the loop may wrap
|
|
// entirely around the sphere (e.g. a loop that defines two revolutions of a
|
|
// candy-cane stripe). Third, the loop may include one or both poles.
|
|
// Note that a small clockwise loop near the equator contains both poles.
|
|
bounder := NewRectBounder()
|
|
for i := 0; i <= len(l.vertices); i++ { // add vertex 0 twice
|
|
bounder.AddPoint(l.Vertex(i))
|
|
}
|
|
b := bounder.RectBound()
|
|
|
|
if l.ContainsPoint(Point{r3.Vector{0, 0, 1}}) {
|
|
b = Rect{r1.Interval{b.Lat.Lo, math.Pi / 2}, s1.FullInterval()}
|
|
}
|
|
// If a loop contains the south pole, then either it wraps entirely
|
|
// around the sphere (full longitude range), or it also contains the
|
|
// north pole in which case b.Lng.IsFull() due to the test above.
|
|
// Either way, we only need to do the south pole containment test if
|
|
// b.Lng.IsFull().
|
|
if b.Lng.IsFull() && l.ContainsPoint(Point{r3.Vector{0, 0, -1}}) {
|
|
b.Lat.Lo = -math.Pi / 2
|
|
}
|
|
l.bound = b
|
|
l.subregionBound = ExpandForSubregions(l.bound)
|
|
}
|
|
|
|
// Validate checks whether this is a valid loop.
|
|
func (l *Loop) Validate() error {
|
|
if err := l.findValidationErrorNoIndex(); err != nil {
|
|
return err
|
|
}
|
|
|
|
// Check for intersections between non-adjacent edges (including at vertices)
|
|
// TODO(roberts): Once shapeutil gets findAnyCrossing uncomment this.
|
|
// return findAnyCrossing(l.index)
|
|
|
|
return nil
|
|
}
|
|
|
|
// findValidationErrorNoIndex reports whether this is not a valid loop, but
|
|
// skips checks that would require a ShapeIndex to be built for the loop. This
|
|
// is primarily used by Polygon to do validation so it doesn't trigger the
|
|
// creation of unneeded ShapeIndices.
|
|
func (l *Loop) findValidationErrorNoIndex() error {
|
|
// All vertices must be unit length.
|
|
for i, v := range l.vertices {
|
|
if !v.IsUnit() {
|
|
return fmt.Errorf("vertex %d is not unit length", i)
|
|
}
|
|
}
|
|
|
|
// Loops must have at least 3 vertices (except for empty and full).
|
|
if len(l.vertices) < 3 {
|
|
if l.isEmptyOrFull() {
|
|
return nil // Skip remaining tests.
|
|
}
|
|
return fmt.Errorf("non-empty, non-full loops must have at least 3 vertices")
|
|
}
|
|
|
|
// Loops are not allowed to have any duplicate vertices or edge crossings.
|
|
// We split this check into two parts. First we check that no edge is
|
|
// degenerate (identical endpoints). Then we check that there are no
|
|
// intersections between non-adjacent edges (including at vertices). The
|
|
// second check needs the ShapeIndex, so it does not fall within the scope
|
|
// of this method.
|
|
for i, v := range l.vertices {
|
|
if v == l.Vertex(i+1) {
|
|
return fmt.Errorf("edge %d is degenerate (duplicate vertex)", i)
|
|
}
|
|
|
|
// Antipodal vertices are not allowed.
|
|
if other := (Point{l.Vertex(i + 1).Mul(-1)}); v == other {
|
|
return fmt.Errorf("vertices %d and %d are antipodal", i,
|
|
(i+1)%len(l.vertices))
|
|
}
|
|
}
|
|
|
|
return nil
|
|
}
|
|
|
|
// Contains reports whether the region contained by this loop is a superset of the
|
|
// region contained by the given other loop.
|
|
func (l *Loop) Contains(o *Loop) bool {
|
|
// For a loop A to contain the loop B, all of the following must
|
|
// be true:
|
|
//
|
|
// (1) There are no edge crossings between A and B except at vertices.
|
|
//
|
|
// (2) At every vertex that is shared between A and B, the local edge
|
|
// ordering implies that A contains B.
|
|
//
|
|
// (3) If there are no shared vertices, then A must contain a vertex of B
|
|
// and B must not contain a vertex of A. (An arbitrary vertex may be
|
|
// chosen in each case.)
|
|
//
|
|
// The second part of (3) is necessary to detect the case of two loops whose
|
|
// union is the entire sphere, i.e. two loops that contains each other's
|
|
// boundaries but not each other's interiors.
|
|
if !l.subregionBound.Contains(o.bound) {
|
|
return false
|
|
}
|
|
|
|
// Special cases to handle either loop being empty or full.
|
|
if l.isEmptyOrFull() || o.isEmptyOrFull() {
|
|
return l.IsFull() || o.IsEmpty()
|
|
}
|
|
|
|
// Check whether there are any edge crossings, and also check the loop
|
|
// relationship at any shared vertices.
|
|
relation := &containsRelation{}
|
|
if hasCrossingRelation(l, o, relation) {
|
|
return false
|
|
}
|
|
|
|
// There are no crossings, and if there are any shared vertices then A
|
|
// contains B locally at each shared vertex.
|
|
if relation.foundSharedVertex {
|
|
return true
|
|
}
|
|
|
|
// Since there are no edge intersections or shared vertices, we just need to
|
|
// test condition (3) above. We can skip this test if we discovered that A
|
|
// contains at least one point of B while checking for edge crossings.
|
|
if !l.ContainsPoint(o.Vertex(0)) {
|
|
return false
|
|
}
|
|
|
|
// We still need to check whether (A union B) is the entire sphere.
|
|
// Normally this check is very cheap due to the bounding box precondition.
|
|
if (o.subregionBound.Contains(l.bound) || o.bound.Union(l.bound).IsFull()) &&
|
|
o.ContainsPoint(l.Vertex(0)) {
|
|
return false
|
|
}
|
|
return true
|
|
}
|
|
|
|
// Intersects reports whether the region contained by this loop intersects the region
|
|
// contained by the other loop.
|
|
func (l *Loop) Intersects(o *Loop) bool {
|
|
// Given two loops, A and B, A.Intersects(B) if and only if !A.Complement().Contains(B).
|
|
//
|
|
// This code is similar to Contains, but is optimized for the case
|
|
// where both loops enclose less than half of the sphere.
|
|
if !l.bound.Intersects(o.bound) {
|
|
return false
|
|
}
|
|
|
|
// Check whether there are any edge crossings, and also check the loop
|
|
// relationship at any shared vertices.
|
|
relation := &intersectsRelation{}
|
|
if hasCrossingRelation(l, o, relation) {
|
|
return true
|
|
}
|
|
if relation.foundSharedVertex {
|
|
return false
|
|
}
|
|
|
|
// Since there are no edge intersections or shared vertices, the loops
|
|
// intersect only if A contains B, B contains A, or the two loops contain
|
|
// each other's boundaries. These checks are usually cheap because of the
|
|
// bounding box preconditions. Note that neither loop is empty (because of
|
|
// the bounding box check above), so it is safe to access vertex(0).
|
|
|
|
// Check whether A contains B, or A and B contain each other's boundaries.
|
|
// (Note that A contains all the vertices of B in either case.)
|
|
if l.subregionBound.Contains(o.bound) || l.bound.Union(o.bound).IsFull() {
|
|
if l.ContainsPoint(o.Vertex(0)) {
|
|
return true
|
|
}
|
|
}
|
|
// Check whether B contains A.
|
|
if o.subregionBound.Contains(l.bound) {
|
|
if o.ContainsPoint(l.Vertex(0)) {
|
|
return true
|
|
}
|
|
}
|
|
return false
|
|
}
|
|
|
|
// Equal reports whether two loops have the same vertices in the same linear order
|
|
// (i.e., cyclic rotations are not allowed).
|
|
func (l *Loop) Equal(other *Loop) bool {
|
|
if len(l.vertices) != len(other.vertices) {
|
|
return false
|
|
}
|
|
|
|
for i, v := range l.vertices {
|
|
if v != other.Vertex(i) {
|
|
return false
|
|
}
|
|
}
|
|
return true
|
|
}
|
|
|
|
// BoundaryEqual reports whether the two loops have the same boundary. This is
|
|
// true if and only if the loops have the same vertices in the same cyclic order
|
|
// (i.e., the vertices may be cyclically rotated). The empty and full loops are
|
|
// considered to have different boundaries.
|
|
func (l *Loop) BoundaryEqual(o *Loop) bool {
|
|
if len(l.vertices) != len(o.vertices) {
|
|
return false
|
|
}
|
|
|
|
// Special case to handle empty or full loops. Since they have the same
|
|
// number of vertices, if one loop is empty/full then so is the other.
|
|
if l.isEmptyOrFull() {
|
|
return l.IsEmpty() == o.IsEmpty()
|
|
}
|
|
|
|
// Loop through the vertices to find the first of ours that matches the
|
|
// starting vertex of the other loop. Use that offset to then 'align' the
|
|
// vertices for comparison.
|
|
for offset, vertex := range l.vertices {
|
|
if vertex == o.Vertex(0) {
|
|
// There is at most one starting offset since loop vertices are unique.
|
|
for i := 0; i < len(l.vertices); i++ {
|
|
if l.Vertex(i+offset) != o.Vertex(i) {
|
|
return false
|
|
}
|
|
}
|
|
return true
|
|
}
|
|
}
|
|
return false
|
|
}
|
|
|
|
// compareBoundary returns +1 if this loop contains the boundary of the other loop,
|
|
// -1 if it excludes the boundary of the other, and 0 if the boundaries of the two
|
|
// loops cross. Shared edges are handled as follows:
|
|
//
|
|
// If XY is a shared edge, define Reversed(XY) to be true if XY
|
|
// appears in opposite directions in both loops.
|
|
// Then this loop contains XY if and only if Reversed(XY) == the other loop is a hole.
|
|
// (Intuitively, this checks whether this loop contains a vanishingly small region
|
|
// extending from the boundary of the other toward the interior of the polygon to
|
|
// which the other belongs.)
|
|
//
|
|
// This function is used for testing containment and intersection of
|
|
// multi-loop polygons. Note that this method is not symmetric, since the
|
|
// result depends on the direction of this loop but not on the direction of
|
|
// the other loop (in the absence of shared edges).
|
|
//
|
|
// This requires that neither loop is empty, and if other loop IsFull, then it must not
|
|
// be a hole.
|
|
func (l *Loop) compareBoundary(o *Loop) int {
|
|
// The bounds must intersect for containment or crossing.
|
|
if !l.bound.Intersects(o.bound) {
|
|
return -1
|
|
}
|
|
|
|
// Full loops are handled as though the loop surrounded the entire sphere.
|
|
if l.IsFull() {
|
|
return 1
|
|
}
|
|
if o.IsFull() {
|
|
return -1
|
|
}
|
|
|
|
// Check whether there are any edge crossings, and also check the loop
|
|
// relationship at any shared vertices.
|
|
relation := newCompareBoundaryRelation(o.IsHole())
|
|
if hasCrossingRelation(l, o, relation) {
|
|
return 0
|
|
}
|
|
if relation.foundSharedVertex {
|
|
if relation.containsEdge {
|
|
return 1
|
|
}
|
|
return -1
|
|
}
|
|
|
|
// There are no edge intersections or shared vertices, so we can check
|
|
// whether A contains an arbitrary vertex of B.
|
|
if l.ContainsPoint(o.Vertex(0)) {
|
|
return 1
|
|
}
|
|
return -1
|
|
}
|
|
|
|
// ContainsOrigin reports true if this loop contains s2.OriginPoint().
|
|
func (l *Loop) ContainsOrigin() bool {
|
|
return l.originInside
|
|
}
|
|
|
|
// ReferencePoint returns the reference point for this loop.
|
|
func (l *Loop) ReferencePoint() ReferencePoint {
|
|
return OriginReferencePoint(l.originInside)
|
|
}
|
|
|
|
// NumEdges returns the number of edges in this shape.
|
|
func (l *Loop) NumEdges() int {
|
|
if l.isEmptyOrFull() {
|
|
return 0
|
|
}
|
|
return len(l.vertices)
|
|
}
|
|
|
|
// Edge returns the endpoints for the given edge index.
|
|
func (l *Loop) Edge(i int) Edge {
|
|
return Edge{l.Vertex(i), l.Vertex(i + 1)}
|
|
}
|
|
|
|
// NumChains reports the number of contiguous edge chains in the Loop.
|
|
func (l *Loop) NumChains() int {
|
|
if l.IsEmpty() {
|
|
return 0
|
|
}
|
|
return 1
|
|
}
|
|
|
|
// Chain returns the i-th edge chain in the Shape.
|
|
func (l *Loop) Chain(chainID int) Chain {
|
|
return Chain{0, l.NumEdges()}
|
|
}
|
|
|
|
// ChainEdge returns the j-th edge of the i-th edge chain.
|
|
func (l *Loop) ChainEdge(chainID, offset int) Edge {
|
|
return Edge{l.Vertex(offset), l.Vertex(offset + 1)}
|
|
}
|
|
|
|
// ChainPosition returns a ChainPosition pair (i, j) such that edgeID is the
|
|
// j-th edge of the Loop.
|
|
func (l *Loop) ChainPosition(edgeID int) ChainPosition {
|
|
return ChainPosition{0, edgeID}
|
|
}
|
|
|
|
// Dimension returns the dimension of the geometry represented by this Loop.
|
|
func (l *Loop) Dimension() int { return 2 }
|
|
|
|
func (l *Loop) typeTag() typeTag { return typeTagNone }
|
|
|
|
func (l *Loop) privateInterface() {}
|
|
|
|
// IsEmpty reports true if this is the special empty loop that contains no points.
|
|
func (l *Loop) IsEmpty() bool {
|
|
return l.isEmptyOrFull() && !l.ContainsOrigin()
|
|
}
|
|
|
|
// IsFull reports true if this is the special full loop that contains all points.
|
|
func (l *Loop) IsFull() bool {
|
|
return l.isEmptyOrFull() && l.ContainsOrigin()
|
|
}
|
|
|
|
// isEmptyOrFull reports true if this loop is either the "empty" or "full" special loops.
|
|
func (l *Loop) isEmptyOrFull() bool {
|
|
return len(l.vertices) == 1
|
|
}
|
|
|
|
// Vertices returns the vertices in the loop.
|
|
func (l *Loop) Vertices() []Point {
|
|
return l.vertices
|
|
}
|
|
|
|
// RectBound returns a tight bounding rectangle. If the loop contains the point,
|
|
// the bound also contains it.
|
|
func (l *Loop) RectBound() Rect {
|
|
return l.bound
|
|
}
|
|
|
|
// CapBound returns a bounding cap that may have more padding than the corresponding
|
|
// RectBound. The bound is conservative such that if the loop contains a point P,
|
|
// the bound also contains it.
|
|
func (l *Loop) CapBound() Cap {
|
|
return l.bound.CapBound()
|
|
}
|
|
|
|
// Vertex returns the vertex for the given index. For convenience, the vertex indices
|
|
// wrap automatically for methods that do index math such as Edge.
|
|
// i.e., Vertex(NumEdges() + n) is the same as Vertex(n).
|
|
func (l *Loop) Vertex(i int) Point {
|
|
return l.vertices[i%len(l.vertices)]
|
|
}
|
|
|
|
// OrientedVertex returns the vertex in reverse order if the loop represents a polygon
|
|
// hole. For example, arguments 0, 1, 2 are mapped to vertices n-1, n-2, n-3, where
|
|
// n == len(vertices). This ensures that the interior of the polygon is always to
|
|
// the left of the vertex chain.
|
|
//
|
|
// This requires: 0 <= i < 2 * len(vertices)
|
|
func (l *Loop) OrientedVertex(i int) Point {
|
|
j := i - len(l.vertices)
|
|
if j < 0 {
|
|
j = i
|
|
}
|
|
if l.IsHole() {
|
|
j = len(l.vertices) - 1 - j
|
|
}
|
|
return l.Vertex(j)
|
|
}
|
|
|
|
// NumVertices returns the number of vertices in this loop.
|
|
func (l *Loop) NumVertices() int {
|
|
return len(l.vertices)
|
|
}
|
|
|
|
// bruteForceContainsPoint reports if the given point is contained by this loop.
|
|
// This method does not use the ShapeIndex, so it is only preferable below a certain
|
|
// size of loop.
|
|
func (l *Loop) bruteForceContainsPoint(p Point) bool {
|
|
origin := OriginPoint()
|
|
inside := l.originInside
|
|
crosser := NewChainEdgeCrosser(origin, p, l.Vertex(0))
|
|
for i := 1; i <= len(l.vertices); i++ { // add vertex 0 twice
|
|
inside = inside != crosser.EdgeOrVertexChainCrossing(l.Vertex(i))
|
|
}
|
|
return inside
|
|
}
|
|
|
|
// ContainsPoint returns true if the loop contains the point.
|
|
func (l *Loop) ContainsPoint(p Point) bool {
|
|
// Empty and full loops don't need a special case, but invalid loops with
|
|
// zero vertices do, so we might as well handle them all at once.
|
|
if len(l.vertices) < 3 {
|
|
return l.originInside
|
|
}
|
|
|
|
// For small loops, and during initial construction, it is faster to just
|
|
// check all the crossing.
|
|
const maxBruteForceVertices = 32
|
|
if len(l.vertices) < maxBruteForceVertices || l.index == nil {
|
|
return l.bruteForceContainsPoint(p)
|
|
}
|
|
|
|
// Otherwise, look up the point in the index.
|
|
it := l.index.Iterator()
|
|
if !it.LocatePoint(p) {
|
|
return false
|
|
}
|
|
return l.iteratorContainsPoint(it, p)
|
|
}
|
|
|
|
// ContainsCell reports whether the given Cell is contained by this Loop.
|
|
func (l *Loop) ContainsCell(target Cell) bool {
|
|
it := l.index.Iterator()
|
|
relation := it.LocateCellID(target.ID())
|
|
|
|
// If "target" is disjoint from all index cells, it is not contained.
|
|
// Similarly, if "target" is subdivided into one or more index cells then it
|
|
// is not contained, since index cells are subdivided only if they (nearly)
|
|
// intersect a sufficient number of edges. (But note that if "target" itself
|
|
// is an index cell then it may be contained, since it could be a cell with
|
|
// no edges in the loop interior.)
|
|
if relation != Indexed {
|
|
return false
|
|
}
|
|
|
|
// Otherwise check if any edges intersect "target".
|
|
if l.boundaryApproxIntersects(it, target) {
|
|
return false
|
|
}
|
|
|
|
// Otherwise check if the loop contains the center of "target".
|
|
return l.iteratorContainsPoint(it, target.Center())
|
|
}
|
|
|
|
// IntersectsCell reports whether this Loop intersects the given cell.
|
|
func (l *Loop) IntersectsCell(target Cell) bool {
|
|
it := l.index.Iterator()
|
|
relation := it.LocateCellID(target.ID())
|
|
|
|
// If target does not overlap any index cell, there is no intersection.
|
|
if relation == Disjoint {
|
|
return false
|
|
}
|
|
// If target is subdivided into one or more index cells, there is an
|
|
// intersection to within the ShapeIndex error bound (see Contains).
|
|
if relation == Subdivided {
|
|
return true
|
|
}
|
|
// If target is an index cell, there is an intersection because index cells
|
|
// are created only if they have at least one edge or they are entirely
|
|
// contained by the loop.
|
|
if it.CellID() == target.id {
|
|
return true
|
|
}
|
|
// Otherwise check if any edges intersect target.
|
|
if l.boundaryApproxIntersects(it, target) {
|
|
return true
|
|
}
|
|
// Otherwise check if the loop contains the center of target.
|
|
return l.iteratorContainsPoint(it, target.Center())
|
|
}
|
|
|
|
// CellUnionBound computes a covering of the Loop.
|
|
func (l *Loop) CellUnionBound() []CellID {
|
|
return l.CapBound().CellUnionBound()
|
|
}
|
|
|
|
// boundaryApproxIntersects reports if the loop's boundary intersects target.
|
|
// It may also return true when the loop boundary does not intersect target but
|
|
// some edge comes within the worst-case error tolerance.
|
|
//
|
|
// This requires that it.Locate(target) returned Indexed.
|
|
func (l *Loop) boundaryApproxIntersects(it *ShapeIndexIterator, target Cell) bool {
|
|
aClipped := it.IndexCell().findByShapeID(0)
|
|
|
|
// If there are no edges, there is no intersection.
|
|
if len(aClipped.edges) == 0 {
|
|
return false
|
|
}
|
|
|
|
// We can save some work if target is the index cell itself.
|
|
if it.CellID() == target.ID() {
|
|
return true
|
|
}
|
|
|
|
// Otherwise check whether any of the edges intersect target.
|
|
maxError := (faceClipErrorUVCoord + intersectsRectErrorUVDist)
|
|
bound := target.BoundUV().ExpandedByMargin(maxError)
|
|
for _, ai := range aClipped.edges {
|
|
v0, v1, ok := ClipToPaddedFace(l.Vertex(ai), l.Vertex(ai+1), target.Face(), maxError)
|
|
if ok && edgeIntersectsRect(v0, v1, bound) {
|
|
return true
|
|
}
|
|
}
|
|
return false
|
|
}
|
|
|
|
// iteratorContainsPoint reports if the iterator that is positioned at the ShapeIndexCell
|
|
// that may contain p, contains the point p.
|
|
func (l *Loop) iteratorContainsPoint(it *ShapeIndexIterator, p Point) bool {
|
|
// Test containment by drawing a line segment from the cell center to the
|
|
// given point and counting edge crossings.
|
|
aClipped := it.IndexCell().findByShapeID(0)
|
|
inside := aClipped.containsCenter
|
|
if len(aClipped.edges) > 0 {
|
|
center := it.Center()
|
|
crosser := NewEdgeCrosser(center, p)
|
|
aiPrev := -2
|
|
for _, ai := range aClipped.edges {
|
|
if ai != aiPrev+1 {
|
|
crosser.RestartAt(l.Vertex(ai))
|
|
}
|
|
aiPrev = ai
|
|
inside = inside != crosser.EdgeOrVertexChainCrossing(l.Vertex(ai+1))
|
|
}
|
|
}
|
|
return inside
|
|
}
|
|
|
|
// RegularLoop creates a loop with the given number of vertices, all
|
|
// located on a circle of the specified radius around the given center.
|
|
func RegularLoop(center Point, radius s1.Angle, numVertices int) *Loop {
|
|
return RegularLoopForFrame(getFrame(center), radius, numVertices)
|
|
}
|
|
|
|
// RegularLoopForFrame creates a loop centered around the z-axis of the given
|
|
// coordinate frame, with the first vertex in the direction of the positive x-axis.
|
|
func RegularLoopForFrame(frame matrix3x3, radius s1.Angle, numVertices int) *Loop {
|
|
return LoopFromPoints(regularPointsForFrame(frame, radius, numVertices))
|
|
}
|
|
|
|
// CanonicalFirstVertex returns a first index and a direction (either +1 or -1)
|
|
// such that the vertex sequence (first, first+dir, ..., first+(n-1)*dir) does
|
|
// not change when the loop vertex order is rotated or inverted. This allows the
|
|
// loop vertices to be traversed in a canonical order. The return values are
|
|
// chosen such that (first, ..., first+n*dir) are in the range [0, 2*n-1] as
|
|
// expected by the Vertex method.
|
|
func (l *Loop) CanonicalFirstVertex() (firstIdx, direction int) {
|
|
firstIdx = 0
|
|
n := len(l.vertices)
|
|
for i := 1; i < n; i++ {
|
|
if l.Vertex(i).Cmp(l.Vertex(firstIdx).Vector) == -1 {
|
|
firstIdx = i
|
|
}
|
|
}
|
|
|
|
// 0 <= firstIdx <= n-1, so (firstIdx+n*dir) <= 2*n-1.
|
|
if l.Vertex(firstIdx+1).Cmp(l.Vertex(firstIdx+n-1).Vector) == -1 {
|
|
return firstIdx, 1
|
|
}
|
|
|
|
// n <= firstIdx <= 2*n-1, so (firstIdx+n*dir) >= 0.
|
|
firstIdx += n
|
|
return firstIdx, -1
|
|
}
|
|
|
|
// TurningAngle returns the sum of the turning angles at each vertex. The return
|
|
// value is positive if the loop is counter-clockwise, negative if the loop is
|
|
// clockwise, and zero if the loop is a great circle. Degenerate and
|
|
// nearly-degenerate loops are handled consistently with Sign. So for example,
|
|
// if a loop has zero area (i.e., it is a very small CCW loop) then the turning
|
|
// angle will always be negative.
|
|
//
|
|
// This quantity is also called the "geodesic curvature" of the loop.
|
|
func (l *Loop) TurningAngle() float64 {
|
|
// For empty and full loops, we return the limit value as the loop area
|
|
// approaches 0 or 4*Pi respectively.
|
|
if l.isEmptyOrFull() {
|
|
if l.ContainsOrigin() {
|
|
return -2 * math.Pi
|
|
}
|
|
return 2 * math.Pi
|
|
}
|
|
|
|
// Don't crash even if the loop is not well-defined.
|
|
if len(l.vertices) < 3 {
|
|
return 0
|
|
}
|
|
|
|
// To ensure that we get the same result when the vertex order is rotated,
|
|
// and that the result is negated when the vertex order is reversed, we need
|
|
// to add up the individual turn angles in a consistent order. (In general,
|
|
// adding up a set of numbers in a different order can change the sum due to
|
|
// rounding errors.)
|
|
//
|
|
// Furthermore, if we just accumulate an ordinary sum then the worst-case
|
|
// error is quadratic in the number of vertices. (This can happen with
|
|
// spiral shapes, where the partial sum of the turning angles can be linear
|
|
// in the number of vertices.) To avoid this we use the Kahan summation
|
|
// algorithm (http://en.wikipedia.org/wiki/Kahan_summation_algorithm).
|
|
n := len(l.vertices)
|
|
i, dir := l.CanonicalFirstVertex()
|
|
sum := TurnAngle(l.Vertex((i+n-dir)%n), l.Vertex(i), l.Vertex((i+dir)%n))
|
|
|
|
compensation := s1.Angle(0)
|
|
for n-1 > 0 {
|
|
i += dir
|
|
angle := TurnAngle(l.Vertex(i-dir), l.Vertex(i), l.Vertex(i+dir))
|
|
oldSum := sum
|
|
angle += compensation
|
|
sum += angle
|
|
compensation = (oldSum - sum) + angle
|
|
n--
|
|
}
|
|
return float64(dir) * float64(sum+compensation)
|
|
}
|
|
|
|
// turningAngleMaxError return the maximum error in TurningAngle. The value is not
|
|
// constant; it depends on the loop.
|
|
func (l *Loop) turningAngleMaxError() float64 {
|
|
// The maximum error can be bounded as follows:
|
|
// 2.24 * dblEpsilon for RobustCrossProd(b, a)
|
|
// 2.24 * dblEpsilon for RobustCrossProd(c, b)
|
|
// 3.25 * dblEpsilon for Angle()
|
|
// 2.00 * dblEpsilon for each addition in the Kahan summation
|
|
// ------------------
|
|
// 9.73 * dblEpsilon
|
|
maxErrorPerVertex := 9.73 * dblEpsilon
|
|
return maxErrorPerVertex * float64(len(l.vertices))
|
|
}
|
|
|
|
// IsHole reports whether this loop represents a hole in its containing polygon.
|
|
func (l *Loop) IsHole() bool { return l.depth&1 != 0 }
|
|
|
|
// Sign returns -1 if this Loop represents a hole in its containing polygon, and +1 otherwise.
|
|
func (l *Loop) Sign() int {
|
|
if l.IsHole() {
|
|
return -1
|
|
}
|
|
return 1
|
|
}
|
|
|
|
// IsNormalized reports whether the loop area is at most 2*pi. Degenerate loops are
|
|
// handled consistently with Sign, i.e., if a loop can be
|
|
// expressed as the union of degenerate or nearly-degenerate CCW triangles,
|
|
// then it will always be considered normalized.
|
|
func (l *Loop) IsNormalized() bool {
|
|
// Optimization: if the longitude span is less than 180 degrees, then the
|
|
// loop covers less than half the sphere and is therefore normalized.
|
|
if l.bound.Lng.Length() < math.Pi {
|
|
return true
|
|
}
|
|
|
|
// We allow some error so that hemispheres are always considered normalized.
|
|
// TODO(roberts): This is no longer required by the Polygon implementation,
|
|
// so alternatively we could create the invariant that a loop is normalized
|
|
// if and only if its complement is not normalized.
|
|
return l.TurningAngle() >= -l.turningAngleMaxError()
|
|
}
|
|
|
|
// Normalize inverts the loop if necessary so that the area enclosed by the loop
|
|
// is at most 2*pi.
|
|
func (l *Loop) Normalize() {
|
|
if !l.IsNormalized() {
|
|
l.Invert()
|
|
}
|
|
}
|
|
|
|
// Invert reverses the order of the loop vertices, effectively complementing the
|
|
// region represented by the loop. For example, the loop ABCD (with edges
|
|
// AB, BC, CD, DA) becomes the loop DCBA (with edges DC, CB, BA, AD).
|
|
// Notice that the last edge is the same in both cases except that its
|
|
// direction has been reversed.
|
|
func (l *Loop) Invert() {
|
|
l.index.Reset()
|
|
if l.isEmptyOrFull() {
|
|
if l.IsFull() {
|
|
l.vertices[0] = emptyLoopPoint
|
|
} else {
|
|
l.vertices[0] = fullLoopPoint
|
|
}
|
|
} else {
|
|
// For non-special loops, reverse the slice of vertices.
|
|
for i := len(l.vertices)/2 - 1; i >= 0; i-- {
|
|
opp := len(l.vertices) - 1 - i
|
|
l.vertices[i], l.vertices[opp] = l.vertices[opp], l.vertices[i]
|
|
}
|
|
}
|
|
|
|
// originInside must be set correctly before building the ShapeIndex.
|
|
l.originInside = !l.originInside
|
|
if l.bound.Lat.Lo > -math.Pi/2 && l.bound.Lat.Hi < math.Pi/2 {
|
|
// The complement of this loop contains both poles.
|
|
l.bound = FullRect()
|
|
l.subregionBound = l.bound
|
|
} else {
|
|
l.initBound()
|
|
}
|
|
l.index.Add(l)
|
|
}
|
|
|
|
// findVertex returns the index of the vertex at the given Point in the range
|
|
// 1..numVertices, and a boolean indicating if a vertex was found.
|
|
func (l *Loop) findVertex(p Point) (index int, ok bool) {
|
|
const notFound = 0
|
|
if len(l.vertices) < 10 {
|
|
// Exhaustive search for loops below a small threshold.
|
|
for i := 1; i <= len(l.vertices); i++ {
|
|
if l.Vertex(i) == p {
|
|
return i, true
|
|
}
|
|
}
|
|
return notFound, false
|
|
}
|
|
|
|
it := l.index.Iterator()
|
|
if !it.LocatePoint(p) {
|
|
return notFound, false
|
|
}
|
|
|
|
aClipped := it.IndexCell().findByShapeID(0)
|
|
for i := aClipped.numEdges() - 1; i >= 0; i-- {
|
|
ai := aClipped.edges[i]
|
|
if l.Vertex(ai) == p {
|
|
if ai == 0 {
|
|
return len(l.vertices), true
|
|
}
|
|
return ai, true
|
|
}
|
|
|
|
if l.Vertex(ai+1) == p {
|
|
return ai + 1, true
|
|
}
|
|
}
|
|
return notFound, false
|
|
}
|
|
|
|
// ContainsNested reports whether the given loops is contained within this loop.
|
|
// This function does not test for edge intersections. The two loops must meet
|
|
// all of the Polygon requirements; for example this implies that their
|
|
// boundaries may not cross or have any shared edges (although they may have
|
|
// shared vertices).
|
|
func (l *Loop) ContainsNested(other *Loop) bool {
|
|
if !l.subregionBound.Contains(other.bound) {
|
|
return false
|
|
}
|
|
|
|
// Special cases to handle either loop being empty or full. Also bail out
|
|
// when B has no vertices to avoid heap overflow on the vertex(1) call
|
|
// below. (This method is called during polygon initialization before the
|
|
// client has an opportunity to call IsValid().)
|
|
if l.isEmptyOrFull() || other.NumVertices() < 2 {
|
|
return l.IsFull() || other.IsEmpty()
|
|
}
|
|
|
|
// We are given that A and B do not share any edges, and that either one
|
|
// loop contains the other or they do not intersect.
|
|
m, ok := l.findVertex(other.Vertex(1))
|
|
if !ok {
|
|
// Since other.vertex(1) is not shared, we can check whether A contains it.
|
|
return l.ContainsPoint(other.Vertex(1))
|
|
}
|
|
|
|
// Check whether the edge order around other.Vertex(1) is compatible with
|
|
// A containing B.
|
|
return WedgeContains(l.Vertex(m-1), l.Vertex(m), l.Vertex(m+1), other.Vertex(0), other.Vertex(2))
|
|
}
|
|
|
|
// surfaceIntegralFloat64 computes the oriented surface integral of some quantity f(x)
|
|
// over the loop interior, given a function f(A,B,C) that returns the
|
|
// corresponding integral over the spherical triangle ABC. Here "oriented
|
|
// surface integral" means:
|
|
//
|
|
// (1) f(A,B,C) must be the integral of f if ABC is counterclockwise,
|
|
// and the integral of -f if ABC is clockwise.
|
|
//
|
|
// (2) The result of this function is *either* the integral of f over the
|
|
// loop interior, or the integral of (-f) over the loop exterior.
|
|
//
|
|
// Note that there are at least two common situations where it easy to work
|
|
// around property (2) above:
|
|
//
|
|
// - If the integral of f over the entire sphere is zero, then it doesn't
|
|
// matter which case is returned because they are always equal.
|
|
//
|
|
// - If f is non-negative, then it is easy to detect when the integral over
|
|
// the loop exterior has been returned, and the integral over the loop
|
|
// interior can be obtained by adding the integral of f over the entire
|
|
// unit sphere (a constant) to the result.
|
|
//
|
|
// Any changes to this method may need corresponding changes to surfaceIntegralPoint as well.
|
|
func (l *Loop) surfaceIntegralFloat64(f func(a, b, c Point) float64) float64 {
|
|
// We sum f over a collection T of oriented triangles, possibly
|
|
// overlapping. Let the sign of a triangle be +1 if it is CCW and -1
|
|
// otherwise, and let the sign of a point x be the sum of the signs of the
|
|
// triangles containing x. Then the collection of triangles T is chosen
|
|
// such that either:
|
|
//
|
|
// (1) Each point in the loop interior has sign +1, and sign 0 otherwise; or
|
|
// (2) Each point in the loop exterior has sign -1, and sign 0 otherwise.
|
|
//
|
|
// The triangles basically consist of a fan from vertex 0 to every loop
|
|
// edge that does not include vertex 0. These triangles will always satisfy
|
|
// either (1) or (2). However, what makes this a bit tricky is that
|
|
// spherical edges become numerically unstable as their length approaches
|
|
// 180 degrees. Of course there is not much we can do if the loop itself
|
|
// contains such edges, but we would like to make sure that all the triangle
|
|
// edges under our control (i.e., the non-loop edges) are stable. For
|
|
// example, consider a loop around the equator consisting of four equally
|
|
// spaced points. This is a well-defined loop, but we cannot just split it
|
|
// into two triangles by connecting vertex 0 to vertex 2.
|
|
//
|
|
// We handle this type of situation by moving the origin of the triangle fan
|
|
// whenever we are about to create an unstable edge. We choose a new
|
|
// location for the origin such that all relevant edges are stable. We also
|
|
// create extra triangles with the appropriate orientation so that the sum
|
|
// of the triangle signs is still correct at every point.
|
|
|
|
// The maximum length of an edge for it to be considered numerically stable.
|
|
// The exact value is fairly arbitrary since it depends on the stability of
|
|
// the function f. The value below is quite conservative but could be
|
|
// reduced further if desired.
|
|
const maxLength = math.Pi - 1e-5
|
|
|
|
var sum float64
|
|
origin := l.Vertex(0)
|
|
for i := 1; i+1 < len(l.vertices); i++ {
|
|
// Let V_i be vertex(i), let O be the current origin, and let length(A,B)
|
|
// be the length of edge (A,B). At the start of each loop iteration, the
|
|
// "leading edge" of the triangle fan is (O,V_i), and we want to extend
|
|
// the triangle fan so that the leading edge is (O,V_i+1).
|
|
//
|
|
// Invariants:
|
|
// 1. length(O,V_i) < maxLength for all (i > 1).
|
|
// 2. Either O == V_0, or O is approximately perpendicular to V_0.
|
|
// 3. "sum" is the oriented integral of f over the area defined by
|
|
// (O, V_0, V_1, ..., V_i).
|
|
if l.Vertex(i+1).Angle(origin.Vector) > maxLength {
|
|
// We are about to create an unstable edge, so choose a new origin O'
|
|
// for the triangle fan.
|
|
oldOrigin := origin
|
|
if origin == l.Vertex(0) {
|
|
// The following point is well-separated from V_i and V_0 (and
|
|
// therefore V_i+1 as well).
|
|
origin = Point{l.Vertex(0).PointCross(l.Vertex(i)).Normalize()}
|
|
} else if l.Vertex(i).Angle(l.Vertex(0).Vector) < maxLength {
|
|
// All edges of the triangle (O, V_0, V_i) are stable, so we can
|
|
// revert to using V_0 as the origin.
|
|
origin = l.Vertex(0)
|
|
} else {
|
|
// (O, V_i+1) and (V_0, V_i) are antipodal pairs, and O and V_0 are
|
|
// perpendicular. Therefore V_0.CrossProd(O) is approximately
|
|
// perpendicular to all of {O, V_0, V_i, V_i+1}, and we can choose
|
|
// this point O' as the new origin.
|
|
origin = Point{l.Vertex(0).Cross(oldOrigin.Vector)}
|
|
|
|
// Advance the edge (V_0,O) to (V_0,O').
|
|
sum += f(l.Vertex(0), oldOrigin, origin)
|
|
}
|
|
// Advance the edge (O,V_i) to (O',V_i).
|
|
sum += f(oldOrigin, l.Vertex(i), origin)
|
|
}
|
|
// Advance the edge (O,V_i) to (O,V_i+1).
|
|
sum += f(origin, l.Vertex(i), l.Vertex(i+1))
|
|
}
|
|
// If the origin is not V_0, we need to sum one more triangle.
|
|
if origin != l.Vertex(0) {
|
|
// Advance the edge (O,V_n-1) to (O,V_0).
|
|
sum += f(origin, l.Vertex(len(l.vertices)-1), l.Vertex(0))
|
|
}
|
|
return sum
|
|
}
|
|
|
|
// surfaceIntegralPoint mirrors the surfaceIntegralFloat64 method but over Points;
|
|
// see that method for commentary. The C++ version uses a templated method.
|
|
// Any changes to this method may need corresponding changes to surfaceIntegralFloat64 as well.
|
|
func (l *Loop) surfaceIntegralPoint(f func(a, b, c Point) Point) Point {
|
|
const maxLength = math.Pi - 1e-5
|
|
var sum r3.Vector
|
|
|
|
origin := l.Vertex(0)
|
|
for i := 1; i+1 < len(l.vertices); i++ {
|
|
if l.Vertex(i+1).Angle(origin.Vector) > maxLength {
|
|
oldOrigin := origin
|
|
if origin == l.Vertex(0) {
|
|
origin = Point{l.Vertex(0).PointCross(l.Vertex(i)).Normalize()}
|
|
} else if l.Vertex(i).Angle(l.Vertex(0).Vector) < maxLength {
|
|
origin = l.Vertex(0)
|
|
} else {
|
|
origin = Point{l.Vertex(0).Cross(oldOrigin.Vector)}
|
|
sum = sum.Add(f(l.Vertex(0), oldOrigin, origin).Vector)
|
|
}
|
|
sum = sum.Add(f(oldOrigin, l.Vertex(i), origin).Vector)
|
|
}
|
|
sum = sum.Add(f(origin, l.Vertex(i), l.Vertex(i+1)).Vector)
|
|
}
|
|
if origin != l.Vertex(0) {
|
|
sum = sum.Add(f(origin, l.Vertex(len(l.vertices)-1), l.Vertex(0)).Vector)
|
|
}
|
|
return Point{sum}
|
|
}
|
|
|
|
// Area returns the area of the loop interior, i.e. the region on the left side of
|
|
// the loop. The return value is between 0 and 4*pi. (Note that the return
|
|
// value is not affected by whether this loop is a "hole" or a "shell".)
|
|
func (l *Loop) Area() float64 {
|
|
// It is surprisingly difficult to compute the area of a loop robustly. The
|
|
// main issues are (1) whether degenerate loops are considered to be CCW or
|
|
// not (i.e., whether their area is close to 0 or 4*pi), and (2) computing
|
|
// the areas of small loops with good relative accuracy.
|
|
//
|
|
// With respect to degeneracies, we would like Area to be consistent
|
|
// with ContainsPoint in that loops that contain many points
|
|
// should have large areas, and loops that contain few points should have
|
|
// small areas. For example, if a degenerate triangle is considered CCW
|
|
// according to s2predicates Sign, then it will contain very few points and
|
|
// its area should be approximately zero. On the other hand if it is
|
|
// considered clockwise, then it will contain virtually all points and so
|
|
// its area should be approximately 4*pi.
|
|
//
|
|
// More precisely, let U be the set of Points for which IsUnitLength
|
|
// is true, let P(U) be the projection of those points onto the mathematical
|
|
// unit sphere, and let V(P(U)) be the Voronoi diagram of the projected
|
|
// points. Then for every loop x, we would like Area to approximately
|
|
// equal the sum of the areas of the Voronoi regions of the points p for
|
|
// which x.ContainsPoint(p) is true.
|
|
//
|
|
// The second issue is that we want to compute the area of small loops
|
|
// accurately. This requires having good relative precision rather than
|
|
// good absolute precision. For example, if the area of a loop is 1e-12 and
|
|
// the error is 1e-15, then the area only has 3 digits of accuracy. (For
|
|
// reference, 1e-12 is about 40 square meters on the surface of the earth.)
|
|
// We would like to have good relative accuracy even for small loops.
|
|
//
|
|
// To achieve these goals, we combine two different methods of computing the
|
|
// area. This first method is based on the Gauss-Bonnet theorem, which says
|
|
// that the area enclosed by the loop equals 2*pi minus the total geodesic
|
|
// curvature of the loop (i.e., the sum of the "turning angles" at all the
|
|
// loop vertices). The big advantage of this method is that as long as we
|
|
// use Sign to compute the turning angle at each vertex, then
|
|
// degeneracies are always handled correctly. In other words, if a
|
|
// degenerate loop is CCW according to the symbolic perturbations used by
|
|
// Sign, then its turning angle will be approximately 2*pi.
|
|
//
|
|
// The disadvantage of the Gauss-Bonnet method is that its absolute error is
|
|
// about 2e-15 times the number of vertices (see turningAngleMaxError).
|
|
// So, it cannot compute the area of small loops accurately.
|
|
//
|
|
// The second method is based on splitting the loop into triangles and
|
|
// summing the area of each triangle. To avoid the difficulty and expense
|
|
// of decomposing the loop into a union of non-overlapping triangles,
|
|
// instead we compute a signed sum over triangles that may overlap (see the
|
|
// comments for surfaceIntegral). The advantage of this method
|
|
// is that the area of each triangle can be computed with much better
|
|
// relative accuracy (using l'Huilier's theorem). The disadvantage is that
|
|
// the result is a signed area: CCW loops may yield a small positive value,
|
|
// while CW loops may yield a small negative value (which is converted to a
|
|
// positive area by adding 4*pi). This means that small errors in computing
|
|
// the signed area may translate into a very large error in the result (if
|
|
// the sign of the sum is incorrect).
|
|
//
|
|
// So, our strategy is to combine these two methods as follows. First we
|
|
// compute the area using the "signed sum over triangles" approach (since it
|
|
// is generally more accurate). We also estimate the maximum error in this
|
|
// result. If the signed area is too close to zero (i.e., zero is within
|
|
// the error bounds), then we double-check the sign of the result using the
|
|
// Gauss-Bonnet method. (In fact we just call IsNormalized, which is
|
|
// based on this method.) If the two methods disagree, we return either 0
|
|
// or 4*pi based on the result of IsNormalized. Otherwise we return the
|
|
// area that we computed originally.
|
|
if l.isEmptyOrFull() {
|
|
if l.ContainsOrigin() {
|
|
return 4 * math.Pi
|
|
}
|
|
return 0
|
|
}
|
|
area := l.surfaceIntegralFloat64(SignedArea)
|
|
|
|
// TODO(roberts): This error estimate is very approximate. There are two
|
|
// issues: (1) SignedArea needs some improvements to ensure that its error
|
|
// is actually never higher than GirardArea, and (2) although the number of
|
|
// triangles in the sum is typically N-2, in theory it could be as high as
|
|
// 2*N for pathological inputs. But in other respects this error bound is
|
|
// very conservative since it assumes that the maximum error is achieved on
|
|
// every triangle.
|
|
maxError := l.turningAngleMaxError()
|
|
|
|
// The signed area should be between approximately -4*pi and 4*pi.
|
|
if area < 0 {
|
|
// We have computed the negative of the area of the loop exterior.
|
|
area += 4 * math.Pi
|
|
}
|
|
|
|
if area > 4*math.Pi {
|
|
area = 4 * math.Pi
|
|
}
|
|
if area < 0 {
|
|
area = 0
|
|
}
|
|
|
|
// If the area is close enough to zero or 4*pi so that the loop orientation
|
|
// is ambiguous, then we compute the loop orientation explicitly.
|
|
if area < maxError && !l.IsNormalized() {
|
|
return 4 * math.Pi
|
|
} else if area > (4*math.Pi-maxError) && l.IsNormalized() {
|
|
return 0
|
|
}
|
|
|
|
return area
|
|
}
|
|
|
|
// Centroid returns the true centroid of the loop multiplied by the area of the
|
|
// loop. The result is not unit length, so you may want to normalize it. Also
|
|
// note that in general, the centroid may not be contained by the loop.
|
|
//
|
|
// We prescale by the loop area for two reasons: (1) it is cheaper to
|
|
// compute this way, and (2) it makes it easier to compute the centroid of
|
|
// more complicated shapes (by splitting them into disjoint regions and
|
|
// adding their centroids).
|
|
//
|
|
// Note that the return value is not affected by whether this loop is a
|
|
// "hole" or a "shell".
|
|
func (l *Loop) Centroid() Point {
|
|
// surfaceIntegralPoint() returns either the integral of position over loop
|
|
// interior, or the negative of the integral of position over the loop
|
|
// exterior. But these two values are the same (!), because the integral of
|
|
// position over the entire sphere is (0, 0, 0).
|
|
return l.surfaceIntegralPoint(TrueCentroid)
|
|
}
|
|
|
|
// Encode encodes the Loop.
|
|
func (l Loop) Encode(w io.Writer) error {
|
|
e := &encoder{w: w}
|
|
l.encode(e)
|
|
return e.err
|
|
}
|
|
|
|
func (l Loop) encode(e *encoder) {
|
|
e.writeInt8(encodingVersion)
|
|
e.writeUint32(uint32(len(l.vertices)))
|
|
for _, v := range l.vertices {
|
|
e.writeFloat64(v.X)
|
|
e.writeFloat64(v.Y)
|
|
e.writeFloat64(v.Z)
|
|
}
|
|
|
|
e.writeBool(l.originInside)
|
|
e.writeInt32(int32(l.depth))
|
|
|
|
// Encode the bound.
|
|
l.bound.encode(e)
|
|
}
|
|
|
|
// Decode decodes a loop.
|
|
func (l *Loop) Decode(r io.Reader) error {
|
|
*l = Loop{}
|
|
d := &decoder{r: asByteReader(r)}
|
|
l.decode(d)
|
|
return d.err
|
|
}
|
|
|
|
func (l *Loop) decode(d *decoder) {
|
|
version := int8(d.readUint8())
|
|
if d.err != nil {
|
|
return
|
|
}
|
|
if version != encodingVersion {
|
|
d.err = fmt.Errorf("cannot decode version %d", version)
|
|
return
|
|
}
|
|
|
|
// Empty loops are explicitly allowed here: a newly created loop has zero vertices
|
|
// and such loops encode and decode properly.
|
|
nvertices := d.readUint32()
|
|
if nvertices > maxEncodedVertices {
|
|
if d.err == nil {
|
|
d.err = fmt.Errorf("too many vertices (%d; max is %d)", nvertices, maxEncodedVertices)
|
|
|
|
}
|
|
return
|
|
}
|
|
l.vertices = make([]Point, nvertices)
|
|
for i := range l.vertices {
|
|
l.vertices[i].X = d.readFloat64()
|
|
l.vertices[i].Y = d.readFloat64()
|
|
l.vertices[i].Z = d.readFloat64()
|
|
}
|
|
l.originInside = d.readBool()
|
|
l.depth = int(d.readUint32())
|
|
l.bound.decode(d)
|
|
l.subregionBound = ExpandForSubregions(l.bound)
|
|
|
|
l.index = NewShapeIndex()
|
|
l.index.Add(l)
|
|
}
|
|
|
|
// Bitmasks to read from properties.
|
|
const (
|
|
originInside = 1 << iota
|
|
boundEncoded
|
|
)
|
|
|
|
func (l *Loop) xyzFaceSiTiVertices() []xyzFaceSiTi {
|
|
ret := make([]xyzFaceSiTi, len(l.vertices))
|
|
for i, v := range l.vertices {
|
|
ret[i].xyz = v
|
|
ret[i].face, ret[i].si, ret[i].ti, ret[i].level = xyzToFaceSiTi(v)
|
|
}
|
|
return ret
|
|
}
|
|
|
|
func (l *Loop) encodeCompressed(e *encoder, snapLevel int, vertices []xyzFaceSiTi) {
|
|
if len(l.vertices) != len(vertices) {
|
|
panic("encodeCompressed: vertices must be the same length as l.vertices")
|
|
}
|
|
if len(vertices) > maxEncodedVertices {
|
|
if e.err == nil {
|
|
e.err = fmt.Errorf("too many vertices (%d; max is %d)", len(vertices), maxEncodedVertices)
|
|
}
|
|
return
|
|
}
|
|
e.writeUvarint(uint64(len(vertices)))
|
|
encodePointsCompressed(e, vertices, snapLevel)
|
|
|
|
props := l.compressedEncodingProperties()
|
|
e.writeUvarint(props)
|
|
e.writeUvarint(uint64(l.depth))
|
|
if props&boundEncoded != 0 {
|
|
l.bound.encode(e)
|
|
}
|
|
}
|
|
|
|
func (l *Loop) compressedEncodingProperties() uint64 {
|
|
var properties uint64
|
|
if l.originInside {
|
|
properties |= originInside
|
|
}
|
|
|
|
// Write whether there is a bound so we can change the threshold later.
|
|
// Recomputing the bound multiplies the decode time taken per vertex
|
|
// by a factor of about 3.5. Without recomputing the bound, decode
|
|
// takes approximately 125 ns / vertex. A loop with 63 vertices
|
|
// encoded without the bound will take ~30us to decode, which is
|
|
// acceptable. At ~3.5 bytes / vertex without the bound, adding
|
|
// the bound will increase the size by <15%, which is also acceptable.
|
|
const minVerticesForBound = 64
|
|
if len(l.vertices) >= minVerticesForBound {
|
|
properties |= boundEncoded
|
|
}
|
|
|
|
return properties
|
|
}
|
|
|
|
func (l *Loop) decodeCompressed(d *decoder, snapLevel int) {
|
|
nvertices := d.readUvarint()
|
|
if d.err != nil {
|
|
return
|
|
}
|
|
if nvertices > maxEncodedVertices {
|
|
d.err = fmt.Errorf("too many vertices (%d; max is %d)", nvertices, maxEncodedVertices)
|
|
return
|
|
}
|
|
l.vertices = make([]Point, nvertices)
|
|
decodePointsCompressed(d, snapLevel, l.vertices)
|
|
properties := d.readUvarint()
|
|
|
|
// Make sure values are valid before using.
|
|
if d.err != nil {
|
|
return
|
|
}
|
|
|
|
l.originInside = (properties & originInside) != 0
|
|
|
|
l.depth = int(d.readUvarint())
|
|
|
|
if (properties & boundEncoded) != 0 {
|
|
l.bound.decode(d)
|
|
if d.err != nil {
|
|
return
|
|
}
|
|
l.subregionBound = ExpandForSubregions(l.bound)
|
|
} else {
|
|
l.initBound()
|
|
}
|
|
|
|
l.index = NewShapeIndex()
|
|
l.index.Add(l)
|
|
}
|
|
|
|
// crossingTarget is an enum representing the possible crossing target cases for relations.
|
|
type crossingTarget int
|
|
|
|
const (
|
|
crossingTargetDontCare crossingTarget = iota
|
|
crossingTargetDontCross
|
|
crossingTargetCross
|
|
)
|
|
|
|
// loopRelation defines the interface for checking a type of relationship between two loops.
|
|
// Some examples of relations are Contains, Intersects, or CompareBoundary.
|
|
type loopRelation interface {
|
|
// Optionally, aCrossingTarget and bCrossingTarget can specify an early-exit
|
|
// condition for the loop relation. If any point P is found such that
|
|
//
|
|
// A.ContainsPoint(P) == aCrossingTarget() &&
|
|
// B.ContainsPoint(P) == bCrossingTarget()
|
|
//
|
|
// then the loop relation is assumed to be the same as if a pair of crossing
|
|
// edges were found. For example, the ContainsPoint relation has
|
|
//
|
|
// aCrossingTarget() == crossingTargetDontCross
|
|
// bCrossingTarget() == crossingTargetCross
|
|
//
|
|
// because if A.ContainsPoint(P) == false and B.ContainsPoint(P) == true
|
|
// for any point P, then it is equivalent to finding an edge crossing (i.e.,
|
|
// since Contains returns false in both cases).
|
|
//
|
|
// Loop relations that do not have an early-exit condition of this form
|
|
// should return crossingTargetDontCare for both crossing targets.
|
|
|
|
// aCrossingTarget reports whether loop A crosses the target point with
|
|
// the given relation type.
|
|
aCrossingTarget() crossingTarget
|
|
// bCrossingTarget reports whether loop B crosses the target point with
|
|
// the given relation type.
|
|
bCrossingTarget() crossingTarget
|
|
|
|
// wedgesCross reports if a shared vertex ab1 and the two associated wedges
|
|
// (a0, ab1, b2) and (b0, ab1, b2) are equivalent to an edge crossing.
|
|
// The loop relation is also allowed to maintain its own internal state, and
|
|
// can return true if it observes any sequence of wedges that are equivalent
|
|
// to an edge crossing.
|
|
wedgesCross(a0, ab1, a2, b0, b2 Point) bool
|
|
}
|
|
|
|
// loopCrosser is a helper type for determining whether two loops cross.
|
|
// It is instantiated twice for each pair of loops to be tested, once for the
|
|
// pair (A,B) and once for the pair (B,A), in order to be able to process
|
|
// edges in either loop nesting order.
|
|
type loopCrosser struct {
|
|
a, b *Loop
|
|
relation loopRelation
|
|
swapped bool
|
|
aCrossingTarget crossingTarget
|
|
bCrossingTarget crossingTarget
|
|
|
|
// state maintained by startEdge and edgeCrossesCell.
|
|
crosser *EdgeCrosser
|
|
aj, bjPrev int
|
|
|
|
// temporary data declared here to avoid repeated memory allocations.
|
|
bQuery *CrossingEdgeQuery
|
|
bCells []*ShapeIndexCell
|
|
}
|
|
|
|
// newLoopCrosser creates a loopCrosser from the given values. If swapped is true,
|
|
// the loops A and B have been swapped. This affects how arguments are passed to
|
|
// the given loop relation, since for example A.Contains(B) is not the same as
|
|
// B.Contains(A).
|
|
func newLoopCrosser(a, b *Loop, relation loopRelation, swapped bool) *loopCrosser {
|
|
l := &loopCrosser{
|
|
a: a,
|
|
b: b,
|
|
relation: relation,
|
|
swapped: swapped,
|
|
aCrossingTarget: relation.aCrossingTarget(),
|
|
bCrossingTarget: relation.bCrossingTarget(),
|
|
bQuery: NewCrossingEdgeQuery(b.index),
|
|
}
|
|
if swapped {
|
|
l.aCrossingTarget, l.bCrossingTarget = l.bCrossingTarget, l.aCrossingTarget
|
|
}
|
|
|
|
return l
|
|
}
|
|
|
|
// startEdge sets the crossers state for checking the given edge of loop A.
|
|
func (l *loopCrosser) startEdge(aj int) {
|
|
l.crosser = NewEdgeCrosser(l.a.Vertex(aj), l.a.Vertex(aj+1))
|
|
l.aj = aj
|
|
l.bjPrev = -2
|
|
}
|
|
|
|
// edgeCrossesCell reports whether the current edge of loop A has any crossings with
|
|
// edges of the index cell of loop B.
|
|
func (l *loopCrosser) edgeCrossesCell(bClipped *clippedShape) bool {
|
|
// Test the current edge of A against all edges of bClipped
|
|
bNumEdges := bClipped.numEdges()
|
|
for j := 0; j < bNumEdges; j++ {
|
|
bj := bClipped.edges[j]
|
|
if bj != l.bjPrev+1 {
|
|
l.crosser.RestartAt(l.b.Vertex(bj))
|
|
}
|
|
l.bjPrev = bj
|
|
if crossing := l.crosser.ChainCrossingSign(l.b.Vertex(bj + 1)); crossing == DoNotCross {
|
|
continue
|
|
} else if crossing == Cross {
|
|
return true
|
|
}
|
|
|
|
// We only need to check each shared vertex once, so we only
|
|
// consider the case where l.aVertex(l.aj+1) == l.b.Vertex(bj+1).
|
|
if l.a.Vertex(l.aj+1) == l.b.Vertex(bj+1) {
|
|
if l.swapped {
|
|
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)) {
|
|
return true
|
|
}
|
|
} else {
|
|
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)) {
|
|
return true
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return false
|
|
}
|
|
|
|
// cellCrossesCell reports whether there are any edge crossings or wedge crossings
|
|
// within the two given cells.
|
|
func (l *loopCrosser) cellCrossesCell(aClipped, bClipped *clippedShape) bool {
|
|
// Test all edges of aClipped against all edges of bClipped.
|
|
for _, edge := range aClipped.edges {
|
|
l.startEdge(edge)
|
|
if l.edgeCrossesCell(bClipped) {
|
|
return true
|
|
}
|
|
}
|
|
|
|
return false
|
|
}
|
|
|
|
// cellCrossesAnySubcell reports whether given an index cell of A, if there are any
|
|
// edge or wedge crossings with any index cell of B contained within bID.
|
|
func (l *loopCrosser) cellCrossesAnySubcell(aClipped *clippedShape, bID CellID) bool {
|
|
// Test all edges of aClipped against all edges of B. The relevant B
|
|
// edges are guaranteed to be children of bID, which lets us find the
|
|
// correct index cells more efficiently.
|
|
bRoot := PaddedCellFromCellID(bID, 0)
|
|
for _, aj := range aClipped.edges {
|
|
// Use an CrossingEdgeQuery starting at bRoot to find the index cells
|
|
// of B that might contain crossing edges.
|
|
l.bCells = l.bQuery.getCells(l.a.Vertex(aj), l.a.Vertex(aj+1), bRoot)
|
|
if len(l.bCells) == 0 {
|
|
continue
|
|
}
|
|
l.startEdge(aj)
|
|
for c := 0; c < len(l.bCells); c++ {
|
|
if l.edgeCrossesCell(l.bCells[c].shapes[0]) {
|
|
return true
|
|
}
|
|
}
|
|
}
|
|
|
|
return false
|
|
}
|
|
|
|
// hasCrossing reports whether given two iterators positioned such that
|
|
// ai.cellID().ContainsCellID(bi.cellID()), there is an edge or wedge crossing
|
|
// anywhere within ai.cellID(). This function advances bi only past ai.cellID().
|
|
func (l *loopCrosser) hasCrossing(ai, bi *rangeIterator) bool {
|
|
// If ai.CellID() intersects many edges of B, then it is faster to use
|
|
// CrossingEdgeQuery to narrow down the candidates. But if it intersects
|
|
// only a few edges, it is faster to check all the crossings directly.
|
|
// We handle this by advancing bi and keeping track of how many edges we
|
|
// would need to test.
|
|
const edgeQueryMinEdges = 20 // Tuned from benchmarks.
|
|
var totalEdges int
|
|
l.bCells = nil
|
|
|
|
for {
|
|
if n := bi.it.IndexCell().shapes[0].numEdges(); n > 0 {
|
|
totalEdges += n
|
|
if totalEdges >= edgeQueryMinEdges {
|
|
// There are too many edges to test them directly, so use CrossingEdgeQuery.
|
|
if l.cellCrossesAnySubcell(ai.it.IndexCell().shapes[0], ai.cellID()) {
|
|
return true
|
|
}
|
|
bi.seekBeyond(ai)
|
|
return false
|
|
}
|
|
l.bCells = append(l.bCells, bi.indexCell())
|
|
}
|
|
bi.next()
|
|
if bi.cellID() > ai.rangeMax {
|
|
break
|
|
}
|
|
}
|
|
|
|
// Test all the edge crossings directly.
|
|
for _, c := range l.bCells {
|
|
if l.cellCrossesCell(ai.it.IndexCell().shapes[0], c.shapes[0]) {
|
|
return true
|
|
}
|
|
}
|
|
|
|
return false
|
|
}
|
|
|
|
// containsCenterMatches reports if the clippedShapes containsCenter boolean corresponds
|
|
// to the crossing target type given. (This is to work around C++ allowing false == 0,
|
|
// true == 1 type implicit conversions and comparisons)
|
|
func containsCenterMatches(a *clippedShape, target crossingTarget) bool {
|
|
return (!a.containsCenter && target == crossingTargetDontCross) ||
|
|
(a.containsCenter && target == crossingTargetCross)
|
|
}
|
|
|
|
// hasCrossingRelation reports whether given two iterators positioned such that
|
|
// ai.cellID().ContainsCellID(bi.cellID()), there is a crossing relationship
|
|
// anywhere within ai.cellID(). Specifically, this method returns true if there
|
|
// is an edge crossing, a wedge crossing, or a point P that matches both relations
|
|
// crossing targets. This function advances both iterators past ai.cellID.
|
|
func (l *loopCrosser) hasCrossingRelation(ai, bi *rangeIterator) bool {
|
|
aClipped := ai.it.IndexCell().shapes[0]
|
|
if aClipped.numEdges() != 0 {
|
|
// The current cell of A has at least one edge, so check for crossings.
|
|
if l.hasCrossing(ai, bi) {
|
|
return true
|
|
}
|
|
ai.next()
|
|
return false
|
|
}
|
|
|
|
if containsCenterMatches(aClipped, l.aCrossingTarget) {
|
|
// The crossing target for A is not satisfied, so we skip over these cells of B.
|
|
bi.seekBeyond(ai)
|
|
ai.next()
|
|
return false
|
|
}
|
|
|
|
// All points within ai.cellID() satisfy the crossing target for A, so it's
|
|
// worth iterating through the cells of B to see whether any cell
|
|
// centers also satisfy the crossing target for B.
|
|
for bi.cellID() <= ai.rangeMax {
|
|
bClipped := bi.it.IndexCell().shapes[0]
|
|
if containsCenterMatches(bClipped, l.bCrossingTarget) {
|
|
return true
|
|
}
|
|
bi.next()
|
|
}
|
|
ai.next()
|
|
return false
|
|
}
|
|
|
|
// hasCrossingRelation checks all edges of loop A for intersection against all edges
|
|
// of loop B and reports if there are any that satisfy the given relation. If there
|
|
// is any shared vertex, the wedges centered at this vertex are sent to the given
|
|
// relation to be tested.
|
|
//
|
|
// If the two loop boundaries cross, this method is guaranteed to return
|
|
// true. It also returns true in certain cases if the loop relationship is
|
|
// equivalent to crossing. For example, if the relation is Contains and a
|
|
// point P is found such that B contains P but A does not contain P, this
|
|
// method will return true to indicate that the result is the same as though
|
|
// a pair of crossing edges were found (since Contains returns false in
|
|
// both cases).
|
|
//
|
|
// See Contains, Intersects and CompareBoundary for the three uses of this function.
|
|
func hasCrossingRelation(a, b *Loop, relation loopRelation) bool {
|
|
// We look for CellID ranges where the indexes of A and B overlap, and
|
|
// then test those edges for crossings.
|
|
ai := newRangeIterator(a.index)
|
|
bi := newRangeIterator(b.index)
|
|
|
|
ab := newLoopCrosser(a, b, relation, false) // Tests edges of A against B
|
|
ba := newLoopCrosser(b, a, relation, true) // Tests edges of B against A
|
|
|
|
for !ai.done() || !bi.done() {
|
|
if ai.rangeMax < bi.rangeMin {
|
|
// The A and B cells don't overlap, and A precedes B.
|
|
ai.seekTo(bi)
|
|
} else if bi.rangeMax < ai.rangeMin {
|
|
// The A and B cells don't overlap, and B precedes A.
|
|
bi.seekTo(ai)
|
|
} else {
|
|
// One cell contains the other. Determine which cell is larger.
|
|
abRelation := int64(ai.it.CellID().lsb() - bi.it.CellID().lsb())
|
|
if abRelation > 0 {
|
|
// A's index cell is larger.
|
|
if ab.hasCrossingRelation(ai, bi) {
|
|
return true
|
|
}
|
|
} else if abRelation < 0 {
|
|
// B's index cell is larger.
|
|
if ba.hasCrossingRelation(bi, ai) {
|
|
return true
|
|
}
|
|
} else {
|
|
// The A and B cells are the same. Since the two cells
|
|
// have the same center point P, check whether P satisfies
|
|
// the crossing targets.
|
|
aClipped := ai.it.IndexCell().shapes[0]
|
|
bClipped := bi.it.IndexCell().shapes[0]
|
|
if containsCenterMatches(aClipped, ab.aCrossingTarget) &&
|
|
containsCenterMatches(bClipped, ab.bCrossingTarget) {
|
|
return true
|
|
}
|
|
// Otherwise test all the edge crossings directly.
|
|
if aClipped.numEdges() > 0 && bClipped.numEdges() > 0 && ab.cellCrossesCell(aClipped, bClipped) {
|
|
return true
|
|
}
|
|
ai.next()
|
|
bi.next()
|
|
}
|
|
}
|
|
}
|
|
return false
|
|
}
|
|
|
|
// containsRelation implements loopRelation for a contains operation. If
|
|
// A.ContainsPoint(P) == false && B.ContainsPoint(P) == true, it is equivalent
|
|
// to having an edge crossing (i.e., Contains returns false).
|
|
type containsRelation struct {
|
|
foundSharedVertex bool
|
|
}
|
|
|
|
func (c *containsRelation) aCrossingTarget() crossingTarget { return crossingTargetDontCross }
|
|
func (c *containsRelation) bCrossingTarget() crossingTarget { return crossingTargetCross }
|
|
func (c *containsRelation) wedgesCross(a0, ab1, a2, b0, b2 Point) bool {
|
|
c.foundSharedVertex = true
|
|
return !WedgeContains(a0, ab1, a2, b0, b2)
|
|
}
|
|
|
|
// intersectsRelation implements loopRelation for an intersects operation. Given
|
|
// two loops, A and B, if A.ContainsPoint(P) == true && B.ContainsPoint(P) == true,
|
|
// it is equivalent to having an edge crossing (i.e., Intersects returns true).
|
|
type intersectsRelation struct {
|
|
foundSharedVertex bool
|
|
}
|
|
|
|
func (i *intersectsRelation) aCrossingTarget() crossingTarget { return crossingTargetCross }
|
|
func (i *intersectsRelation) bCrossingTarget() crossingTarget { return crossingTargetCross }
|
|
func (i *intersectsRelation) wedgesCross(a0, ab1, a2, b0, b2 Point) bool {
|
|
i.foundSharedVertex = true
|
|
return WedgeIntersects(a0, ab1, a2, b0, b2)
|
|
}
|
|
|
|
// compareBoundaryRelation implements loopRelation for comparing boundaries.
|
|
//
|
|
// The compare boundary relation does not have a useful early-exit condition,
|
|
// so we return crossingTargetDontCare for both crossing targets.
|
|
//
|
|
// Aside: A possible early exit condition could be based on the following.
|
|
// If A contains a point of both B and ~B, then A intersects Boundary(B).
|
|
// If ~A contains a point of both B and ~B, then ~A intersects Boundary(B).
|
|
// So if the intersections of {A, ~A} with {B, ~B} are all non-empty,
|
|
// the return value is 0, i.e., Boundary(A) intersects Boundary(B).
|
|
// Unfortunately it isn't worth detecting this situation because by the
|
|
// time we have seen a point in all four intersection regions, we are also
|
|
// guaranteed to have seen at least one pair of crossing edges.
|
|
type compareBoundaryRelation struct {
|
|
reverse bool // True if the other loop should be reversed.
|
|
foundSharedVertex bool // True if any wedge was processed.
|
|
containsEdge bool // True if any edge of the other loop is contained by this loop.
|
|
excludesEdge bool // True if any edge of the other loop is excluded by this loop.
|
|
}
|
|
|
|
func newCompareBoundaryRelation(reverse bool) *compareBoundaryRelation {
|
|
return &compareBoundaryRelation{reverse: reverse}
|
|
}
|
|
|
|
func (c *compareBoundaryRelation) aCrossingTarget() crossingTarget { return crossingTargetDontCare }
|
|
func (c *compareBoundaryRelation) bCrossingTarget() crossingTarget { return crossingTargetDontCare }
|
|
func (c *compareBoundaryRelation) wedgesCross(a0, ab1, a2, b0, b2 Point) bool {
|
|
// Because we don't care about the interior of the other, only its boundary,
|
|
// it is sufficient to check whether this one contains the semiwedge (ab1, b2).
|
|
c.foundSharedVertex = true
|
|
if wedgeContainsSemiwedge(a0, ab1, a2, b2, c.reverse) {
|
|
c.containsEdge = true
|
|
} else {
|
|
c.excludesEdge = true
|
|
}
|
|
return c.containsEdge && c.excludesEdge
|
|
}
|
|
|
|
// wedgeContainsSemiwedge reports whether the wedge (a0, ab1, a2) contains the
|
|
// "semiwedge" defined as any non-empty open set of rays immediately CCW from
|
|
// the edge (ab1, b2). If reverse is true, then substitute clockwise for CCW;
|
|
// this simulates what would happen if the direction of the other loop was reversed.
|
|
func wedgeContainsSemiwedge(a0, ab1, a2, b2 Point, reverse bool) bool {
|
|
if b2 == a0 || b2 == a2 {
|
|
// We have a shared or reversed edge.
|
|
return (b2 == a0) == reverse
|
|
}
|
|
return OrderedCCW(a0, a2, b2, ab1)
|
|
}
|
|
|
|
// containsNonCrossingBoundary reports whether given two loops whose boundaries
|
|
// do not cross (see compareBoundary), if this loop contains the boundary of the
|
|
// other loop. If reverse is true, the boundary of the other loop is reversed
|
|
// first (which only affects the result when there are shared edges). This method
|
|
// is cheaper than compareBoundary because it does not test for edge intersections.
|
|
//
|
|
// This function requires that neither loop is empty, and that if the other is full,
|
|
// then reverse == false.
|
|
func (l *Loop) containsNonCrossingBoundary(other *Loop, reverseOther bool) bool {
|
|
// The bounds must intersect for containment.
|
|
if !l.bound.Intersects(other.bound) {
|
|
return false
|
|
}
|
|
|
|
// Full loops are handled as though the loop surrounded the entire sphere.
|
|
if l.IsFull() {
|
|
return true
|
|
}
|
|
if other.IsFull() {
|
|
return false
|
|
}
|
|
|
|
m, ok := l.findVertex(other.Vertex(0))
|
|
if !ok {
|
|
// Since the other loops vertex 0 is not shared, we can check if this contains it.
|
|
return l.ContainsPoint(other.Vertex(0))
|
|
}
|
|
// Otherwise check whether the edge (b0, b1) is contained by this loop.
|
|
return wedgeContainsSemiwedge(l.Vertex(m-1), l.Vertex(m), l.Vertex(m+1),
|
|
other.Vertex(1), reverseOther)
|
|
}
|
|
|
|
// TODO(roberts): Differences from the C++ version:
|
|
// DistanceToPoint
|
|
// DistanceToBoundary
|
|
// Project
|
|
// ProjectToBoundary
|
|
// BoundaryApproxEqual
|
|
// BoundaryNear
|