s2: Update EdgeTesselator to be more robust and conservative.

Previously the algorithm could not handle points of inflection (that is,
edges that cross the equator in the Plate Carree and Mercator
projections), and was not always conservative even when edges did not
have a point of inflection.

The new algorithm solves these problems by modeling the error as a cubic
polynomial. This is done by evaluating the distance between the edges at
two equally space points rather than at the midpoint. This results is a
small amount of extra tessellation (about 1.5% on average) but is much
more robust.

Signed-off-by: David Symonds <dsymonds@golang.org>

Author: rsned

s2/edge_tessellator.go

@@ -15,18 +15,162 @@
 package s2
 
 import (
-	"math"
-
 	"github.com/golang/geo/r2"
 	"github.com/golang/geo/s1"
 )
 
+// Tessellation is implemented by subdividing the edge until the estimated
+// maximum error is below the given tolerance. Estimating error is a hard
+// problem, especially when the only methods available are point evaluation of
+// the projection and its inverse. (These are the only methods that
+// Projection provides, which makes it easier and less error-prone to
+// implement new projections.)
+//
+// One technique that significantly increases robustness is to treat the
+// geodesic and projected edges as parametric curves rather than geometric ones.
+// Given a spherical edge AB and a projection p:S2->R2, let f(t) be the
+// normalized arc length parametrization of AB and let g(t) be the normalized
+// arc length parameterization of the projected edge p(A)p(B). (In other words,
+// f(0)=A, f(1)=B, g(0)=p(A), g(1)=p(B).)  We now define the geometric error as
+// the maximum distance from the point p^-1(g(t)) to the geodesic edge AB for
+// any t in [0,1], where p^-1 denotes the inverse projection. In other words,
+// the geometric error is the maximum distance from any point on the projected
+// edge (mapped back onto the sphere) to the geodesic edge AB. On the other
+// hand we define the parametric error as the maximum distance between the
+// points f(t) and p^-1(g(t)) for any t in [0,1], i.e. the maximum distance
+// (measured on the sphere) between the geodesic and projected points at the
+// same interpolation fraction t.
+//
+// The easiest way to estimate the parametric error is to simply evaluate both
+// edges at their midpoints and measure the distance between them (the "midpoint
+// method"). This is very fast and works quite well for most edges, however it
+// has one major drawback: it doesn't handle points of inflection (i.e., points
+// where the curvature changes sign). For example, edges in the Mercator and
+// Plate Carree projections always curve towards the equator relative to the
+// corresponding geodesic edge, so in these projections there is a point of
+// inflection whenever the projected edge crosses the equator. The worst case
+// occurs when the edge endpoints have different longitudes but the same
+// absolute latitude, since in that case the error is non-zero but the edges
+// have exactly the same midpoint (on the equator).
+//
+// One solution to this problem is to split the input edges at all inflection
+// points (i.e., along the equator in the case of the Mercator and Plate Carree
+// projections). However for general projections these inflection points can
+// occur anywhere on the sphere (e.g., consider the Transverse Mercator
+// projection). This could be addressed by adding methods to the S2Projection
+// interface to split edges at inflection points but this would make it harder
+// and more error-prone to implement new projections.
+//
+// Another problem with this approach is that the midpoint method sometimes
+// underestimates the true error even when edges do not cross the equator.
+// For the Plate Carree and Mercator projections, the midpoint method can
+// underestimate the error by up to 3%.
+//
+// Both of these problems can be solved as follows. We assume that the error
+// can be modeled as a convex combination of two worst-case functions, one
+// where the error is maximized at the edge midpoint and another where the
+// error is *minimized* (i.e., zero) at the edge midpoint. For example, we
+// could choose these functions as:
+//
+//    E1(x) = 1 - x^2
+//    E2(x) = x * (1 - x^2)
+//
+// where for convenience we use an interpolation parameter "x" in the range
+// [-1, 1] rather than the original "t" in the range [0, 1]. Note that both
+// error functions must have roots at x = {-1, 1} since the error must be zero
+// at the edge endpoints. E1 is simply a parabola whose maximum value is 1
+// attained at x = 0, while E2 is a cubic with an additional root at x = 0,
+// and whose maximum value is 2 * sqrt(3) / 9 attained at x = 1 / sqrt(3).
+//
+// Next, it is convenient to scale these functions so that the both have a
+// maximum value of 1. E1 already satisfies this requirement, and we simply
+// redefine E2 as
+//
+//   E2(x) = x * (1 - x^2) / (2 * sqrt(3) / 9)
+//
+// Now define x0 to be the point where these two functions intersect, i.e. the
+// point in the range (-1, 1) where E1(x0) = E2(x0). This value has the very
+// convenient property that if we evaluate the actual error E(x0), then the
+// maximum error on the entire interval [-1, 1] is bounded by
+//
+//   E(x) <= E(x0) / E1(x0)
+//
+// since whether the error is modeled using E1 or E2, the resulting function
+// has the same maximum value (namely E(x0) / E1(x0)). If it is modeled as
+// some other convex combination of E1 and E2, the maximum value can only
+// decrease.
+//
+// Finally, since E2 is not symmetric about the y-axis, we must also allow for
+// the possibility that the error is a convex combination of E1 and -E2. This
+// can be handled by evaluating the error at E(-x0) as well, and then
+// computing the final error bound as
+//
+//   E(x) <= max(E(x0), E(-x0)) / E1(x0) .
+//
+// Effectively, this method is simply evaluating the error at two points about
+// 1/3 and 2/3 of the way along the edges, and then scaling the maximum of
+// these two errors by a constant factor. Intuitively, the reason this works
+// is that if the two edges cross somewhere in the interior, then at least one
+// of these points will be far from the crossing.
+//
+// The actual algorithm implemented below has some additional refinements.
+// First, edges longer than 90 degrees are always subdivided; this avoids
+// various unusual situations that can happen with very long edges, and there
+// is really no reason to avoid adding vertices to edges that are so long.
+//
+// Second, the error function E1 above needs to be modified to take into
+// account spherical distortions. (It turns out that spherical distortions are
+// beneficial in the case of E2, i.e. they only make its error estimates
+// slightly more conservative.)  To do this, we model E1 as the maximum error
+// in a Plate Carree edge of length 90 degrees or less. This turns out to be
+// an edge from 45:-90 to 45:90 (in lat:lng format). The corresponding error
+// as a function of "x" in the range [-1, 1] can be computed as the distance
+// between the Plate Caree edge point (45, 90 * x) and the geodesic
+// edge point (90 - 45 * abs(x), 90 * sgn(x)). Using the Haversine formula,
+// the corresponding function E1 (normalized to have a maximum value of 1) is:
+//
+//   E1(x) =
+//     asin(sqrt(sin(Pi / 8 * (1 - x)) ^ 2 +
+//               sin(Pi / 4 * (1 - x)) ^ 2 * cos(Pi / 4) * sin(Pi / 4 * x))) /
+//     asin(sqrt((1 - 1 / sqrt(2)) / 2))
+//
+// Note that this function does not need to be evaluated at runtime, it
+// simply affects the calculation of the value x0 where E1(x0) = E2(x0)
+// and the corresponding scaling factor C = 1 / E1(x0).
+//
+// ------------------------------------------------------------------
+//
+// In the case of the Mercator and Plate Carree projections this strategy
+// produces a conservative upper bound (verified using 10 million random
+// edges). Furthermore the bound is nearly tight; the scaling constant is
+// C = 1.19289, whereas the maximum observed value was 1.19254.
+//
+// Compared to the simpler midpoint evaluation method, this strategy requires
+// more function evaluations (currently twice as many, but with a smarter
+// tessellation algorithm it will only be 50% more). It also results in a
+// small amount of additional tessellation (about 1.5%) compared to the
+// midpoint method, but this is due almost entirely to the fact that the
+// midpoint method does not yield conservative error estimates.
+//
+// For random edges with a tolerance of 1 meter, the expected amount of
+// overtessellation is as follows:
+//
+//                   Midpoint Method    Cubic Method
+//   Plate Carree               1.8%            3.0%
+//   Mercator                  15.8%           17.4%
+
 const (
-	// MinTessellationTolerance is the minimum supported tolerance (which
+	// tessellationInterpolationFraction is the fraction at which the two edges
+	// are evaluated in order to measure the error between them. (Edges are
+	// evaluated at two points measured this fraction from either end.)
+	tessellationInterpolationFraction = 0.31215691082248312
+	tessellationScaleFactor           = 0.83829992569888509
+
+	// minTessellationTolerance is the minimum supported tolerance (which
 	// corresponds to a distance less than 1 micrometer on the Earth's
 	// surface, but is still much larger than the expected projection and
 	// interpolation errors).
-	MinTessellationTolerance s1.Angle = 1e-13
+	minTessellationTolerance s1.Angle = 1e-13
 )
 
 // EdgeTessellator converts an edge in a given projection (e.g., Mercator) into
@@ -41,17 +185,18 @@ const (
 //   Projected   | S2 geodesics           | Planar projected edges
 //   Unprojected | Planar projected edges | S2 geodesics
 type EdgeTessellator struct {
-	projection   Projection
-	tolerance    s1.ChordAngle
-	wrapDistance r2.Point
+	projection Projection
+
+	// The given tolerance scaled by a constant fraction so that it can be
+	// compared against the result returned by estimateMaxError.
+	scaledTolerance s1.ChordAngle
 }
 
 // NewEdgeTessellator creates a new edge tessellator for the given projection and tolerance.
 func NewEdgeTessellator(p Projection, tolerance s1.Angle) *EdgeTessellator {
 	return &EdgeTessellator{
-		projection:   p,
-		tolerance:    s1.ChordAngleFromAngle(maxAngle(tolerance, MinTessellationTolerance)),
-		wrapDistance: p.WrapDistance(),
+		projection:      p,
+		scaledTolerance: s1.ChordAngleFromAngle(maxAngle(tolerance, minTessellationTolerance)),
 	}
 }
 
@@ -68,46 +213,25 @@ func (e *EdgeTessellator) AppendProjected(a, b Point, vertices []r2.Point) []r2.
 	if len(vertices) == 0 {
 		vertices = []r2.Point{pa}
 	} else {
-		pa = e.wrapDestination(vertices[len(vertices)-1], pa)
+		pa = e.projection.WrapDestination(vertices[len(vertices)-1], pa)
 	}
 
-	pb := e.wrapDestination(pa, e.projection.Project(b))
+	pb := e.projection.Project(b)
 	return e.appendProjected(pa, a, pb, b, vertices)
 }
 
 // appendProjected splits a geodesic edge AB as necessary and returns the
 // projected vertices appended to the given vertices.
 //
-// The maximum recursion depth is (math.Pi / MinTessellationTolerance) < 45
-func (e *EdgeTessellator) appendProjected(pa r2.Point, a Point, pb r2.Point, b Point, vertices []r2.Point) []r2.Point {
-	// It's impossible to robustly test whether a projected edge is close enough
-	// to a geodesic edge without knowing the details of the projection
-	// function, but the following heuristic works well for a wide range of map
-	// projections. The idea is simply to test whether the midpoint of the
-	// projected edge is close enough to the midpoint of the geodesic edge.
-	//
-	// This measures the distance between the two edges by treating them as
-	// parametric curves rather than geometric ones. The problem with
-	// measuring, say, the minimum distance from the projected midpoint to the
-	// geodesic edge is that this is a lower bound on the value we want, because
-	// the maximum separation between the two curves is generally not attained
-	// at the midpoint of the projected edge. The distance between the curve
-	// midpoints is at least an upper bound on the distance from either midpoint
-	// to opposite curve. It's not necessarily an upper bound on the maximum
-	// distance between the two curves, but it is a powerful requirement because
-	// it demands that the two curves stay parametrically close together. This
-	// turns out to be much more robust with respect for projections with
-	// singularities (e.g., the North and South poles in the rectangular and
-	// Mercator projections) because the curve parameterization speed changes
-	// rapidly near such singularities.
-	mid := Point{a.Add(b.Vector).Normalize()}
-	testMid := e.projection.Unproject(e.projection.Interpolate(0.5, pa, pb))
-
-	if ChordAngleBetweenPoints(mid, testMid) < e.tolerance {
+// The maximum recursion depth is (math.Pi / minTessellationTolerance) < 45
+func (e *EdgeTessellator) appendProjected(pa r2.Point, a Point, pbIn r2.Point, b Point, vertices []r2.Point) []r2.Point {
+	pb := e.projection.WrapDestination(pa, pbIn)
+	if e.estimateMaxError(pa, a, pb, b) <= e.scaledTolerance {
 		return append(vertices, pb)
 	}
 
-	pmid := e.wrapDestination(pa, e.projection.Project(mid))
+	mid := Point{a.Add(b.Vector).Normalize()}
+	pmid := e.projection.WrapDestination(pa, e.projection.Project(mid))
 	vertices = e.appendProjected(pa, a, pmid, mid, vertices)
 	return e.appendProjected(pmid, mid, pb, b, vertices)
 }
@@ -120,7 +244,6 @@ func (e *EdgeTessellator) appendProjected(pa r2.Point, a Point, pb r2.Point, b P
 // (e.g. across the 180 degree meridian) then the first and last vertices may not
 // be exactly the same.
 func (e *EdgeTessellator) AppendUnprojected(pa, pb r2.Point, vertices []Point) []Point {
-	pb2 := e.wrapDestination(pa, pb)
 	a := e.projection.Unproject(pa)
 	b := e.projection.Unproject(pb)
 
@@ -133,35 +256,36 @@ func (e *EdgeTessellator) AppendUnprojected(pa, pb r2.Point, vertices []Point) [
 	// transformed into "0:-175, 0:-181" while the second is transformed into
 	// "0:179, 0:183". The two coordinate pairs for the middle vertex
 	// ("0:-181" and "0:179") may not yield exactly the same S2Point.
-	return e.appendUnprojected(pa, a, pb2, b, vertices)
+	return e.appendUnprojected(pa, a, pb, b, vertices)
 }
 
 // appendUnprojected interpolates a projected edge and appends the corresponding
 // points on the sphere.
-func (e *EdgeTessellator) appendUnprojected(pa r2.Point, a Point, pb r2.Point, b Point, vertices []Point) []Point {
-	pmid := e.projection.Interpolate(0.5, pa, pb)
-	mid := e.projection.Unproject(pmid)
-	testMid := Point{a.Add(b.Vector).Normalize()}
-
-	if ChordAngleBetweenPoints(mid, testMid) < e.tolerance {
+func (e *EdgeTessellator) appendUnprojected(pa r2.Point, a Point, pbIn r2.Point, b Point, vertices []Point) []Point {
+	pb := e.projection.WrapDestination(pa, pbIn)
+	if e.estimateMaxError(pa, a, pb, b) <= e.scaledTolerance {
 		return append(vertices, b)
 	}
 
+	pmid := e.projection.Interpolate(0.5, pa, pb)
+	mid := e.projection.Unproject(pmid)
+
 	vertices = e.appendUnprojected(pa, a, pmid, mid, vertices)
 	return e.appendUnprojected(pmid, mid, pb, b, vertices)
 }
 
-// wrapDestination returns the coordinates of the edge destination wrapped if
-// necessary to obtain the shortest edge.
-func (e *EdgeTessellator) wrapDestination(pa, pb r2.Point) r2.Point {
-	x := pb.X
-	y := pb.Y
-	// The code below ensures that pb is unmodified unless wrapping is required.
-	if e.wrapDistance.X > 0 && math.Abs(x-pa.X) > 0.5*e.wrapDistance.X {
-		x = pa.X + math.Remainder(x-pa.X, e.wrapDistance.X)
-	}
-	if e.wrapDistance.Y > 0 && math.Abs(y-pa.Y) > 0.5*e.wrapDistance.Y {
-		y = pa.Y + math.Remainder(y-pa.Y, e.wrapDistance.Y)
+func (e *EdgeTessellator) estimateMaxError(pa r2.Point, a Point, pb r2.Point, b Point) s1.ChordAngle {
+	// See the algorithm description at the top of this file.
+	// We always tessellate edges longer than 90 degrees on the sphere, since the
+	// approximation below is not robust enough to handle such edges.
+	if a.Dot(b.Vector) < -1e-14 {
+		return s1.InfChordAngle()
 	}
-	return r2.Point{x, y}
+	t1 := tessellationInterpolationFraction
+	t2 := 1 - tessellationInterpolationFraction
+	mid1 := Interpolate(t1, a, b)
+	mid2 := Interpolate(t2, a, b)
+	pmid1 := e.projection.Unproject(e.projection.Interpolate(t1, pa, pb))
+	pmid2 := e.projection.Unproject(e.projection.Interpolate(t2, pa, pb))
+	return maxChordAngle(ChordAngleBetweenPoints(mid1, pmid1), ChordAngleBetweenPoints(mid2, pmid2))
 }

s2/projections.go

@@ -59,6 +59,16 @@ type Projection interface {
 	//
 	// If a given axis does not wrap, its WrapDistance should be set to zero.
 	WrapDistance() r2.Point
+
+	// WrapDestination that wraps the coordinates of B if necessary in order to
+	// obtain the shortest edge AB. For example, suppose that A = [170, 20],
+	// B = [-170, 20], and the projection wraps so that [x, y] == [x + 360, y].
+	// Then this function would return [190, 20] for point B (reducing the edge
+	// length in the "x" direction from 340 to 20).
+	WrapDestination(a, b r2.Point) r2.Point
+
+	// We do not support implementations of this interface outside this package.
+	privateInterface()
 }
 
 // PlateCarreeProjection defines the "plate carree" (square plate) projection,
@@ -127,6 +137,13 @@ func (p *PlateCarreeProjection) WrapDistance() r2.Point {
 	return r2.Point{p.xWrap, 0}
 }
 
+// WrapDestination wraps the points if needed to get the shortest edge.
+func (p *PlateCarreeProjection) WrapDestination(a, b r2.Point) r2.Point {
+	return wrapDestination(a, b, p.WrapDistance)
+}
+
+func (p *PlateCarreeProjection) privateInterface() {}
+
 // MercatorProjection defines the spherical Mercator projection. Google Maps
 // uses this projection together with WGS84 coordinates, in which case it is
 // known as the "Web Mercator" projection (see Wikipedia). This class makes
@@ -201,3 +218,24 @@ func (p *MercatorProjection) Interpolate(f float64, a, b r2.Point) r2.Point {
 func (p *MercatorProjection) WrapDistance() r2.Point {
 	return r2.Point{p.xWrap, 0}
 }
+
+// WrapDestination wraps the points if needed to get the shortest edge.
+func (p *MercatorProjection) WrapDestination(a, b r2.Point) r2.Point {
+	return wrapDestination(a, b, p.WrapDistance)
+}
+
+func (p *MercatorProjection) privateInterface() {}
+
+func wrapDestination(a, b r2.Point, wrapDistance func() r2.Point) r2.Point {
+	wrap := wrapDistance()
+	x := b.X
+	y := b.Y
+	// The code below ensures that "b" is unmodified unless wrapping is required.
+	if wrap.X > 0 && math.Abs(x-a.X) > 0.5*wrap.X {
+		x = a.X + math.Remainder(x-a.X, wrap.X)
+	}
+	if wrap.Y > 0 && math.Abs(y-a.Y) > 0.5*wrap.Y {
+		y = a.Y + math.Remainder(y-a.Y, wrap.Y)
+	}
+	return r2.Point{x, y}
+}

s2/projections_test.go

@@ -57,10 +57,6 @@ func TestPlateCarreeProjectionInterpolate(t *testing.T) {
 	}
 }
 
-func r2pointsApproxEqual(a, b r2.Point) bool {
-	return float64Eq(a.X, b.X) && float64Eq(a.Y, b.Y)
-}
-
 func TestPlateCarreeProjectionProjectUnproject(t *testing.T) {
 	tests := []struct {
 		have Point
@@ -77,7 +73,7 @@ func TestPlateCarreeProjectionProjectUnproject(t *testing.T) {
 	proj := NewPlateCarreeProjection(180)
 
 	for _, test := range tests {
-		if got := proj.Project(test.have); !r2pointsApproxEqual(test.want, got) {
+		if got := proj.Project(test.have); !r2PointsApproxEqual(test.want, got, epsilon) {
 			t.Errorf("proj.Project(%v) = %v, want %v", test.have, got, test.want)
 		}
 		if got := proj.Unproject(test.want); !got.ApproxEqual(test.have) {
@@ -102,7 +98,7 @@ func TestMercatorProjectionProjectUnproject(t *testing.T) {
 	proj := NewMercatorProjection(180)
 
 	for _, test := range tests {
-		if got := proj.Project(test.have); !r2pointsApproxEqual(test.want, got) {
+		if got := proj.Project(test.have); !r2PointsApproxEqual(test.want, got, epsilon) {
 			t.Errorf("proj.Project(%v) = %v, want %v", test.have, got, test.want)
 		}
 		if got := proj.Unproject(test.want); !got.ApproxEqual(test.have) {