freetype/raster: Implement round joins.

R=r, rsc
CC=golang-dev, rog
http://codereview.appspot.com/1746043

example/round/main.go

+// Copyright 2010 The Freetype-Go Authors. All rights reserved.
+// Use of this source code is governed by your choice of either the
+// FreeType License or the GNU General Public License version 2,
+// both of which can be found in the LICENSE file.
+
+// This program visualizes the quadratic approximation to the circle, used to
+// implement round joins when stroking paths. The approximation is used in the
+// stroking code for arcs between 0 and 45 degrees, but is visualized here
+// between 0 and 90 degrees. The discrepancy between the approximation and the
+// true circle is clearly visible at angles above 65 degrees.
+package main
+
+import (
+	"bufio"
+	"fmt"
+	"image"
+	"image/png"
+	"log"
+	"math"
+	"os"
+
+	"freetype-go.googlecode.com/hg/freetype/raster"
+)
+
+func main() {
+	const (
+		n = 17
+		r = 256 * 80
+	)
+	s := raster.Fixed(r * math.Sqrt(2) / 2)
+	t := raster.Fixed(r * math.Tan(math.Pi/8))
+
+	m := image.NewRGBA(800, 600)
+	for y := 0; y < m.Height(); y++ {
+		for x := 0; x < m.Width(); x++ {
+			m.Pixel[y][x] = image.RGBAColor{63, 63, 63, 255}
+		}
+	}
+	mp := raster.NewRGBAPainter(m)
+	mp.SetColor(image.Black)
+	z := raster.NewRasterizer(800, 600)
+
+	for i := 0; i < n; i++ {
+		cx := raster.Fixed(25600 + 51200*(i%4))
+		cy := raster.Fixed(2560 + 32000*(i/4))
+		c := raster.Point{cx, cy}
+		theta := math.Pi * (0.5 + 0.5*float64(i)/(n-1))
+		dx := raster.Fixed(r * math.Cos(theta))
+		dy := raster.Fixed(r * math.Sin(theta))
+		d := raster.Point{dx, dy}
+		// Draw a quarter-circle approximated by two quadratic segments,
+		// with each segment spanning 45 degrees.
+		z.Start(c)
+		z.Add1(c.Add(raster.Point{r, 0}))
+		z.Add2(c.Add(raster.Point{r, t}), c.Add(raster.Point{s, s}))
+		z.Add2(c.Add(raster.Point{t, r}), c.Add(raster.Point{0, r}))
+		// Add another quadratic segment whose angle ranges between 0 and 90 degrees.
+		// For an explanation of the magic constants 22, 150, 181 and 256, read the
+		// comments in the freetype/raster package.
+		dot := 256 * d.Dot(raster.Point{0, r}) / (r * r)
+		multiple := raster.Fixed(150 - 22*(dot-181)/(256-181))
+		z.Add2(c.Add(raster.Point{dx, r + dy}.Mul(multiple)), c.Add(d))
+		// Close the curve.
+		z.Add1(c)
+	}
+	z.Rasterize(mp)
+
+	for i := 0; i < n; i++ {
+		cx := raster.Fixed(25600 + 51200*(i%4))
+		cy := raster.Fixed(2560 + 32000*(i/4))
+		for j := 0; j < n; j++ {
+			theta := math.Pi * float64(j) / (n - 1)
+			dx := raster.Fixed(r * math.Cos(theta))
+			dy := raster.Fixed(r * math.Sin(theta))
+			m.Set(int((cx+dx)/256), int((cy+dy)/256), image.Yellow)
+		}
+	}
+
+	// Save that RGBA image to disk.
+	f, err := os.Open("out.png", os.O_CREAT|os.O_WRONLY, 0600)
+	if err != nil {
+		log.Stderr(err)
+		os.Exit(1)
+	}
+	defer f.Close()
+	b := bufio.NewWriter(f)
+	err = png.Encode(b, m)
+	if err != nil {
+		log.Stderr(err)
+		os.Exit(1)
+	}
+	err = b.Flush()
+	if err != nil {
+		log.Stderr(err)
+		os.Exit(1)
+	}
+	fmt.Println("Wrote out.png OK.")
+}

freetype/raster/geom.go

@@ -13,6 +13,9 @@ import (
 // A Fixed is a 24.8 fixed point number.
 type Fixed int32
 
+// A Fixed64 is a 48.16 fixed point number.
+type Fixed64 int64
+
 // String returns a human-readable representation of a 24.8 fixed point number.
 // For example, the number one-and-a-quarter becomes "1:064".
 func (x Fixed) String() string {
@@ -23,6 +26,16 @@ func (x Fixed) String() string {
 	return fmt.Sprintf("%d:%03d", int32(i), int32(f))
 }
 
+// String returns a human-readable representation of a 48.16 fixed point number.
+// For example, the number one-and-a-quarter becomes "1:00064".
+func (x Fixed64) String() string {
+	i, f := x/65536, x%65536
+	if f < 0 {
+		f = -f
+	}
+	return fmt.Sprintf("%d:%05d", int64(i), int64(f))
+}
+
 // maxAbs returns the maximum of abs(a) and abs(b).
 func maxAbs(a, b Fixed) Fixed {
 	if a < 0 {
@@ -58,6 +71,18 @@ func (p Point) Mul(k Fixed) Point {
 	return Point{p.X * k / 256, p.Y * k / 256}
 }
 
+// Neg returns the vector -p, or equivalently p rotated by 180 degrees.
+func (p Point) Neg() Point {
+	return Point{-p.X, -p.Y}
+}
+
+// Dot returns the dot product p·q.
+func (p Point) Dot(q Point) Fixed64 {
+	px, py := int64(p.X), int64(p.Y)
+	qx, qy := int64(q.X), int64(q.Y)
+	return Fixed64(px*qx + py*qy)
+}
+
 // Len returns the length of the vector p.
 func (p Point) Len() Fixed {
 	// TODO(nigeltao): use fixed point math.
@@ -73,22 +98,64 @@ func (p Point) Norm(length Fixed) Point {
 	if d == 0 {
 		return Point{0, 0}
 	}
-	// TODO(nigeltao): should we check for overflow?
-	return Point{p.X * length / d, p.Y * length / d}
+	s, t := int64(length), int64(d)
+	x := int64(p.X) * s / t
+	y := int64(p.Y) * s / t
+	return Point{Fixed(x), Fixed(y)}
 }
 
-// RotateCW returns the vector p rotated clockwise by 90 degrees.
-// Note that the Y-axis grows downwards, so {1, 0}.RotateCW is {0, 1}.
-func (p Point) RotateCW() Point {
+// Rot45CW returns the vector p rotated clockwise by 45 degrees.
+// Note that the Y-axis grows downwards, so {1, 0}.Rot45CW is {1/√2, 1/√2}.
+func (p Point) Rot45CW() Point {
+	// 181/256 is approximately 1/√2, or sin(π/4).
+	px, py := int64(p.X), int64(p.Y)
+	qx := (+px - py) * 181 / 256
+	qy := (+px + py) * 181 / 256
+	return Point{Fixed(qx), Fixed(qy)}
+}
+
+// Rot90CW returns the vector p rotated clockwise by 90 degrees.
+// Note that the Y-axis grows downwards, so {1, 0}.Rot90CW is {0, 1}.
+func (p Point) Rot90CW() Point {
 	return Point{-p.Y, p.X}
 }
 
-// RotateCCW returns the vector p rotated counter-clockwise by 90 degrees.
-// Note that the Y-axis grows downwards, so {1, 0}.RotateCCW is {0, -1}.
-func (p Point) RotateCCW() Point {
+// Rot135CW returns the vector p rotated clockwise by 135 degrees.
+// Note that the Y-axis grows downwards, so {1, 0}.Rot135CW is {-1/√2, 1/√2}.
+func (p Point) Rot135CW() Point {
+	// 181/256 is approximately 1/√2, or sin(π/4).
+	px, py := int64(p.X), int64(p.Y)
+	qx := (-px - py) * 181 / 256
+	qy := (+px - py) * 181 / 256
+	return Point{Fixed(qx), Fixed(qy)}
+}
+
+// Rot45CCW returns the vector p rotated counter-clockwise by 45 degrees.
+// Note that the Y-axis grows downwards, so {1, 0}.Rot45CCW is {1/√2, -1/√2}.
+func (p Point) Rot45CCW() Point {
+	// 181/256 is approximately 1/√2, or sin(π/4).
+	px, py := int64(p.X), int64(p.Y)
+	qx := (+px + py) * 181 / 256
+	qy := (-px + py) * 181 / 256
+	return Point{Fixed(qx), Fixed(qy)}
+}
+
+// Rot90CCW returns the vector p rotated counter-clockwise by 90 degrees.
+// Note that the Y-axis grows downwards, so {1, 0}.Rot90CCW is {0, -1}.
+func (p Point) Rot90CCW() Point {
 	return Point{p.Y, -p.X}
 }
 
+// Rot135CCW returns the vector p rotated counter-clockwise by 135 degrees.
+// Note that the Y-axis grows downwards, so {1, 0}.Rot135CCW is {-1/√2, -1/√2}.
+func (p Point) Rot135CCW() Point {
+	// 181/256 is approximately 1/√2, or sin(π/4).
+	px, py := int64(p.X), int64(p.Y)
+	qx := (-px + py) * 181 / 256
+	qy := (-px - py) * 181 / 256
+	return Point{Fixed(qx), Fixed(qy)}
+}
+
 // An Adder accumulates points on a curve.
 type Adder interface {
 	// Start starts a new curve at the given point.
@@ -258,10 +325,10 @@ func addCap(p Adder, cap Cap, center, end Point) {
 	switch cap {
 	case RoundCap:
 		// The cubic Bézier approximation to a circle involves the magic number
-		// (sqrt(2) - 1) * 4/3, which is approximately 141 / 256.
+		// (√2 - 1) * 4/3, which is approximately 141/256.
 		const k = 141
 		d := end.Sub(center)
-		e := d.RotateCCW()
+		e := d.Rot90CCW()
 		side := center.Add(e)
 		start := center.Sub(d)
 		d, e = d.Mul(k), e.Mul(k)
@@ -271,7 +338,7 @@ func addCap(p Adder, cap Cap, center, end Point) {
 		p.Add1(end)
 	case SquareCap:
 		d := end.Sub(center)
-		e := d.RotateCCW()
+		e := d.Rot90CCW()
 		side := center.Add(e)
 		p.Add1(side.Sub(d))
 		p.Add1(side.Add(d))
@@ -279,6 +346,110 @@ func addCap(p Adder, cap Cap, center, end Point) {
 	}
 }
 
+func addJoin(lhs, rhs Adder, join Join, a, anorm, bnorm Point) {
+	switch join {
+	case RoundJoin:
+		dot := anorm.Rot90CW().Dot(bnorm)
+		if dot >= 0 {
+			addArc(lhs, a, anorm, bnorm)
+			rhs.Add1(a.Sub(bnorm))
+		} else {
+			lhs.Add1(a.Add(bnorm))
+			addArc(rhs, a, anorm.Neg(), bnorm.Neg())
+		}
+	case BevelJoin:
+		lhs.Add1(a.Add(bnorm))
+		rhs.Add1(a.Sub(bnorm))
+	case MiterJoin:
+		panic("freetype/raster: miter join unimplemented")
+	}
+}
+
+// addArc adds a circular arc from pivot+n0 to pivot+n1 to p. The shorter of
+// the two possible arcs is taken, i.e. the one spanning <= 180 degrees.
+// The two vectors n0 and n1 must be of equal length.
+func addArc(p Adder, pivot, n0, n1 Point) {
+	// r2 is the square of the length of n0.
+	r2 := n0.Dot(n0)
+	if r2 < 4096 {
+		// The arc radius is so small that we collapse to a straight line.
+		p.Add1(pivot.Add(n1))
+		return
+	}
+	// We approximate the arc by 0, 1, 2 or 3 45-degree quadratic segments plus
+	// a final quadratic segment from s to n1. Each 45-degree segment has control
+	// points {1, 0}, {1, tan(π/8)} and {1/√2, 1/√2} suitably scaled, rotated and
+	// translated. tan(π/8) is approximately 106/256.
+	const t = 106
+	var s Point
+	// We determine which octant the angle between n0 and n1 is in via three dot products.
+	// m0, m1 and m2 are n0 rotated clockwise by 45, 90 and 135 degrees.
+	m0 := n0.Rot45CW()
+	m1 := n0.Rot90CW()
+	m2 := m0.Rot90CW()
+	if m1.Dot(n1) >= 0 {
+		if n0.Dot(n1) >= 0 {
+			if m2.Dot(n1) <= 0 {
+				// n1 is between 0 and 45 degrees clockwise of n0.
+				s = n0
+			} else {
+				// n1 is between 45 and 90 degrees clockwise of n0.
+				p.Add2(pivot.Add(n0).Add(m1.Mul(t)), pivot.Add(m0))
+				s = m0
+			}
+		} else {
+			pm1, n0t := pivot.Add(m1), n0.Mul(t)
+			p.Add2(pivot.Add(n0).Add(m1.Mul(t)), pivot.Add(m0))
+			p.Add2(pm1.Add(n0t), pm1)
+			if m0.Dot(n1) >= 0 {
+				// n1 is between 90 and 135 degrees clockwise of n0.
+				s = m1
+			} else {
+				// n1 is between 135 and 180 degrees clockwise of n0.
+				p.Add2(pm1.Sub(n0t), pivot.Add(m2))
+				s = m2
+			}
+		}
+	} else {
+		if n0.Dot(n1) >= 0 {
+			if m0.Dot(n1) >= 0 {
+				// n1 is between 0 and 45 degrees counter-clockwise of n0.
+				s = n0
+			} else {
+				// n1 is between 45 and 90 degrees counter-clockwise of n0.
+				p.Add2(pivot.Add(n0).Sub(m1.Mul(t)), pivot.Sub(m2))
+				s = m2.Neg()
+			}
+		} else {
+			pm1, n0t := pivot.Sub(m1), n0.Mul(t)
+			p.Add2(pivot.Add(n0).Sub(m1.Mul(t)), pivot.Sub(m2))
+			p.Add2(pm1.Add(n0t), pm1)
+			if m2.Dot(n1) <= 0 {
+				// n1 is between 90 and 135 degrees counter-clockwise of n0.
+				s = m1.Neg()
+			} else {
+				// n1 is between 135 and 180 degrees counter-clockwise of n0.
+				p.Add2(pm1.Sub(n0t), pivot.Sub(m0))
+				s = m0.Neg()
+			}
+		}
+	}
+	// The final quadratic segment has two endpoints s and n1 and the middle
+	// control point is a multiple of s.Add(n1), i.e. it is on the angle bisector
+	// of those two points. The multiple ranges between 128/256 and 150/256 as
+	// the angle between s and n1 ranges between 0 and 45 degrees.
+	// When the angle is 0 degrees (i.e. s and n1 are coincident) then s.Add(n1)
+	// is twice s and so the middle control point of the degenerate quadratic
+	// segment should be half s.Add(n1), and half = 128/256.
+	// When the angle is 45 degrees then 150/256 is the ratio of the lengths of
+	// the two vectors {1, tan(π/8)} and {1 + 1/√2, 1/√2}.
+	// d is the normalized dot product between s and n1. Since the angle ranges
+	// between 0 and 45 degrees then d ranges between 256/256 and 181/256.
+	d := 256 * s.Dot(n1) / r2
+	multiple := Fixed(150 - 22*(d-181)/(256-181))
+	p.Add2(pivot.Add(s.Add(n1).Mul(multiple)), pivot.Add(n1))
+}
+
 // stroke adds the stroked Path q to p, where q consists of exactly one curve.
 func stroke(p Adder, q Path, width Fixed, cap Cap, join Join) {
 	// Stroking is implemented by deriving two paths each width/2 apart from q.
@@ -286,27 +457,24 @@ func stroke(p Adder, q Path, width Fixed, cap Cap, join Join) {
 	// path is accumulated in r, and once we've finished adding the LHS to p
 	// we add the RHS in reverse order.
 	r := Path(make([]Fixed, 0, len(q)))
-	var start Point
+	var start, anorm Point
 	a := Point{q[1], q[2]}
 	i := 4
 	for i < len(q) {
 		switch q[i] {
 		case 1:
-			bx, by := q[i+1], q[i+2]
-			delta := Point{bx - a.X, by - a.Y}
-			normal := delta.Norm(width / 2).RotateCCW()
+			b := Point{q[i+1], q[i+2]}
+			bnorm := b.Sub(a).Norm(width / 2).Rot90CCW()
 			if i == 4 {
-				start = Point{a.X + normal.X, a.Y + normal.Y}
+				start = a.Add(bnorm)
 				p.Start(start)
-				r.Start(Point{a.X - normal.X, a.Y - normal.Y})
+				r.Start(a.Sub(bnorm))
 			} else {
-				// TODO(nigeltao): handle joins.
-				p.Add1(Point{a.X + normal.X, a.Y + normal.Y})
-				r.Add1(Point{a.X - normal.X, a.Y - normal.Y})
+				addJoin(p, &r, join, a, anorm, bnorm)
 			}
-			p.Add1(Point{bx + normal.X, by + normal.Y})
-			r.Add1(Point{bx - normal.X, by - normal.Y})
-			a = Point{q[i+1], q[i+2]}
+			p.Add1(b.Add(bnorm))
+			r.Add1(b.Sub(bnorm))
+			a, anorm = b, bnorm
 			i += 4
 		case 2:
 			panic("freetype/raster: stroke unimplemented for quadratic segments")