Cassini projection, some great new tests (#8)

README.md

@@ -64,7 +64,7 @@ Determine how representative of reality a projection is at a point on the sphere
 
 A list of the predefined projections supported by the package.
 
-**Invertability**: all projections in this library have a well defined inverse function. However, only some of them reasonably satisfy the invertability property `x = f_1(f(x))` due to floating-point error in computations, edge-cases resulting in infinities or NaNs, or ambiguity at some type of critical point (commonly the poles or prime meridian). The table below checks of *invertible* for all projections which satisfy the `x = f_1(f(x))` property for all valid spherical locations to a reasonably fine degree of precision, referred to here as **everywhere floating-point invertible**. Oblique transforms of floating-point invertible standard projections do not necessarily share that property.
+**Invertability**: all projections in this library have a well defined inverse function. However, only some of them reasonably satisfy the invertability property `x = f_1(f(x))` due to floating-point error in computations, edge-cases resulting in infinities or NaNs, or ambiguity at some type of critical point (commonly the poles or prime meridian). The table below checks of *invertible* for all projections which satisfy the `x = f_1(f(x))` property for all valid spherical locations **except the poles** to a reasonably fine degree of precision, referred to here as **everywhere floating-point invertible**. Oblique transforms of floating-point invertible standard projections do not necessarily share that property.
 
 Efforts are ongoing to improve the coverage of this property to more projections where possible.
 
@@ -81,7 +81,7 @@ Efforts are ongoing to improve the coverage of this property to more projections
 |Miller|:white_check_mark:|
 |Central|:white_check_mark:|
 |Sinusoidal|:white_check_mark:|
-|HEALPix| |
+|HEALPix|:white_check_mark:|
 |Mollweide| |
 |Homolosine| |
 |Eckert IV| |
@@ -92,6 +92,7 @@ Efforts are ongoing to improve the coverage of this property to more projections
 |Orthographic| |
 |Robinson|:white_check_mark:|
 |Natural Earth|:white_check_mark:|
+|Cassini|:white_check_mark:|
 
 ## Credits
 

bounded_test.go

@@ -81,10 +81,6 @@ func FuzzLambertAzimuthalProjectBounded(f *testing.F) {
 	projectionBoundedFuzz(f, NewLambertAzimuthal())
 }
 
-//func FuzzTransverseMercatorProjectBounded(f *testing.F) {
-//	projectionBoundedFuzz(f, NewObliqueProjection(NewMercator(), 0, math.Pi/2, -math.Pi/2))
-//}
-
 //func FuzzGnomonicProjectBounded(f *testing.F) {
 //	projectionBoundedFuzz(f, NewGnomonic())
 //}
@@ -105,6 +101,14 @@ func FuzzEqualEarthProjectBounded(f *testing.F) {
 	projectionBoundedFuzz(f, NewEqualEarth())
 }
 
+func FuzzCassiniProjectBounded(f *testing.F) {
+	projectionBoundedFuzz(f, NewCassini())
+}
+
+func FuzzTransversePlateCarreeProjectBounded(f *testing.F) {
+	projectionBoundedFuzz(f, NewObliqueProjection(NewPlateCarree(), 0, math.Pi/2, -math.Pi/2))
+}
+
 func projectionBoundedFuzz(f *testing.F, proj Projection) {
 	f.Add(0.0, 0.0)
 	f.Add(0.0, math.Pi)

cylindrical.go

@@ -1,6 +1,8 @@
 package flatsphere
 
-import "math"
+import (
+	"math"
+)
 
 // An early standard cylindrical projection useful for navigation.
 // https://en.wikipedia.org/wiki/Mercator_projection
@@ -201,3 +203,27 @@ func (c Central) Inverse(x float64, y float64) (lat float64, lon float64) {
 func (c Central) PlanarBounds() Bounds {
 	return NewRectangleBounds(2*math.Pi, 2*math.Pi)
 }
+
+// A transverse version of the Plate–Carée projection, implemented directly for efficiency.
+// https://en.wikipedia.org/wiki/Cassini_projection
+type Cassini struct{}
+
+func NewCassini() Cassini {
+	return Cassini{}
+}
+
+func (c Cassini) Project(lat float64, lon float64) (float64, float64) {
+	x := math.Asin(math.Cos(lat) * math.Sin(lon))
+	y := math.Atan2(math.Sin(lat), math.Cos(lat)*math.Cos(lon))
+	return x, y
+}
+
+func (c Cassini) Inverse(x float64, y float64) (float64, float64) {
+	lat := math.Asin(math.Cos(x) * math.Sin(y))
+	lon := math.Atan2(math.Sin(x), math.Cos(x)*math.Cos(y))
+	return lat, lon
+}
+
+func (c Cassini) PlanarBounds() Bounds {
+	return NewRectangleBounds(math.Pi, 2*math.Pi)
+}

healpix.go

@@ -48,7 +48,7 @@ func (h HEALPixStandard) Inverse(x float64, y float64) (float64, float64) {
 		lon := facetX + (x-facetX)/sigma
 		return lat, lon
 	} else {
-		return math.Copysign(math.Pi/2, y), -math.Pi
+		return math.Copysign(math.Pi/2, y), x
 	}
 }
 

invertability_test.go

@@ -53,9 +53,9 @@ func FuzzSinusoidalProjectInverse(f *testing.F) {
 	projectInverseFuzz(f, NewSinusoidal())
 }
 
-//func FuzzHEALPixStandardProjectInverse(f *testing.F) {
-//	projectInverseFuzz(f, NewHEALPixStandard())
-//}
+func FuzzHEALPixStandardProjectInverse(f *testing.F) {
+	projectInverseFuzz(f, NewHEALPixStandard())
+}
 
 //func FuzzMollweideProjectInverse(f *testing.F) {
 //	projectInverseFuzz(f, NewMollweide())
@@ -97,6 +97,10 @@ func FuzzNaturalEarthProjectInverse(f *testing.F) {
 //	projectInverseFuzz(f, NewEqualEarth())
 //}
 
+func FuzzCassiniProjectInverse(f *testing.F) {
+	projectInverseFuzz(f, NewCassini())
+}
+
 func withinTolerance(n1, n2, tolerance float64) bool {
 	if n1 == n2 {
 		return true
@@ -114,16 +118,19 @@ func projectInverseFuzz(f *testing.F, proj Projection) {
 	f.Add(math.Pi/2, math.Pi)
 	f.Add(66.0, 0.0)
 	f.Fuzz(func(t *testing.T, lat float64, lon float64) {
-		if math.Abs(lat) > math.Pi/2 {
-			lat = math.Mod(lat, math.Pi/2)
-		}
-		if math.Abs(lon) > math.Pi {
-			lon = math.Mod(lon, math.Pi)
-		}
+		lat = math.Mod(lat, math.Pi/2)
+		lon = math.Mod(lon, math.Pi)
 		x, y := proj.Project(lat, lon)
 		rlat, rlon := proj.Inverse(x, y)
-		if !withinTolerance(lat, rlat, 0.00001) || !withinTolerance(lon, rlon, 0.00001) {
-			t.Errorf("expected %e,%e, got %e,%e", lat, lon, rlat, rlon)
+
+		if withinTolerance(lat, math.Pi/2, 0.0000001) || withinTolerance(lat, -math.Pi/2, 0.0000001) {
+			if !withinTolerance(lat, rlat, 0.000001) {
+				t.Errorf("expected lat %e, but got %e", lat, rlat)
+			}
+		} else {
+			if !withinTolerance(lat, rlat, 0.00001) || !withinTolerance(lon, rlon, 0.00001) {
+				t.Errorf("expected %e,%e, got %e,%e", lat, lon, rlat, rlon)
+			}
 		}
 	})
 }

oblique.go

@@ -114,8 +114,10 @@ func (o ObliqueProjection) TransformToOblique(latitude float64, longitude float6
 	newLat := math.Asin(preAsin)
 	var newLon float64
 	inner := math.Sin(latitude)/o.cosPoleLat/math.Cos(newLat) - math.Tan(o.poleLat)*math.Tan(newLat)
-	if o.poleLat == -math.Pi/2 {
-		newLon = o.poleLon + rotateLon
+	if o.poleLat == math.Pi/2 {
+		newLon = rotateLon + o.poleLon
+	} else if o.poleLat == -math.Pi/2 {
+		newLon = -rotateLon + o.poleLon + math.Pi
 	} else if math.Abs(inner) > 1 {
 		if (rotateLon == 0 && latitude < -o.poleLat) || (rotateLon != 0 && latitude < o.poleLat) {
 			newLon = o.poleLon + math.Pi

oblique_test.go

@@ -0,0 +1,63 @@
+package flatsphere
+
+import (
+	"math"
+	"testing"
+)
+
+func FuzzObliqueNoOp(f *testing.F) {
+	obliqEq := NewObliqueProjection(NewPlateCarree(), math.Pi/2, 0, 0)
+	cassini := NewPlateCarree()
+
+	f.Add(0.0, 0.0)
+	f.Add(math.Pi/4, math.Pi/4)
+	f.Add(math.Pi/2, 0.0)
+	f.Add(0.0, math.Pi/2)
+	f.Fuzz(func(t *testing.T, lat float64, lon float64) {
+		if math.Abs(lat) > math.Pi/2 {
+			lat = math.Mod(lat, math.Pi/2)
+		}
+		if math.Abs(lon) > math.Pi {
+			lon = math.Mod(lon, math.Pi)
+		}
+
+		xo, yo := obliqEq.Project(lat, lon)
+		xc, yc := cassini.Project(lat, lon)
+
+		if !withinTolerance(xo, xc, 0.000001) || !withinTolerance(yo, yc, 0.000001) {
+			t.Errorf("expected %e,%e, but got %e,%e", xc, yc, xo, yo)
+		}
+	})
+}
+
+func rotatePoint(x, y float64, rad float64) (float64, float64) {
+	rotX := x*math.Cos(rad) - y*math.Sin(rad)
+	rotY := x*math.Sin(rad) + y*math.Cos(rad)
+	return rotX, rotY
+}
+
+func FuzzObliquePlateCarreeVsCassiniProject(f *testing.F) {
+	obliqEq := NewObliqueProjection(NewPlateCarree(), 0, math.Pi/2, -math.Pi/2)
+	cassini := NewCassini()
+
+	f.Add(0.0, 0.0)
+	f.Add(math.Pi/4, math.Pi/4)
+	f.Add(math.Pi/2, 0.0)
+	f.Add(0.0, math.Pi/2)
+	f.Fuzz(func(t *testing.T, lat float64, lon float64) {
+		if math.Abs(lat) > math.Pi/2 {
+			lat = math.Mod(lat, math.Pi/2)
+		}
+		if math.Abs(lon) > math.Pi {
+			lon = math.Mod(lon, math.Pi)
+		}
+
+		xo, yo := obliqEq.Project(lat, lon)
+		xc, yc := cassini.Project(lat, lon)
+		xr, yr := rotatePoint(xc, yc, math.Pi/2)
+
+		if !withinTolerance(xo, xr, 0.000001) || !withinTolerance(yo, yr, 0.000001) {
+			t.Errorf("expected %e,%e, but got %e,%e", xr, yr, xo, yo)
+		}
+	})
+}

project_test.go

@@ -60,6 +60,33 @@ func TestMercatorInverseSanity(t *testing.T) {
 	})
 }
 
+func TestCassiniProjectSanity(t *testing.T) {
+	checkProject(t, "cassini", NewCassini(), []projectTestCase{
+		{0, 0, 0, 0},
+		{0, math.Pi, 0, math.Pi},
+		{0, math.Pi / 2, math.Pi / 2, 0},
+		{math.Pi / 4, 0, 0, math.Pi / 4},
+		{-math.Pi / 4, 0, 0, -math.Pi / 4},
+		{math.Pi / 2, 0, 0, math.Pi / 2},
+		{-math.Pi / 2, 0, 0, -math.Pi / 2},
+		{math.Pi / 2, -math.Pi, 0, math.Pi / 2},
+		{-math.Pi / 2, math.Pi, 0, -math.Pi / 2},
+	})
+}
+
+func TestCassiniInverseSanity(t *testing.T) {
+	checkInverse(t, "invCassini", NewCassini(), []inverseTestCase{
+		{0, 0, 0, 0},
+		{0, math.Pi, 0, math.Pi},
+		{math.Pi / 2, 0, 0, math.Pi / 2},
+		{0, math.Pi / 4, math.Pi / 4, 0},
+		{0, -math.Pi / 4, -math.Pi / 4, 0},
+		{0, math.Pi / 2, math.Pi / 2, 0},
+		{0, -math.Pi / 2, -math.Pi / 2, 0},
+		{math.Pi / 2, math.Pi, 0, math.Pi / 2},
+	})
+}
+
 /*func TestTransverseMercatorProjectSanity(t *testing.T) {
 	xFrom := func(lat, lon float64) float64 {
 		return math.Log((1+math.Sin(lon)*math.Cos(lat))/(1-math.Sin(lon)*math.Cos(lat))) / 2