diff options
Diffstat (limited to 'vendor/github.com/golang/geo/s2/rect_bounder.go')
| -rw-r--r-- | vendor/github.com/golang/geo/s2/rect_bounder.go | 352 |
1 files changed, 352 insertions, 0 deletions
diff --git a/vendor/github.com/golang/geo/s2/rect_bounder.go b/vendor/github.com/golang/geo/s2/rect_bounder.go new file mode 100644 index 000000000..419dea0c1 --- /dev/null +++ b/vendor/github.com/golang/geo/s2/rect_bounder.go @@ -0,0 +1,352 @@ +// Copyright 2017 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 ( + "math" + + "github.com/golang/geo/r1" + "github.com/golang/geo/r3" + "github.com/golang/geo/s1" +) + +// RectBounder is used to compute a bounding rectangle that contains all edges +// defined by a vertex chain (v0, v1, v2, ...). All vertices must be unit length. +// Note that the bounding rectangle of an edge can be larger than the bounding +// rectangle of its endpoints, e.g. consider an edge that passes through the North Pole. +// +// The bounds are calculated conservatively to account for numerical errors +// when points are converted to LatLngs. More precisely, this function +// guarantees the following: +// Let L be a closed edge chain (Loop) such that the interior of the loop does +// not contain either pole. Now if P is any point such that L.ContainsPoint(P), +// then RectBound(L).ContainsPoint(LatLngFromPoint(P)). +type RectBounder struct { + // The previous vertex in the chain. + a Point + // The previous vertex latitude longitude. + aLL LatLng + bound Rect +} + +// NewRectBounder returns a new instance of a RectBounder. +func NewRectBounder() *RectBounder { + return &RectBounder{ + bound: EmptyRect(), + } +} + +// maxErrorForTests returns the maximum error in RectBound provided that the +// result does not include either pole. It is only used for testing purposes +func (r *RectBounder) maxErrorForTests() LatLng { + // The maximum error in the latitude calculation is + // 3.84 * dblEpsilon for the PointCross calculation + // 0.96 * dblEpsilon for the Latitude calculation + // 5 * dblEpsilon added by AddPoint/RectBound to compensate for error + // ----------------- + // 9.80 * dblEpsilon maximum error in result + // + // The maximum error in the longitude calculation is dblEpsilon. RectBound + // does not do any expansion because this isn't necessary in order to + // bound the *rounded* longitudes of contained points. + return LatLng{10 * dblEpsilon * s1.Radian, 1 * dblEpsilon * s1.Radian} +} + +// AddPoint adds the given point to the chain. The Point must be unit length. +func (r *RectBounder) AddPoint(b Point) { + bLL := LatLngFromPoint(b) + + if r.bound.IsEmpty() { + r.a = b + r.aLL = bLL + r.bound = r.bound.AddPoint(bLL) + return + } + + // First compute the cross product N = A x B robustly. This is the normal + // to the great circle through A and B. We don't use RobustSign + // since that method returns an arbitrary vector orthogonal to A if the two + // vectors are proportional, and we want the zero vector in that case. + n := r.a.Sub(b.Vector).Cross(r.a.Add(b.Vector)) // N = 2 * (A x B) + + // The relative error in N gets large as its norm gets very small (i.e., + // when the two points are nearly identical or antipodal). We handle this + // by choosing a maximum allowable error, and if the error is greater than + // this we fall back to a different technique. Since it turns out that + // the other sources of error in converting the normal to a maximum + // latitude add up to at most 1.16 * dblEpsilon, and it is desirable to + // have the total error be a multiple of dblEpsilon, we have chosen to + // limit the maximum error in the normal to be 3.84 * dblEpsilon. + // It is possible to show that the error is less than this when + // + // n.Norm() >= 8 * sqrt(3) / (3.84 - 0.5 - sqrt(3)) * dblEpsilon + // = 1.91346e-15 (about 8.618 * dblEpsilon) + nNorm := n.Norm() + if nNorm < 1.91346e-15 { + // A and B are either nearly identical or nearly antipodal (to within + // 4.309 * dblEpsilon, or about 6 nanometers on the earth's surface). + if r.a.Dot(b.Vector) < 0 { + // The two points are nearly antipodal. The easiest solution is to + // assume that the edge between A and B could go in any direction + // around the sphere. + r.bound = FullRect() + } else { + // The two points are nearly identical (to within 4.309 * dblEpsilon). + // In this case we can just use the bounding rectangle of the points, + // since after the expansion done by GetBound this Rect is + // guaranteed to include the (lat,lng) values of all points along AB. + r.bound = r.bound.Union(RectFromLatLng(r.aLL).AddPoint(bLL)) + } + r.a = b + r.aLL = bLL + return + } + + // Compute the longitude range spanned by AB. + lngAB := s1.EmptyInterval().AddPoint(r.aLL.Lng.Radians()).AddPoint(bLL.Lng.Radians()) + if lngAB.Length() >= math.Pi-2*dblEpsilon { + // The points lie on nearly opposite lines of longitude to within the + // maximum error of the calculation. The easiest solution is to assume + // that AB could go on either side of the pole. + lngAB = s1.FullInterval() + } + + // Next we compute the latitude range spanned by the edge AB. We start + // with the range spanning the two endpoints of the edge: + latAB := r1.IntervalFromPoint(r.aLL.Lat.Radians()).AddPoint(bLL.Lat.Radians()) + + // This is the desired range unless the edge AB crosses the plane + // through N and the Z-axis (which is where the great circle through A + // and B attains its minimum and maximum latitudes). To test whether AB + // crosses this plane, we compute a vector M perpendicular to this + // plane and then project A and B onto it. + m := n.Cross(r3.Vector{0, 0, 1}) + mA := m.Dot(r.a.Vector) + mB := m.Dot(b.Vector) + + // We want to test the signs of "mA" and "mB", so we need to bound + // the error in these calculations. It is possible to show that the + // total error is bounded by + // + // (1 + sqrt(3)) * dblEpsilon * nNorm + 8 * sqrt(3) * (dblEpsilon**2) + // = 6.06638e-16 * nNorm + 6.83174e-31 + + mError := 6.06638e-16*nNorm + 6.83174e-31 + if mA*mB < 0 || math.Abs(mA) <= mError || math.Abs(mB) <= mError { + // Minimum/maximum latitude *may* occur in the edge interior. + // + // The maximum latitude is 90 degrees minus the latitude of N. We + // compute this directly using atan2 in order to get maximum accuracy + // near the poles. + // + // Our goal is compute a bound that contains the computed latitudes of + // all S2Points P that pass the point-in-polygon containment test. + // There are three sources of error we need to consider: + // - the directional error in N (at most 3.84 * dblEpsilon) + // - converting N to a maximum latitude + // - computing the latitude of the test point P + // The latter two sources of error are at most 0.955 * dblEpsilon + // individually, but it is possible to show by a more complex analysis + // that together they can add up to at most 1.16 * dblEpsilon, for a + // total error of 5 * dblEpsilon. + // + // We add 3 * dblEpsilon to the bound here, and GetBound() will pad + // the bound by another 2 * dblEpsilon. + maxLat := math.Min( + math.Atan2(math.Sqrt(n.X*n.X+n.Y*n.Y), math.Abs(n.Z))+3*dblEpsilon, + math.Pi/2) + + // In order to get tight bounds when the two points are close together, + // we also bound the min/max latitude relative to the latitudes of the + // endpoints A and B. First we compute the distance between A and B, + // and then we compute the maximum change in latitude between any two + // points along the great circle that are separated by this distance. + // This gives us a latitude change "budget". Some of this budget must + // be spent getting from A to B; the remainder bounds the round-trip + // distance (in latitude) from A or B to the min or max latitude + // attained along the edge AB. + latBudget := 2 * math.Asin(0.5*(r.a.Sub(b.Vector)).Norm()*math.Sin(maxLat)) + maxDelta := 0.5*(latBudget-latAB.Length()) + dblEpsilon + + // Test whether AB passes through the point of maximum latitude or + // minimum latitude. If the dot product(s) are small enough then the + // result may be ambiguous. + if mA <= mError && mB >= -mError { + latAB.Hi = math.Min(maxLat, latAB.Hi+maxDelta) + } + if mB <= mError && mA >= -mError { + latAB.Lo = math.Max(-maxLat, latAB.Lo-maxDelta) + } + } + r.a = b + r.aLL = bLL + r.bound = r.bound.Union(Rect{latAB, lngAB}) +} + +// RectBound returns the bounding rectangle of the edge chain that connects the +// vertices defined so far. This bound satisfies the guarantee made +// above, i.e. if the edge chain defines a Loop, then the bound contains +// the LatLng coordinates of all Points contained by the loop. +func (r *RectBounder) RectBound() Rect { + return r.bound.expanded(LatLng{s1.Angle(2 * dblEpsilon), 0}).PolarClosure() +} + +// ExpandForSubregions expands a bounding Rect so that it is guaranteed to +// contain the bounds of any subregion whose bounds are computed using +// ComputeRectBound. For example, consider a loop L that defines a square. +// GetBound ensures that if a point P is contained by this square, then +// LatLngFromPoint(P) is contained by the bound. But now consider a diamond +// shaped loop S contained by L. It is possible that GetBound returns a +// *larger* bound for S than it does for L, due to rounding errors. This +// method expands the bound for L so that it is guaranteed to contain the +// bounds of any subregion S. +// +// More precisely, if L is a loop that does not contain either pole, and S +// is a loop such that L.Contains(S), then +// +// ExpandForSubregions(L.RectBound).Contains(S.RectBound). +// +func ExpandForSubregions(bound Rect) Rect { + // Empty bounds don't need expansion. + if bound.IsEmpty() { + return bound + } + + // First we need to check whether the bound B contains any nearly-antipodal + // points (to within 4.309 * dblEpsilon). If so then we need to return + // FullRect, since the subregion might have an edge between two + // such points, and AddPoint returns Full for such edges. Note that + // this can happen even if B is not Full for example, consider a loop + // that defines a 10km strip straddling the equator extending from + // longitudes -100 to +100 degrees. + // + // It is easy to check whether B contains any antipodal points, but checking + // for nearly-antipodal points is trickier. Essentially we consider the + // original bound B and its reflection through the origin B', and then test + // whether the minimum distance between B and B' is less than 4.309 * dblEpsilon. + + // lngGap is a lower bound on the longitudinal distance between B and its + // reflection B'. (2.5 * dblEpsilon is the maximum combined error of the + // endpoint longitude calculations and the Length call.) + lngGap := math.Max(0, math.Pi-bound.Lng.Length()-2.5*dblEpsilon) + + // minAbsLat is the minimum distance from B to the equator (if zero or + // negative, then B straddles the equator). + minAbsLat := math.Max(bound.Lat.Lo, -bound.Lat.Hi) + + // latGapSouth and latGapNorth measure the minimum distance from B to the + // south and north poles respectively. + latGapSouth := math.Pi/2 + bound.Lat.Lo + latGapNorth := math.Pi/2 - bound.Lat.Hi + + if minAbsLat >= 0 { + // The bound B does not straddle the equator. In this case the minimum + // distance is between one endpoint of the latitude edge in B closest to + // the equator and the other endpoint of that edge in B'. The latitude + // distance between these two points is 2*minAbsLat, and the longitude + // distance is lngGap. We could compute the distance exactly using the + // Haversine formula, but then we would need to bound the errors in that + // calculation. Since we only need accuracy when the distance is very + // small (close to 4.309 * dblEpsilon), we substitute the Euclidean + // distance instead. This gives us a right triangle XYZ with two edges of + // length x = 2*minAbsLat and y ~= lngGap. The desired distance is the + // length of the third edge z, and we have + // + // z ~= sqrt(x^2 + y^2) >= (x + y) / sqrt(2) + // + // Therefore the region may contain nearly antipodal points only if + // + // 2*minAbsLat + lngGap < sqrt(2) * 4.309 * dblEpsilon + // ~= 1.354e-15 + // + // Note that because the given bound B is conservative, minAbsLat and + // lngGap are both lower bounds on their true values so we do not need + // to make any adjustments for their errors. + if 2*minAbsLat+lngGap < 1.354e-15 { + return FullRect() + } + } else if lngGap >= math.Pi/2 { + // B spans at most Pi/2 in longitude. The minimum distance is always + // between one corner of B and the diagonally opposite corner of B'. We + // use the same distance approximation that we used above; in this case + // we have an obtuse triangle XYZ with two edges of length x = latGapSouth + // and y = latGapNorth, and angle Z >= Pi/2 between them. We then have + // + // z >= sqrt(x^2 + y^2) >= (x + y) / sqrt(2) + // + // Unlike the case above, latGapSouth and latGapNorth are not lower bounds + // (because of the extra addition operation, and because math.Pi/2 is not + // exactly equal to Pi/2); they can exceed their true values by up to + // 0.75 * dblEpsilon. Putting this all together, the region may contain + // nearly antipodal points only if + // + // latGapSouth + latGapNorth < (sqrt(2) * 4.309 + 1.5) * dblEpsilon + // ~= 1.687e-15 + if latGapSouth+latGapNorth < 1.687e-15 { + return FullRect() + } + } else { + // Otherwise we know that (1) the bound straddles the equator and (2) its + // width in longitude is at least Pi/2. In this case the minimum + // distance can occur either between a corner of B and the diagonally + // opposite corner of B' (as in the case above), or between a corner of B + // and the opposite longitudinal edge reflected in B'. It is sufficient + // to only consider the corner-edge case, since this distance is also a + // lower bound on the corner-corner distance when that case applies. + + // Consider the spherical triangle XYZ where X is a corner of B with + // minimum absolute latitude, Y is the closest pole to X, and Z is the + // point closest to X on the opposite longitudinal edge of B'. This is a + // right triangle (Z = Pi/2), and from the spherical law of sines we have + // + // sin(z) / sin(Z) = sin(y) / sin(Y) + // sin(maxLatGap) / 1 = sin(dMin) / sin(lngGap) + // sin(dMin) = sin(maxLatGap) * sin(lngGap) + // + // where "maxLatGap" = max(latGapSouth, latGapNorth) and "dMin" is the + // desired minimum distance. Now using the facts that sin(t) >= (2/Pi)*t + // for 0 <= t <= Pi/2, that we only need an accurate approximation when + // at least one of "maxLatGap" or lngGap is extremely small (in which + // case sin(t) ~= t), and recalling that "maxLatGap" has an error of up + // to 0.75 * dblEpsilon, we want to test whether + // + // maxLatGap * lngGap < (4.309 + 0.75) * (Pi/2) * dblEpsilon + // ~= 1.765e-15 + if math.Max(latGapSouth, latGapNorth)*lngGap < 1.765e-15 { + return FullRect() + } + } + // Next we need to check whether the subregion might contain any edges that + // span (math.Pi - 2 * dblEpsilon) radians or more in longitude, since AddPoint + // sets the longitude bound to Full in that case. This corresponds to + // testing whether (lngGap <= 0) in lngExpansion below. + + // Otherwise, the maximum latitude error in AddPoint is 4.8 * dblEpsilon. + // In the worst case, the errors when computing the latitude bound for a + // subregion could go in the opposite direction as the errors when computing + // the bound for the original region, so we need to double this value. + // (More analysis shows that it's okay to round down to a multiple of + // dblEpsilon.) + // + // For longitude, we rely on the fact that atan2 is correctly rounded and + // therefore no additional bounds expansion is necessary. + + latExpansion := 9 * dblEpsilon + lngExpansion := 0.0 + if lngGap <= 0 { + lngExpansion = math.Pi + } + return bound.expanded(LatLng{s1.Angle(latExpansion), s1.Angle(lngExpansion)}).PolarClosure() +} |