in s2/rect_bounder.go [68:196]
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})
}