~whereswaldon/gio-x

120523b876108d515e3a6f68ee99eee58d52c934 — Elias Naur 2 years ago 199c05a strokes
stroke: add package for approximating complex strokes

This is a port of Sebastien Binet's stroke approximation algorithms.
They are no longer in Gio core because it is not clear they can be
implemented on a GPU.

Signed-off-by: Elias Naur <mail@eliasnaur.com>
A stroke/approx.go => stroke/approx.go +858 -0
@@ 0,0 1,858 @@
// SPDX-License-Identifier: Unlicense OR MIT

// Most of the algorithms to compute strokes and their offsets have been
// extracted, adapted from (and used as a reference implementation):
//  - github.com/tdewolff/canvas (Licensed under MIT)
//
// These algorithms have been implemented from:
//  Fast, precise flattening of cubic Bézier path and offset curves
//   Thomas F. Hain, et al.
//
// An electronic version is available at:
//  https://seant23.files.wordpress.com/2010/11/fastpreciseflatteningofbeziercurve.pdf
//
// Possible improvements (in term of speed and/or accuracy) on these
// algorithms are:
//
//  - Polar Stroking: New Theory and Methods for Stroking Paths,
//    M. Kilgard
//    https://arxiv.org/pdf/2007.00308.pdf
//
//  - https://raphlinus.github.io/graphics/curves/2019/12/23/flatten-quadbez.html
//    R. Levien

package stroke

import (
	"math"

	"gioui.org/f32"
)

// strokeTolerance is used to reconcile rounding errors arising
// when splitting quads into smaller and smaller segments to approximate
// them into straight lines, and when joining back segments.
//
// The magic value of 0.01 was found by striking a compromise between
// aesthetic looking (curves did look like curves, even after linearization)
// and speed.
const strokeTolerance = 0.01

type quadSegment struct {
	From, Ctrl, To f32.Point
}

type strokeQuad struct {
	Contour int
	Quad    quadSegment
}

type strokeState struct {
	p0, p1 f32.Point // p0 is the start point, p1 the end point.
	n0, n1 f32.Point // n0 is the normal vector at the start point, n1 at the end point.
	r0, r1 float32   // r0 is the curvature at the start point, r1 at the end point.
	ctl    f32.Point // ctl is the control point of the quadratic Bézier segment.
}

type strokeQuads []strokeQuad

func (qs *strokeQuads) setContour(n int) {
	for i := range *qs {
		(*qs)[i].Contour = n
	}
}

func (qs *strokeQuads) pen() f32.Point {
	return (*qs)[len(*qs)-1].Quad.To
}

func (qs *strokeQuads) closed() bool {
	beg := (*qs)[0].Quad.From
	end := (*qs)[len(*qs)-1].Quad.To
	return f32Eq(beg.X, end.X) && f32Eq(beg.Y, end.Y)
}

func (qs *strokeQuads) lineTo(pt f32.Point) {
	end := qs.pen()
	*qs = append(*qs, strokeQuad{
		Quad: quadSegment{
			From: end,
			Ctrl: end.Add(pt).Mul(0.5),
			To:   pt,
		},
	})
}

func (qs *strokeQuads) arc(f1, f2 f32.Point, angle float32) {
	const segments = 16
	pen := qs.pen()
	m := arcTransform(pen, f1.Add(pen), f2.Add(pen), angle, segments)
	for i := 0; i < segments; i++ {
		p0 := qs.pen()
		p1 := m.Transform(p0)
		p2 := m.Transform(p1)
		ctl := p1.Mul(2).Sub(p0.Add(p2).Mul(.5))
		*qs = append(*qs, strokeQuad{
			Quad: quadSegment{
				From: p0, Ctrl: ctl, To: p2,
			},
		})
	}
}

// split splits a slice of quads into slices of quads grouped
// by contours (ie: splitted at move-to boundaries).
func (qs strokeQuads) split() []strokeQuads {
	if len(qs) == 0 {
		return nil
	}

	var (
		c int
		o []strokeQuads
		i = len(o)
	)
	for _, q := range qs {
		if q.Contour != c {
			c = q.Contour
			i = len(o)
			o = append(o, strokeQuads{})
		}
		o[i] = append(o[i], q)
	}

	return o
}

func (qs strokeQuads) stroke(stroke Stroke) strokeQuads {
	if !isSolidLine(stroke.Dashes) {
		qs = qs.dash(stroke.Dashes)
	}

	var (
		o  strokeQuads
		hw = 0.5 * stroke.Width
	)

	for _, ps := range qs.split() {
		rhs, lhs := ps.offset(hw, stroke)
		switch lhs {
		case nil:
			o = o.append(rhs)
		default:
			// Closed path.
			// Inner path should go opposite direction to cancel outer path.
			switch {
			case ps.ccw():
				lhs = lhs.reverse()
				o = o.append(rhs)
				o = o.append(lhs)
			default:
				rhs = rhs.reverse()
				o = o.append(lhs)
				o = o.append(rhs)
			}
		}
	}

	return o
}

// offset returns the right-hand and left-hand sides of the path, offset by
// the half-width hw.
// The stroke handles how segments are joined and ends are capped.
func (qs strokeQuads) offset(hw float32, stroke Stroke) (rhs, lhs strokeQuads) {
	var (
		states []strokeState
		beg    = qs[0].Quad.From
		end    = qs[len(qs)-1].Quad.To
		closed = beg == end
	)
	for i := range qs {
		q := qs[i].Quad

		var (
			n0 = strokePathNorm(q.From, q.Ctrl, q.To, 0, hw)
			n1 = strokePathNorm(q.From, q.Ctrl, q.To, 1, hw)
			r0 = strokePathCurv(q.From, q.Ctrl, q.To, 0)
			r1 = strokePathCurv(q.From, q.Ctrl, q.To, 1)
		)
		states = append(states, strokeState{
			p0:  q.From,
			p1:  q.To,
			n0:  n0,
			n1:  n1,
			r0:  r0,
			r1:  r1,
			ctl: q.Ctrl,
		})
	}

	for i, state := range states {
		rhs = rhs.append(strokeQuadBezier(state, +hw, strokeTolerance))
		lhs = lhs.append(strokeQuadBezier(state, -hw, strokeTolerance))

		// join the current and next segments
		if hasNext := i+1 < len(states); hasNext || closed {
			var next strokeState
			switch {
			case hasNext:
				next = states[i+1]
			case closed:
				next = states[0]
			}
			if state.n1 != next.n0 {
				strokePathJoin(stroke, &rhs, &lhs, hw, state.p1, state.n1, next.n0, state.r1, next.r0)
			}
		}
	}

	if closed {
		rhs.close()
		lhs.close()
		return rhs, lhs
	}

	qbeg := &states[0]
	qend := &states[len(states)-1]

	// Default to counter-clockwise direction.
	lhs = lhs.reverse()
	strokePathCap(stroke, &rhs, hw, qend.p1, qend.n1)

	rhs = rhs.append(lhs)
	strokePathCap(stroke, &rhs, hw, qbeg.p0, qbeg.n0.Mul(-1))

	rhs.close()

	return rhs, nil
}

func (qs *strokeQuads) close() {
	p0 := (*qs)[len(*qs)-1].Quad.To
	p1 := (*qs)[0].Quad.From

	if p1 == p0 {
		return
	}

	*qs = append(*qs, strokeQuad{
		Quad: quadSegment{
			From: p0,
			Ctrl: p0.Add(p1).Mul(0.5),
			To:   p1,
		},
	})
}

// ccw returns whether the path is counter-clockwise.
func (qs strokeQuads) ccw() bool {
	// Use the Shoelace formula:
	//  https://en.wikipedia.org/wiki/Shoelace_formula
	var area float32
	for _, ps := range qs.split() {
		for i := 1; i < len(ps); i++ {
			pi := ps[i].Quad.To
			pj := ps[i-1].Quad.To
			area += (pi.X - pj.X) * (pi.Y + pj.Y)
		}
	}
	return area <= 0.0
}

func (qs strokeQuads) reverse() strokeQuads {
	if len(qs) == 0 {
		return nil
	}

	ps := make(strokeQuads, 0, len(qs))
	for i := range qs {
		q := qs[len(qs)-1-i]
		q.Quad.To, q.Quad.From = q.Quad.From, q.Quad.To
		ps = append(ps, q)
	}

	return ps
}

func (qs strokeQuads) append(ps strokeQuads) strokeQuads {
	switch {
	case len(ps) == 0:
		return qs
	case len(qs) == 0:
		return ps
	}

	// Consolidate quads and smooth out rounding errors.
	// We need to also check for the strokeTolerance to correctly handle
	// join/cap points or on-purpose disjoint quads.
	p0 := qs[len(qs)-1].Quad.To
	p1 := ps[0].Quad.From
	if p0 != p1 && lenPt(p0.Sub(p1)) < strokeTolerance {
		qs = append(qs, strokeQuad{
			Quad: quadSegment{
				From: p0,
				Ctrl: p0.Add(p1).Mul(0.5),
				To:   p1,
			},
		})
	}
	return append(qs, ps...)
}

func (q quadSegment) Transform(t f32.Affine2D) quadSegment {
	q.From = t.Transform(q.From)
	q.Ctrl = t.Transform(q.Ctrl)
	q.To = t.Transform(q.To)
	return q
}

// strokePathNorm returns the normal vector at t.
func strokePathNorm(p0, p1, p2 f32.Point, t, d float32) f32.Point {
	switch t {
	case 0:
		n := p1.Sub(p0)
		if n.X == 0 && n.Y == 0 {
			return f32.Point{}
		}
		n = rot90CW(n)
		return normPt(n, d)
	case 1:
		n := p2.Sub(p1)
		if n.X == 0 && n.Y == 0 {
			return f32.Point{}
		}
		n = rot90CW(n)
		return normPt(n, d)
	}
	panic("impossible")
}

func rot90CW(p f32.Point) f32.Point  { return f32.Pt(+p.Y, -p.X) }
func rot90CCW(p f32.Point) f32.Point { return f32.Pt(-p.Y, +p.X) }

// cosPt returns the cosine of the opening angle between p and q.
func cosPt(p, q f32.Point) float32 {
	np := math.Hypot(float64(p.X), float64(p.Y))
	nq := math.Hypot(float64(q.X), float64(q.Y))
	return dotPt(p, q) / float32(np*nq)
}

func normPt(p f32.Point, l float32) f32.Point {
	d := math.Hypot(float64(p.X), float64(p.Y))
	l64 := float64(l)
	if math.Abs(d-l64) < 1e-10 {
		return f32.Point{}
	}
	n := float32(l64 / d)
	return f32.Point{X: p.X * n, Y: p.Y * n}
}

func lenPt(p f32.Point) float32 {
	return float32(math.Hypot(float64(p.X), float64(p.Y)))
}

func dotPt(p, q f32.Point) float32 {
	return p.X*q.X + p.Y*q.Y
}

func perpDot(p, q f32.Point) float32 {
	return p.X*q.Y - p.Y*q.X
}

// strokePathCurv returns the curvature at t, along the quadratic Bézier
// curve defined by the triplet (beg, ctl, end).
func strokePathCurv(beg, ctl, end f32.Point, t float32) float32 {
	var (
		d1p = quadBezierD1(beg, ctl, end, t)
		d2p = quadBezierD2(beg, ctl, end, t)

		// Negative when bending right, ie: the curve is CW at this point.
		a = float64(perpDot(d1p, d2p))
	)

	// We check early that the segment isn't too line-like and
	// save a costly call to math.Pow that will be discarded by dividing
	// with a too small 'a'.
	if math.Abs(a) < 1e-10 {
		return float32(math.NaN())
	}
	return float32(math.Pow(float64(d1p.X*d1p.X+d1p.Y*d1p.Y), 1.5) / a)
}

// quadBezierSample returns the point on the Bézier curve at t.
//  B(t) = (1-t)^2 P0 + 2(1-t)t P1 + t^2 P2
func quadBezierSample(p0, p1, p2 f32.Point, t float32) f32.Point {
	t1 := 1 - t
	c0 := t1 * t1
	c1 := 2 * t1 * t
	c2 := t * t

	o := p0.Mul(c0)
	o = o.Add(p1.Mul(c1))
	o = o.Add(p2.Mul(c2))
	return o
}

// quadBezierD1 returns the first derivative of the Bézier curve with respect to t.
//  B'(t) = 2(1-t)(P1 - P0) + 2t(P2 - P1)
func quadBezierD1(p0, p1, p2 f32.Point, t float32) f32.Point {
	p10 := p1.Sub(p0).Mul(2 * (1 - t))
	p21 := p2.Sub(p1).Mul(2 * t)

	return p10.Add(p21)
}

// quadBezierD2 returns the second derivative of the Bézier curve with respect to t:
//  B''(t) = 2(P2 - 2P1 + P0)
func quadBezierD2(p0, p1, p2 f32.Point, t float32) f32.Point {
	p := p2.Sub(p1.Mul(2)).Add(p0)
	return p.Mul(2)
}

// quadBezierLen returns the length of the Bézier curve.
// See:
//  https://malczak.linuxpl.com/blog/quadratic-bezier-curve-length/
func quadBezierLen(p0, p1, p2 f32.Point) float32 {
	a := p0.Sub(p1.Mul(2)).Add(p2)
	b := p1.Mul(2).Sub(p0.Mul(2))
	A := float64(4 * dotPt(a, a))
	B := float64(4 * dotPt(a, b))
	C := float64(dotPt(b, b))
	if f64Eq(A, 0.0) {
		// p1 is in the middle between p0 and p2,
		// so it is a straight line from p0 to p2.
		return lenPt(p2.Sub(p0))
	}

	Sabc := 2 * math.Sqrt(A+B+C)
	A2 := math.Sqrt(A)
	A32 := 2 * A * A2
	C2 := 2 * math.Sqrt(C)
	BA := B / A2
	return float32((A32*Sabc + A2*B*(Sabc-C2) + (4*C*A-B*B)*math.Log((2*A2+BA+Sabc)/(BA+C2))) / (4 * A32))
}

func strokeQuadBezier(state strokeState, d, flatness float32) strokeQuads {
	// Gio strokes are only quadratic Bézier curves, w/o any inflection point.
	// So we just have to flatten them.
	var qs strokeQuads
	return flattenQuadBezier(qs, state.p0, state.ctl, state.p1, d, flatness)
}

// flattenQuadBezier splits a Bézier quadratic curve into linear sub-segments,
// themselves also encoded as Bézier (degenerate, flat) quadratic curves.
func flattenQuadBezier(qs strokeQuads, p0, p1, p2 f32.Point, d, flatness float32) strokeQuads {
	var (
		t      float32
		flat64 = float64(flatness)
	)
	for t < 1 {
		s2 := float64((p2.X-p0.X)*(p1.Y-p0.Y) - (p2.Y-p0.Y)*(p1.X-p0.X))
		den := math.Hypot(float64(p1.X-p0.X), float64(p1.Y-p0.Y))
		if s2*den == 0.0 {
			break
		}

		s2 /= den
		t = 2.0 * float32(math.Sqrt(flat64/3.0/math.Abs(s2)))
		if t >= 1.0 {
			break
		}
		var q0, q1, q2 f32.Point
		q0, q1, q2, p0, p1, p2 = quadBezierSplit(p0, p1, p2, t)
		qs.addLine(q0, q1, q2, 0, d)
	}
	qs.addLine(p0, p1, p2, 1, d)
	return qs
}

func (qs *strokeQuads) addLine(p0, ctrl, p1 f32.Point, t, d float32) {

	switch i := len(*qs); i {
	case 0:
		p0 = p0.Add(strokePathNorm(p0, ctrl, p1, 0, d))
	default:
		// Address possible rounding errors and use previous point.
		p0 = (*qs)[i-1].Quad.To
	}

	p1 = p1.Add(strokePathNorm(p0, ctrl, p1, 1, d))

	*qs = append(*qs,
		strokeQuad{
			Quad: quadSegment{
				From: p0,
				Ctrl: p0.Add(p1).Mul(0.5),
				To:   p1,
			},
		},
	)
}

// quadInterp returns the interpolated point at t.
func quadInterp(p, q f32.Point, t float32) f32.Point {
	return f32.Pt(
		(1-t)*p.X+t*q.X,
		(1-t)*p.Y+t*q.Y,
	)
}

// quadBezierSplit returns the pair of triplets (from,ctrl,to) Bézier curve,
// split before (resp. after) the provided parametric t value.
func quadBezierSplit(p0, p1, p2 f32.Point, t float32) (f32.Point, f32.Point, f32.Point, f32.Point, f32.Point, f32.Point) {

	var (
		b0 = p0
		b1 = quadInterp(p0, p1, t)
		b2 = quadBezierSample(p0, p1, p2, t)

		a0 = b2
		a1 = quadInterp(p1, p2, t)
		a2 = p2
	)

	return b0, b1, b2, a0, a1, a2
}

// strokePathJoin joins the two paths rhs and lhs, according to the provided
// stroke operation.
func strokePathJoin(stroke Stroke, rhs, lhs *strokeQuads, hw float32, pivot, n0, n1 f32.Point, r0, r1 float32) {
	if stroke.Miter > 0 {
		strokePathMiterJoin(stroke, rhs, lhs, hw, pivot, n0, n1, r0, r1)
		return
	}
	switch stroke.Join {
	case BevelJoin:
		strokePathBevelJoin(rhs, lhs, hw, pivot, n0, n1, r0, r1)
	case RoundJoin:
		strokePathRoundJoin(rhs, lhs, hw, pivot, n0, n1, r0, r1)
	default:
		panic("impossible")
	}
}

func strokePathBevelJoin(rhs, lhs *strokeQuads, hw float32, pivot, n0, n1 f32.Point, r0, r1 float32) {

	rp := pivot.Add(n1)
	lp := pivot.Sub(n1)

	rhs.lineTo(rp)
	lhs.lineTo(lp)
}

func strokePathRoundJoin(rhs, lhs *strokeQuads, hw float32, pivot, n0, n1 f32.Point, r0, r1 float32) {
	rp := pivot.Add(n1)
	lp := pivot.Sub(n1)
	cw := dotPt(rot90CW(n0), n1) >= 0.0
	switch {
	case cw:
		// Path bends to the right, ie. CW (or 180 degree turn).
		c := pivot.Sub(lhs.pen())
		angle := -math.Acos(float64(cosPt(n0, n1)))
		lhs.arc(c, c, float32(angle))
		lhs.lineTo(lp) // Add a line to accommodate for rounding errors.
		rhs.lineTo(rp)
	default:
		// Path bends to the left, ie. CCW.
		angle := math.Acos(float64(cosPt(n0, n1)))
		c := pivot.Sub(rhs.pen())
		rhs.arc(c, c, float32(angle))
		rhs.lineTo(rp) // Add a line to accommodate for rounding errors.
		lhs.lineTo(lp)
	}
}

func strokePathMiterJoin(stroke Stroke, rhs, lhs *strokeQuads, hw float32, pivot, n0, n1 f32.Point, r0, r1 float32) {
	if n0 == n1.Mul(-1) {
		strokePathBevelJoin(rhs, lhs, hw, pivot, n0, n1, r0, r1)
		return
	}

	// This is to handle nearly linear joints that would be clipped otherwise.
	limit := math.Max(float64(stroke.Miter), 1.001)

	cw := dotPt(rot90CW(n0), n1) >= 0.0
	if cw {
		// hw is used to calculate |R|.
		// When running CW, n0 and n1 point the other way,
		// so the sign of r0 and r1 is negated.
		hw = -hw
	}
	hw64 := float64(hw)

	cos := math.Sqrt(0.5 * (1 + float64(cosPt(n0, n1))))
	d := hw64 / cos
	if math.Abs(limit*hw64) < math.Abs(d) {
		stroke.Miter = 0 // Set miter to zero to disable the miter joint.
		strokePathJoin(stroke, rhs, lhs, hw, pivot, n0, n1, r0, r1)
		return
	}
	mid := pivot.Add(normPt(n0.Add(n1), float32(d)))

	rp := pivot.Add(n1)
	lp := pivot.Sub(n1)
	switch {
	case cw:
		// Path bends to the right, ie. CW.
		lhs.lineTo(mid)
	default:
		// Path bends to the left, ie. CCW.
		rhs.lineTo(mid)
	}
	rhs.lineTo(rp)
	lhs.lineTo(lp)
}

// strokePathCap caps the provided path qs, according to the provided stroke operation.
func strokePathCap(stroke Stroke, qs *strokeQuads, hw float32, pivot, n0 f32.Point) {
	switch stroke.Cap {
	case FlatCap:
		strokePathFlatCap(qs, hw, pivot, n0)
	case SquareCap:
		strokePathSquareCap(qs, hw, pivot, n0)
	case RoundCap:
		strokePathRoundCap(qs, hw, pivot, n0)
	default:
		panic("impossible")
	}
}

// strokePathFlatCap caps the start or end of a path with a flat cap.
func strokePathFlatCap(qs *strokeQuads, hw float32, pivot, n0 f32.Point) {
	end := pivot.Sub(n0)
	qs.lineTo(end)
}

// strokePathSquareCap caps the start or end of a path with a square cap.
func strokePathSquareCap(qs *strokeQuads, hw float32, pivot, n0 f32.Point) {
	var (
		e       = pivot.Add(rot90CCW(n0))
		corner1 = e.Add(n0)
		corner2 = e.Sub(n0)
		end     = pivot.Sub(n0)
	)

	qs.lineTo(corner1)
	qs.lineTo(corner2)
	qs.lineTo(end)
}

// strokePathRoundCap caps the start or end of a path with a round cap.
func strokePathRoundCap(qs *strokeQuads, hw float32, pivot, n0 f32.Point) {
	c := pivot.Sub(qs.pen())
	qs.arc(c, c, math.Pi)
}

// arcTransform computes a transformation that can be used for generating quadratic bézier
// curve approximations for an arc.
//
// The math is extracted from the following paper:
//  "Drawing an elliptical arc using polylines, quadratic or
//   cubic Bezier curves", L. Maisonobe
// An electronic version may be found at:
//  http://spaceroots.org/documents/ellipse/elliptical-arc.pdf
func arcTransform(p, f1, f2 f32.Point, angle float32, segments int) f32.Affine2D {
	c := f32.Point{
		X: 0.5 * (f1.X + f2.X),
		Y: 0.5 * (f1.Y + f2.Y),
	}

	// semi-major axis: 2a = |PF1| + |PF2|
	a := 0.5 * (dist(f1, p) + dist(f2, p))

	// semi-minor axis: c^2 = a^2+b^2 (c: focal distance)
	f := dist(f1, c)
	b := math.Sqrt(a*a - f*f)

	var rx, ry, alpha float64
	switch {
	case a > b:
		rx = a
		ry = b
	default:
		rx = b
		ry = a
	}

	var x float64
	switch {
	case f1 == c || f2 == c:
		// degenerate case of a circle.
		alpha = 0
	default:
		switch {
		case f1.X > c.X:
			x = float64(f1.X - c.X)
			alpha = math.Acos(x / f)
		case f1.X < c.X:
			x = float64(f2.X - c.X)
			alpha = math.Acos(x / f)
		case f1.X == c.X:
			// special case of a "vertical" ellipse.
			alpha = math.Pi / 2
			if f1.Y < c.Y {
				alpha = -alpha
			}
		}
	}

	var (
		θ   = angle / float32(segments)
		ref f32.Affine2D // transform from absolute frame to ellipse-based one
		rot f32.Affine2D // rotation matrix for each segment
		inv f32.Affine2D // transform from ellipse-based frame to absolute one
	)
	ref = ref.Offset(f32.Point{}.Sub(c))
	ref = ref.Rotate(f32.Point{}, float32(-alpha))
	ref = ref.Scale(f32.Point{}, f32.Point{
		X: float32(1 / rx),
		Y: float32(1 / ry),
	})
	inv = ref.Invert()
	rot = rot.Rotate(f32.Point{}, float32(0.5*θ))

	// Instead of invoking math.Sincos for every segment, compute a rotation
	// matrix once and apply for each segment.
	// Before applying the rotation matrix rot, transform the coordinates
	// to a frame centered to the ellipse (and warped into a unit circle), then rotate.
	// Finally, transform back into the original frame.
	return inv.Mul(rot).Mul(ref)
}

func dist(p1, p2 f32.Point) float64 {
	var (
		x1 = float64(p1.X)
		y1 = float64(p1.Y)
		x2 = float64(p2.X)
		y2 = float64(p2.Y)
		dx = x2 - x1
		dy = y2 - y1
	)
	return math.Hypot(dx, dy)
}

func strokePathCommands(stroke Stroke) strokeQuads {
	quads := convertToStrokeQuads(stroke.Path.Segments)
	return quads.stroke(stroke)
}

// convertToStrokeQuads converts segments to quads ready to stroke.
func convertToStrokeQuads(segs []Segment) strokeQuads {
	quads := make(strokeQuads, 0, len(segs))
	var pen f32.Point
	contour := 0
	for _, seg := range segs {
		switch seg.op {
		case segOpMoveTo:
			pen = seg.args[0]
			contour++
		case segOpLineTo:
			var q quadSegment
			q.From, q.To = pen, seg.args[0]
			q.Ctrl = q.From.Add(q.To).Mul(.5)
			quad := strokeQuad{
				Contour: contour,
				Quad:    q,
			}
			quads = append(quads, quad)
			pen = q.To
		case segOpQuadTo:
			var q quadSegment
			q.From, q.Ctrl, q.To = pen, seg.args[0], seg.args[1]
			quad := strokeQuad{
				Contour: contour,
				Quad:    q,
			}
			quads = append(quads, quad)
			pen = q.To
		case segOpCubeTo:
			for _, q := range splitCubic(pen, seg.args[0], seg.args[1], seg.args[2]) {
				quad := strokeQuad{
					Contour: contour,
					Quad:    q,
				}
				quads = append(quads, quad)
			}
			pen = seg.args[2]
		default:
			panic("unsupported segment op")
		}
	}
	return quads
}

func splitCubic(from, ctrl0, ctrl1, to f32.Point) []quadSegment {
	quads := make([]quadSegment, 0, 10)
	// Set the maximum distance proportionally to the longest side
	// of the bounding rectangle.
	hull := f32.Rectangle{
		Min: from,
		Max: ctrl0,
	}.Canon().Add(ctrl1).Add(to)
	l := hull.Dx()
	if h := hull.Dy(); h > l {
		l = h
	}
	approxCubeTo(&quads, 0, l*0.001, from, ctrl0, ctrl1, to)
	return quads
}

// approxCubeTo approximates a cubic Bézier by a series of quadratic
// curves.
func approxCubeTo(quads *[]quadSegment, splits int, maxDist float32, from, ctrl0, ctrl1, to f32.Point) int {
	// The idea is from
	// https://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
	// where a quadratic approximates a cubic by eliminating its t³ term
	// from its polynomial expression anchored at the starting point:
	//
	// P(t) = pen + 3t(ctrl0 - pen) + 3t²(ctrl1 - 2ctrl0 + pen) + t³(to - 3ctrl1 + 3ctrl0 - pen)
	//
	// The control point for the new quadratic Q1 that shares starting point, pen, with P is
	//
	// C1 = (3ctrl0 - pen)/2
	//
	// The reverse cubic anchored at the end point has the polynomial
	//
	// P'(t) = to + 3t(ctrl1 - to) + 3t²(ctrl0 - 2ctrl1 + to) + t³(pen - 3ctrl0 + 3ctrl1 - to)
	//
	// The corresponding quadratic Q2 that shares the end point, to, with P has control
	// point
	//
	// C2 = (3ctrl1 - to)/2
	//
	// The combined quadratic Bézier, Q, shares both start and end points with its cubic
	// and use the midpoint between the two curves Q1 and Q2 as control point:
	//
	// C = (3ctrl0 - pen + 3ctrl1 - to)/4
	c := ctrl0.Mul(3).Sub(from).Add(ctrl1.Mul(3)).Sub(to).Mul(1.0 / 4.0)
	const maxSplits = 32
	if splits >= maxSplits {
		*quads = append(*quads, quadSegment{From: from, Ctrl: c, To: to})
		return splits
	}
	// The maximum distance between the cubic P and its approximation Q given t
	// can be shown to be
	//
	// d = sqrt(3)/36*|to - 3ctrl1 + 3ctrl0 - pen|
	//
	// To save a square root, compare d² with the squared tolerance.
	v := to.Sub(ctrl1.Mul(3)).Add(ctrl0.Mul(3)).Sub(from)
	d2 := (v.X*v.X + v.Y*v.Y) * 3 / (36 * 36)
	if d2 <= maxDist*maxDist {
		*quads = append(*quads, quadSegment{From: from, Ctrl: c, To: to})
		return splits
	}
	// De Casteljau split the curve and approximate the halves.
	t := float32(0.5)
	c0 := from.Add(ctrl0.Sub(from).Mul(t))
	c1 := ctrl0.Add(ctrl1.Sub(ctrl0).Mul(t))
	c2 := ctrl1.Add(to.Sub(ctrl1).Mul(t))
	c01 := c0.Add(c1.Sub(c0).Mul(t))
	c12 := c1.Add(c2.Sub(c1).Mul(t))
	c0112 := c01.Add(c12.Sub(c01).Mul(t))
	splits++
	splits = approxCubeTo(quads, splits, maxDist, from, c0, c01, c0112)
	splits = approxCubeTo(quads, splits, maxDist, c0112, c12, c2, to)
	return splits
}

A stroke/clip_test.go => stroke/clip_test.go +385 -0
@@ 0,0 1,385 @@
// SPDX-License-Identifier: Unlicense OR MIT

package stroke

import (
	"math"
	"testing"

	"golang.org/x/image/colornames"

	"gioui.org/f32"
	"gioui.org/op"
	"gioui.org/op/clip"
	"gioui.org/op/paint"
)

func TestStrokedPathBevelFlat(t *testing.T) {
	run(t, func(o *op.Ops) {
		defer Stroke{
			Path:  fruit,
			Width: 2.5,
			Cap:   FlatCap,
			Join:  BevelJoin,
		}.Op(o).Push(o).Pop()

		paint.Fill(o, red)
	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(10, 50, colornames.Red)
	})
}

func TestStrokedPathBevelRound(t *testing.T) {
	run(t, func(o *op.Ops) {
		defer Stroke{
			Path:  fruit,
			Width: 2.5,
			Cap:   RoundCap,
			Join:  BevelJoin,
		}.Op(o).Push(o).Pop()

		paint.Fill(o, red)
	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(10, 50, colornames.Red)
	})
}

func TestStrokedPathBevelSquare(t *testing.T) {
	run(t, func(o *op.Ops) {
		defer Stroke{
			Path:  fruit,
			Width: 2.5,
			Cap:   SquareCap,
			Join:  BevelJoin,
		}.Op(o).Push(o).Pop()

		paint.Fill(o, red)
	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(10, 50, colornames.Red)
	})
}

func TestStrokedPathFlatMiter(t *testing.T) {
	run(t, func(o *op.Ops) {
		{
			cl := Stroke{
				Path:  zigzag,
				Width: 10,
				Cap:   FlatCap,
				Join:  BevelJoin,
				Miter: 5,
			}.Op(o).Push(o)
			paint.Fill(o, red)
			cl.Pop()
		}
		{
			cl := Stroke{
				Path:  zigzag,
				Width: 2,
				Cap:   FlatCap,
				Join:  BevelJoin,
			}.Op(o).Push(o)
			paint.Fill(o, black)
			cl.Pop()
		}

	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(40, 10, colornames.Black)
		r.expect(40, 12, colornames.Red)
	})
}

func TestStrokedPathFlatMiterInf(t *testing.T) {
	run(t, func(o *op.Ops) {
		{
			cl := Stroke{
				Path:  zigzag,
				Width: 10,
				Cap:   FlatCap,
				Join:  BevelJoin,
				Miter: float32(math.Inf(+1)),
			}.Op(o).Push(o)
			paint.Fill(o, red)
			cl.Pop()
		}
		{
			cl := Stroke{
				Path:  zigzag,
				Width: 2,
				Cap:   FlatCap,
				Join:  BevelJoin,
			}.Op(o).Push(o)
			paint.Fill(o, black)
			cl.Pop()
		}
	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(40, 10, colornames.Black)
		r.expect(40, 12, colornames.Red)
	})
}

func TestStrokedPathZeroWidth(t *testing.T) {
	run(t, func(o *op.Ops) {
		{
			p := new(clip.Path)
			p.Begin(o)
			p.Move(f32.Pt(10, 50))
			p.Line(f32.Pt(50, 0))
			cl := clip.Stroke{
				Path: p.End(),
				Style: clip.StrokeStyle{
					Width: 2,
					Cap:   clip.FlatCap,
					Join:  clip.BevelJoin,
				},
			}.Op().Push(o)

			paint.Fill(o, black)
			cl.Pop()
		}

		{
			p := new(clip.Path)
			p.Begin(o)
			p.Move(f32.Pt(10, 50))
			p.Line(f32.Pt(30, 0))
			cl := clip.Stroke{
				Path: p.End(),
			}.Op().Push(o) // width=0, disable stroke

			paint.Fill(o, red)
			cl.Pop()
		}

	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(10, 50, colornames.Black)
		r.expect(30, 50, colornames.Black)
		r.expect(65, 50, transparent)
	})
}

func TestDashedPathFlatCapEllipse(t *testing.T) {
	run(t, func(o *op.Ops) {
		{
			cl := Stroke{
				Path:  ellipse,
				Width: 10,
				Cap:   FlatCap,
				Join:  BevelJoin,
				Miter: float32(math.Inf(+1)),
				Dashes: Dashes{
					Dashes: []float32{5, 3},
				},
			}.Op(o).Push(o)

			paint.Fill(
				o,
				red,
			)
			cl.Pop()
		}
		{
			cl := Stroke{
				Path:  ellipse,
				Width: 2,
			}.Op(o).Push(o)

			paint.Fill(
				o,
				black,
			)
			cl.Pop()
		}

	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(0, 62, colornames.Red)
		r.expect(0, 65, colornames.Black)
	})
}

func TestDashedPathFlatCapZ(t *testing.T) {
	run(t, func(o *op.Ops) {
		{
			cl := Stroke{
				Path:  zigzag,
				Width: 10,
				Cap:   FlatCap,
				Join:  BevelJoin,
				Miter: float32(math.Inf(+1)),
				Dashes: Dashes{
					Dashes: []float32{5, 3},
				},
			}.Op(o).Push(o)
			paint.Fill(o, red)
			cl.Pop()
		}

		{
			cl := Stroke{
				Path:  zigzag,
				Width: 2,
				Cap:   FlatCap,
				Join:  BevelJoin,
			}.Op(o).Push(o)
			paint.Fill(o, black)
			cl.Pop()
		}
	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(40, 10, colornames.Black)
		r.expect(40, 12, colornames.Red)
		r.expect(46, 12, transparent)
	})
}

func TestDashedPathFlatCapZNoDash(t *testing.T) {
	run(t, func(o *op.Ops) {
		{
			cl := Stroke{
				Path:  zigzag,
				Width: 10,
				Cap:   FlatCap,
				Join:  BevelJoin,
				Miter: float32(math.Inf(+1)),
				Dashes: Dashes{
					Phase: 1,
				},
			}.Op(o).Push(o)
			paint.Fill(o, red)
			cl.Pop()
		}
		{
			cl := Stroke{
				Path:  zigzag,
				Width: 2,
				Cap:   FlatCap,
				Join:  BevelJoin,
			}.Op(o).Push(o)
			paint.Fill(o, black)
			cl.Pop()
		}
	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(40, 10, colornames.Black)
		r.expect(40, 12, colornames.Red)
		r.expect(46, 12, colornames.Red)
	})
}

func TestDashedPathFlatCapZNoPath(t *testing.T) {
	run(t, func(o *op.Ops) {
		{
			cl := Stroke{
				Path:  zigzag,
				Width: 10,
				Cap:   FlatCap,
				Join:  BevelJoin,
				Miter: float32(math.Inf(+1)),
				Dashes: Dashes{
					Dashes: []float32{0},
				},
			}.Op(o).Push(o)
			paint.Fill(o, red)
			cl.Pop()
		}
		{
			cl := Stroke{
				Path:  zigzag,
				Width: 2,
				Cap:   FlatCap,
				Join:  BevelJoin,
			}.Op(o).Push(o)
			paint.Fill(o, black)
			cl.Pop()
		}
	}, func(r result) {
		r.expect(0, 0, transparent)
		r.expect(40, 10, colornames.Black)
		r.expect(40, 12, transparent)
		r.expect(46, 12, transparent)
	})
}

var fruit = Path{
	Segments: []Segment{
		MoveTo(f32.Pt(10, 50)),
		LineTo(f32.Pt(20, 50)),
		QuadTo(f32.Point{X: 20.00035, Y: 48.607147}, f32.Point{X: 20.288229, Y: 47.240997}),
		QuadTo(f32.Point{X: 20.57679, Y: 45.874977}, f32.Point{X: 21.141825, Y: 44.588024}),
		QuadTo(f32.Point{X: 21.707504, Y: 43.301327}, f32.Point{X: 22.527983, Y: 42.143032}),
		QuadTo(f32.Point{X: 23.349041, Y: 40.985104}, f32.Point{X: 24.393435, Y: 39.99998}),
		QuadTo(f32.Point{X: 25.43832, Y: 39.01532}, f32.Point{X: 26.666492, Y: 38.241226}),
		QuadTo(f32.Point{X: 27.89505, Y: 37.467674}, f32.Point{X: 29.259802, Y: 36.934353}),
		QuadTo(f32.Point{X: 30.62482, Y: 36.401638}, f32.Point{X: 32.073708, Y: 36.12959}),
		QuadTo(f32.Point{X: 33.522728, Y: 35.858185}, f32.Point{X: 35.00007, Y: 35.857857}),
		QuadTo(f32.Point{X: 36.47741, Y: 35.858192}, f32.Point{X: 37.92643, Y: 36.129604}),
		QuadTo(f32.Point{X: 39.375313, Y: 36.401665}, f32.Point{X: 40.740334, Y: 36.93438}),
		QuadTo(f32.Point{X: 42.105087, Y: 37.467705}, f32.Point{X: 43.33364, Y: 38.241264}),
		QuadTo(f32.Point{X: 44.561806, Y: 39.01536}, f32.Point{X: 45.60668, Y: 40.00003}),
		QuadTo(f32.Point{X: 46.651073, Y: 40.985157}, f32.Point{X: 47.472122, Y: 42.14309}),
		QuadTo(f32.Point{X: 48.292587, Y: 43.301384}, f32.Point{X: 48.85826, Y: 44.58808}),
		QuadTo(f32.Point{X: 49.42329, Y: 45.87504}, f32.Point{X: 49.711845, Y: 47.241055}),
		QuadTo(f32.Point{X: 49.999718, Y: 48.60721}, f32.Point{X: 50.000053, Y: 50.000053}),
		LineTo(f32.Pt(60, 50)),
		LineTo(f32.Pt(70, 60)),
		QuadTo(f32.Point{X: 75.96515, Y: 60.01108}, f32.Point{X: 81.48046, Y: 62.283623}),
		QuadTo(f32.Point{X: 86.987305, Y: 64.57663}, f32.Point{X: 91.21312, Y: 68.78679}),
		QuadTo(f32.Point{X: 95.423294, Y: 73.01262}, f32.Point{X: 97.71627, Y: 78.519455}),
		QuadTo(f32.Point{X: 99.98879, Y: 84.034775}, f32.Point{X: 99.99987, Y: 89.999916}),
		QuadTo(f32.Point{X: 99.988785, Y: 95.96506}, f32.Point{X: 97.716255, Y: 101.48037}),
		QuadTo(f32.Point{X: 95.42325, Y: 106.9872}, f32.Point{X: 91.21309, Y: 111.21302}),
		QuadTo(f32.Point{X: 86.987274, Y: 115.42317}, f32.Point{X: 81.48043, Y: 117.71617}),
		QuadTo(f32.Point{X: 75.96512, Y: 119.9887}, f32.Point{X: 69.99997, Y: 119.99979}),
		QuadTo(f32.Point{X: 64.03482, Y: 119.9887}, f32.Point{X: 58.51951, Y: 117.71617}),
		QuadTo(f32.Point{X: 53.01267, Y: 115.42317}, f32.Point{X: 48.78685, Y: 111.213005}),
		QuadTo(f32.Point{X: 44.57669, Y: 106.98717}, f32.Point{X: 42.283707, Y: 101.48033}),
		QuadTo(f32.Point{X: 40.011185, Y: 95.96501}, f32.Point{X: 40.000122, Y: 89.99987}),
		QuadTo(f32.Point{X: 40.01121, Y: 84.03473}, f32.Point{X: 42.283745, Y: 78.519424}),
		QuadTo(f32.Point{X: 44.576748, Y: 73.01259}, f32.Point{X: 48.78691, Y: 68.78678}),
		QuadTo(f32.Point{X: 53.012737, Y: 64.57663}, f32.Point{X: 58.51957, Y: 62.283646}),
		QuadTo(f32.Point{X: 64.03488, Y: 60.01113}, f32.Point{X: 70.000015, Y: 60.00006}),
		LineTo(f32.Pt(50, 60)),
		QuadTo(f32.Pt(40, 50), f32.Pt(20, 90)),
	},
}

var zigzag = Path{
	Segments: []Segment{
		MoveTo(f32.Pt(40, 10)),
		LineTo(f32.Pt(90, 10)),
		LineTo(f32.Pt(40, 60)),
		LineTo(f32.Pt(90, 60)),
		QuadTo(f32.Pt(40, 80), f32.Pt(40, 110)),
		LineTo(f32.Pt(90, 110)),
	},
}

var ellipse = Path{
	Segments: []Segment{
		MoveTo(f32.Pt(0, 65)),
		LineTo(f32.Pt(20, 65)),
		QuadTo(f32.Point{X: 20.016611, Y: 57.560127}, f32.Point{X: 23.425419, Y: 50.681286}),
		QuadTo(f32.Point{X: 26.864927, Y: 43.81302}, f32.Point{X: 33.1802, Y: 38.542465}),
		QuadTo(f32.Point{X: 39.51897, Y: 33.291443}, f32.Point{X: 47.779266, Y: 30.431564}),
		QuadTo(f32.Point{X: 56.052277, Y: 27.59721}, f32.Point{X: 65.00003, Y: 27.583397}),
		QuadTo(f32.Point{X: 73.947784, Y: 27.59721}, f32.Point{X: 82.2208, Y: 30.431564}),
		QuadTo(f32.Point{X: 90.4811, Y: 33.291443}, f32.Point{X: 96.81986, Y: 38.542465}),
		QuadTo(f32.Point{X: 103.13513, Y: 43.813015}, f32.Point{X: 106.574646, Y: 50.681282}),
		QuadTo(f32.Point{X: 109.98345, Y: 57.56012}, f32.Point{X: 110.00008, Y: 64.99999}),
		QuadTo(f32.Point{X: 109.98346, Y: 72.439865}, f32.Point{X: 106.57466, Y: 79.3187}),
		QuadTo(f32.Point{X: 103.135155, Y: 86.18697}, f32.Point{X: 96.819885, Y: 91.45753}),
		QuadTo(f32.Point{X: 90.48111, Y: 96.70854}, f32.Point{X: 82.22082, Y: 99.568436}),
		QuadTo(f32.Point{X: 73.9478, Y: 102.40279}, f32.Point{X: 65.000046, Y: 102.4166}),
		QuadTo(f32.Point{X: 56.052288, Y: 102.40279}, f32.Point{X: 47.779274, Y: 99.568436}),
		QuadTo(f32.Point{X: 39.51898, Y: 96.70856}, f32.Point{X: 33.180206, Y: 91.45754}),
		QuadTo(f32.Point{X: 26.86493, Y: 86.18698}, f32.Point{X: 23.425415, Y: 79.318726}),
		QuadTo(f32.Point{X: 20.016602, Y: 72.439896}, f32.Point{X: 19.999983, Y: 65.00001}),
	},
}

A stroke/dash.go => stroke/dash.go +389 -0
@@ 0,0 1,389 @@
// SPDX-License-Identifier: Unlicense OR MIT

// The algorithms to compute dashes have been extracted, adapted from
// (and used as a reference implementation):
//  - github.com/tdewolff/canvas (Licensed under MIT)

package stroke

import (
	"math"
	"sort"

	"gioui.org/f32"
)

func isSolidLine(sty Dashes) bool {
	return sty.Phase == 0 && len(sty.Dashes) == 0
}

func (qs strokeQuads) dash(sty Dashes) strokeQuads {
	sty = dashCanonical(sty)

	switch {
	case len(sty.Dashes) == 0:
		return qs
	case len(sty.Dashes) == 1 && sty.Dashes[0] == 0.0:
		return strokeQuads{}
	}

	if len(sty.Dashes)%2 == 1 {
		// If the dash pattern is of uneven length, dash and space lengths
		// alternate. The following duplicates the pattern so that uneven
		// indices are always spaces.
		sty.Dashes = append(sty.Dashes, sty.Dashes...)
	}

	var (
		i0, pos0 = dashStart(sty)
		out      strokeQuads

		contour int = 1
	)

	for _, ps := range qs.split() {
		var (
			i      = i0
			pos    = pos0
			t      []float64
			length = ps.len()
		)
		for pos+sty.Dashes[i] < length {
			pos += sty.Dashes[i]
			if 0.0 < pos {
				t = append(t, float64(pos))
			}
			i++
			if i == len(sty.Dashes) {
				i = 0
			}
		}

		j0 := 0
		endsInDash := i%2 == 0
		if len(t)%2 == 1 && endsInDash || len(t)%2 == 0 && !endsInDash {
			j0 = 1
		}

		var (
			qd strokeQuads
			pd = ps.splitAt(&contour, t...)
		)
		for j := j0; j < len(pd)-1; j += 2 {
			qd = qd.append(pd[j])
		}
		if endsInDash {
			if ps.closed() {
				qd = pd[len(pd)-1].append(qd)
			} else {
				qd = qd.append(pd[len(pd)-1])
			}
		}
		out = out.append(qd)
		contour++
	}
	return out
}

func dashCanonical(sty Dashes) Dashes {
	var (
		o  = sty
		ds = o.Dashes
	)

	if len(sty.Dashes) == 0 {
		return sty
	}

	// Remove zeros except first and last.
	for i := 1; i < len(ds)-1; i++ {
		if f32Eq(ds[i], 0.0) {
			ds[i-1] += ds[i+1]
			ds = append(ds[:i], ds[i+2:]...)
			i--
		}
	}

	// Remove first zero, collapse with second and last.
	if f32Eq(ds[0], 0.0) {
		if len(ds) < 3 {
			return Dashes{
				Phase:  0.0,
				Dashes: []float32{0.0},
			}
		}
		o.Phase -= ds[1]
		ds[len(ds)-1] += ds[1]
		ds = ds[2:]
	}

	// Remove last zero, collapse with fist and second to last.
	if f32Eq(ds[len(ds)-1], 0.0) {
		if len(ds) < 3 {
			return Dashes{}
		}
		o.Phase += ds[len(ds)-2]
		ds[0] += ds[len(ds)-2]
		ds = ds[:len(ds)-2]
	}

	// If there are zeros or negatives, don't draw dashes.
	for i := 0; i < len(ds); i++ {
		if ds[i] < 0.0 || f32Eq(ds[i], 0.0) {
			return Dashes{
				Phase:  0.0,
				Dashes: []float32{0.0},
			}
		}
	}

	// Remove repeated patterns.
loop:
	for len(ds)%2 == 0 {
		mid := len(ds) / 2
		for i := 0; i < mid; i++ {
			if !f32Eq(ds[i], ds[mid+i]) {
				break loop
			}
		}
		ds = ds[:mid]
	}
	return o
}

func dashStart(sty Dashes) (int, float32) {
	i0 := 0 // i0 is the index into dashes.
	for sty.Dashes[i0] <= sty.Phase {
		sty.Phase -= sty.Dashes[i0]
		i0++
		if i0 == len(sty.Dashes) {
			i0 = 0
		}
	}
	// pos0 may be negative if the offset lands halfway into dash.
	pos0 := -sty.Phase
	if sty.Phase < 0.0 {
		var sum float32
		for _, d := range sty.Dashes {
			sum += d
		}
		pos0 = -(sum + sty.Phase) // handle negative offsets
	}
	return i0, pos0
}

func (qs strokeQuads) len() float32 {
	var sum float32
	for i := range qs {
		q := qs[i].Quad
		sum += quadBezierLen(q.From, q.Ctrl, q.To)
	}
	return sum
}

// splitAt splits the path into separate paths at the specified intervals
// along the path.
// splitAt updates the provided contour counter as it splits the segments.
func (qs strokeQuads) splitAt(contour *int, ts ...float64) []strokeQuads {
	if len(ts) == 0 {
		qs.setContour(*contour)
		return []strokeQuads{qs}
	}

	sort.Float64s(ts)
	if ts[0] == 0 {
		ts = ts[1:]
	}

	var (
		j int     // index into ts
		t float64 // current position along curve
	)

	var oo []strokeQuads
	var oi strokeQuads
	push := func() {
		oo = append(oo, oi)
		oi = nil
	}

	for _, ps := range qs.split() {
		for _, q := range ps {
			if j == len(ts) {
				oi = append(oi, q)
				continue
			}
			speed := func(t float64) float64 {
				return float64(lenPt(quadBezierD1(q.Quad.From, q.Quad.Ctrl, q.Quad.To, float32(t))))
			}
			invL, dt := invSpeedPolynomialChebyshevApprox(20, gaussLegendre7, speed, 0, 1)

			var (
				t0 float64
				r0 = q.Quad.From
				r1 = q.Quad.Ctrl
				r2 = q.Quad.To

				// from keeps track of the start of the 'running' segment.
				from = r0
			)
			for j < len(ts) && t < ts[j] && ts[j] <= t+dt {
				tj := invL(ts[j] - t)
				tsub := (tj - t0) / (1.0 - t0)
				t0 = tj

				var q1 f32.Point
				_, q1, _, r0, r1, r2 = quadBezierSplit(r0, r1, r2, float32(tsub))

				oi = append(oi, strokeQuad{
					Contour: *contour,
					Quad: quadSegment{
						From: from,
						Ctrl: q1,
						To:   r0,
					},
				})
				push()
				(*contour)++

				from = r0
				j++
			}
			if !f64Eq(t0, 1) {
				if len(oi) > 0 {
					r0 = oi.pen()
				}
				oi = append(oi, strokeQuad{
					Contour: *contour,
					Quad: quadSegment{
						From: r0,
						Ctrl: r1,
						To:   r2,
					},
				})
			}
			t += dt
		}
	}
	if len(oi) > 0 {
		push()
		(*contour)++
	}

	return oo
}

func f32Eq(a, b float32) bool {
	const epsilon = 1e-10
	return math.Abs(float64(a-b)) < epsilon
}

func f64Eq(a, b float64) bool {
	const epsilon = 1e-10
	return math.Abs(a-b) < epsilon
}

func invSpeedPolynomialChebyshevApprox(N int, gaussLegendre gaussLegendreFunc, fp func(float64) float64, tmin, tmax float64) (func(float64) float64, float64) {
	// The TODOs below are copied verbatim from tdewolff/canvas:
	//
	// TODO: find better way to determine N. For Arc 10 seems fine, for some
	// Quads 10 is too low, for Cube depending on inflection points is
	// maybe not the best indicator
	//
	// TODO: track efficiency, how many times is fp called?
	// Does a look-up table make more sense?
	fLength := func(t float64) float64 {
		return math.Abs(gaussLegendre(fp, tmin, t))
	}
	totalLength := fLength(tmax)
	t := func(L float64) float64 {
		return bisectionMethod(fLength, L, tmin, tmax)
	}
	return polynomialChebyshevApprox(N, t, 0.0, totalLength, tmin, tmax), totalLength
}

func polynomialChebyshevApprox(N int, f func(float64) float64, xmin, xmax, ymin, ymax float64) func(float64) float64 {
	var (
		invN = 1.0 / float64(N)
		fs   = make([]float64, N)
	)
	for k := 0; k < N; k++ {
		u := math.Cos(math.Pi * (float64(k+1) - 0.5) * invN)
		fs[k] = f(xmin + 0.5*(xmax-xmin)*(u+1))
	}

	c := make([]float64, N)
	for j := 0; j < N; j++ {
		var a float64
		for k := 0; k < N; k++ {
			a += fs[k] * math.Cos(float64(j)*math.Pi*(float64(k+1)-0.5)/float64(N))
		}
		c[j] = 2 * invN * a
	}

	if ymax < ymin {
		ymin, ymax = ymax, ymin
	}
	return func(x float64) float64 {
		x = math.Min(xmax, math.Max(xmin, x))
		u := (x-xmin)/(xmax-xmin)*2 - 1
		var a float64
		for j := 0; j < N; j++ {
			a += c[j] * math.Cos(float64(j)*math.Acos(u))
		}
		y := -0.5*c[0] + a
		if !math.IsNaN(ymin) && !math.IsNaN(ymax) {
			y = math.Min(ymax, math.Max(ymin, y))
		}
		return y
	}
}

// bisectionMethod finds the value x for which f(x) = y in the interval x
// in [xmin, xmax] using the bisection method.
func bisectionMethod(f func(float64) float64, y, xmin, xmax float64) float64 {
	const (
		maxIter   = 100
		tolerance = 0.001 // 0.1%
	)

	var (
		n    = 0
		x    float64
		tolX = math.Abs(xmax-xmin) * tolerance
		tolY = math.Abs(f(xmax)-f(xmin)) * tolerance
	)
	for {
		x = 0.5 * (xmin + xmax)
		if n >= maxIter {
			return x
		}

		dy := f(x) - y
		switch {
		case math.Abs(dy) < tolY, math.Abs(0.5*(xmax-xmin)) < tolX:
			return x
		case dy > 0:
			xmax = x
		default:
			xmin = x
		}
		n++
	}
}

type gaussLegendreFunc func(func(float64) float64, float64, float64) float64

// Gauss-Legendre quadrature integration from a to b with n=7
func gaussLegendre7(f func(float64) float64, a, b float64) float64 {
	c := 0.5 * (b - a)
	d := 0.5 * (a + b)
	Qd1 := f(-0.949108*c + d)
	Qd2 := f(-0.741531*c + d)
	Qd3 := f(-0.405845*c + d)
	Qd4 := f(d)
	Qd5 := f(0.405845*c + d)
	Qd6 := f(0.741531*c + d)
	Qd7 := f(0.949108*c + d)
	return c * (0.129485*(Qd1+Qd7) + 0.279705*(Qd2+Qd6) + 0.381830*(Qd3+Qd5) + 0.417959*Qd4)
}

A stroke/refs/TestDashedPathFlatCapEllipse.png => stroke/refs/TestDashedPathFlatCapEllipse.png +0 -0
A stroke/refs/TestDashedPathFlatCapZ.png => stroke/refs/TestDashedPathFlatCapZ.png +0 -0
A stroke/refs/TestDashedPathFlatCapZNoDash.png => stroke/refs/TestDashedPathFlatCapZNoDash.png +0 -0
A stroke/refs/TestDashedPathFlatCapZNoPath.png => stroke/refs/TestDashedPathFlatCapZNoPath.png +0 -0
A stroke/refs/TestStrokedPathBevelFlat.png => stroke/refs/TestStrokedPathBevelFlat.png +0 -0
A stroke/refs/TestStrokedPathBevelRound.png => stroke/refs/TestStrokedPathBevelRound.png +0 -0
A stroke/refs/TestStrokedPathBevelSquare.png => stroke/refs/TestStrokedPathBevelSquare.png +0 -0
A stroke/refs/TestStrokedPathFlatMiter.png => stroke/refs/TestStrokedPathFlatMiter.png +0 -0
A stroke/refs/TestStrokedPathFlatMiterInf.png => stroke/refs/TestStrokedPathFlatMiterInf.png +0 -0
A stroke/refs/TestStrokedPathZeroWidth.png => stroke/refs/TestStrokedPathZeroWidth.png +0 -0
A stroke/stroke.go => stroke/stroke.go +140 -0
@@ 0,0 1,140 @@
// SPDX-License-Identifier: Unlicense OR MIT

// Package stroke converts complex strokes to gioui.org/op/clip operations.
package stroke

import (
	"gioui.org/f32"
	"gioui.org/op"
	"gioui.org/op/clip"
)

// Path defines the shape of a Stroke.
type Path struct {
	Segments []Segment
}

type Segment struct {
	// op is the operator.
	op segmentOp
	// args is up to three (x, y) coordinates.
	args [3]f32.Point
}

// Dashes defines the dash pattern of a Stroke.
type Dashes struct {
	Phase  float32
	Dashes []float32
}

// Stroke defines a stroke.
type Stroke struct {
	Path  Path
	Width float32 // Width of the stroked path.

	// Miter is the limit to apply to a miter joint.
	// The zero Miter disables the miter joint; setting Miter to +∞
	// unconditionally enables the miter joint.
	Miter float32
	Cap   StrokeCap  // Cap describes the head or tail of a stroked path.
	Join  StrokeJoin // Join describes how stroked paths are collated.

	Dashes Dashes
}

type segmentOp uint8

const (
	segOpMoveTo segmentOp = iota
	segOpLineTo
	segOpQuadTo
	segOpCubeTo
)

// StrokeCap describes the head or tail of a stroked path.
type StrokeCap uint8

const (
	// RoundCap caps stroked paths with a round cap, joining the right-hand and
	// left-hand sides of a stroked path with a half disc of diameter the
	// stroked path's width.
	RoundCap StrokeCap = iota

	// FlatCap caps stroked paths with a flat cap, joining the right-hand
	// and left-hand sides of a stroked path with a straight line.
	FlatCap

	// SquareCap caps stroked paths with a square cap, joining the right-hand
	// and left-hand sides of a stroked path with a half square of length
	// the stroked path's width.
	SquareCap
)

// StrokeJoin describes how stroked paths are collated.
type StrokeJoin uint8

const (
	// RoundJoin joins path segments with a round segment.
	RoundJoin StrokeJoin = iota

	// BevelJoin joins path segments with sharp bevels.
	BevelJoin
)

func MoveTo(p f32.Point) Segment {
	s := Segment{
		op: segOpMoveTo,
	}
	s.args[0] = p
	return s
}

func LineTo(p f32.Point) Segment {
	s := Segment{
		op: segOpLineTo,
	}
	s.args[0] = p
	return s
}

func QuadTo(ctrl, end f32.Point) Segment {
	s := Segment{
		op: segOpQuadTo,
	}
	s.args[0] = ctrl
	s.args[1] = end
	return s
}

func CubeTo(ctrl0, ctrl1, end f32.Point) Segment {
	s := Segment{
		op: segOpCubeTo,
	}
	s.args[0] = ctrl0
	s.args[1] = ctrl1
	s.args[2] = end
	return s
}

// Op returns a clip operation that approximates stroke.
func (s Stroke) Op(ops *op.Ops) clip.Op {
	if len(s.Path.Segments) == 0 {
		return clip.Op{}
	}

	// Approximate and output path data.
	var outline clip.Path
	outline.Begin(ops)
	quads := strokePathCommands(s)
	pen := f32.Pt(0, 0)
	for _, quad := range quads {
		q := quad.Quad
		if q.From != pen {
			pen = q.From
			outline.MoveTo(pen)
		}
		outline.QuadTo(q.Ctrl, q.To)
		pen = q.To
	}
	return clip.Outline{Path: outline.End()}.Op()
}

A stroke/util_test.go => stroke/util_test.go +205 -0
@@ 0,0 1,205 @@
// SPDX-License-Identifier: Unlicense OR MIT

package stroke

import (
	"bytes"
	"flag"
	"fmt"
	"image"
	"image/color"
	"image/png"
	"io/ioutil"
	"path/filepath"
	"strconv"
	"testing"

	"golang.org/x/image/colornames"

	"gioui.org/gpu/headless"
	"gioui.org/internal/f32color"
	"gioui.org/op"
)

var (
	dumpImages = flag.Bool("saveimages", false, "save test images")
)

var (
	red         = f32color.RGBAToNRGBA(colornames.Red)
	black       = f32color.RGBAToNRGBA(colornames.Black)
	transparent = color.RGBA{}
)

func drawImage(t *testing.T, size int, ops *op.Ops, draw func(o *op.Ops)) (im *image.RGBA, err error) {
	sz := image.Point{X: size, Y: size}
	w := newWindow(t, sz.X, sz.Y)
	defer w.Release()
	draw(ops)
	if err := w.Frame(ops); err != nil {
		return nil, err
	}
	return w.Screenshot()
}

func run(t *testing.T, f func(o *op.Ops), c func(r result)) {
	// Draw a few times and check that it is correct each time, to
	// ensure any caching effects still generate the correct images.
	var img *image.RGBA
	var err error
	ops := new(op.Ops)
	for i := 0; i < 3; i++ {
		ops.Reset()
		img, err = drawImage(t, 128, ops, f)
		if err != nil {
			t.Error("error rendering:", err)
			return
		}
		// Check for a reference image and make sure it is identical.
		if !verifyRef(t, img, 0) {
			name := fmt.Sprintf("%s-%d-bad.png", t.Name(), i)
			saveImage(t, name, img)
		}
		c(result{t: t, img: img})
	}
}

func verifyRef(t *testing.T, img *image.RGBA, frame int) (ok bool) {
	// ensure identical to ref data
	path := filepath.Join("refs", t.Name()+".png")
	if frame != 0 {
		path = filepath.Join("refs", t.Name()+"_"+strconv.Itoa(frame)+".png")
	}
	b, err := ioutil.ReadFile(path)
	if err != nil {
		t.Error("could not open ref:", err)
		return
	}
	r, err := png.Decode(bytes.NewReader(b))
	if err != nil {
		t.Error("could not decode ref:", err)
		return
	}
	if img.Bounds() != r.Bounds() {
		t.Errorf("reference image is %v, expected %v", r.Bounds(), img.Bounds())
		return false
	}
	var ref *image.RGBA
	switch r := r.(type) {
	case *image.RGBA:
		ref = r
	case *image.NRGBA:
		ref = image.NewRGBA(r.Bounds())
		bnd := r.Bounds()
		for x := bnd.Min.X; x < bnd.Max.X; x++ {
			for y := bnd.Min.Y; y < bnd.Max.Y; y++ {
				ref.SetRGBA(x, y, f32color.NRGBAToRGBA(r.NRGBAAt(x, y)))
			}
		}
	default:
		t.Fatalf("reference image is a %T, expected *image.NRGBA or *image.RGBA", r)
	}
	bnd := img.Bounds()
	for x := bnd.Min.X; x < bnd.Max.X; x++ {
		for y := bnd.Min.Y; y < bnd.Max.Y; y++ {
			exp := ref.RGBAAt(x, y)
			got := img.RGBAAt(x, y)
			if !colorsClose(exp, got) {
				t.Error("not equal to ref at", x, y, " ", got, exp)
				return false
			}
		}
	}
	return true
}

func colorsClose(c1, c2 color.RGBA) bool {
	const delta = 0.01 // magic value obtained from experimentation.
	return yiqEqApprox(c1, c2, delta)
}

// yiqEqApprox compares the colors of 2 pixels, in the NTSC YIQ color space,
// as described in:
//
//   Measuring perceived color difference using YIQ NTSC
//   transmission color space in mobile applications.
//   Yuriy Kotsarenko, Fernando Ramos.
//
// An electronic version is available at:
//
// - http://www.progmat.uaem.mx:8080/artVol2Num2/Articulo3Vol2Num2.pdf
func yiqEqApprox(c1, c2 color.RGBA, d2 float64) bool {
	const max = 35215.0 // difference between 2 maximally different pixels.

	var (
		r1 = float64(c1.R)
		g1 = float64(c1.G)
		b1 = float64(c1.B)

		r2 = float64(c2.R)
		g2 = float64(c2.G)
		b2 = float64(c2.B)

		y1 = r1*0.29889531 + g1*0.58662247 + b1*0.11448223
		i1 = r1*0.59597799 - g1*0.27417610 - b1*0.32180189
		q1 = r1*0.21147017 - g1*0.52261711 + b1*0.31114694

		y2 = r2*0.29889531 + g2*0.58662247 + b2*0.11448223
		i2 = r2*0.59597799 - g2*0.27417610 - b2*0.32180189
		q2 = r2*0.21147017 - g2*0.52261711 + b2*0.31114694

		y = y1 - y2
		i = i1 - i2
		q = q1 - q2

		diff = 0.5053*y*y + 0.299*i*i + 0.1957*q*q
	)
	return diff <= max*d2
}

func (r result) expect(x, y int, col color.RGBA) {
	r.t.Helper()
	if r.img == nil {
		return
	}
	c := r.img.RGBAAt(x, y)
	if !colorsClose(c, col) {
		r.t.Error("expected ", col, " at ", "(", x, ",", y, ") but got ", c)
	}
}

type result struct {
	t   *testing.T
	img *image.RGBA
}

func saveImage(t testing.TB, file string, img *image.RGBA) {
	if !*dumpImages {
		return
	}
	// Only NRGBA images are losslessly encoded by png.Encode.
	nrgba := image.NewNRGBA(img.Bounds())
	bnd := img.Bounds()
	for x := bnd.Min.X; x < bnd.Max.X; x++ {
		for y := bnd.Min.Y; y < bnd.Max.Y; y++ {
			nrgba.SetNRGBA(x, y, f32color.RGBAToNRGBA(img.RGBAAt(x, y)))
		}
	}
	var buf bytes.Buffer
	if err := png.Encode(&buf, nrgba); err != nil {
		t.Error(err)
		return
	}
	if err := ioutil.WriteFile(file, buf.Bytes(), 0666); err != nil {
		t.Error(err)
		return
	}
}

func newWindow(t testing.TB, width, height int) *headless.Window {
	w, err := headless.NewWindow(width, height)
	if err != nil {
		t.Skipf("failed to create headless window, skipping: %v", err)
	}
	return w
}