stroke.mx raw

   1  // SPDX-License-Identifier: Unlicense OR MIT
   2  
   3  // Stroke-to-fill conversion. Quadratic Bezier offset curves, round joins/caps.
   4  // Ported from gioui.org/internal/stroke.
   5  //
   6  // Algorithms from:
   7  //   Fast, precise flattening of cubic Bezier path and offset curves
   8  //   Thomas F. Hain, et al.
   9  
  10  package gio
  11  
  12  import "math"
  13  
  14  type StrokeStyle struct {
  15  	Width float32
  16  }
  17  
  18  // strokeTolerance reconciles rounding errors when splitting/joining quads.
  19  const strokeTolerance = 0.01
  20  
  21  type QuadSegment struct {
  22  	From, Ctrl, To Point
  23  }
  24  
  25  type StrokeQuad struct {
  26  	Contour uint32
  27  	Quad    QuadSegment
  28  }
  29  
  30  type strokeState struct {
  31  	p0, p1 Point   // start, end
  32  	n0, n1 Point   // normals at start, end
  33  	r0, r1 float32 // curvature at start, end
  34  	ctl    Point   // control point
  35  }
  36  
  37  type StrokeQuads []StrokeQuad
  38  
  39  func (qs *StrokeQuads) pen() Point {
  40  	return (*qs)[len(*qs)-1].Quad.To
  41  }
  42  
  43  func (qs *StrokeQuads) lineTo(pt Point) {
  44  	end := qs.pen()
  45  	*qs = append(*qs, StrokeQuad{
  46  		Quad: QuadSegment{
  47  			From: end,
  48  			Ctrl: end.Add(pt).Mul(0.5),
  49  			To:   pt,
  50  		},
  51  	})
  52  }
  53  
  54  func (qs *StrokeQuads) arc(f1, f2 Point, angle float32) {
  55  	pen := qs.pen()
  56  	m, segments := ArcTransform(pen, f1.Add(pen), f2.Add(pen), angle)
  57  	for i := 0; i < segments; i++ {
  58  		p0 := qs.pen()
  59  		p1 := m.Transform(p0)
  60  		p2 := m.Transform(p1)
  61  		ctl := p1.Mul(2).Sub(p0.Add(p2).Mul(.5))
  62  		*qs = append(*qs, StrokeQuad{
  63  			Quad: QuadSegment{
  64  				From: p0, Ctrl: ctl, To: p2,
  65  			},
  66  		})
  67  	}
  68  }
  69  
  70  // split splits quads into slices grouped by contour.
  71  func (qs StrokeQuads) split() []StrokeQuads {
  72  	if len(qs) == 0 {
  73  		return nil
  74  	}
  75  
  76  	var (
  77  		c uint32
  78  		o []StrokeQuads
  79  		i = len(o)
  80  	)
  81  	for _, q := range qs {
  82  		if q.Contour != c {
  83  			c = q.Contour
  84  			i = len(o)
  85  			o = append(o, StrokeQuads{})
  86  		}
  87  		o[i] = append(o[i], q)
  88  	}
  89  
  90  	return o
  91  }
  92  
  93  func (qs StrokeQuads) stroke(style StrokeStyle) StrokeQuads {
  94  	var (
  95  		o  StrokeQuads
  96  		hw = 0.5 * style.Width
  97  	)
  98  
  99  	for _, ps := range qs.split() {
 100  		rhs, lhs := ps.offset(hw, style)
 101  		switch lhs {
 102  		case nil:
 103  			o = o.appendQuads(rhs)
 104  		default:
 105  			// Closed path. Inner path goes opposite direction.
 106  			switch {
 107  			case ps.ccw():
 108  				lhs = lhs.reverse()
 109  				o = o.appendQuads(rhs)
 110  				o = o.appendQuads(lhs)
 111  			default:
 112  				rhs = rhs.reverse()
 113  				o = o.appendQuads(lhs)
 114  				o = o.appendQuads(rhs)
 115  			}
 116  		}
 117  	}
 118  
 119  	return o
 120  }
 121  
 122  // offset returns right-hand and left-hand sides of the path, offset by hw.
 123  func (qs StrokeQuads) offset(hw float32, style StrokeStyle) (rhs, lhs StrokeQuads) {
 124  	var (
 125  		states []strokeState
 126  		beg    = qs[0].Quad.From
 127  		end    = qs[len(qs)-1].Quad.To
 128  		closed = beg == end
 129  	)
 130  	for i := range qs {
 131  		q := qs[i].Quad
 132  
 133  		var (
 134  			n0 = strokePathNorm(q.From, q.Ctrl, q.To, 0, hw)
 135  			n1 = strokePathNorm(q.From, q.Ctrl, q.To, 1, hw)
 136  			r0 = strokePathCurv(q.From, q.Ctrl, q.To, 0)
 137  			r1 = strokePathCurv(q.From, q.Ctrl, q.To, 1)
 138  		)
 139  		states = append(states, strokeState{
 140  			p0:  q.From,
 141  			p1:  q.To,
 142  			n0:  n0,
 143  			n1:  n1,
 144  			r0:  r0,
 145  			r1:  r1,
 146  			ctl: q.Ctrl,
 147  		})
 148  	}
 149  
 150  	for i, state := range states {
 151  		rhs = rhs.appendQuads(strokeQuadBezier(state, +hw, strokeTolerance))
 152  		lhs = lhs.appendQuads(strokeQuadBezier(state, -hw, strokeTolerance))
 153  
 154  		if hasNext := i+1 < len(states); hasNext || closed {
 155  			var next strokeState
 156  			switch {
 157  			case hasNext:
 158  				next = states[i+1]
 159  			case closed:
 160  				next = states[0]
 161  			}
 162  			if state.n1 != next.n0 {
 163  				strokePathRoundJoin(&rhs, &lhs, hw, state.p1, state.n1, next.n0, state.r1, next.r0)
 164  			}
 165  		}
 166  	}
 167  
 168  	if closed {
 169  		rhs.close()
 170  		lhs.close()
 171  		return rhs, lhs
 172  	}
 173  
 174  	qbeg := &states[0]
 175  	qend := &states[len(states)-1]
 176  
 177  	lhs = lhs.reverse()
 178  	strokePathCap(style, &rhs, hw, qend.p1, qend.n1)
 179  
 180  	rhs = rhs.appendQuads(lhs)
 181  	strokePathCap(style, &rhs, hw, qbeg.p0, qbeg.n0.Mul(-1))
 182  
 183  	rhs.close()
 184  
 185  	return rhs, nil
 186  }
 187  
 188  func (qs *StrokeQuads) close() {
 189  	p0 := (*qs)[len(*qs)-1].Quad.To
 190  	p1 := (*qs)[0].Quad.From
 191  
 192  	if p1 == p0 {
 193  		return
 194  	}
 195  
 196  	*qs = append(*qs, StrokeQuad{
 197  		Quad: QuadSegment{
 198  			From: p0,
 199  			Ctrl: p0.Add(p1).Mul(0.5),
 200  			To:   p1,
 201  		},
 202  	})
 203  }
 204  
 205  // ccw returns whether the path is counter-clockwise (Shoelace formula).
 206  func (qs StrokeQuads) ccw() bool {
 207  	var area float32
 208  	for _, ps := range qs.split() {
 209  		for i := 1; i < len(ps); i++ {
 210  			pi := ps[i].Quad.To
 211  			pj := ps[i-1].Quad.To
 212  			area += (pi.X - pj.X) * (pi.Y + pj.Y)
 213  		}
 214  	}
 215  	return area <= 0.0
 216  }
 217  
 218  func (qs StrokeQuads) reverse() StrokeQuads {
 219  	if len(qs) == 0 {
 220  		return nil
 221  	}
 222  
 223  	ps := (make)(StrokeQuads, 0, len(qs))
 224  	for i := range qs {
 225  		q := qs[len(qs)-1-i]
 226  		q.Quad.To, q.Quad.From = q.Quad.From, q.Quad.To
 227  		ps = append(ps, q)
 228  	}
 229  
 230  	return ps
 231  }
 232  
 233  // appendQuads joins two quad sequences, smoothing rounding errors at the seam.
 234  // Named appendQuads to avoid collision with builtin append.
 235  func (qs StrokeQuads) appendQuads(ps StrokeQuads) StrokeQuads {
 236  	switch {
 237  	case len(ps) == 0:
 238  		return qs
 239  	case len(qs) == 0:
 240  		return ps
 241  	}
 242  
 243  	p0 := qs[len(qs)-1].Quad.To
 244  	p1 := ps[0].Quad.From
 245  	if p0 != p1 && lenPt(p0.Sub(p1)) < strokeTolerance {
 246  		qs = append(qs, StrokeQuad{
 247  			Quad: QuadSegment{
 248  				From: p0,
 249  				Ctrl: p0.Add(p1).Mul(0.5),
 250  				To:   p1,
 251  			},
 252  		})
 253  	}
 254  	return append(qs, ps...)
 255  }
 256  
 257  func (q QuadSegment) Transform(t Affine2D) QuadSegment {
 258  	q.From = t.Transform(q.From)
 259  	q.Ctrl = t.Transform(q.Ctrl)
 260  	q.To = t.Transform(q.To)
 261  	return q
 262  }
 263  
 264  // strokePathNorm returns the normal vector at t.
 265  func strokePathNorm(p0, p1, p2 Point, t, d float32) Point {
 266  	switch t {
 267  	case 0:
 268  		n := p1.Sub(p0)
 269  		if n.X == 0 && n.Y == 0 {
 270  			return Point{}
 271  		}
 272  		n = rot90CW(n)
 273  		return normPt(n, d)
 274  	case 1:
 275  		n := p2.Sub(p1)
 276  		if n.X == 0 && n.Y == 0 {
 277  			return Point{}
 278  		}
 279  		n = rot90CW(n)
 280  		return normPt(n, d)
 281  	}
 282  	panic("impossible")
 283  }
 284  
 285  func rot90CW(p Point) Point { return Pt(+p.Y, -p.X) }
 286  
 287  func normPt(p Point, l float32) Point {
 288  	if (p.X == 0 && p.Y == 0) || l == 0 {
 289  		return Point{}
 290  	}
 291  	isVerticalUnit := p.X == 0 && (p.Y == l || p.Y == -l)
 292  	isHorizontalUnit := p.Y == 0 && (p.X == l || p.X == -l)
 293  	if isVerticalUnit || isHorizontalUnit {
 294  		if math.Signbit(float64(l)) {
 295  			return Point{X: -p.X, Y: -p.Y}
 296  		} else {
 297  			return Point{X: p.X, Y: p.Y}
 298  		}
 299  	}
 300  	d := math.Hypot(float64(p.X), float64(p.Y))
 301  	l64 := float64(l)
 302  	if math.Abs(d-l64) < 1e-10 {
 303  		if math.Signbit(float64(l)) {
 304  			return Point{X: -p.X, Y: -p.Y}
 305  		} else {
 306  			return Point{X: p.X, Y: p.Y}
 307  		}
 308  	}
 309  	n := float32(l64 / d)
 310  	return Point{X: p.X * n, Y: p.Y * n}
 311  }
 312  
 313  func lenPt(p Point) float32 {
 314  	return float32(math.Hypot(float64(p.X), float64(p.Y)))
 315  }
 316  
 317  func perpDot(p, q Point) float32 {
 318  	return p.X*q.Y - p.Y*q.X
 319  }
 320  
 321  func angleBetween(n0, n1 Point) float64 {
 322  	return math.Atan2(float64(n1.Y), float64(n1.X)) -
 323  		math.Atan2(float64(n0.Y), float64(n0.X))
 324  }
 325  
 326  // strokePathCurv returns the curvature at t along the quadratic Bezier (beg, ctl, end).
 327  func strokePathCurv(beg, ctl, end Point, t float32) float32 {
 328  	var (
 329  		d1p = quadBezierD1(beg, ctl, end, t)
 330  		d2p = quadBezierD2(beg, ctl, end, t)
 331  		a   = float64(perpDot(d1p, d2p))
 332  	)
 333  
 334  	if math.Abs(a) < 1e-10 {
 335  		return float32(math.NaN())
 336  	}
 337  	return float32(math.Pow(float64(d1p.X*d1p.X+d1p.Y*d1p.Y), 1.5) / a)
 338  }
 339  
 340  // quadBezierSample: B(t) = (1-t)^2 P0 + 2(1-t)t P1 + t^2 P2
 341  func quadBezierSample(p0, p1, p2 Point, t float32) Point {
 342  	t1 := 1 - t
 343  	c0 := t1 * t1
 344  	c1 := 2 * t1 * t
 345  	c2 := t * t
 346  
 347  	o := p0.Mul(c0)
 348  	o = o.Add(p1.Mul(c1))
 349  	o = o.Add(p2.Mul(c2))
 350  	return o
 351  }
 352  
 353  // quadBezierD1: B'(t) = 2(1-t)(P1 - P0) + 2t(P2 - P1)
 354  func quadBezierD1(p0, p1, p2 Point, t float32) Point {
 355  	p10 := p1.Sub(p0).Mul(2 * (1 - t))
 356  	p21 := p2.Sub(p1).Mul(2 * t)
 357  	return p10.Add(p21)
 358  }
 359  
 360  // quadBezierD2: B''(t) = 2(P2 - 2P1 + P0)
 361  func quadBezierD2(p0, p1, p2 Point, t float32) Point {
 362  	p := p2.Sub(p1.Mul(2)).Add(p0)
 363  	return p.Mul(2)
 364  }
 365  
 366  func strokeQuadBezier(state strokeState, d, flatness float32) StrokeQuads {
 367  	var qs StrokeQuads
 368  	return flattenQuadBezier(qs, state.p0, state.ctl, state.p1, d, flatness)
 369  }
 370  
 371  // flattenQuadBezier splits a quadratic Bezier into flat sub-segments.
 372  func flattenQuadBezier(qs StrokeQuads, p0, p1, p2 Point, d, flatness float32) StrokeQuads {
 373  	var (
 374  		t      float32
 375  		flat64 = float64(flatness)
 376  	)
 377  	for t < 1 {
 378  		s2 := float64((p2.X-p0.X)*(p1.Y-p0.Y) - (p2.Y-p0.Y)*(p1.X-p0.X))
 379  		den := math.Hypot(float64(p1.X-p0.X), float64(p1.Y-p0.Y))
 380  		if s2*den == 0.0 {
 381  			break
 382  		}
 383  
 384  		s2 /= den
 385  		t = 2.0 * float32(math.Sqrt(flat64/3.0/math.Abs(s2)))
 386  		if t >= 1.0 {
 387  			break
 388  		}
 389  		var q0, q1, q2 Point
 390  		q0, q1, q2, p0, p1, p2 = quadBezierSplit(p0, p1, p2, t)
 391  		qs.addLine(q0, q1, q2, 0, d)
 392  	}
 393  	qs.addLine(p0, p1, p2, 1, d)
 394  	return qs
 395  }
 396  
 397  func (qs *StrokeQuads) addLine(p0, ctrl, p1 Point, t, d float32) {
 398  	switch i := len(*qs); i {
 399  	case 0:
 400  		p0 = p0.Add(strokePathNorm(p0, ctrl, p1, 0, d))
 401  	default:
 402  		p0 = (*qs)[i-1].Quad.To
 403  	}
 404  
 405  	p1 = p1.Add(strokePathNorm(p0, ctrl, p1, 1, d))
 406  
 407  	*qs = append(*qs,
 408  		StrokeQuad{
 409  			Quad: QuadSegment{
 410  				From: p0,
 411  				Ctrl: p0.Add(p1).Mul(0.5),
 412  				To:   p1,
 413  			},
 414  		},
 415  	)
 416  }
 417  
 418  // quadInterp returns the interpolated point at t.
 419  func quadInterp(p, q Point, t float32) Point {
 420  	return Pt(
 421  		(1-t)*p.X+t*q.X,
 422  		(1-t)*p.Y+t*q.Y,
 423  	)
 424  }
 425  
 426  // quadBezierSplit returns (before, after) triplets split at t.
 427  func quadBezierSplit(p0, p1, p2 Point, t float32) (Point, Point, Point, Point, Point, Point) {
 428  	var (
 429  		b0 = p0
 430  		b1 = quadInterp(p0, p1, t)
 431  		b2 = quadBezierSample(p0, p1, p2, t)
 432  
 433  		a0 = b2
 434  		a1 = quadInterp(p1, p2, t)
 435  		a2 = p2
 436  	)
 437  
 438  	return b0, b1, b2, a0, a1, a2
 439  }
 440  
 441  // strokePathRoundJoin creates a round join between rhs and lhs.
 442  func strokePathRoundJoin(rhs, lhs *StrokeQuads, hw float32, pivot, n0, n1 Point, r0, r1 float32) {
 443  	rp := pivot.Add(n1)
 444  	lp := pivot.Sub(n1)
 445  	angle := angleBetween(n0, n1)
 446  	switch {
 447  	case angle <= 0:
 448  		// CW bend or 180 degree turn.
 449  		c := pivot.Sub(lhs.pen())
 450  		lhs.arc(c, c, float32(angle))
 451  		lhs.lineTo(lp)
 452  		rhs.lineTo(rp)
 453  	default:
 454  		// CCW bend.
 455  		c := pivot.Sub(rhs.pen())
 456  		rhs.arc(c, c, float32(angle))
 457  		rhs.lineTo(rp)
 458  		lhs.lineTo(lp)
 459  	}
 460  }
 461  
 462  // strokePathCap caps the path with a round cap.
 463  func strokePathCap(style StrokeStyle, qs *StrokeQuads, hw float32, pivot, n0 Point) {
 464  	strokePathRoundCap(qs, hw, pivot, n0)
 465  }
 466  
 467  func strokePathRoundCap(qs *StrokeQuads, hw float32, pivot, n0 Point) {
 468  	c := pivot.Sub(qs.pen())
 469  	qs.arc(c, c, math.Pi)
 470  }
 471  
 472  // ArcTransform computes a transform for generating quadratic Bezier arc approximations.
 473  //
 474  // Math from: "Drawing an elliptical arc using polylines, quadratic or
 475  // cubic Bezier curves", L. Maisonobe
 476  func ArcTransform(p, f1, f2 Point, angle float32) (transform Affine2D, segments int) {
 477  	const segmentsPerCircle = 16
 478  	const anglePerSegment = 2 * math.Pi / segmentsPerCircle
 479  
 480  	s := angle / anglePerSegment
 481  	if s < 0 {
 482  		s = -s
 483  	}
 484  	segments = int(math.Ceil(float64(s)))
 485  	if segments <= 0 {
 486  		segments = 1
 487  	}
 488  
 489  	var rx, ry, alpha float64
 490  	if f1 == f2 {
 491  		rx = dist(f1, p)
 492  		ry = rx
 493  	} else {
 494  		a := 0.5 * (dist(f1, p) + dist(f2, p))
 495  		c := dist(f1, f2) * 0.5
 496  		b := math.Sqrt(a*a - c*c)
 497  		switch {
 498  		case a > b:
 499  			rx = a
 500  			ry = b
 501  		default:
 502  			rx = b
 503  			ry = a
 504  		}
 505  		if f1.X == f2.X {
 506  			alpha = math.Pi / 2
 507  			if f1.Y < f2.Y {
 508  				alpha = -alpha
 509  			}
 510  		} else {
 511  			x := float64(f1.X-f2.X) * 0.5
 512  			if x < 0 {
 513  				x = -x
 514  			}
 515  			alpha = math.Acos(x / c)
 516  		}
 517  	}
 518  
 519  	th := angle / float32(segments)
 520  	ref := AffineId()
 521  	rot := AffineId()
 522  	inv := AffineId()
 523  	center := Point{
 524  		X: 0.5 * (f1.X + f2.X),
 525  		Y: 0.5 * (f1.Y + f2.Y),
 526  	}
 527  	ref = ref.Offset(Point{}.Sub(center))
 528  	ref = ref.Rotate(Point{}, float32(-alpha))
 529  	ref = ref.Scale(Point{}, Point{
 530  		X: float32(1 / rx),
 531  		Y: float32(1 / ry),
 532  	})
 533  	inv = ref.Invert()
 534  	rot = rot.Rotate(Point{}, 0.5*th)
 535  
 536  	return inv.Mul(rot).Mul(ref), segments
 537  }
 538  
 539  func dist(p1, p2 Point) float64 {
 540  	var (
 541  		x1 = float64(p1.X)
 542  		y1 = float64(p1.Y)
 543  		x2 = float64(p2.X)
 544  		y2 = float64(p2.Y)
 545  		dx = x2 - x1
 546  		dy = y2 - y1
 547  	)
 548  	return math.Hypot(dx, dy)
 549  }
 550  
 551  func StrokePathCommands(style StrokeStyle, scene []byte) StrokeQuads {
 552  	quads := decodeToStrokeQuads(scene)
 553  	return quads.stroke(style)
 554  }
 555  
 556  // decodeToStrokeQuads decodes scene commands to quads ready to stroke.
 557  func decodeToStrokeQuads(pathData []byte) StrokeQuads {
 558  	quads := (make)(StrokeQuads, 0, 2*len(pathData)/(CommandSize+4))
 559  	scratch := []QuadSegment{:0:10}
 560  	for len(pathData) >= CommandSize+4 {
 561  		contour := leUint32(pathData)
 562  		cmd := DecodeCommand(pathData[4:])
 563  		switch cmd.Op() {
 564  		case OpLine:
 565  			var q QuadSegment
 566  			q.From, q.To = DecodeLine(cmd)
 567  			q.Ctrl = q.From.Add(q.To).Mul(.5)
 568  			quad := StrokeQuad{
 569  				Contour: contour,
 570  				Quad:    q,
 571  			}
 572  			quads = append(quads, quad)
 573  		case OpGap:
 574  			// Ignore gaps for strokes.
 575  		case OpQuad:
 576  			var q QuadSegment
 577  			q.From, q.Ctrl, q.To = DecodeQuad(cmd)
 578  			quad := StrokeQuad{
 579  				Contour: contour,
 580  				Quad:    q,
 581  			}
 582  			quads = append(quads, quad)
 583  		case OpCubic:
 584  			from, ctrl0, ctrl1, to := DecodeCubic(cmd)
 585  			scratch = SplitCubic(from, ctrl0, ctrl1, to, scratch[:0])
 586  			for _, q := range scratch {
 587  				quad := StrokeQuad{
 588  					Contour: contour,
 589  					Quad:    q,
 590  				}
 591  				quads = append(quads, quad)
 592  			}
 593  		default:
 594  			panic("unsupported scene command")
 595  		}
 596  		pathData = pathData[CommandSize+4:]
 597  	}
 598  	return quads
 599  }
 600  
 601  // leUint32 reads a little-endian uint32 from b.
 602  func leUint32(b []byte) uint32 {
 603  	return uint32(b[0]) | uint32(b[1])<<8 | uint32(b[2])<<16 | uint32(b[3])<<24
 604  }
 605  
 606  func SplitCubic(from, ctrl0, ctrl1, to Point, quads []QuadSegment) []QuadSegment {
 607  	hull := Rectangle{
 608  		Min: from,
 609  		Max: ctrl0,
 610  	}.Canon().Union(Rectangle{
 611  		Min: ctrl1,
 612  		Max: to,
 613  	}.Canon())
 614  	l := hull.Dx()
 615  	if h := hull.Dy(); h > l {
 616  		l = h
 617  	}
 618  	maxDist := l * 0.001
 619  	approxCubeTo(&quads, 0, maxDist*maxDist, from, ctrl0, ctrl1, to)
 620  	return quads
 621  }
 622  
 623  // approxCubeTo approximates a cubic Bezier by a series of quadratic curves.
 624  func approxCubeTo(quads *[]QuadSegment, splits int, maxDistSq float32, from, ctrl0, ctrl1, to Point) int {
 625  	// Quadratic approximation: eliminate the t^3 term.
 626  	// C = (3*ctrl0 - from + 3*ctrl1 - to) / 4
 627  	q0 := ctrl0.Mul(3).Sub(from)
 628  	q1 := ctrl1.Mul(3).Sub(to)
 629  	c := q0.Add(q1).Mul(1.0 / 4.0)
 630  	const maxSplits = 32
 631  	if splits >= maxSplits {
 632  		*quads = append(*quads, QuadSegment{From: from, Ctrl: c, To: to})
 633  		return splits
 634  	}
 635  	// d = sqrt(3)/36 * |q0 - q1|; compare d^2 with tolerance^2.
 636  	v := q0.Sub(q1)
 637  	d2 := (v.X*v.X + v.Y*v.Y) * 3 / (36 * 36)
 638  	if d2 <= maxDistSq {
 639  		*quads = append(*quads, QuadSegment{From: from, Ctrl: c, To: to})
 640  		return splits
 641  	}
 642  	// De Casteljau split at t=0.5.
 643  	t := float32(0.5)
 644  	c0 := from.Add(ctrl0.Sub(from).Mul(t))
 645  	c1 := ctrl0.Add(ctrl1.Sub(ctrl0).Mul(t))
 646  	c2 := ctrl1.Add(to.Sub(ctrl1).Mul(t))
 647  	c01 := c0.Add(c1.Sub(c0).Mul(t))
 648  	c12 := c1.Add(c2.Sub(c1).Mul(t))
 649  	c0112 := c01.Add(c12.Sub(c01).Mul(t))
 650  	splits++
 651  	splits = approxCubeTo(quads, splits, maxDistSq, from, c0, c01, c0112)
 652  	splits = approxCubeTo(quads, splits, maxDistSq, c0112, c12, c2, to)
 653  	return splits
 654  }
 655