/*
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

// This file contains a collection of methods for:
//
//   (1) Robustly clipping geodesic edges to the faces of the S2 biunit cube
//       (see s2stuv), and
//
//   (2) Robustly clipping 2D edges against 2D rectangles.
//
// These functions can be used to efficiently find the set of CellIDs that
// are intersected by a geodesic edge (e.g., see CrossingEdgeQuery).

import (
	"math"

	"github.com/golang/geo/r1"
	"github.com/golang/geo/r2"
	"github.com/golang/geo/r3"
)

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
)

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

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

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

// 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, interpolateFloat64(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, interpolateFloat64(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
}

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

// TODO(roberts): Differences from C++:
// type FaceSegment
// FaceSegments
