vector: add a fixed point math implementation.

name                      old time/op  new time/op  delta

GlyphAlpha16Over-8        4.48µs ± 1%  3.56µs ± 0%  -20.70%   (p=0.000 n=9+10)
GlyphAlpha16Src-8         4.17µs ± 0%  3.38µs ± 1%  -19.09%  (p=0.000 n=10+10)
GlyphAlpha32Over-8        9.03µs ± 0%  6.74µs ± 0%  -25.33%   (p=0.000 n=9+10)
GlyphAlpha32Src-8         7.46µs ± 1%  5.98µs ± 0%  -19.80%   (p=0.000 n=10+9)
GlyphAlpha64Over-8        21.3µs ± 0%  16.4µs ± 0%  -22.84%  (p=0.000 n=10+10)
GlyphAlpha64Src-8         16.2µs ± 1%  13.1µs ± 0%  -19.33%  (p=0.000 n=10+10)
GlyphAlpha128Over-8       59.8µs ± 0%  47.2µs ± 0%  -21.11%    (p=0.000 n=9+9)
GlyphAlpha128Src-8        41.3µs ± 1%  33.0µs ± 0%  -20.26%   (p=0.000 n=9+10)
GlyphAlpha256Over-8        197µs ± 0%   158µs ± 0%  -19.44%   (p=0.000 n=9+10)
GlyphAlpha256Src-8         124µs ± 0%    98µs ± 0%  -21.17%    (p=0.000 n=9+9)

GlyphAlphaLoose16Over-8   4.73µs ± 0%  3.97µs ± 1%  -16.06%  (p=0.000 n=10+10)
GlyphAlphaLoose16Src-8    4.41µs ± 0%  3.64µs ± 1%  -17.50%  (p=0.000 n=10+10)
GlyphAlphaLoose32Over-8   9.62µs ± 0%  8.47µs ± 0%  -11.95%  (p=0.000 n=10+10)
GlyphAlphaLoose32Src-8    8.25µs ± 0%  7.19µs ± 0%  -12.88%    (p=0.000 n=9+9)
GlyphAlphaLoose64Over-8   25.6µs ± 0%  22.2µs ± 0%  -13.01%    (p=0.000 n=9+9)
GlyphAlphaLoose64Src-8    20.2µs ± 0%  17.2µs ± 1%  -14.98%  (p=0.000 n=10+10)
GlyphAlphaLoose128Over-8  83.4µs ± 1%  68.2µs ± 0%  -18.27%  (p=0.000 n=10+10)
GlyphAlphaLoose128Src-8   59.8µs ± 0%  47.4µs ± 0%  -20.77%   (p=0.000 n=10+9)
GlyphAlphaLoose256Over-8   273µs ± 1%   239µs ± 0%  -12.52%   (p=0.000 n=10+9)
GlyphAlphaLoose256Src-8    187µs ± 0%   155µs ± 1%  -16.91%   (p=0.000 n=9+10)

GlyphRGBA16Over-8         5.99µs ± 0%  5.24µs ± 1%  -12.60%   (p=0.000 n=9+10)
GlyphRGBA16Src-8          5.48µs ± 0%  4.68µs ± 0%  -14.68%   (p=0.000 n=9+10)
GlyphRGBA32Over-8         14.6µs ± 0%  13.5µs ± 0%   -7.60%    (p=0.000 n=9+9)
GlyphRGBA32Src-8          12.6µs ± 0%  11.4µs ± 0%   -9.62%    (p=0.000 n=9+9)
GlyphRGBA64Over-8         44.8µs ± 0%  42.2µs ± 0%   -5.69%    (p=0.000 n=9+9)
GlyphRGBA64Src-8          36.6µs ± 1%  33.5µs ± 1%   -8.55%    (p=0.000 n=9+9)
GlyphRGBA128Over-8         162µs ± 0%   148µs ± 1%   -8.85%   (p=0.000 n=10+9)
GlyphRGBA128Src-8          129µs ± 1%   114µs ± 0%  -11.61%   (p=0.000 n=9+10)
GlyphRGBA256Over-8         588µs ± 0%   573µs ± 0%   -2.53%   (p=0.000 n=9+10)
GlyphRGBA256Src-8          455µs ± 0%   426µs ± 1%   -6.51%   (p=0.000 n=9+10)

GlyphNRGBA16Over-8        27.0µs ± 4%  26.3µs ± 2%   -2.65%   (p=0.001 n=9+10)
GlyphNRGBA16Src-8         19.4µs ± 3%  18.6µs ± 1%   -4.35%   (p=0.000 n=9+10)
GlyphNRGBA32Over-8        97.4µs ± 3%  96.8µs ± 2%     ~      (p=0.447 n=9+10)
GlyphNRGBA32Src-8         66.6µs ± 3%  64.5µs ± 1%   -3.21%   (p=0.000 n=10+9)
GlyphNRGBA64Over-8         372µs ± 3%   368µs ± 1%     ~     (p=0.105 n=10+10)
GlyphNRGBA64Src-8          235µs ± 1%   234µs ± 1%     ~       (p=0.130 n=8+8)
GlyphNRGBA128Over-8       1.45ms ± 2%  1.48ms ± 3%   +2.06%    (p=0.014 n=9+9)
GlyphNRGBA128Src-8         926µs ± 3%   937µs ± 1%     ~      (p=0.113 n=10+9)
GlyphNRGBA256Over-8       5.76ms ± 2%  5.90ms ± 3%   +2.29%   (p=0.001 n=9+10)
GlyphNRGBA256Src-8        3.59ms ± 1%  3.86ms ± 1%   +7.46%   (p=0.000 n=9+10)

Change-Id: I72f25193b5be4e57af09e9eea4eee50545a34cbf
Reviewed-on: https://go-review.googlesource.com/29972
Reviewed-by: David Crawshaw <crawshaw@golang.org>

vector/raster_fixed.go

@@ -0,0 +1,252 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package vector
+
+// This file contains a fixed point math implementation of the vector
+// graphics rasterizer.
+
+import (
+	"golang.org/x/image/math/f32"
+)
+
+const (
+	// ϕ is the number of binary digits after the fixed point.
+	//
+	// For example, if ϕ == 10 (and int1ϕ is based on the int32 type) then we
+	// are using 22.10 fixed point math.
+	//
+	// When changing this number, also change the assembly code (search for ϕ
+	// in the .s files).
+	ϕ = 10
+
+	one          int1ϕ = 1 << ϕ
+	oneAndAHalf  int1ϕ = 1<<ϕ + 1<<(ϕ-1)
+	oneMinusIota int1ϕ = 1<<ϕ - 1 // Used for rounding up.
+)
+
+// int1ϕ is a signed fixed-point number with 1*ϕ binary digits after the fixed
+// point.
+type int1ϕ int32
+
+// int2ϕ is a signed fixed-point number with 2*ϕ binary digits after the fixed
+// point.
+//
+// The Rasterizer's bufU32 field, nominally of type []uint32 (since that slice
+// is also used by other code), can be thought of as a []int2ϕ during the
+// fixedLineTo method. Lines of code that are actually like:
+//	buf[i] += uint32(etc) // buf has type []uint32.
+// can be thought of as
+//	buf[i] += int2ϕ(etc)  // buf has type []int2ϕ.
+type int2ϕ int32
+
+func fixedMax(x, y int1ϕ) int1ϕ {
+	if x > y {
+		return x
+	}
+	return y
+}
+
+func fixedMin(x, y int1ϕ) int1ϕ {
+	if x < y {
+		return x
+	}
+	return y
+}
+
+func fixedFloor(x int1ϕ) int32 { return int32(x >> ϕ) }
+func fixedCeil(x int1ϕ) int32  { return int32((x + oneMinusIota) >> ϕ) }
+
+func (z *Rasterizer) fixedLineTo(b f32.Vec2) {
+	a := z.pen
+	z.pen = b
+	dir := int1ϕ(1)
+	if a[1] > b[1] {
+		dir, a, b = -1, b, a
+	}
+	// Horizontal line segments yield no change in coverage. Almost horizontal
+	// segments would yield some change, in ideal math, but the computation
+	// further below, involving 1 / (b[1] - a[1]), is unstable in fixed point
+	// math, so we treat the segment as if it was perfectly horizontal.
+	if b[1]-a[1] <= 0.000001 {
+		return
+	}
+	dxdy := (b[0] - a[0]) / (b[1] - a[1])
+
+	ay := int1ϕ(a[1] * float32(one))
+	by := int1ϕ(b[1] * float32(one))
+
+	x := int1ϕ(a[0] * float32(one))
+	y := fixedFloor(ay)
+	yMax := fixedCeil(by)
+	if yMax > int32(z.size.Y) {
+		yMax = int32(z.size.Y)
+	}
+	width := int32(z.size.X)
+
+	for ; y < yMax; y++ {
+		dy := fixedMin(int1ϕ(y+1)<<ϕ, by) - fixedMax(int1ϕ(y)<<ϕ, ay)
+		xNext := x + int1ϕ(float32(dy)*dxdy)
+		if y < 0 {
+			x = xNext
+			continue
+		}
+		buf := z.bufU32[y*width:]
+		d := dy * dir
+		x0, x1 := x, xNext
+		if x > xNext {
+			x0, x1 = x1, x0
+		}
+		x0i := fixedFloor(x0)
+		x0Floor := int1ϕ(x0i) << ϕ
+		x1i := fixedCeil(x1)
+		x1Ceil := int1ϕ(x1i) << ϕ
+
+		if x1i <= x0i+1 {
+			xmf := (x+xNext)>>1 - x0Floor
+			if i := clamp(x0i+0, width); i < uint(len(buf)) {
+				buf[i] += uint32(d * (one - xmf))
+			}
+			if i := clamp(x0i+1, width); i < uint(len(buf)) {
+				buf[i] += uint32(d * xmf)
+			}
+		} else {
+			oneOverS := x1 - x0
+			twoOverS := 2 * oneOverS
+			x0f := x0 - x0Floor
+			oneMinusX0f := one - x0f
+			oneMinusX0fSquared := oneMinusX0f * oneMinusX0f
+			x1f := x1 - x1Ceil + one
+			x1fSquared := x1f * x1f
+
+			// These next two variables are unused, as rounding errors are
+			// minimized when we delay the division by oneOverS for as long as
+			// possible. These lines of code (and the "In ideal math" comments
+			// below) are commented out instead of deleted in order to aid the
+			// comparison with the floating point version of the rasterizer.
+			//
+			// a0 := ((oneMinusX0f * oneMinusX0f) >> 1) / oneOverS
+			// am := ((x1f * x1f) >> 1) / oneOverS
+
+			if i := clamp(x0i, width); i < uint(len(buf)) {
+				// In ideal math: buf[i] += uint32(d * a0)
+				D := oneMinusX0fSquared
+				D *= d
+				D /= twoOverS
+				buf[i] += uint32(D)
+			}
+
+			if x1i == x0i+2 {
+				if i := clamp(x0i+1, width); i < uint(len(buf)) {
+					// In ideal math: buf[i] += uint32(d * (one - a0 - am))
+					D := twoOverS<<ϕ - oneMinusX0fSquared - x1fSquared
+					D *= d
+					D /= twoOverS
+					buf[i] += uint32(D)
+				}
+			} else {
+				// This is commented out for the same reason as a0 and am.
+				//
+				// a1 := ((oneAndAHalf - x0f) << ϕ) / oneOverS
+
+				if i := clamp(x0i+1, width); i < uint(len(buf)) {
+					// In ideal math: buf[i] += uint32(d * (a1 - a0))
+					//
+					// Convert to int64 to avoid overflow. Without that,
+					// TestRasterizePolygon fails.
+					D := int64((oneAndAHalf-x0f)<<(ϕ+1) - oneMinusX0fSquared)
+					D *= int64(d)
+					D /= int64(twoOverS)
+					buf[i] += uint32(D)
+				}
+				dTimesS := uint32((d << (2 * ϕ)) / oneOverS)
+				for xi := x0i + 2; xi < x1i-1; xi++ {
+					if i := clamp(xi, width); i < uint(len(buf)) {
+						buf[i] += dTimesS
+					}
+				}
+
+				// This is commented out for the same reason as a0 and am.
+				//
+				// a2 := a1 + (int1ϕ(x1i-x0i-3)<<(2*ϕ))/oneOverS
+
+				if i := clamp(x1i-1, width); i < uint(len(buf)) {
+					// In ideal math: buf[i] += uint32(d * (one - a2 - am))
+					//
+					// Convert to int64 to avoid overflow. Without that,
+					// TestRasterizePolygon fails.
+					D := int64(twoOverS << ϕ)
+					D -= int64((oneAndAHalf - x0f) << (ϕ + 1))
+					D -= int64((x1i - x0i - 3) << (2*ϕ + 1))
+					D -= int64(x1fSquared)
+					D *= int64(d)
+					D /= int64(twoOverS)
+					buf[i] += uint32(D)
+				}
+			}
+
+			if i := clamp(x1i, width); i < uint(len(buf)) {
+				// In ideal math: buf[i] += uint32(d * am)
+				D := x1fSquared
+				D *= d
+				D /= twoOverS
+				buf[i] += uint32(D)
+			}
+		}
+
+		x = xNext
+	}
+}
+
+func fixedAccumulateOpSrc(dst []uint8, src []uint32) {
+	acc := int2ϕ(0)
+	for i, v := range src {
+		acc += int2ϕ(v)
+		a := acc
+		if a < 0 {
+			a = -a
+		}
+		a >>= 2*ϕ - 8
+		if a > 0xff {
+			a = 0xff
+		}
+		dst[i] = uint8(a)
+	}
+}
+
+func fixedAccumulateOpOver(dst []uint8, src []uint32) {
+	acc := int2ϕ(0)
+	for i, v := range src {
+		acc += int2ϕ(v)
+		a := acc
+		if a < 0 {
+			a = -a
+		}
+		a >>= 2*ϕ - 16
+		if a > 0xffff {
+			a = 0xffff
+		}
+		// This algorithm comes from the standard library's image/draw package.
+		dstA := uint32(dst[i]) * 0x101
+		maskA := uint32(a)
+		outA := dstA*(0xffff-maskA)/0xffff + maskA
+		dst[i] = uint8(outA >> 8)
+	}
+}
+
+func fixedAccumulateMask(buf []uint32) {
+	acc := int2ϕ(0)
+	for i, v := range buf {
+		acc += int2ϕ(v)
+		a := acc
+		if a < 0 {
+			a = -a
+		}
+		a >>= 2*ϕ - 16
+		if a > 0xffff {
+			a = 0xffff
+		}
+		buf[i] = uint32(a)
+	}
+}

vector/vector.go

@@ -25,6 +25,26 @@ import (
 	"golang.org/x/image/math/f32"
 )
 
+// floatingPointMathThreshold is the width or hight above which the rasterizer
+// chooses to used floating point math instead of fixed point math.
+//
+// Both implementations of line segmentation rasterization (see raster_fixed.go
+// and raster_floating.go) implement the same algorithm (in ideal, infinite
+// precision math) but they perform differently in practice. The fixed point
+// math version is roughtly 1.25x faster (on GOARCH=amd64) on the benchmarks,
+// but at sufficiently large scales, the computations will overflow and hence
+// show rendering artifacts. The floating point math version has more
+// consistent quality over larger scales, but it is significantly slower.
+//
+// This constant determines when to use the faster implementation and when to
+// use the better quality implementation.
+//
+// The rationale for this particular value is that TestRasterizePolygon in
+// vector_test.go checks the rendering quality of polygon edges at various
+// angles, inscribed in a circle of diameter 2048. It may be that a higher
+// value would still produce acceptable quality, but 2048 seems to work.
+const floatingPointMathThreshold = 2048
+
 func midPoint(p, q f32.Vec2) f32.Vec2 {
 	return f32.Vec2{
 		(p[0] + q[0]) * 0.5,
@@ -52,10 +72,9 @@ func clamp(i, width int32) uint {
 // NewRasterizer returns a new Rasterizer whose rendered mask image is bounded
 // by the given width and height.
 func NewRasterizer(w, h int) *Rasterizer {
-	return &Rasterizer{
+	z := &Rasterizer{}
-		bufF32: make([]float32, w*h),
+	z.Reset(w, h)
-		size:   image.Point{w, h},
+	return z
-	}
 }
 
 // Raster is a 2-D vector graphics rasterizer.
@@ -77,11 +96,11 @@ type Rasterizer struct {
 	//	bufU32[i] = math.Float32bits(x + math.Float32frombits(bufU32[i]))
 	//
 	// See golang.org/issue/17220 for some discussion.
-	//
-	// TODO: use bufU32 in the fixed point math implementation.
 	bufF32 []float32
 	bufU32 []uint32
 
+	useFloatingPointMath bool
+
 	size  image.Point
 	first f32.Vec2
 	pen   f32.Vec2
@@ -99,18 +118,33 @@ type Rasterizer struct {
 //
 // This includes setting z.DrawOp to draw.Over.
 func (z *Rasterizer) Reset(w, h int) {
-	if n := w * h; n > cap(z.bufF32) {
-		z.bufF32 = make([]float32, n)
-	} else {
-		z.bufF32 = z.bufF32[:n]
-		for i := range z.bufF32 {
-			z.bufF32[i] = 0
-		}
-	}
 	z.size = image.Point{w, h}
 	z.first = f32.Vec2{}
 	z.pen = f32.Vec2{}
 	z.DrawOp = draw.Over
+
+	z.useFloatingPointMath = w > floatingPointMathThreshold || h > floatingPointMathThreshold
+
+	// Make z.bufF32 or z.bufU32 large enough to hold w*h samples.
+	if z.useFloatingPointMath {
+		if n := w * h; n > cap(z.bufF32) {
+			z.bufF32 = make([]float32, n)
+		} else {
+			z.bufF32 = z.bufF32[:n]
+			for i := range z.bufF32 {
+				z.bufF32[i] = 0
+			}
+		}
+	} else {
+		if n := w * h; n > cap(z.bufU32) {
+			z.bufU32 = make([]uint32, n)
+		} else {
+			z.bufU32 = z.bufU32[:n]
+			for i := range z.bufU32 {
+				z.bufU32[i] = 0
+			}
+		}
+	}
 }
 
 // Size returns the width and height passed to NewRasterizer or Reset.
@@ -147,8 +181,11 @@ func (z *Rasterizer) MoveTo(a f32.Vec2) {
 //
 // The coordinates are allowed to be out of the Rasterizer's bounds.
 func (z *Rasterizer) LineTo(b f32.Vec2) {
-	// TODO: add a fixed point math implementation.
+	if z.useFloatingPointMath {
-	z.floatingLineTo(b)
+		z.floatingLineTo(b)
+	} else {
+		z.fixedLineTo(b)
+	}
 }
 
 // QuadTo adds a quadratic Bézier segment, from the pen via b to c, and moves
@@ -258,22 +295,29 @@ func (z *Rasterizer) Draw(dst draw.Image, r image.Rectangle, src image.Image, sp
 }
 
 func (z *Rasterizer) accumulateMask() {
-	if n := z.size.X * z.size.Y; n > cap(z.bufU32) {
+	if z.useFloatingPointMath {
-		z.bufU32 = make([]uint32, n)
+		if n := z.size.X * z.size.Y; n > cap(z.bufU32) {
+			z.bufU32 = make([]uint32, n)
+		} else {
+			z.bufU32 = z.bufU32[:n]
+		}
+		floatingAccumulateMask(z.bufU32, z.bufF32)
 	} else {
-		z.bufU32 = z.bufU32[:n]
+		fixedAccumulateMask(z.bufU32)
 	}
-	floatingAccumulateMask(z.bufU32, z.bufF32)
 }
 
 func (z *Rasterizer) rasterizeDstAlphaSrcOpaqueOpOver(dst *image.Alpha, r image.Rectangle) {
 	// TODO: add SIMD implementations.
-	// TODO: add a fixed point math implementation.
 	// TODO: non-zero vs even-odd winding?
 	if r == dst.Bounds() && r == z.Bounds() {
 		// We bypass the z.accumulateMask step and convert straight from
-		// z.bufF32 to dst.Pix.
+		// z.bufF32 or z.bufU32 to dst.Pix.
-		floatingAccumulateOpOver(dst.Pix, z.bufF32)
+		if z.useFloatingPointMath {
+			floatingAccumulateOpOver(dst.Pix, z.bufF32)
+		} else {
+			fixedAccumulateOpOver(dst.Pix, z.bufU32)
+		}
 		return
 	}
 
@@ -294,12 +338,15 @@ func (z *Rasterizer) rasterizeDstAlphaSrcOpaqueOpOver(dst *image.Alpha, r image.
 
 func (z *Rasterizer) rasterizeDstAlphaSrcOpaqueOpSrc(dst *image.Alpha, r image.Rectangle) {
 	// TODO: add SIMD implementations.
-	// TODO: add a fixed point math implementation.
 	// TODO: non-zero vs even-odd winding?
 	if r == dst.Bounds() && r == z.Bounds() {
 		// We bypass the z.accumulateMask step and convert straight from
-		// z.bufF32 to dst.Pix.
+		// z.bufF32 or z.bufU32 to dst.Pix.
-		floatingAccumulateOpSrc(dst.Pix, z.bufF32)
+		if z.useFloatingPointMath {
+			floatingAccumulateOpSrc(dst.Pix, z.bufF32)
+		} else {
+			fixedAccumulateOpSrc(dst.Pix, z.bufU32)
+		}
 		return
 	}