summaryrefslogtreecommitdiff
path: root/vendor/github.com/golang/geo/s2/polyline.go
diff options
context:
space:
mode:
authorLibravatar Tobi Smethurst <31960611+tsmethurst@users.noreply.github.com>2021-08-12 21:03:24 +0200
committerLibravatar GitHub <noreply@github.com>2021-08-12 21:03:24 +0200
commit98263a7de64269898a2f81207e38943b5c8e8653 (patch)
tree743c90f109a6c5d27832d1dcef2388d939f0f77a /vendor/github.com/golang/geo/s2/polyline.go
parentText duplication fix (#137) (diff)
downloadgotosocial-98263a7de64269898a2f81207e38943b5c8e8653.tar.xz
Grand test fixup (#138)
* start fixing up tests * fix up tests + automate with drone * fiddle with linting * messing about with drone.yml * some more fiddling * hmmm * add cache * add vendor directory * verbose * ci updates * update some little things * update sig
Diffstat (limited to 'vendor/github.com/golang/geo/s2/polyline.go')
-rw-r--r--vendor/github.com/golang/geo/s2/polyline.go589
1 files changed, 589 insertions, 0 deletions
diff --git a/vendor/github.com/golang/geo/s2/polyline.go b/vendor/github.com/golang/geo/s2/polyline.go
new file mode 100644
index 000000000..517968342
--- /dev/null
+++ b/vendor/github.com/golang/geo/s2/polyline.go
@@ -0,0 +1,589 @@
+// Copyright 2016 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/s1"
+)
+
+// Polyline represents a sequence of zero or more vertices connected by
+// straight edges (geodesics). Edges of length 0 and 180 degrees are not
+// allowed, i.e. adjacent vertices should not be identical or antipodal.
+type Polyline []Point
+
+// PolylineFromLatLngs creates a new Polyline from the given LatLngs.
+func PolylineFromLatLngs(points []LatLng) *Polyline {
+ p := make(Polyline, len(points))
+ for k, v := range points {
+ p[k] = PointFromLatLng(v)
+ }
+ return &p
+}
+
+// Reverse reverses the order of the Polyline vertices.
+func (p *Polyline) Reverse() {
+ for i := 0; i < len(*p)/2; i++ {
+ (*p)[i], (*p)[len(*p)-i-1] = (*p)[len(*p)-i-1], (*p)[i]
+ }
+}
+
+// Length returns the length of this Polyline.
+func (p *Polyline) Length() s1.Angle {
+ var length s1.Angle
+
+ for i := 1; i < len(*p); i++ {
+ length += (*p)[i-1].Distance((*p)[i])
+ }
+ return length
+}
+
+// Centroid returns the true centroid of the polyline multiplied by the length of the
+// polyline. The result is not unit length, so you may wish to normalize it.
+//
+// Scaling by the Polyline length makes it easy to compute the centroid
+// of several Polylines (by simply adding up their centroids).
+func (p *Polyline) Centroid() Point {
+ var centroid Point
+ for i := 1; i < len(*p); i++ {
+ // The centroid (multiplied by length) is a vector toward the midpoint
+ // of the edge, whose length is twice the sin of half the angle between
+ // the two vertices. Defining theta to be this angle, we have:
+ vSum := (*p)[i-1].Add((*p)[i].Vector) // Length == 2*cos(theta)
+ vDiff := (*p)[i-1].Sub((*p)[i].Vector) // Length == 2*sin(theta)
+
+ // Length == 2*sin(theta)
+ centroid = Point{centroid.Add(vSum.Mul(math.Sqrt(vDiff.Norm2() / vSum.Norm2())))}
+ }
+ return centroid
+}
+
+// Equal reports whether the given Polyline is exactly the same as this one.
+func (p *Polyline) Equal(b *Polyline) bool {
+ if len(*p) != len(*b) {
+ return false
+ }
+ for i, v := range *p {
+ if v != (*b)[i] {
+ return false
+ }
+ }
+
+ return true
+}
+
+// ApproxEqual reports whether two polylines have the same number of vertices,
+// and corresponding vertex pairs are separated by no more the standard margin.
+func (p *Polyline) ApproxEqual(o *Polyline) bool {
+ return p.approxEqual(o, s1.Angle(epsilon))
+}
+
+// approxEqual reports whether two polylines are equal within the given margin.
+func (p *Polyline) approxEqual(o *Polyline, maxError s1.Angle) bool {
+ if len(*p) != len(*o) {
+ return false
+ }
+ for offset, val := range *p {
+ if !val.approxEqual((*o)[offset], maxError) {
+ return false
+ }
+ }
+ return true
+}
+
+// CapBound returns the bounding Cap for this Polyline.
+func (p *Polyline) CapBound() Cap {
+ return p.RectBound().CapBound()
+}
+
+// RectBound returns the bounding Rect for this Polyline.
+func (p *Polyline) RectBound() Rect {
+ rb := NewRectBounder()
+ for _, v := range *p {
+ rb.AddPoint(v)
+ }
+ return rb.RectBound()
+}
+
+// ContainsCell reports whether this Polyline contains the given Cell. Always returns false
+// because "containment" is not numerically well-defined except at the Polyline vertices.
+func (p *Polyline) ContainsCell(cell Cell) bool {
+ return false
+}
+
+// IntersectsCell reports whether this Polyline intersects the given Cell.
+func (p *Polyline) IntersectsCell(cell Cell) bool {
+ if len(*p) == 0 {
+ return false
+ }
+
+ // We only need to check whether the cell contains vertex 0 for correctness,
+ // but these tests are cheap compared to edge crossings so we might as well
+ // check all the vertices.
+ for _, v := range *p {
+ if cell.ContainsPoint(v) {
+ return true
+ }
+ }
+
+ cellVertices := []Point{
+ cell.Vertex(0),
+ cell.Vertex(1),
+ cell.Vertex(2),
+ cell.Vertex(3),
+ }
+
+ for j := 0; j < 4; j++ {
+ crosser := NewChainEdgeCrosser(cellVertices[j], cellVertices[(j+1)&3], (*p)[0])
+ for i := 1; i < len(*p); i++ {
+ if crosser.ChainCrossingSign((*p)[i]) != DoNotCross {
+ // There is a proper crossing, or two vertices were the same.
+ return true
+ }
+ }
+ }
+ return false
+}
+
+// ContainsPoint returns false since Polylines are not closed.
+func (p *Polyline) ContainsPoint(point Point) bool {
+ return false
+}
+
+// CellUnionBound computes a covering of the Polyline.
+func (p *Polyline) CellUnionBound() []CellID {
+ return p.CapBound().CellUnionBound()
+}
+
+// NumEdges returns the number of edges in this shape.
+func (p *Polyline) NumEdges() int {
+ if len(*p) == 0 {
+ return 0
+ }
+ return len(*p) - 1
+}
+
+// Edge returns endpoints for the given edge index.
+func (p *Polyline) Edge(i int) Edge {
+ return Edge{(*p)[i], (*p)[i+1]}
+}
+
+// ReferencePoint returns the default reference point with negative containment because Polylines are not closed.
+func (p *Polyline) ReferencePoint() ReferencePoint {
+ return OriginReferencePoint(false)
+}
+
+// NumChains reports the number of contiguous edge chains in this Polyline.
+func (p *Polyline) NumChains() int {
+ return minInt(1, p.NumEdges())
+}
+
+// Chain returns the i-th edge Chain in the Shape.
+func (p *Polyline) Chain(chainID int) Chain {
+ return Chain{0, p.NumEdges()}
+}
+
+// ChainEdge returns the j-th edge of the i-th edge Chain.
+func (p *Polyline) ChainEdge(chainID, offset int) Edge {
+ return Edge{(*p)[offset], (*p)[offset+1]}
+}
+
+// ChainPosition returns a pair (i, j) such that edgeID is the j-th edge
+func (p *Polyline) ChainPosition(edgeID int) ChainPosition {
+ return ChainPosition{0, edgeID}
+}
+
+// Dimension returns the dimension of the geometry represented by this Polyline.
+func (p *Polyline) Dimension() int { return 1 }
+
+// IsEmpty reports whether this shape contains no points.
+func (p *Polyline) IsEmpty() bool { return defaultShapeIsEmpty(p) }
+
+// IsFull reports whether this shape contains all points on the sphere.
+func (p *Polyline) IsFull() bool { return defaultShapeIsFull(p) }
+
+func (p *Polyline) typeTag() typeTag { return typeTagPolyline }
+
+func (p *Polyline) privateInterface() {}
+
+// findEndVertex reports the maximal end index such that the line segment between
+// the start index and this one such that the line segment between these two
+// vertices passes within the given tolerance of all interior vertices, in order.
+func findEndVertex(p Polyline, tolerance s1.Angle, index int) int {
+ // The basic idea is to keep track of the "pie wedge" of angles
+ // from the starting vertex such that a ray from the starting
+ // vertex at that angle will pass through the discs of radius
+ // tolerance centered around all vertices processed so far.
+ //
+ // First we define a coordinate frame for the tangent and normal
+ // spaces at the starting vertex. Essentially this means picking
+ // three orthonormal vectors X,Y,Z such that X and Y span the
+ // tangent plane at the starting vertex, and Z is up. We use
+ // the coordinate frame to define a mapping from 3D direction
+ // vectors to a one-dimensional ray angle in the range (-π,
+ // π]. The angle of a direction vector is computed by
+ // transforming it into the X,Y,Z basis, and then calculating
+ // atan2(y,x). This mapping allows us to represent a wedge of
+ // angles as a 1D interval. Since the interval wraps around, we
+ // represent it as an Interval, i.e. an interval on the unit
+ // circle.
+ origin := p[index]
+ frame := getFrame(origin)
+
+ // As we go along, we keep track of the current wedge of angles
+ // and the distance to the last vertex (which must be
+ // non-decreasing).
+ currentWedge := s1.FullInterval()
+ var lastDistance s1.Angle
+
+ for index++; index < len(p); index++ {
+ candidate := p[index]
+ distance := origin.Distance(candidate)
+
+ // We don't allow simplification to create edges longer than
+ // 90 degrees, to avoid numeric instability as lengths
+ // approach 180 degrees. We do need to allow for original
+ // edges longer than 90 degrees, though.
+ if distance > math.Pi/2 && lastDistance > 0 {
+ break
+ }
+
+ // Vertices must be in increasing order along the ray, except
+ // for the initial disc around the origin.
+ if distance < lastDistance && lastDistance > tolerance {
+ break
+ }
+
+ lastDistance = distance
+
+ // Points that are within the tolerance distance of the origin
+ // do not constrain the ray direction, so we can ignore them.
+ if distance <= tolerance {
+ continue
+ }
+
+ // If the current wedge of angles does not contain the angle
+ // to this vertex, then stop right now. Note that the wedge
+ // of possible ray angles is not necessarily empty yet, but we
+ // can't continue unless we are willing to backtrack to the
+ // last vertex that was contained within the wedge (since we
+ // don't create new vertices). This would be more complicated
+ // and also make the worst-case running time more than linear.
+ direction := toFrame(frame, candidate)
+ center := math.Atan2(direction.Y, direction.X)
+ if !currentWedge.Contains(center) {
+ break
+ }
+
+ // To determine how this vertex constrains the possible ray
+ // angles, consider the triangle ABC where A is the origin, B
+ // is the candidate vertex, and C is one of the two tangent
+ // points between A and the spherical cap of radius
+ // tolerance centered at B. Then from the spherical law of
+ // sines, sin(a)/sin(A) = sin(c)/sin(C), where a and c are
+ // the lengths of the edges opposite A and C. In our case C
+ // is a 90 degree angle, therefore A = asin(sin(a) / sin(c)).
+ // Angle A is the half-angle of the allowable wedge.
+ halfAngle := math.Asin(math.Sin(tolerance.Radians()) / math.Sin(distance.Radians()))
+ target := s1.IntervalFromPointPair(center, center).Expanded(halfAngle)
+ currentWedge = currentWedge.Intersection(target)
+ }
+
+ // We break out of the loop when we reach a vertex index that
+ // can't be included in the line segment, so back up by one
+ // vertex.
+ return index - 1
+}
+
+// SubsampleVertices returns a subsequence of vertex indices such that the
+// polyline connecting these vertices is never further than the given tolerance from
+// the original polyline. Provided the first and last vertices are distinct,
+// they are always preserved; if they are not, the subsequence may contain
+// only a single index.
+//
+// Some useful properties of the algorithm:
+//
+// - It runs in linear time.
+//
+// - The output always represents a valid polyline. In particular, adjacent
+// output vertices are never identical or antipodal.
+//
+// - The method is not optimal, but it tends to produce 2-3% fewer
+// vertices than the Douglas-Peucker algorithm with the same tolerance.
+//
+// - The output is parametrically equivalent to the original polyline to
+// within the given tolerance. For example, if a polyline backtracks on
+// itself and then proceeds onwards, the backtracking will be preserved
+// (to within the given tolerance). This is different than the
+// Douglas-Peucker algorithm which only guarantees geometric equivalence.
+func (p *Polyline) SubsampleVertices(tolerance s1.Angle) []int {
+ var result []int
+
+ if len(*p) < 1 {
+ return result
+ }
+
+ result = append(result, 0)
+ clampedTolerance := s1.Angle(math.Max(tolerance.Radians(), 0))
+
+ for index := 0; index+1 < len(*p); {
+ nextIndex := findEndVertex(*p, clampedTolerance, index)
+ // Don't create duplicate adjacent vertices.
+ if (*p)[nextIndex] != (*p)[index] {
+ result = append(result, nextIndex)
+ }
+ index = nextIndex
+ }
+
+ return result
+}
+
+// Encode encodes the Polyline.
+func (p Polyline) Encode(w io.Writer) error {
+ e := &encoder{w: w}
+ p.encode(e)
+ return e.err
+}
+
+func (p Polyline) encode(e *encoder) {
+ e.writeInt8(encodingVersion)
+ e.writeUint32(uint32(len(p)))
+ for _, v := range p {
+ e.writeFloat64(v.X)
+ e.writeFloat64(v.Y)
+ e.writeFloat64(v.Z)
+ }
+}
+
+// Decode decodes the polyline.
+func (p *Polyline) Decode(r io.Reader) error {
+ d := decoder{r: asByteReader(r)}
+ p.decode(d)
+ return d.err
+}
+
+func (p *Polyline) decode(d decoder) {
+ version := d.readInt8()
+ if d.err != nil {
+ return
+ }
+ if int(version) != int(encodingVersion) {
+ d.err = fmt.Errorf("can't decode version %d; my version: %d", version, encodingVersion)
+ return
+ }
+ nvertices := d.readUint32()
+ if d.err != nil {
+ return
+ }
+ if nvertices > maxEncodedVertices {
+ d.err = fmt.Errorf("too many vertices (%d; max is %d)", nvertices, maxEncodedVertices)
+ return
+ }
+ *p = make([]Point, nvertices)
+ for i := range *p {
+ (*p)[i].X = d.readFloat64()
+ (*p)[i].Y = d.readFloat64()
+ (*p)[i].Z = d.readFloat64()
+ }
+}
+
+// Project returns a point on the polyline that is closest to the given point,
+// and the index of the next vertex after the projected point. The
+// value of that index is always in the range [1, len(polyline)].
+// The polyline must not be empty.
+func (p *Polyline) Project(point Point) (Point, int) {
+ if len(*p) == 1 {
+ // If there is only one vertex, it is always closest to any given point.
+ return (*p)[0], 1
+ }
+
+ // Initial value larger than any possible distance on the unit sphere.
+ minDist := 10 * s1.Radian
+ minIndex := -1
+
+ // Find the line segment in the polyline that is closest to the point given.
+ for i := 1; i < len(*p); i++ {
+ if dist := DistanceFromSegment(point, (*p)[i-1], (*p)[i]); dist < minDist {
+ minDist = dist
+ minIndex = i
+ }
+ }
+
+ // Compute the point on the segment found that is closest to the point given.
+ closest := Project(point, (*p)[minIndex-1], (*p)[minIndex])
+ if closest == (*p)[minIndex] {
+ minIndex++
+ }
+
+ return closest, minIndex
+}
+
+// IsOnRight reports whether the point given is on the right hand side of the
+// polyline, using a naive definition of "right-hand-sideness" where the point
+// is on the RHS of the polyline iff the point is on the RHS of the line segment
+// in the polyline which it is closest to.
+// The polyline must have at least 2 vertices.
+func (p *Polyline) IsOnRight(point Point) bool {
+ // If the closest point C is an interior vertex of the polyline, let B and D
+ // be the previous and next vertices. The given point P is on the right of
+ // the polyline (locally) if B, P, D are ordered CCW around vertex C.
+ closest, next := p.Project(point)
+ if closest == (*p)[next-1] && next > 1 && next < len(*p) {
+ if point == (*p)[next-1] {
+ // Polyline vertices are not on the RHS.
+ return false
+ }
+ return OrderedCCW((*p)[next-2], point, (*p)[next], (*p)[next-1])
+ }
+ // Otherwise, the closest point C is incident to exactly one polyline edge.
+ // We test the point P against that edge.
+ if next == len(*p) {
+ next--
+ }
+ return Sign(point, (*p)[next], (*p)[next-1])
+}
+
+// Validate checks whether this is a valid polyline or not.
+func (p *Polyline) Validate() error {
+ // All vertices must be unit length.
+ for i, pt := range *p {
+ if !pt.IsUnit() {
+ return fmt.Errorf("vertex %d is not unit length", i)
+ }
+ }
+
+ // Adjacent vertices must not be identical or antipodal.
+ for i := 1; i < len(*p); i++ {
+ prev, cur := (*p)[i-1], (*p)[i]
+ if prev == cur {
+ return fmt.Errorf("vertices %d and %d are identical", i-1, i)
+ }
+ if prev == (Point{cur.Mul(-1)}) {
+ return fmt.Errorf("vertices %d and %d are antipodal", i-1, i)
+ }
+ }
+
+ return nil
+}
+
+// Intersects reports whether this polyline intersects the given polyline. If
+// the polylines share a vertex they are considered to be intersecting. When a
+// polyline endpoint is the only intersection with the other polyline, the
+// function may return true or false arbitrarily.
+//
+// The running time is quadratic in the number of vertices.
+func (p *Polyline) Intersects(o *Polyline) bool {
+ if len(*p) == 0 || len(*o) == 0 {
+ return false
+ }
+
+ if !p.RectBound().Intersects(o.RectBound()) {
+ return false
+ }
+
+ // TODO(roberts): Use ShapeIndex here.
+ for i := 1; i < len(*p); i++ {
+ crosser := NewChainEdgeCrosser((*p)[i-1], (*p)[i], (*o)[0])
+ for j := 1; j < len(*o); j++ {
+ if crosser.ChainCrossingSign((*o)[j]) != DoNotCross {
+ return true
+ }
+ }
+ }
+ return false
+}
+
+// Interpolate returns the point whose distance from vertex 0 along the polyline is
+// the given fraction of the polyline's total length, and the index of
+// the next vertex after the interpolated point P. Fractions less than zero
+// or greater than one are clamped. The return value is unit length. The cost of
+// this function is currently linear in the number of vertices.
+//
+// This method allows the caller to easily construct a given suffix of the
+// polyline by concatenating P with the polyline vertices starting at that next
+// vertex. Note that P is guaranteed to be different than the point at the next
+// vertex, so this will never result in a duplicate vertex.
+//
+// The polyline must not be empty. Note that if fraction >= 1.0, then the next
+// vertex will be set to len(p) (indicating that no vertices from the polyline
+// need to be appended). The value of the next vertex is always between 1 and
+// len(p).
+//
+// This method can also be used to construct a prefix of the polyline, by
+// taking the polyline vertices up to next vertex-1 and appending the
+// returned point P if it is different from the last vertex (since in this
+// case there is no guarantee of distinctness).
+func (p *Polyline) Interpolate(fraction float64) (Point, int) {
+ // We intentionally let the (fraction >= 1) case fall through, since
+ // we need to handle it in the loop below in any case because of
+ // possible roundoff errors.
+ if fraction <= 0 {
+ return (*p)[0], 1
+ }
+ target := s1.Angle(fraction) * p.Length()
+
+ for i := 1; i < len(*p); i++ {
+ length := (*p)[i-1].Distance((*p)[i])
+ if target < length {
+ // This interpolates with respect to arc length rather than
+ // straight-line distance, and produces a unit-length result.
+ result := InterpolateAtDistance(target, (*p)[i-1], (*p)[i])
+
+ // It is possible that (result == vertex(i)) due to rounding errors.
+ if result == (*p)[i] {
+ return result, i + 1
+ }
+ return result, i
+ }
+ target -= length
+ }
+
+ return (*p)[len(*p)-1], len(*p)
+}
+
+// Uninterpolate is the inverse operation of Interpolate. Given a point on the
+// polyline, it returns the ratio of the distance to the point from the
+// beginning of the polyline over the length of the polyline. The return
+// value is always betwen 0 and 1 inclusive.
+//
+// The polyline should not be empty. If it has fewer than 2 vertices, the
+// return value is zero.
+func (p *Polyline) Uninterpolate(point Point, nextVertex int) float64 {
+ if len(*p) < 2 {
+ return 0
+ }
+
+ var sum s1.Angle
+ for i := 1; i < nextVertex; i++ {
+ sum += (*p)[i-1].Distance((*p)[i])
+ }
+ lengthToPoint := sum + (*p)[nextVertex-1].Distance(point)
+ for i := nextVertex; i < len(*p); i++ {
+ sum += (*p)[i-1].Distance((*p)[i])
+ }
+ // The ratio can be greater than 1.0 due to rounding errors or because the
+ // point is not exactly on the polyline.
+ return minFloat64(1.0, float64(lengthToPoint/sum))
+}
+
+// TODO(roberts): Differences from C++.
+// NearlyCoversPolyline
+// InitToSnapped
+// InitToSimplified
+// SnapLevel
+// encode/decode compressed