dash.go raw

   1  // SPDX-License-Identifier: Unlicense OR MIT
   2  
   3  // The algorithms to compute dashes have been extracted, adapted from
   4  // (and used as a reference implementation):
   5  //  - github.com/tdewolff/canvas (Licensed under MIT)
   6  
   7  package stroke
   8  
   9  import (
  10  	"math"
  11  	"sort"
  12  
  13  	"github.com/p9c/p9/pkg/gel/gio/f32"
  14  )
  15  
  16  type DashOp struct {
  17  	Phase  float32
  18  	Dashes []float32
  19  }
  20  
  21  func IsSolidLine(sty DashOp) bool {
  22  	return sty.Phase == 0 && len(sty.Dashes) == 0
  23  }
  24  
  25  func (qs StrokeQuads) dash(sty DashOp) StrokeQuads {
  26  	sty = dashCanonical(sty)
  27  
  28  	switch {
  29  	case len(sty.Dashes) == 0:
  30  		return qs
  31  	case len(sty.Dashes) == 1 && sty.Dashes[0] == 0.0:
  32  		return StrokeQuads{}
  33  	}
  34  
  35  	if len(sty.Dashes)%2 == 1 {
  36  		// If the dash pattern is of uneven length, dash and space lengths
  37  		// alternate. The following duplicates the pattern so that uneven
  38  		// indices are always spaces.
  39  		sty.Dashes = append(sty.Dashes, sty.Dashes...)
  40  	}
  41  
  42  	var (
  43  		i0, pos0 = dashStart(sty)
  44  		out      StrokeQuads
  45  
  46  		contour uint32 = 1
  47  	)
  48  
  49  	for _, ps := range qs.split() {
  50  		var (
  51  			i      = i0
  52  			pos    = pos0
  53  			t      []float64
  54  			length = ps.len()
  55  		)
  56  		for pos+sty.Dashes[i] < length {
  57  			pos += sty.Dashes[i]
  58  			if 0.0 < pos {
  59  				t = append(t, float64(pos))
  60  			}
  61  			i++
  62  			if i == len(sty.Dashes) {
  63  				i = 0
  64  			}
  65  		}
  66  
  67  		j0 := 0
  68  		endsInDash := i%2 == 0
  69  		if len(t)%2 == 1 && endsInDash || len(t)%2 == 0 && !endsInDash {
  70  			j0 = 1
  71  		}
  72  
  73  		var (
  74  			qd StrokeQuads
  75  			pd = ps.splitAt(&contour, t...)
  76  		)
  77  		for j := j0; j < len(pd)-1; j += 2 {
  78  			qd = qd.append(pd[j])
  79  		}
  80  		if endsInDash {
  81  			if ps.closed() {
  82  				qd = pd[len(pd)-1].append(qd)
  83  			} else {
  84  				qd = qd.append(pd[len(pd)-1])
  85  			}
  86  		}
  87  		out = out.append(qd)
  88  		contour++
  89  	}
  90  	return out
  91  }
  92  
  93  func dashCanonical(sty DashOp) DashOp {
  94  	var (
  95  		o  = sty
  96  		ds = o.Dashes
  97  	)
  98  
  99  	if len(sty.Dashes) == 0 {
 100  		return sty
 101  	}
 102  
 103  	// Remove zeros except first and last.
 104  	for i := 1; i < len(ds)-1; i++ {
 105  		if f32Eq(ds[i], 0.0) {
 106  			ds[i-1] += ds[i+1]
 107  			ds = append(ds[:i], ds[i+2:]...)
 108  			i--
 109  		}
 110  	}
 111  
 112  	// Remove first zero, collapse with second and last.
 113  	if f32Eq(ds[0], 0.0) {
 114  		if len(ds) < 3 {
 115  			return DashOp{
 116  				Phase:  0.0,
 117  				Dashes: []float32{0.0},
 118  			}
 119  		}
 120  		o.Phase -= ds[1]
 121  		ds[len(ds)-1] += ds[1]
 122  		ds = ds[2:]
 123  	}
 124  
 125  	// Remove last zero, collapse with fist and second to last.
 126  	if f32Eq(ds[len(ds)-1], 0.0) {
 127  		if len(ds) < 3 {
 128  			return DashOp{}
 129  		}
 130  		o.Phase += ds[len(ds)-2]
 131  		ds[0] += ds[len(ds)-2]
 132  		ds = ds[:len(ds)-2]
 133  	}
 134  
 135  	// If there are zeros or negatives, don't draw dashes.
 136  	for i := 0; i < len(ds); i++ {
 137  		if ds[i] < 0.0 || f32Eq(ds[i], 0.0) {
 138  			return DashOp{
 139  				Phase:  0.0,
 140  				Dashes: []float32{0.0},
 141  			}
 142  		}
 143  	}
 144  
 145  	// Remove repeated patterns.
 146  loop:
 147  	for len(ds)%2 == 0 {
 148  		mid := len(ds) / 2
 149  		for i := 0; i < mid; i++ {
 150  			if !f32Eq(ds[i], ds[mid+i]) {
 151  				break loop
 152  			}
 153  		}
 154  		ds = ds[:mid]
 155  	}
 156  	return o
 157  }
 158  
 159  func dashStart(sty DashOp) (int, float32) {
 160  	i0 := 0 // i0 is the index into dashes.
 161  	for sty.Dashes[i0] <= sty.Phase {
 162  		sty.Phase -= sty.Dashes[i0]
 163  		i0++
 164  		if i0 == len(sty.Dashes) {
 165  			i0 = 0
 166  		}
 167  	}
 168  	// pos0 may be negative if the offset lands halfway into dash.
 169  	pos0 := -sty.Phase
 170  	if sty.Phase < 0.0 {
 171  		var sum float32
 172  		for _, d := range sty.Dashes {
 173  			sum += d
 174  		}
 175  		pos0 = -(sum + sty.Phase) // handle negative offsets
 176  	}
 177  	return i0, pos0
 178  }
 179  
 180  func (qs StrokeQuads) len() float32 {
 181  	var sum float32
 182  	for i := range qs {
 183  		q := qs[i].Quad
 184  		sum += quadBezierLen(q.From, q.Ctrl, q.To)
 185  	}
 186  	return sum
 187  }
 188  
 189  // splitAt splits the path into separate paths at the specified intervals
 190  // along the path.
 191  // splitAt updates the provided contour counter as it splits the segments.
 192  func (qs StrokeQuads) splitAt(contour *uint32, ts ...float64) []StrokeQuads {
 193  	if len(ts) == 0 {
 194  		qs.setContour(*contour)
 195  		return []StrokeQuads{qs}
 196  	}
 197  
 198  	sort.Float64s(ts)
 199  	if ts[0] == 0 {
 200  		ts = ts[1:]
 201  	}
 202  
 203  	var (
 204  		j int     // index into ts
 205  		t float64 // current position along curve
 206  	)
 207  
 208  	var oo []StrokeQuads
 209  	var oi StrokeQuads
 210  	push := func() {
 211  		oo = append(oo, oi)
 212  		oi = nil
 213  	}
 214  
 215  	for _, ps := range qs.split() {
 216  		for _, q := range ps {
 217  			if j == len(ts) {
 218  				oi = append(oi, q)
 219  				continue
 220  			}
 221  			speed := func(t float64) float64 {
 222  				return float64(lenPt(quadBezierD1(q.Quad.From, q.Quad.Ctrl, q.Quad.To, float32(t))))
 223  			}
 224  			invL, dt := invSpeedPolynomialChebyshevApprox(20, gaussLegendre7, speed, 0, 1)
 225  
 226  			var (
 227  				t0 float64
 228  				r0 = q.Quad.From
 229  				r1 = q.Quad.Ctrl
 230  				r2 = q.Quad.To
 231  
 232  				// from keeps track of the start of the 'running' segment.
 233  				from = r0
 234  			)
 235  			for j < len(ts) && t < ts[j] && ts[j] <= t+dt {
 236  				tj := invL(ts[j] - t)
 237  				tsub := (tj - t0) / (1.0 - t0)
 238  				t0 = tj
 239  
 240  				var q1 f32.Point
 241  				_, q1, _, r0, r1, r2 = quadBezierSplit(r0, r1, r2, float32(tsub))
 242  
 243  				oi = append(oi, StrokeQuad{
 244  					Contour: *contour,
 245  					Quad: QuadSegment{
 246  						From: from,
 247  						Ctrl: q1,
 248  						To:   r0,
 249  					},
 250  				})
 251  				push()
 252  				(*contour)++
 253  
 254  				from = r0
 255  				j++
 256  			}
 257  			if !f64Eq(t0, 1) {
 258  				if len(oi) > 0 {
 259  					r0 = oi.pen()
 260  				}
 261  				oi = append(oi, StrokeQuad{
 262  					Contour: *contour,
 263  					Quad: QuadSegment{
 264  						From: r0,
 265  						Ctrl: r1,
 266  						To:   r2,
 267  					},
 268  				})
 269  			}
 270  			t += dt
 271  		}
 272  	}
 273  	if len(oi) > 0 {
 274  		push()
 275  		(*contour)++
 276  	}
 277  
 278  	return oo
 279  }
 280  
 281  func f32Eq(a, b float32) bool {
 282  	const epsilon = 1e-10
 283  	return math.Abs(float64(a-b)) < epsilon
 284  }
 285  
 286  func f64Eq(a, b float64) bool {
 287  	const epsilon = 1e-10
 288  	return math.Abs(a-b) < epsilon
 289  }
 290  
 291  func invSpeedPolynomialChebyshevApprox(N int, gaussLegendre gaussLegendreFunc, fp func(float64) float64, tmin, tmax float64) (func(float64) float64, float64) {
 292  	// The TODOs below are copied verbatim from tdewolff/canvas:
 293  	//
 294  	// TODO: find better way to determine N. For Arc 10 seems fine, for some
 295  	// Quads 10 is too low, for Cube depending on inflection points is
 296  	// maybe not the best indicator
 297  	//
 298  	// TODO: track efficiency, how many times is fp called?
 299  	// Does a look-up table make more sense?
 300  	fLength := func(t float64) float64 {
 301  		return math.Abs(gaussLegendre(fp, tmin, t))
 302  	}
 303  	totalLength := fLength(tmax)
 304  	t := func(L float64) float64 {
 305  		return bisectionMethod(fLength, L, tmin, tmax)
 306  	}
 307  	return polynomialChebyshevApprox(N, t, 0.0, totalLength, tmin, tmax), totalLength
 308  }
 309  
 310  func polynomialChebyshevApprox(N int, f func(float64) float64, xmin, xmax, ymin, ymax float64) func(float64) float64 {
 311  	var (
 312  		invN = 1.0 / float64(N)
 313  		fs   = make([]float64, N)
 314  	)
 315  	for k := 0; k < N; k++ {
 316  		u := math.Cos(math.Pi * (float64(k+1) - 0.5) * invN)
 317  		fs[k] = f(xmin + 0.5*(xmax-xmin)*(u+1))
 318  	}
 319  
 320  	c := make([]float64, N)
 321  	for j := 0; j < N; j++ {
 322  		var a float64
 323  		for k := 0; k < N; k++ {
 324  			a += fs[k] * math.Cos(float64(j)*math.Pi*(float64(k+1)-0.5)/float64(N))
 325  		}
 326  		c[j] = 2 * invN * a
 327  	}
 328  
 329  	if ymax < ymin {
 330  		ymin, ymax = ymax, ymin
 331  	}
 332  	return func(x float64) float64 {
 333  		x = math.Min(xmax, math.Max(xmin, x))
 334  		u := (x-xmin)/(xmax-xmin)*2 - 1
 335  		var a float64
 336  		for j := 0; j < N; j++ {
 337  			a += c[j] * math.Cos(float64(j)*math.Acos(u))
 338  		}
 339  		y := -0.5*c[0] + a
 340  		if !math.IsNaN(ymin) && !math.IsNaN(ymax) {
 341  			y = math.Min(ymax, math.Max(ymin, y))
 342  		}
 343  		return y
 344  	}
 345  }
 346  
 347  // bisectionMethod finds the value x for which f(x) = y in the interval x
 348  // in [xmin, xmax] using the bisection method.
 349  func bisectionMethod(f func(float64) float64, y, xmin, xmax float64) float64 {
 350  	const (
 351  		maxIter   = 100
 352  		tolerance = 0.001 // 0.1%
 353  	)
 354  
 355  	var (
 356  		n    = 0
 357  		x    float64
 358  		tolX = math.Abs(xmax-xmin) * tolerance
 359  		tolY = math.Abs(f(xmax)-f(xmin)) * tolerance
 360  	)
 361  	for {
 362  		x = 0.5 * (xmin + xmax)
 363  		if n >= maxIter {
 364  			return x
 365  		}
 366  
 367  		dy := f(x) - y
 368  		switch {
 369  		case math.Abs(dy) < tolY, math.Abs(0.5*(xmax-xmin)) < tolX:
 370  			return x
 371  		case dy > 0:
 372  			xmax = x
 373  		default:
 374  			xmin = x
 375  		}
 376  		n++
 377  	}
 378  }
 379  
 380  type gaussLegendreFunc func(func(float64) float64, float64, float64) float64
 381  
 382  // Gauss-Legendre quadrature integration from a to b with n=7
 383  func gaussLegendre7(f func(float64) float64, a, b float64) float64 {
 384  	c := 0.5 * (b - a)
 385  	d := 0.5 * (a + b)
 386  	Qd1 := f(-0.949108*c + d)
 387  	Qd2 := f(-0.741531*c + d)
 388  	Qd3 := f(-0.405845*c + d)
 389  	Qd4 := f(d)
 390  	Qd5 := f(0.405845*c + d)
 391  	Qd6 := f(0.741531*c + d)
 392  	Qd7 := f(0.949108*c + d)
 393  	return c * (0.129485*(Qd1+Qd7) + 0.279705*(Qd2+Qd6) + 0.381830*(Qd3+Qd5) + 0.417959*Qd4)
 394  }
 395