diff options
Diffstat (limited to 'vendor/github.com/golang/geo/s2/polyline.go')
| -rw-r--r-- | vendor/github.com/golang/geo/s2/polyline.go | 589 |
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 |