1 // Copyright 2018 Google Inc. All rights reserved. 2 // 3 // Licensed under the Apache License, Version 2.0 (the "License"); 4 // you may not use this file except in compliance with the License. 5 // You may obtain a copy of the License at 6 // 7 // http://www.apache.org/licenses/LICENSE-2.0 8 // 9 // Unless required by applicable law or agreed to in writing, software 10 // distributed under the License is distributed on an "AS IS" BASIS, 11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 // See the License for the specific language governing permissions and 13 // limitations under the License. 14 15 package s2 16 17 import ( 18 "github.com/golang/geo/r2" 19 "github.com/golang/geo/s1" 20 ) 21 22 // Tessellation is implemented by subdividing the edge until the estimated 23 // maximum error is below the given tolerance. Estimating error is a hard 24 // problem, especially when the only methods available are point evaluation of 25 // the projection and its inverse. (These are the only methods that 26 // Projection provides, which makes it easier and less error-prone to 27 // implement new projections.) 28 // 29 // One technique that significantly increases robustness is to treat the 30 // geodesic and projected edges as parametric curves rather than geometric ones. 31 // Given a spherical edge AB and a projection p:S2->R2, let f(t) be the 32 // normalized arc length parametrization of AB and let g(t) be the normalized 33 // arc length parameterization of the projected edge p(A)p(B). (In other words, 34 // f(0)=A, f(1)=B, g(0)=p(A), g(1)=p(B).) We now define the geometric error as 35 // the maximum distance from the point p^-1(g(t)) to the geodesic edge AB for 36 // any t in [0,1], where p^-1 denotes the inverse projection. In other words, 37 // the geometric error is the maximum distance from any point on the projected 38 // edge (mapped back onto the sphere) to the geodesic edge AB. On the other 39 // hand we define the parametric error as the maximum distance between the 40 // points f(t) and p^-1(g(t)) for any t in [0,1], i.e. the maximum distance 41 // (measured on the sphere) between the geodesic and projected points at the 42 // same interpolation fraction t. 43 // 44 // The easiest way to estimate the parametric error is to simply evaluate both 45 // edges at their midpoints and measure the distance between them (the "midpoint 46 // method"). This is very fast and works quite well for most edges, however it 47 // has one major drawback: it doesn't handle points of inflection (i.e., points 48 // where the curvature changes sign). For example, edges in the Mercator and 49 // Plate Carree projections always curve towards the equator relative to the 50 // corresponding geodesic edge, so in these projections there is a point of 51 // inflection whenever the projected edge crosses the equator. The worst case 52 // occurs when the edge endpoints have different longitudes but the same 53 // absolute latitude, since in that case the error is non-zero but the edges 54 // have exactly the same midpoint (on the equator). 55 // 56 // One solution to this problem is to split the input edges at all inflection 57 // points (i.e., along the equator in the case of the Mercator and Plate Carree 58 // projections). However for general projections these inflection points can 59 // occur anywhere on the sphere (e.g., consider the Transverse Mercator 60 // projection). This could be addressed by adding methods to the S2Projection 61 // interface to split edges at inflection points but this would make it harder 62 // and more error-prone to implement new projections. 63 // 64 // Another problem with this approach is that the midpoint method sometimes 65 // underestimates the true error even when edges do not cross the equator. 66 // For the Plate Carree and Mercator projections, the midpoint method can 67 // underestimate the error by up to 3%. 68 // 69 // Both of these problems can be solved as follows. We assume that the error 70 // can be modeled as a convex combination of two worst-case functions, one 71 // where the error is maximized at the edge midpoint and another where the 72 // error is *minimized* (i.e., zero) at the edge midpoint. For example, we 73 // could choose these functions as: 74 // 75 // E1(x) = 1 - x^2 76 // E2(x) = x * (1 - x^2) 77 // 78 // where for convenience we use an interpolation parameter "x" in the range 79 // [-1, 1] rather than the original "t" in the range [0, 1]. Note that both 80 // error functions must have roots at x = {-1, 1} since the error must be zero 81 // at the edge endpoints. E1 is simply a parabola whose maximum value is 1 82 // attained at x = 0, while E2 is a cubic with an additional root at x = 0, 83 // and whose maximum value is 2 * sqrt(3) / 9 attained at x = 1 / sqrt(3). 84 // 85 // Next, it is convenient to scale these functions so that the both have a 86 // maximum value of 1. E1 already satisfies this requirement, and we simply 87 // redefine E2 as 88 // 89 // E2(x) = x * (1 - x^2) / (2 * sqrt(3) / 9) 90 // 91 // Now define x0 to be the point where these two functions intersect, i.e. the 92 // point in the range (-1, 1) where E1(x0) = E2(x0). This value has the very 93 // convenient property that if we evaluate the actual error E(x0), then the 94 // maximum error on the entire interval [-1, 1] is bounded by 95 // 96 // E(x) <= E(x0) / E1(x0) 97 // 98 // since whether the error is modeled using E1 or E2, the resulting function 99 // has the same maximum value (namely E(x0) / E1(x0)). If it is modeled as 100 // some other convex combination of E1 and E2, the maximum value can only 101 // decrease. 102 // 103 // Finally, since E2 is not symmetric about the y-axis, we must also allow for 104 // the possibility that the error is a convex combination of E1 and -E2. This 105 // can be handled by evaluating the error at E(-x0) as well, and then 106 // computing the final error bound as 107 // 108 // E(x) <= max(E(x0), E(-x0)) / E1(x0) . 109 // 110 // Effectively, this method is simply evaluating the error at two points about 111 // 1/3 and 2/3 of the way along the edges, and then scaling the maximum of 112 // these two errors by a constant factor. Intuitively, the reason this works 113 // is that if the two edges cross somewhere in the interior, then at least one 114 // of these points will be far from the crossing. 115 // 116 // The actual algorithm implemented below has some additional refinements. 117 // First, edges longer than 90 degrees are always subdivided; this avoids 118 // various unusual situations that can happen with very long edges, and there 119 // is really no reason to avoid adding vertices to edges that are so long. 120 // 121 // Second, the error function E1 above needs to be modified to take into 122 // account spherical distortions. (It turns out that spherical distortions are 123 // beneficial in the case of E2, i.e. they only make its error estimates 124 // slightly more conservative.) To do this, we model E1 as the maximum error 125 // in a Plate Carree edge of length 90 degrees or less. This turns out to be 126 // an edge from 45:-90 to 45:90 (in lat:lng format). The corresponding error 127 // as a function of "x" in the range [-1, 1] can be computed as the distance 128 // between the Plate Caree edge point (45, 90 * x) and the geodesic 129 // edge point (90 - 45 * abs(x), 90 * sgn(x)). Using the Haversine formula, 130 // the corresponding function E1 (normalized to have a maximum value of 1) is: 131 // 132 // E1(x) = 133 // asin(sqrt(sin(Pi / 8 * (1 - x)) ^ 2 + 134 // sin(Pi / 4 * (1 - x)) ^ 2 * cos(Pi / 4) * sin(Pi / 4 * x))) / 135 // asin(sqrt((1 - 1 / sqrt(2)) / 2)) 136 // 137 // Note that this function does not need to be evaluated at runtime, it 138 // simply affects the calculation of the value x0 where E1(x0) = E2(x0) 139 // and the corresponding scaling factor C = 1 / E1(x0). 140 // 141 // ------------------------------------------------------------------ 142 // 143 // In the case of the Mercator and Plate Carree projections this strategy 144 // produces a conservative upper bound (verified using 10 million random 145 // edges). Furthermore the bound is nearly tight; the scaling constant is 146 // C = 1.19289, whereas the maximum observed value was 1.19254. 147 // 148 // Compared to the simpler midpoint evaluation method, this strategy requires 149 // more function evaluations (currently twice as many, but with a smarter 150 // tessellation algorithm it will only be 50% more). It also results in a 151 // small amount of additional tessellation (about 1.5%) compared to the 152 // midpoint method, but this is due almost entirely to the fact that the 153 // midpoint method does not yield conservative error estimates. 154 // 155 // For random edges with a tolerance of 1 meter, the expected amount of 156 // overtessellation is as follows: 157 // 158 // Midpoint Method Cubic Method 159 // Plate Carree 1.8% 3.0% 160 // Mercator 15.8% 17.4% 161 162 const ( 163 // tessellationInterpolationFraction is the fraction at which the two edges 164 // are evaluated in order to measure the error between them. (Edges are 165 // evaluated at two points measured this fraction from either end.) 166 tessellationInterpolationFraction = 0.31215691082248312 167 tessellationScaleFactor = 0.83829992569888509 168 169 // minTessellationTolerance is the minimum supported tolerance (which 170 // corresponds to a distance less than 1 micrometer on the Earth's 171 // surface, but is still much larger than the expected projection and 172 // interpolation errors). 173 minTessellationTolerance s1.Angle = 1e-13 174 ) 175 176 // EdgeTessellator converts an edge in a given projection (e.g., Mercator) into 177 // a chain of spherical geodesic edges such that the maximum distance between 178 // the original edge and the geodesic edge chain is at most the requested 179 // tolerance. Similarly, it can convert a spherical geodesic edge into a chain 180 // of edges in a given 2D projection such that the maximum distance between the 181 // geodesic edge and the chain of projected edges is at most the requested tolerance. 182 // 183 // Method | Input | Output 184 // ------------|------------------------|----------------------- 185 // Projected | S2 geodesics | Planar projected edges 186 // Unprojected | Planar projected edges | S2 geodesics 187 type EdgeTessellator struct { 188 projection Projection 189 190 // The given tolerance scaled by a constant fraction so that it can be 191 // compared against the result returned by estimateMaxError. 192 scaledTolerance s1.ChordAngle 193 } 194 195 // NewEdgeTessellator creates a new edge tessellator for the given projection and tolerance. 196 func NewEdgeTessellator(p Projection, tolerance s1.Angle) *EdgeTessellator { 197 return &EdgeTessellator{ 198 projection: p, 199 scaledTolerance: s1.ChordAngleFromAngle(maxAngle(tolerance, minTessellationTolerance)), 200 } 201 } 202 203 // AppendProjected converts the spherical geodesic edge AB to a chain of planar edges 204 // in the given projection and returns the corresponding vertices. 205 // 206 // If the given projection has one or more coordinate axes that wrap, then 207 // every vertex's coordinates will be as close as possible to the previous 208 // vertex's coordinates. Note that this may yield vertices whose 209 // coordinates are outside the usual range. For example, tessellating the 210 // edge (0:170, 0:-170) (in lat:lng notation) yields (0:170, 0:190). 211 func (e *EdgeTessellator) AppendProjected(a, b Point, vertices []r2.Point) []r2.Point { 212 pa := e.projection.Project(a) 213 if len(vertices) == 0 { 214 vertices = []r2.Point{pa} 215 } else { 216 pa = e.projection.WrapDestination(vertices[len(vertices)-1], pa) 217 } 218 219 pb := e.projection.Project(b) 220 return e.appendProjected(pa, a, pb, b, vertices) 221 } 222 223 // appendProjected splits a geodesic edge AB as necessary and returns the 224 // projected vertices appended to the given vertices. 225 // 226 // The maximum recursion depth is (math.Pi / minTessellationTolerance) < 45 227 func (e *EdgeTessellator) appendProjected(pa r2.Point, a Point, pbIn r2.Point, b Point, vertices []r2.Point) []r2.Point { 228 pb := e.projection.WrapDestination(pa, pbIn) 229 if e.estimateMaxError(pa, a, pb, b) <= e.scaledTolerance { 230 return append(vertices, pb) 231 } 232 233 mid := Point{a.Add(b.Vector).Normalize()} 234 pmid := e.projection.WrapDestination(pa, e.projection.Project(mid)) 235 vertices = e.appendProjected(pa, a, pmid, mid, vertices) 236 return e.appendProjected(pmid, mid, pb, b, vertices) 237 } 238 239 // AppendUnprojected converts the planar edge AB in the given projection to a chain of 240 // spherical geodesic edges and returns the vertices. 241 // 242 // Note that to construct a Loop, you must eliminate the duplicate first and last 243 // vertex. Note also that if the given projection involves coordinate wrapping 244 // (e.g. across the 180 degree meridian) then the first and last vertices may not 245 // be exactly the same. 246 func (e *EdgeTessellator) AppendUnprojected(pa, pb r2.Point, vertices []Point) []Point { 247 a := e.projection.Unproject(pa) 248 b := e.projection.Unproject(pb) 249 250 if len(vertices) == 0 { 251 vertices = []Point{a} 252 } 253 254 // Note that coordinate wrapping can create a small amount of error. For 255 // example in the edge chain "0:-175, 0:179, 0:-177", the first edge is 256 // transformed into "0:-175, 0:-181" while the second is transformed into 257 // "0:179, 0:183". The two coordinate pairs for the middle vertex 258 // ("0:-181" and "0:179") may not yield exactly the same S2Point. 259 return e.appendUnprojected(pa, a, pb, b, vertices) 260 } 261 262 // appendUnprojected interpolates a projected edge and appends the corresponding 263 // points on the sphere. 264 func (e *EdgeTessellator) appendUnprojected(pa r2.Point, a Point, pbIn r2.Point, b Point, vertices []Point) []Point { 265 pb := e.projection.WrapDestination(pa, pbIn) 266 if e.estimateMaxError(pa, a, pb, b) <= e.scaledTolerance { 267 return append(vertices, b) 268 } 269 270 pmid := e.projection.Interpolate(0.5, pa, pb) 271 mid := e.projection.Unproject(pmid) 272 273 vertices = e.appendUnprojected(pa, a, pmid, mid, vertices) 274 return e.appendUnprojected(pmid, mid, pb, b, vertices) 275 } 276 277 func (e *EdgeTessellator) estimateMaxError(pa r2.Point, a Point, pb r2.Point, b Point) s1.ChordAngle { 278 // See the algorithm description at the top of this file. 279 // We always tessellate edges longer than 90 degrees on the sphere, since the 280 // approximation below is not robust enough to handle such edges. 281 if a.Dot(b.Vector) < -1e-14 { 282 return s1.InfChordAngle() 283 } 284 t1 := tessellationInterpolationFraction 285 t2 := 1 - tessellationInterpolationFraction 286 mid1 := Interpolate(t1, a, b) 287 mid2 := Interpolate(t2, a, b) 288 pmid1 := e.projection.Unproject(e.projection.Interpolate(t1, pa, pb)) 289 pmid2 := e.projection.Unproject(e.projection.Interpolate(t2, pa, pb)) 290 return maxChordAngle(ChordAngleBetweenPoints(mid1, pmid1), ChordAngleBetweenPoints(mid2, pmid2)) 291 } 292