stroke.go raw
1 // SPDX-License-Identifier: Unlicense OR MIT
2
3 // Most of the algorithms to compute strokes and their offsets have been
4 // extracted, adapted from (and used as a reference implementation):
5 // - github.com/tdewolff/canvas (Licensed under MIT)
6 //
7 // These algorithms have been implemented from:
8 // Fast, precise flattening of cubic Bézier path and offset curves
9 // Thomas F. Hain, et al.
10 //
11 // An electronic version is available at:
12 // https://seant23.files.wordpress.com/2010/11/fastpreciseflatteningofbeziercurve.pdf
13 //
14 // Possible improvements (in term of speed and/or accuracy) on these
15 // algorithms are:
16 //
17 // - Polar Stroking: New Theory and Methods for Stroking Paths,
18 // M. Kilgard
19 // https://arxiv.org/pdf/2007.00308.pdf
20 //
21 // - https://raphlinus.github.io/graphics/curves/2019/12/23/flatten-quadbez.html
22 // R. Levien
23
24 // Package stroke implements conversion of strokes to filled outlines. It is used as a
25 // fallback for stroke configurations not natively supported by the renderer.
26 package stroke
27
28 import (
29 "encoding/binary"
30 "math"
31
32 "github.com/p9c/p9/pkg/gel/gio/f32"
33 "github.com/p9c/p9/pkg/gel/gio/internal/ops"
34 "github.com/p9c/p9/pkg/gel/gio/internal/scene"
35 )
36
37 // The following are copies of types from op/clip to avoid a circular import of
38 // that package.
39 // TODO: when the old renderer is gone, this package can be merged with
40 // op/clip, eliminating the duplicate types.
41 type StrokeStyle struct {
42 Width float32
43 Miter float32
44 Cap StrokeCap
45 Join StrokeJoin
46 }
47
48 type StrokeCap uint8
49
50 const (
51 RoundCap StrokeCap = iota
52 FlatCap
53 SquareCap
54 )
55
56 type StrokeJoin uint8
57
58 const (
59 RoundJoin StrokeJoin = iota
60 BevelJoin
61 )
62
63 // strokeTolerance is used to reconcile rounding errors arising
64 // when splitting quads into smaller and smaller segments to approximate
65 // them into straight lines, and when joining back segments.
66 //
67 // The magic value of 0.01 was found by striking a compromise between
68 // aesthetic looking (curves did look like curves, even after linearization)
69 // and speed.
70 const strokeTolerance = 0.01
71
72 type QuadSegment struct {
73 From, Ctrl, To f32.Point
74 }
75
76 type StrokeQuad struct {
77 Contour uint32
78 Quad QuadSegment
79 }
80
81 type strokeState struct {
82 p0, p1 f32.Point // p0 is the start point, p1 the end point.
83 n0, n1 f32.Point // n0 is the normal vector at the start point, n1 at the end point.
84 r0, r1 float32 // r0 is the curvature at the start point, r1 at the end point.
85 ctl f32.Point // ctl is the control point of the quadratic Bézier segment.
86 }
87
88 type StrokeQuads []StrokeQuad
89
90 func (qs *StrokeQuads) setContour(n uint32) {
91 for i := range *qs {
92 (*qs)[i].Contour = n
93 }
94 }
95
96 func (qs *StrokeQuads) pen() f32.Point {
97 return (*qs)[len(*qs)-1].Quad.To
98 }
99
100 func (qs *StrokeQuads) closed() bool {
101 beg := (*qs)[0].Quad.From
102 end := (*qs)[len(*qs)-1].Quad.To
103 return f32Eq(beg.X, end.X) && f32Eq(beg.Y, end.Y)
104 }
105
106 func (qs *StrokeQuads) lineTo(pt f32.Point) {
107 end := qs.pen()
108 *qs = append(*qs, StrokeQuad{
109 Quad: QuadSegment{
110 From: end,
111 Ctrl: end.Add(pt).Mul(0.5),
112 To: pt,
113 },
114 })
115 }
116
117 func (qs *StrokeQuads) arc(f1, f2 f32.Point, angle float32) {
118 const segments = 16
119 pen := qs.pen()
120 m := ArcTransform(pen, f1.Add(pen), f2.Add(pen), angle, segments)
121 for i := 0; i < segments; i++ {
122 p0 := qs.pen()
123 p1 := m.Transform(p0)
124 p2 := m.Transform(p1)
125 ctl := p1.Mul(2).Sub(p0.Add(p2).Mul(.5))
126 *qs = append(*qs, StrokeQuad{
127 Quad: QuadSegment{
128 From: p0, Ctrl: ctl, To: p2,
129 },
130 })
131 }
132 }
133
134 // split splits a slice of quads into slices of quads grouped
135 // by contours (ie: splitted at move-to boundaries).
136 func (qs StrokeQuads) split() []StrokeQuads {
137 if len(qs) == 0 {
138 return nil
139 }
140
141 var (
142 c uint32
143 o []StrokeQuads
144 i = len(o)
145 )
146 for _, q := range qs {
147 if q.Contour != c {
148 c = q.Contour
149 i = len(o)
150 o = append(o, StrokeQuads{})
151 }
152 o[i] = append(o[i], q)
153 }
154
155 return o
156 }
157
158 func (qs StrokeQuads) stroke(stroke StrokeStyle, dashes DashOp) StrokeQuads {
159 if !IsSolidLine(dashes) {
160 qs = qs.dash(dashes)
161 }
162
163 var (
164 o StrokeQuads
165 hw = 0.5 * stroke.Width
166 )
167
168 for _, ps := range qs.split() {
169 rhs, lhs := ps.offset(hw, stroke)
170 switch lhs {
171 case nil:
172 o = o.append(rhs)
173 default:
174 // Closed path.
175 // Inner path should go opposite direction to cancel outer path.
176 switch {
177 case ps.ccw():
178 lhs = lhs.reverse()
179 o = o.append(rhs)
180 o = o.append(lhs)
181 default:
182 rhs = rhs.reverse()
183 o = o.append(lhs)
184 o = o.append(rhs)
185 }
186 }
187 }
188
189 return o
190 }
191
192 // offset returns the right-hand and left-hand sides of the path, offset by
193 // the half-width hw.
194 // The stroke handles how segments are joined and ends are capped.
195 func (qs StrokeQuads) offset(hw float32, stroke StrokeStyle) (rhs, lhs StrokeQuads) {
196 var (
197 states []strokeState
198 beg = qs[0].Quad.From
199 end = qs[len(qs)-1].Quad.To
200 closed = beg == end
201 )
202 for i := range qs {
203 q := qs[i].Quad
204
205 var (
206 n0 = strokePathNorm(q.From, q.Ctrl, q.To, 0, hw)
207 n1 = strokePathNorm(q.From, q.Ctrl, q.To, 1, hw)
208 r0 = strokePathCurv(q.From, q.Ctrl, q.To, 0)
209 r1 = strokePathCurv(q.From, q.Ctrl, q.To, 1)
210 )
211 states = append(states, strokeState{
212 p0: q.From,
213 p1: q.To,
214 n0: n0,
215 n1: n1,
216 r0: r0,
217 r1: r1,
218 ctl: q.Ctrl,
219 })
220 }
221
222 for i, state := range states {
223 rhs = rhs.append(strokeQuadBezier(state, +hw, strokeTolerance))
224 lhs = lhs.append(strokeQuadBezier(state, -hw, strokeTolerance))
225
226 // join the current and next segments
227 if hasNext := i+1 < len(states); hasNext || closed {
228 var next strokeState
229 switch {
230 case hasNext:
231 next = states[i+1]
232 case closed:
233 next = states[0]
234 }
235 if state.n1 != next.n0 {
236 strokePathJoin(stroke, &rhs, &lhs, hw, state.p1, state.n1, next.n0, state.r1, next.r0)
237 }
238 }
239 }
240
241 if closed {
242 rhs.close()
243 lhs.close()
244 return rhs, lhs
245 }
246
247 qbeg := &states[0]
248 qend := &states[len(states)-1]
249
250 // Default to counter-clockwise direction.
251 lhs = lhs.reverse()
252 strokePathCap(stroke, &rhs, hw, qend.p1, qend.n1)
253
254 rhs = rhs.append(lhs)
255 strokePathCap(stroke, &rhs, hw, qbeg.p0, qbeg.n0.Mul(-1))
256
257 rhs.close()
258
259 return rhs, nil
260 }
261
262 func (qs *StrokeQuads) close() {
263 p0 := (*qs)[len(*qs)-1].Quad.To
264 p1 := (*qs)[0].Quad.From
265
266 if p1 == p0 {
267 return
268 }
269
270 *qs = append(*qs, StrokeQuad{
271 Quad: QuadSegment{
272 From: p0,
273 Ctrl: p0.Add(p1).Mul(0.5),
274 To: p1,
275 },
276 })
277 }
278
279 // ccw returns whether the path is counter-clockwise.
280 func (qs StrokeQuads) ccw() bool {
281 // Use the Shoelace formula:
282 // https://en.wikipedia.org/wiki/Shoelace_formula
283 var area float32
284 for _, ps := range qs.split() {
285 for i := 1; i < len(ps); i++ {
286 pi := ps[i].Quad.To
287 pj := ps[i-1].Quad.To
288 area += (pi.X - pj.X) * (pi.Y + pj.Y)
289 }
290 }
291 return area <= 0.0
292 }
293
294 func (qs StrokeQuads) reverse() StrokeQuads {
295 if len(qs) == 0 {
296 return nil
297 }
298
299 ps := make(StrokeQuads, 0, len(qs))
300 for i := range qs {
301 q := qs[len(qs)-1-i]
302 q.Quad.To, q.Quad.From = q.Quad.From, q.Quad.To
303 ps = append(ps, q)
304 }
305
306 return ps
307 }
308
309 func (qs StrokeQuads) append(ps StrokeQuads) StrokeQuads {
310 switch {
311 case len(ps) == 0:
312 return qs
313 case len(qs) == 0:
314 return ps
315 }
316
317 // Consolidate quads and smooth out rounding errors.
318 // We need to also check for the strokeTolerance to correctly handle
319 // join/cap points or on-purpose disjoint quads.
320 p0 := qs[len(qs)-1].Quad.To
321 p1 := ps[0].Quad.From
322 if p0 != p1 && lenPt(p0.Sub(p1)) < strokeTolerance {
323 qs = append(qs, StrokeQuad{
324 Quad: QuadSegment{
325 From: p0,
326 Ctrl: p0.Add(p1).Mul(0.5),
327 To: p1,
328 },
329 })
330 }
331 return append(qs, ps...)
332 }
333
334 func (q QuadSegment) Transform(t f32.Affine2D) QuadSegment {
335 q.From = t.Transform(q.From)
336 q.Ctrl = t.Transform(q.Ctrl)
337 q.To = t.Transform(q.To)
338 return q
339 }
340
341 // strokePathNorm returns the normal vector at t.
342 func strokePathNorm(p0, p1, p2 f32.Point, t, d float32) f32.Point {
343 switch t {
344 case 0:
345 n := p1.Sub(p0)
346 if n.X == 0 && n.Y == 0 {
347 return f32.Point{}
348 }
349 n = rot90CW(n)
350 return normPt(n, d)
351 case 1:
352 n := p2.Sub(p1)
353 if n.X == 0 && n.Y == 0 {
354 return f32.Point{}
355 }
356 n = rot90CW(n)
357 return normPt(n, d)
358 }
359 panic("impossible")
360 }
361
362 func rot90CW(p f32.Point) f32.Point { return f32.Pt(+p.Y, -p.X) }
363 func rot90CCW(p f32.Point) f32.Point { return f32.Pt(-p.Y, +p.X) }
364
365 // cosPt returns the cosine of the opening angle between p and q.
366 func cosPt(p, q f32.Point) float32 {
367 np := math.Hypot(float64(p.X), float64(p.Y))
368 nq := math.Hypot(float64(q.X), float64(q.Y))
369 return dotPt(p, q) / float32(np*nq)
370 }
371
372 func normPt(p f32.Point, l float32) f32.Point {
373 d := math.Hypot(float64(p.X), float64(p.Y))
374 l64 := float64(l)
375 if math.Abs(d-l64) < 1e-10 {
376 return f32.Point{}
377 }
378 n := float32(l64 / d)
379 return f32.Point{X: p.X * n, Y: p.Y * n}
380 }
381
382 func lenPt(p f32.Point) float32 {
383 return float32(math.Hypot(float64(p.X), float64(p.Y)))
384 }
385
386 func dotPt(p, q f32.Point) float32 {
387 return p.X*q.X + p.Y*q.Y
388 }
389
390 func perpDot(p, q f32.Point) float32 {
391 return p.X*q.Y - p.Y*q.X
392 }
393
394 // strokePathCurv returns the curvature at t, along the quadratic Bézier
395 // curve defined by the triplet (beg, ctl, end).
396 func strokePathCurv(beg, ctl, end f32.Point, t float32) float32 {
397 var (
398 d1p = quadBezierD1(beg, ctl, end, t)
399 d2p = quadBezierD2(beg, ctl, end, t)
400
401 // Negative when bending right, ie: the curve is CW at this point.
402 a = float64(perpDot(d1p, d2p))
403 )
404
405 // We check early that the segment isn't too line-like and
406 // save a costly call to math.Pow that will be discarded by dividing
407 // with a too small 'a'.
408 if math.Abs(a) < 1e-10 {
409 return float32(math.NaN())
410 }
411 return float32(math.Pow(float64(d1p.X*d1p.X+d1p.Y*d1p.Y), 1.5) / a)
412 }
413
414 // quadBezierSample returns the point on the Bézier curve at t.
415 // B(t) = (1-t)^2 P0 + 2(1-t)t P1 + t^2 P2
416 func quadBezierSample(p0, p1, p2 f32.Point, t float32) f32.Point {
417 t1 := 1 - t
418 c0 := t1 * t1
419 c1 := 2 * t1 * t
420 c2 := t * t
421
422 o := p0.Mul(c0)
423 o = o.Add(p1.Mul(c1))
424 o = o.Add(p2.Mul(c2))
425 return o
426 }
427
428 // quadBezierD1 returns the first derivative of the Bézier curve with respect to t.
429 // B'(t) = 2(1-t)(P1 - P0) + 2t(P2 - P1)
430 func quadBezierD1(p0, p1, p2 f32.Point, t float32) f32.Point {
431 p10 := p1.Sub(p0).Mul(2 * (1 - t))
432 p21 := p2.Sub(p1).Mul(2 * t)
433
434 return p10.Add(p21)
435 }
436
437 // quadBezierD2 returns the second derivative of the Bézier curve with respect to t:
438 // B''(t) = 2(P2 - 2P1 + P0)
439 func quadBezierD2(p0, p1, p2 f32.Point, t float32) f32.Point {
440 p := p2.Sub(p1.Mul(2)).Add(p0)
441 return p.Mul(2)
442 }
443
444 // quadBezierLen returns the length of the Bézier curve.
445 // See:
446 // https://malczak.linuxpl.com/blog/quadratic-bezier-curve-length/
447 func quadBezierLen(p0, p1, p2 f32.Point) float32 {
448 a := p0.Sub(p1.Mul(2)).Add(p2)
449 b := p1.Mul(2).Sub(p0.Mul(2))
450 A := float64(4 * dotPt(a, a))
451 B := float64(4 * dotPt(a, b))
452 C := float64(dotPt(b, b))
453 if f64Eq(A, 0.0) {
454 // p1 is in the middle between p0 and p2,
455 // so it is a straight line from p0 to p2.
456 return lenPt(p2.Sub(p0))
457 }
458
459 Sabc := 2 * math.Sqrt(A+B+C)
460 A2 := math.Sqrt(A)
461 A32 := 2 * A * A2
462 C2 := 2 * math.Sqrt(C)
463 BA := B / A2
464 return float32((A32*Sabc + A2*B*(Sabc-C2) + (4*C*A-B*B)*math.Log((2*A2+BA+Sabc)/(BA+C2))) / (4 * A32))
465 }
466
467 func strokeQuadBezier(state strokeState, d, flatness float32) StrokeQuads {
468 // Gio strokes are only quadratic Bézier curves, w/o any inflection point.
469 // So we just have to flatten them.
470 var qs StrokeQuads
471 return flattenQuadBezier(qs, state.p0, state.ctl, state.p1, d, flatness)
472 }
473
474 // flattenQuadBezier splits a Bézier quadratic curve into linear sub-segments,
475 // themselves also encoded as Bézier (degenerate, flat) quadratic curves.
476 func flattenQuadBezier(qs StrokeQuads, p0, p1, p2 f32.Point, d, flatness float32) StrokeQuads {
477 var (
478 t float32
479 flat64 = float64(flatness)
480 )
481 for t < 1 {
482 s2 := float64((p2.X-p0.X)*(p1.Y-p0.Y) - (p2.Y-p0.Y)*(p1.X-p0.X))
483 den := math.Hypot(float64(p1.X-p0.X), float64(p1.Y-p0.Y))
484 if s2*den == 0.0 {
485 break
486 }
487
488 s2 /= den
489 t = 2.0 * float32(math.Sqrt(flat64/3.0/math.Abs(s2)))
490 if t >= 1.0 {
491 break
492 }
493 var q0, q1, q2 f32.Point
494 q0, q1, q2, p0, p1, p2 = quadBezierSplit(p0, p1, p2, t)
495 qs.addLine(q0, q1, q2, 0, d)
496 }
497 qs.addLine(p0, p1, p2, 1, d)
498 return qs
499 }
500
501 func (qs *StrokeQuads) addLine(p0, ctrl, p1 f32.Point, t, d float32) {
502
503 switch i := len(*qs); i {
504 case 0:
505 p0 = p0.Add(strokePathNorm(p0, ctrl, p1, 0, d))
506 default:
507 // Address possible rounding errors and use previous point.
508 p0 = (*qs)[i-1].Quad.To
509 }
510
511 p1 = p1.Add(strokePathNorm(p0, ctrl, p1, 1, d))
512
513 *qs = append(*qs,
514 StrokeQuad{
515 Quad: QuadSegment{
516 From: p0,
517 Ctrl: p0.Add(p1).Mul(0.5),
518 To: p1,
519 },
520 },
521 )
522 }
523
524 // quadInterp returns the interpolated point at t.
525 func quadInterp(p, q f32.Point, t float32) f32.Point {
526 return f32.Pt(
527 (1-t)*p.X+t*q.X,
528 (1-t)*p.Y+t*q.Y,
529 )
530 }
531
532 // quadBezierSplit returns the pair of triplets (from,ctrl,to) Bézier curve,
533 // split before (resp. after) the provided parametric t value.
534 func quadBezierSplit(p0, p1, p2 f32.Point, t float32) (f32.Point, f32.Point, f32.Point, f32.Point, f32.Point, f32.Point) {
535
536 var (
537 b0 = p0
538 b1 = quadInterp(p0, p1, t)
539 b2 = quadBezierSample(p0, p1, p2, t)
540
541 a0 = b2
542 a1 = quadInterp(p1, p2, t)
543 a2 = p2
544 )
545
546 return b0, b1, b2, a0, a1, a2
547 }
548
549 // strokePathJoin joins the two paths rhs and lhs, according to the provided
550 // stroke operation.
551 func strokePathJoin(stroke StrokeStyle, rhs, lhs *StrokeQuads, hw float32, pivot, n0, n1 f32.Point, r0, r1 float32) {
552 if stroke.Miter > 0 {
553 strokePathMiterJoin(stroke, rhs, lhs, hw, pivot, n0, n1, r0, r1)
554 return
555 }
556 switch stroke.Join {
557 case BevelJoin:
558 strokePathBevelJoin(rhs, lhs, hw, pivot, n0, n1, r0, r1)
559 case RoundJoin:
560 strokePathRoundJoin(rhs, lhs, hw, pivot, n0, n1, r0, r1)
561 default:
562 panic("impossible")
563 }
564 }
565
566 func strokePathBevelJoin(rhs, lhs *StrokeQuads, hw float32, pivot, n0, n1 f32.Point, r0, r1 float32) {
567
568 rp := pivot.Add(n1)
569 lp := pivot.Sub(n1)
570
571 rhs.lineTo(rp)
572 lhs.lineTo(lp)
573 }
574
575 func strokePathRoundJoin(rhs, lhs *StrokeQuads, hw float32, pivot, n0, n1 f32.Point, r0, r1 float32) {
576 rp := pivot.Add(n1)
577 lp := pivot.Sub(n1)
578 cw := dotPt(rot90CW(n0), n1) >= 0.0
579 switch {
580 case cw:
581 // Path bends to the right, ie. CW (or 180 degree turn).
582 c := pivot.Sub(lhs.pen())
583 angle := -math.Acos(float64(cosPt(n0, n1)))
584 lhs.arc(c, c, float32(angle))
585 lhs.lineTo(lp) // Add a line to accommodate for rounding errors.
586 rhs.lineTo(rp)
587 default:
588 // Path bends to the left, ie. CCW.
589 angle := math.Acos(float64(cosPt(n0, n1)))
590 c := pivot.Sub(rhs.pen())
591 rhs.arc(c, c, float32(angle))
592 rhs.lineTo(rp) // Add a line to accommodate for rounding errors.
593 lhs.lineTo(lp)
594 }
595 }
596
597 func strokePathMiterJoin(stroke StrokeStyle, rhs, lhs *StrokeQuads, hw float32, pivot, n0, n1 f32.Point, r0, r1 float32) {
598 if n0 == n1.Mul(-1) {
599 strokePathBevelJoin(rhs, lhs, hw, pivot, n0, n1, r0, r1)
600 return
601 }
602
603 // This is to handle nearly linear joints that would be clipped otherwise.
604 limit := math.Max(float64(stroke.Miter), 1.001)
605
606 cw := dotPt(rot90CW(n0), n1) >= 0.0
607 if cw {
608 // hw is used to calculate |R|.
609 // When running CW, n0 and n1 point the other way,
610 // so the sign of r0 and r1 is negated.
611 hw = -hw
612 }
613 hw64 := float64(hw)
614
615 cos := math.Sqrt(0.5 * (1 + float64(cosPt(n0, n1))))
616 d := hw64 / cos
617 if math.Abs(limit*hw64) < math.Abs(d) {
618 stroke.Miter = 0 // Set miter to zero to disable the miter joint.
619 strokePathJoin(stroke, rhs, lhs, hw, pivot, n0, n1, r0, r1)
620 return
621 }
622 mid := pivot.Add(normPt(n0.Add(n1), float32(d)))
623
624 rp := pivot.Add(n1)
625 lp := pivot.Sub(n1)
626 switch {
627 case cw:
628 // Path bends to the right, ie. CW.
629 lhs.lineTo(mid)
630 default:
631 // Path bends to the left, ie. CCW.
632 rhs.lineTo(mid)
633 }
634 rhs.lineTo(rp)
635 lhs.lineTo(lp)
636 }
637
638 // strokePathCap caps the provided path qs, according to the provided stroke operation.
639 func strokePathCap(stroke StrokeStyle, qs *StrokeQuads, hw float32, pivot, n0 f32.Point) {
640 switch stroke.Cap {
641 case FlatCap:
642 strokePathFlatCap(qs, hw, pivot, n0)
643 case SquareCap:
644 strokePathSquareCap(qs, hw, pivot, n0)
645 case RoundCap:
646 strokePathRoundCap(qs, hw, pivot, n0)
647 default:
648 panic("impossible")
649 }
650 }
651
652 // strokePathFlatCap caps the start or end of a path with a flat cap.
653 func strokePathFlatCap(qs *StrokeQuads, hw float32, pivot, n0 f32.Point) {
654 end := pivot.Sub(n0)
655 qs.lineTo(end)
656 }
657
658 // strokePathSquareCap caps the start or end of a path with a square cap.
659 func strokePathSquareCap(qs *StrokeQuads, hw float32, pivot, n0 f32.Point) {
660 var (
661 e = pivot.Add(rot90CCW(n0))
662 corner1 = e.Add(n0)
663 corner2 = e.Sub(n0)
664 end = pivot.Sub(n0)
665 )
666
667 qs.lineTo(corner1)
668 qs.lineTo(corner2)
669 qs.lineTo(end)
670 }
671
672 // strokePathRoundCap caps the start or end of a path with a round cap.
673 func strokePathRoundCap(qs *StrokeQuads, hw float32, pivot, n0 f32.Point) {
674 c := pivot.Sub(qs.pen())
675 qs.arc(c, c, math.Pi)
676 }
677
678 // ArcTransform computes a transformation that can be used for generating quadratic bézier
679 // curve approximations for an arc.
680 //
681 // The math is extracted from the following paper:
682 // "Drawing an elliptical arc using polylines, quadratic or
683 // cubic Bezier curves", L. Maisonobe
684 // An electronic version may be found at:
685 // http://spaceroots.org/documents/ellipse/elliptical-arc.pdf
686 func ArcTransform(p, f1, f2 f32.Point, angle float32, segments int) f32.Affine2D {
687 c := f32.Point{
688 X: 0.5 * (f1.X + f2.X),
689 Y: 0.5 * (f1.Y + f2.Y),
690 }
691
692 // semi-major axis: 2a = |PF1| + |PF2|
693 a := 0.5 * (dist(f1, p) + dist(f2, p))
694
695 // semi-minor axis: c^2 = a^2+b^2 (c: focal distance)
696 f := dist(f1, c)
697 b := math.Sqrt(a*a - f*f)
698
699 var rx, ry, alpha, start float64
700 switch {
701 case a > b:
702 rx = a
703 ry = b
704 default:
705 rx = b
706 ry = a
707 }
708
709 var x float64
710 switch {
711 case f1 == c || f2 == c:
712 // degenerate case of a circle.
713 alpha = 0
714 default:
715 switch {
716 case f1.X > c.X:
717 x = float64(f1.X - c.X)
718 alpha = math.Acos(x / f)
719 case f1.X < c.X:
720 x = float64(f2.X - c.X)
721 alpha = math.Acos(x / f)
722 case f1.X == c.X:
723 // special case of a "vertical" ellipse.
724 alpha = math.Pi / 2
725 if f1.Y < c.Y {
726 alpha = -alpha
727 }
728 }
729 }
730
731 start = math.Acos(float64(p.X-c.X) / dist(c, p))
732 if c.Y > p.Y {
733 start = -start
734 }
735 start -= alpha
736
737 var (
738 θ = angle / float32(segments)
739 ref f32.Affine2D // transform from absolute frame to ellipse-based one
740 rot f32.Affine2D // rotation matrix for each segment
741 inv f32.Affine2D // transform from ellipse-based frame to absolute one
742 )
743 ref = ref.Offset(f32.Point{}.Sub(c))
744 ref = ref.Rotate(f32.Point{}, float32(-alpha))
745 ref = ref.Scale(f32.Point{}, f32.Point{
746 X: float32(1 / rx),
747 Y: float32(1 / ry),
748 })
749 inv = ref.Invert()
750 rot = rot.Rotate(f32.Point{}, float32(0.5*θ))
751
752 // Instead of invoking math.Sincos for every segment, compute a rotation
753 // matrix once and apply for each segment.
754 // Before applying the rotation matrix rot, transform the coordinates
755 // to a frame centered to the ellipse (and warped into a unit circle), then rotate.
756 // Finally, transform back into the original frame.
757 return inv.Mul(rot).Mul(ref)
758 }
759
760 func dist(p1, p2 f32.Point) float64 {
761 var (
762 x1 = float64(p1.X)
763 y1 = float64(p1.Y)
764 x2 = float64(p2.X)
765 y2 = float64(p2.Y)
766 dx = x2 - x1
767 dy = y2 - y1
768 )
769 return math.Hypot(dx, dy)
770 }
771
772 func StrokePathCommands(style StrokeStyle, dashes DashOp, scene []byte) StrokeQuads {
773 quads := decodeToStrokeQuads(scene)
774 return quads.stroke(style, dashes)
775 }
776
777 // decodeToStrokeQuads decodes scene commands to quads ready to stroke.
778 func decodeToStrokeQuads(pathData []byte) StrokeQuads {
779 quads := make(StrokeQuads, 0, 2*len(pathData)/(scene.CommandSize+4))
780 for len(pathData) >= scene.CommandSize+4 {
781 contour := binary.LittleEndian.Uint32(pathData)
782 cmd := ops.DecodeCommand(pathData[4:])
783 switch cmd.Op() {
784 case scene.OpLine:
785 var q QuadSegment
786 q.From, q.To = scene.DecodeLine(cmd)
787 q.Ctrl = q.From.Add(q.To).Mul(.5)
788 quad := StrokeQuad{
789 Contour: contour,
790 Quad: q,
791 }
792 quads = append(quads, quad)
793 case scene.OpQuad:
794 var q QuadSegment
795 q.From, q.Ctrl, q.To = scene.DecodeQuad(cmd)
796 quad := StrokeQuad{
797 Contour: contour,
798 Quad: q,
799 }
800 quads = append(quads, quad)
801 case scene.OpCubic:
802 for _, q := range SplitCubic(scene.DecodeCubic(cmd)) {
803 quad := StrokeQuad{
804 Contour: contour,
805 Quad: q,
806 }
807 quads = append(quads, quad)
808 }
809 default:
810 panic("unsupported scene command")
811 }
812 pathData = pathData[scene.CommandSize+4:]
813 }
814 return quads
815 }
816
817 func SplitCubic(from, ctrl0, ctrl1, to f32.Point) []QuadSegment {
818 quads := make([]QuadSegment, 0, 10)
819 // Set the maximum distance proportionally to the longest side
820 // of the bounding rectangle.
821 hull := f32.Rectangle{
822 Min: from,
823 Max: ctrl0,
824 }.Canon().Add(ctrl1).Add(to)
825 l := hull.Dx()
826 if h := hull.Dy(); h > l {
827 l = h
828 }
829 approxCubeTo(&quads, 0, l*0.001, from, ctrl0, ctrl1, to)
830 return quads
831 }
832
833 // approxCubeTo approximates a cubic Bézier by a series of quadratic
834 // curves.
835 func approxCubeTo(quads *[]QuadSegment, splits int, maxDist float32, from, ctrl0, ctrl1, to f32.Point) int {
836 // The idea is from
837 // https://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
838 // where a quadratic approximates a cubic by eliminating its t³ term
839 // from its polynomial expression anchored at the starting point:
840 //
841 // P(t) = pen + 3t(ctrl0 - pen) + 3t²(ctrl1 - 2ctrl0 + pen) + t³(to - 3ctrl1 + 3ctrl0 - pen)
842 //
843 // The control point for the new quadratic Q1 that shares starting point, pen, with P is
844 //
845 // C1 = (3ctrl0 - pen)/2
846 //
847 // The reverse cubic anchored at the end point has the polynomial
848 //
849 // P'(t) = to + 3t(ctrl1 - to) + 3t²(ctrl0 - 2ctrl1 + to) + t³(pen - 3ctrl0 + 3ctrl1 - to)
850 //
851 // The corresponding quadratic Q2 that shares the end point, to, with P has control
852 // point
853 //
854 // C2 = (3ctrl1 - to)/2
855 //
856 // The combined quadratic Bézier, Q, shares both start and end points with its cubic
857 // and use the midpoint between the two curves Q1 and Q2 as control point:
858 //
859 // C = (3ctrl0 - pen + 3ctrl1 - to)/4
860 c := ctrl0.Mul(3).Sub(from).Add(ctrl1.Mul(3)).Sub(to).Mul(1.0 / 4.0)
861 const maxSplits = 32
862 if splits >= maxSplits {
863 *quads = append(*quads, QuadSegment{From: from, Ctrl: c, To: to})
864 return splits
865 }
866 // The maximum distance between the cubic P and its approximation Q given t
867 // can be shown to be
868 //
869 // d = sqrt(3)/36*|to - 3ctrl1 + 3ctrl0 - pen|
870 //
871 // To save a square root, compare d² with the squared tolerance.
872 v := to.Sub(ctrl1.Mul(3)).Add(ctrl0.Mul(3)).Sub(from)
873 d2 := (v.X*v.X + v.Y*v.Y) * 3 / (36 * 36)
874 if d2 <= maxDist*maxDist {
875 *quads = append(*quads, QuadSegment{From: from, Ctrl: c, To: to})
876 return splits
877 }
878 // De Casteljau split the curve and approximate the halves.
879 t := float32(0.5)
880 c0 := from.Add(ctrl0.Sub(from).Mul(t))
881 c1 := ctrl0.Add(ctrl1.Sub(ctrl0).Mul(t))
882 c2 := ctrl1.Add(to.Sub(ctrl1).Mul(t))
883 c01 := c0.Add(c1.Sub(c0).Mul(t))
884 c12 := c1.Add(c2.Sub(c1).Mul(t))
885 c0112 := c01.Add(c12.Sub(c01).Mul(t))
886 splits++
887 splits = approxCubeTo(quads, splits, maxDist, from, c0, c01, c0112)
888 splits = approxCubeTo(quads, splits, maxDist, c0112, c12, c2, to)
889 return splits
890 }
891