WIP: splines

geom/float2.go

@@ -0,0 +1,63 @@
+package geom
+
+import "fmt"
+
+// Float2 is an element of float64^2
+type Float2 struct {
+	X, Y float64
+}
+
+// String returns a string representation of p like "(3.0,4.0,5.0)".
+func (p Float2) String() string {
+	return fmt.Sprintf("(%f,%f)", p.X, p.Y)
+}
+
+// Add performs vector addition.
+func (p Float2) Add(q Float2) Float2 {
+	return Float2{p.X + q.X, p.Y + q.Y}
+}
+
+// Sub performs vector subtraction.
+func (p Float2) Sub(q Float2) Float2 {
+	return p.Add(q.Neg())
+}
+
+// CMul performs componentwise multiplication.
+func (p Float2) CMul(q Float2) Float2 {
+	return Float2{p.X * q.X, p.Y * q.Y}
+}
+
+// Mul performs scalar multiplication.
+func (p Float2) Mul(k float64) Float2 {
+	return Float2{p.X * k, p.Y * k}
+}
+
+// CDiv performs componentwise division.
+func (p Float2) CDiv(q Float2) Float2 {
+	return Float2{p.X / q.X, p.Y / q.Y}
+}
+
+// Div performs scalar division by k.
+func (p Float2) Div(k float64) Float2 {
+	return Float2{p.X / k, p.Y / k}
+}
+
+// Neg returns the vector pointing in the opposite direction.
+func (p Float2) Neg() Float2 {
+	return Float2{-p.X, -p.Y}
+}
+
+// Coord returns the components of the vector.
+func (p Float2) Coord() (x, y float64) {
+	return p.X, p.Y
+}
+
+// Sign returns a sign vector.
+func (p Float2) Sign() Float2 {
+	return Float2{FSign(p.X), FSign(p.Y)}
+}
+
+// Dot returns the dot product of the two vectors.
+func (p Float2) Dot(q Float2) float64 {
+	return p.X*q.X + p.Y*q.Y
+}

geom/float3.go

@@ -0,0 +1,63 @@
+package geom
+
+import "fmt"
+
+// Float3 is an element of float64^3.
+type Float3 struct {
+	X, Y, Z float64
+}
+
+// String returns a string representation of p like "(3.0,4.0,5.0)".
+func (p Float3) String() string {
+	return fmt.Sprintf("(%f,%f,%f)", p.X, p.Y, p.Z)
+}
+
+// Add performs vector addition.
+func (p Float3) Add(q Float3) Float3 {
+	return Float3{p.X + q.X, p.Y + q.Y, p.Z + q.Z}
+}
+
+// Sub performs vector subtraction.
+func (p Float3) Sub(q Float3) Float3 {
+	return p.Add(q.Neg())
+}
+
+// CMul performs componentwise multiplication.
+func (p Float3) CMul(q Float3) Float3 {
+	return Float3{p.X * q.X, p.Y * q.Y, p.Z * q.Z}
+}
+
+// Mul performs scalar multiplication.
+func (p Float3) Mul(k float64) Float3 {
+	return Float3{p.X * k, p.Y * k, p.Z * k}
+}
+
+// CDiv performs componentwise division.
+func (p Float3) CDiv(q Float3) Float3 {
+	return Float3{p.X / q.X, p.Y / q.Y, p.Z / q.Z}
+}
+
+// Div performs scalar division by k.
+func (p Float3) Div(k float64) Float3 {
+	return Float3{p.X / k, p.Y / k, p.Z / k}
+}
+
+// Neg returns the vector pointing in the opposite direction.
+func (p Float3) Neg() Float3 {
+	return Float3{-p.X, -p.Y, -p.Z}
+}
+
+// Coord returns the components of the vector.
+func (p Float3) Coord() (x, y, z float64) {
+	return p.X, p.Y, p.Z
+}
+
+// Sign returns a sign vector.
+func (p Float3) Sign() Float3 {
+	return Float3{FSign(p.X), FSign(p.Y), FSign(p.Z)}
+}
+
+// Dot returns the dot product of the two vectors.
+func (p Float3) Dot(q Float3) float64 {
+	return p.X*q.X + p.Y*q.Y + p.Z*q.Z
+}

geom/int3.go

@@ -17,7 +17,6 @@ limitations under the License.
 package geom
 
 import (
-	"fmt"
 	"image"
 	"strconv"
 )
@@ -96,85 +95,3 @@ func (p Int3) Sign() Int3 {
 func (p Int3) Dot(q Int3) int {
 	return p.X*q.X + p.Y*q.Y + p.Z*q.Z
 }
-
-// Sign returns the sign of the int (-1, 0, or 1).
-func Sign(m int) int {
-	if m == 0 {
-		return 0
-	}
-	if m < 0 {
-		return -1
-	}
-	return 1
-}
-
-// FSign returns the sign of the float64 (-1, 0, or 1).
-func FSign(m float64) float64 {
-	if m == 0 {
-		return 0
-	}
-	if m < 0 {
-		return -1
-	}
-	return 1
-}
-
-// Float3 is an element of float64^3.
-type Float3 struct {
-	X, Y, Z float64
-}
-
-// String returns a string representation of p like "(3.0,4.0,5.0)".
-func (p Float3) String() string {
-	return fmt.Sprintf("(%f,%f,%f)", p.X, p.Y, p.Z)
-}
-
-// Add performs vector addition.
-func (p Float3) Add(q Float3) Float3 {
-	return Float3{p.X + q.X, p.Y + q.Y, p.Z + q.Z}
-}
-
-// Sub performs vector subtraction.
-func (p Float3) Sub(q Float3) Float3 {
-	return p.Add(q.Neg())
-}
-
-// CMul performs componentwise multiplication.
-func (p Float3) CMul(q Float3) Float3 {
-	return Float3{p.X * q.X, p.Y * q.Y, p.Z * q.Z}
-}
-
-// Mul performs scalar multiplication.
-func (p Float3) Mul(k float64) Float3 {
-	return Float3{p.X * k, p.Y * k, p.Z * k}
-}
-
-// CDiv performs componentwise division.
-func (p Float3) CDiv(q Float3) Float3 {
-	return Float3{p.X / q.X, p.Y / q.Y, p.Z / q.Z}
-}
-
-// Div performs scalar division by k.
-func (p Float3) Div(k float64) Float3 {
-	return Float3{p.X / k, p.Y / k, p.Z / k}
-}
-
-// Neg returns the vector pointing in the opposite direction.
-func (p Float3) Neg() Float3 {
-	return Float3{-p.X, -p.Y, -p.Z}
-}
-
-// Coord returns the components of the vector.
-func (p Float3) Coord() (x, y, z float64) {
-	return p.X, p.Y, p.Z
-}
-
-// Sign returns a sign vector.
-func (p Float3) Sign() Float3 {
-	return Float3{FSign(p.X), FSign(p.Y), FSign(p.Z)}
-}
-
-// Dot returns the dot product of the two vectors.
-func (p Float3) Dot(q Float3) float64 {
-	return p.X*q.X + p.Y*q.Y + p.Z*q.Z
-}

geom/misc.go

@@ -39,3 +39,27 @@ func CFloat(p image.Point) (x, y float64) {
 func Dot(p, q image.Point) int {
 	return p.X*q.X + p.Y*q.Y
 }
+
+// ---------- Some other helpers ----------
+
+// FSign returns the sign of the float64 (-1, 0, or 1).
+func FSign(m float64) float64 {
+	if m == 0 {
+		return 0
+	}
+	if m < 0 {
+		return -1
+	}
+	return 1
+}
+
+// Sign returns the sign of the int (-1, 0, or 1).
+func Sign(m int) int {
+	if m == 0 {
+		return 0
+	}
+	if m < 0 {
+		return -1
+	}
+	return 1
+}

geom/spline.go

@@ -0,0 +1,100 @@
+package geom
+
+import (
+	"errors"
+	"fmt"
+	"sort"
+)
+
+// LinearSpline implements a linear spline.
+type LinearSpline struct {
+	Points []Float2
+
+	deriv []float64 // slope of segment between points i and i+1
+}
+
+// Prepare sorts the points and computes internal information.
+func (s *LinearSpline) Prepare() error {
+	if len(s.Points) < 1 {
+		return errors.New("spline needs at least 1 point")
+	}
+
+	// Ensure Points is sorted.
+	sort.Slice(s.Points, func(i, j int) bool {
+		return s.Points[i].X < s.Points[j].X
+	})
+
+	// Check for points with equal X and compute derivatives
+	s.deriv = make([]float64, len(s.Points)-1)
+	for i := range s.Points[1:] {
+		if s.Points[i].X == s.Points[i+1].X {
+			return fmt.Errorf("spline value defined twice [%v, %v]", s.Points[i], s.Points[i+1])
+		}
+		s.deriv[i] = (s.Points[i+1].Y - s.Points[i].Y) / (s.Points[i+1].X - s.Points[i].X)
+	}
+	return nil
+}
+
+// Interpolate, given x, returns y where (x,y) is a point on the spline.
+// If x is outside the spline, it extrapolates from either the first or
+// last segments of the spline.
+func (s *LinearSpline) Interpolate(x float64) float64 {
+	N := len(s.Points)
+	if N == 1 {
+		return s.Points[0].Y
+	}
+	if x < s.Points[1].X {
+		// Comes before the end of the first segment
+		return s.Points[0].Y + (x-s.Points[0].X)*s.deriv[0]
+	}
+	if x > s.Points[N-2].X {
+		// Comes after the start of the last segment
+		return s.Points[N-1].Y + (x-s.Points[N-1].X)*s.deriv[N-2]
+	}
+	// Somewhere in the middle
+	i := sort.Search(N, func(n int) bool {
+		return x <= s.Points[n].X
+	})
+	if x == s.Points[i].X {
+		// Hit a point exactly
+		return s.Points[i].Y
+	}
+	// In the interval between point i and point i+1
+	return s.Points[i].Y + (x-s.Points[i].X)*s.deriv[i]
+}
+
+// CubicSpline implements a cubic spline.
+type CubicSpline struct {
+	// Points on the spline
+	Points []Float2
+
+	deriv, deriv2 []float64
+}
+
+// Prepare
+func (s *CubicSpline) Prepare() error {
+	if len(s.Points) < 1 {
+		return errors.New("spline needs at least 1 point")
+	}
+	// Ensure Points is sorted.
+	sort.Slice(s.Points, func(i, j int) bool {
+		return s.Points[i].X < s.Points[j].X
+	})
+	// Check for points with equal X.
+	for i := range s.Points[1:] {
+		if s.Points[i].X == s.Points[i+1].X {
+			return fmt.Errorf("spline value defined twice [%v, %v]", s.Points[i], s.Points[i+1])
+		}
+	}
+	// TODO: compute deriv and deriv2
+	return nil
+}
+
+func (s *CubicSpline) Interpolate(x float64) float64 {
+	N := len(s.Points)
+	if N == 1 {
+		return s.Points[0].Y
+	}
+	// TODO
+	return 0
+}

geom/spline_test.go

@@ -0,0 +1,79 @@
+package geom
+
+import "testing"
+
+func TestLinearSplineNoPoints(t *testing.T) {
+	s := &LinearSpline{}
+	if err := s.Prepare(); err == nil {
+		t.Errorf("s.Prepare() = %v, want error", err)
+	}
+}
+
+func TestLinearSplineEqualXPoints(t *testing.T) {
+	s := &LinearSpline{
+		Points: []Float2{{-5, 1}, {-2, 7}, {-2, -3}, {0, 2}, {3, -2}},
+	}
+	if err := s.Prepare(); err == nil {
+		t.Errorf("s.Prepare() = %v, want error", err)
+	}
+}
+
+func TestLinearSplineOnePoint(t *testing.T) {
+	s := &LinearSpline{
+		Points: []Float2{{-2, -3}},
+	}
+	if err := s.Prepare(); err != nil {
+		t.Errorf("s.Prepare() = %v, want nil", err)
+	}
+	for _, x := range []float64{-5, -4, -2, 0, 1, 7} {
+		if got, want := s.Interpolate(x), float64(-3); got != want {
+			t.Errorf("s.Interpolate(%v) = %v, want %v", x, got, want)
+		}
+	}
+}
+
+func TestLinearSpline(t *testing.T) {
+	s := &LinearSpline{
+		Points: []Float2{{-7, -2}, {-5, 1}, {-3, 0}, {-2, -3}, {0, 2}, {1, -5}, {3, -2}, {4, 4}},
+	}
+	if err := s.Prepare(); err != nil {
+		t.Errorf("s.Prepare() = %v, want nil", err)
+	}
+	tests := []struct {
+		x, want float64
+	}{
+		{x: -8, want: -3.5},
+		{x: -7.5, want: -2.75},
+		{x: -7, want: -2},
+		{x: -6.5, want: -1.25},
+		{x: -6, want: -0.5},
+		{x: -5.5, want: 0.25},
+		{x: -5, want: 1},
+		{x: -4.5, want: 4.5},
+		{x: -4, want: 3},
+		{x: -3.5, want: 1.5},
+		{x: -3, want: 0},
+		{x: -2.5, want: -4.25},
+		{x: -2, want: -3},
+		{x: -1.5, want: 12.5},
+		{x: -1, want: 9},
+		{x: -0.5, want: 5.5},
+		{x: 0, want: 2},
+		{x: 0.5, want: -5.75},
+		{x: 1, want: -5},
+		{x: 1.5, want: -11},
+		{x: 2, want: -8},
+		{x: 2.5, want: -5},
+		{x: 3, want: -2},
+		{x: 3.5, want: 1},
+		{x: 4, want: 4},
+		{x: 4.5, want: 7},
+		{x: 5, want: 10},
+		{x: 5.5, want: 13},
+	}
+	for _, test := range tests {
+		if got := s.Interpolate(test.x); got != test.want {
+			t.Errorf("s.Interpolate(%v) = %v, want %v", test.x, got, test.want)
+		}
+	}
+}