/*
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()
}
