/
polyline.go
392 lines (339 loc) · 12 KB
/
polyline.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
/*
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
}
// Equals reports whether the given Polyline is exactly the same as this one.
func (p *Polyline) Equals(b *Polyline) bool {
if len(*p) != len(*b) {
return false
}
for i, v := range *p {
if v != (*b)[i] {
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]}
}
// HasInterior returns false as Polylines are not closed.
func (p *Polyline) HasInterior() bool {
return false
}
// 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 min(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() dimension { return polylineGeometry }
// 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()
}
}
// TODO(roberts): Differences from C++.
// IsValid
// Suffix
// Interpolate/UnInterpolate
// Project
// IsPointOnRight
// Intersects(Polyline)
// Reverse
// ApproxEqual
// NearlyCoversPolyline