Add methods for testing Points normalizability. (#127)

Add IsNormalizable and EnsureNormalizable and some basic tests for both.

Author: rsned

s2/point.go

@@ -317,3 +317,48 @@ func Rotate(p, axis Point, angle s1.Angle) Point {
 func (p Point) stableAngle(o Point) s1.Angle {
 	return s1.Angle(2 * math.Atan2(p.Sub(o.Vector).Norm(), p.Add(o.Vector).Norm()))
 }
+
+// IsNormalizable reports if the given Point's magnitude is large enough such that the
+// angle to another vector of the same magnitude can be measured using Angle()
+// without loss of precision due to floating-point underflow.  (This requirement
+// is also sufficient to ensure that Normalize() can be called without risk of
+// precision loss.)
+func (p Point) IsNormalizable() bool {
+	// Let ab = RobustCrossProd(a, b) and cd = RobustCrossProd(cd).  In order for
+	// ab.Angle(cd) to not lose precision, the squared magnitudes of ab and cd
+	// must each be at least 2**-484.  This ensures that the sum of the squared
+	// magnitudes of ab.CrossProd(cd) and ab.DotProd(cd) is at least 2**-968,
+	// which ensures that any denormalized terms in these two calculations do
+	// not affect the accuracy of the result (since all denormalized numbers are
+	// smaller than 2**-1022, which is less than dblError * 2**-968).
+	//
+	// The fastest way to ensure this is to test whether the largest component of
+	// the result has a magnitude of at least 2**-242.
+	return maxFloat64(math.Abs(p.X), math.Abs(p.Y), math.Abs(p.Z)) >= math.Ldexp(1, -242)
+}
+
+// EnsureNormalizable scales a vector as necessary to ensure that the result can
+// be normalized without loss of precision due to floating-point underflow.
+//
+// This requires p != (0, 0, 0)
+func (p Point) EnsureNormalizable() Point {
+	// TODO(rsned): Zero vector isn't normalizable, and we don't have DCHECK in Go.
+	// What is the appropriate return value in this case? Is it {NaN, NaN, NaN}?
+	if p == (Point{r3.Vector{0, 0, 0}}) {
+		return p
+	}
+	if !p.IsNormalizable() {
+		// We can't just scale by a fixed factor because the smallest representable
+		// double is 2**-1074, so if we multiplied by 2**(1074 - 242) then the
+		// result might be so large that we couldn't square it without overflow.
+		//
+		// Note that we must scale by a power of two to avoid rounding errors.
+		// The code below scales "p" such that the largest component is
+		// in the range [1, 2).
+		pMax := maxFloat64(math.Abs(p.X), math.Abs(p.Y), math.Abs(p.Z))
+
+		// This avoids signed overflow for any value of Ilogb().
+		return Point{p.Mul(math.Ldexp(2, -1-math.Ilogb(pMax)))}
+	}
+	return p
+}

s2/point_test.go

@@ -22,6 +22,12 @@ import (
 	"github.com/golang/geo/s1"
 )
 
+// pointNear reports if each component of the two points is within the given epsilon.
+// This is similar to Point/Vector.ApproxEqual but with a user supplied epsilon.
+func pointNear(a, b Point, ε float64) bool {
+	return math.Abs(a.X-b.X) < ε && math.Abs(a.Y-b.Y) < ε && math.Abs(a.Z-b.Z) < ε
+}
+
 func TestOriginPoint(t *testing.T) {
 	if math.Abs(OriginPoint().Norm()-1) > 1e-15 {
 		t.Errorf("Origin point norm = %v, want 1", OriginPoint().Norm())
@@ -341,3 +347,117 @@ func TestPointRotate(t *testing.T) {
 		}
 	}
 }
+
+func TestPointIsNormalizable(t *testing.T) {
+	tests := []struct {
+		have Point
+		want bool
+	}{
+		{
+			// 0,0,0 is not normalizeable.
+			have: Point{r3.Vector{X: 0, Y: 0, Z: 0}},
+			want: false,
+		},
+		{
+			have: Point{r3.Vector{X: 1, Y: 1, Z: 1}},
+			want: true,
+		},
+
+		// The approximate cutoff is ~1.4149498560666738e-73
+		{
+			have: Point{r3.Vector{X: 1, Y: 0, Z: 0}},
+			want: true,
+		},
+		{
+			// Only one too small component.
+			have: Point{r3.Vector{X: 1e-75, Y: 1, Z: 1}},
+			want: true,
+		},
+		{
+			// All three components exact boundary case.
+			have: Point{r3.Vector{
+				X: math.Ldexp(1, -242),
+				Y: math.Ldexp(1, -242),
+				Z: math.Ldexp(1, -242)}},
+			want: true,
+		},
+		{
+			// All three components too small.
+			have: Point{r3.Vector{
+				X: math.Ldexp(1, -243),
+				Y: math.Ldexp(1, -243),
+				Z: math.Ldexp(1, -243)}},
+			want: false,
+		},
+	}
+
+	for _, test := range tests {
+		if got := test.have.IsNormalizable(); got != test.want {
+			t.Errorf("%+v.IsNormalizable() = %t, want %t", test.have, got, test.want)
+		}
+	}
+}
+
+func TestPointEnsureNormalizable(t *testing.T) {
+	tests := []struct {
+		have Point
+		want Point
+	}{
+		{
+			// 0,0,0 is not normalizeable.
+			have: Point{r3.Vector{X: 0, Y: 0, Z: 0}},
+			want: Point{r3.Vector{X: 0, Y: 0, Z: 0}},
+		},
+		{
+			have: Point{r3.Vector{X: 1, Y: 0, Z: 0}},
+			want: Point{r3.Vector{X: 1, Y: 0, Z: 0}},
+		},
+		{
+			// All three components exact border for still normalizeable.
+			have: Point{r3.Vector{
+				X: math.Ldexp(1, -242),
+				Y: math.Ldexp(1, -242),
+				Z: math.Ldexp(1, -242),
+			}},
+			want: Point{r3.Vector{
+				X: math.Ldexp(1, -242),
+				Y: math.Ldexp(1, -242),
+				Z: math.Ldexp(1, -242),
+			}},
+		},
+		{
+			// All three components too small but the same.
+			have: Point{r3.Vector{
+				X: math.Ldexp(1, -243),
+				Y: math.Ldexp(1, -243),
+				Z: math.Ldexp(1, -243),
+			}},
+			want: Point{r3.Vector{
+				X: 1,
+				Y: 1,
+				Z: 1,
+			}},
+		},
+		{
+			// All three components too small but different.
+			have: Point{r3.Vector{
+				X: math.Ldexp(1, -243),
+				Y: math.Ldexp(1, -486),
+				Z: math.Ldexp(1, -729),
+			}},
+			want: Point{r3.Vector{
+				X: 1,
+				Y: 0,
+				Z: 0,
+			}},
+		},
+	}
+
+	for _, test := range tests {
+		got := test.have.EnsureNormalizable()
+		if !pointNear(got, test.want, 1e-50) {
+			t.Errorf("%+v.EnsureNormalizable() = %+v, want %+v",
+				test.have, got, test.want)
+		}
+	}
+}