From d62384658de5ac0969733f7c590a2f27759cf836 Mon Sep 17 00:00:00 2001 From: mineracks <134782215+mineracks@users.noreply.github.com> Date: Thu, 28 May 2026 18:29:29 +1000 Subject: [PATCH] Lift bezier/ + bspline/ from Gangleri42 fork First batch of code lifts: hardware-agnostic curve math packages. - bezier/: cubic + quadratic Bezier curve representation + tessellation. Lifted verbatim from Gangleri42/seedhammer @ 0a3c63ef path bezier/. Stdlib only. - bspline/: B-spline curves + an LP-based optimiser that finds the smallest engrave-stroke representation of a path. Lifted from Gangleri42/seedhammer @ 0a3c63ef path bspline/. Depends on bezier/ + gonum (BSD-3, brought in transitively). Import paths rewritten from `seedhammer.com/X` to `github.com/mineracks/seedhammer-v1-companion/X`. No other changes. go build ./... and go test ./bezier/... ./bspline/... both clean. Doc.go stubs updated to reflect LIFTED status with the source commit pinned. Deferred to next commit: - backup/ (needs bc/fountain + bc/ur transitively for UR codes) - font/{comfortaa,poppins,constant}/ (needs font/bitmap runtime) Co-Authored-By: Claude Opus 4.7 (1M context) --- bezier/bezier.go | 489 ++++++++++++++++++++++++++++++++ bezier/bezier_test.go | 168 +++++++++++ bezier/doc.go | 9 +- bspline/bspline.go | 212 ++++++++++++++ bspline/bspline_test.go | 472 +++++++++++++++++++++++++++++++ bspline/doc.go | 12 +- bspline/optimize.go | 535 ++++++++++++++++++++++++++++++++++++ bspline/testdata/curves.png | Bin 0 -> 655 bytes go.mod | 4 +- go.sum | 2 + 10 files changed, 1894 insertions(+), 9 deletions(-) create mode 100644 bezier/bezier.go create mode 100644 bezier/bezier_test.go create mode 100644 bspline/bspline.go create mode 100644 bspline/bspline_test.go create mode 100644 bspline/optimize.go create mode 100644 bspline/testdata/curves.png create mode 100644 go.sum diff --git a/bezier/bezier.go b/bezier/bezier.go new file mode 100644 index 0000000..6e5caa7 --- /dev/null +++ b/bezier/bezier.go @@ -0,0 +1,489 @@ +package bezier + +import ( + "math" +) + +type Cubic struct { + C0, C1, C2, C3 Point +} + +type cubic16 struct { + C0, C1, C2, C3 point16 +} + +// Interpolator interpolates a [Cubic] curve. +type Interpolator struct { + int fraction + c cubic16 + // off tracks the offset from the current + // curve segment to c, which is offset to + // fill out its limited precision. + off Point + + // The remaining part of the segment. + rem struct { + c Cubic + steps uint + } +} + +const ( + logOne = 16 + // 1.0 in fixed-point for [Bezier.Interpolate16]. + one = 1 << logOne +) + +// fraction implements integer interpolation of the interval ]0;1]. +type fraction struct { + // lim is the interval end. + lim uint32 + // n is the number of ticks. + n uint32 + // 2*d/m2 is the fractional step rate. + d, m2 uint32 + // frac2 accumulates the fractional steps, + // multiplied by 2. + frac2 uint32 + // p1 is 1-p, where p is the interpolation parameter. + p1 uint32 +} + +// Construct a new fraction that steps through ]0;n] in d steps. +func newFraction(n, d uint32) fraction { + f := fraction{ + lim: n, + n: d, + // Round to nearest step by initializing the fraction + // to 1/2*steps (=> 2*frac = steps). + frac2: d, + } + if f.n > 0 { + f.p1 = n + f.d = n / d + f.m2 = (n % d) * 2 + } + return f +} + +func (f *fraction) Value() uint32 { + return f.lim - f.p1 +} + +func (f *fraction) Step() bool { + if f.p1 == 0 { + return false + } + f.p1 -= f.d + f.frac2 += f.m2 + if n2 := 2 * f.n; f.frac2 >= n2 { + f.frac2 -= n2 + f.p1-- + } + return true +} + +func (in *Interpolator) Segment(c Cubic, steps uint) { + // Track discontinuities. + end := in.c.C3.P().Add(in.off) + diff := c.C0.Sub(end) + in.off = in.off.Add(diff) + in.rem.c, in.rem.steps = c, steps +} + +// Step the interpolator or return false if there are +// no more steps. +func (in *Interpolator) Step() bool { + for !in.int.Step() { + if in.rem.steps == 0 { + return false + } + in.advance() + } + return true +} + +// Position of the bézier interpolation. +func (in *Interpolator) Position() Point { + p := in.int.Value() + return in.c.Interpolate16(p).P().Add(in.off) +} + +func (in *Interpolator) advance() { + c, steps := in.rem.c, in.rem.steps + if steps == 0 { + panic("overstepping") + } + c1, c2, t1 := splitBezier16(c, steps) + t2 := steps - t1 + in.rem.c = c2 + in.rem.steps = t2 + // off is the offset caused by centering of the + // 16-bit segment c16. + off := in.c.C3.P().Sub(c1.C0.P()) + in.off = in.off.Add(off) + in.c = c1 + in.int = newFraction(one, uint32(t1)) +} + +// Split a part off c that fit in a [bezier16]. +func splitBezier16(c Cubic, ticks uint) (cubic16, Cubic, uint) { + // Firt, compute the split that makes the duration + // fit 16 bits. + splitDiv := (ticks + one - 1) / one + // Then split until the bezier control points fit. + for { + d := uint32(splitDiv) + p := (one + d - 1) / d + c1, c2 := c.split16(p) + // Center the bezier to maximize use of the available + // precision. + if c16, ok := clampBezier16(centerBezier(c1)); ok { + return c16, c2, ticks / splitDiv + } + splitDiv *= 2 + } +} + +func centerBezier(c Cubic) Cubic { + b := Bounds{ + Min: Point{ + X: min(c.C0.X, c.C1.X, c.C2.X, c.C3.X), + Y: min(c.C0.Y, c.C1.Y, c.C2.Y, c.C3.Y), + }, + Max: Point{ + X: max(c.C0.X, c.C1.X, c.C2.X, c.C3.X), + Y: max(c.C0.Y, c.C1.Y, c.C2.Y, c.C3.Y), + }, + } + center := b.Max.Add(b.Min).Div(2) + return c.Sub(center) +} + +func clampBezier16(c Cubic) (cubic16, bool) { + c16 := cubic16{ + C0: clampP16(c.C0), + C1: clampP16(c.C1), + C2: clampP16(c.C2), + C3: clampP16(c.C3), + } + if c16.B() != c { + return cubic16{}, false + } + return c16, true +} + +func clampP16(p Point) point16 { + return point16{ + X: clamp16(p.X), + Y: clamp16(p.Y), + } +} + +func clamp16(v int) int16 { + return int16(min(max(v, math.MinInt16), math.MaxInt16)) +} + +func (b *cubic16) B() Cubic { + return Cubic{ + C0: b.C0.P(), + C1: b.C1.P(), + C2: b.C2.P(), + C3: b.C3.P(), + } +} + +func (p point16) P() Point { + return Point{ + X: int(p.X), + Y: int(p.Y), + } +} + +// Sample the curve at p. p must be in the interval [0;one]. +// The result is exact. +func (s *cubic16) Interpolate16(p16 uint32) point16 { + // p fit in 16 bits, except for the extremes. + switch { + case p16 == 0: + return s.C0 + case p16 == one: + return s.C3 + } + p := uint16(p16) + p1 := uint16(one - p16) + q0 := interpolate16(p, p1, s.C0, s.C1) + q1 := interpolate16(p, p1, s.C1, s.C2) + q2 := interpolate16(p, p1, s.C2, s.C3) + r0 := interpolate32(p, p1, q0, q1) + r1 := interpolate32(p, p1, q1, q2) + x := interpolate64(p, p1, r0, r1) + + const ( + prec = 3 * logOne + rounding = 1 << (prec - 1) + ) + return point16{ + X: int16((x.X + rounding) >> prec), + Y: int16((x.Y + rounding) >> prec), + } +} + +// split16 splits a bezier at p16, in the iterval [0;one]. +func (s *Cubic) split16(p16 uint32) (Cubic, Cubic) { + switch { + case p16 == 0: + return Cubic{}, *s + case p16 == one: + return *s, Cubic{} + } + p := uint16(p16) + p1 := uint16(one - p16) + q0 := interpolate32(p, p1, s.C0, s.C1) + q1 := interpolate32(p, p1, s.C1, s.C2) + q2 := interpolate32(p, p1, s.C2, s.C3) + r064 := interpolate64(p, p1, q0, q1) + r164 := interpolate64(p, p1, q1, q2) + // We're out of bits - round and shift down once. + r0 := roundP64(r064) + r1 := roundP64(r164) + // Round twice. + x64 := interpolate64(p, p1, r0, r1) + x := Point{ + X: int((x64.X + 1<<(2*logOne-1)) >> (2 * logOne)), + Y: int((x64.Y + 1<<(2*logOne-1)) >> (2 * logOne)), + } + c11 := roundP64(q0).Point() + c12 := roundP64(r0).Point() + c21 := roundP64(r1).Point() + c22 := roundP64(q2).Point() + return Cubic{s.C0, c11, c12, x}, + Cubic{x, c21, c22, s.C3} +} + +// roundP64 performs a 16-bit rounding shifting. +func roundP64(p Point64) Point64 { + return Point64{ + X: (p.X + 1<<(logOne-1)) >> logOne, + Y: (p.Y + 1<<(logOne-1)) >> logOne, + } +} + +type Point struct { + X, Y int +} + +type point16 struct { + X, Y int16 +} + +type Point64 struct { + X, Y int64 +} + +func P64(p Point) Point64 { + return Point64{ + X: int64(p.X), + Y: int64(p.Y), + } +} + +func (p Point64) Mul(s int) Point64 { + return Point64{ + X: p.X * int64(s), + Y: p.Y * int64(s), + } +} + +func (p Point64) Div(s int) Point64 { + return Point64{ + X: p.X / int64(s), + Y: p.Y / int64(s), + } +} + +func (p Point64) Add(p2 Point64) Point64 { + return Point64{ + X: p.X + p2.X, + Y: p.Y + p2.Y, + } +} + +func (p Point64) Point() Point { + return Point{ + X: int(p.X), + Y: int(p.Y), + } +} + +func Pt(x, y int) Point { + return Point{ + X: x, + Y: y, + } +} + +func (p Point) Mul16(s uint16) Point { + return Point{ + X: p.X * int(s), + Y: p.Y * int(s), + } +} + +func (p Point) Mul(s int) Point { + return Point{ + X: p.X * s, + Y: p.Y * s, + } +} + +func (p Point) Div(s int) Point { + return Point{ + X: p.X / s, + Y: p.Y / s, + } +} + +func (p Point) Add(p2 Point) Point { + return Point{ + X: p.X + p2.X, + Y: p.Y + p2.Y, + } +} + +func (p Point) Sub(p2 Point) Point { + return Point{ + X: p.X - p2.X, + Y: p.Y - p2.Y, + } +} + +func interpolate16(p, p1 uint16, a, b point16) Point { + b1 := b.P().Mul16(p) + a1 := a.P().Mul16(p1) + s := b1.Add(a1) + return s +} + +func interpolate32(p, p1 uint16, a, b Point) Point64 { + b1 := P64(b).Mul(int(p)) + a1 := P64(a).Mul(int(p1)) + s := b1.Add(a1) + return s +} + +func interpolate64(p, p1 uint16, a, b Point64) Point64 { + b1 := b.Mul(int(p)) + a1 := a.Mul(int(p1)) + s := b1.Add(a1) + return s +} + +func (s *Cubic) Add(off Point) Cubic { + return Cubic{ + s.C0.Add(off), + s.C1.Add(off), + s.C2.Add(off), + s.C3.Add(off), + } +} + +func (s *Cubic) Sub(off Point) Cubic { + return Cubic{ + s.C0.Sub(off), + s.C1.Sub(off), + s.C2.Sub(off), + s.C3.Sub(off), + } +} + +// Bounds is like [image.Rectangle] with its upper +// bound inclusive. +type Bounds struct { + Min, Max Point +} + +func (b Bounds) In(b2 Bounds) bool { + return inBounds(b.Min, b2) && inBounds(b.Max, b2) +} + +func (b Bounds) Union(b2 Bounds) Bounds { + return Bounds{ + Min: Point{ + X: min(b.Min.X, b2.Min.X), + Y: min(b.Min.Y, b2.Min.Y), + }, + Max: Point{ + X: max(b.Max.X, b2.Max.X), + Y: max(b.Max.Y, b2.Max.Y), + }, + } +} + +func (b Bounds) Empty() bool { + return b.Max.X < b.Min.X || b.Max.Y < b.Min.X +} + +func (b Bounds) Dx() int { + return int(b.Max.X - b.Min.X) +} + +func (b Bounds) Dy() int { + return int(b.Max.Y - b.Min.Y) +} + +func inBounds(p Point, b Bounds) bool { + return b.Min.X <= p.X && p.X <= b.Max.X && + b.Min.Y <= p.Y && p.Y <= b.Max.Y +} + +// Sample samples enough points on b that chords between samples +// are close to spacing apart. The samples are appended to points. +func Sample(points []Point, b Cubic, spacing int) []Point { + // Estimate the curve length by many small samples. + const samplingRate = 200 + var totalDist int + var first Point + if len(points) > 0 { + first = points[len(points)-1] + } + prev := first + in := new(Interpolator) + in.Segment(b, samplingRate) + for in.Step() { + s := in.Position() + totalDist += dist(prev, s) + prev = s + } + // Given the total distance, compute the number of samples. + nsamples := (totalDist + spacing - 1) / spacing + // Sample short segments in the middle. + nsamples = max(nsamples, 2) + adjSpacing := (totalDist + nsamples - 1) / nsamples + prev = first + var d int + pi := 0 + // Sample inner points. + in = new(Interpolator) + in.Segment(b, samplingRate) + for range nsamples - 1 { + var s Point + for d < adjSpacing { + pi++ + in.Step() + s = in.Position() + d += dist(prev, s) + prev = s + } + points = append(points, s) + d -= adjSpacing + } + // Force endpoint to align with curve segment end. + points = append(points, b.C3) + return points +} + +func dist(a, b Point) int { + d := b.Sub(a) + return int(math.Round(math.Hypot(float64(d.X), float64(d.Y)))) +} diff --git a/bezier/bezier_test.go b/bezier/bezier_test.go new file mode 100644 index 0000000..c1dccca --- /dev/null +++ b/bezier/bezier_test.go @@ -0,0 +1,168 @@ +package bezier + +import ( + "bytes" + "fmt" + "image" + "image/png" + "math" + "os" + "testing" +) + +func TestBezierAccuracy(t *testing.T) { + curves := []struct { + steps uint + b Cubic + }{ + { + 671998, + Cubic{Point{213333, 170666}, Point{170666, 192000}, Point{277333, 192000}, Point{319999, 170666}}, + }, + { + 0, + Cubic{Point{0, 0}, Point{0, 0}, Point{0, 0}, Point{0, -71}}, + }, + { + one, + Cubic{Point{math.MinInt16, math.MinInt16}, Point{0, 0}, Point{0, 0}, Point{math.MaxInt16, math.MaxInt16}}, + }, + } + for _, test := range curves { + verifyBezierSmoothness(t, test.b, test.steps) + } +} + +func FuzzInterpolatorAccuracy(f *testing.F) { + f.Fuzz(func(t *testing.T, c0x, c0y, c1x, c1y, c2x, c2y, c3x, c3y int, ticks int32) { + b := Cubic{ + C0: Pt(c0x, c0y), + C1: Pt(c1x, c1y), + C2: Pt(c2x, c2y), + C3: Pt(c3x, c3y), + } + verifyBezierSmoothness(t, b, uint(ticks)) + }) +} + +func verifyBezierSmoothness(t *testing.T, b Cubic, ticks uint) { + t.Helper() + var dir [2]int + pos := b.C0 + if ticks == 0 { + pos = b.C3 + } + errors := 10 + step := 0 + reportErr := func(f string, args ...any) { + t.Helper() + if errors == 0 { + return + } + err := fmt.Errorf(f, args...) + t.Errorf("step %d/%d of %v: %v", step, ticks, b, err) + errors-- + } + nturns := [2]int{} + in := new(Interpolator) + in.Segment(b, ticks) + for in.Step() { + step++ + newPos := in.Position() + step := newPos.Sub(pos) + pos = newPos + step.X = max(min(step.X, 1), -1) + step.Y = max(min(step.Y, 1), -1) + + for axis, step := range []int{step.X, step.Y} { + if step == 0 { + continue + } + // Detect turns. + if dir[axis] != 0 && step != dir[axis] { + nturns[axis]++ + } + dir[axis] = step + } + } + for axis, nturns := range nturns { + // A third degree bézier has at most 2 turns. + if nturns > 2 { + reportErr("axis %d: %d turns, expected <= 2", axis, nturns) + } + } + if end := b.C3; pos != end { + reportErr("ended in %v, expected %v", pos, end) + } +} + +// quadBezier represents a second degree Bézier curve. +type quadBezier struct { + C0, C1, C2 Point +} + +func velocityCurve(b Cubic) quadBezier { + // Derived quadratic velocity curve: + // + // 3{(1−t)²(C1−C0)+2t(1−t)(C2−C1)+t²(C3−C2)} + return quadBezier{ + C0: b.C1.Sub(b.C0).Mul(3), + C1: b.C2.Sub(b.C1).Mul(3), + C2: b.C3.Sub(b.C2).Mul(3), + } +} + +func bezierMaxVelocity(b Cubic) uint { + c := velocityCurve(b) + v0, v1, v2 := c.C0, c.C1, c.C2 + return uint(max( + v0.X, v1.X, v2.X, + -v0.X, -v1.X, -v2.X, + v0.Y, v1.Y, v2.Y, + -v0.Y, -v1.Y, -v2.Y, + )) +} + +func compareImages(imgPath string, update bool, img image.Image) error { + if update { + var buf bytes.Buffer + if err := png.Encode(&buf, img); err != nil { + return err + } + return os.WriteFile(imgPath, buf.Bytes(), 0o640) + } + f, err := os.Open(imgPath) + if err != nil { + return err + } + want, _, err := image.Decode(f) + if err != nil { + return err + } + f.Close() + if w, g := want.Bounds().Size(), img.Bounds().Size(); w != g { + return fmt.Errorf("golden image bounds mismatch: got %v, want %v", g, w) + } + mismatches := 0 + pixels := 0 + width, height := want.Bounds().Dx(), want.Bounds().Dy() + gotOff := img.Bounds().Min + for y := range height { + for x := range width { + wanta, _, _, _ := want.At(x, y).RGBA() + want := wanta != 0 + gota, _, _, _ := img.At(gotOff.X+x, gotOff.Y+y).RGBA() + got := gota != 0 + if want { + pixels++ + } + if got != want { + mismatches++ + } + } + } + if mismatches > 0 { + return fmt.Errorf("%d/%d pixels golden image mismatches", mismatches, pixels) + } + return nil +} diff --git a/bezier/doc.go b/bezier/doc.go index 51eabdf..3cc93df 100644 --- a/bezier/doc.go +++ b/bezier/doc.go @@ -1,10 +1,11 @@ -// Package bezier provides quadratic and cubic Bezier curve tessellation -// for the engrave pipeline. +// Package bezier provides quadratic and cubic Bezier curve representations +// plus tessellation helpers for the engrave pipeline. // // Used to convert OpenType glyph QuadTo/CubeTo segments and SH1E SVG path // Q/C commands into linear MoveTo/LineTo sequences at the engraver's // resolution. // -// Status: STUB — to be lifted from upstream at v1.3.0 + cross-checked -// against Gangleri42's fork (the math is hardware-agnostic and unchanged). +// LIFTED from https://github.com/Gangleri42/seedhammer/tree/seedhammer-features/bezier +// at commit 0a3c63efb125d17d8ec86ce739ecd058c8747cfe. Hardware-agnostic +// math; identical to upstream apart from the Go import path. package bezier diff --git a/bspline/bspline.go b/bspline/bspline.go new file mode 100644 index 0000000..69c7d22 --- /dev/null +++ b/bspline/bspline.go @@ -0,0 +1,212 @@ +package bspline + +import ( + "iter" + "math" + + "github.com/mineracks/seedhammer-v1-companion/bezier" +) + +type segmentBounds struct { + Knots [4]Knot +} + +type Segment struct { + // The previous segment end point. + prev bezier.Point + // Knots is a sliding window of knots. + Knots [3]Knot +} + +// Curve is an iterator over the knots of a b-spline. +type Curve = iter.Seq[Knot] + +type Knot struct { + Ctrl bezier.Point + T uint + Engrave bool +} + +func (s *segmentBounds) Knot(k Knot) (Bounds, bool) { + c0, c1, c2, c3 := s.Knots[0].Ctrl, s.Knots[1].Ctrl, s.Knots[2].Ctrl, s.Knots[3].Ctrl + segKnot := s.Knots[2] + copy(s.Knots[:3], s.Knots[1:]) + s.Knots[3] = k + return Bounds{ + Min: bezier.Point{ + X: min(c0.X, c1.X, c2.X, c3.X), + Y: min(c0.Y, c1.Y, c2.Y, c3.Y), + }, + Max: bezier.Point{ + X: max(c0.X, c1.X, c2.X, c3.X), + Y: max(c0.Y, c1.Y, c2.Y, c3.Y), + }, + }, segKnot.Engrave +} + +func (s *Segment) Knot(k Knot) (bezier.Cubic, uint, bool) { + // Extract Bézier curve through Böhm's algorithm. + // See "An Introduction to B-Spline Curves" by Thomas W. Sederberg. + dt2, dt3, dt4, dt5 := s.Knots[0].T, s.Knots[1].T, s.Knots[2].T, k.T + P234, P345, P456 := s.Knots[0].Ctrl, s.Knots[1].Ctrl, s.Knots[2].Ctrl + engrave := s.Knots[1].Engrave + + // Shift the knot window. + copy(s.Knots[:2], s.Knots[1:]) + s.Knots[2] = k + // The start point of this segment is the end point of the previous. + P333 := s.prev + + if dt3 == 0 { + s.prev = P345 + // Zero duration curve. + return bezier.Cubic{ + C0: P234, + C1: P234, + C2: P345, + C3: P345, + }, 0, engrave + } + d1 := int(dt4 + dt3 + dt2) + P334 := bezier.P64(P234).Mul(int(dt4 + dt3)).Add(bezier.P64(P345).Mul(int(dt2))).Div(d1).Point() + P344 := bezier.P64(P234).Mul(int(dt4)).Add(bezier.P64(P345).Mul(int(dt3 + dt2))).Div(d1).Point() + d3 := int(dt5 + dt4 + dt3) + P445 := bezier.P64(P345).Mul(int(dt5 + dt4)).Add(bezier.P64(P456).Mul(int(dt3))).Div(d3).Point() + d5 := int(dt4 + dt3) + P444 := bezier.P64(P344).Mul(int(dt4)).Add(bezier.P64(P445).Mul(int(dt3))).Div(d5).Point() + s.prev = P444 + return bezier.Cubic{ + C0: P333, + C1: P334, + C2: P344, + C3: P444, + }, dt3, engrave +} + +// Kinematics track the derivative values of a B-Spline. +type Kinematics struct { + Velocity bezier.Point + Acceleration bezier.Point + Jerk bezier.Point + + // Sliding window of knots and state. + knots [3]uint + p [2]bezier.Point + v bezier.Point + a bezier.Point +} + +// Knot shifts the spline one knot and updates the kinematics. +func (k *Kinematics) Knot(t uint, ctrl bezier.Point, scale uint) { + copy(k.knots[:], k.knots[1:]) + k.knots[2] = t + v := k.derive(k.p[0], k.p[1], 3, scale) + copy(k.p[:1], k.p[1:]) + k.p[1] = ctrl + a := k.derive(k.v, v, 2, scale) + k.v = v + j := k.derive(k.a, a, 1, scale) + k.a = a + k.Velocity = v + k.Acceleration = a + k.Jerk = j +} + +func (k *Kinematics) Max() (v, a, j uint) { + vv, aa, jj := k.Velocity, k.Acceleration, k.Jerk + return uint(max(vv.X, -vv.X, vv.Y, -vv.Y)), + uint(max(aa.X, -aa.X, aa.Y, -aa.Y)), + uint(max(jj.X, -jj.X, jj.Y, -jj.Y)) +} + +func (k *Kinematics) derive(p0, p1 bezier.Point, degree int, scale uint) bezier.Point { + t := uint(0) + knots := k.knots[:degree] + for _, k := range knots { + t += k + } + if t == 0 { + return bezier.Point{} + } + return bezier.P64(p1.Sub(p0)).Mul(degree * int(scale)).Div(int(t)).Point() +} + +// ComputeKinematics compute the absolute kinematic maximums of the +// clamped B-spline. +func ComputeKinematics(spline []Knot, scale uint) (v, a, j uint) { + var deriv Kinematics + var maxv, maxa, maxj uint + for _, p := range spline { + deriv.Knot(p.T, p.Ctrl, scale) + v, a, j := deriv.Max() + maxv = max(maxv, v) + maxa = max(maxa, a) + maxj = max(maxj, j) + } + return maxv, maxa, maxj +} + +type Attributes struct { + Bounds Bounds + Duration uint +} + +// Bounds is like [image.Rectangle] with its upper +// bound inclusive. +type Bounds struct { + Min, Max bezier.Point +} + +func (b Bounds) In(b2 Bounds) bool { + return inBounds(b.Min, b2) && inBounds(b.Max, b2) +} + +func (b Bounds) Union(b2 Bounds) Bounds { + return Bounds{ + Min: bezier.Point{ + X: min(b.Min.X, b2.Min.X), + Y: min(b.Min.Y, b2.Min.Y), + }, + Max: bezier.Point{ + X: max(b.Max.X, b2.Max.X), + Y: max(b.Max.Y, b2.Max.Y), + }, + } +} + +func (b Bounds) Empty() bool { + return b.Max.X < b.Min.X || b.Max.Y < b.Min.X +} + +func (b Bounds) Dx() int { + return int(b.Max.X - b.Min.X) +} + +func (b Bounds) Dy() int { + return int(b.Max.Y - b.Min.Y) +} + +func inBounds(p bezier.Point, b Bounds) bool { + return b.Min.X <= p.X && p.X <= b.Max.X && + b.Min.Y <= p.Y && p.Y <= b.Max.Y +} + +func Measure(spline Curve) Attributes { + inf := Bounds{ + Min: bezier.Point{X: math.MaxInt32, Y: math.MaxInt32}, + Max: bezier.Point{X: math.MinInt32, Y: math.MinInt32}, + } + bounds := inf + var t uint + var bsb segmentBounds + for c := range spline { + t += c.T + if bb, engrave := bsb.Knot(c); engrave { + bounds = bounds.Union(bb) + } + } + return Attributes{ + Bounds: bounds, + Duration: t, + } +} diff --git a/bspline/bspline_test.go b/bspline/bspline_test.go new file mode 100644 index 0000000..7db61c2 --- /dev/null +++ b/bspline/bspline_test.go @@ -0,0 +1,472 @@ +package bspline + +import ( + "bytes" + "flag" + "fmt" + "image" + "image/color" + "image/draw" + "image/png" + "math" + "os" + "path/filepath" + "slices" + "testing" + + "github.com/mineracks/seedhammer-v1-companion/bezier" +) + +var ( + update = flag.Bool("update", false, "update golden files") + dump = flag.String("dump", "", "dump SVG files to directory") +) + +func TestInterpolator(t *testing.T) { + splines := [][]bezier.Point{ + { + bezier.Pt(1, 1), + bezier.Pt(1, 1), + bezier.Pt(1, 1), + bezier.Pt(1, 1), + }, + { + bezier.Pt(10, 10), + bezier.Pt(10, 10), + bezier.Pt(50, 20), + bezier.Pt(50, 20), + }, + { + bezier.Pt(20, 40), + bezier.Pt(50, 100), + bezier.Pt(140, -50), + bezier.Pt(100, 30), + }, + { + bezier.Pt(50, 100), + bezier.Pt(200, 30), + bezier.Pt(40, 30), + bezier.Pt(140, 100), + }, + } + img := image.NewAlpha(image.Rectangle{Max: image.Pt(150, 150)}) + for _, s := range splines { + rasterizeUniformBSpline(img, s) + } + p := filepath.Join("testdata", "curves.png") + if err := compareImages(p, *update, img); err != nil { + t.Error(err) + } +} + +func TestMeasureBSpline(t *testing.T) { + tests := []struct { + spline []Knot + bounds Bounds + }{ + { + []Knot{ + {Ctrl: bezier.Pt(0, 0)}, + {Ctrl: bezier.Pt(0, 0)}, + {Ctrl: bezier.Pt(0, 0)}, + {T: 1, Ctrl: bezier.Pt(1, 0)}, + {T: 1, Ctrl: bezier.Pt(100, 10)}, + {T: 1, Ctrl: bezier.Pt(10, 30)}, + {T: 1, Ctrl: bezier.Pt(60, 30), Engrave: true}, + {T: 1, Ctrl: bezier.Pt(50, 10)}, + {T: 1, Ctrl: bezier.Pt(0, 0)}, + {Ctrl: bezier.Pt(0, 0)}, + {Ctrl: bezier.Pt(0, 0)}, + }, + Bounds{Min: bezier.Pt(10, 10), Max: bezier.Pt(100, 30)}, + }, + } + for _, test := range tests { + attrs := Measure(slices.Values(test.spline)) + if got := attrs.Bounds; got != test.bounds { + t.Errorf("measured spline bounds to %+v, want %+v", got, test.bounds) + } + } +} + +func TestKinematics(t *testing.T) { + tests := []struct { + spline []Knot + velocity []bezier.Point + acceleration []bezier.Point + jerk []bezier.Point + }{ + { + []Knot{ + {Ctrl: bezier.Pt(0, 100), T: 0}, + {Ctrl: bezier.Pt(0, 100), T: 0}, + {Ctrl: bezier.Pt(0, 100), T: 0}, + {Ctrl: bezier.Pt(0, 100), T: 1}, + {Ctrl: bezier.Pt(80, 100), T: 1}, + {Ctrl: bezier.Pt(80, 100), T: 1}, + {Ctrl: bezier.Pt(80, 100), T: 0}, + {Ctrl: bezier.Pt(80, 100), T: 0}, + }, + []bezier.Point{{X: 0, Y: 0}, {X: 0, Y: 0}, {X: 80, Y: 0}, {X: 0, Y: 0}, {X: 0, Y: 0}}, + []bezier.Point{{X: 0, Y: 0}, {X: 80, Y: 0}, {X: -80, Y: 0}, {X: 0, Y: 0}}, + []bezier.Point{{X: 80, Y: 0}, {X: -160, Y: 0}, {X: 80, Y: 0}}, + }, + { + []Knot{ + {Ctrl: bezier.Pt(10, 20), T: 0}, + {Ctrl: bezier.Pt(10, 20), T: 0}, + {Ctrl: bezier.Pt(10, 20), T: 0}, + {Ctrl: bezier.Pt(5, 30), T: 2}, + {Ctrl: bezier.Pt(50, 25), T: 3}, + {Ctrl: bezier.Pt(200, -100), T: 1}, + {Ctrl: bezier.Pt(-200, -200), T: 4}, + {Ctrl: bezier.Pt(80, 100), T: 2}, + {Ctrl: bezier.Pt(80, 100), T: 0}, + {Ctrl: bezier.Pt(80, 100), T: 0}, + }, + []bezier.Point{{0, 0}, {-3, 6}, {22, -2}, {56, -46}, {-171, -42}, {140, 150}, {0, 0}}, + []bezier.Point{{-3, 6}, {10, -3}, {17, -22}, {-90, 1}, {103, 64}, {-140, -150}}, + []bezier.Point{{6, -4}, {2, -6}, {-107, 23}, {48, 15}, {-121, -107}}, + }, + } + for _, test := range tests { + vel, accel, jerk := extractKinematics(t, test.spline) + if !slices.Equal(vel, test.velocity) { + t.Errorf("got velocity %v, expected %v", vel, test.velocity) + } + if !slices.Equal(accel, test.acceleration) { + t.Errorf("got acceleration %v, expected %v", accel, test.acceleration) + } + if !slices.Equal(jerk, test.jerk) { + t.Errorf("got jerk %v, expected %v", jerk, test.jerk) + } + } +} + +func TestInterpolatePoints(t *testing.T) { + tests := [][]bezier.Point{ + {{0, 0}, {1000, 1000}}, + {{100, 100}, {400, 100}, {700, 100}, {1000, 100}, {1000, 400}, {1000, 700}, {1000, 1000}, {700, 1000}, {400, 1000}, {100, 1000}, {100, 700}, {100, 400}, {100, 100}}, + {{X: 2900, Y: -2300}, {X: 2605, Y: -2506}, {X: 2269, Y: -2637}, {X: 1909, Y: -2676}, {X: 1550, Y: -2607}, {X: 1244, Y: -2417}, {X: 1100, Y: -2100}, {X: 1231, Y: -1709}, {X: 1585, Y: -1501}, {X: 2000, Y: -1400}, {X: 2342, Y: -1335}, {X: 2665, Y: -1231}, {X: 2916, Y: -1025}, {X: 3000, Y: -700}, {X: 2846, Y: -328}, {X: 2509, Y: -123}, {X: 2111, Y: -47}, {X: 1714, Y: -73}, {X: 1337, Y: -189}, {X: 1000, Y: -400}}, + } + for i, waypoints := range tests { + uspline, v, a, j, err := interpolatePoints(waypoints) + if err != nil { + t.Fatal(err) + } + cuspline := expandUniformBSpline(uspline) + if dir := *dump; dir != "" { + path := filepath.Join(dir, fmt.Sprintf("dump%d.svg", i)) + if err := dumpBSplineToSVG(path, cuspline, waypoints); err != nil { + t.Fatal(err) + } + } + var seg Segment + // Test that the spline passes through every line segment + // spanned by waypoints. + for i, k := range cuspline { + b, _, _ := seg.Knot(k) + if i < 6 { + // Skip starting points. + continue + } + p := b.C0 + w0, w1 := waypoints[i-6], waypoints[i-5] + if !closeToSegment(p, w0, w1) { + t.Errorf("spline passes through %v, which is not on the line segment %v-%v", p, w0, w1) + } + } + // Test that the optimization goal matches the kinematics. + wantv, wanta, wantj := extractKinematics(t, cuspline) + if !pointsNearlyEqual(v, wantv) { + t.Errorf("got velocities\n%v, expected\n%v", v, wantv) + } + if !pointsNearlyEqual(a, wanta) { + t.Errorf("got accelerations\n%v, expected\n%v", a, wanta) + } + if !pointsNearlyEqual(j, wantj) { + t.Errorf("got jerks\n%v, expected\n%v", j, wantj) + } + } +} + +func TestExprConst(t *testing.T) { + tests := []struct { + name string + want float64 + expr expr + }{ + {"0", 0, constExpr(0)}, + {"1", 1, constExpr(1)}, + {"0*10", 0, constExpr(0).Mul(10)}, + {"1*10", 10, constExpr(1).Mul(10)}, + {"10/2", 2, constExpr(10).Div(5)}, + {"10+2", 12, constExpr(10).Add(constExpr(2))}, + {"10-10", 0, constExpr(10).Sub(constExpr(10))}, + {"[1,1]-[0,1]", 1, varExpr(0).Add(varExpr(1)).Sub(varExpr(1))}, + {"[10,0,10,1]-[1,0,10,1]", 9, constExpr(10).Add(varExpr(2).Mul(10)).Add(varExpr(3)). + Sub(constExpr(1).Add(varExpr(2).Mul(10)).Add(varExpr(3)))}, + } + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + if got := test.expr; got.Const() != test.want { + t.Errorf("%s == %s", test.name, got) + } + }) + } +} + +func TestExprVar(t *testing.T) { + tests := []struct { + name string + want []float64 + expr expr + }{ + {"[1]", []float64{1}, varExpr(0)}, + {"[0,1]", []float64{0, 1}, varExpr(1)}, + {"[1,0]+[0,1]", []float64{1, 1}, varExpr(0).Add(varExpr(1))}, + {"[10,0]+[0,1]", []float64{10, 1}, constExpr(10).Add(varExpr(1))}, + {"[10,0,10,1]+[1,0,10,1]", []float64{11, 0, 20, 2}, constExpr(10).Add(varExpr(2).Mul(10)).Add(varExpr(3)). + Add(constExpr(1).Add(varExpr(2).Mul(10)).Add(varExpr(3)))}, + } + for _, test := range tests { + t.Run(test.name, func(t *testing.T) { + r := test.expr.Explode() + if got := test.expr; !slices.Equal(r, test.want) { + t.Errorf("%s == %s", test.name, got) + } + }) + } +} + +// closeToSegment checks if point p lies on the segment between a and b. +func closeToSegment(p, a, b bezier.Point) bool { + const slack = 2 + + // Compute distance from p to the line passing a and b. + nom := (b.Y-a.Y)*p.X - (b.X-a.X)*p.Y + b.X*a.Y - b.Y*a.X + nom = max(nom, -nom) + denom := math.Hypot(float64(b.Y-a.Y), float64(b.X-a.X)) + dist := int(math.Round(float64(nom) / denom)) + if dist > slack { + return false + } + + // p must also lie in the bounding box with a and b as corners. + return p.X >= min(a.X, b.X)-slack && p.X <= max(a.X, b.X)+slack && + p.Y >= min(a.Y, b.Y)-slack && p.Y <= max(a.Y, b.Y)+slack +} + +// expandUniformBSpline expands a uniform B-spline to its non-uniform +// representation. +func expandUniformBSpline(bspline []bezier.Point) []Knot { + knots := make([]Knot, 0, len(bspline)) + for i, c := range bspline { + k := Knot{ + Ctrl: c, + } + if i > 2 && i < len(bspline)-2 { + k.T = 1 + } + knots = append(knots, k) + } + return knots +} + +func uniformBSpline(uspline []bezier.Point) []Knot { + s := uspline[0] + e := uspline[len(uspline)-1] + // Clamp. + clamped := append([]bezier.Point{s, s}, uspline...) + clamped = append(clamped, e, e) + spline := expandUniformBSpline(clamped) + maxv, _, _ := ComputeKinematics(spline, 1) + for i := range spline[3 : len(spline)-2] { + spline[i+3].T = maxv + } + return spline +} + +func rasterizeUniformBSpline(img draw.Image, uspline []bezier.Point) { + spline := uniformBSpline(uspline) + var in bezier.Interpolator + dot := func(pos bezier.Point) { + img.Set(int(pos.X), int(pos.Y), color.Alpha{A: 255}) + } + dot(uspline[0]) + var seg Segment + for _, k := range spline { + c, ticks, _ := seg.Knot(k) + in.Segment(c, uint(ticks)) + for in.Step() { + dot(in.Position()) + } + } +} + +func compareImages(imgPath string, update bool, img image.Image) error { + if update { + var buf bytes.Buffer + if err := png.Encode(&buf, img); err != nil { + return err + } + return os.WriteFile(imgPath, buf.Bytes(), 0o640) + } + f, err := os.Open(imgPath) + if err != nil { + return err + } + want, _, err := image.Decode(f) + if err != nil { + return err + } + f.Close() + if w, g := want.Bounds().Size(), img.Bounds().Size(); w != g { + return fmt.Errorf("golden image bounds mismatch: got %v, want %v", g, w) + } + mismatches := 0 + pixels := 0 + width, height := want.Bounds().Dx(), want.Bounds().Dy() + gotOff := img.Bounds().Min + for y := range height { + for x := range width { + wanta, _, _, _ := want.At(x, y).RGBA() + want := wanta != 0 + gota, _, _, _ := img.At(gotOff.X+x, gotOff.Y+y).RGBA() + got := gota != 0 + if want { + pixels++ + } + if got != want { + mismatches++ + } + } + } + if mismatches > 0 { + return fmt.Errorf("%d/%d pixels golden image mismatches", mismatches, pixels) + } + return nil +} + +func dumpBSplineToSVG(filename string, spline []Knot, waypoints []bezier.Point) (err error) { + f, err := os.Create(filename) + if err != nil { + return err + } + defer func() { + if cerr := f.Close(); err == nil { + err = cerr + } + }() + + minX, minY, maxX, maxY := math.MaxInt32, math.MaxInt32, math.MinInt32, math.MinInt32 + for _, p := range spline { + minX = min(minX, p.Ctrl.X) + minY = min(minY, p.Ctrl.Y) + maxX = max(maxX, p.Ctrl.X) + maxY = max(maxY, p.Ctrl.Y) + } + for _, p := range waypoints { + minX = min(minX, p.X) + minY = min(minY, p.Y) + maxX = max(maxX, p.X) + maxY = max(maxY, p.Y) + } + const margin = 20 + fmt.Fprintf(f, "\n", + minX-margin, minY-margin, (maxX-minX)+2*margin, (maxY-minY)+2*margin) + + fmt.Fprintln(f, ``) + + fmt.Fprintf(f, ``, minX-margin, minY-margin, (maxX-minX)+2*margin, (maxY-minY)+2*margin) + + fmt.Fprint(f, ``) + + fmt.Fprint(f, ``) + + for _, p := range spline { + fmt.Fprintf(f, ``+"\n", p.Ctrl.X, p.Ctrl.Y) + } + for _, p := range waypoints { + fmt.Fprintf(f, ``+"\n", p.X, p.Y) + } + + fmt.Fprintln(f, "") + return nil +} + +func extractKinematics(t *testing.T, spline []Knot) (v, a, j []bezier.Point) { + var kin Kinematics + var vel, accel, jerk []bezier.Point + var maxv, maxa, maxj uint + for i, k := range spline { + kin.Knot(k.T, k.Ctrl, 1) + v := kin.Velocity + absv := uint(max(v.X, -v.X, v.Y, -v.Y)) + maxv = max(maxv, absv) + a := kin.Acceleration + absa := uint(max(a.X, -a.X, a.Y, -a.Y)) + maxa = max(maxa, absa) + j := kin.Jerk + absj := uint(max(j.X, -j.X, j.Y, -j.Y)) + maxj = max(maxj, absj) + mv, ma, mj := kin.Max() + if mv != absv || ma != absa || mj != absj { + t.Errorf("got absolute kinematics (%d, %d, %d), expected (%d, %d, %d)", mv, ma, mj, absv, absa, absj) + } + if i >= 3 { + vel = append(vel, v) + if i >= 4 { + accel = append(accel, a) + if i >= 5 { + jerk = append(jerk, j) + } + } + } + } + gotv, gota, gotj := ComputeKinematics(spline, 1) + if gotv != uint(maxv) || gota != uint(maxa) || gotj != uint(maxj) { + t.Errorf("got kinematics (v, a, j) = (%d, %d, %d), expected (%d, %d, %d)", v, a, j, maxv, maxa, maxj) + } + return vel, accel, jerk +} + +func pointsNearlyEqual(a, b []bezier.Point) bool { + const slack = 4 + if len(a) != len(b) { + return false + } + for i, pa := range a { + pb := b[i] + dx, dy := pa.X-pb.X, pa.Y-pb.Y + dist := dx*dx + dy*dy + if dist > slack*slack { + return false + } + } + return true +} diff --git a/bspline/doc.go b/bspline/doc.go index 9d3f488..84d30f7 100644 --- a/bspline/doc.go +++ b/bspline/doc.go @@ -1,7 +1,11 @@ -// Package bspline provides B-spline tessellation for the engrave pipeline. +// Package bspline provides B-spline curve tessellation + optimisation for +// the engrave pipeline. // -// Used by the OpenType font rasteriser for the few glyphs that use -// B-spline rather than Bezier segments. +// Used by the OpenType font rasteriser for glyphs that use B-spline rather +// than Bezier segments, and by the optimiser that finds the smallest +// engrave-stroke representation of a given path. // -// Status: STUB — to be lifted from upstream at v1.3.0. +// LIFTED from https://github.com/Gangleri42/seedhammer/tree/seedhammer-features/bspline +// at commit 0a3c63efb125d17d8ec86ce739ecd058c8747cfe. Hardware-agnostic. +// Depends on bezier/ + the gonum linear-programming library. package bspline diff --git a/bspline/optimize.go b/bspline/optimize.go new file mode 100644 index 0000000..0f162d5 --- /dev/null +++ b/bspline/optimize.go @@ -0,0 +1,535 @@ +package bspline + +import ( + "fmt" + "math" + + "gonum.org/v1/gonum/mat" + "gonum.org/v1/gonum/optimize/convex/lp" + "github.com/mineracks/seedhammer-v1-companion/bezier" +) + +// InterpolatePoints takes a clamped cubic uniform b-spline and +// returns a similar spline that minimizes maximum speed, +// acceleration and jerk. The returned spline has a control point +// in each control point interval of the input spline. +// As a special case, a zero-length input B-spline is returned as is. +// +// For example, the clamped uniform B-spline with n+4 control +// points, +// +// P0P0P0 - P1 - P2 - P3 - … - P{n-1}P{n-1}P{n-1} +// +// is turned into another clamped uniform B-spline n+5 control +// points +// +// Q0Q0Q0 - Q1 - Q2 - Q3 - Q4 - … - Q{n}Q{n}Q{n} +// +// where +// +// Q0 = P0, +// Q{n} = P{n-1}, and +// Q{i} for i∈[1,n-1] is on the line segment P{i-1} - P{i}. +func InterpolatePoints(pts []bezier.Point) (knots []bezier.Point, err error) { + knots, _, _, _, err = interpolatePoints(pts) + return +} + +func interpolatePoints(pts []bezier.Point) (knots, v, a, j []bezier.Point, err error) { + // Find the placement of the knots of the smooth B-spline, + // whose control points are on the line segments of the original. + // The placement should minimize the maximum of the kinetic + // properties velocity, acceleration and jerk. + // + // Because of the strong convex hull property of B-splines, + // and because the derivative of a B-spline is another B-spline, + // it suffices to minimize the kinetic properties at the control + // points of the B-splines that corresponds to the kinetic properties. + // + // Construct a linear program that minimizes J >= 0, which is the + // maximum of all properties of all control points of the input + // B-spline. + // + // The program is on the form + // + // minimize cᵀ x + // such that Ax = b + // x >= 0 . + // + // To confine the position of internal control points to line + // segments of the input spline, define each point Q{i} in terms of + // a scalar weight, w{i}∈]0,1[ such that + // + // Q{i} = (1-w{i})P{i}+w{i}P{i+1} + // = w{i}(P{i+1} - P{i}) + P{i} + // + // The weight interval is open to ensure unique control points. + // Any duplicate control point would violate the C² continuity + // of the cubic B-spline. + // + // To clamp the output spline, it must contain 3 duplicate control + // points at each end. Output control point lie between input points, + // so the input spline is extended by one point at each end to ensure 4 + // identical start and end points. + // Since there is an output control point per interval, the output spline + // contains one control point more than the input. + + const naxis = 2 + // nsegs is the number of segments, which is the number of + // variable control points. + nsegs := len(pts) - 1 + + // Variables begin at offset 1; offset 0 is for the + // constants. + const varOff = 1 + // nctrl is the number of control point variables, one per segment + // and axis, one λ per segment, and the goal variable J. + nctrl := nsegs*(naxis+1) + 1 + ctrl := func(axis, i int) expr { + idx := varOff + axis*nsegs + i + return varExpr(idx) + } + λ := func(i int) expr { + λOff := varOff + nsegs*naxis + return varExpr(λOff + i) + } + // The goal variable is last. + jOff := varOff + nctrl - 1 + J := varExpr(jOff) + // Separating equality constraints seems to help LP + // robustness. + var eqs, ineqs []expr + // addConstraint adds a constraint on the form + // + // expr op 0 where op is '=' or '≤'. + // + // It returns an expr for the slack variable, if any. + addConstraint := func(cons expr, op rune) expr { + if cons.IsZero() { + return expr{} + } + if op == '≤' { + slackIdx := jOff + 1 + len(ineqs) + slack := varExpr(slackIdx) + cons = cons.Add(slack) + ineqs = append(ineqs, cons) + return slack + } + eqs = append(eqs, cons) + return expr{} + } + // Add an inequality constraint on the form + // + // |e| <= J + // + // It returns an equivalent expression, using + // the goal and a slack variable. + addKinematic := func(e expr, scale float64) expr { + if e.IsZero() { + return expr{} + } + // The absolute function expands to an inequality constraint + // for each direction, on the standard form + // + // ±e - J <= 0 + // + s := addConstraint(e.Mul(+scale).Sub(J), '≤') + addConstraint(e.Mul(-scale).Sub(J), '≤') + return J.Sub(s).Div(scale) + } + + // Acceleration and jerk scale non-linearly (squared and cubed) with + // time scaling. Linear programs don't support non-linear scaling, so + // instead scale the contributions to the goal. + const ( + vScale = 1 + aScale = 1 + jScale = 10 + ) + derive := func(knots [4]uint, p0, p1 expr, degree int) expr { + t := uint(0) + for _, k := range knots[1 : degree+1] { + t += k + } + res := expr{} + if t != 0 { + res = p1.Sub(p0).Mul(float64(degree) / float64(t)) + } + return res + } + type state struct { + knots [4]uint + λ [2]expr + s [2]float64 + p [3]expr + v, a expr + } + var vExprs, aExprs, jExprs []expr + nknots := 0 + knot := func(last state, t uint, p expr, s float64, λ expr) state { + copy(last.knots[:], last.knots[1:]) + last.knots[3] = t + // Construct geometric constraints such that each knot of the B-spline + // hit a line segment. The start and end control points do not need + // constraints because they're fixed. + // + // A uniform B-spline at knot t{i} evaluates to + // + // BS(t{i}) = B{i-1}*C{i-1} + B{i}*C{i} + B{i+1}*C{i+1} + // + // Where B{i} is the ith basis function of third degree. + // + // Knot t{i} must lie in the line segment between S{i} and S{i+1}: + // + // BS(t{i}) = (1-λ{i})*S{i} + λ{i}*S{i+1}, λ{i}∈]0;1[ + // = S{i} + (S{i+1}-S{i})*λ{i} + // + // The λ intervals are open to avoid overlapping control points. + // Overlapping control points reduce the continuity of B-splines. + // + // The combination leads to the constraint, + // + // B{i-1}*C{i-1} + B{i}*C{i} + B{i+1}*C{i+1} = S{i} + (S{i+1}-S{i})*λ{i} + // + // On standard form: + // + // B{i-1}*C{i-1} + B{i}*C{i} + B{i+1}*C{i+1} - (S{i+1}-S{i})*λ{i} - S{i} = 0 + B := bsplineBasis(last.knots) + s0, s1 := last.s[0], last.s[1] + ctrl := expr{} + for j, b := range B { + ctrl = ctrl.Add(last.p[j].Mul(b)) + } + cons := ctrl.Add(last.λ[0].Mul(-(s1 - s0))).Sub(constExpr(s0)) + addConstraint(cons, '=') + + // Shift control point and waypoint into state. + copy(last.s[:], last.s[1:]) + last.s[1] = s + copy(last.p[:], last.p[1:]) + last.p[2] = p + copy(last.λ[:], last.λ[1:]) + last.λ[1] = λ + v := derive(last.knots, last.p[0], last.p[1], 3) + a := derive(last.knots, last.v, v, 2) + last.v = v + j := derive(last.knots, last.a, a, 1) + last.a = a + + // Add kinetic constraints. + if nknots >= 3 { + vExprs = append(vExprs, addKinematic(v, vScale)) + if nknots >= 4 { + aExprs = append(aExprs, addKinematic(a, aScale)) + if nknots >= 5 { + jExprs = append(jExprs, addKinematic(j, jScale)) + } + } + } + nknots++ + return last + } + // Improve the linear program conditioning by normalizing the input points + // to the [1;2] range. The [0;1] range is left for control points because + // the standard linear program model forces control coordinates to be + // non-negative. The alternative is to split control coordinates in two + // variables per axis. + var minPt, ptScale [naxis]float64 + const ptOffset = 1 + // For each axis, construct symbolic control points and apply + // constraints on them. + for axis := range naxis { + nknots = 0 + mi, ma := math.Inf(+1), math.Inf(-1) + for _, p := range pts { + pf := [naxis]float64{float64(p.X), float64(p.Y)}[axis] + mi, ma = min(mi, pf), max(ma, pf) + } + scale := max(1, ma-mi) + s := func(i int) float64 { + p := pts[i] + pf := [naxis]float64{float64(p.X), float64(p.Y)}[axis] + return (pf-mi)/scale + ptOffset + } + minPt[axis] = mi + ptScale[axis] = scale + // state is a sliding window of previous control points and + // kinetic values. + var last state + // Fixed start control points. + start := s(0) + for range 3 { + last = knot(last, 0, constExpr(start), start, expr{}) + } + // Variable inner control points. + for i := range nsegs { + last = knot(last, 1, ctrl(axis, i), s(i), λ(i)) + } + // Fixed end points. + end := s(nsegs) + for i := range 3 { + t := uint(1) + if i > 0 { + t = 0 + } + last = knot(last, t, constExpr(end), end, expr{}) + } + } + + // Constrain the weights by expressing w{i}∈]0,1[ in + // terms of a small constant: + // + // 0 <= w{i} <= 1-ε + // + // Since the first constraint is implied in an LP, there + // is only the upper limit: + // + // w{i} <= 1-ε + // + const ε = .05 + for i := range nsegs { + cons := λ(i).Add(constExpr(-(1 - ε))) + addConstraint(cons, '≤') + } + // As an edge case, the first weight must not be + // zero, otherwise its point may overlap the previous, + // fixed point. In standard form: + // + // -w{0} <= -ε + cons := λ(0).Mul(-1).Add(constExpr(ε)) + addConstraint(cons, '≤') + // Total number of constraints. + ncons := len(ineqs) + len(eqs) + // Total number of variables. + nvars := nctrl + len(ineqs) + A := mat.NewDense(ncons, nvars, nil) + b := make([]float64, ncons) + // Copy constraints to A. + for row, c := range ineqs { + cons := c.Explode() + for j, v := range cons[1:] { + A.Set(row, j, v) + } + b[row] = -cons[0] + } + for i, c := range eqs { + row := len(ineqs) + i + cons := c.Explode() + for j, v := range cons[1:] { + A.Set(row, j, v) + } + b[row] = -cons[0] + } + // Minimize the goal, J. + c := make([]float64, nvars) + copy(c, J.Explode()[1:]) + _, x, err := lp.Simplex(c, A, b, 1e-6, nil) + if err != nil { + return nil, nil, nil, nil, err + } + // eval evaluates an expression on the LP result. + eval := func(e expr) float64 { + coeffs := e.Explode() + v := coeffs[0] + for i, coeff := range coeffs[1:] { + v += coeff * x[i] + } + return v + } + // Recover internal control points. + toPoint := func(v [2]float64) bezier.Point { + return bezier.Point{ + X: int(math.Round(v[0])), + Y: int(math.Round(v[1])), + } + } + ctrls := make([]bezier.Point, 0, nsegs+2) + // Clamping start point. + s, e := pts[0], pts[len(pts)-1] + ctrls = append(ctrls, s, s, s) + for i := range nsegs { + var c [naxis]float64 + for axis := range naxis { + // Undo robustness transformation. + v := eval(ctrl(axis, i)) + scale, mi := ptScale[axis], minPt[axis] + c[axis] = (v-ptOffset)*scale + mi + } + ctrls = append(ctrls, toPoint(c)) + } + // Clamping end points. + ctrls = append(ctrls, e, e, e) + + // Recover kinematic values. + recoverKin := func(exprs []expr) []bezier.Point { + var v []bezier.Point + nvals := len(exprs) / 2 + for i := range nvals { + var c [naxis]float64 + for axis := range c { + scale := ptScale[axis] + e := exprs[axis*nvals+i] + c[axis] = eval(e) * scale + } + v = append(v, toPoint(c)) + } + return v + } + v = recoverKin(vExprs) + a = recoverKin(aExprs) + j = recoverKin(jExprs) + return ctrls, v, a, j, nil +} + +// bsplineBasis computes the coefficients of the B-spline control +// points at the start of the segment. +func bsplineBasis(knots [4]uint) [3]float64 { + dt1, dt2, dt3, dt4 := float64(knots[0]), float64(knots[1]), float64(knots[2]), float64(knots[3]) + if dt3 == 0 { + // Zero duration curve. + return [...]float64{0, 1, 0} + } + d1 := dt4 + dt3 + dt2 + p334c2, p334c3 := (dt4+dt3)/d1, dt2/d1 + d2 := dt3 + dt2 + dt1 + p323c1, p323c2 := dt3/d2, (dt2+dt1)/d2 + d4 := dt3 + dt2 + c1, c2, c3 := p323c1*dt3/d4, (p323c2*dt3+p334c2*dt2)/d4, p334c3*dt2/d4 + return [...]float64{c1, c2, c3} +} + +// constExpr returns an expression with the constant coefficient +// set to c. +func constExpr(c float64) expr { + return expr{ + c0: c, + } +} + +// varExpr returns an expression with the ith coefficient set +// to 1. +func varExpr(i int) expr { + if i == 0 { + return constExpr(1) + } + s := expr{ + zeros: i - 1, + c: make([]float64, 1), + } + s.c[0] = 1 + return s.normalize() +} + +// expr is a value represented by a vector of coefficients +// c{i} for implicit variables v{i}: +// +// e = c{0} + c{1}v{0}...c{n}v{n-1} +type expr struct { + // c0 stores c{0}. + c0 float64 + // c stores the extra coefficients. The layout is + // + // c[i] = c{zeros+i} + c []float64 + // zeros store the number of implied zero + // coefficients between c0 and c[0]. + zeros int +} + +func (s expr) IsZero() bool { + return s.c0 == 0 && len(s.c) == 0 +} + +func (s expr) Explode() []float64 { + r := make([]float64, s.numCoeffs()) + r[0] = s.c0 + copy(r[1+s.zeros:], s.c) + return r +} + +func (s expr) numCoeffs() int { + return 1 + s.zeros + len(s.c) +} + +func (s expr) String() string { + coeffs := s.Explode() + if len(coeffs) == 1 { + return fmt.Sprintf("%g", coeffs[0]) + } + return fmt.Sprintf("%g", coeffs) +} + +func (s expr) Const() float64 { + if len(s.c) > 0 { + panic("non-const expression") + } + return s.c0 +} + +func (s expr) copy() expr { + c := expr{ + zeros: s.zeros, + c0: s.c0, + c: make([]float64, len(s.c)), + } + copy(c.c, s.c) + return c +} + +func (s expr) Mul(f float64) expr { + sf := s.copy() + sf.c0 *= f + for i := range sf.c { + sf.c[i] *= f + } + return sf.normalize() +} + +func (s expr) Div(f float64) expr { + sf := s.copy() + sf.c0 /= f + for i := range sf.c { + sf.c[i] /= f + } + return sf.normalize() +} + +func (s expr) Sub(s2 expr) expr { + return s.Add(s2.Mul(-1)) +} + +func (s expr) Add(s2 expr) expr { + cmin := min(s.zeros, s2.zeros) + cmax := max(s.zeros+len(s.c), s2.zeros+len(s2.c)) + r := expr{ + c0: s.c0 + s2.c0, + zeros: cmin, + c: make([]float64, cmax-cmin), + } + copy(r.c[s.zeros-cmin:], s.c) + for i, c := range s2.c { + r.c[s2.zeros-cmin+i] += c + } + return r.normalize() +} + +// normalize left-adjusts the coefficients. +func (s expr) normalize() expr { + for len(s.c) > 0 { + n := len(s.c) + if s.c[n-1] != 0 { + break + } + s.c = s.c[:n-1] + } + for i, c := range s.c { + if c == 0 { + continue + } + copy(s.c, s.c[i:]) + s.c = s.c[:len(s.c)-i] + s.zeros += i + break + } + return s +} diff --git a/bspline/testdata/curves.png b/bspline/testdata/curves.png new file mode 100644 index 0000000000000000000000000000000000000000..3abb5c0bed2c6a8f6a87344c2a93e3ba1ed319fc GIT binary patch literal 655 zcmeAS@N?(olHy`uVBq!ia0vp^(?FPm1x)VDWME(l^K@|xskrra?%ATl3OsCo|LObc z+nQvmP3xIpckta-V+k?N&&iiw-(vqOz_lRPpqT zrXy3mMK_)Kd`u+t@$`GEf+iRISs2eXtv5V&%HC@mckFw0sU>8kclFx-nCV%blY6v% zGanT$Re2j5I{SP1zD?hko6b>8^ZdPj)mx?Xy&Jw?isZKnft%W+=@GLNl`9yo9U>Uc z7n&fqmLXJh)iL*guYngP@0VYAN@Q)5_L-|JTRHp1uWKDQW^ib&Jojs^ z)D`ZzslwMo zGlk!t^4-3o&nx!LIorHE?H#wT9SF%e)0MaU)k_B1z^l(wr@Y%y{q4}>Z)fLERW6*C zywat_Rr~&8$3u4YB3qWZ@0}cZka_){_1rZt>*x0QGOg-a*`cvnFGR9c;nW)LP`B6) zqt{MXuPl);xn4WhW#MVlJ&RA9<~mHzcIwibE>yBJWx?sCjInyG?60egUfP6QKQmEh zQOEjFUfV}OX-~wBrC)d5cfJvJeV4X+hn4Xz6FoJBEb-W%(uZr#FB76!q=zd)uDCblt6tz$AnVWc^@#R{YIG SPk-Y#kP=T