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, 0 insertions, 589 deletions
diff --git a/vendor/github.com/golang/geo/s2/polyline.go b/vendor/github.com/golang/geo/s2/polyline.go deleted file mode 100644 index 517968342..000000000 --- a/vendor/github.com/golang/geo/s2/polyline.go +++ /dev/null @@ -1,589 +0,0 @@ -// 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 |