1 // Copyright 2017 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 "math" 19 20 "github.com/golang/geo/r1" 21 "github.com/golang/geo/r3" 22 "github.com/golang/geo/s1" 23 ) 24 25 // RectBounder is used to compute a bounding rectangle that contains all edges 26 // defined by a vertex chain (v0, v1, v2, ...). All vertices must be unit length. 27 // Note that the bounding rectangle of an edge can be larger than the bounding 28 // rectangle of its endpoints, e.g. consider an edge that passes through the North Pole. 29 // 30 // The bounds are calculated conservatively to account for numerical errors 31 // when points are converted to LatLngs. More precisely, this function 32 // guarantees the following: 33 // Let L be a closed edge chain (Loop) such that the interior of the loop does 34 // not contain either pole. Now if P is any point such that L.ContainsPoint(P), 35 // then RectBound(L).ContainsPoint(LatLngFromPoint(P)). 36 type RectBounder struct { 37 // The previous vertex in the chain. 38 a Point 39 // The previous vertex latitude longitude. 40 aLL LatLng 41 bound Rect 42 } 43 44 // NewRectBounder returns a new instance of a RectBounder. 45 func NewRectBounder() *RectBounder { 46 return &RectBounder{ 47 bound: EmptyRect(), 48 } 49 } 50 51 // maxErrorForTests returns the maximum error in RectBound provided that the 52 // result does not include either pole. It is only used for testing purposes 53 func (r *RectBounder) maxErrorForTests() LatLng { 54 // The maximum error in the latitude calculation is 55 // 3.84 * dblEpsilon for the PointCross calculation 56 // 0.96 * dblEpsilon for the Latitude calculation 57 // 5 * dblEpsilon added by AddPoint/RectBound to compensate for error 58 // ----------------- 59 // 9.80 * dblEpsilon maximum error in result 60 // 61 // The maximum error in the longitude calculation is dblEpsilon. RectBound 62 // does not do any expansion because this isn't necessary in order to 63 // bound the *rounded* longitudes of contained points. 64 return LatLng{10 * dblEpsilon * s1.Radian, 1 * dblEpsilon * s1.Radian} 65 } 66 67 // AddPoint adds the given point to the chain. The Point must be unit length. 68 func (r *RectBounder) AddPoint(b Point) { 69 bLL := LatLngFromPoint(b) 70 71 if r.bound.IsEmpty() { 72 r.a = b 73 r.aLL = bLL 74 r.bound = r.bound.AddPoint(bLL) 75 return 76 } 77 78 // First compute the cross product N = A x B robustly. This is the normal 79 // to the great circle through A and B. We don't use RobustSign 80 // since that method returns an arbitrary vector orthogonal to A if the two 81 // vectors are proportional, and we want the zero vector in that case. 82 n := r.a.Sub(b.Vector).Cross(r.a.Add(b.Vector)) // N = 2 * (A x B) 83 84 // The relative error in N gets large as its norm gets very small (i.e., 85 // when the two points are nearly identical or antipodal). We handle this 86 // by choosing a maximum allowable error, and if the error is greater than 87 // this we fall back to a different technique. Since it turns out that 88 // the other sources of error in converting the normal to a maximum 89 // latitude add up to at most 1.16 * dblEpsilon, and it is desirable to 90 // have the total error be a multiple of dblEpsilon, we have chosen to 91 // limit the maximum error in the normal to be 3.84 * dblEpsilon. 92 // It is possible to show that the error is less than this when 93 // 94 // n.Norm() >= 8 * sqrt(3) / (3.84 - 0.5 - sqrt(3)) * dblEpsilon 95 // = 1.91346e-15 (about 8.618 * dblEpsilon) 96 nNorm := n.Norm() 97 if nNorm < 1.91346e-15 { 98 // A and B are either nearly identical or nearly antipodal (to within 99 // 4.309 * dblEpsilon, or about 6 nanometers on the earth's surface). 100 if r.a.Dot(b.Vector) < 0 { 101 // The two points are nearly antipodal. The easiest solution is to 102 // assume that the edge between A and B could go in any direction 103 // around the sphere. 104 r.bound = FullRect() 105 } else { 106 // The two points are nearly identical (to within 4.309 * dblEpsilon). 107 // In this case we can just use the bounding rectangle of the points, 108 // since after the expansion done by GetBound this Rect is 109 // guaranteed to include the (lat,lng) values of all points along AB. 110 r.bound = r.bound.Union(RectFromLatLng(r.aLL).AddPoint(bLL)) 111 } 112 r.a = b 113 r.aLL = bLL 114 return 115 } 116 117 // Compute the longitude range spanned by AB. 118 lngAB := s1.EmptyInterval().AddPoint(r.aLL.Lng.Radians()).AddPoint(bLL.Lng.Radians()) 119 if lngAB.Length() >= math.Pi-2*dblEpsilon { 120 // The points lie on nearly opposite lines of longitude to within the 121 // maximum error of the calculation. The easiest solution is to assume 122 // that AB could go on either side of the pole. 123 lngAB = s1.FullInterval() 124 } 125 126 // Next we compute the latitude range spanned by the edge AB. We start 127 // with the range spanning the two endpoints of the edge: 128 latAB := r1.IntervalFromPoint(r.aLL.Lat.Radians()).AddPoint(bLL.Lat.Radians()) 129 130 // This is the desired range unless the edge AB crosses the plane 131 // through N and the Z-axis (which is where the great circle through A 132 // and B attains its minimum and maximum latitudes). To test whether AB 133 // crosses this plane, we compute a vector M perpendicular to this 134 // plane and then project A and B onto it. 135 m := n.Cross(r3.Vector{0, 0, 1}) 136 mA := m.Dot(r.a.Vector) 137 mB := m.Dot(b.Vector) 138 139 // We want to test the signs of "mA" and "mB", so we need to bound 140 // the error in these calculations. It is possible to show that the 141 // total error is bounded by 142 // 143 // (1 + sqrt(3)) * dblEpsilon * nNorm + 8 * sqrt(3) * (dblEpsilon**2) 144 // = 6.06638e-16 * nNorm + 6.83174e-31 145 146 mError := 6.06638e-16*nNorm + 6.83174e-31 147 if mA*mB < 0 || math.Abs(mA) <= mError || math.Abs(mB) <= mError { 148 // Minimum/maximum latitude *may* occur in the edge interior. 149 // 150 // The maximum latitude is 90 degrees minus the latitude of N. We 151 // compute this directly using atan2 in order to get maximum accuracy 152 // near the poles. 153 // 154 // Our goal is compute a bound that contains the computed latitudes of 155 // all S2Points P that pass the point-in-polygon containment test. 156 // There are three sources of error we need to consider: 157 // - the directional error in N (at most 3.84 * dblEpsilon) 158 // - converting N to a maximum latitude 159 // - computing the latitude of the test point P 160 // The latter two sources of error are at most 0.955 * dblEpsilon 161 // individually, but it is possible to show by a more complex analysis 162 // that together they can add up to at most 1.16 * dblEpsilon, for a 163 // total error of 5 * dblEpsilon. 164 // 165 // We add 3 * dblEpsilon to the bound here, and GetBound() will pad 166 // the bound by another 2 * dblEpsilon. 167 maxLat := math.Min( 168 math.Atan2(math.Sqrt(n.X*n.X+n.Y*n.Y), math.Abs(n.Z))+3*dblEpsilon, 169 math.Pi/2) 170 171 // In order to get tight bounds when the two points are close together, 172 // we also bound the min/max latitude relative to the latitudes of the 173 // endpoints A and B. First we compute the distance between A and B, 174 // and then we compute the maximum change in latitude between any two 175 // points along the great circle that are separated by this distance. 176 // This gives us a latitude change "budget". Some of this budget must 177 // be spent getting from A to B; the remainder bounds the round-trip 178 // distance (in latitude) from A or B to the min or max latitude 179 // attained along the edge AB. 180 latBudget := 2 * math.Asin(0.5*(r.a.Sub(b.Vector)).Norm()*math.Sin(maxLat)) 181 maxDelta := 0.5*(latBudget-latAB.Length()) + dblEpsilon 182 183 // Test whether AB passes through the point of maximum latitude or 184 // minimum latitude. If the dot product(s) are small enough then the 185 // result may be ambiguous. 186 if mA <= mError && mB >= -mError { 187 latAB.Hi = math.Min(maxLat, latAB.Hi+maxDelta) 188 } 189 if mB <= mError && mA >= -mError { 190 latAB.Lo = math.Max(-maxLat, latAB.Lo-maxDelta) 191 } 192 } 193 r.a = b 194 r.aLL = bLL 195 r.bound = r.bound.Union(Rect{latAB, lngAB}) 196 } 197 198 // RectBound returns the bounding rectangle of the edge chain that connects the 199 // vertices defined so far. This bound satisfies the guarantee made 200 // above, i.e. if the edge chain defines a Loop, then the bound contains 201 // the LatLng coordinates of all Points contained by the loop. 202 func (r *RectBounder) RectBound() Rect { 203 return r.bound.expanded(LatLng{s1.Angle(2 * dblEpsilon), 0}).PolarClosure() 204 } 205 206 // ExpandForSubregions expands a bounding Rect so that it is guaranteed to 207 // contain the bounds of any subregion whose bounds are computed using 208 // ComputeRectBound. For example, consider a loop L that defines a square. 209 // GetBound ensures that if a point P is contained by this square, then 210 // LatLngFromPoint(P) is contained by the bound. But now consider a diamond 211 // shaped loop S contained by L. It is possible that GetBound returns a 212 // *larger* bound for S than it does for L, due to rounding errors. This 213 // method expands the bound for L so that it is guaranteed to contain the 214 // bounds of any subregion S. 215 // 216 // More precisely, if L is a loop that does not contain either pole, and S 217 // is a loop such that L.Contains(S), then 218 // 219 // ExpandForSubregions(L.RectBound).Contains(S.RectBound). 220 func ExpandForSubregions(bound Rect) Rect { 221 // Empty bounds don't need expansion. 222 if bound.IsEmpty() { 223 return bound 224 } 225 226 // First we need to check whether the bound B contains any nearly-antipodal 227 // points (to within 4.309 * dblEpsilon). If so then we need to return 228 // FullRect, since the subregion might have an edge between two 229 // such points, and AddPoint returns Full for such edges. Note that 230 // this can happen even if B is not Full for example, consider a loop 231 // that defines a 10km strip straddling the equator extending from 232 // longitudes -100 to +100 degrees. 233 // 234 // It is easy to check whether B contains any antipodal points, but checking 235 // for nearly-antipodal points is trickier. Essentially we consider the 236 // original bound B and its reflection through the origin B', and then test 237 // whether the minimum distance between B and B' is less than 4.309 * dblEpsilon. 238 239 // lngGap is a lower bound on the longitudinal distance between B and its 240 // reflection B'. (2.5 * dblEpsilon is the maximum combined error of the 241 // endpoint longitude calculations and the Length call.) 242 lngGap := math.Max(0, math.Pi-bound.Lng.Length()-2.5*dblEpsilon) 243 244 // minAbsLat is the minimum distance from B to the equator (if zero or 245 // negative, then B straddles the equator). 246 minAbsLat := math.Max(bound.Lat.Lo, -bound.Lat.Hi) 247 248 // latGapSouth and latGapNorth measure the minimum distance from B to the 249 // south and north poles respectively. 250 latGapSouth := math.Pi/2 + bound.Lat.Lo 251 latGapNorth := math.Pi/2 - bound.Lat.Hi 252 253 if minAbsLat >= 0 { 254 // The bound B does not straddle the equator. In this case the minimum 255 // distance is between one endpoint of the latitude edge in B closest to 256 // the equator and the other endpoint of that edge in B'. The latitude 257 // distance between these two points is 2*minAbsLat, and the longitude 258 // distance is lngGap. We could compute the distance exactly using the 259 // Haversine formula, but then we would need to bound the errors in that 260 // calculation. Since we only need accuracy when the distance is very 261 // small (close to 4.309 * dblEpsilon), we substitute the Euclidean 262 // distance instead. This gives us a right triangle XYZ with two edges of 263 // length x = 2*minAbsLat and y ~= lngGap. The desired distance is the 264 // length of the third edge z, and we have 265 // 266 // z ~= sqrt(x^2 + y^2) >= (x + y) / sqrt(2) 267 // 268 // Therefore the region may contain nearly antipodal points only if 269 // 270 // 2*minAbsLat + lngGap < sqrt(2) * 4.309 * dblEpsilon 271 // ~= 1.354e-15 272 // 273 // Note that because the given bound B is conservative, minAbsLat and 274 // lngGap are both lower bounds on their true values so we do not need 275 // to make any adjustments for their errors. 276 if 2*minAbsLat+lngGap < 1.354e-15 { 277 return FullRect() 278 } 279 } else if lngGap >= math.Pi/2 { 280 // B spans at most Pi/2 in longitude. The minimum distance is always 281 // between one corner of B and the diagonally opposite corner of B'. We 282 // use the same distance approximation that we used above; in this case 283 // we have an obtuse triangle XYZ with two edges of length x = latGapSouth 284 // and y = latGapNorth, and angle Z >= Pi/2 between them. We then have 285 // 286 // z >= sqrt(x^2 + y^2) >= (x + y) / sqrt(2) 287 // 288 // Unlike the case above, latGapSouth and latGapNorth are not lower bounds 289 // (because of the extra addition operation, and because math.Pi/2 is not 290 // exactly equal to Pi/2); they can exceed their true values by up to 291 // 0.75 * dblEpsilon. Putting this all together, the region may contain 292 // nearly antipodal points only if 293 // 294 // latGapSouth + latGapNorth < (sqrt(2) * 4.309 + 1.5) * dblEpsilon 295 // ~= 1.687e-15 296 if latGapSouth+latGapNorth < 1.687e-15 { 297 return FullRect() 298 } 299 } else { 300 // Otherwise we know that (1) the bound straddles the equator and (2) its 301 // width in longitude is at least Pi/2. In this case the minimum 302 // distance can occur either between a corner of B and the diagonally 303 // opposite corner of B' (as in the case above), or between a corner of B 304 // and the opposite longitudinal edge reflected in B'. It is sufficient 305 // to only consider the corner-edge case, since this distance is also a 306 // lower bound on the corner-corner distance when that case applies. 307 308 // Consider the spherical triangle XYZ where X is a corner of B with 309 // minimum absolute latitude, Y is the closest pole to X, and Z is the 310 // point closest to X on the opposite longitudinal edge of B'. This is a 311 // right triangle (Z = Pi/2), and from the spherical law of sines we have 312 // 313 // sin(z) / sin(Z) = sin(y) / sin(Y) 314 // sin(maxLatGap) / 1 = sin(dMin) / sin(lngGap) 315 // sin(dMin) = sin(maxLatGap) * sin(lngGap) 316 // 317 // where "maxLatGap" = max(latGapSouth, latGapNorth) and "dMin" is the 318 // desired minimum distance. Now using the facts that sin(t) >= (2/Pi)*t 319 // for 0 <= t <= Pi/2, that we only need an accurate approximation when 320 // at least one of "maxLatGap" or lngGap is extremely small (in which 321 // case sin(t) ~= t), and recalling that "maxLatGap" has an error of up 322 // to 0.75 * dblEpsilon, we want to test whether 323 // 324 // maxLatGap * lngGap < (4.309 + 0.75) * (Pi/2) * dblEpsilon 325 // ~= 1.765e-15 326 if math.Max(latGapSouth, latGapNorth)*lngGap < 1.765e-15 { 327 return FullRect() 328 } 329 } 330 // Next we need to check whether the subregion might contain any edges that 331 // span (math.Pi - 2 * dblEpsilon) radians or more in longitude, since AddPoint 332 // sets the longitude bound to Full in that case. This corresponds to 333 // testing whether (lngGap <= 0) in lngExpansion below. 334 335 // Otherwise, the maximum latitude error in AddPoint is 4.8 * dblEpsilon. 336 // In the worst case, the errors when computing the latitude bound for a 337 // subregion could go in the opposite direction as the errors when computing 338 // the bound for the original region, so we need to double this value. 339 // (More analysis shows that it's okay to round down to a multiple of 340 // dblEpsilon.) 341 // 342 // For longitude, we rely on the fact that atan2 is correctly rounded and 343 // therefore no additional bounds expansion is necessary. 344 345 latExpansion := 9 * dblEpsilon 346 lngExpansion := 0.0 347 if lngGap <= 0 { 348 lngExpansion = math.Pi 349 } 350 return bound.expanded(LatLng{s1.Angle(latExpansion), s1.Angle(lngExpansion)}).PolarClosure() 351 } 352