Initial commit.

internal/helper/clamp.go

@@ -0,0 +1,63 @@
+package helper
+
+import (
+	"math"
+)
+
+// 0 if NaN or -Inf, 1 otherwise.
+func oneIfFinite(x float64) float64 {
+	if x > math.Inf(-1) {
+		return 1
+	}
+
+	// NaN or -Inf
+	return 0
+}
+
+// x if >= 0, 0 otherwise (including NaN).
+func zeroOrMore(x float64) float64 {
+	if x >= 0 {
+		return x
+	}
+
+	return 0
+}
+
+// the builtin max function doesn't ignore NaNs like I'd prefer, so do it ourselves.
+func max3(a, b, c float64) float64 {
+	result := a
+
+	if result != result || b > result {
+		result = b
+	}
+
+	if result != result || c > result {
+		result = c
+	}
+
+	return result
+}
+
+func ClampRGB(r, g, b float64) (_, _, _ float64) {
+	// if any components are greater than 1, scale them down back into a legal range and
+	// fade to white based to how for out of range they are.
+	if m := max3(r, g, b); m > 1 {
+		m2 := m * m
+
+		if math.IsInf(m2, 1) {
+			// This would be white if all components were sensible finite values,
+			// although we will return zeros for any that were NaN or -Inf.
+			return oneIfFinite(r), oneIfFinite(g), oneIfFinite(b)
+		}
+
+		c := 1 - 1/m
+		r = c + r/m2
+		g = c + g/m2
+		b = c + b/m2
+	}
+
+	// make sure no components are NaN or less than zero.
+	// note that we do this last so that the fade to white logic has a chance to bring
+	// components back into legal ranges.
+	return zeroOrMore(r), zeroOrMore(g), zeroOrMore(b)
+}

internal/helper/clamp_test.go

@@ -0,0 +1,75 @@
+package helper
+
+import (
+	"math"
+	"testing"
+)
+
+func TestClampRGB(t *testing.T) {
+	t.Run("values in legal ranges should be unmodified", func(t *testing.T) {
+		const steps = 8
+
+		for rI := range steps {
+			for gI := range steps {
+				for bI := range steps {
+					want := [3]float64{
+						float64(rI) * (1 / (steps - 1)),
+						float64(gI) * (1 / (steps - 1)),
+						float64(bI) * (1 / (steps - 1)),
+					}
+
+					if got := collect3(ClampRGB(want[0], want[1], want[2])); got != want {
+						t.Errorf("Clamp(%v) = %v, expected values to be unmodified", want, got)
+						return
+					}
+				}
+			}
+		}
+	})
+
+	tests := []struct {
+		name   string
+		values [3]float64
+		want   [3]float64
+	}{
+
+		{
+			// any component being infinity should result in white.
+			name:   "+inf",
+			values: [3]float64{-1, .5, math.Inf(1)},
+			want:   [3]float64{1, 1, 1},
+		}, {
+			// ... except the case where any other component was NaN or -infinity.
+			name:   "+inf, -inf, NaN",
+			values: [3]float64{math.Inf(-1), math.NaN(), math.Inf(1)},
+			want:   [3]float64{0, 0, 1},
+		}, {
+			// colors that are too bright should be scaled back to the legal range, and then
+			// interpolate to white by 1-1/max_value.
+			name:   "normalize over-bright colours and fade them to white",
+			values: [3]float64{1, 2, 3},
+			want:   [3]float64{1./3*(1./3) + (1 - 1./3), 2./3*(1./3) + (1 - 1./3), 1}, // note that
+		}, {
+			name:   "negative values should be clamped to 0",
+			values: [3]float64{-1, math.Inf(-1), .5},
+			want:   [3]float64{0, 0, .5},
+		}, {
+			name:   "except where the logic for over-bright colors would bring them back to the legal range",
+			values: [3]float64{-1, 0, 3},
+			want:   [3]float64{-1./3*(1./3) + (1 - 1./3), 1 - 1./3, 1},
+		},
+	}
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			for _, order := range permuteOrder3 {
+				values := permute3(tt.values, order)
+				want := permute3(tt.want, order)
+
+				if got := collect3(ClampRGB(values[0], values[1], values[2])); !EqFloat64SliceFuzzy(got[:], want[:]) {
+					t.Errorf("ClampRGB(%v) = %v, want %v", values, got, want)
+					return
+				}
+			}
+		})
+	}
+}

internal/helper/cmp.go

@@ -0,0 +1,52 @@
+package helper
+
+import (
+	"math"
+)
+
+// EqFloat64Fuzzy returns true if two floats aren't meaningfully distinct from each other.
+//
+// NaNs aren't considered distinct (meaning this function will return true if both inputs are NaN).
+func EqFloat64Fuzzy(a, b float64) bool {
+	// if either input is NaN...
+	if math.IsNaN(a) || math.IsNaN(b) {
+		// return true if they'be both NaN (think SQL's "IS NOT DISTINCT FROM")
+		// otherwise was was NaN and the other was not, so return false.
+		return math.IsNaN(a) == math.IsNaN(b)
+	}
+
+	// if either input is infinity...
+	if math.IsInf(a, 0) || math.IsInf(b, 0) {
+		// return true if they're the same value (both +infinity or both -infinity)
+		// false otherwise (infinity vs a finite number, or an infinity with the opposite sign)
+		return a == b
+	}
+
+	const epsilon = 1e-9
+
+	absA, absB, absDiff := math.Abs(a), math.Abs(b), math.Abs(a-b)
+
+	// For numbers close to zero, use absolute epsilon
+	if min(absA, absB, absDiff) < math.SmallestNonzeroFloat64 {
+		return absDiff < epsilon
+	}
+
+	return absDiff < epsilon*max(absA, absB)
+}
+
+// EqFloat64SliceFuzzy returns true if two lists of floats aren't meaningfully distinct from each other.
+//
+// Returns false if the lists are of different lengths, [EqFloat64Fuzzy] returns false for any pair of floats.
+func EqFloat64SliceFuzzy(a, b []float64) bool {
+	if len(a) != len(b) {
+		return false
+	}
+
+	for i := range a {
+		if !EqFloat64Fuzzy(a[i], b[i]) {
+			return false
+		}
+	}
+
+	return true
+}

internal/helper/cmp_test.go

@@ -0,0 +1,85 @@
+package helper
+
+import (
+	"math"
+	"testing"
+)
+
+func TestEqFloat64Fuzzy(t *testing.T) {
+	tests := []struct {
+		name string
+		a, b float64
+		want bool
+	}{
+		{"exactly equal", 1, 1, true},
+		{"nearly equal", 1, math.Nextafter(1, math.Inf(1)), true},
+		{"zero equal to itself", +0., -0., true},
+		{"zero not equal to non-zero", 0., 1e-9, false},
+		{"definitely not equal", 1, 1 + 1e-9, false},
+		{"infinity equal to itself", math.Inf(1), math.Inf(1), true},
+		{"infinity not equal to a finite value", math.Inf(1), 1, false},
+		{"NaN equal to itself", math.NaN(), math.NaN(), true},
+		{"NaN not equal to a finite value", math.NaN(), 1, false},
+		{"NaN not equal to infinity", math.NaN(), math.Inf(1), false},
+
+		// these are actual numbers encountered that should be equal, but failed the test at some point.
+		// keeping them around as test cases.
+		{"testcase1", 0.0015172579307272023, 0.001517257930727202, true},
+		{"testcase2", 0, -8.131516293641283e-20, true},
+		{"testcase3", 0, -5.0186702124817295e-20, true},
+		{"testcase4", 0, -6.776263578034403e-21, true},
+	}
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			if got := EqFloat64Fuzzy(tt.a, tt.b); got != tt.want {
+				t.Errorf("EqFloat64Fuzzy(%v, %v) = %v, want %v", tt.a, tt.b, got, tt.want)
+			}
+
+			// swapping the arguments shouldn't change the outcome.
+			if got := EqFloat64Fuzzy(tt.b, tt.a); got != tt.want {
+				t.Errorf("EqFloat64Fuzzy(%v, %v) = %v, want %v", tt.b, tt.a, got, tt.want)
+			}
+
+			// negating the arguments shouldn't change the outcome either
+			if got := EqFloat64Fuzzy(-tt.a, -tt.b); got != tt.want {
+				t.Errorf("EqFloat64Fuzzy(%v, %v) = %v, want %v", -tt.a, -tt.b, got, tt.want)
+			}
+
+			if got := EqFloat64Fuzzy(tt.b, tt.a); got != tt.want {
+				t.Errorf("EqFloat64Fuzzy(%v, %v) = %v, want %v", -tt.b, -tt.a, got, tt.want)
+			}
+		})
+	}
+}
+
+func TestEqFloat64SliceFuzzy(t *testing.T) {
+	tests := []struct {
+		name string
+		a, b []float64
+		want bool
+	}{
+		{
+			"equivalent float slices",
+			[]float64{1, 2, 3, math.NaN()},
+			[]float64{1, 2, math.Nextafter(3, math.Inf(1)), math.NaN()},
+			true,
+		}, {
+			"dissimilar float slices",
+			[]float64{1, 2, 4},
+			[]float64{1, 3, 4},
+			false,
+		}, {
+			"different lengths",
+			[]float64{1, 2},
+			[]float64{1, 2, 3},
+			false,
+		},
+	}
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			if got := EqFloat64SliceFuzzy(tt.a, tt.b); got != tt.want {
+				t.Errorf("EqFloat64SliceFuzzy(%v, %v) = %v, want %v", tt.a, tt.b, got, tt.want)
+			}
+		})
+	}
+}

internal/helper/collect.go

@@ -0,0 +1,33 @@
+package helper
+
+// we have several places where a function returns multiple values,
+// and collecting them into an array so that we can treat them as a single value
+// is convenient.
+
+func collect3[T any](a, b, c T) [3]T {
+	return [3]T{a, b, c}
+}
+
+func collect4[T any](a, b, c, d T) [4]T {
+	return [4]T{a, b, c, d}
+}
+
+// it's also convenient to permute these things, for
+// tests where the order shouldn't matter.
+
+var permuteOrder3 = [][3]int{
+	{0, 1, 2},
+	{0, 2, 1},
+	{1, 0, 2},
+	{1, 2, 0},
+	{2, 0, 1},
+	{2, 1, 0},
+}
+
+func permute3[T any](in [3]T, order [3]int) (out [3]T) {
+	for i, j := range order {
+		out[i] = in[j]
+	}
+
+	return
+}

internal/helper/distance.go

@@ -0,0 +1,46 @@
+package helper
+
+import (
+	"image/color"
+	"math"
+	"slices"
+)
+
+func TestDistance[T tester[T], C color.Color](t T, alpha bool, midpoint func(c0, c1 C) C, f func(c0, c1 C) float64, m color.Model) {
+	colors := slices.Collect(EnumColor[C](alpha, false, m))
+
+	for i, c0 := range colors {
+		// a colour should have a distance of zero to itself.
+		if d := f(c0, c0); !EqFloat64Fuzzy(d, 0) {
+			t.Errorf("Distance(%#+v, %#+v) = %f, want 0", c0, c0, d)
+			return
+		}
+
+		for j := i + 1; j < len(colors); j++ {
+			c1 := colors[j]
+			d, d2 := f(c0, c1), f(c1, c0)
+
+			switch {
+			case math.IsNaN(d) || math.IsInf(d, 0):
+				t.Errorf("Distance(%#+v, %#+v) = %f, want finite", c0, c1, d)
+				return
+
+			case d < 0 || EqFloat64Fuzzy(d, 0):
+				t.Errorf("Distance(%#+v, %#+v) = %f, want > 0", c0, c1, d)
+				return
+
+			case !EqFloat64Fuzzy(d, d2):
+				t.Errorf("Distance(%#+v, %#+v) != Distance(%#+v, %#+v), want %f == %f", c1, c0, c0, c1, d, d2)
+				return
+			}
+
+			// traveling from c0 to c1 via mid can't possibly be
+			// shorter than traveling from c0 to c1 directly.
+			mid := midpoint(c0, c1)
+			if cumulative := f(c0, mid) + f(mid, c1); !(d < cumulative || EqFloat64Fuzzy(d, cumulative)) {
+				t.Errorf("Distance(%#+v, %#+v)+Distance(%#+v, %#+v) < Distance(%#+v, %#+v), want %f >= %f", c0, mid, mid, c1, c0, c1, cumulative, d)
+				return
+			}
+		}
+	}
+}

internal/helper/distance_test.go

@@ -0,0 +1,116 @@
+package helper
+
+import (
+	"cmp"
+	"image/color"
+	"math"
+	"testing"
+)
+
+func TestTestDistance(t *testing.T) {
+	mt := mockTester{t: t}
+
+	midpoint := func(c0, c1 color.RGBA) color.RGBA {
+		return color.RGBA{
+			uint8((uint16(c0.R) + uint16(c1.R)) / 2),
+			uint8((uint16(c0.G) + uint16(c1.G)) / 2),
+			uint8((uint16(c0.B) + uint16(c1.B)) / 2),
+			uint8((uint16(c0.A) + uint16(c1.A)) / 2),
+		}
+	}
+
+	mt.run("non-zero distance for identical colours", func(t *mockTest) {
+		TestDistance(t, true, midpoint, func(c0, c1 color.RGBA) float64 {
+			return 1
+		}, color.RGBAModel)
+	})
+
+	mt.run("NaN distance", func(t *mockTest) {
+		TestDistance(t, true, midpoint, func(c0, c1 color.RGBA) float64 {
+			if c0 == c1 {
+				return 0
+			}
+
+			return math.NaN()
+		}, color.RGBAModel)
+	})
+
+	mt.run("negative distance", func(t *mockTest) {
+		TestDistance(t, true, midpoint, func(c0, c1 color.RGBA) float64 {
+			if c0 == c1 {
+				return 0
+			}
+
+			return -1
+		}, color.RGBAModel)
+	})
+
+	mt.run("asymmetric distance", func(t *mockTest) {
+		TestDistance(t, true, midpoint, func(c0, c1 color.RGBA) float64 {
+			if c0 == c1 {
+				return 0
+			}
+
+			if cmp.Or(int(c0.R)-int(c1.R), int(c0.G)-int(c1.G), int(c0.B)-int(c1.B), int(c0.A)-int(c1.A)) > 0 {
+				return 1
+			}
+
+			return 2
+		}, color.RGBAModel)
+	})
+
+	mt.run("triangle inequality", func(t *mockTest) {
+		TestDistance(t, true, midpoint, func(c0, c1 color.RGBA) float64 {
+			dR := int(c0.R) - int(c1.R)
+			dG := int(c0.G) - int(c1.G)
+			dB := int(c0.B) - int(c1.B)
+			dA := int(c0.A) - int(c1.A)
+
+			d2 := float64(dR*dR + dG*dG + dB*dB + dA*dA)
+
+			return d2
+		}, color.RGBAModel)
+	})
+
+	mt.run("euclidean distance", func(t *mockTest) {
+		TestDistance(t, true, midpoint, func(c0, c1 color.RGBA) float64 {
+			dR := int(c0.R) - int(c1.R)
+			dG := int(c0.G) - int(c1.G)
+			dB := int(c0.B) - int(c1.B)
+			dA := int(c0.A) - int(c1.A)
+
+			d2 := float64(dR*dR + dG*dG + dB*dB + dA*dA)
+
+			return math.Sqrt(d2)
+		}, color.RGBAModel)
+	})
+
+	mt.expectError(
+		"non-zero distance for identical colours",
+		`Distance(color.RGBA{R:0x0, G:0x0, B:0x0, A:0x0}, color.RGBA{R:0x0, G:0x0, B:0x0, A:0x0}) = 1.000000, want 0`,
+	)
+
+	mt.expectError(
+		"NaN distance",
+		`Distance(color.RGBA{R:0x0, G:0x0, B:0x0, A:0x0}, color.RGBA{R:0x0, G:0x0, B:0x0, A:0x55}) = NaN, want finite`,
+	)
+
+	mt.expectError(
+		"negative distance",
+		`Distance(color.RGBA{R:0x0, G:0x0, B:0x0, A:0x0}, color.RGBA{R:0x0, G:0x0, B:0x0, A:0x55}) = -1.000000, want > 0`,
+	)
+
+	mt.expectError(
+		"asymmetric distance",
+		`Distance(color.RGBA{R:0x0, G:0x0, B:0x0, A:0x55}, color.RGBA{R:0x0, G:0x0, B:0x0, A:0x0}) != Distance(color.RGBA{R:0x0, G:0x0, B:0x0, A:0x0}, color.RGBA{R:0x0, G:0x0, B:0x0, A:0x55}), want 2.000000 == 1.000000`,
+	)
+
+	mt.expectError(
+		"triangle inequality",
+		`Distance(color.RGBA{R:0x0, G:0x0, B:0x0, A:0x0}, color.RGBA{R:0x0, G:0x0, B:0x0, A:0x2a})+Distance(color.RGBA{R:0x0, G:0x0, B:0x0, A:0x2a}, color.RGBA{R:0x0, G:0x0, B:0x0, A:0x55}) < Distance(color.RGBA{R:0x0, G:0x0, B:0x0, A:0x0}, color.RGBA{R:0x0, G:0x0, B:0x0, A:0x55}), want 3613.000000 >= 7225.000000`,
+	)
+
+	mt.expectSuccess("euclidean distance")
+
+	mt.expectAllHandled()
+}

internal/helper/enum.go

@@ -0,0 +1,60 @@
+package helper
+
+import (
+	"image/color"
+	"iter"
+)
+
+// Enum iterates over a sparse sample of the RGBA colour space.
+//
+// If alpha is true, the colours will include transparency,
+// otherwise the returned colours will be fully opaque.
+//
+// If slow is false, an even smaller number of samples will be returned
+// making this suitable for use in a nested loop.
+//
+// alpha=true slow=true: 87481 samples.
+// alpha=false slow=true: 140608 samples.
+// alpha=true slow=false: 649 samples.
+// alpha=false slow=false: 216 samples.
+func Enum(alpha, slow bool) iter.Seq[color.RGBA64] {
+	var aStart, aStep, cDiv uint32
+
+	switch {
+	case alpha && slow:
+		aStart, aStep, cDiv = 0, 0xffff/15, 17
+	case alpha: // alpha && !slow
+		aStart, aStep, cDiv = 0, 0xffff/3, 5
+	case slow: // !alpha && slow
+		aStart, aStep, cDiv = 0xffff, 1, 51
+	default: // !alpha && !slow
+		aStart, aStep, cDiv = 0xffff, 1, 5
+	}
+
+	return func(yield func(color.RGBA64) bool) {
+		for a := aStart; a <= 0xffff; a += aStep {
+			cStep := max(1, a/cDiv)
+
+			for b := uint32(0); b <= a; b += cStep {
+				for g := uint32(0); g <= a; g += cStep {
+					for r := uint32(0); r <= a; r += cStep {
+						if !yield(color.RGBA64{uint16(r), uint16(g), uint16(b), uint16(a)}) {
+							return
+						}
+					}
+				}
+			}
+		}
+	}
+}
+
+// EnumColor is identical to [Enum], but invokes a [color.Model] to return a concrete colour type.
+func EnumColor[C color.Color](alpha, slow bool, m color.Model) iter.Seq[C] {
+	return func(yield func(C) bool) {
+		for rgba := range Enum(alpha, slow) {
+			if !yield(m.Convert(rgba).(C)) {
+				return
+			}
+		}
+	}
+}

internal/helper/enum_test.go

@@ -0,0 +1,131 @@
+package helper
+
+import (
+	"cmp"
+	"fmt"
+	"image/color"
+	"iter"
+	"slices"
+	"testing"
+)
+
+func cmpRGBA64(a, b color.RGBA64) int {
+	return cmp.Or(
+		cmp.Compare(a.R, b.R),
+		cmp.Compare(a.G, b.G),
+		cmp.Compare(a.B, b.B),
+		cmp.Compare(a.A, b.A),
+	)
+}
+
+func eqRGBA64(a, b color.RGBA64) bool {
+	return cmpRGBA64(a, b) == 0
+}
+
+func TestEnum(t *testing.T) {
+	tests := []struct {
+		alpha, slow   bool
+		expectedCount int
+	}{
+		{true, true, 87481},
+		{false, true, 140608},
+		{true, false, 649},
+		{false, false, 216},
+	}
+	for _, tt := range tests {
+		t.Run(fmt.Sprintf("alpha=%v, slow=%v", tt.alpha, tt.slow), func(t *testing.T) {
+			t.Run("sequence meets expected criteria", func(t *testing.T) {
+				list := slices.Collect(Enum(tt.alpha, tt.slow))
+				gotCount := len(list)
+				if gotCount != tt.expectedCount {
+					t.Errorf("Enum(%v, %v) returned %d items, wanted %d", tt.alpha, tt.slow, gotCount, tt.expectedCount)
+				}
+
+				slices.SortFunc(list, cmpRGBA64)
+				list = slices.CompactFunc(list, eqRGBA64)
+
+				if len(list) != gotCount {
+					t.Errorf("Enum(%v, %v) returned %d duplicate items", tt.alpha, tt.slow, gotCount-len(list))
+				}
+
+				listHasAlpha := false
+				for _, c := range list {
+					if c.A != 0xffff {
+						if !tt.alpha {
+							t.Errorf("Enum(%v, %v) returned non-opaque color: %v", tt.alpha, tt.slow, c)
+						}
+
+						listHasAlpha = true
+						break
+					}
+				}
+
+				if !listHasAlpha && tt.alpha {
+					t.Errorf("Enum(%v, %v) didn't return non-opaque colors", tt.alpha, tt.slow)
+				}
+			})
+
+			t.Run("cancel", func(t *testing.T) {
+				// make sure cancelling the iteration doesn't panic.
+				// But mostly, we want that sweet, sweet test coverage.
+				next, stop := iter.Pull(Enum(tt.alpha, tt.slow))
+				// need to invoke next to actually have the generated be started.
+				if _, ok := next(); !ok {
+					t.Error("iteration stopped before we could cancel it")
+				}
+				stop()
+			})
+		})
+	}
+
+}
+
+func TestEnumColor(t *testing.T) {
+	tests := []struct {
+		alpha, slow bool
+	}{
+		{true, true},
+		{false, true},
+		{true, false},
+		{false, false},
+	}
+	for _, tt := range tests {
+		t.Run(fmt.Sprintf("alpha=%v, slow=%v", tt.alpha, tt.slow), func(t *testing.T) {
+			t.Run("sequence equivalence", func(t *testing.T) {
+				nextRGBA64, stop1 := iter.Pull(Enum(tt.alpha, tt.slow))
+				defer stop1()
+				nextNRGBA, stop2 := iter.Pull(EnumColor[color.NRGBA](tt.alpha, tt.slow, color.NRGBAModel))
+				defer stop2()
+
+				for i := 0; ; i++ {
+					rgba64, gotRgba64 := nextRGBA64()
+					nrgba, gotNrgba := nextNRGBA()
+
+					if gotRgba64 != gotNrgba {
+						t.Errorf("one sequence terminated at i=%d: gotRgba64=%v, gotNrgba=%v", i, gotRgba64, gotNrgba)
+						return
+					}
+
+					if !gotRgba64 {
+						break
+					}
+
+					if wantNrgba := color.NRGBAModel.Convert(rgba64).(color.NRGBA); nrgba != wantNrgba {
+						t.Errorf("i=%d: got %#+v, expected %#+v", i, nrgba, wantNrgba)
+					}
+				}
+			})
+
+			t.Run("cancel", func(t *testing.T) {
+				// make sure cancelling the iteration doesn't panic.
+				// But mostly, we want that sweet, sweet test coverage.
+				next, stop := iter.Pull(EnumColor[color.NRGBA](tt.alpha, tt.slow, color.NRGBAModel))
+				// need to invoke next to actually have the generated be started.
+				if _, ok := next(); !ok {
+					t.Error("iteration stopped before we could cancel it")
+				}
+				stop()
+			})
+		})
+	}
+}

internal/helper/gamma.go

@@ -0,0 +1,28 @@
+package helper
+
+import "math"
+
+// Linearize converts an sRGB component in the range [0, 0xffff] to a linearRGB component in the range [0, 1].
+func Linearize(c uint32) float64 {
+	l := float64(c) / 0xffff
+
+	if l <= 0.039285714285714285714285714285714 {
+		return l / 12.923210180787861094641554898407
+	}
+
+	return math.Pow((l+0.055)/1.055, 2.4)
+}
+
+// Delinearize converts a linearRGB component in the range [0, 1] to an sRGB component in the range [0, 0xffff].
+func Delinearize(l float64) uint32 {
+	switch {
+	case l <= 0:
+		return 0
+	case l <= 0.0030399346397784299969770436366690:
+		return uint32(l*846922.57919793247683733430026710 + 0.5)
+	case l >= 1:
+		return 0xffff
+	default:
+		return uint32(69139.425*math.Pow(l, 1/2.4) - 3603.925)
+	}
+}

internal/helper/gamma_test.go

@@ -0,0 +1,84 @@
+package helper
+
+import (
+	"fmt"
+	"math"
+	"testing"
+)
+
+func TestLinearize(t *testing.T) {
+	tests := []struct {
+		value uint32
+		want  float64
+	}{
+		// the minimum and maximum legal values should map to 0 and 1, respectively.
+		{0x0000, 0.0},
+		{0xffff, 1.0},
+
+		// check what would be 0x01 in 8-bit sRGB.
+		// this is below the point where it the function switches from linear to exponential.
+		{0x0101, 0x1.3e312a36f1977p-12},
+
+		// check the midpoint of the gamma curve.
+		{0x7fff, 0x1.b6577fc57aa37p-03},
+
+		// We do support values beyond the maximum legal value, although you probably shouldn't depend on it
+		// as the converse function will map these to the maximum legal value, creating an asymmetry.
+		{0x10000, 0x1.000246626b604p+00},
+		{math.MaxUint32, 0x1.2912a0c535107p+38},
+	}
+	for _, tt := range tests {
+		t.Run(fmt.Sprintf("0x%04x", tt.value), func(t *testing.T) {
+			if got := Linearize(tt.value); !EqFloat64Fuzzy(got, tt.want) {
+				t.Errorf("Linearize(0x%04x) = %x: want %x", tt.value, got, tt.want)
+			}
+		})
+	}
+
+	t.Run("monotonically increasing", func(t *testing.T) {
+		for i, prev := uint32(1), Linearize(0); i < 0x10000; i++ {
+			got := Linearize(i)
+			if got <= prev {
+				t.Errorf("Linearize(0x%04x) = %x; want > %x", i, got, prev)
+			}
+			prev = got
+		}
+	})
+}
+
+func TestDelinearize(t *testing.T) {
+	tests := []struct {
+		value float64
+		want  uint32
+	}{
+		// make sure clamping to legal values is happening.
+		{-1, 0x0000},
+		{2, 0xffff},
+
+		// again with the next values below 0 and above 1.
+		{math.Nextafter(0, math.Inf(-1)), 0},
+		{math.Nextafter(1, math.Inf(1)), 0xffff},
+
+		// and lastly, the non-finites.
+		{math.Inf(-1), 0x0000},
+		{math.Inf(1), 0xffff},
+		{math.NaN(), 0x0000},
+	}
+	for _, tt := range tests {
+		t.Run(fmt.Sprintf("%x", tt.value), func(t *testing.T) {
+			if got := Delinearize(tt.value); got != tt.want {
+				t.Errorf("Delinearize(%x) = 0x%04x: want 0x%04x", tt.value, got, tt.want)
+			}
+		})
+	}
+
+	t.Run("lossless conversion of legal values", func(t *testing.T) {
+		for c := uint32(0); c < 0x10000; c++ {
+			got := Delinearize(Linearize(c))
+			if got != c {
+				t.Errorf("Delinearize(Linearize(0x%04x)) != 0x%04x", c, got)
+				return
+			}
+		}
+	})
+}

internal/helper/model.go

@@ -0,0 +1,92 @@
+package helper
+
+import (
+	"image/color"
+)
+
+func Model[C color.Color](fromColor func(color.Color) C) color.Model {
+	return color.ModelFunc(func(c color.Color) color.Color {
+		return fromColor(c)
+	})
+}
+
+// Interface that the colours used in this package are expected to implement.
+type Color interface {
+	comparable
+	color.Color
+	NRGBA() (r, g, b, a uint32)
+	NLRGBA() (r, g, b, a float64)
+	NXYZA() (x, y, z, a float64)
+	NOkLabA() (lightness, chromaA, chromaB, a float64)
+}
+
+type ConvertTest[C Color] struct {
+	Name string
+	In   color.Color
+	Out  C
+}
+
+func TestModel[T tester[T], C Color](t T, alpha bool, m color.Model, eq func(c0, c1 C) bool, extra []ConvertTest[C]) {
+	t.Run("legal colours", func(t T) {
+		for wantRGBA := range Enum(alpha, true) {
+			_gotC := m.Convert(wantRGBA)
+			gotC, ok := _gotC.(C)
+
+			if !ok {
+				t.Errorf("model.Convert(%#+v) returned %T, expected %T", wantRGBA, _gotC, gotC)
+				return
+			}
+
+			r, g, b, a := gotC.RGBA()
+			gotRGBA := color.RGBA64{uint16(r), uint16(g), uint16(b), uint16(a)}
+
+			if gotRGBA != wantRGBA {
+				t.Errorf("%#+v.RGBA() = %v, want %v", gotC, gotRGBA, wantRGBA)
+				return
+			}
+
+			wantNRGBA := color.NRGBA64Model.Convert(wantRGBA)
+			r, g, b, a = gotC.NRGBA()
+			gotNRGBA := color.NRGBA64{uint16(r), uint16(g), uint16(b), uint16(a)}
+			if gotNRGBA != wantNRGBA {
+				t.Errorf("%#+v.NRGBA() = %v, want %v", gotC, gotNRGBA, wantNRGBA)
+				return
+			}
+
+			wantNLRGBA := collect4(NRGBAtoNLRGBA(r, g, b, a))
+			if gotNLRGBA := collect4(gotC.NLRGBA()); !EqFloat64SliceFuzzy(gotNLRGBA[:], wantNLRGBA[:]) {
+				t.Errorf("%#+v.NLRGBA() = %v, want %v", gotC, gotNLRGBA, wantNLRGBA)
+				return
+			}
+
+			var wantNXYZA [4]float64
+			wantNXYZA[0], wantNXYZA[1], wantNXYZA[2] = LRGBtoXYZ(wantNLRGBA[0], wantNLRGBA[1], wantNLRGBA[2])
+			wantNXYZA[3] = wantNLRGBA[3]
+
+			if gotNXYZA := collect4(gotC.NXYZA()); !EqFloat64SliceFuzzy(gotNXYZA[:], wantNXYZA[:]) {
+				t.Errorf("%#+v.NXYZA() = %v want %v", gotC, gotNXYZA, wantNXYZA)
+				return
+			}
+
+			var wantNOkLabA [4]float64
+			wantNOkLabA[0], wantNOkLabA[1], wantNOkLabA[2] = LMStoOkLab(LRGBtoLMS(wantNLRGBA[0], wantNLRGBA[1], wantNLRGBA[2]))
+			wantNOkLabA[3] = wantNLRGBA[3]
+
+			if gotNOkLabA := collect4(gotC.NOkLabA()); !EqFloat64SliceFuzzy(gotNOkLabA[:], wantNOkLabA[:]) {
+				t.Errorf("%#+v.NOkLabA()[:3] = %v want %v", gotC, gotNOkLabA[:], wantNOkLabA[:])
+				return
+			}
+		}
+	})
+
+	for _, tt := range extra {
+		t.Run(tt.Name, func(t T) {
+			gotC := m.Convert(tt.In).(C)
+
+			if !eq(gotC, tt.Out) {
+				t.Errorf("model.Convert(%#+v) = %#+v, want %#+v", tt.In, gotC, tt.Out)
+				return
+			}
+		})
+	}
+}

internal/helper/model_test.go

@@ -0,0 +1,168 @@
+package helper
+
+import (
+	"image/color"
+	"testing"
+)
+
+func eq[T comparable](a, b T) bool {
+	return a == b
+}
+
+type nrgba64 struct {
+	color.NRGBA64
+}
+
+func (c *nrgba64) set(v color.NRGBA64) {
+	c.NRGBA64 = v
+}
+
+func (c nrgba64) NRGBA() (_, _, _, _ uint32) {
+	return uint32(c.R), uint32(c.G), uint32(c.B), uint32(c.A)
+}
+
+func (c nrgba64) NLRGBA() (_, _, _, _ float64) {
+	return NRGBAtoNLRGBA(c.NRGBA())
+}
+
+func (c nrgba64) NXYZA() (_, _, _, _ float64) {
+	r, g, b, a := c.NLRGBA()
+	x, y, z := LRGBtoXYZ(r, g, b)
+	return x, y, z, a
+}
+
+func (c nrgba64) NOkLabA() (_, _, _, _ float64) {
+	r, g, b, a := c.NLRGBA()
+	lightness, chromaA, chromaB := LMStoOkLab(LRGBtoLMS(r, g, b))
+	return lightness, chromaA, chromaB, a
+}
+
+func convert[C Color, P interface {
+	*C
+	set(color.NRGBA64)
+}](c color.Color) C {
+	var result C
+	P(&result).set(color.NRGBA64Model.Convert(c).(color.NRGBA64))
+	return result
+}
+
+type nrgba64BadRGBA struct {
+	nrgba64
+}
+
+func (nrgba64BadRGBA) RGBA() (_, _, _, _ uint32) {
+	return 1, 2, 3, 4
+}
+
+type nrgba64BadNRGBA struct {
+	nrgba64
+}
+
+func (nrgba64BadNRGBA) NRGBA() (_, _, _, _ uint32) {
+	return 1, 2, 3, 4
+}
+
+type nrgba64BadNLRGBA struct {
+	nrgba64
+}
+
+func (nrgba64BadNLRGBA) NLRGBA() (_, _, _, _ float64) {
+	return 1, 2, 3, 4
+}
+
+type nrgba64BadNXYZA struct {
+	nrgba64
+}
+
+func (nrgba64BadNXYZA) NXYZA() (_, _, _, _ float64) {
+	return 1, 2, 3, 4
+}
+
+type nrgba64BadNOkLabA struct {
+	nrgba64
+}
+
+func (nrgba64BadNOkLabA) NOkLabA() (_, _, _, _ float64) {
+	return 1, 2, 3, 4
+}
+
+func TestTestModel(t *testing.T) {
+	mt := mockTester{t: t}
+
+	mt.run("wrong colour type", func(t *mockTest) {
+		TestModel(t, true, color.RGBAModel, eq[nrgba64], nil)
+	})
+
+	mt.run("bad RGBA", func(t *mockTest) {
+		TestModel(t, false, Model(convert[nrgba64BadRGBA]), eq[nrgba64BadRGBA], nil)
+	})
+
+	mt.run("bad NRGBA", func(t *mockTest) {
+		TestModel(t, false, Model(convert[nrgba64BadNRGBA]), eq[nrgba64BadNRGBA], nil)
+	})
+
+	mt.run("bad NLRGBA", func(t *mockTest) {
+		TestModel(t, false, Model(convert[nrgba64BadNLRGBA]), eq[nrgba64BadNLRGBA], nil)
+	})
+
+	mt.run("bad NXYZA", func(t *mockTest) {
+		TestModel(t, false, Model(convert[nrgba64BadNXYZA]), eq[nrgba64BadNXYZA], nil)
+	})
+
+	mt.run("bad NOkLabA", func(t *mockTest) {
+		TestModel(t, false, Model(convert[nrgba64BadNOkLabA]), eq[nrgba64BadNOkLabA], nil)
+	})
+
+	mt.run("working model", func(t *mockTest) {
+		TestModel(t, true, Model(convert[nrgba64]), eq[nrgba64], []ConvertTest[nrgba64]{
+			{"good", color.NRGBA64{0x0123, 0x4567, 0x89ab, 0xcdef}, nrgba64{color.NRGBA64{0x0123, 0x4567, 0x89ab, 0xcdef}}},
+			{"bad", color.NRGBA64{0xcafe, 0xf00d, 0x54ac, 0xce55}, nrgba64{color.NRGBA64{0x0123, 0x4567, 0x89ab, 0xcdef}}},
+		})
+	})
+
+	mt.expectFailedChildren("wrong colour type")
+	mt.expectError(
+		"wrong colour type/legal colours",
+		`model.Convert(color.RGBA64{R:0x0, G:0x0, B:0x0, A:0x0}) returned color.RGBA, expected helper.nrgba64`,
+	)
+
+	mt.expectFailedChildren("bad RGBA")
+	mt.expectError(
+		"bad RGBA/legal colours",
+		`helper.nrgba64BadRGBA{nrgba64:helper.nrgba64{NRGBA64:color.NRGBA64{R:0x0, G:0x0, B:0x0, A:0xffff}}}.RGBA() = {1 2 3 4}, want {0 0 0 65535}`,
+	)
+
+	mt.expectFailedChildren("bad NRGBA")
+	mt.expectError(
+		"bad NRGBA/legal colours",
+		`helper.nrgba64BadNRGBA{nrgba64:helper.nrgba64{NRGBA64:color.NRGBA64{R:0x0, G:0x0, B:0x0, A:0xffff}}}.NRGBA() = {1 2 3 4}, want {0 0 0 65535}`,
+	)
+
+	mt.expectFailedChildren("bad NLRGBA")
+	mt.expectError(
+		"bad NLRGBA/legal colours",
+		`helper.nrgba64BadNLRGBA{nrgba64:helper.nrgba64{NRGBA64:color.NRGBA64{R:0x0, G:0x0, B:0x0, A:0xffff}}}.NLRGBA() = [1 2 3 4], want [0 0 0 1]`,
+	)
+
+	mt.expectFailedChildren("bad NXYZA")
+	mt.expectError(
+		"bad NXYZA/legal colours",
+		`helper.nrgba64BadNXYZA{nrgba64:helper.nrgba64{NRGBA64:color.NRGBA64{R:0x0, G:0x0, B:0x0, A:0xffff}}}.NXYZA() = [1 2 3 4] want [0 0 0 1]`,
+	)
+
+	mt.expectFailedChildren("bad NOkLabA")
+	mt.expectError(
+		"bad NOkLabA/legal colours",
+		`helper.nrgba64BadNOkLabA{nrgba64:helper.nrgba64{NRGBA64:color.NRGBA64{R:0x0, G:0x0, B:0x0, A:0xffff}}}.NOkLabA()[:3] = [1 2 3 4] want [0 0 0 1]`,
+	)
+
+	mt.expectFailedChildren("working model")
+	mt.expectSuccess("working model/legal colours")
+	mt.expectSuccess("working model/good")
+	mt.expectError(
+		"working model/bad",
+		`model.Convert(color.NRGBA64{R:0xcafe, G:0xf00d, B:0x54ac, A:0xce55}) = helper.nrgba64{NRGBA64:color.NRGBA64{R:0xcafe, G:0xf00d, B:0x54ac, A:0xce55}}, want helper.nrgba64{NRGBA64:color.NRGBA64{R:0x123, G:0x4567, B:0x89ab, A:0xcdef}}`,
+	)
+
+	mt.expectAllHandled()
+}

internal/helper/oklab.go

@@ -0,0 +1,103 @@
+package helper
+
+import (
+	"image/color"
+	"math"
+)
+
+func XYZtoLMS(x, y, z float64) (_, _, _ float64) {
+	// https://bottosson.github.io/posts/oklab/
+	//
+	// OkLab first defines a transform from xyz to an intermediate space:
+	//
+	// [+0.8189330101, +0.3618667424, -0.1288597137]
+	// [+0.0329845436, +0.9293118715, +0.0361456387]
+	// [+0.0482003018, +0.2643662691, +0.633851707 ]
+	//
+	// Inverse:
+	//
+	// [+1.2270138511035210261251539010893,   -0.55779998065182223833890733780747,  +0.28125614896646780760667886762980 ]
+	// [-0.040580178423280593980748617561551, +1.1122568696168301049956590765194,   -0.071676678665601200581102747142872]
+	// [-0.076381284505706892872271894590358, -0.42148197841801273056818761141308,  +1.5861632204407947575338479416771  ]
+
+	return +0.8189330101*x + 0.3618667424*y - 0.1288597137*z,
+		0.0329845436*x + 0.9293118715*y + 0.0361456387*z,
+		0.0482003018*x + 0.2643662691*y + 0.633851707*z
+}
+
+func LMStoXYZ(l, m, s float64) (_, _, _ float64) {
+	return 1.2270138511035210261251539010893*l - 0.55779998065182223833890733780747*m + 0.28125614896646780760667886762980*s,
+		-0.040580178423280593980748617561551*l + 1.1122568696168301049956590765194*m - 0.071676678665601200581102747142872*s,
+		-0.076381284505706892872271894590358*l - 0.42148197841801273056818761141308*m + 1.5861632204407947575338479416771*s
+}
+
+func LRGBtoLMS(r, g, b float64) (_, _, _ float64) {
+	// we can combine the LRGB to D65 CIE XYZ transform and the XYZ to LMS transform above to go straight from LRGB to LMS:
+	//
+	// [+0.41217385032507, +0.5362974607032, +0.05146302925248]
+	// [+0.21187214048845, +0.6807476834212, +0.10740645682645]
+	// [+0.08831541121808, +0.2818663070584, +0.63026344660742]
+	//
+	// Inverse:
+	//
+	// [+4.0767584135565013494237930518854,    -3.3072279873944731418619352916485,  +0.23072145994488563247301883404900]
+	// [-1.2681810851624033989047813181437,    +2.6092932102856398573991970933594,  -0.34111211654775355696796160418220]
+	// [-0.0040984077180314400491332639337372, -0.70350366010241732765095902557887, +1.7068604529788013559365593912662]
+
+	return +0.41217385032507*r + 0.5362974607032*g + 0.05146302925248*b,
+		0.21187214048845*r + 0.6807476834212*g + 0.10740645682645*b,
+		0.08831541121808*r + 0.2818663070584*g + 0.63026344660742*b
+}
+
+func LMStoLRGB(l, m, s float64) (_, _, _ float64) {
+	return +4.0767584135565013494237930518854*l - 3.3072279873944731418619352916485*m + 0.23072145994488563247301883404900*s,
+		-1.2681810851624033989047813181437*l + 2.6092932102856398573991970933594*m - 0.34111211654775355696796160418220*s,
+		-0.0040984077180314400491332639337372*l - 0.70350366010241732765095902557887*m + 1.7068604529788013559365593912662*s
+}
+
+func LMStoOkLab(l, m, s float64) (_, _, _ float64) {
+	// After a non-linear transformation (cube root), OkLab applies a second matrix:
+	//
+	// [+0.2104542553, +0.793617785,  -0.0040720468]
+	// [+1.9779984951, -2.428592205,  +0.4505937099]
+	// [+0.0259040371, +0.7827717662, -0.808675766 ]
+	//
+	// Inverse:
+	//
+	// [+0.99999999845051981426207542502031, +0.39633779217376785682345989261573,  +0.21580375806075880342314146183004 ]
+	// [+1.0000000088817607767160752456705,  -0.10556134232365634941095687705472,  -0.063854174771705903405254198817796]
+	// [+1.0000000546724109177012928651534,  -0.089484182094965759689052745863391, -1.2914855378640917399489287529148  ]
+
+	lP, mP, sP := math.Cbrt(l), math.Cbrt(m), math.Cbrt(s)
+
+	return 0.2104542553*lP + 0.793617785*mP - 0.0040720468*sP,
+		1.9779984951*lP - 2.428592205*mP + 0.4505937099*sP,
+		0.0259040371*lP + 0.7827717662*mP - 0.808675766*sP
+}
+
+func cube(v float64) float64 {
+	return v * v * v
+}
+
+func OkLabToLMS(l, a, b float64) (_, _, _ float64) {
+	return cube(0.99999999845051981426207542502031*l + 0.39633779217376785682345989261573*a + 0.21580375806075880342314146183004*b),
+		cube(1.0000000088817607767160752456705*l - 0.10556134232365634941095687705472*a - 0.063854174771705903405254198817796*b),
+		cube(1.0000000546724109177012928651534*l - 0.089484182094965759689052745863391*a - 1.2914855378640917399489287529148*b)
+}
+
+func ColorToNOkLabA(c color.Color) (lightness, chromaA, chromaB, a float64) {
+	switch c := c.(type) {
+	case interface {
+		NOkLabA() (lightness, chromaA, chromaB, a float64)
+	}:
+		return c.NOkLabA()
+	case interface{ NXYZA() (x, y, z, a float64) }:
+		x, y, z, a := c.NXYZA()
+		lightness, chromaA, chromaB = LMStoOkLab(XYZtoLMS(x, y, z))
+		return lightness, chromaA, chromaB, a
+	default:
+		r, g, b, a := ColorToNLRGBA(c)
+		lightness, chromaA, chromaB = LMStoOkLab(LRGBtoLMS(r, g, b))
+		return lightness, chromaA, chromaB, a
+	}
+}

internal/helper/oklab_test.go

@@ -0,0 +1,136 @@
+package helper
+
+import (
+	"fmt"
+	"image/color"
+	"math"
+	"testing"
+)
+
+func TestLMSToXYZ(t *testing.T) {
+	for c := range Enum(false, true) {
+		want := collect3(LRGBtoXYZ(RGBtoLRGB(uint32(c.R), uint32(c.G), uint32(c.B))))
+
+		if got := collect3(LMStoXYZ(XYZtoLMS(want[0], want[1], want[2]))); !EqFloat64SliceFuzzy(want[:], got[:]) {
+			t.Errorf("LMStoXYZ(XYZtoLMS(%v)) = %v, want unchanged", want, got)
+			return
+		}
+	}
+}
+
+func TestLMSToLRGB(t *testing.T) {
+	for c := range Enum(false, true) {
+		want := collect3(RGBtoLRGB(uint32(c.R), uint32(c.G), uint32(c.B)))
+
+		l, m, s := LRGBtoLMS(want[0], want[1], want[2])
+
+		// test via the optimized LMStoLRGB function.
+		if got := collect3(LMStoLRGB(l, m, s)); !EqFloat64SliceFuzzy(want[:], got[:]) {
+			t.Errorf("LMStoLRGB(LRGBtoLMS(%v)) = %v, want unchanged", want, got)
+			return
+		}
+
+		// make sure this is equivalent to going through the XYZ colourspace.
+		if got := collect3(XYZtoLRGB(LMStoXYZ(l, m, s))); !EqFloat64SliceFuzzy(want[:], got[:]) {
+			t.Errorf("XYZtoLRGB(LMStoXYZ(LRGBtoLMS(%v))) = %v, want unchanged", want, got)
+			return
+		}
+	}
+}
+
+func TestOKLabToLMS(t *testing.T) {
+	for c := range Enum(false, true) {
+		want := collect3(LRGBtoLMS(RGBtoLRGB(uint32(c.R), uint32(c.G), uint32(c.B))))
+		if got := collect3(OkLabToLMS(LMStoOkLab(want[0], want[1], want[2]))); !EqFloat64SliceFuzzy(want[:], got[:]) {
+			t.Errorf("OkLabToLMS(LMStoOKLab(%v)) = %v, want unchanged", want, got)
+			return
+		}
+	}
+}
+
+func TestOkLabExamplePairs(t *testing.T) {
+	// The page https://bottosson.github.io/posts/oklab/ lists example XYZ and OkLab pairs,
+	// with the results rounded to three decimal places.
+	examples := []struct{ xyz, lab [3]float64 }{
+		{[3]float64{0.950, 1.000, 1.089}, [3]float64{1.000, 0.000, 0.000}},
+		{[3]float64{1.000, 0.000, 0.000}, [3]float64{0.450, 1.236, -0.019}},
+		{[3]float64{0.000, 1.000, 0.000}, [3]float64{0.922, -0.671, 0.263}},
+		{[3]float64{0.000, 0.000, 1.000}, [3]float64{0.153, -1.415, -0.449}},
+	}
+
+	round := func(x float64) float64 {
+		return math.Round(x*1000) / 1000
+	}
+
+	round3 := func(a, b, c float64) [3]float64 {
+		return [3]float64{round(a), round(b), round(c)}
+	}
+
+	for i, e := range examples {
+		if gotLab := round3(LMStoOkLab(XYZtoLMS(e.xyz[0], e.xyz[1], e.xyz[2]))); gotLab != e.lab {
+			t.Errorf("pair %d: computed lab=%v, want=%v", i+1, gotLab, e.lab)
+		}
+
+		// note that the example table isn't suitable fo testing OkLab to XYZ conversion due to
+		// the errors introduced by rounding.
+		//
+		// we are depending the round trip conversions being correct, which is verified by TestLMSToXYZ and TestOKLabToLMS.
+	}
+}
+
+type testNXYZAColor [4]float64
+
+func (c testNXYZAColor) RGBA() (_, _, _, _ uint32) {
+	panic("should not be called")
+}
+
+func (c testNXYZAColor) NXYZA() (_, _, _, _ float64) {
+	return c[0], c[1], c[2], c[3]
+}
+
+type testNOkLabAColor [4]float64
+
+func (c testNOkLabAColor) RGBA() (_, _, _, _ uint32) {
+	panic("should not be called")
+}
+
+func (c testNOkLabAColor) NOkLabA() (_, _, _, _ float64) {
+	return c[0], c[1], c[2], c[3]
+}
+
+func TestColorToNOkLabA(t *testing.T) {
+	tests := []struct {
+		input color.Color
+		want  [4]float64
+	}{
+		{
+			// test special NRGBA handling.
+			color.NRGBA{0x01, 0x23, 0x45, 0x67},
+			[4]float64{0.25462381167525894, -0.02293028913883799, -0.07098467472369072, float64(0x6767) / 0xffff},
+		}, {
+			// test special NRGBA64 handling.
+			color.NRGBA64{0x0123, 0x4567, 0x89ab, 0},
+			[4]float64{0.39601873251000413, -0.03369278598612779, -0.12401844116020128, 0},
+		}, {
+			// test a colour that can return its linear NRGBA values directly.
+			testNLRGBA{color.NRGBA64{0x0123, 0x4567, 0x89ab, 0xcdef}},
+			[4]float64{0.39601873251000413, -0.03369278598612779, -0.12401844116020128, float64(0xcdef) / 0xffff},
+		}, {
+			// test conversion of the values from a a colour that can return NXYZA values directly.
+			testNXYZAColor{0.95, 1., 1.089, .5},
+			// these were from the canonical test pairs, these values are 1, 0, 0 when rounded to the nearest thousandth.
+			[4]float64{0.9999686754143632, -0.0002580058168537569, -0.00011499756458199784, .5},
+		}, {
+			// test that we get the values from a colour that can return NOkLabA directly.
+			testNOkLabAColor{math.Inf(1), math.NaN(), math.Inf(-1), -1},
+			[4]float64{math.Inf(1), math.NaN(), math.Inf(-1), -1},
+		},
+	}
+	for _, tt := range tests {
+		t.Run(fmt.Sprintf("%#+v", tt.input), func(t *testing.T) {
+			if got := collect4(ColorToNOkLabA(tt.input)); !EqFloat64SliceFuzzy(got[:], tt.want[:]) {
+				t.Errorf("ColorToNOkLabA(%#+v) = %v, want %v", tt.input, got, tt.want)
+			}
+		})
+	}
+}

internal/helper/rgb.go

@@ -0,0 +1,74 @@
+package helper
+
+import (
+	"image/color"
+)
+
+func RGBtoLRGB(r, g, b uint32) (_, _, _ float64) {
+	return Linearize(r), Linearize(g), Linearize(b)
+}
+
+func NRGBAtoNLRGBA(r, g, b, a uint32) (_, _, _, _ float64) {
+	_r, _g, _b := RGBtoLRGB(r, g, b)
+	return _r, _g, _b, float64(a) / 0xffff
+}
+
+func RGBAtoNLRGBA(r, g, b, a uint32) (_, _, _, _ float64) {
+	switch a {
+	case 0:
+		return 0, 0, 0, 0
+	case 0xffff:
+		_r, _g, _b := RGBtoLRGB(r, g, b)
+		return _r, _g, _b, 1
+	default:
+		// note that we round up here, as the inverse function, ToRGBA, rounds down.
+		return NRGBAtoNLRGBA((r*0xffff+a-1)/a, (g*0xffff+a-1)/a, (b*0xffff+a-1)/a, a)
+	}
+}
+
+func ColorToNLRGBA(c color.Color) (_, _, _, _ float64) {
+	switch c := c.(type) {
+	case color.NRGBA:
+		return NRGBAtoNLRGBA(uint32(c.R)*0x101, uint32(c.G)*0x101, uint32(c.B)*0x101, uint32(c.A)*0x101)
+	case color.NRGBA64:
+		return NRGBAtoNLRGBA(uint32(c.R), uint32(c.G), uint32(c.B), uint32(c.A))
+	case interface{ NLRGBA() (r, g, b, a float64) }:
+		return c.NLRGBA()
+	case interface{ NXYZA() (x, y, z, a float64) }:
+		x, y, z, a := c.NXYZA()
+		r, g, b := XYZtoLRGB(x, y, z)
+		return r, g, b, a
+	case interface {
+		NOkLabA() (lightness, chromaA, chromaB, a float64)
+	}:
+		lightness, chromaA, chromaB, a := c.NOkLabA()
+		r, g, b := LMStoLRGB(OkLabToLMS(lightness, chromaA, chromaB))
+		return r, g, b, a
+
+	case interface{ NRGBA() (r, g, b, a uint32) }:
+		return NRGBAtoNLRGBA(c.NRGBA())
+	default:
+		return RGBAtoNLRGBA(c.RGBA())
+	}
+}
+
+func LRGBtoRGB(r, g, b float64) (_, _, _ uint32) {
+	r, g, b = ClampRGB(r, g, b)
+	return Delinearize(r), Delinearize(g), Delinearize(b)
+}
+
+func NLRGBAtoNRGBA(r, g, b, a float64) (_, _, _, _ uint32) {
+	_r, _g, _b := LRGBtoRGB(r, g, b)
+	return _r, _g, _b, uint32(min(1, max(0, a))*0xffff + .5)
+}
+
+func NLRGBAtoRGBA(r, g, b, a float64) (_, _, _, _ uint32) {
+	switch _r, _g, _b, _a := NLRGBAtoNRGBA(r, g, b, a); _a {
+	case 0:
+		return 0, 0, 0, 0
+	case 0xffff:
+		return _r, _g, _b, 0xffff
+	default:
+		return _r * _a / 0xffff, _g * _a / 0xffff, _b * _a / 0xffff, _a
+	}
+}

internal/helper/rgb_test.go

@@ -0,0 +1,123 @@
+package helper
+
+import (
+	"fmt"
+	"image/color"
+	"testing"
+)
+
+type testNRGBA struct {
+	color.NRGBA64
+}
+
+func (c testNRGBA) NRGBA() (_, _, _, _ uint32) {
+	return uint32(c.R), uint32(c.G), uint32(c.B), uint32(c.A)
+}
+
+type testNLRGBA struct {
+	color.NRGBA64
+}
+
+func (c testNLRGBA) NLRGBA() (_, _, _, _ float64) {
+	return Linearize(uint32(c.R)), Linearize(uint32(c.G)), Linearize(uint32(c.B)), float64(c.A) / 0xffff
+}
+
+type testNXYZA struct {
+	color.NRGBA64
+}
+
+func (c testNXYZA) NXYZA() (_, _, _, _ float64) {
+	x, y, z := LRGBtoXYZ(Linearize(uint32(c.R)), Linearize(uint32(c.G)), Linearize(uint32(c.B)))
+	return x, y, z, float64(c.A) / 0xffff
+}
+
+type testNOkLabA struct {
+	color.NRGBA64
+}
+
+func (c testNOkLabA) NOkLabA() (_, _, _, _ float64) {
+	l, a, b := LMStoOkLab(LRGBtoLMS(Linearize(uint32(c.R)), Linearize(uint32(c.G)), Linearize(uint32(c.B))))
+	return l, a, b, float64(c.A) / 0xffff
+}
+
+func TestColorToNLRGBA(t *testing.T) {
+	tests := []struct {
+		input color.Color
+		want  [4]float64
+	}{
+		{
+			// test special NRGBA handling.
+			color.NRGBA{0x01, 0x23, 0x45, 0x67},
+			[4]float64{Linearize(0x0101), Linearize(0x2323), Linearize(0x4545), float64(0x6767) / 0xffff},
+		}, {
+			// test special NRGBA64 handling.
+			color.NRGBA64{0x0123, 0x4567, 0x89ab, 0xcdef},
+			[4]float64{Linearize(0x0123), Linearize(0x4567), Linearize(0x89ab), float64(0xcdef) / 0xffff},
+		}, {
+			// test a colour that can returns is NRGBA values directly.
+			testNRGBA{color.NRGBA64{0x0123, 0x4567, 0x89ab, 0xcdef}},
+			[4]float64{Linearize(0x0123), Linearize(0x4567), Linearize(0x89ab), float64(0xcdef) / 0xffff},
+		}, {
+			// test a colour that can return its NLRGBA values directly.
+			testNLRGBA{color.NRGBA64{0x0123, 0x4567, 0x89ab, 0xcdef}},
+			[4]float64{Linearize(0x0123), Linearize(0x4567), Linearize(0x89ab), float64(0xcdef) / 0xffff},
+		}, {
+			// test a colour that can returns its NXYZA values directly.
+			testNXYZA{color.NRGBA64{0x0123, 0x4567, 0x89ab, 0xcdef}},
+			[4]float64{Linearize(0x0123), Linearize(0x4567), Linearize(0x89ab), float64(0xcdef) / 0xffff},
+		}, {
+			// test a colour that can returns its NOkLabA values directly.
+			testNOkLabA{color.NRGBA64{0x0123, 0x4567, 0x89ab, 0xcdef}},
+			[4]float64{Linearize(0x0123), Linearize(0x4567), Linearize(0x89ab), float64(0xcdef) / 0xffff},
+		}, {
+			// the FromRGBA codepath with partial transparency
+			color.RGBA64{0x0123, 0x4567, 0x89ab, 0xcdef},
+			[4]float64{
+				Linearize((0x0123*0xffff + 0xcdee - 1) / 0xcdef),
+				Linearize((0x4567*0xffff + 0xcdee - 1) / 0xcdef),
+				Linearize((0x89ab*0xffff + 0xcdee - 1) / 0xcdef),
+				float64(0xcdef) / 0xffff,
+			},
+		}, {
+			// the FromRGBA codepath with full transparency
+			color.RGBA64{0x0000, 0x0000, 0x0000, 0x0000},
+			[4]float64{0, 0, 0, 0},
+		}, {
+			// the FromRGBA codepath with full opacity
+			color.RGBA64{0x0123, 0x4567, 0x89ab, 0xffff},
+			[4]float64{Linearize(0x0123), Linearize(0x4567), Linearize(0x89ab), 1},
+		},
+	}
+	for _, tt := range tests {
+		t.Run(fmt.Sprintf("%#+v", tt.input), func(t *testing.T) {
+			if got := collect4(ColorToNLRGBA(tt.input)); !EqFloat64SliceFuzzy(got[:], tt.want[:]) {
+				t.Errorf("ColorToNLRGBA(%#+v) = %v, want %v", tt.input, got, tt.want)
+			}
+		})
+	}
+}
+
+func TestNLRGBAtoRGBA(t *testing.T) {
+	tests := []struct {
+		input [4]float64
+		want  [4]uint32
+	}{
+		{
+			[4]float64{0, .5, 1, 0},
+			[4]uint32{0, 0, 0, 0},
+		}, {
+			[4]float64{.25, .5, .75, 1},
+			[4]uint32{Delinearize(.25), Delinearize(.5), Delinearize(.75), 0xffff},
+		}, {
+			[4]float64{.25, .5, .75, .5},
+			[4]uint32{Delinearize(.25) / 2, Delinearize(.5) / 2, Delinearize(.75) / 2, 0x8000},
+		},
+	}
+	for _, tt := range tests {
+		t.Run(fmt.Sprintf("%v", tt.input), func(t *testing.T) {
+			if got := collect4(NLRGBAtoRGBA(tt.input[0], tt.input[1], tt.input[2], tt.input[3])); got != tt.want {
+				t.Errorf("NLRGBAtoRGBA(%v) = %v, want %v", tt.input, got, tt.want)
+			}
+		})
+	}
+}

internal/helper/test.go

@@ -0,0 +1,6 @@
+package helper
+
+type tester[T any] interface {
+	Errorf(string, ...any)
+	Run(name string, f func(T)) bool
+}

internal/helper/test_test.go

@@ -0,0 +1,263 @@
+package helper
+
+import (
+	"fmt"
+	"strings"
+	"sync"
+	"testing"
+)
+
+type testStatus struct {
+	parent *testStatus
+
+	m           sync.Mutex
+	errors      []string
+	childFailed bool
+	panicValue  any
+	handled     bool
+}
+
+func (s *testStatus) setChildFailed() {
+	// the status object is allowed to be nil, to simplify
+	// marking a child as failed when it doesn't have a parent.
+	if s == nil {
+		return
+	}
+
+	// propagate to potential parents first, to avoid
+	// having multiple locks simultaneous.
+	s.parent.setChildFailed()
+
+	s.m.Lock()
+	defer s.m.Unlock()
+	s.childFailed = true
+}
+
+func (s *testStatus) addError(msg string) {
+	s.parent.setChildFailed()
+
+	s.m.Lock()
+	defer s.m.Unlock()
+	s.errors = append(s.errors, msg)
+}
+
+func (s *testStatus) hasError(text string) bool {
+	s.m.Lock()
+	defer s.m.Unlock()
+
+	for _, msg := range s.errors {
+		if strings.Contains(msg, text) {
+			return true
+		}
+	}
+	return false
+}
+
+func (s *testStatus) setPanic(v any) {
+	s.parent.setChildFailed()
+
+	s.m.Lock()
+	defer s.m.Unlock()
+	s.panicValue = v
+}
+
+func (s *testStatus) getPanic() any {
+	s.m.Lock()
+	defer s.m.Unlock()
+	return s.panicValue
+}
+
+func (s *testStatus) setHandled() {
+	s.m.Lock()
+	defer s.m.Unlock()
+	s.handled = true
+}
+
+func (s *testStatus) hasFailed() bool {
+	s.m.Lock()
+	defer s.m.Unlock()
+	return s.childFailed || s.panicValue != nil || len(s.errors) > 0
+}
+
+func (s *testStatus) hasFailedChildren() bool {
+	s.m.Lock()
+	defer s.m.Unlock()
+	return s.childFailed
+}
+
+func (s *testStatus) wasHandled() bool {
+	s.m.Lock()
+	defer s.m.Unlock()
+	return s.handled
+}
+
+func (s *testStatus) log(t *testing.T, name string) {
+	s.m.Lock()
+	defer s.m.Unlock()
+
+	success := true
+
+	for _, msg := range s.errors {
+		t.Logf("%s: %s", name, msg)
+		success = false
+	}
+
+	if s.panicValue != nil {
+		t.Logf("%s: panic: %v", name, s.panicValue)
+		success = false
+	}
+
+	if s.childFailed {
+		t.Logf("%s: has failed children", name)
+		success = false
+	}
+
+	if success {
+		t.Logf("%s: success", name)
+	}
+}
+
+type mockTest struct {
+	*mockTester
+	*testStatus
+	name string
+}
+
+func (m *mockTest) run(f func(*mockTest)) (success bool) {
+	defer func() {
+		if r := recover(); r != nil {
+			m.setPanic(r)
+		}
+		success = !m.hasFailed()
+	}()
+
+	f(m)
+	return
+}
+
+func (m *mockTest) Errorf(f string, args ...any) {
+	m.addError(fmt.Sprintf(f, args...))
+}
+
+func (m *mockTest) Run(name string, f func(*mockTest)) bool {
+	child := &mockTest{
+		mockTester: m.mockTester,
+		testStatus: m.create(m.testStatus, m.name+"/"+name),
+		name:       m.name + "/" + name,
+	}
+	return child.run(f)
+}
+
+type mockTester struct {
+	t       *testing.T
+	m       sync.Mutex
+	results map[string]*testStatus
+}
+
+func (m *mockTester) create(parent *testStatus, name string) *testStatus {
+	m.m.Lock()
+	defer m.m.Unlock()
+
+	if _, ok := m.results[name]; ok {
+		m.t.Fatalf("%s: test already exists", name)
+		return nil
+	}
+
+	if m.results == nil {
+		m.results = make(map[string]*testStatus)
+	}
+
+	s := &testStatus{parent: parent}
+	m.results[name] = s
+	return s
+}
+
+func (m *mockTester) get(name string) *testStatus {
+	m.m.Lock()
+	defer m.m.Unlock()
+
+	if s, ok := m.results[name]; ok {
+		return s
+	}
+
+	m.t.Fatalf("%s: test doesn't exist", name)
+	return nil
+}
+
+// run the test.
+func (m *mockTester) run(name string, f func(*mockTest)) bool {
+	t := &mockTest{
+		mockTester: m,
+		testStatus: m.create(nil, name),
+		name:       name,
+	}
+	return t.run(f)
+}
+
+// panic if the named test doesn't exist or doesn't have an error containing the given text.
+func (m *mockTester) expectError(name string, text string) {
+	if s := m.get(name); s != nil {
+		s.setHandled()
+
+		if !s.hasError(text) {
+			m.t.Errorf("%s: doesn't contain error message: %s", name, text)
+			s.log(m.t, name)
+		}
+	}
+}
+
+func (m *mockTester) expectFailedChildren(name string) {
+	if s := m.get(name); s != nil {
+		s.setHandled()
+
+		if !s.hasFailedChildren() {
+			m.t.Errorf("%s: doesn't have failed children", name)
+			s.log(m.t, name)
+		}
+	}
+}
+
+// panic if the named test doesn't exist or didn't panic (and returns the panic value)
+func (m *mockTester) expectPanic(name string) any {
+	if s := m.get(name); s != nil {
+		s.setHandled()
+
+		if r := s.getPanic(); r != nil {
+			return r
+		}
+
+		m.t.Errorf("%s: didn't panic", name)
+		s.log(m.t, name)
+	}
+
+	return nil
+}
+
+// panic if the named test doesn't exist or has failed.
+func (m *mockTester) expectSuccess(name string) {
+	if s := m.get(name); s != nil {
+		s.setHandled()
+
+		if s.hasFailed() {
+			m.t.Errorf("%s: failed", name)
+			s.log(m.t, name)
+		}
+	}
+}
+
+func (m *mockTester) expectAllHandled() {
+	m.m.Lock()
+	defer m.m.Unlock()
+
+	if len(m.results) == 0 {
+		m.t.Errorf("no tests were run")
+		return
+	}
+
+	for name, s := range m.results {
+		if !s.wasHandled() {
+			m.t.Errorf("%s: not handled", name)
+			s.log(m.t, name)
+		}
+	}
+}

internal/helper/xyz.go

@@ -0,0 +1,27 @@
+package helper
+
+func LRGBtoXYZ(r, g, b float64) (_, _, _ float64) {
+	// https://en.wikipedia.org/wiki/SRGB#Correspondence_to_CIE_XYZ_stimulus
+	//
+	// Wikipedia lists this matrix for converting from linear sRGB to D65 CIE XYZ, so
+	// I'm considering it canonical:
+	//
+	// [+0.4124, +0.3576, +0.1805]
+	// [+0.2126, +0.7152, +0.0722]
+	// [+0.0193, +0.1192, +0.9505]
+	//
+	// The inverse:
+	// [+3.2406254773200531456132481428905,   -1.5372079722103185962799221761846,  -0.49862859869824785916021137156360 ]
+	// [-0.96893071472931930204316125127115,  +1.8757560608852411526964057125165,  +0.041517523842953942971183706902422]
+	// [+0.055710120445510610303218445022341, -0.20402105059848668752573283843409, +1.0569959422543882942447416955375  ]
+
+	return 0.4124*r + 0.3576*g + 0.1805*b,
+		0.2126*r + 0.7152*g + 0.0722*b,
+		0.0193*r + 0.1192*g + 0.9505*b
+}
+
+func XYZtoLRGB(x, y, z float64) (_, _, _ float64) {
+	return 3.2406254773200531456132481428905*x - 1.5372079722103185962799221761846*y - 0.49862859869824785916021137156360*z,
+		-0.96893071472931930204316125127115*x + 1.8757560608852411526964057125165*y + 0.041517523842953942971183706902422*z,
+		0.055710120445510610303218445022341*x - 0.20402105059848668752573283843409*y + 1.0569959422543882942447416955375*z
+}

internal/helper/xyz_test.go

@@ -0,0 +1,15 @@
+package helper
+
+import (
+	"testing"
+)
+
+func TestXYZtoLRGB(t *testing.T) {
+	for c := range Enum(false, true) {
+		want := collect3(RGBtoLRGB(uint32(c.R), uint32(c.G), uint32(c.B)))
+		if got := collect3(XYZtoLRGB(LRGBtoXYZ(want[0], want[1], want[2]))); !EqFloat64SliceFuzzy(want[:], got[:]) {
+			t.Errorf("XYZtoLRGB(LRGBtoXYZ(%v)) = %v, want unchanged", want, got)
+			return
+		}
+	}
+}