/
edge_tessellator.go
291 lines (272 loc) · 14 KB
/
edge_tessellator.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
// Copyright 2018 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 (
"github.com/golang/geo/r2"
"github.com/golang/geo/s1"
)
// Tessellation is implemented by subdividing the edge until the estimated
// maximum error is below the given tolerance. Estimating error is a hard
// problem, especially when the only methods available are point evaluation of
// the projection and its inverse. (These are the only methods that
// Projection provides, which makes it easier and less error-prone to
// implement new projections.)
//
// One technique that significantly increases robustness is to treat the
// geodesic and projected edges as parametric curves rather than geometric ones.
// Given a spherical edge AB and a projection p:S2->R2, let f(t) be the
// normalized arc length parametrization of AB and let g(t) be the normalized
// arc length parameterization of the projected edge p(A)p(B). (In other words,
// f(0)=A, f(1)=B, g(0)=p(A), g(1)=p(B).) We now define the geometric error as
// the maximum distance from the point p^-1(g(t)) to the geodesic edge AB for
// any t in [0,1], where p^-1 denotes the inverse projection. In other words,
// the geometric error is the maximum distance from any point on the projected
// edge (mapped back onto the sphere) to the geodesic edge AB. On the other
// hand we define the parametric error as the maximum distance between the
// points f(t) and p^-1(g(t)) for any t in [0,1], i.e. the maximum distance
// (measured on the sphere) between the geodesic and projected points at the
// same interpolation fraction t.
//
// The easiest way to estimate the parametric error is to simply evaluate both
// edges at their midpoints and measure the distance between them (the "midpoint
// method"). This is very fast and works quite well for most edges, however it
// has one major drawback: it doesn't handle points of inflection (i.e., points
// where the curvature changes sign). For example, edges in the Mercator and
// Plate Carree projections always curve towards the equator relative to the
// corresponding geodesic edge, so in these projections there is a point of
// inflection whenever the projected edge crosses the equator. The worst case
// occurs when the edge endpoints have different longitudes but the same
// absolute latitude, since in that case the error is non-zero but the edges
// have exactly the same midpoint (on the equator).
//
// One solution to this problem is to split the input edges at all inflection
// points (i.e., along the equator in the case of the Mercator and Plate Carree
// projections). However for general projections these inflection points can
// occur anywhere on the sphere (e.g., consider the Transverse Mercator
// projection). This could be addressed by adding methods to the S2Projection
// interface to split edges at inflection points but this would make it harder
// and more error-prone to implement new projections.
//
// Another problem with this approach is that the midpoint method sometimes
// underestimates the true error even when edges do not cross the equator.
// For the Plate Carree and Mercator projections, the midpoint method can
// underestimate the error by up to 3%.
//
// Both of these problems can be solved as follows. We assume that the error
// can be modeled as a convex combination of two worst-case functions, one
// where the error is maximized at the edge midpoint and another where the
// error is *minimized* (i.e., zero) at the edge midpoint. For example, we
// could choose these functions as:
//
// E1(x) = 1 - x^2
// E2(x) = x * (1 - x^2)
//
// where for convenience we use an interpolation parameter "x" in the range
// [-1, 1] rather than the original "t" in the range [0, 1]. Note that both
// error functions must have roots at x = {-1, 1} since the error must be zero
// at the edge endpoints. E1 is simply a parabola whose maximum value is 1
// attained at x = 0, while E2 is a cubic with an additional root at x = 0,
// and whose maximum value is 2 * sqrt(3) / 9 attained at x = 1 / sqrt(3).
//
// Next, it is convenient to scale these functions so that the both have a
// maximum value of 1. E1 already satisfies this requirement, and we simply
// redefine E2 as
//
// E2(x) = x * (1 - x^2) / (2 * sqrt(3) / 9)
//
// Now define x0 to be the point where these two functions intersect, i.e. the
// point in the range (-1, 1) where E1(x0) = E2(x0). This value has the very
// convenient property that if we evaluate the actual error E(x0), then the
// maximum error on the entire interval [-1, 1] is bounded by
//
// E(x) <= E(x0) / E1(x0)
//
// since whether the error is modeled using E1 or E2, the resulting function
// has the same maximum value (namely E(x0) / E1(x0)). If it is modeled as
// some other convex combination of E1 and E2, the maximum value can only
// decrease.
//
// Finally, since E2 is not symmetric about the y-axis, we must also allow for
// the possibility that the error is a convex combination of E1 and -E2. This
// can be handled by evaluating the error at E(-x0) as well, and then
// computing the final error bound as
//
// E(x) <= max(E(x0), E(-x0)) / E1(x0) .
//
// Effectively, this method is simply evaluating the error at two points about
// 1/3 and 2/3 of the way along the edges, and then scaling the maximum of
// these two errors by a constant factor. Intuitively, the reason this works
// is that if the two edges cross somewhere in the interior, then at least one
// of these points will be far from the crossing.
//
// The actual algorithm implemented below has some additional refinements.
// First, edges longer than 90 degrees are always subdivided; this avoids
// various unusual situations that can happen with very long edges, and there
// is really no reason to avoid adding vertices to edges that are so long.
//
// Second, the error function E1 above needs to be modified to take into
// account spherical distortions. (It turns out that spherical distortions are
// beneficial in the case of E2, i.e. they only make its error estimates
// slightly more conservative.) To do this, we model E1 as the maximum error
// in a Plate Carree edge of length 90 degrees or less. This turns out to be
// an edge from 45:-90 to 45:90 (in lat:lng format). The corresponding error
// as a function of "x" in the range [-1, 1] can be computed as the distance
// between the Plate Caree edge point (45, 90 * x) and the geodesic
// edge point (90 - 45 * abs(x), 90 * sgn(x)). Using the Haversine formula,
// the corresponding function E1 (normalized to have a maximum value of 1) is:
//
// E1(x) =
// asin(sqrt(sin(Pi / 8 * (1 - x)) ^ 2 +
// sin(Pi / 4 * (1 - x)) ^ 2 * cos(Pi / 4) * sin(Pi / 4 * x))) /
// asin(sqrt((1 - 1 / sqrt(2)) / 2))
//
// Note that this function does not need to be evaluated at runtime, it
// simply affects the calculation of the value x0 where E1(x0) = E2(x0)
// and the corresponding scaling factor C = 1 / E1(x0).
//
// ------------------------------------------------------------------
//
// In the case of the Mercator and Plate Carree projections this strategy
// produces a conservative upper bound (verified using 10 million random
// edges). Furthermore the bound is nearly tight; the scaling constant is
// C = 1.19289, whereas the maximum observed value was 1.19254.
//
// Compared to the simpler midpoint evaluation method, this strategy requires
// more function evaluations (currently twice as many, but with a smarter
// tessellation algorithm it will only be 50% more). It also results in a
// small amount of additional tessellation (about 1.5%) compared to the
// midpoint method, but this is due almost entirely to the fact that the
// midpoint method does not yield conservative error estimates.
//
// For random edges with a tolerance of 1 meter, the expected amount of
// overtessellation is as follows:
//
// Midpoint Method Cubic Method
// Plate Carree 1.8% 3.0%
// Mercator 15.8% 17.4%
const (
// tessellationInterpolationFraction is the fraction at which the two edges
// are evaluated in order to measure the error between them. (Edges are
// evaluated at two points measured this fraction from either end.)
tessellationInterpolationFraction = 0.31215691082248312
tessellationScaleFactor = 0.83829992569888509
// minTessellationTolerance is the minimum supported tolerance (which
// corresponds to a distance less than 1 micrometer on the Earth's
// surface, but is still much larger than the expected projection and
// interpolation errors).
minTessellationTolerance s1.Angle = 1e-13
)
// EdgeTessellator converts an edge in a given projection (e.g., Mercator) into
// a chain of spherical geodesic edges such that the maximum distance between
// the original edge and the geodesic edge chain is at most the requested
// tolerance. Similarly, it can convert a spherical geodesic edge into a chain
// of edges in a given 2D projection such that the maximum distance between the
// geodesic edge and the chain of projected edges is at most the requested tolerance.
//
// Method | Input | Output
// ------------|------------------------|-----------------------
// Projected | S2 geodesics | Planar projected edges
// Unprojected | Planar projected edges | S2 geodesics
type EdgeTessellator struct {
projection Projection
// The given tolerance scaled by a constant fraction so that it can be
// compared against the result returned by estimateMaxError.
scaledTolerance s1.ChordAngle
}
// NewEdgeTessellator creates a new edge tessellator for the given projection and tolerance.
func NewEdgeTessellator(p Projection, tolerance s1.Angle) *EdgeTessellator {
return &EdgeTessellator{
projection: p,
scaledTolerance: s1.ChordAngleFromAngle(maxAngle(tolerance, minTessellationTolerance)),
}
}
// AppendProjected converts the spherical geodesic edge AB to a chain of planar edges
// in the given projection and returns the corresponding vertices.
//
// If the given projection has one or more coordinate axes that wrap, then
// every vertex's coordinates will be as close as possible to the previous
// vertex's coordinates. Note that this may yield vertices whose
// coordinates are outside the usual range. For example, tessellating the
// edge (0:170, 0:-170) (in lat:lng notation) yields (0:170, 0:190).
func (e *EdgeTessellator) AppendProjected(a, b Point, vertices []r2.Point) []r2.Point {
pa := e.projection.Project(a)
if len(vertices) == 0 {
vertices = []r2.Point{pa}
} else {
pa = e.projection.WrapDestination(vertices[len(vertices)-1], pa)
}
pb := e.projection.Project(b)
return e.appendProjected(pa, a, pb, b, vertices)
}
// appendProjected splits a geodesic edge AB as necessary and returns the
// projected vertices appended to the given vertices.
//
// The maximum recursion depth is (math.Pi / minTessellationTolerance) < 45
func (e *EdgeTessellator) appendProjected(pa r2.Point, a Point, pbIn r2.Point, b Point, vertices []r2.Point) []r2.Point {
pb := e.projection.WrapDestination(pa, pbIn)
if e.estimateMaxError(pa, a, pb, b) <= e.scaledTolerance {
return append(vertices, pb)
}
mid := Point{a.Add(b.Vector).Normalize()}
pmid := e.projection.WrapDestination(pa, e.projection.Project(mid))
vertices = e.appendProjected(pa, a, pmid, mid, vertices)
return e.appendProjected(pmid, mid, pb, b, vertices)
}
// AppendUnprojected converts the planar edge AB in the given projection to a chain of
// spherical geodesic edges and returns the vertices.
//
// Note that to construct a Loop, you must eliminate the duplicate first and last
// vertex. Note also that if the given projection involves coordinate wrapping
// (e.g. across the 180 degree meridian) then the first and last vertices may not
// be exactly the same.
func (e *EdgeTessellator) AppendUnprojected(pa, pb r2.Point, vertices []Point) []Point {
a := e.projection.Unproject(pa)
b := e.projection.Unproject(pb)
if len(vertices) == 0 {
vertices = []Point{a}
}
// Note that coordinate wrapping can create a small amount of error. For
// example in the edge chain "0:-175, 0:179, 0:-177", the first edge is
// transformed into "0:-175, 0:-181" while the second is transformed into
// "0:179, 0:183". The two coordinate pairs for the middle vertex
// ("0:-181" and "0:179") may not yield exactly the same S2Point.
return e.appendUnprojected(pa, a, pb, b, vertices)
}
// appendUnprojected interpolates a projected edge and appends the corresponding
// points on the sphere.
func (e *EdgeTessellator) appendUnprojected(pa r2.Point, a Point, pbIn r2.Point, b Point, vertices []Point) []Point {
pb := e.projection.WrapDestination(pa, pbIn)
if e.estimateMaxError(pa, a, pb, b) <= e.scaledTolerance {
return append(vertices, b)
}
pmid := e.projection.Interpolate(0.5, pa, pb)
mid := e.projection.Unproject(pmid)
vertices = e.appendUnprojected(pa, a, pmid, mid, vertices)
return e.appendUnprojected(pmid, mid, pb, b, vertices)
}
func (e *EdgeTessellator) estimateMaxError(pa r2.Point, a Point, pb r2.Point, b Point) s1.ChordAngle {
// See the algorithm description at the top of this file.
// We always tessellate edges longer than 90 degrees on the sphere, since the
// approximation below is not robust enough to handle such edges.
if a.Dot(b.Vector) < -1e-14 {
return s1.InfChordAngle()
}
t1 := tessellationInterpolationFraction
t2 := 1 - tessellationInterpolationFraction
mid1 := Interpolate(t1, a, b)
mid2 := Interpolate(t2, a, b)
pmid1 := e.projection.Unproject(e.projection.Interpolate(t1, pa, pb))
pmid2 := e.projection.Unproject(e.projection.Interpolate(t2, pa, pb))
return maxChordAngle(ChordAngleBetweenPoints(mid1, pmid1), ChordAngleBetweenPoints(mid2, pmid2))
}