vector: add CubeTo.

Change-Id: Ia64b5b077172303981538de99cc14c5bcc99a3e4
Reviewed-on: https://go-review.googlesource.com/28489
Reviewed-by: David Crawshaw <crawshaw@golang.org>
diff --git a/vector/vector.go b/vector/vector.go
index 3723a84..218e764 100644
--- a/vector/vector.go
+++ b/vector/vector.go
@@ -134,35 +134,67 @@
 //
 // The coordinates are allowed to be out of the Rasterizer's bounds.
 func (z *Rasterizer) QuadTo(b, c f32.Vec2) {
-	// We make a linear approximation to the curve.
-	// http://lists.nongnu.org/archive/html/freetype-devel/2016-08/msg00080.html
-	// gives the rationale for this evenly spaced heuristic instead of a
-	// recursive de Casteljau approach:
-	//
-	// The reason for the subdivision by n is that I expect the "flatness"
-	// computation to be semi-expensive (it's done once rather than on each
-	// potential subdivision) and also because you'll often get fewer
-	// subdivisions. Taking a circular arc as a simplifying assumption (ie a
-	// spherical cow), where I get n, a recursive approach would get 2^⌈lg n⌉,
-	// which, if I haven't made any horrible mistakes, is expected to be 33%
-	// more in the limit.
 	a := z.pen
-	devx := a[0] - 2*b[0] + c[0]
-	devy := a[1] - 2*b[1] + c[1]
-	devsq := devx*devx + devy*devy
+	devsq := devSquared(a, b, c)
 	if devsq >= 0.333 {
 		const tol = 3
 		n := 1 + int(math.Sqrt(math.Sqrt(tol*float64(devsq))))
 		t, nInv := float32(0), 1/float32(n)
 		for i := 0; i < n-1; i++ {
 			t += nInv
-			z.LineTo(lerp(t, lerp(t, a, b), lerp(t, b, c)))
+			ab := lerp(t, a, b)
+			bc := lerp(t, b, c)
+			z.LineTo(lerp(t, ab, bc))
 		}
 	}
 	z.LineTo(c)
 }
 
-// TODO: CubeTo for cubic Béziers.
+// CubeTo adds a cubic Bézier segment, from the pen via b and c to d, and moves
+// the pen to d.
+//
+// The coordinates are allowed to be out of the Rasterizer's bounds.
+func (z *Rasterizer) CubeTo(b, c, d f32.Vec2) {
+	a := z.pen
+	devsq := devSquared(a, b, d)
+	if devsqAlt := devSquared(a, c, d); devsq < devsqAlt {
+		devsq = devsqAlt
+	}
+	if devsq >= 0.333 {
+		const tol = 3
+		n := 1 + int(math.Sqrt(math.Sqrt(tol*float64(devsq))))
+		t, nInv := float32(0), 1/float32(n)
+		for i := 0; i < n-1; i++ {
+			t += nInv
+			ab := lerp(t, a, b)
+			bc := lerp(t, b, c)
+			cd := lerp(t, c, d)
+			abc := lerp(t, ab, bc)
+			bcd := lerp(t, bc, cd)
+			z.LineTo(lerp(t, abc, bcd))
+		}
+	}
+	z.LineTo(d)
+}
+
+// devSquared returns a measure of how curvy the sequnce a to b to c is. It
+// determines how many line segments will approximate a Bézier curve segment.
+//
+// http://lists.nongnu.org/archive/html/freetype-devel/2016-08/msg00080.html
+// gives the rationale for this evenly spaced heuristic instead of a recursive
+// de Casteljau approach:
+//
+// The reason for the subdivision by n is that I expect the "flatness"
+// computation to be semi-expensive (it's done once rather than on each
+// potential subdivision) and also because you'll often get fewer subdivisions.
+// Taking a circular arc as a simplifying assumption (ie a spherical cow),
+// where I get n, a recursive approach would get 2^⌈lg n⌉, which, if I haven't
+// made any horrible mistakes, is expected to be 33% more in the limit.
+func devSquared(a, b, c f32.Vec2) float32 {
+	devx := a[0] - 2*b[0] + c[0]
+	devy := a[1] - 2*b[1] + c[1]
+	return devx*devx + devy*devy
+}
 
 // Draw implements the Drawer interface from the standard library's image/draw
 // package.
diff --git a/vector/vector_test.go b/vector/vector_test.go
index 9b8b642..3aa8c19 100644
--- a/vector/vector_test.go
+++ b/vector/vector_test.go
@@ -9,35 +9,52 @@
 import (
 	"image"
 	"image/draw"
+	"image/png"
+	"os"
 	"testing"
 
 	"golang.org/x/image/math/f32"
 )
 
+// encodePNG is useful for manually debugging the tests.
+func encodePNG(dstFilename string, src image.Image) error {
+	f, err := os.Create(dstFilename)
+	if err != nil {
+		return err
+	}
+	encErr := png.Encode(f, src)
+	closeErr := f.Close()
+	if encErr != nil {
+		return encErr
+	}
+	return closeErr
+}
+
 func TestBasicPath(t *testing.T) {
 	want := []byte{
 		0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
 		0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-		0x00, 0x00, 0xd4, 0xdd, 0xc5, 0xab, 0x63, 0x12, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-		0x00, 0x00, 0x9f, 0xff, 0xff, 0xff, 0xff, 0xfb, 0xb3, 0x20, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
-		0x00, 0x00, 0x60, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xf1, 0x5b, 0x00, 0x00, 0x00, 0x00, 0x00,
-		0x00, 0x00, 0x1f, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x5b, 0x00, 0x00, 0x00, 0x00,
-		0x00, 0x00, 0x00, 0xdf, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xf1, 0x20, 0x00, 0x00, 0x00,
-		0x00, 0x00, 0x00, 0x9f, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xb3, 0x00, 0x00, 0x00,
-		0x00, 0x00, 0x00, 0x5f, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfb, 0x12, 0x00, 0x00,
-		0x00, 0x00, 0x00, 0x1f, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x63, 0x00, 0x00,
-		0x00, 0x00, 0x00, 0x00, 0xdf, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xab, 0x00, 0x00,
-		0x00, 0x00, 0x00, 0x00, 0x9f, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xc5, 0x00, 0x00,
-		0x00, 0x00, 0x00, 0x00, 0x5f, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xdd, 0x00, 0x00,
-		0x00, 0x00, 0x00, 0x00, 0x1f, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xf4, 0x00, 0x00,
+		0x00, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xe3, 0xaa, 0x3e, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0x00, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfa, 0x5f, 0x00, 0x00, 0x00, 0x00,
+		0x00, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfc, 0x24, 0x00, 0x00, 0x00,
+		0x00, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xa1, 0x00, 0x00, 0x00,
+		0x00, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfc, 0x14, 0x00, 0x00,
+		0x00, 0x00, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x4a, 0x00, 0x00,
+		0x00, 0x00, 0xcc, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x81, 0x00, 0x00,
+		0x00, 0x00, 0x66, 0xff, 0xff, 0xff, 0xff, 0xff, 0xef, 0xe4, 0xff, 0xff, 0xff, 0xb6, 0x00, 0x00,
+		0x00, 0x00, 0x0c, 0xf2, 0xff, 0xff, 0xfe, 0x9e, 0x15, 0x00, 0x15, 0x96, 0xff, 0xce, 0x00, 0x00,
+		0x00, 0x00, 0x00, 0x88, 0xfc, 0xe3, 0x43, 0x00, 0x00, 0x00, 0x00, 0x06, 0xcd, 0xdc, 0x00, 0x00,
+		0x00, 0x00, 0x00, 0x00, 0x10, 0x0f, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x25, 0xde, 0x00, 0x00,
+		0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x56, 0x00, 0x00,
 		0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
 		0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
 	}
 
 	z := NewRasterizer(16, 16)
 	z.MoveTo(f32.Vec2{2, 2})
+	z.LineTo(f32.Vec2{8, 2})
 	z.QuadTo(f32.Vec2{14, 2}, f32.Vec2{14, 14})
-	z.LineTo(f32.Vec2{5, 14})
+	z.CubeTo(f32.Vec2{8, 2}, f32.Vec2{5, 20}, f32.Vec2{2, 8})
 	z.ClosePath()
 
 	dst := image.NewAlpha(z.Bounds())