More projections, more tests

README.md

@@ -1,12 +1,13 @@
 # flatsphere
 
 [![GoReportCard](https://goreportcard.com/badge/owlpinetech/flatsphere)](https://goreportcard.com/report/github.com/owlpinetech/flatsphere)
+![test](https://github.com/owlpinetech/flatsphere/actions/workflows/go.yml/badge.svg)
 
-A [Go][] library for converting between spherical and planar coordinates.
+A [Go](https://go.dev/) library for converting between spherical and planar coordinates.
 
 ## Prerequisites
 
-- **[Go][]**: any one of the **three latest major** [releases][go-releases].
+- **[Go](https://go.dev/)**: any one of the **three latest major** [releases](https://go.dev/doc/devel/release).
 
 ## Install
 

analysis.go

@@ -0,0 +1,32 @@
+package flatsphere
+
+import "math"
+
+// Try to find a root of the given target function. Returns NaN if the process fails.
+func newtonsMethod(
+	initialGuess float64,
+	targetFunc func(float64) float64,
+	derivative func(float64) float64,
+	tolerance float64,
+	epsilon float64,
+	maxIterations int,
+) float64 {
+	currentGuess := initialGuess
+	for i := 0; i < maxIterations; i++ {
+		// get function/derivative values for current guess
+		y := targetFunc(currentGuess)
+		yPrime := derivative(currentGuess)
+
+		// if denominator is too small error out
+		if math.Abs(yPrime) < epsilon {
+			break
+		}
+
+		nextGuess := currentGuess - y/yPrime
+		if math.Abs(nextGuess-currentGuess) < tolerance {
+			return nextGuess
+		}
+		currentGuess = nextGuess
+	}
+	return math.NaN()
+}

analysis_test.go

@@ -0,0 +1,26 @@
+package flatsphere
+
+import (
+	"math"
+	"testing"
+)
+
+func FuzzNewtonsMethodSqrt(f *testing.F) {
+	f.Add(612.0)
+	f.Add(4.0)
+	f.Add(100.0)
+	f.Add(2.0)
+	f.Add(123.0)
+	f.Skip(0.0)
+	f.Fuzz(func(t *testing.T, sqrd float64) {
+		f := func(x float64) float64 { return x*x - sqrd }
+		d := func(x float64) float64 { return 2 * x }
+		approx := newtonsMethod(sqrd/2, f, d, 1e-6, 1e-12, 100)
+		native := math.Sqrt(sqrd)
+		if math.IsNaN(approx) && !math.IsNaN(native) {
+			t.Errorf("expected %e, got %e", native, approx)
+		} else if !math.IsNaN(approx) && !withinTolerance(approx, native, 0.00001) {
+			t.Errorf("expected %e, got %e", native, approx)
+		}
+	})
+}

cylindrical.go

@@ -18,7 +18,7 @@ func (m Mercator) Inverse(x float64, y float64) (lat float64, lon float64) {
 	return math.Atan(math.Sinh(y)), x
 }
 
-func (m Mercator) Bounds() Bounds {
+func (m Mercator) PlanarBounds() Bounds {
 	return NewRectangleBounds(2*math.Pi, 2*math.Pi)
 }
 
@@ -40,7 +40,7 @@ func (p PlateCarree) Inverse(x float64, y float64) (lat float64, lon float64) {
 	return y, x
 }
 
-func (p PlateCarree) Bounds() Bounds {
+func (p PlateCarree) PlanarBounds() Bounds {
 	return NewRectangleBounds(2*math.Pi, math.Pi)
 }
 
@@ -65,7 +65,7 @@ func (e Equirectangular) Inverse(x float64, y float64) (lat float64, lon float64
 	return y * math.Cos(e.Parallel), x
 }
 
-func (e Equirectangular) Bounds() Bounds {
+func (e Equirectangular) PlanarBounds() Bounds {
 	return NewRectangleBounds(2*math.Pi, math.Pi/math.Cos(e.Parallel))
 }
 
@@ -90,14 +90,14 @@ func (l CylindricalEqualArea) Parallel() float64 {
 }
 
 func (l CylindricalEqualArea) Project(lat float64, lon float64) (x float64, y float64) {
-	return lon, math.Sin(lat) * l.Bounds().YMax
+	return lon, math.Sin(lat) * l.PlanarBounds().YMax
 }
 
 func (l CylindricalEqualArea) Inverse(x float64, y float64) (lat float64, lon float64) {
-	return math.Asin(y / l.Bounds().YMax), x
+	return math.Asin(y / l.PlanarBounds().YMax), x
 }
 
-func (l CylindricalEqualArea) Bounds() Bounds {
+func (l CylindricalEqualArea) PlanarBounds() Bounds {
 	return NewRectangleBounds(2*math.Pi, 2/l.Stretch)
 }
 
@@ -157,7 +157,7 @@ func (g GallStereographic) Inverse(x float64, y float64) (lat float64, lon float
 	return 2 * math.Atan(y/(1+math.Sqrt(2))), x
 }
 
-func (g GallStereographic) Bounds() Bounds {
+func (g GallStereographic) PlanarBounds() Bounds {
 	return NewRectangleBounds(2*math.Pi, 1.5*math.Pi)
 }
 
@@ -177,7 +177,7 @@ func (m Miller) Inverse(x float64, y float64) (lat float64, lon float64) {
 	return math.Atan(math.Sinh(y*0.8)) / 0.8, x
 }
 
-func (m Miller) Bounds() Bounds {
+func (m Miller) PlanarBounds() Bounds {
 	return NewRectangleBounds(2*math.Pi, 2.5*math.Log(math.Tan(9*math.Pi/20)))
 }
 
@@ -198,6 +198,6 @@ func (c Central) Inverse(x float64, y float64) (lat float64, lon float64) {
 	return math.Atan(y), x
 }
 
-func (c Central) Bounds() Bounds {
+func (c Central) PlanarBounds() Bounds {
 	return NewRectangleBounds(2*math.Pi, 2*math.Pi)
 }

cylindrical_test.go

@@ -21,8 +21,8 @@ func TestEqualAreaStretch(t *testing.T) {
 	for _, tc := range testCases {
 		t.Run(tc.name, func(t *testing.T) {
 			proj := NewCylindricalEqualArea(tc.lat)
-			if !withinTolerance(proj.Stretch, tc.stretch, 0.000001) {
-				t.Errorf("expected streth factor %f, got %f", tc.stretch, proj.Stretch)
+			if !withinTolerance(proj.Stretch, tc.stretch, 0.00001) {
+				t.Errorf("expected stretch factor %e, got %e", tc.stretch, proj.Stretch)
 			}
 		})
 	}

healpix.go

@@ -6,54 +6,53 @@ import (
 
 // An equal area projection combining Lambert cylindrical with interrupted Collignon for the polar regions.
 // https://en.wikipedia.org/wiki/HEALPix
-type HEALPix struct{}
+// See also: "Mapping on the HEALPix Grid", https://arxiv.org/pdf/astro-ph/0412607.pdf
+type HEALPixStandard struct{}
 
-func NewHEALPix() HEALPix {
-	return HEALPix{}
+func NewHEALPixStandard() HEALPixStandard {
+	return HEALPixStandard{}
 }
 
-func (h HEALPix) Project(lat float64, lon float64) (x float64, y float64) {
-	colatitude := lat + math.Pi/2
-	z := math.Cos(colatitude)
+func (h HEALPixStandard) Project(lat float64, lon float64) (x float64, y float64) {
+	z := math.Sin(lat)
 	if math.Abs(z) <= 2.0/3.0 {
 		// equatiorial region
 		return lon, 3 * (math.Pi / 8) * z
 	} else {
 		// polar region
-		facetX := math.Mod(lon, math.Pi/2)
-		sigma := 2 - math.Sqrt(3*(1-math.Abs(z)))
-		if z < 0 {
-			sigma = -sigma
-		}
-		y := (math.Pi / 4) * sigma
-		x := lon - (math.Abs(sigma)-1)*(facetX-math.Pi/4)
+		sigma := math.Sqrt(3 * (1 - math.Abs(z)))
+		y := (math.Pi / 4) * (2 - sigma)
+		y = math.Copysign(y, lat)
+
+		facetX := (math.Pi / 4) * (2*math.Floor(2+(2*lon)/math.Pi) - 3)
+		x := facetX + sigma*(lon-facetX)
 		return x, y
 	}
 }
 
-func (h HEALPix) Inverse(x float64, y float64) (lat float64, lon float64) {
+func (h HEALPixStandard) Inverse(x float64, y float64) (float64, float64) {
 	absY := math.Abs(y)
-	//if absY >= math.Pi/2 {
-	//	panic(fmt.Sprintf("flatsphere: domain error in projection coordinate y dimension, %v too big", y))
-	//}
-
 	if absY <= math.Pi/4 {
 		// equatorial region
 		z := (8 / (3 * math.Pi)) * y
-		colat := math.Acos(z)
-		return math.Pi/2 - colat, x
-	} else {
+		lat := math.Asin(z)
+		return lat, x
+	} else if absY < math.Pi/2 {
 		// polar region
-		tt := math.Mod(x, math.Pi/2)
-		lng := x - ((absY-math.Pi/4)/(absY-math.Pi/2))*(tt-math.Pi/4)
-		zz := 2 - 4*absY/math.Pi
-		z := (1 - 1.0/3.0*(zz*zz)) * (y / absY)
-		colat := math.Acos(z)
-		return math.Pi/2 - colat, lng
+		sigma := 2 - math.Abs(4*y)/math.Pi
+		z := 1 - (sigma * sigma / 3)
+		lat := math.Asin(z)
+		lat = math.Copysign(lat, y)
+
+		facetX := (math.Pi / 4) * (2*math.Floor(2+(2*x)/math.Pi) - 3)
+		lon := facetX + (x-facetX)/sigma
+		return lat, lon
+	} else {
+		return math.Copysign(math.Pi/2, y), -math.Pi
 	}
 }
 
-func (h HEALPix) Bounds() Bounds {
+func (h HEALPixStandard) PlanarBounds() Bounds {
 	return Bounds{
 		XMin: 0,
 		XMax: 2 * math.Pi,

invertability_test.go

@@ -53,30 +53,41 @@ func FuzzSinusoidal(f *testing.F) {
 	projectInverseFuzz(f, NewSinusoidal())
 }
 
-func FuzzHEALPix(f *testing.F) {
-	projectInverseFuzz(f, NewHEALPix())
-}
+//func FuzzHEALPixStandard(f *testing.F) {
+//	projectInverseFuzz(f, NewHEALPixStandard())
+//}
+
+//func FuzzMollweide(f *testing.F) {
+//	projectInverseFuzz(f, NewMollweide())
+//}
 
 func withinTolerance(n1, n2, tolerance float64) bool {
 	if n1 == n2 {
 		return true
 	}
 	diff := math.Abs(n1 - n2)
-	if n2 == 0 {
-		return diff < tolerance
-	}
-	return (diff / math.Abs(n2)) < tolerance
+	return diff < tolerance
 }
 
 func projectInverseFuzz(f *testing.F, proj Projection) {
 	f.Add(0.0, 0.0)
+	f.Add(0.0, math.Pi)
 	f.Add(math.Pi/2, math.Pi/4)
+	f.Add(math.Pi/2, 0.0)
 	f.Add(-math.Pi/2, -math.Pi/4)
+	f.Add(math.Pi/2, math.Pi)
+	f.Add(66.0, 0.0)
 	f.Fuzz(func(t *testing.T, lat float64, lon float64) {
+		if lat > math.Pi/2 {
+			lat = math.Mod(lat, math.Pi/2)
+		}
+		if lon > math.Pi {
+			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 %f,%f, got %f,%f", lat, lon, rlat, rlon)
+			t.Errorf("expected %e,%e, got %e,%e", lat, lon, rlat, rlon)
 		}
 	})
 }

projection.go

@@ -1,10 +1,13 @@
 package flatsphere
 
+// A package of functionality for converting from spherical locations to planar coordinates, or from
+// planar coordinates into spherical locations. Contains information specifying some characteristics
+// of the planar space mapped to by the projection functions, which can differ between projections.
 type Projection interface {
 	// Convert a location on the sphere (in radians) to a coordinate on the plane.
 	Project(lat float64, lon float64) (x float64, y float64)
 	// Convert a coordinate on the plane to a location in radians on the sphere.
 	Inverse(x float64, y float64) (lat float64, lon float64)
 	// Retrieve the planar bounds of the projection.
-	Bounds() Bounds
+	PlanarBounds() Bounds
 }

pseudocylindrical.go

@@ -18,6 +18,86 @@ func (s Sinusoidal) Inverse(x float64, y float64) (lat float64, lon float64) {
 	return y, x / math.Cos(y)
 }
 
-func (s Sinusoidal) Bounds() Bounds {
+func (s Sinusoidal) PlanarBounds() Bounds {
+	return NewRectangleBounds(2*math.Pi, math.Pi)
+}
+
+// An equal-area pseudocylindrical map commonly used for maps of the celestial sphere.
+// https://en.wikipedia.org/wiki/Mollweide_projection
+type Mollweide struct{}
+
+func NewMollweide() Mollweide {
+	return Mollweide{}
+}
+
+func (m Mollweide) Project(lat float64, lon float64) (x float64, y float64) {
+	f := func(t float64) float64 { return 2*t + math.Sin(2*t) }
+	d := func(t float64) float64 { return 2 + 2*math.Cos(2*t) }
+	theta := newtonsMethod(math.Pi*math.Sin(lat), f, d, 1e-9, 1e-15, 125)
+	if math.IsNaN(theta) {
+		theta = math.Copysign(math.Pi/2, lat)
+	}
+	return lon / math.Pi * 2 * math.Cos(theta), math.Sin(theta)
+}
+
+func (m Mollweide) Inverse(x float64, y float64) (lat float64, lon float64) {
+	theta := math.Asin(y)
+	lat = math.Asin((2*theta + math.Sin(2*theta)) / math.Pi)
+	lon = x / math.Cos(theta) * math.Pi / 2
+	return lat, lon
+}
+
+func (m Mollweide) PlanarBounds() Bounds {
+	return NewRectangleBounds(2*math.Pi, math.Pi)
+}
+
+// An equal-area projection combining Sinusoidal and Mollweide at different hemispheres.
+// While most presentations of this projection use an interrupted form, this type is an
+// uninterrupted version of Homolosine.
+// https://en.wikipedia.org/wiki/Goode_homolosine_projection
+type Homolosine struct {
+	m     Mollweide
+	s     Sinusoidal
+	phiH  float64
+	scale float64
+	yH    float64
+}
+
+func NewHomolosine() Homolosine {
+	m := NewMollweide()
+	_, y := m.Project(0.71098, 0)
+	return Homolosine{
+		m,
+		NewSinusoidal(),
+		0.71098,
+		math.Sqrt2,
+		y * math.Sqrt2,
+	}
+}
+
+func (h Homolosine) Project(lat float64, lon float64) (x float64, y float64) {
+	if math.Abs(lat) <= h.phiH {
+		return h.s.Project(lat, lon)
+	} else {
+		x, y := h.m.Project(lat, lon)
+		if lat > 0 {
+			return x * h.scale, y*h.scale + h.phiH - h.yH
+		} else {
+			return x * h.scale, y*h.scale - h.phiH + h.yH
+		}
+	}
+}
+
+func (h Homolosine) Inverse(x float64, y float64) (lat float64, lon float64) {
+	if math.Abs(y) <= h.phiH {
+		return h.s.Inverse(x, y)
+	} else if y > 0 {
+		return h.m.Inverse(x/h.scale, (y-h.phiH+h.yH)/h.scale)
+	} else {
+		return h.m.Inverse(x/h.scale, (y+h.phiH-h.yH)/h.scale)
+	}
+}
+
+func (h Homolosine) PlanarBounds() Bounds {
 	return NewRectangleBounds(2*math.Pi, math.Pi)
 }