From 87013da1482c3e5452b2611b9d9dc944ca7d334a Mon Sep 17 00:00:00 2001
From: Nigel Tao <nigeltao@golang.org>
Date: Fri, 13 Mar 2015 17:44:34 +1100
Subject: [PATCH] draw: implement NearestNeighbor and ApproxBiLinear Transform.

Change-Id: I70a5e3703dea436354e9591fce7b704ec749c2d1
Reviewed-on: https://go-review.googlesource.com/7541
Reviewed-by: Rob Pike <r@golang.org>
---
 draw/example_test.go                |  11 +-
 draw/gen.go                         | 153 ++++-
 draw/impl.go                        | 888 ++++++++++++++++++++++++++--
 draw/scale.go                       |  69 +++
 draw/scale_test.go                  | 143 +++--
 testdata/go-turns-two-rotate-ab.png | Bin 0 -> 8761 bytes
 testdata/go-turns-two-rotate-nn.png | Bin 0 -> 4993 bytes
 7 files changed, 1135 insertions(+), 129 deletions(-)
 create mode 100644 testdata/go-turns-two-rotate-ab.png
 create mode 100644 testdata/go-turns-two-rotate-nn.png

diff --git a/draw/example_test.go b/draw/example_test.go
index 02e4c5d..f8545ad 100644
--- a/draw/example_test.go
+++ b/draw/example_test.go
@@ -9,7 +9,6 @@ import (
 	"image"
 	"image/png"
 	"log"
-	"math"
 	"os"
 
 	"golang.org/x/image/draw"
@@ -34,19 +33,17 @@ func ExampleDraw() {
 		draw.ApproxBiLinear,
 		draw.CatmullRom,
 	}
-	c, s := math.Cos(math.Pi/3), math.Sin(math.Pi/3)
+	const cos60, sin60 = 0.5, 0.866025404
 	t := &f64.Aff3{
-		+2 * c, -2 * s, 100,
-		+2 * s, +2 * c, 100,
+		+2 * cos60, -2 * sin60, 100,
+		+2 * sin60, +2 * cos60, 100,
 	}
 
 	draw.Copy(dst, image.Point{20, 30}, src, sr, nil)
 	for i, q := range qs {
 		q.Scale(dst, image.Rect(200+10*i, 100*i, 600+10*i, 150+100*i), src, sr, nil)
 	}
-	// TODO: delete the "_ = t" and uncomment this when Transform is implemented.
-	// draw.NearestNeighbor.Transform(dst, t, src, sr, nil)
-	_ = t
+	draw.NearestNeighbor.Transform(dst, t, src, sr, nil)
 
 	// Change false to true to write the resultant image to disk.
 	if false {
diff --git a/draw/gen.go b/draw/gen.go
index bc11b3e..587b969 100644
--- a/draw/gen.go
+++ b/draw/gen.go
@@ -27,12 +27,13 @@ func main() {
 		"package draw\n\nimport (\n" +
 		"\"image\"\n" +
 		"\"image/color\"\n" +
+		"\"math\"\n" +
 		"\n" +
 		"\"golang.org/x/image/math/f64\"\n" +
 		")\n")
 
-	gen(w, "nnInterpolator", codeNNScaleLeaf)
-	gen(w, "ablInterpolator", codeABLScaleLeaf)
+	gen(w, "nnInterpolator", codeNNScaleLeaf, codeNNTransformLeaf)
+	gen(w, "ablInterpolator", codeABLScaleLeaf, codeABLTransformLeaf)
 	genKernel(w)
 
 	if *debug {
@@ -90,14 +91,16 @@ type data struct {
 	receiver string
 }
 
-func gen(w *bytes.Buffer, receiver string, code string) {
+func gen(w *bytes.Buffer, receiver string, codes ...string) {
 	expn(w, codeRoot, &data{receiver: receiver})
-	for _, t := range dsTypes {
-		expn(w, code, &data{
-			dType:    t.dType,
-			sType:    t.sType,
-			receiver: receiver,
-		})
+	for _, code := range codes {
+		for _, t := range dsTypes {
+			expn(w, code, &data{
+				dType:    t.dType,
+				sType:    t.sType,
+				receiver: receiver,
+			})
+		}
 	}
 }
 
@@ -227,7 +230,7 @@ func expnDollar(prefix, dollar, suffix string, d *data) string {
 				"dstColorRGBA64.G = uint16(%sg)\n"+
 				"dstColorRGBA64.B = uint16(%sb)\n"+
 				"dstColorRGBA64.A = uint16(%sa)\n"+
-				"dst.Set(dr.Min.X+int(%s), dr.Min.Y+int(%s), dstColor)",
+				"dst.Set(%s, %s, dstColor)",
 				args[2], args[2], args[2], args[2],
 				args[0], args[1],
 			)
@@ -236,8 +239,7 @@ func expnDollar(prefix, dollar, suffix string, d *data) string {
 				"dst.Pix[d+0] = uint8(uint32(%sr) >> 8)\n"+
 				"dst.Pix[d+1] = uint8(uint32(%sg) >> 8)\n"+
 				"dst.Pix[d+2] = uint8(uint32(%sb) >> 8)\n"+
-				"dst.Pix[d+3] = uint8(uint32(%sa) >> 8)\n"+
-				"d += 4",
+				"dst.Pix[d+3] = uint8(uint32(%sa) >> 8)",
 				args[2], args[2], args[2], args[2],
 			)
 		}
@@ -256,7 +258,7 @@ func expnDollar(prefix, dollar, suffix string, d *data) string {
 				"dstColorRGBA64.G = ftou(%sg * %s)\n"+
 				"dstColorRGBA64.B = ftou(%sb * %s)\n"+
 				"dstColorRGBA64.A = ftou(%sa * %s)\n"+
-				"dst.Set(dr.Min.X+int(%s), dr.Min.Y+int(%s), dstColor)",
+				"dst.Set(%s, %s, dstColor)",
 				args[2], args[3], args[2], args[3], args[2], args[3], args[2], args[3],
 				args[0], args[1],
 			)
@@ -292,14 +294,14 @@ func expnDollar(prefix, dollar, suffix string, d *data) string {
 			log.Fatalf("bad sType %q", d.sType)
 		case "image.Image", "*image.Gray", "*image.NRGBA", "*image.Uniform", "*image.YCbCr": // TODO: separate code for concrete types.
 			fmt.Fprintf(buf, "%sr%s, %sg%s, %sb%s, %sa%s := "+
-				"src.At(sr.Min.X + int(%s), sr.Min.Y+int(%s)).RGBA()\n",
+				"src.At(%s, %s).RGBA()\n",
 				lhs, tmp, lhs, tmp, lhs, tmp, lhs, tmp,
 				args[0], args[1],
 			)
 		case "*image.RGBA":
 			// TODO: there's no need to multiply by 0x101 if the next thing
 			// we're going to do is shift right by 8.
-			fmt.Fprintf(buf, "%si := src.PixOffset(sr.Min.X + int(%s), sr.Min.Y+int(%s))\n"+
+			fmt.Fprintf(buf, "%si := src.PixOffset(%s, %s)\n"+
 				"%sr%s := uint32(src.Pix[%si+0]) * 0x101\n"+
 				"%sg%s := uint32(src.Pix[%si+1]) * 0x101\n"+
 				"%sb%s := uint32(src.Pix[%si+2]) * 0x101\n"+
@@ -327,6 +329,12 @@ func expnDollar(prefix, dollar, suffix string, d *data) string {
 
 		return strings.TrimSpace(buf.String())
 
+	case "tweakDx":
+		if d.dType == "*image.RGBA" {
+			return strings.Replace(suffix, "dx++", "dx, d = dx+1, d+4", 1)
+		}
+		return suffix
+
 	case "tweakDy":
 		if d.dType == "*image.RGBA" {
 			return strings.Replace(suffix, "for dy, s", "for _, s", 1)
@@ -428,8 +436,15 @@ const (
 			}
 		}
 
-		func (z $receiver) Transform(dst Image, m *f64.Aff3, src image.Image, sr image.Rectangle, opts *Options) {
-			panic("unimplemented")
+		func (z $receiver) Transform(dst Image, s2d *f64.Aff3, src image.Image, sr image.Rectangle, opts *Options) {
+			dr := transformRect(s2d, &sr)
+			// adr is the affected destination pixels, relative to dr.Min.
+			adr := dst.Bounds().Intersect(dr).Sub(dr.Min)
+			if adr.Empty() || sr.Empty() {
+				return
+			}
+			d2s := invert(s2d)
+			z.transform_Image_Image(dst, dr, adr, &d2s, src, sr)
 		}
 	`
 
@@ -443,10 +458,30 @@ const (
 			for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 				sy := (2*uint64(dy) + 1) * sh / dh2
 				$preInner
-				for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+				$tweakDx for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
 					sx := (2*uint64(dx) + 1) * sw / dw2
-					p := $srcu[sx, sy]
-					$outputu[dx, dy, p]
+					p := $srcu[sr.Min.X + int(sx), sr.Min.Y + int(sy)]
+					$outputu[dr.Min.X + int(dx), dr.Min.Y + int(dy), p]
+				}
+			}
+		}
+	`
+
+	codeNNTransformLeaf = `
+		func (nnInterpolator) transform_$dTypeRN_$sTypeRN(dst $dType, dr, adr image.Rectangle, d2s *f64.Aff3, src $sType, sr image.Rectangle) {
+			$preOuter
+			for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+				dyf := float64(dr.Min.Y + int(dy)) + 0.5
+				$preInner
+				$tweakDx for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+					dxf := float64(dr.Min.X + int(dx)) + 0.5
+					sx0 := int(math.Floor(d2s[0]*dxf + d2s[1]*dyf + d2s[2]))
+					sy0 := int(math.Floor(d2s[3]*dxf + d2s[4]*dyf + d2s[5]))
+					if !(image.Point{sx0, sy0}).In(sr) {
+						continue
+					}
+					p := $srcu[sx0, sy0]
+					$outputu[dr.Min.X + int(dx), dr.Min.Y + int(dy), p]
 				}
 			}
 		}
@@ -458,9 +493,14 @@ const (
 			sh := int32(sr.Dy())
 			yscale := float64(sh) / float64(dr.Dy())
 			xscale := float64(sw) / float64(dr.Dx())
+			swMinus1, shMinus1 := sw - 1, sh - 1
 			$preOuter
+
 			for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 				sy := (float64(dy)+0.5)*yscale - 0.5
+				// If sy < 0, we will clamp sy0 to 0 anyway, so it doesn't matter if
+				// we say int32(sy) instead of int32(math.Floor(sy)). Similarly for
+				// sx, below.
 				sy0 := int32(sy)
 				yFrac0 := sy - float64(sy0)
 				yFrac1 := 1 - yFrac0
@@ -468,12 +508,13 @@ const (
 				if sy < 0 {
 					sy0, sy1 = 0, 0
 					yFrac0, yFrac1 = 0, 1
-				} else if sy1 >= sh {
-					sy1 = sy0
+				} else if sy1 > shMinus1 {
+					sy0, sy1 = shMinus1, shMinus1
 					yFrac0, yFrac1 = 1, 0
 				}
 				$preInner
-				for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+
+				$tweakDx for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
 					sx := (float64(dx)+0.5)*xscale - 0.5
 					sx0 := int32(sx)
 					xFrac0 := sx - float64(sx0)
@@ -482,10 +523,66 @@ const (
 					if sx < 0 {
 						sx0, sx1 = 0, 0
 						xFrac0, xFrac1 = 0, 1
-					} else if sx1 >= sw {
-						sx1 = sx0
+					} else if sx1 > swMinus1 {
+						sx0, sx1 = swMinus1, swMinus1
 						xFrac0, xFrac1 = 1, 0
 					}
+
+					s00 := $srcf[sr.Min.X + int(sx0), sr.Min.Y + int(sy0)]
+					s10 := $srcf[sr.Min.X + int(sx1), sr.Min.Y + int(sy0)]
+					$blend[xFrac1, s00, xFrac0, s10]
+					s01 := $srcf[sr.Min.X + int(sx0), sr.Min.Y + int(sy1)]
+					s11 := $srcf[sr.Min.X + int(sx1), sr.Min.Y + int(sy1)]
+					$blend[xFrac1, s01, xFrac0, s11]
+					$blend[yFrac1, s10, yFrac0, s11]
+					$outputu[dr.Min.X + int(dx), dr.Min.Y + int(dy), s11]
+				}
+			}
+		}
+	`
+
+	codeABLTransformLeaf = `
+		func (ablInterpolator) transform_$dTypeRN_$sTypeRN(dst $dType, dr, adr image.Rectangle, d2s *f64.Aff3, src $sType, sr image.Rectangle) {
+			$preOuter
+			for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+				dyf := float64(dr.Min.Y + int(dy)) + 0.5
+				$preInner
+				$tweakDx for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+					dxf := float64(dr.Min.X + int(dx)) + 0.5
+					sx := d2s[0]*dxf + d2s[1]*dyf + d2s[2]
+					sy := d2s[3]*dxf + d2s[4]*dyf + d2s[5]
+					if !(image.Point{int(math.Floor(sx)), int(math.Floor(sy))}).In(sr) {
+						continue
+					}
+
+					sx -= 0.5
+					sxf := math.Floor(sx)
+					xFrac0 := sx - sxf
+					xFrac1 := 1 - xFrac0
+					sx0 := int(sxf)
+					sx1 := sx0 + 1
+					if sx0 < sr.Min.X {
+						sx0, sx1 = sr.Min.X, sr.Min.X
+						xFrac0, xFrac1 = 0, 1
+					} else if sx1 >= sr.Max.X {
+						sx0, sx1 = sr.Max.X-1, sr.Max.X-1
+						xFrac0, xFrac1 = 1, 0
+					}
+
+					sy -= 0.5
+					syf := math.Floor(sy)
+					yFrac0 := sy - syf
+					yFrac1 := 1 - yFrac0
+					sy0 := int(syf)
+					sy1 := sy0 + 1
+					if sy0 < sr.Min.Y {
+						sy0, sy1 = sr.Min.Y, sr.Min.Y
+						yFrac0, yFrac1 = 0, 1
+					} else if sy1 >= sr.Max.Y {
+						sy0, sy1 = sr.Max.Y-1, sr.Max.Y-1
+						yFrac0, yFrac1 = 1, 0
+					}
+
 					s00 := $srcf[sx0, sy0]
 					s10 := $srcf[sx1, sy0]
 					$blend[xFrac1, s00, xFrac0, s10]
@@ -493,7 +590,7 @@ const (
 					s11 := $srcf[sx1, sy1]
 					$blend[xFrac1, s01, xFrac0, s11]
 					$blend[yFrac1, s10, yFrac0, s11]
-					$outputu[dx, dy, s11]
+					$outputu[dr.Min.X + int(dx), dr.Min.Y + int(dy), s11]
 				}
 			}
 		}
@@ -540,7 +637,7 @@ const (
 				for _, s := range z.horizontal.sources {
 					var pr, pg, pb, pa float64
 					for _, c := range z.horizontal.contribs[s.i:s.j] {
-						p += $srcf[c.coord, y] * c.weight
+						p += $srcf[sr.Min.X + int(c.coord), sr.Min.Y + int(y)] * c.weight
 					}
 					tmp[t] = [4]float64{
 						pr * s.invTotalWeightFFFF,
@@ -568,7 +665,7 @@ const (
 						pb += p[2] * c.weight
 						pa += p[3] * c.weight
 					}
-					$outputf[dx, adr.Min.Y+dy, p, s.invTotalWeight]
+					$outputf[dr.Min.X + int(dx), dr.Min.Y + int(adr.Min.Y + dy), p, s.invTotalWeight]
 				}
 			}
 		}
diff --git a/draw/impl.go b/draw/impl.go
index ff9f988..cba9349 100644
--- a/draw/impl.go
+++ b/draw/impl.go
@@ -5,6 +5,7 @@ package draw
 import (
 	"image"
 	"image/color"
+	"math"
 
 	"golang.org/x/image/math/f64"
 )
@@ -46,8 +47,15 @@ func (z nnInterpolator) Scale(dst Image, dr image.Rectangle, src image.Image, sr
 	}
 }
 
-func (z nnInterpolator) Transform(dst Image, m *f64.Aff3, src image.Image, sr image.Rectangle, opts *Options) {
-	panic("unimplemented")
+func (z nnInterpolator) Transform(dst Image, s2d *f64.Aff3, src image.Image, sr image.Rectangle, opts *Options) {
+	dr := transformRect(s2d, &sr)
+	// adr is the affected destination pixels, relative to dr.Min.
+	adr := dst.Bounds().Intersect(dr).Sub(dr.Min)
+	if adr.Empty() || sr.Empty() {
+		return
+	}
+	d2s := invert(s2d)
+	z.transform_Image_Image(dst, dr, adr, &d2s, src, sr)
 }
 
 func (nnInterpolator) scale_RGBA_Gray(dst *image.RGBA, dr, adr image.Rectangle, src *image.Gray, sr image.Rectangle) {
@@ -58,14 +66,13 @@ func (nnInterpolator) scale_RGBA_Gray(dst *image.RGBA, dr, adr image.Rectangle,
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (2*uint64(dy) + 1) * sh / dh2
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (2*uint64(dx) + 1) * sw / dw2
 			pr, pg, pb, pa := src.At(sr.Min.X+int(sx), sr.Min.Y+int(sy)).RGBA()
 			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
 			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
 			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
 			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
-			d += 4
 		}
 	}
 }
@@ -78,14 +85,13 @@ func (nnInterpolator) scale_RGBA_NRGBA(dst *image.RGBA, dr, adr image.Rectangle,
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (2*uint64(dy) + 1) * sh / dh2
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (2*uint64(dx) + 1) * sw / dw2
 			pr, pg, pb, pa := src.At(sr.Min.X+int(sx), sr.Min.Y+int(sy)).RGBA()
 			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
 			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
 			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
 			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
-			d += 4
 		}
 	}
 }
@@ -98,7 +104,7 @@ func (nnInterpolator) scale_RGBA_RGBA(dst *image.RGBA, dr, adr image.Rectangle,
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (2*uint64(dy) + 1) * sh / dh2
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (2*uint64(dx) + 1) * sw / dw2
 			pi := src.PixOffset(sr.Min.X+int(sx), sr.Min.Y+int(sy))
 			pr := uint32(src.Pix[pi+0]) * 0x101
@@ -109,7 +115,6 @@ func (nnInterpolator) scale_RGBA_RGBA(dst *image.RGBA, dr, adr image.Rectangle,
 			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
 			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
 			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
-			d += 4
 		}
 	}
 }
@@ -122,14 +127,13 @@ func (nnInterpolator) scale_RGBA_Uniform(dst *image.RGBA, dr, adr image.Rectangl
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (2*uint64(dy) + 1) * sh / dh2
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (2*uint64(dx) + 1) * sw / dw2
 			pr, pg, pb, pa := src.At(sr.Min.X+int(sx), sr.Min.Y+int(sy)).RGBA()
 			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
 			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
 			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
 			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
-			d += 4
 		}
 	}
 }
@@ -142,14 +146,13 @@ func (nnInterpolator) scale_RGBA_YCbCr(dst *image.RGBA, dr, adr image.Rectangle,
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (2*uint64(dy) + 1) * sh / dh2
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (2*uint64(dx) + 1) * sw / dw2
 			pr, pg, pb, pa := src.At(sr.Min.X+int(sx), sr.Min.Y+int(sy)).RGBA()
 			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
 			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
 			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
 			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
-			d += 4
 		}
 	}
 }
@@ -162,14 +165,13 @@ func (nnInterpolator) scale_RGBA_Image(dst *image.RGBA, dr, adr image.Rectangle,
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (2*uint64(dy) + 1) * sh / dh2
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (2*uint64(dx) + 1) * sw / dw2
 			pr, pg, pb, pa := src.At(sr.Min.X+int(sx), sr.Min.Y+int(sy)).RGBA()
 			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
 			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
 			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
 			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
-			d += 4
 		}
 	}
 }
@@ -195,6 +197,152 @@ func (nnInterpolator) scale_Image_Image(dst Image, dr, adr image.Rectangle, src
 	}
 }
 
+func (nnInterpolator) transform_RGBA_Gray(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.Gray, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx0 := int(math.Floor(d2s[0]*dxf + d2s[1]*dyf + d2s[2]))
+			sy0 := int(math.Floor(d2s[3]*dxf + d2s[4]*dyf + d2s[5]))
+			if !(image.Point{sx0, sy0}).In(sr) {
+				continue
+			}
+			pr, pg, pb, pa := src.At(sx0, sy0).RGBA()
+			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
+			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
+			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
+			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
+		}
+	}
+}
+
+func (nnInterpolator) transform_RGBA_NRGBA(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.NRGBA, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx0 := int(math.Floor(d2s[0]*dxf + d2s[1]*dyf + d2s[2]))
+			sy0 := int(math.Floor(d2s[3]*dxf + d2s[4]*dyf + d2s[5]))
+			if !(image.Point{sx0, sy0}).In(sr) {
+				continue
+			}
+			pr, pg, pb, pa := src.At(sx0, sy0).RGBA()
+			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
+			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
+			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
+			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
+		}
+	}
+}
+
+func (nnInterpolator) transform_RGBA_RGBA(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.RGBA, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx0 := int(math.Floor(d2s[0]*dxf + d2s[1]*dyf + d2s[2]))
+			sy0 := int(math.Floor(d2s[3]*dxf + d2s[4]*dyf + d2s[5]))
+			if !(image.Point{sx0, sy0}).In(sr) {
+				continue
+			}
+			pi := src.PixOffset(sx0, sy0)
+			pr := uint32(src.Pix[pi+0]) * 0x101
+			pg := uint32(src.Pix[pi+1]) * 0x101
+			pb := uint32(src.Pix[pi+2]) * 0x101
+			pa := uint32(src.Pix[pi+3]) * 0x101
+			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
+			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
+			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
+			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
+		}
+	}
+}
+
+func (nnInterpolator) transform_RGBA_Uniform(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.Uniform, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx0 := int(math.Floor(d2s[0]*dxf + d2s[1]*dyf + d2s[2]))
+			sy0 := int(math.Floor(d2s[3]*dxf + d2s[4]*dyf + d2s[5]))
+			if !(image.Point{sx0, sy0}).In(sr) {
+				continue
+			}
+			pr, pg, pb, pa := src.At(sx0, sy0).RGBA()
+			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
+			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
+			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
+			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
+		}
+	}
+}
+
+func (nnInterpolator) transform_RGBA_YCbCr(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.YCbCr, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx0 := int(math.Floor(d2s[0]*dxf + d2s[1]*dyf + d2s[2]))
+			sy0 := int(math.Floor(d2s[3]*dxf + d2s[4]*dyf + d2s[5]))
+			if !(image.Point{sx0, sy0}).In(sr) {
+				continue
+			}
+			pr, pg, pb, pa := src.At(sx0, sy0).RGBA()
+			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
+			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
+			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
+			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
+		}
+	}
+}
+
+func (nnInterpolator) transform_RGBA_Image(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src image.Image, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx0 := int(math.Floor(d2s[0]*dxf + d2s[1]*dyf + d2s[2]))
+			sy0 := int(math.Floor(d2s[3]*dxf + d2s[4]*dyf + d2s[5]))
+			if !(image.Point{sx0, sy0}).In(sr) {
+				continue
+			}
+			pr, pg, pb, pa := src.At(sx0, sy0).RGBA()
+			dst.Pix[d+0] = uint8(uint32(pr) >> 8)
+			dst.Pix[d+1] = uint8(uint32(pg) >> 8)
+			dst.Pix[d+2] = uint8(uint32(pb) >> 8)
+			dst.Pix[d+3] = uint8(uint32(pa) >> 8)
+		}
+	}
+}
+
+func (nnInterpolator) transform_Image_Image(dst Image, dr, adr image.Rectangle, d2s *f64.Aff3, src image.Image, sr image.Rectangle) {
+	dstColorRGBA64 := &color.RGBA64{}
+	dstColor := color.Color(dstColorRGBA64)
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx0 := int(math.Floor(d2s[0]*dxf + d2s[1]*dyf + d2s[2]))
+			sy0 := int(math.Floor(d2s[3]*dxf + d2s[4]*dyf + d2s[5]))
+			if !(image.Point{sx0, sy0}).In(sr) {
+				continue
+			}
+			pr, pg, pb, pa := src.At(sx0, sy0).RGBA()
+			dstColorRGBA64.R = uint16(pr)
+			dstColorRGBA64.G = uint16(pg)
+			dstColorRGBA64.B = uint16(pb)
+			dstColorRGBA64.A = uint16(pa)
+			dst.Set(dr.Min.X+int(dx), dr.Min.Y+int(dy), dstColor)
+		}
+	}
+}
+
 func (z ablInterpolator) Scale(dst Image, dr image.Rectangle, src image.Image, sr image.Rectangle, opts *Options) {
 	// adr is the affected destination pixels, relative to dr.Min.
 	adr := dst.Bounds().Intersect(dr).Sub(dr.Min)
@@ -232,8 +380,15 @@ func (z ablInterpolator) Scale(dst Image, dr image.Rectangle, src image.Image, s
 	}
 }
 
-func (z ablInterpolator) Transform(dst Image, m *f64.Aff3, src image.Image, sr image.Rectangle, opts *Options) {
-	panic("unimplemented")
+func (z ablInterpolator) Transform(dst Image, s2d *f64.Aff3, src image.Image, sr image.Rectangle, opts *Options) {
+	dr := transformRect(s2d, &sr)
+	// adr is the affected destination pixels, relative to dr.Min.
+	adr := dst.Bounds().Intersect(dr).Sub(dr.Min)
+	if adr.Empty() || sr.Empty() {
+		return
+	}
+	d2s := invert(s2d)
+	z.transform_Image_Image(dst, dr, adr, &d2s, src, sr)
 }
 
 func (ablInterpolator) scale_RGBA_Gray(dst *image.RGBA, dr, adr image.Rectangle, src *image.Gray, sr image.Rectangle) {
@@ -241,8 +396,13 @@ func (ablInterpolator) scale_RGBA_Gray(dst *image.RGBA, dr, adr image.Rectangle,
 	sh := int32(sr.Dy())
 	yscale := float64(sh) / float64(dr.Dy())
 	xscale := float64(sw) / float64(dr.Dx())
+	swMinus1, shMinus1 := sw-1, sh-1
+
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (float64(dy)+0.5)*yscale - 0.5
+		// If sy < 0, we will clamp sy0 to 0 anyway, so it doesn't matter if
+		// we say int32(sy) instead of int32(math.Floor(sy)). Similarly for
+		// sx, below.
 		sy0 := int32(sy)
 		yFrac0 := sy - float64(sy0)
 		yFrac1 := 1 - yFrac0
@@ -250,12 +410,13 @@ func (ablInterpolator) scale_RGBA_Gray(dst *image.RGBA, dr, adr image.Rectangle,
 		if sy < 0 {
 			sy0, sy1 = 0, 0
 			yFrac0, yFrac1 = 0, 1
-		} else if sy1 >= sh {
-			sy1 = sy0
+		} else if sy1 > shMinus1 {
+			sy0, sy1 = shMinus1, shMinus1
 			yFrac0, yFrac1 = 1, 0
 		}
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (float64(dx)+0.5)*xscale - 0.5
 			sx0 := int32(sx)
 			xFrac0 := sx - float64(sx0)
@@ -264,10 +425,11 @@ func (ablInterpolator) scale_RGBA_Gray(dst *image.RGBA, dr, adr image.Rectangle,
 			if sx < 0 {
 				sx0, sx1 = 0, 0
 				xFrac0, xFrac1 = 0, 1
-			} else if sx1 >= sw {
-				sx1 = sx0
+			} else if sx1 > swMinus1 {
+				sx0, sx1 = swMinus1, swMinus1
 				xFrac0, xFrac1 = 1, 0
 			}
+
 			s00ru, s00gu, s00bu, s00au := src.At(sr.Min.X+int(sx0), sr.Min.Y+int(sy0)).RGBA()
 			s00r := float64(s00ru)
 			s00g := float64(s00gu)
@@ -304,7 +466,6 @@ func (ablInterpolator) scale_RGBA_Gray(dst *image.RGBA, dr, adr image.Rectangle,
 			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
 			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
 			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
-			d += 4
 		}
 	}
 }
@@ -314,8 +475,13 @@ func (ablInterpolator) scale_RGBA_NRGBA(dst *image.RGBA, dr, adr image.Rectangle
 	sh := int32(sr.Dy())
 	yscale := float64(sh) / float64(dr.Dy())
 	xscale := float64(sw) / float64(dr.Dx())
+	swMinus1, shMinus1 := sw-1, sh-1
+
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (float64(dy)+0.5)*yscale - 0.5
+		// If sy < 0, we will clamp sy0 to 0 anyway, so it doesn't matter if
+		// we say int32(sy) instead of int32(math.Floor(sy)). Similarly for
+		// sx, below.
 		sy0 := int32(sy)
 		yFrac0 := sy - float64(sy0)
 		yFrac1 := 1 - yFrac0
@@ -323,12 +489,13 @@ func (ablInterpolator) scale_RGBA_NRGBA(dst *image.RGBA, dr, adr image.Rectangle
 		if sy < 0 {
 			sy0, sy1 = 0, 0
 			yFrac0, yFrac1 = 0, 1
-		} else if sy1 >= sh {
-			sy1 = sy0
+		} else if sy1 > shMinus1 {
+			sy0, sy1 = shMinus1, shMinus1
 			yFrac0, yFrac1 = 1, 0
 		}
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (float64(dx)+0.5)*xscale - 0.5
 			sx0 := int32(sx)
 			xFrac0 := sx - float64(sx0)
@@ -337,10 +504,11 @@ func (ablInterpolator) scale_RGBA_NRGBA(dst *image.RGBA, dr, adr image.Rectangle
 			if sx < 0 {
 				sx0, sx1 = 0, 0
 				xFrac0, xFrac1 = 0, 1
-			} else if sx1 >= sw {
-				sx1 = sx0
+			} else if sx1 > swMinus1 {
+				sx0, sx1 = swMinus1, swMinus1
 				xFrac0, xFrac1 = 1, 0
 			}
+
 			s00ru, s00gu, s00bu, s00au := src.At(sr.Min.X+int(sx0), sr.Min.Y+int(sy0)).RGBA()
 			s00r := float64(s00ru)
 			s00g := float64(s00gu)
@@ -377,7 +545,6 @@ func (ablInterpolator) scale_RGBA_NRGBA(dst *image.RGBA, dr, adr image.Rectangle
 			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
 			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
 			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
-			d += 4
 		}
 	}
 }
@@ -387,8 +554,13 @@ func (ablInterpolator) scale_RGBA_RGBA(dst *image.RGBA, dr, adr image.Rectangle,
 	sh := int32(sr.Dy())
 	yscale := float64(sh) / float64(dr.Dy())
 	xscale := float64(sw) / float64(dr.Dx())
+	swMinus1, shMinus1 := sw-1, sh-1
+
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (float64(dy)+0.5)*yscale - 0.5
+		// If sy < 0, we will clamp sy0 to 0 anyway, so it doesn't matter if
+		// we say int32(sy) instead of int32(math.Floor(sy)). Similarly for
+		// sx, below.
 		sy0 := int32(sy)
 		yFrac0 := sy - float64(sy0)
 		yFrac1 := 1 - yFrac0
@@ -396,12 +568,13 @@ func (ablInterpolator) scale_RGBA_RGBA(dst *image.RGBA, dr, adr image.Rectangle,
 		if sy < 0 {
 			sy0, sy1 = 0, 0
 			yFrac0, yFrac1 = 0, 1
-		} else if sy1 >= sh {
-			sy1 = sy0
+		} else if sy1 > shMinus1 {
+			sy0, sy1 = shMinus1, shMinus1
 			yFrac0, yFrac1 = 1, 0
 		}
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (float64(dx)+0.5)*xscale - 0.5
 			sx0 := int32(sx)
 			xFrac0 := sx - float64(sx0)
@@ -410,10 +583,11 @@ func (ablInterpolator) scale_RGBA_RGBA(dst *image.RGBA, dr, adr image.Rectangle,
 			if sx < 0 {
 				sx0, sx1 = 0, 0
 				xFrac0, xFrac1 = 0, 1
-			} else if sx1 >= sw {
-				sx1 = sx0
+			} else if sx1 > swMinus1 {
+				sx0, sx1 = swMinus1, swMinus1
 				xFrac0, xFrac1 = 1, 0
 			}
+
 			s00i := src.PixOffset(sr.Min.X+int(sx0), sr.Min.Y+int(sy0))
 			s00ru := uint32(src.Pix[s00i+0]) * 0x101
 			s00gu := uint32(src.Pix[s00i+1]) * 0x101
@@ -466,7 +640,6 @@ func (ablInterpolator) scale_RGBA_RGBA(dst *image.RGBA, dr, adr image.Rectangle,
 			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
 			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
 			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
-			d += 4
 		}
 	}
 }
@@ -476,8 +649,13 @@ func (ablInterpolator) scale_RGBA_Uniform(dst *image.RGBA, dr, adr image.Rectang
 	sh := int32(sr.Dy())
 	yscale := float64(sh) / float64(dr.Dy())
 	xscale := float64(sw) / float64(dr.Dx())
+	swMinus1, shMinus1 := sw-1, sh-1
+
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (float64(dy)+0.5)*yscale - 0.5
+		// If sy < 0, we will clamp sy0 to 0 anyway, so it doesn't matter if
+		// we say int32(sy) instead of int32(math.Floor(sy)). Similarly for
+		// sx, below.
 		sy0 := int32(sy)
 		yFrac0 := sy - float64(sy0)
 		yFrac1 := 1 - yFrac0
@@ -485,12 +663,13 @@ func (ablInterpolator) scale_RGBA_Uniform(dst *image.RGBA, dr, adr image.Rectang
 		if sy < 0 {
 			sy0, sy1 = 0, 0
 			yFrac0, yFrac1 = 0, 1
-		} else if sy1 >= sh {
-			sy1 = sy0
+		} else if sy1 > shMinus1 {
+			sy0, sy1 = shMinus1, shMinus1
 			yFrac0, yFrac1 = 1, 0
 		}
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (float64(dx)+0.5)*xscale - 0.5
 			sx0 := int32(sx)
 			xFrac0 := sx - float64(sx0)
@@ -499,10 +678,11 @@ func (ablInterpolator) scale_RGBA_Uniform(dst *image.RGBA, dr, adr image.Rectang
 			if sx < 0 {
 				sx0, sx1 = 0, 0
 				xFrac0, xFrac1 = 0, 1
-			} else if sx1 >= sw {
-				sx1 = sx0
+			} else if sx1 > swMinus1 {
+				sx0, sx1 = swMinus1, swMinus1
 				xFrac0, xFrac1 = 1, 0
 			}
+
 			s00ru, s00gu, s00bu, s00au := src.At(sr.Min.X+int(sx0), sr.Min.Y+int(sy0)).RGBA()
 			s00r := float64(s00ru)
 			s00g := float64(s00gu)
@@ -539,7 +719,6 @@ func (ablInterpolator) scale_RGBA_Uniform(dst *image.RGBA, dr, adr image.Rectang
 			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
 			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
 			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
-			d += 4
 		}
 	}
 }
@@ -549,8 +728,13 @@ func (ablInterpolator) scale_RGBA_YCbCr(dst *image.RGBA, dr, adr image.Rectangle
 	sh := int32(sr.Dy())
 	yscale := float64(sh) / float64(dr.Dy())
 	xscale := float64(sw) / float64(dr.Dx())
+	swMinus1, shMinus1 := sw-1, sh-1
+
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (float64(dy)+0.5)*yscale - 0.5
+		// If sy < 0, we will clamp sy0 to 0 anyway, so it doesn't matter if
+		// we say int32(sy) instead of int32(math.Floor(sy)). Similarly for
+		// sx, below.
 		sy0 := int32(sy)
 		yFrac0 := sy - float64(sy0)
 		yFrac1 := 1 - yFrac0
@@ -558,12 +742,13 @@ func (ablInterpolator) scale_RGBA_YCbCr(dst *image.RGBA, dr, adr image.Rectangle
 		if sy < 0 {
 			sy0, sy1 = 0, 0
 			yFrac0, yFrac1 = 0, 1
-		} else if sy1 >= sh {
-			sy1 = sy0
+		} else if sy1 > shMinus1 {
+			sy0, sy1 = shMinus1, shMinus1
 			yFrac0, yFrac1 = 1, 0
 		}
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (float64(dx)+0.5)*xscale - 0.5
 			sx0 := int32(sx)
 			xFrac0 := sx - float64(sx0)
@@ -572,10 +757,11 @@ func (ablInterpolator) scale_RGBA_YCbCr(dst *image.RGBA, dr, adr image.Rectangle
 			if sx < 0 {
 				sx0, sx1 = 0, 0
 				xFrac0, xFrac1 = 0, 1
-			} else if sx1 >= sw {
-				sx1 = sx0
+			} else if sx1 > swMinus1 {
+				sx0, sx1 = swMinus1, swMinus1
 				xFrac0, xFrac1 = 1, 0
 			}
+
 			s00ru, s00gu, s00bu, s00au := src.At(sr.Min.X+int(sx0), sr.Min.Y+int(sy0)).RGBA()
 			s00r := float64(s00ru)
 			s00g := float64(s00gu)
@@ -612,7 +798,6 @@ func (ablInterpolator) scale_RGBA_YCbCr(dst *image.RGBA, dr, adr image.Rectangle
 			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
 			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
 			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
-			d += 4
 		}
 	}
 }
@@ -622,8 +807,13 @@ func (ablInterpolator) scale_RGBA_Image(dst *image.RGBA, dr, adr image.Rectangle
 	sh := int32(sr.Dy())
 	yscale := float64(sh) / float64(dr.Dy())
 	xscale := float64(sw) / float64(dr.Dx())
+	swMinus1, shMinus1 := sw-1, sh-1
+
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (float64(dy)+0.5)*yscale - 0.5
+		// If sy < 0, we will clamp sy0 to 0 anyway, so it doesn't matter if
+		// we say int32(sy) instead of int32(math.Floor(sy)). Similarly for
+		// sx, below.
 		sy0 := int32(sy)
 		yFrac0 := sy - float64(sy0)
 		yFrac1 := 1 - yFrac0
@@ -631,12 +821,13 @@ func (ablInterpolator) scale_RGBA_Image(dst *image.RGBA, dr, adr image.Rectangle
 		if sy < 0 {
 			sy0, sy1 = 0, 0
 			yFrac0, yFrac1 = 0, 1
-		} else if sy1 >= sh {
-			sy1 = sy0
+		} else if sy1 > shMinus1 {
+			sy0, sy1 = shMinus1, shMinus1
 			yFrac0, yFrac1 = 1, 0
 		}
 		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
-		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
 			sx := (float64(dx)+0.5)*xscale - 0.5
 			sx0 := int32(sx)
 			xFrac0 := sx - float64(sx0)
@@ -645,10 +836,11 @@ func (ablInterpolator) scale_RGBA_Image(dst *image.RGBA, dr, adr image.Rectangle
 			if sx < 0 {
 				sx0, sx1 = 0, 0
 				xFrac0, xFrac1 = 0, 1
-			} else if sx1 >= sw {
-				sx1 = sx0
+			} else if sx1 > swMinus1 {
+				sx0, sx1 = swMinus1, swMinus1
 				xFrac0, xFrac1 = 1, 0
 			}
+
 			s00ru, s00gu, s00bu, s00au := src.At(sr.Min.X+int(sx0), sr.Min.Y+int(sy0)).RGBA()
 			s00r := float64(s00ru)
 			s00g := float64(s00gu)
@@ -685,7 +877,6 @@ func (ablInterpolator) scale_RGBA_Image(dst *image.RGBA, dr, adr image.Rectangle
 			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
 			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
 			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
-			d += 4
 		}
 	}
 }
@@ -695,10 +886,15 @@ func (ablInterpolator) scale_Image_Image(dst Image, dr, adr image.Rectangle, src
 	sh := int32(sr.Dy())
 	yscale := float64(sh) / float64(dr.Dy())
 	xscale := float64(sw) / float64(dr.Dx())
+	swMinus1, shMinus1 := sw-1, sh-1
 	dstColorRGBA64 := &color.RGBA64{}
 	dstColor := color.Color(dstColorRGBA64)
+
 	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
 		sy := (float64(dy)+0.5)*yscale - 0.5
+		// If sy < 0, we will clamp sy0 to 0 anyway, so it doesn't matter if
+		// we say int32(sy) instead of int32(math.Floor(sy)). Similarly for
+		// sx, below.
 		sy0 := int32(sy)
 		yFrac0 := sy - float64(sy0)
 		yFrac1 := 1 - yFrac0
@@ -706,10 +902,11 @@ func (ablInterpolator) scale_Image_Image(dst Image, dr, adr image.Rectangle, src
 		if sy < 0 {
 			sy0, sy1 = 0, 0
 			yFrac0, yFrac1 = 0, 1
-		} else if sy1 >= sh {
-			sy1 = sy0
+		} else if sy1 > shMinus1 {
+			sy0, sy1 = shMinus1, shMinus1
 			yFrac0, yFrac1 = 1, 0
 		}
+
 		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
 			sx := (float64(dx)+0.5)*xscale - 0.5
 			sx0 := int32(sx)
@@ -719,10 +916,11 @@ func (ablInterpolator) scale_Image_Image(dst Image, dr, adr image.Rectangle, src
 			if sx < 0 {
 				sx0, sx1 = 0, 0
 				xFrac0, xFrac1 = 0, 1
-			} else if sx1 >= sw {
-				sx1 = sx0
+			} else if sx1 > swMinus1 {
+				sx0, sx1 = swMinus1, swMinus1
 				xFrac0, xFrac1 = 1, 0
 			}
+
 			s00ru, s00gu, s00bu, s00au := src.At(sr.Min.X+int(sx0), sr.Min.Y+int(sy0)).RGBA()
 			s00r := float64(s00ru)
 			s00g := float64(s00gu)
@@ -764,6 +962,584 @@ func (ablInterpolator) scale_Image_Image(dst Image, dr, adr image.Rectangle, src
 	}
 }
 
+func (ablInterpolator) transform_RGBA_Gray(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.Gray, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx := d2s[0]*dxf + d2s[1]*dyf + d2s[2]
+			sy := d2s[3]*dxf + d2s[4]*dyf + d2s[5]
+			if !(image.Point{int(math.Floor(sx)), int(math.Floor(sy))}).In(sr) {
+				continue
+			}
+
+			sx -= 0.5
+			sxf := math.Floor(sx)
+			xFrac0 := sx - sxf
+			xFrac1 := 1 - xFrac0
+			sx0 := int(sxf)
+			sx1 := sx0 + 1
+			if sx0 < sr.Min.X {
+				sx0, sx1 = sr.Min.X, sr.Min.X
+				xFrac0, xFrac1 = 0, 1
+			} else if sx1 >= sr.Max.X {
+				sx0, sx1 = sr.Max.X-1, sr.Max.X-1
+				xFrac0, xFrac1 = 1, 0
+			}
+
+			sy -= 0.5
+			syf := math.Floor(sy)
+			yFrac0 := sy - syf
+			yFrac1 := 1 - yFrac0
+			sy0 := int(syf)
+			sy1 := sy0 + 1
+			if sy0 < sr.Min.Y {
+				sy0, sy1 = sr.Min.Y, sr.Min.Y
+				yFrac0, yFrac1 = 0, 1
+			} else if sy1 >= sr.Max.Y {
+				sy0, sy1 = sr.Max.Y-1, sr.Max.Y-1
+				yFrac0, yFrac1 = 1, 0
+			}
+
+			s00ru, s00gu, s00bu, s00au := src.At(sx0, sy0).RGBA()
+			s00r := float64(s00ru)
+			s00g := float64(s00gu)
+			s00b := float64(s00bu)
+			s00a := float64(s00au)
+			s10ru, s10gu, s10bu, s10au := src.At(sx1, sy0).RGBA()
+			s10r := float64(s10ru)
+			s10g := float64(s10gu)
+			s10b := float64(s10bu)
+			s10a := float64(s10au)
+			s10r = xFrac1*s00r + xFrac0*s10r
+			s10g = xFrac1*s00g + xFrac0*s10g
+			s10b = xFrac1*s00b + xFrac0*s10b
+			s10a = xFrac1*s00a + xFrac0*s10a
+			s01ru, s01gu, s01bu, s01au := src.At(sx0, sy1).RGBA()
+			s01r := float64(s01ru)
+			s01g := float64(s01gu)
+			s01b := float64(s01bu)
+			s01a := float64(s01au)
+			s11ru, s11gu, s11bu, s11au := src.At(sx1, sy1).RGBA()
+			s11r := float64(s11ru)
+			s11g := float64(s11gu)
+			s11b := float64(s11bu)
+			s11a := float64(s11au)
+			s11r = xFrac1*s01r + xFrac0*s11r
+			s11g = xFrac1*s01g + xFrac0*s11g
+			s11b = xFrac1*s01b + xFrac0*s11b
+			s11a = xFrac1*s01a + xFrac0*s11a
+			s11r = yFrac1*s10r + yFrac0*s11r
+			s11g = yFrac1*s10g + yFrac0*s11g
+			s11b = yFrac1*s10b + yFrac0*s11b
+			s11a = yFrac1*s10a + yFrac0*s11a
+			dst.Pix[d+0] = uint8(uint32(s11r) >> 8)
+			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
+			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
+			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
+		}
+	}
+}
+
+func (ablInterpolator) transform_RGBA_NRGBA(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.NRGBA, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx := d2s[0]*dxf + d2s[1]*dyf + d2s[2]
+			sy := d2s[3]*dxf + d2s[4]*dyf + d2s[5]
+			if !(image.Point{int(math.Floor(sx)), int(math.Floor(sy))}).In(sr) {
+				continue
+			}
+
+			sx -= 0.5
+			sxf := math.Floor(sx)
+			xFrac0 := sx - sxf
+			xFrac1 := 1 - xFrac0
+			sx0 := int(sxf)
+			sx1 := sx0 + 1
+			if sx0 < sr.Min.X {
+				sx0, sx1 = sr.Min.X, sr.Min.X
+				xFrac0, xFrac1 = 0, 1
+			} else if sx1 >= sr.Max.X {
+				sx0, sx1 = sr.Max.X-1, sr.Max.X-1
+				xFrac0, xFrac1 = 1, 0
+			}
+
+			sy -= 0.5
+			syf := math.Floor(sy)
+			yFrac0 := sy - syf
+			yFrac1 := 1 - yFrac0
+			sy0 := int(syf)
+			sy1 := sy0 + 1
+			if sy0 < sr.Min.Y {
+				sy0, sy1 = sr.Min.Y, sr.Min.Y
+				yFrac0, yFrac1 = 0, 1
+			} else if sy1 >= sr.Max.Y {
+				sy0, sy1 = sr.Max.Y-1, sr.Max.Y-1
+				yFrac0, yFrac1 = 1, 0
+			}
+
+			s00ru, s00gu, s00bu, s00au := src.At(sx0, sy0).RGBA()
+			s00r := float64(s00ru)
+			s00g := float64(s00gu)
+			s00b := float64(s00bu)
+			s00a := float64(s00au)
+			s10ru, s10gu, s10bu, s10au := src.At(sx1, sy0).RGBA()
+			s10r := float64(s10ru)
+			s10g := float64(s10gu)
+			s10b := float64(s10bu)
+			s10a := float64(s10au)
+			s10r = xFrac1*s00r + xFrac0*s10r
+			s10g = xFrac1*s00g + xFrac0*s10g
+			s10b = xFrac1*s00b + xFrac0*s10b
+			s10a = xFrac1*s00a + xFrac0*s10a
+			s01ru, s01gu, s01bu, s01au := src.At(sx0, sy1).RGBA()
+			s01r := float64(s01ru)
+			s01g := float64(s01gu)
+			s01b := float64(s01bu)
+			s01a := float64(s01au)
+			s11ru, s11gu, s11bu, s11au := src.At(sx1, sy1).RGBA()
+			s11r := float64(s11ru)
+			s11g := float64(s11gu)
+			s11b := float64(s11bu)
+			s11a := float64(s11au)
+			s11r = xFrac1*s01r + xFrac0*s11r
+			s11g = xFrac1*s01g + xFrac0*s11g
+			s11b = xFrac1*s01b + xFrac0*s11b
+			s11a = xFrac1*s01a + xFrac0*s11a
+			s11r = yFrac1*s10r + yFrac0*s11r
+			s11g = yFrac1*s10g + yFrac0*s11g
+			s11b = yFrac1*s10b + yFrac0*s11b
+			s11a = yFrac1*s10a + yFrac0*s11a
+			dst.Pix[d+0] = uint8(uint32(s11r) >> 8)
+			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
+			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
+			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
+		}
+	}
+}
+
+func (ablInterpolator) transform_RGBA_RGBA(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.RGBA, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx := d2s[0]*dxf + d2s[1]*dyf + d2s[2]
+			sy := d2s[3]*dxf + d2s[4]*dyf + d2s[5]
+			if !(image.Point{int(math.Floor(sx)), int(math.Floor(sy))}).In(sr) {
+				continue
+			}
+
+			sx -= 0.5
+			sxf := math.Floor(sx)
+			xFrac0 := sx - sxf
+			xFrac1 := 1 - xFrac0
+			sx0 := int(sxf)
+			sx1 := sx0 + 1
+			if sx0 < sr.Min.X {
+				sx0, sx1 = sr.Min.X, sr.Min.X
+				xFrac0, xFrac1 = 0, 1
+			} else if sx1 >= sr.Max.X {
+				sx0, sx1 = sr.Max.X-1, sr.Max.X-1
+				xFrac0, xFrac1 = 1, 0
+			}
+
+			sy -= 0.5
+			syf := math.Floor(sy)
+			yFrac0 := sy - syf
+			yFrac1 := 1 - yFrac0
+			sy0 := int(syf)
+			sy1 := sy0 + 1
+			if sy0 < sr.Min.Y {
+				sy0, sy1 = sr.Min.Y, sr.Min.Y
+				yFrac0, yFrac1 = 0, 1
+			} else if sy1 >= sr.Max.Y {
+				sy0, sy1 = sr.Max.Y-1, sr.Max.Y-1
+				yFrac0, yFrac1 = 1, 0
+			}
+
+			s00i := src.PixOffset(sx0, sy0)
+			s00ru := uint32(src.Pix[s00i+0]) * 0x101
+			s00gu := uint32(src.Pix[s00i+1]) * 0x101
+			s00bu := uint32(src.Pix[s00i+2]) * 0x101
+			s00au := uint32(src.Pix[s00i+3]) * 0x101
+			s00r := float64(s00ru)
+			s00g := float64(s00gu)
+			s00b := float64(s00bu)
+			s00a := float64(s00au)
+			s10i := src.PixOffset(sx1, sy0)
+			s10ru := uint32(src.Pix[s10i+0]) * 0x101
+			s10gu := uint32(src.Pix[s10i+1]) * 0x101
+			s10bu := uint32(src.Pix[s10i+2]) * 0x101
+			s10au := uint32(src.Pix[s10i+3]) * 0x101
+			s10r := float64(s10ru)
+			s10g := float64(s10gu)
+			s10b := float64(s10bu)
+			s10a := float64(s10au)
+			s10r = xFrac1*s00r + xFrac0*s10r
+			s10g = xFrac1*s00g + xFrac0*s10g
+			s10b = xFrac1*s00b + xFrac0*s10b
+			s10a = xFrac1*s00a + xFrac0*s10a
+			s01i := src.PixOffset(sx0, sy1)
+			s01ru := uint32(src.Pix[s01i+0]) * 0x101
+			s01gu := uint32(src.Pix[s01i+1]) * 0x101
+			s01bu := uint32(src.Pix[s01i+2]) * 0x101
+			s01au := uint32(src.Pix[s01i+3]) * 0x101
+			s01r := float64(s01ru)
+			s01g := float64(s01gu)
+			s01b := float64(s01bu)
+			s01a := float64(s01au)
+			s11i := src.PixOffset(sx1, sy1)
+			s11ru := uint32(src.Pix[s11i+0]) * 0x101
+			s11gu := uint32(src.Pix[s11i+1]) * 0x101
+			s11bu := uint32(src.Pix[s11i+2]) * 0x101
+			s11au := uint32(src.Pix[s11i+3]) * 0x101
+			s11r := float64(s11ru)
+			s11g := float64(s11gu)
+			s11b := float64(s11bu)
+			s11a := float64(s11au)
+			s11r = xFrac1*s01r + xFrac0*s11r
+			s11g = xFrac1*s01g + xFrac0*s11g
+			s11b = xFrac1*s01b + xFrac0*s11b
+			s11a = xFrac1*s01a + xFrac0*s11a
+			s11r = yFrac1*s10r + yFrac0*s11r
+			s11g = yFrac1*s10g + yFrac0*s11g
+			s11b = yFrac1*s10b + yFrac0*s11b
+			s11a = yFrac1*s10a + yFrac0*s11a
+			dst.Pix[d+0] = uint8(uint32(s11r) >> 8)
+			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
+			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
+			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
+		}
+	}
+}
+
+func (ablInterpolator) transform_RGBA_Uniform(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.Uniform, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx := d2s[0]*dxf + d2s[1]*dyf + d2s[2]
+			sy := d2s[3]*dxf + d2s[4]*dyf + d2s[5]
+			if !(image.Point{int(math.Floor(sx)), int(math.Floor(sy))}).In(sr) {
+				continue
+			}
+
+			sx -= 0.5
+			sxf := math.Floor(sx)
+			xFrac0 := sx - sxf
+			xFrac1 := 1 - xFrac0
+			sx0 := int(sxf)
+			sx1 := sx0 + 1
+			if sx0 < sr.Min.X {
+				sx0, sx1 = sr.Min.X, sr.Min.X
+				xFrac0, xFrac1 = 0, 1
+			} else if sx1 >= sr.Max.X {
+				sx0, sx1 = sr.Max.X-1, sr.Max.X-1
+				xFrac0, xFrac1 = 1, 0
+			}
+
+			sy -= 0.5
+			syf := math.Floor(sy)
+			yFrac0 := sy - syf
+			yFrac1 := 1 - yFrac0
+			sy0 := int(syf)
+			sy1 := sy0 + 1
+			if sy0 < sr.Min.Y {
+				sy0, sy1 = sr.Min.Y, sr.Min.Y
+				yFrac0, yFrac1 = 0, 1
+			} else if sy1 >= sr.Max.Y {
+				sy0, sy1 = sr.Max.Y-1, sr.Max.Y-1
+				yFrac0, yFrac1 = 1, 0
+			}
+
+			s00ru, s00gu, s00bu, s00au := src.At(sx0, sy0).RGBA()
+			s00r := float64(s00ru)
+			s00g := float64(s00gu)
+			s00b := float64(s00bu)
+			s00a := float64(s00au)
+			s10ru, s10gu, s10bu, s10au := src.At(sx1, sy0).RGBA()
+			s10r := float64(s10ru)
+			s10g := float64(s10gu)
+			s10b := float64(s10bu)
+			s10a := float64(s10au)
+			s10r = xFrac1*s00r + xFrac0*s10r
+			s10g = xFrac1*s00g + xFrac0*s10g
+			s10b = xFrac1*s00b + xFrac0*s10b
+			s10a = xFrac1*s00a + xFrac0*s10a
+			s01ru, s01gu, s01bu, s01au := src.At(sx0, sy1).RGBA()
+			s01r := float64(s01ru)
+			s01g := float64(s01gu)
+			s01b := float64(s01bu)
+			s01a := float64(s01au)
+			s11ru, s11gu, s11bu, s11au := src.At(sx1, sy1).RGBA()
+			s11r := float64(s11ru)
+			s11g := float64(s11gu)
+			s11b := float64(s11bu)
+			s11a := float64(s11au)
+			s11r = xFrac1*s01r + xFrac0*s11r
+			s11g = xFrac1*s01g + xFrac0*s11g
+			s11b = xFrac1*s01b + xFrac0*s11b
+			s11a = xFrac1*s01a + xFrac0*s11a
+			s11r = yFrac1*s10r + yFrac0*s11r
+			s11g = yFrac1*s10g + yFrac0*s11g
+			s11b = yFrac1*s10b + yFrac0*s11b
+			s11a = yFrac1*s10a + yFrac0*s11a
+			dst.Pix[d+0] = uint8(uint32(s11r) >> 8)
+			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
+			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
+			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
+		}
+	}
+}
+
+func (ablInterpolator) transform_RGBA_YCbCr(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src *image.YCbCr, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx := d2s[0]*dxf + d2s[1]*dyf + d2s[2]
+			sy := d2s[3]*dxf + d2s[4]*dyf + d2s[5]
+			if !(image.Point{int(math.Floor(sx)), int(math.Floor(sy))}).In(sr) {
+				continue
+			}
+
+			sx -= 0.5
+			sxf := math.Floor(sx)
+			xFrac0 := sx - sxf
+			xFrac1 := 1 - xFrac0
+			sx0 := int(sxf)
+			sx1 := sx0 + 1
+			if sx0 < sr.Min.X {
+				sx0, sx1 = sr.Min.X, sr.Min.X
+				xFrac0, xFrac1 = 0, 1
+			} else if sx1 >= sr.Max.X {
+				sx0, sx1 = sr.Max.X-1, sr.Max.X-1
+				xFrac0, xFrac1 = 1, 0
+			}
+
+			sy -= 0.5
+			syf := math.Floor(sy)
+			yFrac0 := sy - syf
+			yFrac1 := 1 - yFrac0
+			sy0 := int(syf)
+			sy1 := sy0 + 1
+			if sy0 < sr.Min.Y {
+				sy0, sy1 = sr.Min.Y, sr.Min.Y
+				yFrac0, yFrac1 = 0, 1
+			} else if sy1 >= sr.Max.Y {
+				sy0, sy1 = sr.Max.Y-1, sr.Max.Y-1
+				yFrac0, yFrac1 = 1, 0
+			}
+
+			s00ru, s00gu, s00bu, s00au := src.At(sx0, sy0).RGBA()
+			s00r := float64(s00ru)
+			s00g := float64(s00gu)
+			s00b := float64(s00bu)
+			s00a := float64(s00au)
+			s10ru, s10gu, s10bu, s10au := src.At(sx1, sy0).RGBA()
+			s10r := float64(s10ru)
+			s10g := float64(s10gu)
+			s10b := float64(s10bu)
+			s10a := float64(s10au)
+			s10r = xFrac1*s00r + xFrac0*s10r
+			s10g = xFrac1*s00g + xFrac0*s10g
+			s10b = xFrac1*s00b + xFrac0*s10b
+			s10a = xFrac1*s00a + xFrac0*s10a
+			s01ru, s01gu, s01bu, s01au := src.At(sx0, sy1).RGBA()
+			s01r := float64(s01ru)
+			s01g := float64(s01gu)
+			s01b := float64(s01bu)
+			s01a := float64(s01au)
+			s11ru, s11gu, s11bu, s11au := src.At(sx1, sy1).RGBA()
+			s11r := float64(s11ru)
+			s11g := float64(s11gu)
+			s11b := float64(s11bu)
+			s11a := float64(s11au)
+			s11r = xFrac1*s01r + xFrac0*s11r
+			s11g = xFrac1*s01g + xFrac0*s11g
+			s11b = xFrac1*s01b + xFrac0*s11b
+			s11a = xFrac1*s01a + xFrac0*s11a
+			s11r = yFrac1*s10r + yFrac0*s11r
+			s11g = yFrac1*s10g + yFrac0*s11g
+			s11b = yFrac1*s10b + yFrac0*s11b
+			s11a = yFrac1*s10a + yFrac0*s11a
+			dst.Pix[d+0] = uint8(uint32(s11r) >> 8)
+			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
+			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
+			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
+		}
+	}
+}
+
+func (ablInterpolator) transform_RGBA_Image(dst *image.RGBA, dr, adr image.Rectangle, d2s *f64.Aff3, src image.Image, sr image.Rectangle) {
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		d := dst.PixOffset(dr.Min.X+adr.Min.X, dr.Min.Y+int(dy))
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx, d = dx+1, d+4 {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx := d2s[0]*dxf + d2s[1]*dyf + d2s[2]
+			sy := d2s[3]*dxf + d2s[4]*dyf + d2s[5]
+			if !(image.Point{int(math.Floor(sx)), int(math.Floor(sy))}).In(sr) {
+				continue
+			}
+
+			sx -= 0.5
+			sxf := math.Floor(sx)
+			xFrac0 := sx - sxf
+			xFrac1 := 1 - xFrac0
+			sx0 := int(sxf)
+			sx1 := sx0 + 1
+			if sx0 < sr.Min.X {
+				sx0, sx1 = sr.Min.X, sr.Min.X
+				xFrac0, xFrac1 = 0, 1
+			} else if sx1 >= sr.Max.X {
+				sx0, sx1 = sr.Max.X-1, sr.Max.X-1
+				xFrac0, xFrac1 = 1, 0
+			}
+
+			sy -= 0.5
+			syf := math.Floor(sy)
+			yFrac0 := sy - syf
+			yFrac1 := 1 - yFrac0
+			sy0 := int(syf)
+			sy1 := sy0 + 1
+			if sy0 < sr.Min.Y {
+				sy0, sy1 = sr.Min.Y, sr.Min.Y
+				yFrac0, yFrac1 = 0, 1
+			} else if sy1 >= sr.Max.Y {
+				sy0, sy1 = sr.Max.Y-1, sr.Max.Y-1
+				yFrac0, yFrac1 = 1, 0
+			}
+
+			s00ru, s00gu, s00bu, s00au := src.At(sx0, sy0).RGBA()
+			s00r := float64(s00ru)
+			s00g := float64(s00gu)
+			s00b := float64(s00bu)
+			s00a := float64(s00au)
+			s10ru, s10gu, s10bu, s10au := src.At(sx1, sy0).RGBA()
+			s10r := float64(s10ru)
+			s10g := float64(s10gu)
+			s10b := float64(s10bu)
+			s10a := float64(s10au)
+			s10r = xFrac1*s00r + xFrac0*s10r
+			s10g = xFrac1*s00g + xFrac0*s10g
+			s10b = xFrac1*s00b + xFrac0*s10b
+			s10a = xFrac1*s00a + xFrac0*s10a
+			s01ru, s01gu, s01bu, s01au := src.At(sx0, sy1).RGBA()
+			s01r := float64(s01ru)
+			s01g := float64(s01gu)
+			s01b := float64(s01bu)
+			s01a := float64(s01au)
+			s11ru, s11gu, s11bu, s11au := src.At(sx1, sy1).RGBA()
+			s11r := float64(s11ru)
+			s11g := float64(s11gu)
+			s11b := float64(s11bu)
+			s11a := float64(s11au)
+			s11r = xFrac1*s01r + xFrac0*s11r
+			s11g = xFrac1*s01g + xFrac0*s11g
+			s11b = xFrac1*s01b + xFrac0*s11b
+			s11a = xFrac1*s01a + xFrac0*s11a
+			s11r = yFrac1*s10r + yFrac0*s11r
+			s11g = yFrac1*s10g + yFrac0*s11g
+			s11b = yFrac1*s10b + yFrac0*s11b
+			s11a = yFrac1*s10a + yFrac0*s11a
+			dst.Pix[d+0] = uint8(uint32(s11r) >> 8)
+			dst.Pix[d+1] = uint8(uint32(s11g) >> 8)
+			dst.Pix[d+2] = uint8(uint32(s11b) >> 8)
+			dst.Pix[d+3] = uint8(uint32(s11a) >> 8)
+		}
+	}
+}
+
+func (ablInterpolator) transform_Image_Image(dst Image, dr, adr image.Rectangle, d2s *f64.Aff3, src image.Image, sr image.Rectangle) {
+	dstColorRGBA64 := &color.RGBA64{}
+	dstColor := color.Color(dstColorRGBA64)
+	for dy := int32(adr.Min.Y); dy < int32(adr.Max.Y); dy++ {
+		dyf := float64(dr.Min.Y+int(dy)) + 0.5
+		for dx := int32(adr.Min.X); dx < int32(adr.Max.X); dx++ {
+			dxf := float64(dr.Min.X+int(dx)) + 0.5
+			sx := d2s[0]*dxf + d2s[1]*dyf + d2s[2]
+			sy := d2s[3]*dxf + d2s[4]*dyf + d2s[5]
+			if !(image.Point{int(math.Floor(sx)), int(math.Floor(sy))}).In(sr) {
+				continue
+			}
+
+			sx -= 0.5
+			sxf := math.Floor(sx)
+			xFrac0 := sx - sxf
+			xFrac1 := 1 - xFrac0
+			sx0 := int(sxf)
+			sx1 := sx0 + 1
+			if sx0 < sr.Min.X {
+				sx0, sx1 = sr.Min.X, sr.Min.X
+				xFrac0, xFrac1 = 0, 1
+			} else if sx1 >= sr.Max.X {
+				sx0, sx1 = sr.Max.X-1, sr.Max.X-1
+				xFrac0, xFrac1 = 1, 0
+			}
+
+			sy -= 0.5
+			syf := math.Floor(sy)
+			yFrac0 := sy - syf
+			yFrac1 := 1 - yFrac0
+			sy0 := int(syf)
+			sy1 := sy0 + 1
+			if sy0 < sr.Min.Y {
+				sy0, sy1 = sr.Min.Y, sr.Min.Y
+				yFrac0, yFrac1 = 0, 1
+			} else if sy1 >= sr.Max.Y {
+				sy0, sy1 = sr.Max.Y-1, sr.Max.Y-1
+				yFrac0, yFrac1 = 1, 0
+			}
+
+			s00ru, s00gu, s00bu, s00au := src.At(sx0, sy0).RGBA()
+			s00r := float64(s00ru)
+			s00g := float64(s00gu)
+			s00b := float64(s00bu)
+			s00a := float64(s00au)
+			s10ru, s10gu, s10bu, s10au := src.At(sx1, sy0).RGBA()
+			s10r := float64(s10ru)
+			s10g := float64(s10gu)
+			s10b := float64(s10bu)
+			s10a := float64(s10au)
+			s10r = xFrac1*s00r + xFrac0*s10r
+			s10g = xFrac1*s00g + xFrac0*s10g
+			s10b = xFrac1*s00b + xFrac0*s10b
+			s10a = xFrac1*s00a + xFrac0*s10a
+			s01ru, s01gu, s01bu, s01au := src.At(sx0, sy1).RGBA()
+			s01r := float64(s01ru)
+			s01g := float64(s01gu)
+			s01b := float64(s01bu)
+			s01a := float64(s01au)
+			s11ru, s11gu, s11bu, s11au := src.At(sx1, sy1).RGBA()
+			s11r := float64(s11ru)
+			s11g := float64(s11gu)
+			s11b := float64(s11bu)
+			s11a := float64(s11au)
+			s11r = xFrac1*s01r + xFrac0*s11r
+			s11g = xFrac1*s01g + xFrac0*s11g
+			s11b = xFrac1*s01b + xFrac0*s11b
+			s11a = xFrac1*s01a + xFrac0*s11a
+			s11r = yFrac1*s10r + yFrac0*s11r
+			s11g = yFrac1*s10g + yFrac0*s11g
+			s11b = yFrac1*s10b + yFrac0*s11b
+			s11a = yFrac1*s10a + yFrac0*s11a
+			dstColorRGBA64.R = uint16(s11r)
+			dstColorRGBA64.G = uint16(s11g)
+			dstColorRGBA64.B = uint16(s11b)
+			dstColorRGBA64.A = uint16(s11a)
+			dst.Set(dr.Min.X+int(dx), dr.Min.Y+int(dy), dstColor)
+		}
+	}
+}
+
 func (z *kernelScaler) Scale(dst Image, dr image.Rectangle, src image.Image, sr image.Rectangle, opts *Options) {
 	if z.dw != int32(dr.Dx()) || z.dh != int32(dr.Dy()) || z.sw != int32(sr.Dx()) || z.sh != int32(sr.Dy()) {
 		z.kernel.Scale(dst, dr, src, sr, opts)
diff --git a/draw/scale.go b/draw/scale.go
index 8fc1644..4f94086 100644
--- a/draw/scale.go
+++ b/draw/scale.go
@@ -249,3 +249,72 @@ func ftou(f float64) uint16 {
 	}
 	return 0
 }
+
+// invert returns the inverse of m.
+//
+// TODO: move this into the f64 package, once we work out the convention for
+// matrix methods in that package: do they modify the receiver, take a dst
+// pointer argument, or return a new value?
+func invert(m *f64.Aff3) f64.Aff3 {
+	m00 := +m[3*1+1]
+	m01 := -m[3*0+1]
+	m02 := +m[3*1+2]*m[3*0+1] - m[3*1+1]*m[3*0+2]
+	m10 := -m[3*1+0]
+	m11 := +m[3*0+0]
+	m12 := +m[3*1+0]*m[3*0+2] - m[3*1+2]*m[3*0+0]
+
+	det := m00*m11 - m10*m01
+
+	return f64.Aff3{
+		m00 / det,
+		m01 / det,
+		m02 / det,
+		m10 / det,
+		m11 / det,
+		m12 / det,
+	}
+}
+
+// transformRect returns a rectangle dr that contains sr transformed by s2d.
+func transformRect(s2d *f64.Aff3, sr *image.Rectangle) (dr image.Rectangle) {
+	ps := [...]image.Point{
+		{sr.Min.X, sr.Min.Y},
+		{sr.Max.X, sr.Min.Y},
+		{sr.Min.X, sr.Max.Y},
+		{sr.Max.X, sr.Max.Y},
+	}
+	for i, p := range ps {
+		sxf := float64(p.X)
+		syf := float64(p.Y)
+		dx := int(math.Floor(s2d[0]*sxf + s2d[1]*syf + s2d[2]))
+		dy := int(math.Floor(s2d[3]*sxf + s2d[4]*syf + s2d[5]))
+
+		// The +1 adjustments below are because an image.Rectangle is inclusive
+		// on the low end but exclusive on the high end.
+
+		if i == 0 {
+			dr = image.Rectangle{
+				Min: image.Point{dx + 0, dy + 0},
+				Max: image.Point{dx + 1, dy + 1},
+			}
+			continue
+		}
+
+		if dr.Min.X > dx {
+			dr.Min.X = dx
+		}
+		dx++
+		if dr.Max.X < dx {
+			dr.Max.X = dx
+		}
+
+		if dr.Min.Y > dy {
+			dr.Min.Y = dy
+		}
+		dy++
+		if dr.Max.Y < dy {
+			dr.Max.Y = dy
+		}
+	}
+	return dr
+}
diff --git a/draw/scale_test.go b/draw/scale_test.go
index 30a82c6..c15810e 100644
--- a/draw/scale_test.go
+++ b/draw/scale_test.go
@@ -16,18 +16,28 @@ import (
 	"reflect"
 	"testing"
 
+	"golang.org/x/image/math/f64"
+
 	_ "image/jpeg"
 )
 
-var genScaleFiles = flag.Bool("gen_scale_files", false, "whether to generate the TestScaleXxx golden files.")
+var genGoldenFiles = flag.Bool("gen_golden_files", false, "whether to generate the TestXxx golden files.")
 
-// testScale tests that scaling the source image gives the exact destination
-// image. This is to ensure that any refactoring or optimization of the scaling
-// code doesn't change the scaling behavior. Changing the actual algorithm or
-// kernel used by any particular quality setting will obviously change the
-// resultant pixels. In such a case, use the gen_scale_files flag to regenerate
-// the golden files.
-func testScale(t *testing.T, w int, h int, direction, srcFilename string) {
+var transformMatrix = func() *f64.Aff3 {
+	const scale, cos30, sin30 = 3.75, 0.866025404, 0.5
+	return &f64.Aff3{
+		+scale * cos30, -scale * sin30, 40,
+		+scale * sin30, +scale * cos30, 10,
+	}
+}()
+
+// testInterp tests that interpolating the source image gives the exact
+// destination image. This is to ensure that any refactoring or optimization of
+// the interpolation code doesn't change the behavior. Changing the actual
+// algorithm or kernel used by any particular quality setting will obviously
+// change the resultant pixels. In such a case, use the gen_golden_files flag
+// to regenerate the golden files.
+func testInterp(t *testing.T, w int, h int, direction, srcFilename string) {
 	f, err := os.Open("../testdata/go-turns-two-" + srcFilename)
 	if err != nil {
 		t.Fatalf("Open: %v", err)
@@ -44,12 +54,21 @@ func testScale(t *testing.T, w int, h int, direction, srcFilename string) {
 		"cr": CatmullRom,
 	}
 	for name, q := range testCases {
-		gotFilename := fmt.Sprintf("../testdata/go-turns-two-%s-%s.png", direction, name)
+		goldenFilename := fmt.Sprintf("../testdata/go-turns-two-%s-%s.png", direction, name)
 
 		got := image.NewRGBA(image.Rect(0, 0, w, h))
-		q.Scale(got, got.Bounds(), src, src.Bounds(), nil)
-		if *genScaleFiles {
-			g, err := os.Create(gotFilename)
+		if direction == "rotate" {
+			if name == "bl" || name == "cr" {
+				// TODO: implement Kernel.Transform.
+				continue
+			}
+			q.Transform(got, transformMatrix, src, src.Bounds(), nil)
+		} else {
+			q.Scale(got, got.Bounds(), src, src.Bounds(), nil)
+		}
+
+		if *genGoldenFiles {
+			g, err := os.Create(goldenFilename)
 			if err != nil {
 				t.Errorf("Create: %v", err)
 				continue
@@ -62,27 +81,35 @@ func testScale(t *testing.T, w int, h int, direction, srcFilename string) {
 			continue
 		}
 
-		g, err := os.Open(gotFilename)
+		g, err := os.Open(goldenFilename)
 		if err != nil {
 			t.Errorf("Open: %v", err)
 			continue
 		}
 		defer g.Close()
-		want, err := png.Decode(g)
+		wantRaw, err := png.Decode(g)
 		if err != nil {
 			t.Errorf("Decode: %v", err)
 			continue
 		}
+		// convert wantRaw to RGBA.
+		want, ok := wantRaw.(*image.RGBA)
+		if !ok {
+			b := wantRaw.Bounds()
+			want = image.NewRGBA(b)
+			Draw(want, b, wantRaw, b.Min, Src)
+		}
 
 		if !reflect.DeepEqual(got, want) {
-			t.Errorf("%s: actual image differs from golden image", gotFilename)
+			t.Errorf("%s: actual image differs from golden image", goldenFilename)
 			continue
 		}
 	}
 }
 
-func TestScaleDown(t *testing.T) { testScale(t, 100, 100, "down", "280x360.jpeg") }
-func TestScaleUp(t *testing.T)   { testScale(t, 75, 100, "up", "14x18.png") }
+func TestScaleDown(t *testing.T) { testInterp(t, 100, 100, "down", "280x360.jpeg") }
+func TestScaleUp(t *testing.T)   { testInterp(t, 75, 100, "up", "14x18.png") }
+func TestTransform(t *testing.T) { testInterp(t, 100, 100, "rotate", "14x18.png") }
 
 func fillPix(r *rand.Rand, pixs ...[]byte) {
 	for _, pix := range pixs {
@@ -92,7 +119,7 @@ func fillPix(r *rand.Rand, pixs ...[]byte) {
 	}
 }
 
-func TestScaleClipCommute(t *testing.T) {
+func TestInterpClipCommute(t *testing.T) {
 	src := image.NewNRGBA(image.Rect(0, 0, 20, 20))
 	fillPix(rand.New(rand.NewSource(0)), src.Pix)
 
@@ -103,28 +130,46 @@ func TestScaleClipCommute(t *testing.T) {
 		ApproxBiLinear,
 		CatmullRom,
 	}
-	for _, q := range qs {
-		dst0 := image.NewRGBA(image.Rect(1, 1, 10, 10))
-		dst1 := image.NewRGBA(image.Rect(1, 1, 10, 10))
-		for i := range dst0.Pix {
-			dst0.Pix[i] = uint8(i / 4)
-			dst1.Pix[i] = uint8(i / 4)
-		}
+	for _, transform := range []bool{false, true} {
+		for _, q := range qs {
+			if transform && q == CatmullRom {
+				// TODO: implement Kernel.Transform.
+				continue
+			}
 
-		// Scale then clip.
-		q.Scale(dst0, outer, src, src.Bounds(), nil)
-		dst0 = dst0.SubImage(inner).(*image.RGBA)
+			dst0 := image.NewRGBA(image.Rect(1, 1, 10, 10))
+			dst1 := image.NewRGBA(image.Rect(1, 1, 10, 10))
+			for i := range dst0.Pix {
+				dst0.Pix[i] = uint8(i / 4)
+				dst1.Pix[i] = uint8(i / 4)
+			}
 
-		// Clip then scale.
-		dst1 = dst1.SubImage(inner).(*image.RGBA)
-		q.Scale(dst1, outer, src, src.Bounds(), nil)
+			var interp func(dst *image.RGBA)
+			if transform {
+				interp = func(dst *image.RGBA) {
+					q.Transform(dst, transformMatrix, src, src.Bounds(), nil)
+				}
+			} else {
+				interp = func(dst *image.RGBA) {
+					q.Scale(dst, outer, src, src.Bounds(), nil)
+				}
+			}
 
-	loop:
-		for y := inner.Min.Y; y < inner.Max.Y; y++ {
-			for x := inner.Min.X; x < inner.Max.X; x++ {
-				if c0, c1 := dst0.RGBAAt(x, y), dst1.RGBAAt(x, y); c0 != c1 {
-					t.Errorf("q=%T: at (%d, %d): c0=%v, c1=%v", q, x, y, c0, c1)
-					break loop
+			// Interpolate then clip.
+			interp(dst0)
+			dst0 = dst0.SubImage(inner).(*image.RGBA)
+
+			// Clip then interpolate.
+			dst1 = dst1.SubImage(inner).(*image.RGBA)
+			interp(dst1)
+
+		loop:
+			for y := inner.Min.Y; y < inner.Max.Y; y++ {
+				for x := inner.Min.X; x < inner.Max.X; x++ {
+					if c0, c1 := dst0.RGBAAt(x, y), dst1.RGBAAt(x, y); c0 != c1 {
+						t.Errorf("q=%T: at (%d, %d): c0=%v, c1=%v", q, x, y, c0, c1)
+						break loop
+					}
 				}
 			}
 		}
@@ -184,7 +229,7 @@ func TestSrcTranslationInvariance(t *testing.T) {
 				t.Errorf("pix differ for delta=%v, q=%T", delta, q)
 			}
 
-			// TODO: Transform.
+			// TODO: Transform, once Kernel.Transform is implemented.
 		}
 	}
 }
@@ -250,6 +295,8 @@ func TestFastPaths(t *testing.T) {
 					if !bytes.Equal(dst0.Pix, dst1.Pix) {
 						t.Errorf("pix differ for dr=%v, src=%T, sr=%v, q=%T", dr, src, sr, q)
 					}
+
+					// TODO: Transform, once Kernel.Transform is implemented.
 				}
 			}
 		}
@@ -331,6 +378,20 @@ func benchScale(b *testing.B, srcf func(image.Rectangle) (image.Image, error), w
 	}
 }
 
+func benchTform(b *testing.B, srcf func(image.Rectangle) (image.Image, error), w int, h int, q Interpolator) {
+	dst := image.NewRGBA(image.Rect(0, 0, w, h))
+	src, err := srcf(image.Rect(0, 0, 1024, 768))
+	if err != nil {
+		b.Fatal(err)
+	}
+	sr := src.Bounds()
+
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		q.Transform(dst, transformMatrix, src, sr, nil)
+	}
+}
+
 func BenchmarkScaleLargeDownNN(b *testing.B) { benchScale(b, srcYCbCrLarge, 200, 150, NearestNeighbor) }
 func BenchmarkScaleLargeDownAB(b *testing.B) { benchScale(b, srcYCbCrLarge, 200, 150, ApproxBiLinear) }
 func BenchmarkScaleLargeDownBL(b *testing.B) { benchScale(b, srcYCbCrLarge, 200, 150, BiLinear) }
@@ -351,3 +412,9 @@ func BenchmarkScaleSrcNRGBA(b *testing.B)   { benchScale(b, srcNRGBA, 200, 150,
 func BenchmarkScaleSrcRGBA(b *testing.B)    { benchScale(b, srcRGBA, 200, 150, ApproxBiLinear) }
 func BenchmarkScaleSrcUniform(b *testing.B) { benchScale(b, srcUniform, 200, 150, ApproxBiLinear) }
 func BenchmarkScaleSrcYCbCr(b *testing.B)   { benchScale(b, srcYCbCr, 200, 150, ApproxBiLinear) }
+
+func BenchmarkTformSrcGray(b *testing.B)    { benchTform(b, srcGray, 200, 150, ApproxBiLinear) }
+func BenchmarkTformSrcNRGBA(b *testing.B)   { benchTform(b, srcNRGBA, 200, 150, ApproxBiLinear) }
+func BenchmarkTformSrcRGBA(b *testing.B)    { benchTform(b, srcRGBA, 200, 150, ApproxBiLinear) }
+func BenchmarkTformSrcUniform(b *testing.B) { benchTform(b, srcUniform, 200, 150, ApproxBiLinear) }
+func BenchmarkTformSrcYCbCr(b *testing.B)   { benchTform(b, srcYCbCr, 200, 150, ApproxBiLinear) }
diff --git a/testdata/go-turns-two-rotate-ab.png b/testdata/go-turns-two-rotate-ab.png
new file mode 100644
index 0000000000000000000000000000000000000000..b04ab3ccab2f37ce43aa7b4589020447fa8337db
GIT binary patch
literal 8761
zcmV-9BF5c`P)<h;3K|Lk000e1NJLTq003kF003kN1^@s6aN?Cz001HYNkl<Zc%1Bg
z+mEH!mDk$)yPjKBb@kotjO`2{?3uuDiNHWIjKlyE%z!jfG*ZMw5uktsNDv|Ngb*)C
zh-aSQAK;ajyoE$0QQ~_X+sQ*NPU7x%*Zo|+`}b{rzqR)%JHg}F?&_-3R&Q(fIaPHo
z-?x5yt@T@%y=~>u*r@RckPRD;0NJqd2#^gMj{w=Q@d%I&8;=0lu<;0x4I7UD*|6~l
zkPRCf56Iq+Uvw9J4d*$`MPI=Cjb8%87`gh?U)t#FJ&Uo?8e&Y4e|XsY$x9!5ZQ0lr
zclgtn9Hq~H{y!9qm8(zw#f`o8`xzT~Q^!Ah*{y*X1}8&!PK8me0U+TAx%!FU-pG6S
zs3X4U3wW>LO%;F`M8W9s&#9&uCvj}V0u+b!wD*hG9AAFr)^B|*fZ4d4SiJR$*AC34
z9i;7@0zzQguGi<rT1nGbCgXIiC^~m<KK!6;oQ7CP7WyD{M=%CyIsyW407g-@0A)~3
ziQ`D>wwLKRkv!vVI61oWeO01MHy>Uzwow;e-1(jgw1ZX@RVQULK<I=1yeu186pd6(
zCp~^1jS|_OX0nAj!Sit*OB{Jq=Z6C@e_^BU?}Lnu(vShFV60J7RW-eob&tmZz;v=K
zS^!cjC^dfG-5#l#;*z3`8!$&%Bxz!GmC+BK+<nz;-}ux+1<J+%a{T6tPJIPde*?2r
z!&I%xP6`0CDqEc%ltz*ylA`V;TD20vpIA)#6+=f4U{W}JSTqG?aqs)i0HW!&hr_>Y
z1R(oAdd{u2qwNh$6jwcb2?Si#jm%e#2LzM;WE>}wr73`jRD+TlGY0#L8~`$Ew16po
zCz8U9QIl4CuR9_#zWU5Y1G5o;_&9_7*5J8SGemH#C8NSp7pvxr8nQhbNtWj_$}<^{
zVrWAoHVQxvRv)2MAO`@b12Y3qqtpN>(a$aR0SpitU4D8afw^!%?)>n(j=>*=G(hmw
z3L8aWQBDDg@I2C=W0_4YocRb&d?eFp4o9Bi_g0;G67zjlO&t9Rr5)k@)DpNHfFcTv
z4auVF@C}#mZXPHX4hVxknv?-RtN{TG8jMF+;!zYU7%3#XwM9mi%XBgVyBbS^ScE3Z
zgM#<t!|J*QGwT!-YDorvQYQkR<R?KHXZ&XT{^$ov$+Am-eq(^S&<pv+8{bud46yD!
z)&}QbVg*En+rXI@l^iWf=^@}eqnQFIW0))eW1}Q!kk^u<c`5T{72=g%HH;2EON`{4
z00RsefEbzNR)6pFWj^qnY=9ZNFb%o=gXg{HU<6Lv@4SZe#)i~MBn|?#aS9M%dO&Bh
zEDgRRgZ3l}h=m4Gi|Fu|wGm(fU*l;rKR<=gC!+{}Koj8n`B>m}1)fhwRwkMCW(?P{
zfa53A+qr&mjO@YyQK2(E0=boqk*cTJf-?rC&0wA=q}U*?05BOqNU+k=6oSnXJd$C>
zr+@%9sZJU>Svfge405vcQGExhMQ1*O(6?u?Omk=poqS?cOX}824N%6J*A|*N`m_DF
zUUgS)T}=FPA%NWd;R`;T1}LlA$-F{j&b*bw3eKS?Gb=M_0TYQCOqVsFwUMaf&Jg^0
zFk*d-0Llqyao59{4`5=5U2Y#Bs$bR`)jND<iuYZ+l*pwm5CZ@}v$ZvbHV%<y0ou+p
zF`75^X8`8To3H3K-hASG)BG2sA+%Er+5p6H;p7N&0bqLgk|_l14Cf%15RG#sLv)@o
z0TgSE27nP9eVnAAx>(tKPfrzt-+<w3>Q>&`t7QMUR$sCU?Wo%f)H#uzX{2B%y;hSM
zDFFTv0{NAaod?DNnBO@UV9p(oJ8wL%*<7y7Ie<7RogBc-94w&?Sf!mKZS8a-qL64i
z#>!5ydVOrc70<!v<YO%YN2fplGyvIw2-H<0Wl_uCK_UC|R&HG#=>p*QG7lC<@Jn#&
zuxesqD@h1a4J?f&Hx+c@_jn)2@%GPO(fH-s;}14JcJ6=>4c^$i^kAF-lw&*=4Pt|~
zN6gT{_S2NKdT}X)oIp4`040MEhu#N~6j)uAvV^%R@a6%pj{r(h*0KcCn!~#vpA=wI
zsS0)fq=SQxKnfCQBJ(ORcL*wz#HS<s!Kg3t8pz++aIFN!g#4Ju!R=Sw<;NbJe>qn`
zZol>&U2<hhFhqgSSs+Ge2I1l_VQ%UMObgFPwwK%{m?Xpis-_RIA*&02v4pQ!)NtYr
zm{lWh?-vTh34k~}TEYblvO9%_Z6}CbQjKRwsYwRXG*fWEp5Z^JolNOVU|hINL(ttR
zEmnv^-R0=+tF*at>BiZo`OgiIw_g3Wb3y+Xm7fIV2y;<#W%tsyyv+_&6aWmVOhIkS
zsK<=O^pUjCkcNyXHCih<Su{G20m#w3mOFbZxp!E|!Eq_~_7-yM=8jyyl*!dABiY_U
zKsiojdz$)`{e%H#EMtJd$w;1$G=C%qK0+Dr+-Om(7ayyT?M0#a;^=Ouznmd3=LU$X
zzmQE?jTc~P@LZH-9W7c2wv`I8!O&ynK~T|9OYGn~FmsbikN36E9D-6KCMm(n=7>Fx
z=an1+jC)6o>>c6xNu?Td_3~J5Tpr2e*T-`6+Egy>jJ23GgQo1vQuw7@{nALm@SxBK
z@%uHU4|Ydq&x$4KGkq6vV1bORP1F1Vn1kEjgBd#;-*PU1-2UGGCBxSYE&(VIcNDFG
zP!Wf~R5dVR$IBKbr_{cr3T9@gVO|<6+ZJGSSmis3g2G8EgtsBM>$Pggpzr`d9E0l4
zi%yGL*8#}oOSxRx$>oX1W^#E4j16iq11MY2l03!dfIu+U#W*BSpamUPiJB>|$qi^T
zO)0CDJ+slqJN||-QO9`6Di!BqyZ7^#)dF1qjr;2zI~PDaV=)Ru+4)&vvOKGku>hPk
z<ST7eOSvkgDnM0Bn5tn6Ay61JC3N;!nONGrF}~iej#sq)L7`RnkgK~RXhe=+Fjcd}
zXni)0<>r+!G$g}Jqy{_dp#>3zYzKAxtgE=euY7)X+4NeuWD2M2eUQoT@ySduy3~Ub
zM{5f0<nEg<x*JdY<&O-?0|(^pD_>V0PZr}r5EZ_r$)WyDc+T`ms&{MQygT@Vr8EG9
z>5VPK4ih})h)E0sL^yfJ*HS>tRi#G%zKz#i+x0AJUDnOOtYYKi3%+H0nt1KN8fjY!
zv4dt{K@6RmRCVeMSSw$pBv~j4UH~Lj07O66^D`gE#MTFu_d6aqAf%DB$MqVB!PhNV
zKw1OSF^{Fep3w^6fcH=XSW#NQ1g&VAI$uf5JEy7nC~``R#!TB&b2=8QKx@;Hex@@F
z9F4v%;0z!o1dcFGL#wmzJAAIE&EmXS;&m+m#tBUzDD<jq@cFm_aq84rD!#m(>+{4K
z%{Vi2Mi%ziQ-5W1P1eEd-*j9(Jp|6lsG89@Zd=n8UNfL60EjnbjXAB6C+avOU9nx~
zBYQ^U2-qcsSQGH#uB-7gUW>>(jaaNo5tjfw$eUVXBA;rHHOddQ5yUQ}%x<7&v;e1u
z`6|_aP+Mrb8Zf$5#ex`u%cDXcYtV7BFA9BUmZ{mwaeYJTdzM&6aGR_lM?d_IIyvh|
z9a-G4Dj=pk)|%1)kh-^O&eqK$E+qJvSK^I!aHti8v0M~#C5O!h7}SawuPK5lp@$kn
z4JKDISs<<`D(0s|{Y#v5l34_q@IeL|(KW1-^y<(JwW6;<6rcsGCBCcC7-c{x+O^bd
zr}GFj`T(&_+xpztpZUa&GPM!}VKR!9789r^|J<kV8<GbGh(@;7YlZ<B5grrR`nWSH
zC_3^6rluws;@OP`WX{|1zU!6YteEFoh%84iUGNDl%*AAan53-4U^*6~n1D<X_~p4=
zN>Z7?loi#WDRu=IG`m_wa$-C6Az8xfBOskvfD}VBSdD2~0u0l`vM@8vcpNnOILHE2
z`tJT*cDFKh^i9)i8ivlE8u|36p50v9d-&S7oMSfvmyL}jfvxBv<s2)!Uk$pV=tIb=
zhR{*%SL2{|oaT$!i>hkn2xe&;t9%6~%~?Q8_5ohpmk{&<3*57%Y-=i+X`zvs<XGkO
zQ!CQy9)Aaupf(H@D0CTu!&&TUFZ7{jR+-Km(R-(L7&vYOrF7_QNvZ*C2@^&)!pbGJ
zlB|sWt|l<i*)`=q5I{8ZYP>1xVe1`eD{rLGMg-p;syB+2W%@?zxHJV$X{gb<hy*Lf
z-<L%t#j=Bm=%j#&ID`qyq!Qbe(v=IX1dhR&Mi5jCAy;i87Fyw$5m><oxmcPa0~ioE
z0O4#3m>W>xjEoCTUs)Tqxwq=$7{-=}Zz{GTfmkeo?#>Jm|8yi<i1PE)n<^H<KJ&Le
zySb72y%)Zw+A$cKBql;2Z+c>ws*{ybs#%9dI0FETBR>fu_$45;9HP#X%%_2AD&eqC
zz?>w5nTi?+A2uJ~W$`_ESIG|iONQ6qW*VRZuqD19j}?CGiUJdW6wngP7FyB(knx!G
z80IMgQ-gEs)Q1gnaI}=eqm`U2$h2AoYa2jpfhEikU`{Zbw)g(=49A%R(n32T=c%K{
z7!Y_5dtAUq?G;A$=SEAc4V<e9=@jGpeqb?R+yPDr+t7#FgfcScl=d$1o;frn2ir+g
z@2F|VYqkijoIti3=nW-=TY@|Qkfs1OHGLC*E)|SEXh}Q41Sl)U6f6akhA}_|M?XBO
z<OmGw=ooWeK`T<tg%OBxo`iIYk-zu*zkj9)#xn)vN8kFYBOq1h%@Kn@GNuL;^bkTj
zfRGNCR*sfVSC;Qu_g#ls7Hds{Gl51G5O581(F{I2ynus0f`%k;zzIOkjjw++SpF4$
z*MiM<_|k5_k=kPQ!c0~8zSW|aV}4dgu&gX0*rv+7KdBOFVY*nycNR2Ri(UW-4C|Pw
z7^d(Sm3FR>O-(>57;KV`)8t0O*=OO;43O8p{=XHB0a9L8Y$Cu)YkeWAVlj#Wtzrnr
z;%0mg6)e4^q>K>2n8`Y@uA#$9Mnul+YYhrp6{VcOxw8Oe2hWP4fm(<O6l`69X;lDm
zkUJTK9Ou#kkjkX8!1wH(bh0;x231Zow^uKXKo&B9oyqk&krP}8KYj3hM=-CEEC~q9
z0*e+J*?<6K8VqJ>r123QInyPd|KoqSx#{--G}S>U6aqUxR>5{tODbpsERCGhgZ?dM
z@qrBr?7w2Pn=;GFpo}I1<^~aW&8#j=T%F@{V0Ox4Df=y*0a1pdg_DipH*8aTP?}Oq
zK0=^|SR=|*OI9s45<ndyz&mLr$M=>{m8$Y{0qvk!Bys4(Ye{PpDX}o-SWt)1NP>`S
z&I({eK9NZO_oXlX;R7jBoe3Z>ef57irXU;`%7Qg=IB&#JswHIj4Qy@6zp?NXC)zm@
zXD}PM`4Lua+kl#OR-<Rbx?dXrQBW`uZ8upD^U}$|63hzGbviVXfT2Y#)_2pwA7Pa*
zV6x~}K$H6jrV{|e{3|F23jn_GTO7*JQGc*VwA+M!tcP~93Npk1)4*)udYA;=-5F^P
zD~fz+mtZ9QgWG$0CV)^o4D3G1$fBgl5(Q*Y`J7Y-b5KJ|ik?AT3XMo*1h$fmN1Ebp
z2wZElJWEaKBNmSEz5og=3_t-9NU@Ru%;LD0+h9N`gF0wSJTUFzZ>4i20njASS*!tQ
zi8nz>rtONdgq9X>7y7sb8a;+7iZxHD>0X$x5@9O*542=D%hk_WpJ>n^Z~#I0^1u7s
z=9X}u|L_0Jxy4eN(G);fwbk??EJC&L5_KPIL{2B0XipV8*AwOO3Dek!$VPCcEmm#c
z4q6`pBti?8f{u<;*nIkv6yyN2N_=XC(evScF9rNV22-QiR`IONfWlH^7>w4=GpICr
z=bEVzG)Hqp^RrQ+*<$*JNJ$o$8Ph$QTk}9LP^E@2jboz+yOTV+pEtkor7t?uXq5>7
zSUDLonIkBh*AZ4&qJXq6!HP@78Swp4;%Ou3;|Pz&tAsje3M`TkmpCR7d+jP=n5O+$
z08oX#E^Nzs+kw?ai50Z0>Z~Rz2SerrYRU;;5+{jLTGrIrKE!5jt$i4pC8v3{!LUWb
zcTQ5vW`X8%rxxq+QK)(R=db+9=60Q!pWPORP`E+xRBd5pjo^HfOH$Y^X^jU&g@Ugc
zgVyGtjNA`M#=_h024(<54O_R|cFLBt!N4%M9JHNO|HBH6CL<N7v6yP2A|_Gr`b<vd
zh0Jjsbi%^`;X@;h#7e7e=%6=tT}T=<_|r$N^BbR2bM#>Bu#XiC0>KI+)9lm^YR!Wl
zA6Y|Q{H@>S7AG-PEvtb}FjjlbNsVZuTpF7Lj70lzmB+tp4xLYb^U=XR4eQC>I&|*N
zSM%2WSJteX9hs!eG+tnnw(c$A<{sxw6)RYzum&v7Xfw9;4cT#|8$~1CPZD-rMV@{~
z?4IZEZKTh)L35d6XU=Ooq}n<D>wo%Z4;+z?tRaq^HtE4gV3OPkZ2AOr6|{3U5rf|&
z4A25nZ!?)@MycaSr#<a5q;~Xeqn(Ema#_RRf#F)V#-~%9@ii3oFdFB}!Fis^Y-dZ;
z7p%!M-EILqQ!HIzL1aZL@W=#~Bw;g<RbU3^lbQUT=0UADRgB76<|QE1W@^!Y`H%l}
zb2}HGyZ*Rypq|5ks23)t1pwUwbksm&#!dytI75;gP8r&f`i5tgPPxU56q&xlZ_%OI
z=r!H$6O4W6i&{6;y^})m^(()rvttDipT%tnlc}$H8)`u9w;~L_rsWw@07OErjDOJ2
zH43|k=Fyx+;9e_yX9}PD(<IJL;J<EsBtV>&<PbzT22xm&5#oplUUkxYFcMJAd@Ob{
zkx4ut3U@NejCOZO0AdXWeWlaoN$NYOwybN@Yfw}t@Jj{L_|n|E6(`3ECPC0c1;Tr9
zQR*Wk9K~S+V{)f}_}@pNS9sl^KmfEw2$OXWE}#y`{)aDr@!T4&KQtiEU3pA9#<Vz7
zwJhtvW>*kuq2*Lf7LvhCX1Xt90wGQkfRZvaa9X>Nb>A$rz|lIr>-~UMy@%nnTadpR
z#IIvMvIfB%B%7oAKGA?q9j^9KfiknvO+*W0K4t9pE3VH;(NiK=ag6<p!rz#N-`z^J
zk9c=G)8{mO_5lWg?=?OYAkW=+!WWN?3*CIUK}L0MC1$eYyEcL&&W1`5RzWjBb7P~K
zSoIs?)9n(p!{^(rZ!lnQt@x80RKQ`sTQDm6(jnM<vMFm*@n@n0v&YGlzG2T-7aFCV
znk$+;e@cX$o&>@6AWiS%Z&v}rm2LKa<@%^`2^P%1{#T#h;Oz)jH-GB}f2E|M(nn!x
z*Z8{_(YRSi1{0KlMdk1#u{!1ehV#I<wRj@08EbKw-h*w!zH`;P-~#+B{J&&SjBQAp
zHTBUC3Wze;%4icEUz;S}+{MAEYfQriH@zh|a<IN@JDCS%C(|-7wd3<&_~UcmwDW-h
zdGXdK{rOei4WEJmek0L%SCVNC*AeT!zKlt~wvYs#H)f4m8`?Z;(%{c`$_8OOHFa8z
zioAjEgR$4ocP5+yWK=thd#CokF8OjbGiD$t?4D=45S_d>IpBJ@zmPdAk^r*3$v9CB
zqiLnCPVz{uQ9HI%)ezq@dNF6EePBSsHay&9po*lbEScuHSmt~%Wh)oA#Y}50ZC^<Q
z7G@Bs_l@%(pmyB*8S8w2=7MmIz9ax+y<hRoDDQNjOL)D6=W6B{uRyG^vLxR+q+6t|
z-w~B*Q=91CLNYzJma>&p&qj%Tn-RMYsU??j?NdB5E%S$;`=ia>()Q}Dr*ti^8xoB7
zm11Ek;@L>S=-Ipsbx%Qi>xr*r)OPU6LT<jX{@gLUrk#A?XR%|%o%R5ULSO<gr!4Co
ztuO%;<2+O2>qriI?rDb^jo06S_l)!;2dAJw31oL-<@!z{H?B^!@tVLfFUFJ%<Bgq3
ztoIn%i#z%70|8=uc0mPb!(yk;NVjGZAvlY-M`8g+x1C5eHB#Zm>oGI6kyLH1B@9wq
z7bvc#Fit_y-RbWt=dRuEpt>=Hzsaeo@@8xOK0QuD4}oynnaGoxXrdU5(oT}u5tK20
zcL~?DI~KWudE;8HDV3JBW{fl(KtP)N?BD;#o4X75wa1?F*`Bpl7!O9oG{ZIknN1~{
zPE}Lv7&dxl0ZIfwVkzjF5sBAyv>hX{M9ZY>(ybp(eFH`B50YE^mUo%6aH8BA6cbF9
z2_DVpeslm4`_6O<-zAJ$L>u)oXwEJ);_^ts1VE-$jb$^og!@`|E?wk!UqD`e?31TQ
zP_FH~_r0KmTL3`8nE0e$nrN<?yIJ$CSZV}U{lQ8#q^X#*VnHwyGn&A;j%1y$)yRDC
zfOyIfl8FHrrM{sL%B{~i2!N;`QQF(9*`sE3@Dp_<nnB?2X)@WMJxe?gS4M+e$(a!V
zAgrSdqG=RrEx~x^cmLo*A^89%%XpKa!;3_db9pAwDEAZ3HR4soG%S;5bu0#wacsX?
z0u-1nNk%YVjV2XG)0sCxVF#_#&7MBE({H6oBPa~|NV9tv0z@A-thHqbIWWxtIPct<
zx@VHI!I!ccxdrvn3;+Py0_}@kt6eDChXqi~;M-RJ>}USb1?$KA0OGa7kKiM{FLBx|
z<<k+)X}Z`uQcq_M+xy#6+T)e9^Gcd!CFP>gDx8Ze%?@*KMqQQ?rvNocG`B^{N}(~p
zQw=dbi#v3K3YVFIkSPzjs8AbOn=kql%(n5lELI}9A(X8`9k)T?XC5YSvWoOOXVX|U
zljcs7ofmDZ!uV4GdFJ+y^$K6V`IJgx?QEkE?HaB)E_?RZ&KJ-a&`k&WJJd42dq-*j
zQlAu3F8~U6O+~qM@d&oPlqf6ph$Wv!PsDO7EZa<(k{?2nQM^7f)7ph?WyY9OK)kB1
zcg-5Dj5HyQ(uW;SZRa=ol&oSF?CBEOo8y+_P6pVj$(RhQR}1i)zxUwEvL9yrT9ft6
zouBB1zH$ADH8+K$hW#zM1Ez$45d}n=YsAXqifsV`z&Kn=b+VFj4sDTCTANE5vm84{
z06CT%fTURtElG97TI+eSnl1W-byk4i*n3wtgaZeB8rNt&Lu^PMTBD~qXh1mlKy<`d
zgYU!V9WQE`fq6|Qy^O|i^zbEYRv6gQ6JJaa_+H}!2Zm3-{iX+o`7F(0!44nK)hYXw
zd9PKzYVSY__>O9^km`6TZPhCvjm@NimbB?a`pHyUIQJ~ewfP_+AfckC+GnFvN9#gw
zsa+=<ozI4`80DD$U~7)$wgzZOJ*a6T7&K=*3eXIm9uinq1&hXaL2|DOQ!51TnSXG;
z`KtFDA3EFn^xZej8#kZyX*<%{6W@V8Kx^8fL@>x+tXBFG&U`-CV=ii8T|Jc=fVlCD
z-ET6Cr&8r(up*c;=DHGJ>!y}aW4NQA9q3wubs=Iu<$Oj-fO;@zHG9jW{<w>mZKZnX
zOdOq-X*)*q_`F0v(`z1m6b~v$hRx9qj*nCidHT*zy??p&Nv$9a?5leFP+%7zT<c?s
zrj`s^G-^9pRV!^w>tL>;tx`J;$-vBbA_D+vVZH`7B9JTq1lMt<(10~4>kT~t5td(-
zB}F=Uk`VWj)N4Sv^4qX~V<4bbw?wc5pJ|dT<`!!<Y+L!vpL}U!w1eZmnjD^a%RBy8
zKk=JRwPjiOE)qn%8CLL^TMs(Fb-4m)JQ<)T|L<GxhX_nQ(macU>FU`$$D^0}1Rn#f
zi&HVso`_8i>$Z3i?G}-+W47%Oi0Wx^#&5Caehr@O%TjLVw?1r*W6g7An%go;i@RFn
z@BWK_vr%AhoJphHr+)q;^YUX)y7eB97yyjd0cD6GOccusE4={4Ef>e4jaL9fY5)Qa
z0qqh;O~LG$ummV3W1&1#t+BzMaKE8X$NMIQ3O|=R35wOKV^?$QTIbUuuyfNak})@o
z!Z(c~r+*vHaFfkhjx+1+{hiz1j6HYxM%Y^jK`{NTkv&(Yl_j?uvloiM)TNXJq!|nQ
zho~V=Gm-!XAUQ<;?DVlP_OUjsMe6LCZ)$Z{QzOif^F)Ud!sf|7G+iAV?bvYF`dDN)
z2E8r8uF_H_Y}ideC}van+*khn#%Tw~nePwz^ub-d(eLbB)i{m~oR_(sjj0VhAK|Yn
zJm%1*V&E>N6@*s<1OS4#1fCw6f;h!&Y3E`(_Iv$!#2<@EdJG`T!SMV0q}Hn%U=R@J
zwIleC0YB4WPoh>P(9j78zy#XXw}qIZJnP#&v+=+-eg3uMJq705vt6|fBWeji(FwE0
z{f0SQ0lbbBZ3B&|r6tRP>53+x#IXB$o5`q8^hBV9=Pc?*bKj4(IO!yhds#WWuc*Q1
z;5SIc{pNRW^OYUrWM(><6s?TeXH3momGV#i=l|G9?cjLe2M0X6I5gh`DC?7eCYUKf
z!DEIMoDnEKcRT~7k!4wGx5o%tGlpZI3=ljb4|X&r(3U7q)sHmG7Via_9x6%uQo5>=
zVnq{3ve3&6pSK0RJYz(TNGOAo>5ELh@Qtr+EEpW;a@NV;1Sov``t*`pe<=mVJ@+hI
z<d~yi(xzFemE7EQG6fJ52sE?4$ixzCoUxmI1pqw|Sy#veKij61r~xQNC9XnL4sDx4
zp~iUc$il(1Id16B<96?Boa+H*zq>qAV7?ZBp=QvxjCn8(pJOr}VHBR|&-k15bq>Z4
zDPtT2oz-F~2?EcwSxQ>XWeS22b5;PQtYik0nDT#VAq3KCcDMt}`<HirvC-PWaUq9b
z{$2p){{R@qd}B^-L*`Cc(p#oVK$lbag&Y)iVw|?2av$R2XfDZOB`HA3mP?r;>Q8C_
z15=k^W?*0J<<-r)pa40T7C+NQo39-l7kbqA-v%(>fPe7`Tq85_Ek?T0zicIfPe>OK
z_(?ALbgT^wOkAc*m@LeES*m$TO3sxetCemmAmwJ8Xv*3MH{sytf4iJ-ymoM0_}TBD
z4gTd@J6H77PsS|?cZ;%5*!s?Oo5bQU%X^0!?dJd`2T(JB5mnF#Y6adyvS94skCC;$
zcbEIlBe$7B;Cqb=2V_0|#__%a^W7^qoo<zeme|9&nBJFiYdCA>o>(vG0YtVcpdDaT
z0K$gpcs64c&Q^*1i3}0qoqxVqJ`~!)v9TBW!j0c>CLEf~<HzkZXKRS|YPD^v=T~uG
zdW*QCgvRiso+0c(Gro&j+wuPB_;_RQ@x8`IU+hb_o^+@G`&E`(L;a&SzG<6>Ci6(Z
z7TTjV54IZlO@XrRBL0J;qm91DUmXv{#l8$+{0U*9H_Q6YV3umx#6Tjj-fBQ|#ls6a
z{Qd779zGQJ_wMmfUGU4-A9tsqh}RNrjtD=y(D>ktKfLU-hX)VU{c${W7yR<ITh6@;
zly&dc5PE>uWs_$Q_aC}@<9H}9n&YJ_H%{|or<;2MD6F6V^~1e~^4=cGi+=Zb;qrAS
zZhiEz2j$ZbDHtEC@%-*J$9d7$alhl@`tHwo{MEyIkF?^^IQ#Lj(U6bJcm&9XjYoiN
j*mwlUhK>Ih00960W3{zdYQ>Vg00000NkvXXu0mjf#f^h<

literal 0
HcmV?d00001

diff --git a/testdata/go-turns-two-rotate-nn.png b/testdata/go-turns-two-rotate-nn.png
new file mode 100644
index 0000000000000000000000000000000000000000..da93978fb8192bdef63e253f78c223388bbd9d5a
GIT binary patch
literal 4993
zcmV-{6MpQ8P)<h;3K|Lk000e1NJLTq003kF003kN1^@s6aN?Cz000w6Nkl<Zc%1B<
zd2k%nnZRGqDUD`E*5T1@*|M-@@PUovlx++?z`&9<BxJE+Lm=!zg{`D2%PxTmZ0(WC
z{*%ogY;B5cQk9}OVo0q`2(jZ31M-b5+sL|X%ZCojx}}jdGn#v5=X+1D`@PrQv9U)p
zEth^(U-fjidZwp;{k^N-$b@Dw6thAM#;g#7F)PGi%nC6WvqB8UtPq1SE5u;T3NaY7
zLJWrCg$!RhmCW`zWD+UcC-9dLWz4NvW$5!w2ZqlOmf_AXZux96>>W4W)0B+H63PaN
zKh?PlRvPx)Zw5ntDUd;5%Tfd(NsA3HJymxY@;lrphA*<L78CLY;@mbW0~os8mW0Uk
zc}*9F6T)T(i`O5+p@qegL}UG}w}ceK39(sCI4DC`qul=Y_Mg)|7NpE{F@`FCC;OYZ
zi}`u`2&(kBxWD<+=NIMI+aFwpvx7xKN<*NaY`LMn-}PX)3=!!a^2x7R8NmDdTF}v^
ze%)ieg5a5eVeX)ji;ei=IUIJp4%s2nc4>%PUr<@Zt*^hcM0JeSY{pBcXJ)j<43QyJ
z1~54p*8H2pr~8rD7smmJGJwHxGq+A<0REx2<m5=Z;WNe%LPSOac2rwkA+LiB7LZa|
zIiK5BJHH$+K?(#bg9W4*SxnZ2TyAgBy*MC+6+_u}{~Ej$rLY`17sWv_nfMaRs8X_C
zuChXM?ItvtbaVS#nR??yYY$#hA(YbcGBiHaEw6WvL{*=&VTp}fH=C7v(OyuO_3J>%
z>WPTiF&HVtOOI1&y0!`xd3jkWw+=E`K#VA4GD^PY1Iw*=$;zNaVhQA$dwbT!WbFh7
z&Yes~<9yGX3ZV=Q4RGsstyn0(?jH)^r6}ZAzdMVUYZh7XlFDHDtCe<cJ;DxRrpTt4
ztkxk`vjqnr$^e$HSfP7tnI{jeSX!&Q|9FS?FpGCQfZqLV42K^5-KVowzhvc3fb0;r
zy?H2%?pr-qUN0@<`xI8>SG|wNqh!b6SU}DW%Im)&>%V$+nY;~zuv8RT(7M~5^16SL
zFYd~!S7h`WWxbUJfDB;&iDA4HWdNSSIXL8X4_{QBzi?Ui{h#|u3l1$Az|nTKY)XH3
zQxf7DUm<N|JVW~0l}?~3Rk>0bEM}7%J@N27dA*|0!tIZTe8}wdaN7<aSF)qNH6UjP
zW$lV$`SsWDCM`KCgtD?WAH`!KW{cWGS@oUMdd;$$A&GcGwGCtdJNNY9&?D|wT3@UT
zVB^v}6!)ozUML?No>0k<Ki`srLkqjjjK=y~rZrnKDum4r;O7To$d_sb`7~NWl&<D3
zZXL2iew~uSilJ=#{-4PHDJg7rP@Dx-I4F=EEM|*c_nbh=bcxBBklyAa$wVxoyC0;y
z@j*AYZEY3Fpd<AcaA@(HeP>Z;&kzoP%AkOhZ#}Sp+xE&cO6v`z-1C*in(xhK3uFh2
zBuU6S(3*tM{Zco8cHT}SRLUDiwC{!g+0}x+_t-MceT|2k@p9Gb1)BS*6qaB7_#wQc
zQYy>bn#Z~dN|Db?`#t<M?i`rJ%er|k&0|4ILFv-8)GZmQ)FC^>ZM#}YkwBjazVnAJ
zZog^rDth}z>$q(;%Otmcc7$&bys=w(kFyuOcu9p&9^bkMFORm6i@&PYp*mlmmmkc$
z(98E(mQ<^yhQs7WE3Lk{xnmhAvMfoyZBJzYo_woH2IZ;$Jc<JlWdw8WsJcq~ZU^cX
z7N8$Ivy$6JWdQ4zm1-UroE%3MdoB(@lmQ%X=SL8z3}B%9=4M(>yO4p?_4ouq3Ir>I
zGQXw-hkR10>|kl`9z!QChBV(Fq(HFQ0X(>2K3<-^IDkX**n$um9aSzk$Y22}P3QR5
znlC_Z43GlB%18$pZ}!eF<THBru?c0LAV!2rf$R|Z>G!K~=&|nheBFH@g%v}2Z|`Zm
z6s53ChGIA<hAy(kJC1Yd+t5!cMUWzjq4b<j^|+z9m?|<}`rA5gKS-&K@nyU@p=1;m
z^BwimM4-2K5C?Fjw~O0%@!S`vrRgLNEkMdEKiki3yRAYgGDR^r%CuGLtc*~Qq&XEr
zd8?7!4<kdQ0>&}~Uiv0~*FRFn*Cm}Mb)}5_WGLj53*_^j{cFBi>z2sp`I1Tp5QMO_
z_3;hR_x{ktt-t!+b(|2OX7QSvzMdH=vT|}$7tPHVMN|e*>?yz@;)+k8U%cOf0}y2Z
z!GnLqA>&uZrMA@ikkj9(dtNRnghP~(Us#RSEv>+zhc(p?>FQHDsp*gbOp6e)sG@vz
zu>*$|7g{eMcfR&M0HnOPhrhV(hY!fFIVl3a-MJ43K(m8od6j;>1yXkZQM*$Luu{^&
zv~&!kesU*Uo#X<Ctb9?mzcqx)Q&mJ?P%VH;NsjR?qo-GsdxuKd@vypcxZX#~|05D9
zq8gg`${3`$;z6|isg3fsM1WN0R7!V1dmPJVCPVa*xF7x73+uUk>B<fu?HMwBs$QA7
zo)Xj6F{=Bw<+(&0WCzQHcO3bo92|PIA34FT-?vg3{av3theHHp0PnwDkC#*iOJZ~g
z2j!oCdVt&i-N#nxo>Mx=NF8#qtxJCWlU@8p@J|F#>AW)CV?hed4$AU5Q<9}T5*qyS
zXaG3_%4ji_K>;bB{GONE*`qYUL5e7YGB6ax0c_s1C9Mskv<hiE@KJKer<AE3eab}u
zDXbU~kOBeOA@b&_Fb+MA9BfDb{(8OUK9ItSp%hw3hYqBOVkmhw5|n2jUBvCDQqoD<
zgn$hCnw7$`v&Dvz(GU)xV^Cf9Km5>AZrv#<?*&r!p1P!a4pa(c2g}R9JcEAt-v@E%
z@lUV)j(@Lr&d9GtF_b5_e=Y6ZhP1sQD+KtUfpn5K-A?KOkRfv0vL$Hm0p$YUR?ZJ`
zAJ}(Hb&LmHkzc>EIUgN4(8%rote<?|liQRLV-UiU=W?RgpHsx|J>JXh|K9VPIU&HB
z6-u$Xpu#mRGJv$aY2)wTO*-=VP7)PD8Hn@yK?a3N5rw$TWKi;>J)}rsWl%PlJGgZ!
zgmQX-Z%Bd+3Y8)XX&ff!QBzGC3akvut1mubY$5bDh!JIUhn3kRb_=NsSRsHdauv08
zEA<DJ!P2njLvH<ihjJmGT{5YX0rVtE=gvDwhD=x?fH^l8-M5Gz8>2!fJGSxDk)}%q
zaMN`ND+74HK|&E^)&2W$%&sW|xTAUw4lOKEA6mb4v%JpA0DgVOflgiFg-{tlK}m^9
z24#Oo7zh47f{d;AEJS1D`m<8f%}sR-Ppw>;jO*;9Qu@RE1zvksCEg485*jsSiYf=G
z^sk1=#Q_-<D&>5~Tih|;0C_JerMO%@z*tyDik`84GB^p@!SaRIjx4Dw5{sz^38|Ey
z{l|00o^@jW*Z;x+d~{Z0h}*7xR{}0m1VxAS=Tqg|wRML=F|_0FpTwa@GRc26CsW6P
z7%E2#i-p{$kR6#u`$rVl0b~ceZmK_uyyFQJsIHN>^`lDi0)EDvinVuoRPPfE1W+^_
z#$n5`qlsaD$RytyMY(n}3Puz;F=j6A(pAflWX(g}oyw4WQ5mhXBpi~I#P7{wu}R30
zmy0792%)06e51o;GNHSc6sJB{`K;;A4*a?4dZS~N_u^&AA@a;Xmcd|2#Bfl?rYc&B
zw_n|gq7=GpI6}fRhqe~dLc3?P{*XpNp%Yon%7q0PEJYe4&EMQgvZK&x!I23WzzqpG
zKF5=cbzVUb0F^-jDehAKOKxE1&tfS-^{6WlkHv8Sc8S0EQZk7uigR&j0a6yN;ES2U
zB4q#*q_AQrwp0S!QdtNAP#M|arZa>}VY&1-|A?0$MHItQ=vFHA$Ryvp2Pwgr{(J$D
zLbF4p$m7PLiPeH6lm5*PQYwr2*4nE-{7%NE+<xg+7Gwvz-jh>|N*X`LuLr-fRr8n!
zR=@yYQvR(*EUa97DrG=;_9RK<PNq^E8nxY~d&uM6&kmBQqywY;<Op|ZK7t&xMAdy7
zD#a@N&ix-ABkxJ2^p0pxu9z;qRE9KGEKEvy+LOfq5EEI<k43F5QD$X0tlD>cEEYk<
z#V#BGmBA7<@kXL?rIH6Rb&KZeo>$B|)i1L0_us;57n)b$#-T?hvje!HA`7HIKz3+R
z+5BhSV{>C(<T5LDLrj=nN~N$^twaV$!4cPP4sTqm{D3Da6ianAPE=Q=|9l0c(CnZ*
z{p|Gh$o07HA`7xZ+?GE<nwnSLq)mVMTQAD{K#F4-Jh{`JkVd6=g+5p`6h!W+q1B0i
z&YUOOQZtd$?;lVef>9}dZn&s<D3ua5sYerP%cdsB^{0PW{d&3X49P)!F*5E|9yEa%
zRtC%HzyJ;)8RVaJ1cQ?(C!t?PgBXWJ$c8Xe*(9{5-JELgMMW+<4q$1OdLGxy+xcf3
zKlr-RvH~&L$PVDTKD4D`0hb*E8V`}mUwQ#AKMs(NA;}H}S7Pcpdx7#>!eEfE$$b$%
z31eaYPK`=7UPy|ZBoW#<Ova)>hRA$p0EI2elf@UGxP#jtj;M>k^$$Io6(KiNhZH*b
zN|4F`w&q>Jp+|HovqTwy+Zw^4hZNSI;jv(@9ft_W0PHcc6lDNq`4$|MY{>wwDWtu;
zhP!PeQBrh)6j2Pta{nf@abXAtK(m7aQrs5)qot(sh!&*SBHF)_l~5c9K(m7~nrilL
zTw07n3y@Oa;~(c1<|@}-6q619;+??Ch)d)p6+$U^>;b$CddNip86qdo^T#-x{PZ7R
zR9%UN0%UY`GD#}hScHG3ZV8b)-fUMom!m>Iq^pbn0nlrY@(m2HYaS;As985XPXe!(
z>%JY#NeM~u7>Y0ssSF@@h|GZ!yI}iwDf3|qmTy+c0K!uj(q2SHi>VC2AqCOk<sp>o
z(Qd89r)Hmz`+2GR>RaX4b@R3RPIKY`lt_J)$w`KFvMmG9EwaR-isDFRibK-w;LTY%
z56x|F(ETuP^QWkO_iwpvBLm9lZdka(S%O0YH~^X*ECu>q40N{3hpeFq-TeuF<{m&@
zi)-EyRE#A28w@n{9~kV9+$V1WAuK!3@Rjg+VSHrcHy-DN0K?s6W@7PT(w~6rP@VLo
zA2~~v=NVFr&yJn@T2;p^D^+fymv<~eNeBM`t9I*-tY-&s%}V|DzDp8>-v021q(~BF
z?E3jz$oA<6IP^HQ^KEXO$^b^rUBDp%G64So>6D5xfT}7F4hky+=<?}5V$L>S`mQ)r
zGA~K9gA$!;F;%|yjJyqmu$*k}Mmam4m)Ff<zLH+D>3&WKV2vv^Qpi;#%V2R>Ql~RV
z8Bspm#%I*<*aW)G#Y@}%u+sFcE^!+sJAfNHkuyU55*3=T4F&4N10evFA_}oZ$;Ear
zzEj@D3IW1qzMtU`MruKfmJonSVP#OZJoWs{7D53rq6{G9_o35O%JX$9gOXP;Cslzw
z!fm5M0FWZeNCcEhIqD`OyP}YjUE?TXReI!92Fq(ZYPfYH$pF%FGR^a9m0wHw`HKA4
z8pAAD2qh(kl>x*|q~ljovs~WB3IQUaAj&IO{u2Yp;P$7)h(cc2pmb!5DvAv!14zp=
z_2oSmbidfsE+-BOCnJUV;D@^JJ9w!Xhn5T=(5F9RrKhDC1+R3;>#Pi5^osruIkQRH
z?~vn#B}(sJV}8^BS@;G4```A=Yq!YR!Lno_83h>ePyJ{n#8=X%J35hl{&H^Hg|qzd
z`-?mH{g1!$Bg16}koMn)tQi}WAhh;H5(%^|q6{92|N1a~>=TsucrUk&$^eR;ew7SL
zxT_flAj){|>M7*=-#4%OJ7kk-ze|p}i*Uw(43RxK{LIC8fbTfo@xr%N$1GezriUIn
zL^_m^9W22%@&nG1g?G!_KnTm<b#@psJAi5b1=ZS#QOz$81=E9}2o60wSNPT$l>uzn
zP>Ms3;){pzawf?DGJ0|coCQk#BJ8uL%CFr_akgn2Ki6xhqM8!|G`_3+$fhAwB4-E7
zy^#@ASgwpB7+pw4|Fy=wey;?fB?I_|O?&b+@U^<RXh~8L#>xP`nlq^Wya;JYjrH~5
zFuDvND^FuKr=-C0Qzm05AjM|SMeOBWvt_mXno4OD53E0qHqw`_N?FECBdh;&5yaRv
z{>wB_@>vrOEiS#Y2i={}pVhe0*Ndui`SJLd!Xd+F2axrD6TANEs02aRZFYX%ZspN-
z)c51JG>-)-SWGHASmw<kollS=iZMEdyK5RfmE5gwWggL%|I4>=cCcLiyiI=H5+TiI
zDunX*__*P|ryhnqQ4qt*z)ZEmLW`A$Q>+jm&*eaU%)_*qBm*$?$u?C9kC7vMpJmXi
zo|ytd04ilx$TWb(MARL~MCzgsPbdSAR0iel(UBSP{a9wy2{)A}^-Dr{#y^=<4`j{?
z$p9cllmR5B#>zrdKV_JSGJqL*(!i-=W#&RcW8iZ2*wBppZZq<vSwM^^WTwivr8rUO
zNkU}%JZ_e(_tewHY*x%NE!;90a;q>a#9+({F&MK#495Ql00960eIlK+OHVA#00000
LNkvXXu0mjfRReuA

literal 0
HcmV?d00001