/*
Copyright 2015 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/r2"
	"github.com/golang/geo/r3"
	"github.com/golang/geo/s1"
)

const (
	// edgeClipErrorUVCoord is the maximum error in a u- or v-coordinate
	// compared to the exact result, assuming that the points A and B are in
	// the rectangle [-1,1]x[1,1] or slightly outside it (by 1e-10 or less).
	edgeClipErrorUVCoord = 2.25 * dblEpsilon

	// edgeClipErrorUVDist is the maximum distance from a clipped point to
	// the corresponding exact result. It is equal to the error in a single
	// coordinate because at most one coordinate is subject to error.
	edgeClipErrorUVDist = 2.25 * dblEpsilon

	// faceClipErrorRadians is the maximum angle between a returned vertex
	// and the nearest point on the exact edge AB. It is equal to the
	// maximum directional error in PointCross, plus the error when
	// projecting points onto a cube face.
	faceClipErrorRadians = 3 * dblEpsilon

	// faceClipErrorDist is the same angle expressed as a maximum distance
	// in (u,v)-space. In other words, a returned vertex is at most this far
	// from the exact edge AB projected into (u,v)-space.
	faceClipErrorUVDist = 9 * dblEpsilon

	// faceClipErrorUVCoord is the maximum angle between a returned vertex
	// and the nearest point on the exact edge AB expressed as the maximum error
	// in an individual u- or v-coordinate. In other words, for each
	// returned vertex there is a point on the exact edge AB whose u- and
	// v-coordinates differ from the vertex by at most this amount.
	faceClipErrorUVCoord = 9.0 * (1.0 / math.Sqrt2) * dblEpsilon

	// intersectsRectErrorUVDist is the maximum error when computing if a point
	// intersects with a given Rect. If some point of AB is inside the
	// rectangle by at least this distance, the result is guaranteed to be true;
	// if all points of AB are outside the rectangle by at least this distance,
	// the result is guaranteed to be false. This bound assumes that rect is
	// a subset of the rectangle [-1,1]x[-1,1] or extends slightly outside it
	// (e.g., by 1e-10 or less).
	intersectsRectErrorUVDist = 3 * math.Sqrt2 * dblEpsilon

	// intersectionError can be set somewhat arbitrarily, because the algorithm
	// uses more precision if necessary in order to achieve the specified error.
	// The only strict requirement is that intersectionError >= dblEpsilon
	// radians. However, using a larger error tolerance makes the algorithm more
	// efficient because it reduces the number of cases where exact arithmetic is
	// needed.
	intersectionError = s1.Angle(4 * dblEpsilon)

	// intersectionMergeRadius is used to ensure that intersection points that
	// are supposed to be coincident are merged back together into a single
	// vertex. This is required in order for various polygon operations (union,
	// intersection, etc) to work correctly. It is twice the intersection error
	// because two coincident intersection points might have errors in
	// opposite directions.
	intersectionMergeRadius = 2 * intersectionError
)

// SimpleCrossing reports whether edge AB crosses CD at a point that is interior
// to both edges. Properties:
//
//  (1) SimpleCrossing(b,a,c,d) == SimpleCrossing(a,b,c,d)
//  (2) SimpleCrossing(c,d,a,b) == SimpleCrossing(a,b,c,d)
func SimpleCrossing(a, b, c, d Point) bool {
	// We compute the equivalent of Sign for triangles ACB, CBD, BDA,
	// and DAC. All of these triangles need to have the same orientation
	// (CW or CCW) for an intersection to exist.

	ab := a.Vector.Cross(b.Vector)
	acb := -(ab.Dot(c.Vector))
	bda := ab.Dot(d.Vector)
	if acb*bda <= 0 {
		return false
	}

	cd := c.Vector.Cross(d.Vector)
	cbd := -(cd.Dot(b.Vector))
	dac := cd.Dot(a.Vector)
	return (acb*cbd > 0) && (acb*dac > 0)
}

// VertexCrossing reports whether two edges "cross" in such a way that point-in-polygon
// containment tests can be implemented by counting the number of edge crossings.
//
// Given two edges AB and CD where at least two vertices are identical
// (i.e. CrossingSign(a,b,c,d) == 0), the basic rule is that a "crossing"
// occurs if AB is encountered after CD during a CCW sweep around the shared
// vertex starting from a fixed reference point.
//
// Note that according to this rule, if AB crosses CD then in general CD
// does not cross AB.  However, this leads to the correct result when
// counting polygon edge crossings.  For example, suppose that A,B,C are
// three consecutive vertices of a CCW polygon.  If we now consider the edge
// crossings of a segment BP as P sweeps around B, the crossing number
// changes parity exactly when BP crosses BA or BC.
//
// Useful properties of VertexCrossing (VC):
//
//  (1) VC(a,a,c,d) == VC(a,b,c,c) == false
//  (2) VC(a,b,a,b) == VC(a,b,b,a) == true
//  (3) VC(a,b,c,d) == VC(a,b,d,c) == VC(b,a,c,d) == VC(b,a,d,c)
//  (3) If exactly one of a,b equals one of c,d, then exactly one of
//      VC(a,b,c,d) and VC(c,d,a,b) is true
//
// It is an error to call this method with 4 distinct vertices.
func VertexCrossing(a, b, c, d Point) bool {
	// If A == B or C == D there is no intersection. We need to check this
	// case first in case 3 or more input points are identical.
	if a.ApproxEqual(b) || c.ApproxEqual(d) {
		return false
	}

	// If any other pair of vertices is equal, there is a crossing if and only
	// if OrderedCCW indicates that the edge AB is further CCW around the
	// shared vertex O (either A or B) than the edge CD, starting from an
	// arbitrary fixed reference point.
	switch {
	case a.ApproxEqual(d):
		return OrderedCCW(Point{a.Ortho()}, c, b, a)
	case b.ApproxEqual(c):
		return OrderedCCW(Point{b.Ortho()}, d, a, b)
	case a.ApproxEqual(c):
		return OrderedCCW(Point{a.Ortho()}, d, b, a)
	case b.ApproxEqual(d):
		return OrderedCCW(Point{b.Ortho()}, c, a, b)
	}

	return false
}

// DistanceFraction returns the distance ratio of the point X along an edge AB.
// If X is on the line segment AB, this is the fraction T such
// that X == Interpolate(T, A, B).
//
// This requires that A and B are distinct.
func DistanceFraction(x, a, b Point) float64 {
	d0 := x.Angle(a.Vector)
	d1 := x.Angle(b.Vector)
	return float64(d0 / (d0 + d1))
}

// Interpolate returns the point X along the line segment AB whose distance from A
// is the given fraction "t" of the distance AB. Does NOT require that "t" be
// between 0 and 1. Note that all distances are measured on the surface of
// the sphere, so this is more complicated than just computing (1-t)*a + t*b
// and normalizing the result.
func Interpolate(t float64, a, b Point) Point {
	if t == 0 {
		return a
	}
	if t == 1 {
		return b
	}
	ab := a.Angle(b.Vector)
	return InterpolateAtDistance(s1.Angle(t)*ab, a, b)
}

// InterpolateAtDistance returns the point X along the line segment AB whose
// distance from A is the angle ax.
func InterpolateAtDistance(ax s1.Angle, a, b Point) Point {
	aRad := ax.Radians()

	// Use PointCross to compute the tangent vector at A towards B. The
	// result is always perpendicular to A, even if A=B or A=-B, but it is not
	// necessarily unit length. (We effectively normalize it below.)
	normal := a.PointCross(b)
	tangent := normal.Vector.Cross(a.Vector)

	// Now compute the appropriate linear combination of A and "tangent". With
	// infinite precision the result would always be unit length, but we
	// normalize it anyway to ensure that the error is within acceptable bounds.
	// (Otherwise errors can build up when the result of one interpolation is
	// fed into another interpolation.)
	return Point{(a.Mul(math.Cos(aRad)).Add(tangent.Mul(math.Sin(aRad) / tangent.Norm()))).Normalize()}
}

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

// EdgeCrosser allows edges to be efficiently tested for intersection with a
// given fixed edge AB. It is especially efficient when testing for
// intersection with an edge chain connecting vertices v0, v1, v2, ...
type EdgeCrosser struct {
	a   Point
	b   Point
	aXb Point

	// To reduce the number of calls to expensiveSign, we compute an
	// outward-facing tangent at A and B if necessary. If the plane
	// perpendicular to one of these tangents separates AB from CD (i.e., one
	// edge on each side) then there is no intersection.
	aTangent Point // Outward-facing tangent at A.
	bTangent Point // Outward-facing tangent at B.

	// The fields below are updated for each vertex in the chain.
	c   Point     // Previous vertex in the vertex chain.
	acb Direction // The orientation of triangle ACB.
}

// NewEdgeCrosser returns an EdgeCrosser with the fixed edge AB.
func NewEdgeCrosser(a, b Point) *EdgeCrosser {
	norm := a.PointCross(b)
	return &EdgeCrosser{
		a:        a,
		b:        b,
		aXb:      Point{a.Cross(b.Vector)},
		aTangent: Point{a.Cross(norm.Vector)},
		bTangent: Point{norm.Cross(b.Vector)},
	}
}

// A Crossing indicates how edges cross.
type Crossing int

const (
	// Cross means the edges cross.
	Cross Crossing = iota
	// MaybeCross means two vertices from different edges are the same.
	MaybeCross
	// DoNotCross means the edges do not cross.
	DoNotCross
)

// CrossingSign reports whether the edge AB intersects the edge CD. If any two
// vertices from different edges are the same, returns MaybeCross. If either edge
// is degenerate (A == B or C == D), returns either DoNotCross or MaybeCross.
//
// Properties of CrossingSign:
//
//  (1) CrossingSign(b,a,c,d) == CrossingSign(a,b,c,d)
//  (2) CrossingSign(c,d,a,b) == CrossingSign(a,b,c,d)
//  (3) CrossingSign(a,b,c,d) == MaybeCross if a==c, a==d, b==c, b==d
//  (3) CrossingSign(a,b,c,d) == DoNotCross or MaybeCross if a==b or c==d
//
// Note that if you want to check an edge against a chain of other edges,
// it is slightly more efficient to use the single-argument version
// ChainCrossingSign below.
func (e *EdgeCrosser) CrossingSign(c, d Point) Crossing {
	if c != e.c {
		e.RestartAt(c)
	}
	return e.ChainCrossingSign(d)
}

// EdgeOrVertexCrossing reports whether if CrossingSign(c, d) > 0, or AB and
// CD share a vertex and VertexCrossing(a, b, c, d) is true.
//
// This method extends the concept of a "crossing" to the case where AB
// and CD have a vertex in common. The two edges may or may not cross,
// according to the rules defined in VertexCrossing above. The rules
// are designed so that point containment tests can be implemented simply
// by counting edge crossings. Similarly, determining whether one edge
// chain crosses another edge chain can be implemented by counting.
func (e *EdgeCrosser) EdgeOrVertexCrossing(c, d Point) bool {
	if c != e.c {
		e.RestartAt(c)
	}
	return e.EdgeOrVertexChainCrossing(d)
}

// NewChainEdgeCrosser is a convenience constructor that uses AB as the fixed edge,
// and C as the first vertex of the vertex chain (equivalent to calling RestartAt(c)).
//
// You don't need to use this or any of the chain functions unless you're trying to
// squeeze out every last drop of performance. Essentially all you are saving is a test
// whether the first vertex of the current edge is the same as the second vertex of the
// previous edge.
func NewChainEdgeCrosser(a, b, c Point) *EdgeCrosser {
	e := NewEdgeCrosser(a, b)
	e.RestartAt(c)
	return e
}

// RestartAt sets the current point of the edge crosser to be c.
// Call this method when your chain 'jumps' to a new place.
// The argument must point to a value that persists until the next call.
func (e *EdgeCrosser) RestartAt(c Point) {
	e.c = c
	e.acb = -triageSign(e.a, e.b, e.c)
}

// ChainCrossingSign is like CrossingSign, but uses the last vertex passed to one of
// the crossing methods (or RestartAt) as the first vertex of the current edge.
func (e *EdgeCrosser) ChainCrossingSign(d Point) Crossing {
	// For there to be an edge crossing, the triangles ACB, CBD, BDA, DAC must
	// all be oriented the same way (CW or CCW). We keep the orientation of ACB
	// as part of our state. When each new point D arrives, we compute the
	// orientation of BDA and check whether it matches ACB. This checks whether
	// the points C and D are on opposite sides of the great circle through AB.

	// Recall that triageSign is invariant with respect to rotating its
	// arguments, i.e. ABD has the same orientation as BDA.
	bda := triageSign(e.a, e.b, d)
	if e.acb == -bda && bda != Indeterminate {
		// The most common case -- triangles have opposite orientations. Save the
		// current vertex D as the next vertex C, and also save the orientation of
		// the new triangle ACB (which is opposite to the current triangle BDA).
		e.c = d
		e.acb = -bda
		return DoNotCross
	}
	return e.crossingSign(d, bda)
}

// EdgeOrVertexChainCrossing is like EdgeOrVertexCrossing, but uses the last vertex
// passed to one of the crossing methods (or RestartAt) as the first vertex of the current edge.
func (e *EdgeCrosser) EdgeOrVertexChainCrossing(d Point) bool {
	// We need to copy e.c since it is clobbered by ChainCrossingSign.
	c := e.c
	switch e.ChainCrossingSign(d) {
	case DoNotCross:
		return false
	case Cross:
		return true
	}
	return VertexCrossing(e.a, e.b, c, d)
}

// crossingSign handle the slow path of CrossingSign.
func (e *EdgeCrosser) crossingSign(d Point, bda Direction) Crossing {
	// Compute the actual result, and then save the current vertex D as the next
	// vertex C, and save the orientation of the next triangle ACB (which is
	// opposite to the current triangle BDA).
	defer func() {
		e.c = d
		e.acb = -bda
	}()

	// At this point, a very common situation is that A,B,C,D are four points on
	// a line such that AB does not overlap CD. (For example, this happens when
	// a line or curve is sampled finely, or when geometry is constructed by
	// computing the union of S2CellIds.) Most of the time, we can determine
	// that AB and CD do not intersect using the two outward-facing
	// tangents at A and B (parallel to AB) and testing whether AB and CD are on
	// opposite sides of the plane perpendicular to one of these tangents. This
	// is moderately expensive but still much cheaper than expensiveSign.

	// The error in RobustCrossProd is insignificant. The maximum error in
	// the call to CrossProd (i.e., the maximum norm of the error vector) is
	// (0.5 + 1/sqrt(3)) * dblEpsilon. The maximum error in each call to
	// DotProd below is dblEpsilon. (There is also a small relative error
	// term that is insignificant because we are comparing the result against a
	// constant that is very close to zero.)
	maxError := (1.5 + 1/math.Sqrt(3)) * dblEpsilon
	if (e.c.Dot(e.aTangent.Vector) > maxError && d.Dot(e.aTangent.Vector) > maxError) || (e.c.Dot(e.bTangent.Vector) > maxError && d.Dot(e.bTangent.Vector) > maxError) {
		return DoNotCross
	}

	// Otherwise, eliminate the cases where any two vertices are equal. (These
	// cases could be handled in the code below, but since ExpensiveSign lives
	// up to its name we would rather avoid calling it if possible.)
	//
	// These are the cases where two vertices from different edges are equal.
	if e.a == e.c || e.a == d || e.b == e.c || e.b == d {
		return MaybeCross
	}

	// These are the cases where an input edge is degenerate. (Note that in
	// most cases, if CD is degenerate then this method is not even called
	// because acb and bda have different signs. That's why this method is
	// documented to return either MaybeCross or DoNotCross when an input
	// edge is degenerate.)
	if e.a == e.b || e.c == d {
		return MaybeCross
	}

	// Otherwise it's time to break out the big guns.
	if e.acb == Indeterminate {
		e.acb = -expensiveSign(e.a, e.b, e.c)
	}
	if bda == Indeterminate {
		bda = expensiveSign(e.a, e.b, d)
	}

	if bda != e.acb {
		return DoNotCross
	}

	cbd := -RobustSign(e.c, d, e.b)
	if cbd != e.acb {
		return DoNotCross
	}
	dac := RobustSign(e.c, d, e.a)
	if dac != e.acb {
		return DoNotCross
	}
	return Cross
}

// pointUVW represents a Point in (u,v,w) coordinate space of a cube face.
type pointUVW Point

// intersectsFace reports whether a given directed line L intersects the cube face F.
// The line L is defined by its normal N in the (u,v,w) coordinates of F.
func (p pointUVW) intersectsFace() bool {
	// L intersects the [-1,1]x[-1,1] square in (u,v) if and only if the dot
	// products of N with the four corner vertices (-1,-1,1), (1,-1,1), (1,1,1),
	// and (-1,1,1) do not all have the same sign. This is true exactly when
	// |Nu| + |Nv| >= |Nw|. The code below evaluates this expression exactly.
	u := math.Abs(p.X)
	v := math.Abs(p.Y)
	w := math.Abs(p.Z)

	// We only need to consider the cases where u or v is the smallest value,
	// since if w is the smallest then both expressions below will have a
	// positive LHS and a negative RHS.
	return (v >= w-u) && (u >= w-v)
}

// intersectsOppositeEdges reports whether a directed line L intersects two
// opposite edges of a cube face F. This includs the case where L passes
// exactly through a corner vertex of F. The directed line L is defined
// by its normal N in the (u,v,w) coordinates of F.
func (p pointUVW) intersectsOppositeEdges() bool {
	// The line L intersects opposite edges of the [-1,1]x[-1,1] (u,v) square if
	// and only exactly two of the corner vertices lie on each side of L. This
	// is true exactly when ||Nu| - |Nv|| >= |Nw|. The code below evaluates this
	// expression exactly.
	u := math.Abs(p.X)
	v := math.Abs(p.Y)
	w := math.Abs(p.Z)

	// If w is the smallest, the following line returns an exact result.
	if math.Abs(u-v) != w {
		return math.Abs(u-v) >= w
	}

	// Otherwise u - v = w exactly, or w is not the smallest value. In either
	// case the following returns the correct result.
	if u >= v {
		return u-w >= v
	}
	return v-w >= u
}

// axis represents the possible results of exitAxis.
type axis int

const (
	axisU axis = iota
	axisV
)

// exitAxis reports which axis the directed line L exits the cube face F on.
// The directed line L is represented by its CCW normal N in the (u,v,w) coordinates
// of F. It returns axisU if L exits through the u=-1 or u=+1 edge, and axisV if L exits
// through the v=-1 or v=+1 edge. Either result is acceptable if L exits exactly
// through a corner vertex of the cube face.
func (p pointUVW) exitAxis() axis {
	if p.intersectsOppositeEdges() {
		// The line passes through through opposite edges of the face.
		// It exits through the v=+1 or v=-1 edge if the u-component of N has a
		// larger absolute magnitude than the v-component.
		if math.Abs(p.X) >= math.Abs(p.Y) {
			return axisV
		}
		return axisU
	}

	// The line passes through through two adjacent edges of the face.
	// It exits the v=+1 or v=-1 edge if an even number of the components of N
	// are negative. We test this using signbit() rather than multiplication
	// to avoid the possibility of underflow.
	var x, y, z int
	if math.Signbit(p.X) {
		x = 1
	}
	if math.Signbit(p.Y) {
		y = 1
	}
	if math.Signbit(p.Z) {
		z = 1
	}

	if x^y^z == 0 {
		return axisV
	}
	return axisU
}

// exitPoint returns the UV coordinates of the point where a directed line L (represented
// by the CCW normal of this point), exits the cube face this point is derived from along
// the given axis.
func (p pointUVW) exitPoint(a axis) r2.Point {
	if a == axisU {
		u := -1.0
		if p.Y > 0 {
			u = 1.0
		}
		return r2.Point{u, (-u*p.X - p.Z) / p.Y}
	}

	v := -1.0
	if p.X < 0 {
		v = 1.0
	}
	return r2.Point{(-v*p.Y - p.Z) / p.X, v}
}

// clipDestination returns a score which is used to indicate if the clipped edge AB
// on the given face intersects the face at all. This function returns the score for
// the given endpoint, which is an integer ranging from 0 to 3. If the sum of the scores
// from both of the endpoints is 3 or more, then edge AB does not intersect this face.
//
// First, it clips the line segment AB to find the clipped destination B' on a given
// face. (The face is specified implicitly by expressing *all arguments* in the (u,v,w)
// coordinates of that face.) Second, it partially computes whether the segment AB
// intersects this face at all. The actual condition is fairly complicated, but it
// turns out that it can be expressed as a "score" that can be computed independently
// when clipping the two endpoints A and B.
func clipDestination(a, b, scaledN, aTan, bTan pointUVW, scaleUV float64) (r2.Point, int) {
	var uv r2.Point

	// Optimization: if B is within the safe region of the face, use it.
	maxSafeUVCoord := 1 - faceClipErrorUVCoord
	if b.Z > 0 {
		uv = r2.Point{b.X / b.Z, b.Y / b.Z}
		if math.Max(math.Abs(uv.X), math.Abs(uv.Y)) <= maxSafeUVCoord {
			return uv, 0
		}
	}

	// Otherwise find the point B' where the line AB exits the face.
	uv = scaledN.exitPoint(scaledN.exitAxis()).Mul(scaleUV)

	p := pointUVW(Point{r3.Vector{uv.X, uv.Y, 1.0}})

	// Determine if the exit point B' is contained within the segment. We do this
	// by computing the dot products with two inward-facing tangent vectors at A
	// and B. If either dot product is negative, we say that B' is on the "wrong
	// side" of that point. As the point B' moves around the great circle AB past
	// the segment endpoint B, it is initially on the wrong side of B only; as it
	// moves further it is on the wrong side of both endpoints; and then it is on
	// the wrong side of A only. If the exit point B' is on the wrong side of
	// either endpoint, we can't use it; instead the segment is clipped at the
	// original endpoint B.
	//
	// We reject the segment if the sum of the scores of the two endpoints is 3
	// or more. Here is what that rule encodes:
	//  - If B' is on the wrong side of A, then the other clipped endpoint A'
	//    must be in the interior of AB (otherwise AB' would go the wrong way
	//    around the circle). There is a similar rule for A'.
	//  - If B' is on the wrong side of either endpoint (and therefore we must
	//    use the original endpoint B instead), then it must be possible to
	//    project B onto this face (i.e., its w-coordinate must be positive).
	//    This rule is only necessary to handle certain zero-length edges (A=B).
	score := 0
	if p.Sub(a.Vector).Dot(aTan.Vector) < 0 {
		score = 2 // B' is on wrong side of A.
	} else if p.Sub(b.Vector).Dot(bTan.Vector) < 0 {
		score = 1 // B' is on wrong side of B.
	}

	if score > 0 { // B' is not in the interior of AB.
		if b.Z <= 0 {
			score = 3 // B cannot be projected onto this face.
		} else {
			uv = r2.Point{b.X / b.Z, b.Y / b.Z}
		}
	}

	return uv, score
}

// ClipToFace returns the (u,v) coordinates for the portion of the edge AB that
// intersects the given face, or false if the edge AB does not intersect.
// This method guarantees that the clipped vertices lie within the [-1,1]x[-1,1]
// cube face rectangle and are within faceClipErrorUVDist of the line AB, but
// the results may differ from those produced by faceSegments.
func ClipToFace(a, b Point, face int) (aUV, bUV r2.Point, intersects bool) {
	return ClipToPaddedFace(a, b, face, 0.0)
}

// ClipToPaddedFace returns the (u,v) coordinates for the portion of the edge AB that
// intersects the given face, but rather than clipping to the square [-1,1]x[-1,1]
// in (u,v) space, this method clips to [-R,R]x[-R,R] where R=(1+padding).
// Padding must be non-negative.
func ClipToPaddedFace(a, b Point, f int, padding float64) (aUV, bUV r2.Point, intersects bool) {
	// Fast path: both endpoints are on the given face.
	if face(a.Vector) == f && face(b.Vector) == f {
		au, av := validFaceXYZToUV(f, a.Vector)
		bu, bv := validFaceXYZToUV(f, b.Vector)
		return r2.Point{au, av}, r2.Point{bu, bv}, true
	}

	// Convert everything into the (u,v,w) coordinates of the given face. Note
	// that the cross product *must* be computed in the original (x,y,z)
	// coordinate system because PointCross (unlike the mathematical cross
	// product) can produce different results in different coordinate systems
	// when one argument is a linear multiple of the other, due to the use of
	// symbolic perturbations.
	normUVW := pointUVW(faceXYZtoUVW(f, a.PointCross(b)))
	aUVW := pointUVW(faceXYZtoUVW(f, a))
	bUVW := pointUVW(faceXYZtoUVW(f, b))

	// Padding is handled by scaling the u- and v-components of the normal.
	// Letting R=1+padding, this means that when we compute the dot product of
	// the normal with a cube face vertex (such as (-1,-1,1)), we will actually
	// compute the dot product with the scaled vertex (-R,-R,1). This allows
	// methods such as intersectsFace, exitAxis, etc, to handle padding
	// with no further modifications.
	scaleUV := 1 + padding
	scaledN := pointUVW{r3.Vector{X: scaleUV * normUVW.X, Y: scaleUV * normUVW.Y, Z: normUVW.Z}}
	if !scaledN.intersectsFace() {
		return aUV, bUV, false
	}

	// TODO(roberts): This is a workaround for extremely small vectors where some
	// loss of precision can occur in Normalize causing underflow. When PointCross
	// is updated to work around this, this can be removed.
	if math.Max(math.Abs(normUVW.X), math.Max(math.Abs(normUVW.Y), math.Abs(normUVW.Z))) < math.Ldexp(1, -511) {
		normUVW = pointUVW{normUVW.Mul(math.Ldexp(1, 563))}
	}

	normUVW = pointUVW{normUVW.Normalize()}

	aTan := pointUVW{normUVW.Cross(aUVW.Vector)}
	bTan := pointUVW{bUVW.Cross(normUVW.Vector)}

	// As described in clipDestination, if the sum of the scores from clipping the two
	// endpoints is 3 or more, then the segment does not intersect this face.
	aUV, aScore := clipDestination(bUVW, aUVW, pointUVW{scaledN.Mul(-1)}, bTan, aTan, scaleUV)
	bUV, bScore := clipDestination(aUVW, bUVW, scaledN, aTan, bTan, scaleUV)

	return aUV, bUV, aScore+bScore < 3
}

// interpolateDouble returns a value with the same combination of a1 and b1 as the
// given value x is of a and b. This function makes the following guarantees:
//  - If x == a, then x1 = a1 (exactly).
//  - If x == b, then x1 = b1 (exactly).
//  - If a <= x <= b, then a1 <= x1 <= b1 (even if a1 == b1).
// This requires a != b.
func interpolateDouble(x, a, b, a1, b1 float64) float64 {
	// To get results that are accurate near both A and B, we interpolate
	// starting from the closer of the two points.
	if math.Abs(a-x) <= math.Abs(b-x) {
		return a1 + (b1-a1)*(x-a)/(b-a)
	}
	return b1 + (a1-b1)*(x-b)/(a-b)
}

// updateEndpoint returns the interval with the specified endpoint updated to
// the given value. If the value lies beyond the opposite endpoint, nothing is
// changed and false is returned.
func updateEndpoint(bound r1.Interval, highEndpoint bool, value float64) (r1.Interval, bool) {
	if !highEndpoint {
		if bound.Hi < value {
			return bound, false
		}
		if bound.Lo < value {
			bound.Lo = value
		}
		return bound, true
	}

	if bound.Lo > value {
		return bound, false
	}
	if bound.Hi > value {
		bound.Hi = value
	}
	return bound, true
}

// clipBoundAxis returns the clipped versions of the bounding intervals for the given
// axes for the line segment from (a0,a1) to (b0,b1) so that neither extends beyond the
// given clip interval. negSlope is a precomputed helper variable that indicates which
// diagonal of the bounding box is spanned by AB; it is false if AB has positive slope,
// and true if AB has negative slope. If the clipping interval doesn't overlap the bounds,
// false is returned.
func clipBoundAxis(a0, b0 float64, bound0 r1.Interval, a1, b1 float64, bound1 r1.Interval,
	negSlope bool, clip r1.Interval) (bound0c, bound1c r1.Interval, updated bool) {

	if bound0.Lo < clip.Lo {
		// If the upper bound is below the clips lower bound, there is nothing to do.
		if bound0.Hi < clip.Lo {
			return bound0, bound1, false
		}
		// narrow the intervals lower bound to the clip bound.
		bound0.Lo = clip.Lo
		if bound1, updated = updateEndpoint(bound1, negSlope, interpolateDouble(clip.Lo, a0, b0, a1, b1)); !updated {
			return bound0, bound1, false
		}
	}

	if bound0.Hi > clip.Hi {
		// If the lower bound is above the clips upper bound, there is nothing to do.
		if bound0.Lo > clip.Hi {
			return bound0, bound1, false
		}
		// narrow the intervals upper bound to the clip bound.
		bound0.Hi = clip.Hi
		if bound1, updated = updateEndpoint(bound1, !negSlope, interpolateDouble(clip.Hi, a0, b0, a1, b1)); !updated {
			return bound0, bound1, false
		}
	}
	return bound0, bound1, true
}

// edgeIntersectsRect reports whether the edge defined by AB intersects the
// given closed rectangle to within the error bound.
func edgeIntersectsRect(a, b r2.Point, r r2.Rect) bool {
	// First check whether the bounds of a Rect around AB intersects the given rect.
	if !r.Intersects(r2.RectFromPoints(a, b)) {
		return false
	}

	// Otherwise AB intersects the rect if and only if all four vertices of rect
	// do not lie on the same side of the extended line AB. We test this by finding
	// the two vertices of rect with minimum and maximum projections onto the normal
	// of AB, and computing their dot products with the edge normal.
	n := b.Sub(a).Ortho()

	i := 0
	if n.X >= 0 {
		i = 1
	}
	j := 0
	if n.Y >= 0 {
		j = 1
	}

	max := n.Dot(r.VertexIJ(i, j).Sub(a))
	min := n.Dot(r.VertexIJ(1-i, 1-j).Sub(a))

	return (max >= 0) && (min <= 0)
}

// clippedEdgeBound returns the bounding rectangle of the portion of the edge defined
// by AB intersected by clip. The resulting bound may be empty. This is a convenience
// function built on top of clipEdgeBound.
func clippedEdgeBound(a, b r2.Point, clip r2.Rect) r2.Rect {
	bound := r2.RectFromPoints(a, b)
	if b1, intersects := clipEdgeBound(a, b, clip, bound); intersects {
		return b1
	}
	return r2.EmptyRect()
}

// clipEdgeBound clips an edge AB to sequence of rectangles efficiently.
// It represents the clipped edges by their bounding boxes rather than as a pair of
// endpoints. Specifically, let A'B' be some portion of an edge AB, and let bound be
// a tight bound of A'B'. This function returns the bound that is a tight bound
// of A'B' intersected with a given rectangle. If A'B' does not intersect clip,
// it returns false and the original bound.
func clipEdgeBound(a, b r2.Point, clip, bound r2.Rect) (r2.Rect, bool) {
	// negSlope indicates which diagonal of the bounding box is spanned by AB: it
	// is false if AB has positive slope, and true if AB has negative slope. This is
	// used to determine which interval endpoints need to be updated each time
	// the edge is clipped.
	negSlope := (a.X > b.X) != (a.Y > b.Y)

	b0x, b0y, up1 := clipBoundAxis(a.X, b.X, bound.X, a.Y, b.Y, bound.Y, negSlope, clip.X)
	if !up1 {
		return bound, false
	}
	b1y, b1x, up2 := clipBoundAxis(a.Y, b.Y, b0y, a.X, b.X, b0x, negSlope, clip.Y)
	if !up2 {
		return r2.Rect{b0x, b0y}, false
	}
	return r2.Rect{X: b1x, Y: b1y}, true
}

// ClipEdge returns the portion of the edge defined by AB that is contained by the
// given rectangle. If there is no intersection, false is returned and aClip and bClip
// are undefined.
func ClipEdge(a, b r2.Point, clip r2.Rect) (aClip, bClip r2.Point, intersects bool) {
	// Compute the bounding rectangle of AB, clip it, and then extract the new
	// endpoints from the clipped bound.
	bound := r2.RectFromPoints(a, b)
	if bound, intersects = clipEdgeBound(a, b, clip, bound); !intersects {
		return aClip, bClip, false
	}
	ai := 0
	if a.X > b.X {
		ai = 1
	}
	aj := 0
	if a.Y > b.Y {
		aj = 1
	}

	return bound.VertexIJ(ai, aj), bound.VertexIJ(1-ai, 1-aj), true
}

// ClosestPoint returns the point along the edge AB that is closest to the point X.
// The fractional distance of this point along the edge AB can be obtained
// using DistanceFraction.
//
// This requires that all points are unit length.
func ClosestPoint(x, a, b Point) Point {
	aXb := a.PointCross(b)
	// Find the closest point to X along the great circle through AB.
	p := x.Sub(aXb.Mul(x.Dot(aXb.Vector) / aXb.Vector.Norm2()))

	// If this point is on the edge AB, then it's the closest point.
	if Sign(aXb, a, Point{p}) && Sign(Point{p}, b, aXb) {
		return Point{p.Normalize()}
	}

	// Otherwise, the closest point is either A or B.
	if x.Sub(a.Vector).Norm2() <= x.Sub(b.Vector).Norm2() {
		return a
	}
	return b
}

// minUpdateDistanceMaxError returns the maximum error in the result of
// MaybeUpdateMinDistance (and the associated functions such as
// MaybeUpdateMinInteriorDistance, IsDistanceLess, etc), assuming that all
// input points are normalized to within the bounds guaranteed by r3.Vector's
// Normalize. The error can be added or subtracted from an s1.ChordAngle
// using its Expanded method.
func minUpdateDistanceMaxError(dist s1.ChordAngle) float64 {
	// There are two cases for the maximum error in MaybeUpdateMinDistance(),
	// depending on whether the closest point is interior to the edge.
	return math.Max(minUpdateInteriorDistanceMaxError(dist), dist.MaxPointError())
}

// minUpdateInteriorDistanceMaxError returns the maximum error in the result of
// MaybeUpdateMinInteriorDistance, assuming that all input points are normalized
// to within the bounds guaranteed by Point's Normalize. The error can be added
// or subtracted from an s1.ChordAngle using its Expanded method.
func minUpdateInteriorDistanceMaxError(dist s1.ChordAngle) float64 {
	// This bound includes all source of error, assuming that the input points
	// are normalized. a and b are components of chord length that are
	// perpendicular and parallel to a plane containing the edge respectively.
	b := 0.5 * float64(dist) * float64(dist)
	a := float64(dist) * math.Sqrt(1-0.5*b)
	return ((2.5+2*math.Sqrt(3)+8.5*a)*a +
		(2+2*math.Sqrt(3)/3+6.5*(1-b))*b +
		(23+16/math.Sqrt(3))*dblEpsilon) * dblEpsilon
}

// DistanceFromSegment returns the distance of point X from line segment AB.
// The points are expected to be normalized.
func DistanceFromSegment(x, a, b Point) s1.Angle {
	var c s1.ChordAngle
	d, _ := updateMinDistance(x, a, b, c, true)
	return d.Angle()
}

// IsDistanceLess reports whether the distance from X to the edge AB is less
// than limit. This method is faster than DistanceFromSegment(). If you want to
// compare against a fixed s1.Angle, you should convert it to an s1.ChordAngle
// once and save the value, since this conversion is relatively expensive.
func IsDistanceLess(x, a, b Point, limit s1.ChordAngle) bool {
	_, less := MaybeUpdateMinDistance(x, a, b, limit)
	return less
}

// MaybeUpdateMinDistance checks if the distance from X to the edge AB is less
// then minDist, and if returns the possibly updated value and whether it was updated.
// The case A == B is handled correctly.
//
// Use this method when you want to compute many distances and keep track of
// the minimum. It is significantly faster than using DistanceFromSegment
// because (1) using s1.ChordAngle is much faster than s1.Angle, and (2) it
// can save a lot of work by not actually computing the distance when it is
// obviously larger than the current minimum.
func MaybeUpdateMinDistance(x, a, b Point, minDist s1.ChordAngle) (s1.ChordAngle, bool) {
	return updateMinDistance(x, a, b, minDist, false)
}

// IsInteriorDistanceLess reports whether the minimum distance from X to the
// edge AB is attained at an interior point of AB (i.e., not an endpoint), and
// that distance is less than limit.
func IsInteriorDistanceLess(x, a, b Point, limit s1.ChordAngle) bool {
	_, less := MaybeUpdateMinInteriorDistance(x, a, b, limit)
	return less
}

// MaybeUpdateMinInteriorDistance reports whether the minimum distance from X to AB
// is attained at an interior point of AB (i.e., not an endpoint), and that distance
// is less than minDist. If so, the value of minDist is updated and returns true.
// Otherwise it is unchanged and returns false.
func MaybeUpdateMinInteriorDistance(x, a, b Point, minDist s1.ChordAngle) (s1.ChordAngle, bool) {
	return interiorDist(x, a, b, minDist, false)
}

// updateMinDistance computes the distance from a point X to a line segment AB,
// and if either the distance was less than the given minDist, or alwaysUpdate is
// true, the value and whether it was updated are returned.
func updateMinDistance(x, a, b Point, minDist s1.ChordAngle, alwaysUpdate bool) (s1.ChordAngle, bool) {
	if d, ok := interiorDist(x, a, b, minDist, alwaysUpdate); ok {
		// Minimum distance is attained along the edge interior.
		return d, true
	}

	// Otherwise the minimum distance is to one of the endpoints.
	xa2, xb2 := (x.Sub(a.Vector)).Norm2(), x.Sub(b.Vector).Norm2()
	dist := s1.ChordAngle(math.Min(xa2, xb2))
	if !alwaysUpdate && dist >= minDist {
		return minDist, false
	}
	return dist, true
}

// interiorDist returns the shortest distance from point x to edge ab, assuming
// that the closest point to X is interior to AB. If the closest point is not
// interior to AB, interiorDist returns (minDist, false). If alwaysUpdate is set to
// false, the distance is only updated when the value exceeds certain the given minDist.
func interiorDist(x, a, b Point, minDist s1.ChordAngle, alwaysUpdate bool) (s1.ChordAngle, bool) {
	// Chord distance of x to both end points a and b.
	xa2, xb2 := (x.Sub(a.Vector)).Norm2(), x.Sub(b.Vector).Norm2()

	// The closest point on AB could either be one of the two vertices (the
	// vertex case) or in the interior (the interior case). Let C = A x B.
	// If X is in the spherical wedge extending from A to B around the axis
	// through C, then we are in the interior case. Otherwise we are in the
	// vertex case.
	//
	// Check whether we might be in the interior case. For this to be true, XAB
	// and XBA must both be acute angles. Checking this condition exactly is
	// expensive, so instead we consider the planar triangle ABX (which passes
	// through the sphere's interior). The planar angles XAB and XBA are always
	// less than the corresponding spherical angles, so if we are in the
	// interior case then both of these angles must be acute.
	//
	// We check this by computing the squared edge lengths of the planar
	// triangle ABX, and testing acuteness using the law of cosines:
	//
	//   max(XA^2, XB^2) < min(XA^2, XB^2) + AB^2
	if math.Max(xa2, xb2) >= math.Min(xa2, xb2)+(a.Sub(b.Vector)).Norm2() {
		return minDist, false
	}

	// The minimum distance might be to a point on the edge interior. Let R
	// be closest point to X that lies on the great circle through AB. Rather
	// than computing the geodesic distance along the surface of the sphere,
	// instead we compute the "chord length" through the sphere's interior.
	//
	// The squared chord length XR^2 can be expressed as XQ^2 + QR^2, where Q
	// is the point X projected onto the plane through the great circle AB.
	// The distance XQ^2 can be written as (X.C)^2 / |C|^2 where C = A x B.
	// We ignore the QR^2 term and instead use XQ^2 as a lower bound, since it
	// is faster and the corresponding distance on the Earth's surface is
	// accurate to within 1% for distances up to about 1800km.
	c := a.PointCross(b)
	c2 := c.Norm2()
	xDotC := x.Dot(c.Vector)
	xDotC2 := xDotC * xDotC
	if !alwaysUpdate && xDotC2 >= c2*float64(minDist) {
		// The closest point on the great circle AB is too far away.
		return minDist, false
	}

	// Otherwise we do the exact, more expensive test for the interior case.
	// This test is very likely to succeed because of the conservative planar
	// test we did initially.
	cx := c.Cross(x.Vector)
	if a.Dot(cx) >= 0 || b.Dot(cx) <= 0 {
		return minDist, false
	}

	// Compute the squared chord length XR^2 = XQ^2 + QR^2 (see above).
	// This calculation has good accuracy for all chord lengths since it
	// is based on both the dot product and cross product (rather than
	// deriving one from the other). However, note that the chord length
	// representation itself loses accuracy as the angle approaches π.
	qr := 1 - math.Sqrt(cx.Norm2()/c2)
	dist := s1.ChordAngle((xDotC2 / c2) + (qr * qr))

	if !alwaysUpdate && dist >= minDist {
		return minDist, false
	}

	return dist, true
}

// WedgeRel enumerates the possible relation between two wedges A and B.
type WedgeRel int

// Define the different possible relationships between two wedges.
const (
	WedgeEquals              WedgeRel = iota // A and B are equal.
	WedgeProperlyContains                    // A is a strict superset of B.
	WedgeIsProperlyContained                 // A is a strict subset of B.
	WedgeProperlyOverlaps                    // A-B, B-A, and A intersect B are non-empty.
	WedgeIsDisjoint                          // A and B are disjoint.
)

// WedgeRelation reports the relation between two non-empty wedges
// A=(a0, ab1, a2) and B=(b0, ab1, b2).
func WedgeRelation(a0, ab1, a2, b0, b2 Point) WedgeRel {
	// There are 6 possible edge orderings at a shared vertex (all
	// of these orderings are circular, i.e. abcd == bcda):
	//
	//  (1) a2 b2 b0 a0: A contains B
	//  (2) a2 a0 b0 b2: B contains A
	//  (3) a2 a0 b2 b0: A and B are disjoint
	//  (4) a2 b0 a0 b2: A and B intersect in one wedge
	//  (5) a2 b2 a0 b0: A and B intersect in one wedge
	//  (6) a2 b0 b2 a0: A and B intersect in two wedges
	//
	// We do not distinguish between 4, 5, and 6.
	// We pay extra attention when some of the edges overlap.  When edges
	// overlap, several of these orderings can be satisfied, and we take
	// the most specific.
	if a0 == b0 && a2 == b2 {
		return WedgeEquals
	}

	// Cases 1, 2, 5, and 6
	if OrderedCCW(a0, a2, b2, ab1) {
		// The cases with this vertex ordering are 1, 5, and 6,
		if OrderedCCW(b2, b0, a0, ab1) {
			return WedgeProperlyContains
		}

		// We are in case 5 or 6, or case 2 if a2 == b2.
		if a2 == b2 {
			return WedgeIsProperlyContained
		}
		return WedgeProperlyOverlaps

	}
	// We are in case 2, 3, or 4.
	if OrderedCCW(a0, b0, b2, ab1) {
		return WedgeIsProperlyContained
	}

	if OrderedCCW(a0, b0, a2, ab1) {
		return WedgeIsDisjoint
	}
	return WedgeProperlyOverlaps
}

// WedgeContains reports whether non-empty wedge A=(a0, ab1, a2) contains B=(b0, ab1, b2).
// Equivalent to WedgeRelation == WedgeProperlyContains || WedgeEquals.
func WedgeContains(a0, ab1, a2, b0, b2 Point) bool {
	// For A to contain B (where each loop interior is defined to be its left
	// side), the CCW edge order around ab1 must be a2 b2 b0 a0.  We split
	// this test into two parts that test three vertices each.
	return OrderedCCW(a2, b2, b0, ab1) && OrderedCCW(b0, a0, a2, ab1)
}

// WedgeIntersects reports whether non-empty wedge A=(a0, ab1, a2) intersects B=(b0, ab1, b2).
// Equivalent to WedgeRelation == WedgeIsDisjoint
func WedgeIntersects(a0, ab1, a2, b0, b2 Point) bool {
	// For A not to intersect B (where each loop interior is defined to be
	// its left side), the CCW edge order around ab1 must be a0 b2 b0 a2.
	// Note that it's important to write these conditions as negatives
	// (!OrderedCCW(a,b,c,o) rather than Ordered(c,b,a,o)) to get correct
	// results when two vertices are the same.
	return !(OrderedCCW(a0, b2, b0, ab1) && OrderedCCW(b0, a2, a0, ab1))
}

// TODO(roberts): Differences from C++
//  LongitudePruner
//  FaceSegments
//  PointFromExact
//  IntersectionExact
//  intersectionExactError