Adding upstream version 1.20.14.upstream/1.20.14upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

10 files changed, 3781 insertions, 0 deletions

diff --git a/src/image/jpeg/dct_test.go b/src/image/jpeg/dct_test.go
new file mode 100644
index 0000000..ed5b73d
--- /dev/null
+++ b/src/image/jpeg/dct_test.go
@@ -0,0 +1,299 @@
+// Copyright 2012 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 jpeg
+
+import (
+	"fmt"
+	"math"
+	"math/rand"
+	"strings"
+	"testing"
+)
+
+func benchmarkDCT(b *testing.B, f func(*block)) {
+	b.StopTimer()
+	blocks := make([]block, 0, b.N*len(testBlocks))
+	for i := 0; i < b.N; i++ {
+		blocks = append(blocks, testBlocks[:]...)
+	}
+	b.StartTimer()
+	for i := range blocks {
+		f(&blocks[i])
+	}
+}
+
+func BenchmarkFDCT(b *testing.B) {
+	benchmarkDCT(b, fdct)
+}
+
+func BenchmarkIDCT(b *testing.B) {
+	benchmarkDCT(b, idct)
+}
+
+func TestDCT(t *testing.T) {
+	blocks := make([]block, len(testBlocks))
+	copy(blocks, testBlocks[:])
+
+	// Append some randomly generated blocks of varying sparseness.
+	r := rand.New(rand.NewSource(123))
+	for i := 0; i < 100; i++ {
+		b := block{}
+		n := r.Int() % 64
+		for j := 0; j < n; j++ {
+			b[r.Int()%len(b)] = r.Int31() % 256
+		}
+		blocks = append(blocks, b)
+	}
+
+	// Check that the FDCT and IDCT functions are inverses, after a scale and
+	// level shift. Scaling reduces the rounding errors in the conversion from
+	// floats to ints.
+	for i, b := range blocks {
+		got, want := b, b
+		for j := range got {
+			got[j] = (got[j] - 128) * 8
+		}
+		slowFDCT(&got)
+		slowIDCT(&got)
+		for j := range got {
+			got[j] = got[j]/8 + 128
+		}
+		if differ(&got, &want) {
+			t.Errorf("i=%d: IDCT(FDCT)\nsrc\n%s\ngot\n%s\nwant\n%s\n", i, &b, &got, &want)
+		}
+	}
+
+	// Check that the optimized and slow FDCT implementations agree.
+	// The fdct function already does a scale and level shift.
+	for i, b := range blocks {
+		got, want := b, b
+		fdct(&got)
+		for j := range want {
+			want[j] = (want[j] - 128) * 8
+		}
+		slowFDCT(&want)
+		if differ(&got, &want) {
+			t.Errorf("i=%d: FDCT\nsrc\n%s\ngot\n%s\nwant\n%s\n", i, &b, &got, &want)
+		}
+	}
+
+	// Check that the optimized and slow IDCT implementations agree.
+	for i, b := range blocks {
+		got, want := b, b
+		idct(&got)
+		slowIDCT(&want)
+		if differ(&got, &want) {
+			t.Errorf("i=%d: IDCT\nsrc\n%s\ngot\n%s\nwant\n%s\n", i, &b, &got, &want)
+		}
+	}
+}
+
+// differ reports whether any pair-wise elements in b0 and b1 differ by 2 or
+// more. That tolerance is because there isn't a single definitive decoding of
+// a given JPEG image, even before the YCbCr to RGB conversion; implementations
+// can have different IDCT rounding errors.
+func differ(b0, b1 *block) bool {
+	for i := range b0 {
+		delta := b0[i] - b1[i]
+		if delta < -2 || +2 < delta {
+			return true
+		}
+	}
+	return false
+}
+
+// alpha returns 1 if i is 0 and returns √2 otherwise.
+func alpha(i int) float64 {
+	if i == 0 {
+		return 1
+	}
+	return math.Sqrt2
+}
+
+var cosines [32]float64 // cosines[k] = cos(π/2 * k/8)
+
+func init() {
+	for k := range cosines {
+		cosines[k] = math.Cos(math.Pi * float64(k) / 16)
+	}
+}
+
+// slowFDCT performs the 8*8 2-dimensional forward discrete cosine transform:
+//
+//	dst[u,v] = (1/8) * Σ_x Σ_y alpha(u) * alpha(v) * src[x,y] *
+//		cos((π/2) * (2*x + 1) * u / 8) *
+//		cos((π/2) * (2*y + 1) * v / 8)
+//
+// x and y are in pixel space, and u and v are in transform space.
+//
+// b acts as both dst and src.
+func slowFDCT(b *block) {
+	var dst [blockSize]float64
+	for v := 0; v < 8; v++ {
+		for u := 0; u < 8; u++ {
+			sum := 0.0
+			for y := 0; y < 8; y++ {
+				for x := 0; x < 8; x++ {
+					sum += alpha(u) * alpha(v) * float64(b[8*y+x]) *
+						cosines[((2*x+1)*u)%32] *
+						cosines[((2*y+1)*v)%32]
+				}
+			}
+			dst[8*v+u] = sum / 8
+		}
+	}
+	// Convert from float64 to int32.
+	for i := range dst {
+		b[i] = int32(dst[i] + 0.5)
+	}
+}
+
+// slowIDCT performs the 8*8 2-dimensional inverse discrete cosine transform:
+//
+//	dst[x,y] = (1/8) * Σ_u Σ_v alpha(u) * alpha(v) * src[u,v] *
+//		cos((π/2) * (2*x + 1) * u / 8) *
+//		cos((π/2) * (2*y + 1) * v / 8)
+//
+// x and y are in pixel space, and u and v are in transform space.
+//
+// b acts as both dst and src.
+func slowIDCT(b *block) {
+	var dst [blockSize]float64
+	for y := 0; y < 8; y++ {
+		for x := 0; x < 8; x++ {
+			sum := 0.0
+			for v := 0; v < 8; v++ {
+				for u := 0; u < 8; u++ {
+					sum += alpha(u) * alpha(v) * float64(b[8*v+u]) *
+						cosines[((2*x+1)*u)%32] *
+						cosines[((2*y+1)*v)%32]
+				}
+			}
+			dst[8*y+x] = sum / 8
+		}
+	}
+	// Convert from float64 to int32.
+	for i := range dst {
+		b[i] = int32(dst[i] + 0.5)
+	}
+}
+
+func (b *block) String() string {
+	s := &strings.Builder{}
+	fmt.Fprintf(s, "{\n")
+	for y := 0; y < 8; y++ {
+		fmt.Fprintf(s, "\t")
+		for x := 0; x < 8; x++ {
+			fmt.Fprintf(s, "0x%04x, ", uint16(b[8*y+x]))
+		}
+		fmt.Fprintln(s)
+	}
+	fmt.Fprintf(s, "}")
+	return s.String()
+}
+
+// testBlocks are the first 10 pre-IDCT blocks from ../testdata/video-001.jpeg.
+var testBlocks = [10]block{
+	{
+		0x7f, 0xf6, 0x01, 0x07, 0xff, 0x00, 0x00, 0x00,
+		0xf5, 0x01, 0xfa, 0x01, 0xfe, 0x00, 0x01, 0x00,
+		0x05, 0x05, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0x01, 0xff, 0xf8, 0x00, 0x01, 0xff, 0x00, 0x00,
+		0x00, 0x01, 0x00, 0x01, 0x00, 0xff, 0xff, 0x00,
+		0xff, 0x0c, 0x00, 0x00, 0x00, 0x00, 0xff, 0x01,
+		0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00,
+		0x01, 0x00, 0x00, 0x01, 0xff, 0x01, 0x00, 0xfe,
+	},
+	{
+		0x29, 0x07, 0x00, 0xfc, 0x01, 0x01, 0x00, 0x00,
+		0x07, 0x00, 0x03, 0x00, 0x01, 0x00, 0xff, 0xff,
+		0xff, 0xfd, 0xff, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0x00, 0x00, 0x04, 0x00, 0xff, 0x01, 0x00, 0x00,
+		0x01, 0x00, 0x01, 0xff, 0x00, 0x00, 0x00, 0x00,
+		0x01, 0xfa, 0x01, 0x00, 0x01, 0x00, 0x01, 0xff,
+		0x00, 0x00, 0xff, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0x00, 0x00, 0x00, 0xff, 0x00, 0xff, 0x00, 0x02,
+	},
+	{
+		0xc5, 0xfa, 0x01, 0x00, 0x00, 0x01, 0x00, 0xff,
+		0x02, 0xff, 0x01, 0x00, 0x01, 0x00, 0xff, 0x00,
+		0xff, 0xff, 0x00, 0xff, 0x01, 0x00, 0x00, 0x00,
+		0xff, 0x00, 0x01, 0x00, 0x00, 0x00, 0xff, 0x00,
+		0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xff,
+		0x00, 0xff, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0xff, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+	},
+	{
+		0x86, 0x05, 0x00, 0x02, 0x00, 0x00, 0x01, 0x00,
+		0xf2, 0x06, 0x00, 0x00, 0x01, 0x02, 0x00, 0x00,
+		0xf6, 0xfa, 0xf9, 0x00, 0xff, 0x01, 0x00, 0x00,
+		0xf9, 0x00, 0x00, 0xff, 0x00, 0x00, 0x00, 0x00,
+		0x00, 0xff, 0x00, 0xff, 0xff, 0xff, 0x00, 0x00,
+		0xff, 0x00, 0x00, 0x01, 0x00, 0xff, 0x01, 0x00,
+		0x00, 0x00, 0x00, 0xff, 0x00, 0x00, 0x00, 0x01,
+		0x00, 0x01, 0xff, 0x01, 0x00, 0xff, 0x00, 0x00,
+	},
+	{
+		0x24, 0xfe, 0x00, 0xff, 0x00, 0xff, 0xff, 0x00,
+		0x08, 0xfd, 0x00, 0x01, 0x01, 0x00, 0x01, 0x00,
+		0x06, 0x03, 0x03, 0xff, 0x00, 0x00, 0x00, 0x00,
+		0x04, 0xff, 0x00, 0x00, 0x00, 0x00, 0x00, 0xff,
+		0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x01,
+		0x01, 0x00, 0x01, 0xff, 0x00, 0x01, 0x00, 0x00,
+		0x01, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0xff, 0x01,
+	},
+	{
+		0xcd, 0xff, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01,
+		0x03, 0xff, 0x00, 0x00, 0x00, 0x00, 0x00, 0xff,
+		0x01, 0x01, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00,
+		0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0x01, 0x00, 0x00, 0x00, 0x00, 0x01, 0x01, 0x00,
+		0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00,
+		0x00, 0xff, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0x00, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0xff,
+	},
+	{
+		0x81, 0xfe, 0x05, 0xff, 0x01, 0xff, 0x01, 0x00,
+		0xef, 0xf9, 0x00, 0xf9, 0x00, 0xff, 0x00, 0xff,
+		0x05, 0xf9, 0x00, 0xf8, 0x01, 0xff, 0x01, 0xff,
+		0x00, 0xff, 0x07, 0x00, 0x01, 0x00, 0x00, 0x00,
+		0x01, 0x00, 0x01, 0x01, 0x00, 0x00, 0x00, 0x00,
+		0x01, 0x00, 0x00, 0x00, 0xff, 0xff, 0x00, 0x01,
+		0xff, 0x01, 0x01, 0x00, 0xff, 0x00, 0x00, 0x00,
+		0x01, 0x01, 0x00, 0xff, 0x00, 0x00, 0x00, 0xff,
+	},
+	{
+		0x28, 0x00, 0xfe, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0x0b, 0x02, 0x01, 0x03, 0x00, 0xff, 0x00, 0x01,
+		0xfe, 0x02, 0x01, 0x03, 0xff, 0x00, 0x00, 0x00,
+		0x01, 0x00, 0xfd, 0x00, 0x01, 0x00, 0xff, 0x00,
+		0x01, 0xff, 0x00, 0xff, 0x01, 0x00, 0x00, 0x00,
+		0x00, 0x00, 0x00, 0xff, 0x01, 0x01, 0x00, 0xff,
+		0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+		0xff, 0xff, 0x00, 0x00, 0x00, 0xff, 0x00, 0x01,
+	},
+	{
+		0xdf, 0xf9, 0xfe, 0x00, 0x03, 0x01, 0xff, 0xff,
+		0x04, 0x01, 0x00, 0x01, 0x00, 0x00, 0x00, 0x00,
+		0xff, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x01,
+		0x00, 0x00, 0xfe, 0x01, 0x00, 0x00, 0x00, 0x00,
+		0x00, 0x00, 0xff, 0x01, 0x00, 0x00, 0x00, 0x01,
+		0xff, 0x00, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00,
+		0x00, 0xff, 0x00, 0xff, 0x01, 0x00, 0x00, 0x01,
+		0xff, 0xff, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00,
+	},
+	{
+		0x88, 0xfd, 0x00, 0x00, 0xff, 0x00, 0x01, 0xff,
+		0xe1, 0x06, 0x06, 0x01, 0xff, 0x00, 0x01, 0x00,
+		0x08, 0x00, 0xfa, 0x00, 0xff, 0xff, 0xff, 0xff,
+		0x08, 0x01, 0x00, 0xff, 0x01, 0xff, 0x00, 0x00,
+		0xf5, 0xff, 0x00, 0x01, 0xff, 0x01, 0x01, 0x00,
+		0xff, 0xff, 0x01, 0xff, 0x01, 0x00, 0x01, 0x00,
+		0x00, 0x01, 0x01, 0xff, 0x00, 0xff, 0x00, 0x01,
+		0x02, 0x00, 0x00, 0xff, 0xff, 0x00, 0xff, 0x00,
+	},
+}
diff --git a/src/image/jpeg/fdct.go b/src/image/jpeg/fdct.go
new file mode 100644
index 0000000..c7a973e
--- /dev/null
+++ b/src/image/jpeg/fdct.go
@@ -0,0 +1,192 @@
+// Copyright 2011 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 jpeg
+
+// This file implements a Forward Discrete Cosine Transformation.
+
+/*
+It is based on the code in jfdctint.c from the Independent JPEG Group,
+found at http://www.ijg.org/files/jpegsrc.v8c.tar.gz.
+
+The "LEGAL ISSUES" section of the README in that archive says:
+
+In plain English:
+
+1. We don't promise that this software works.  (But if you find any bugs,
+   please let us know!)
+2. You can use this software for whatever you want.  You don't have to pay us.
+3. You may not pretend that you wrote this software.  If you use it in a
+   program, you must acknowledge somewhere in your documentation that
+   you've used the IJG code.
+
+In legalese:
+
+The authors make NO WARRANTY or representation, either express or implied,
+with respect to this software, its quality, accuracy, merchantability, or
+fitness for a particular purpose.  This software is provided "AS IS", and you,
+its user, assume the entire risk as to its quality and accuracy.
+
+This software is copyright (C) 1991-2011, Thomas G. Lane, Guido Vollbeding.
+All Rights Reserved except as specified below.
+
+Permission is hereby granted to use, copy, modify, and distribute this
+software (or portions thereof) for any purpose, without fee, subject to these
+conditions:
+(1) If any part of the source code for this software is distributed, then this
+README file must be included, with this copyright and no-warranty notice
+unaltered; and any additions, deletions, or changes to the original files
+must be clearly indicated in accompanying documentation.
+(2) If only executable code is distributed, then the accompanying
+documentation must state that "this software is based in part on the work of
+the Independent JPEG Group".
+(3) Permission for use of this software is granted only if the user accepts
+full responsibility for any undesirable consequences; the authors accept
+NO LIABILITY for damages of any kind.
+
+These conditions apply to any software derived from or based on the IJG code,
+not just to the unmodified library.  If you use our work, you ought to
+acknowledge us.
+
+Permission is NOT granted for the use of any IJG author's name or company name
+in advertising or publicity relating to this software or products derived from
+it.  This software may be referred to only as "the Independent JPEG Group's
+software".
+
+We specifically permit and encourage the use of this software as the basis of
+commercial products, provided that all warranty or liability claims are
+assumed by the product vendor.
+*/
+
+// Trigonometric constants in 13-bit fixed point format.
+const (
+	fix_0_298631336 = 2446
+	fix_0_390180644 = 3196
+	fix_0_541196100 = 4433
+	fix_0_765366865 = 6270
+	fix_0_899976223 = 7373
+	fix_1_175875602 = 9633
+	fix_1_501321110 = 12299
+	fix_1_847759065 = 15137
+	fix_1_961570560 = 16069
+	fix_2_053119869 = 16819
+	fix_2_562915447 = 20995
+	fix_3_072711026 = 25172
+)
+
+const (
+	constBits     = 13
+	pass1Bits     = 2
+	centerJSample = 128
+)
+
+// fdct performs a forward DCT on an 8x8 block of coefficients, including a
+// level shift.
+func fdct(b *block) {
+	// Pass 1: process rows.
+	for y := 0; y < 8; y++ {
+		y8 := y * 8
+		s := b[y8 : y8+8 : y8+8] // Small cap improves performance, see https://golang.org/issue/27857
+		x0 := s[0]
+		x1 := s[1]
+		x2 := s[2]
+		x3 := s[3]
+		x4 := s[4]
+		x5 := s[5]
+		x6 := s[6]
+		x7 := s[7]
+
+		tmp0 := x0 + x7
+		tmp1 := x1 + x6
+		tmp2 := x2 + x5
+		tmp3 := x3 + x4
+
+		tmp10 := tmp0 + tmp3
+		tmp12 := tmp0 - tmp3
+		tmp11 := tmp1 + tmp2
+		tmp13 := tmp1 - tmp2
+
+		tmp0 = x0 - x7
+		tmp1 = x1 - x6
+		tmp2 = x2 - x5
+		tmp3 = x3 - x4
+
+		s[0] = (tmp10 + tmp11 - 8*centerJSample) << pass1Bits
+		s[4] = (tmp10 - tmp11) << pass1Bits
+		z1 := (tmp12 + tmp13) * fix_0_541196100
+		z1 += 1 << (constBits - pass1Bits - 1)
+		s[2] = (z1 + tmp12*fix_0_765366865) >> (constBits - pass1Bits)
+		s[6] = (z1 - tmp13*fix_1_847759065) >> (constBits - pass1Bits)
+
+		tmp10 = tmp0 + tmp3
+		tmp11 = tmp1 + tmp2
+		tmp12 = tmp0 + tmp2
+		tmp13 = tmp1 + tmp3
+		z1 = (tmp12 + tmp13) * fix_1_175875602
+		z1 += 1 << (constBits - pass1Bits - 1)
+		tmp0 *= fix_1_501321110
+		tmp1 *= fix_3_072711026
+		tmp2 *= fix_2_053119869
+		tmp3 *= fix_0_298631336
+		tmp10 *= -fix_0_899976223
+		tmp11 *= -fix_2_562915447
+		tmp12 *= -fix_0_390180644
+		tmp13 *= -fix_1_961570560
+
+		tmp12 += z1
+		tmp13 += z1
+		s[1] = (tmp0 + tmp10 + tmp12) >> (constBits - pass1Bits)
+		s[3] = (tmp1 + tmp11 + tmp13) >> (constBits - pass1Bits)
+		s[5] = (tmp2 + tmp11 + tmp12) >> (constBits - pass1Bits)
+		s[7] = (tmp3 + tmp10 + tmp13) >> (constBits - pass1Bits)
+	}
+	// Pass 2: process columns.
+	// We remove pass1Bits scaling, but leave results scaled up by an overall factor of 8.
+	for x := 0; x < 8; x++ {
+		tmp0 := b[0*8+x] + b[7*8+x]
+		tmp1 := b[1*8+x] + b[6*8+x]
+		tmp2 := b[2*8+x] + b[5*8+x]
+		tmp3 := b[3*8+x] + b[4*8+x]
+
+		tmp10 := tmp0 + tmp3 + 1<<(pass1Bits-1)
+		tmp12 := tmp0 - tmp3
+		tmp11 := tmp1 + tmp2
+		tmp13 := tmp1 - tmp2
+
+		tmp0 = b[0*8+x] - b[7*8+x]
+		tmp1 = b[1*8+x] - b[6*8+x]
+		tmp2 = b[2*8+x] - b[5*8+x]
+		tmp3 = b[3*8+x] - b[4*8+x]
+
+		b[0*8+x] = (tmp10 + tmp11) >> pass1Bits
+		b[4*8+x] = (tmp10 - tmp11) >> pass1Bits
+
+		z1 := (tmp12 + tmp13) * fix_0_541196100
+		z1 += 1 << (constBits + pass1Bits - 1)
+		b[2*8+x] = (z1 + tmp12*fix_0_765366865) >> (constBits + pass1Bits)
+		b[6*8+x] = (z1 - tmp13*fix_1_847759065) >> (constBits + pass1Bits)
+
+		tmp10 = tmp0 + tmp3
+		tmp11 = tmp1 + tmp2
+		tmp12 = tmp0 + tmp2
+		tmp13 = tmp1 + tmp3
+		z1 = (tmp12 + tmp13) * fix_1_175875602
+		z1 += 1 << (constBits + pass1Bits - 1)
+		tmp0 *= fix_1_501321110
+		tmp1 *= fix_3_072711026
+		tmp2 *= fix_2_053119869
+		tmp3 *= fix_0_298631336
+		tmp10 *= -fix_0_899976223
+		tmp11 *= -fix_2_562915447
+		tmp12 *= -fix_0_390180644
+		tmp13 *= -fix_1_961570560
+
+		tmp12 += z1
+		tmp13 += z1
+		b[1*8+x] = (tmp0 + tmp10 + tmp12) >> (constBits + pass1Bits)
+		b[3*8+x] = (tmp1 + tmp11 + tmp13) >> (constBits + pass1Bits)
+		b[5*8+x] = (tmp2 + tmp11 + tmp12) >> (constBits + pass1Bits)
+		b[7*8+x] = (tmp3 + tmp10 + tmp13) >> (constBits + pass1Bits)
+	}
+}
diff --git a/src/image/jpeg/fuzz_test.go b/src/image/jpeg/fuzz_test.go
new file mode 100644
index 0000000..bd534a9
--- /dev/null
+++ b/src/image/jpeg/fuzz_test.go
@@ -0,0 +1,65 @@
+// Copyright 2021 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 jpeg
+
+import (
+	"bytes"
+	"image"
+	"os"
+	"path/filepath"
+	"strings"
+	"testing"
+)
+
+func FuzzDecode(f *testing.F) {
+	if testing.Short() {
+		f.Skip("Skipping in short mode")
+	}
+
+	testdata, err := os.ReadDir("../testdata")
+	if err != nil {
+		f.Fatalf("failed to read testdata directory: %s", err)
+	}
+	for _, de := range testdata {
+		if de.IsDir() || !strings.HasSuffix(de.Name(), ".jpeg") {
+			continue
+		}
+		b, err := os.ReadFile(filepath.Join("../testdata", de.Name()))
+		if err != nil {
+			f.Fatalf("failed to read testdata: %s", err)
+		}
+		f.Add(b)
+	}
+
+	f.Fuzz(func(t *testing.T, b []byte) {
+		cfg, _, err := image.DecodeConfig(bytes.NewReader(b))
+		if err != nil {
+			return
+		}
+		if cfg.Width*cfg.Height > 1e6 {
+			return
+		}
+		img, typ, err := image.Decode(bytes.NewReader(b))
+		if err != nil || typ != "jpeg" {
+			return
+		}
+		for q := 1; q <= 100; q++ {
+			var w bytes.Buffer
+			err := Encode(&w, img, &Options{Quality: q})
+			if err != nil {
+				t.Fatalf("failed to encode valid image: %s", err)
+			}
+			img1, err := Decode(&w)
+			if err != nil {
+				t.Fatalf("failed to decode roundtripped image: %s", err)
+			}
+			got := img1.Bounds()
+			want := img.Bounds()
+			if !got.Eq(want) {
+				t.Fatalf("roundtripped image bounds have changed, got: %s, want: %s", got, want)
+			}
+		}
+	})
+}
diff --git a/src/image/jpeg/huffman.go b/src/image/jpeg/huffman.go
new file mode 100644
index 0000000..95aaf71
--- /dev/null
+++ b/src/image/jpeg/huffman.go
@@ -0,0 +1,247 @@
+// Copyright 2009 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 jpeg
+
+import (
+	"io"
+)
+
+// maxCodeLength is the maximum (inclusive) number of bits in a Huffman code.
+const maxCodeLength = 16
+
+// maxNCodes is the maximum (inclusive) number of codes in a Huffman tree.
+const maxNCodes = 256
+
+// lutSize is the log-2 size of the Huffman decoder's look-up table.
+const lutSize = 8
+
+// huffman is a Huffman decoder, specified in section C.
+type huffman struct {
+	// length is the number of codes in the tree.
+	nCodes int32
+	// lut is the look-up table for the next lutSize bits in the bit-stream.
+	// The high 8 bits of the uint16 are the encoded value. The low 8 bits
+	// are 1 plus the code length, or 0 if the value is too large to fit in
+	// lutSize bits.
+	lut [1 << lutSize]uint16
+	// vals are the decoded values, sorted by their encoding.
+	vals [maxNCodes]uint8
+	// minCodes[i] is the minimum code of length i, or -1 if there are no
+	// codes of that length.
+	minCodes [maxCodeLength]int32
+	// maxCodes[i] is the maximum code of length i, or -1 if there are no
+	// codes of that length.
+	maxCodes [maxCodeLength]int32
+	// valsIndices[i] is the index into vals of minCodes[i].
+	valsIndices [maxCodeLength]int32
+}
+
+// errShortHuffmanData means that an unexpected EOF occurred while decoding
+// Huffman data.
+var errShortHuffmanData = FormatError("short Huffman data")
+
+// ensureNBits reads bytes from the byte buffer to ensure that d.bits.n is at
+// least n. For best performance (avoiding function calls inside hot loops),
+// the caller is the one responsible for first checking that d.bits.n < n.
+func (d *decoder) ensureNBits(n int32) error {
+	for {
+		c, err := d.readByteStuffedByte()
+		if err != nil {
+			if err == io.EOF {
+				return errShortHuffmanData
+			}
+			return err
+		}
+		d.bits.a = d.bits.a<<8 | uint32(c)
+		d.bits.n += 8
+		if d.bits.m == 0 {
+			d.bits.m = 1 << 7
+		} else {
+			d.bits.m <<= 8
+		}
+		if d.bits.n >= n {
+			break
+		}
+	}
+	return nil
+}
+
+// receiveExtend is the composition of RECEIVE and EXTEND, specified in section
+// F.2.2.1.
+func (d *decoder) receiveExtend(t uint8) (int32, error) {
+	if d.bits.n < int32(t) {
+		if err := d.ensureNBits(int32(t)); err != nil {
+			return 0, err
+		}
+	}
+	d.bits.n -= int32(t)
+	d.bits.m >>= t
+	s := int32(1) << t
+	x := int32(d.bits.a>>uint8(d.bits.n)) & (s - 1)
+	if x < s>>1 {
+		x += ((-1) << t) + 1
+	}
+	return x, nil
+}
+
+// processDHT processes a Define Huffman Table marker, and initializes a huffman
+// struct from its contents. Specified in section B.2.4.2.
+func (d *decoder) processDHT(n int) error {
+	for n > 0 {
+		if n < 17 {
+			return FormatError("DHT has wrong length")
+		}
+		if err := d.readFull(d.tmp[:17]); err != nil {
+			return err
+		}
+		tc := d.tmp[0] >> 4
+		if tc > maxTc {
+			return FormatError("bad Tc value")
+		}
+		th := d.tmp[0] & 0x0f
+		// The baseline th <= 1 restriction is specified in table B.5.
+		if th > maxTh || (d.baseline && th > 1) {
+			return FormatError("bad Th value")
+		}
+		h := &d.huff[tc][th]
+
+		// Read nCodes and h.vals (and derive h.nCodes).
+		// nCodes[i] is the number of codes with code length i.
+		// h.nCodes is the total number of codes.
+		h.nCodes = 0
+		var nCodes [maxCodeLength]int32
+		for i := range nCodes {
+			nCodes[i] = int32(d.tmp[i+1])
+			h.nCodes += nCodes[i]
+		}
+		if h.nCodes == 0 {
+			return FormatError("Huffman table has zero length")
+		}
+		if h.nCodes > maxNCodes {
+			return FormatError("Huffman table has excessive length")
+		}
+		n -= int(h.nCodes) + 17
+		if n < 0 {
+			return FormatError("DHT has wrong length")
+		}
+		if err := d.readFull(h.vals[:h.nCodes]); err != nil {
+			return err
+		}
+
+		// Derive the look-up table.
+		for i := range h.lut {
+			h.lut[i] = 0
+		}
+		var x, code uint32
+		for i := uint32(0); i < lutSize; i++ {
+			code <<= 1
+			for j := int32(0); j < nCodes[i]; j++ {
+				// The codeLength is 1+i, so shift code by 8-(1+i) to
+				// calculate the high bits for every 8-bit sequence
+				// whose codeLength's high bits matches code.
+				// The high 8 bits of lutValue are the encoded value.
+				// The low 8 bits are 1 plus the codeLength.
+				base := uint8(code << (7 - i))
+				lutValue := uint16(h.vals[x])<<8 | uint16(2+i)
+				for k := uint8(0); k < 1<<(7-i); k++ {
+					h.lut[base|k] = lutValue
+				}
+				code++
+				x++
+			}
+		}
+
+		// Derive minCodes, maxCodes, and valsIndices.
+		var c, index int32
+		for i, n := range nCodes {
+			if n == 0 {
+				h.minCodes[i] = -1
+				h.maxCodes[i] = -1
+				h.valsIndices[i] = -1
+			} else {
+				h.minCodes[i] = c
+				h.maxCodes[i] = c + n - 1
+				h.valsIndices[i] = index
+				c += n
+				index += n
+			}
+			c <<= 1
+		}
+	}
+	return nil
+}
+
+// decodeHuffman returns the next Huffman-coded value from the bit-stream,
+// decoded according to h.
+func (d *decoder) decodeHuffman(h *huffman) (uint8, error) {
+	if h.nCodes == 0 {
+		return 0, FormatError("uninitialized Huffman table")
+	}
+
+	if d.bits.n < 8 {
+		if err := d.ensureNBits(8); err != nil {
+			if err != errMissingFF00 && err != errShortHuffmanData {
+				return 0, err
+			}
+			// There are no more bytes of data in this segment, but we may still
+			// be able to read the next symbol out of the previously read bits.
+			// First, undo the readByte that the ensureNBits call made.
+			if d.bytes.nUnreadable != 0 {
+				d.unreadByteStuffedByte()
+			}
+			goto slowPath
+		}
+	}
+	if v := h.lut[(d.bits.a>>uint32(d.bits.n-lutSize))&0xff]; v != 0 {
+		n := (v & 0xff) - 1
+		d.bits.n -= int32(n)
+		d.bits.m >>= n
+		return uint8(v >> 8), nil
+	}
+
+slowPath:
+	for i, code := 0, int32(0); i < maxCodeLength; i++ {
+		if d.bits.n == 0 {
+			if err := d.ensureNBits(1); err != nil {
+				return 0, err
+			}
+		}
+		if d.bits.a&d.bits.m != 0 {
+			code |= 1
+		}
+		d.bits.n--
+		d.bits.m >>= 1
+		if code <= h.maxCodes[i] {
+			return h.vals[h.valsIndices[i]+code-h.minCodes[i]], nil
+		}
+		code <<= 1
+	}
+	return 0, FormatError("bad Huffman code")
+}
+
+func (d *decoder) decodeBit() (bool, error) {
+	if d.bits.n == 0 {
+		if err := d.ensureNBits(1); err != nil {
+			return false, err
+		}
+	}
+	ret := d.bits.a&d.bits.m != 0
+	d.bits.n--
+	d.bits.m >>= 1
+	return ret, nil
+}
+
+func (d *decoder) decodeBits(n int32) (uint32, error) {
+	if d.bits.n < n {
+		if err := d.ensureNBits(n); err != nil {
+			return 0, err
+		}
+	}
+	ret := d.bits.a >> uint32(d.bits.n-n)
+	ret &= (1 << uint32(n)) - 1
+	d.bits.n -= n
+	d.bits.m >>= uint32(n)
+	return ret, nil
+}
diff --git a/src/image/jpeg/idct.go b/src/image/jpeg/idct.go
new file mode 100644
index 0000000..a3957c8
--- /dev/null
+++ b/src/image/jpeg/idct.go
@@ -0,0 +1,194 @@
+// Copyright 2009 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 jpeg
+
+// This is a Go translation of idct.c from
+//
+// http://standards.iso.org/ittf/PubliclyAvailableStandards/ISO_IEC_13818-4_2004_Conformance_Testing/Video/verifier/mpeg2decode_960109.tar.gz
+//
+// which carries the following notice:
+
+/* Copyright (C) 1996, MPEG Software Simulation Group. All Rights Reserved. */
+
+/*
+ * Disclaimer of Warranty
+ *
+ * These software programs are available to the user without any license fee or
+ * royalty on an "as is" basis.  The MPEG Software Simulation Group disclaims
+ * any and all warranties, whether express, implied, or statuary, including any
+ * implied warranties or merchantability or of fitness for a particular
+ * purpose.  In no event shall the copyright-holder be liable for any
+ * incidental, punitive, or consequential damages of any kind whatsoever
+ * arising from the use of these programs.
+ *
+ * This disclaimer of warranty extends to the user of these programs and user's
+ * customers, employees, agents, transferees, successors, and assigns.
+ *
+ * The MPEG Software Simulation Group does not represent or warrant that the
+ * programs furnished hereunder are free of infringement of any third-party
+ * patents.
+ *
+ * Commercial implementations of MPEG-1 and MPEG-2 video, including shareware,
+ * are subject to royalty fees to patent holders.  Many of these patents are
+ * general enough such that they are unavoidable regardless of implementation
+ * design.
+ *
+ */
+
+const blockSize = 64 // A DCT block is 8x8.
+
+type block [blockSize]int32
+
+const (
+	w1 = 2841 // 2048*sqrt(2)*cos(1*pi/16)
+	w2 = 2676 // 2048*sqrt(2)*cos(2*pi/16)
+	w3 = 2408 // 2048*sqrt(2)*cos(3*pi/16)
+	w5 = 1609 // 2048*sqrt(2)*cos(5*pi/16)
+	w6 = 1108 // 2048*sqrt(2)*cos(6*pi/16)
+	w7 = 565  // 2048*sqrt(2)*cos(7*pi/16)
+
+	w1pw7 = w1 + w7
+	w1mw7 = w1 - w7
+	w2pw6 = w2 + w6
+	w2mw6 = w2 - w6
+	w3pw5 = w3 + w5
+	w3mw5 = w3 - w5
+
+	r2 = 181 // 256/sqrt(2)
+)
+
+// idct performs a 2-D Inverse Discrete Cosine Transformation.
+//
+// The input coefficients should already have been multiplied by the
+// appropriate quantization table. We use fixed-point computation, with the
+// number of bits for the fractional component varying over the intermediate
+// stages.
+//
+// For more on the actual algorithm, see Z. Wang, "Fast algorithms for the
+// discrete W transform and for the discrete Fourier transform", IEEE Trans. on
+// ASSP, Vol. ASSP- 32, pp. 803-816, Aug. 1984.
+func idct(src *block) {
+	// Horizontal 1-D IDCT.
+	for y := 0; y < 8; y++ {
+		y8 := y * 8
+		s := src[y8 : y8+8 : y8+8] // Small cap improves performance, see https://golang.org/issue/27857
+		// If all the AC components are zero, then the IDCT is trivial.
+		if s[1] == 0 && s[2] == 0 && s[3] == 0 &&
+			s[4] == 0 && s[5] == 0 && s[6] == 0 && s[7] == 0 {
+			dc := s[0] << 3
+			s[0] = dc
+			s[1] = dc
+			s[2] = dc
+			s[3] = dc
+			s[4] = dc
+			s[5] = dc
+			s[6] = dc
+			s[7] = dc
+			continue
+		}
+
+		// Prescale.
+		x0 := (s[0] << 11) + 128
+		x1 := s[4] << 11
+		x2 := s[6]
+		x3 := s[2]
+		x4 := s[1]
+		x5 := s[7]
+		x6 := s[5]
+		x7 := s[3]
+
+		// Stage 1.
+		x8 := w7 * (x4 + x5)
+		x4 = x8 + w1mw7*x4
+		x5 = x8 - w1pw7*x5
+		x8 = w3 * (x6 + x7)
+		x6 = x8 - w3mw5*x6
+		x7 = x8 - w3pw5*x7
+
+		// Stage 2.
+		x8 = x0 + x1
+		x0 -= x1
+		x1 = w6 * (x3 + x2)
+		x2 = x1 - w2pw6*x2
+		x3 = x1 + w2mw6*x3
+		x1 = x4 + x6
+		x4 -= x6
+		x6 = x5 + x7
+		x5 -= x7
+
+		// Stage 3.
+		x7 = x8 + x3
+		x8 -= x3
+		x3 = x0 + x2
+		x0 -= x2
+		x2 = (r2*(x4+x5) + 128) >> 8
+		x4 = (r2*(x4-x5) + 128) >> 8
+
+		// Stage 4.
+		s[0] = (x7 + x1) >> 8
+		s[1] = (x3 + x2) >> 8
+		s[2] = (x0 + x4) >> 8
+		s[3] = (x8 + x6) >> 8
+		s[4] = (x8 - x6) >> 8
+		s[5] = (x0 - x4) >> 8
+		s[6] = (x3 - x2) >> 8
+		s[7] = (x7 - x1) >> 8
+	}
+
+	// Vertical 1-D IDCT.
+	for x := 0; x < 8; x++ {
+		// Similar to the horizontal 1-D IDCT case, if all the AC components are zero, then the IDCT is trivial.
+		// However, after performing the horizontal 1-D IDCT, there are typically non-zero AC components, so
+		// we do not bother to check for the all-zero case.
+		s := src[x : x+57 : x+57] // Small cap improves performance, see https://golang.org/issue/27857
+
+		// Prescale.
+		y0 := (s[8*0] << 8) + 8192
+		y1 := s[8*4] << 8
+		y2 := s[8*6]
+		y3 := s[8*2]
+		y4 := s[8*1]
+		y5 := s[8*7]
+		y6 := s[8*5]
+		y7 := s[8*3]
+
+		// Stage 1.
+		y8 := w7*(y4+y5) + 4
+		y4 = (y8 + w1mw7*y4) >> 3
+		y5 = (y8 - w1pw7*y5) >> 3
+		y8 = w3*(y6+y7) + 4
+		y6 = (y8 - w3mw5*y6) >> 3
+		y7 = (y8 - w3pw5*y7) >> 3
+
+		// Stage 2.
+		y8 = y0 + y1
+		y0 -= y1
+		y1 = w6*(y3+y2) + 4
+		y2 = (y1 - w2pw6*y2) >> 3
+		y3 = (y1 + w2mw6*y3) >> 3
+		y1 = y4 + y6
+		y4 -= y6
+		y6 = y5 + y7
+		y5 -= y7
+
+		// Stage 3.
+		y7 = y8 + y3
+		y8 -= y3
+		y3 = y0 + y2
+		y0 -= y2
+		y2 = (r2*(y4+y5) + 128) >> 8
+		y4 = (r2*(y4-y5) + 128) >> 8
+
+		// Stage 4.
+		s[8*0] = (y7 + y1) >> 14
+		s[8*1] = (y3 + y2) >> 14
+		s[8*2] = (y0 + y4) >> 14
+		s[8*3] = (y8 + y6) >> 14
+		s[8*4] = (y8 - y6) >> 14
+		s[8*5] = (y0 - y4) >> 14
+		s[8*6] = (y3 - y2) >> 14
+		s[8*7] = (y7 - y1) >> 14
+	}
+}
diff --git a/src/image/jpeg/reader.go b/src/image/jpeg/reader.go
new file mode 100644
index 0000000..b340723
--- /dev/null
+++ b/src/image/jpeg/reader.go
@@ -0,0 +1,815 @@
+// Copyright 2009 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 jpeg implements a JPEG image decoder and encoder.
+//
+// JPEG is defined in ITU-T T.81: https://www.w3.org/Graphics/JPEG/itu-t81.pdf.
+package jpeg
+
+import (
+	"image"
+	"image/color"
+	"image/internal/imageutil"
+	"io"
+)
+
+// A FormatError reports that the input is not a valid JPEG.
+type FormatError string
+
+func (e FormatError) Error() string { return "invalid JPEG format: " + string(e) }
+
+// An UnsupportedError reports that the input uses a valid but unimplemented JPEG feature.
+type UnsupportedError string
+
+func (e UnsupportedError) Error() string { return "unsupported JPEG feature: " + string(e) }
+
+var errUnsupportedSubsamplingRatio = UnsupportedError("luma/chroma subsampling ratio")
+
+// Component specification, specified in section B.2.2.
+type component struct {
+	h  int   // Horizontal sampling factor.
+	v  int   // Vertical sampling factor.
+	c  uint8 // Component identifier.
+	tq uint8 // Quantization table destination selector.
+}
+
+const (
+	dcTable = 0
+	acTable = 1
+	maxTc   = 1
+	maxTh   = 3
+	maxTq   = 3
+
+	maxComponents = 4
+)
+
+const (
+	sof0Marker = 0xc0 // Start Of Frame (Baseline Sequential).
+	sof1Marker = 0xc1 // Start Of Frame (Extended Sequential).
+	sof2Marker = 0xc2 // Start Of Frame (Progressive).
+	dhtMarker  = 0xc4 // Define Huffman Table.
+	rst0Marker = 0xd0 // ReSTart (0).
+	rst7Marker = 0xd7 // ReSTart (7).
+	soiMarker  = 0xd8 // Start Of Image.
+	eoiMarker  = 0xd9 // End Of Image.
+	sosMarker  = 0xda // Start Of Scan.
+	dqtMarker  = 0xdb // Define Quantization Table.
+	driMarker  = 0xdd // Define Restart Interval.
+	comMarker  = 0xfe // COMment.
+	// "APPlication specific" markers aren't part of the JPEG spec per se,
+	// but in practice, their use is described at
+	// https://www.sno.phy.queensu.ca/~phil/exiftool/TagNames/JPEG.html
+	app0Marker  = 0xe0
+	app14Marker = 0xee
+	app15Marker = 0xef
+)
+
+// See https://www.sno.phy.queensu.ca/~phil/exiftool/TagNames/JPEG.html#Adobe
+const (
+	adobeTransformUnknown = 0
+	adobeTransformYCbCr   = 1
+	adobeTransformYCbCrK  = 2
+)
+
+// unzig maps from the zig-zag ordering to the natural ordering. For example,
+// unzig[3] is the column and row of the fourth element in zig-zag order. The
+// value is 16, which means first column (16%8 == 0) and third row (16/8 == 2).
+var unzig = [blockSize]int{
+	0, 1, 8, 16, 9, 2, 3, 10,
+	17, 24, 32, 25, 18, 11, 4, 5,
+	12, 19, 26, 33, 40, 48, 41, 34,
+	27, 20, 13, 6, 7, 14, 21, 28,
+	35, 42, 49, 56, 57, 50, 43, 36,
+	29, 22, 15, 23, 30, 37, 44, 51,
+	58, 59, 52, 45, 38, 31, 39, 46,
+	53, 60, 61, 54, 47, 55, 62, 63,
+}
+
+// Deprecated: Reader is not used by the image/jpeg package and should
+// not be used by others. It is kept for compatibility.
+type Reader interface {
+	io.ByteReader
+	io.Reader
+}
+
+// bits holds the unprocessed bits that have been taken from the byte-stream.
+// The n least significant bits of a form the unread bits, to be read in MSB to
+// LSB order.
+type bits struct {
+	a uint32 // accumulator.
+	m uint32 // mask. m==1<<(n-1) when n>0, with m==0 when n==0.
+	n int32  // the number of unread bits in a.
+}
+
+type decoder struct {
+	r    io.Reader
+	bits bits
+	// bytes is a byte buffer, similar to a bufio.Reader, except that it
+	// has to be able to unread more than 1 byte, due to byte stuffing.
+	// Byte stuffing is specified in section F.1.2.3.
+	bytes struct {
+		// buf[i:j] are the buffered bytes read from the underlying
+		// io.Reader that haven't yet been passed further on.
+		buf  [4096]byte
+		i, j int
+		// nUnreadable is the number of bytes to back up i after
+		// overshooting. It can be 0, 1 or 2.
+		nUnreadable int
+	}
+	width, height int
+
+	img1        *image.Gray
+	img3        *image.YCbCr
+	blackPix    []byte
+	blackStride int
+
+	ri    int // Restart Interval.
+	nComp int
+
+	// As per section 4.5, there are four modes of operation (selected by the
+	// SOF? markers): sequential DCT, progressive DCT, lossless and
+	// hierarchical, although this implementation does not support the latter
+	// two non-DCT modes. Sequential DCT is further split into baseline and
+	// extended, as per section 4.11.
+	baseline    bool
+	progressive bool
+
+	jfif                bool
+	adobeTransformValid bool
+	adobeTransform      uint8
+	eobRun              uint16 // End-of-Band run, specified in section G.1.2.2.
+
+	comp       [maxComponents]component
+	progCoeffs [maxComponents][]block // Saved state between progressive-mode scans.
+	huff       [maxTc + 1][maxTh + 1]huffman
+	quant      [maxTq + 1]block // Quantization tables, in zig-zag order.
+	tmp        [2 * blockSize]byte
+}
+
+// fill fills up the d.bytes.buf buffer from the underlying io.Reader. It
+// should only be called when there are no unread bytes in d.bytes.
+func (d *decoder) fill() error {
+	if d.bytes.i != d.bytes.j {
+		panic("jpeg: fill called when unread bytes exist")
+	}
+	// Move the last 2 bytes to the start of the buffer, in case we need
+	// to call unreadByteStuffedByte.
+	if d.bytes.j > 2 {
+		d.bytes.buf[0] = d.bytes.buf[d.bytes.j-2]
+		d.bytes.buf[1] = d.bytes.buf[d.bytes.j-1]
+		d.bytes.i, d.bytes.j = 2, 2
+	}
+	// Fill in the rest of the buffer.
+	n, err := d.r.Read(d.bytes.buf[d.bytes.j:])
+	d.bytes.j += n
+	if n > 0 {
+		err = nil
+	}
+	return err
+}
+
+// unreadByteStuffedByte undoes the most recent readByteStuffedByte call,
+// giving a byte of data back from d.bits to d.bytes. The Huffman look-up table
+// requires at least 8 bits for look-up, which means that Huffman decoding can
+// sometimes overshoot and read one or two too many bytes. Two-byte overshoot
+// can happen when expecting to read a 0xff 0x00 byte-stuffed byte.
+func (d *decoder) unreadByteStuffedByte() {
+	d.bytes.i -= d.bytes.nUnreadable
+	d.bytes.nUnreadable = 0
+	if d.bits.n >= 8 {
+		d.bits.a >>= 8
+		d.bits.n -= 8
+		d.bits.m >>= 8
+	}
+}
+
+// readByte returns the next byte, whether buffered or not buffered. It does
+// not care about byte stuffing.
+func (d *decoder) readByte() (x byte, err error) {
+	for d.bytes.i == d.bytes.j {
+		if err = d.fill(); err != nil {
+			return 0, err
+		}
+	}
+	x = d.bytes.buf[d.bytes.i]
+	d.bytes.i++
+	d.bytes.nUnreadable = 0
+	return x, nil
+}
+
+// errMissingFF00 means that readByteStuffedByte encountered an 0xff byte (a
+// marker byte) that wasn't the expected byte-stuffed sequence 0xff, 0x00.
+var errMissingFF00 = FormatError("missing 0xff00 sequence")
+
+// readByteStuffedByte is like readByte but is for byte-stuffed Huffman data.
+func (d *decoder) readByteStuffedByte() (x byte, err error) {
+	// Take the fast path if d.bytes.buf contains at least two bytes.
+	if d.bytes.i+2 <= d.bytes.j {
+		x = d.bytes.buf[d.bytes.i]
+		d.bytes.i++
+		d.bytes.nUnreadable = 1
+		if x != 0xff {
+			return x, err
+		}
+		if d.bytes.buf[d.bytes.i] != 0x00 {
+			return 0, errMissingFF00
+		}
+		d.bytes.i++
+		d.bytes.nUnreadable = 2
+		return 0xff, nil
+	}
+
+	d.bytes.nUnreadable = 0
+
+	x, err = d.readByte()
+	if err != nil {
+		return 0, err
+	}
+	d.bytes.nUnreadable = 1
+	if x != 0xff {
+		return x, nil
+	}
+
+	x, err = d.readByte()
+	if err != nil {
+		return 0, err
+	}
+	d.bytes.nUnreadable = 2
+	if x != 0x00 {
+		return 0, errMissingFF00
+	}
+	return 0xff, nil
+}
+
+// readFull reads exactly len(p) bytes into p. It does not care about byte
+// stuffing.
+func (d *decoder) readFull(p []byte) error {
+	// Unread the overshot bytes, if any.
+	if d.bytes.nUnreadable != 0 {
+		if d.bits.n >= 8 {
+			d.unreadByteStuffedByte()
+		}
+		d.bytes.nUnreadable = 0
+	}
+
+	for {
+		n := copy(p, d.bytes.buf[d.bytes.i:d.bytes.j])
+		p = p[n:]
+		d.bytes.i += n
+		if len(p) == 0 {
+			break
+		}
+		if err := d.fill(); err != nil {
+			if err == io.EOF {
+				err = io.ErrUnexpectedEOF
+			}
+			return err
+		}
+	}
+	return nil
+}
+
+// ignore ignores the next n bytes.
+func (d *decoder) ignore(n int) error {
+	// Unread the overshot bytes, if any.
+	if d.bytes.nUnreadable != 0 {
+		if d.bits.n >= 8 {
+			d.unreadByteStuffedByte()
+		}
+		d.bytes.nUnreadable = 0
+	}
+
+	for {
+		m := d.bytes.j - d.bytes.i
+		if m > n {
+			m = n
+		}
+		d.bytes.i += m
+		n -= m
+		if n == 0 {
+			break
+		}
+		if err := d.fill(); err != nil {
+			if err == io.EOF {
+				err = io.ErrUnexpectedEOF
+			}
+			return err
+		}
+	}
+	return nil
+}
+
+// Specified in section B.2.2.
+func (d *decoder) processSOF(n int) error {
+	if d.nComp != 0 {
+		return FormatError("multiple SOF markers")
+	}
+	switch n {
+	case 6 + 3*1: // Grayscale image.
+		d.nComp = 1
+	case 6 + 3*3: // YCbCr or RGB image.
+		d.nComp = 3
+	case 6 + 3*4: // YCbCrK or CMYK image.
+		d.nComp = 4
+	default:
+		return UnsupportedError("number of components")
+	}
+	if err := d.readFull(d.tmp[:n]); err != nil {
+		return err
+	}
+	// We only support 8-bit precision.
+	if d.tmp[0] != 8 {
+		return UnsupportedError("precision")
+	}
+	d.height = int(d.tmp[1])<<8 + int(d.tmp[2])
+	d.width = int(d.tmp[3])<<8 + int(d.tmp[4])
+	if int(d.tmp[5]) != d.nComp {
+		return FormatError("SOF has wrong length")
+	}
+
+	for i := 0; i < d.nComp; i++ {
+		d.comp[i].c = d.tmp[6+3*i]
+		// Section B.2.2 states that "the value of C_i shall be different from
+		// the values of C_1 through C_(i-1)".
+		for j := 0; j < i; j++ {
+			if d.comp[i].c == d.comp[j].c {
+				return FormatError("repeated component identifier")
+			}
+		}
+
+		d.comp[i].tq = d.tmp[8+3*i]
+		if d.comp[i].tq > maxTq {
+			return FormatError("bad Tq value")
+		}
+
+		hv := d.tmp[7+3*i]
+		h, v := int(hv>>4), int(hv&0x0f)
+		if h < 1 || 4 < h || v < 1 || 4 < v {
+			return FormatError("luma/chroma subsampling ratio")
+		}
+		if h == 3 || v == 3 {
+			return errUnsupportedSubsamplingRatio
+		}
+		switch d.nComp {
+		case 1:
+			// If a JPEG image has only one component, section A.2 says "this data
+			// is non-interleaved by definition" and section A.2.2 says "[in this
+			// case...] the order of data units within a scan shall be left-to-right
+			// and top-to-bottom... regardless of the values of H_1 and V_1". Section
+			// 4.8.2 also says "[for non-interleaved data], the MCU is defined to be
+			// one data unit". Similarly, section A.1.1 explains that it is the ratio
+			// of H_i to max_j(H_j) that matters, and similarly for V. For grayscale
+			// images, H_1 is the maximum H_j for all components j, so that ratio is
+			// always 1. The component's (h, v) is effectively always (1, 1): even if
+			// the nominal (h, v) is (2, 1), a 20x5 image is encoded in three 8x8
+			// MCUs, not two 16x8 MCUs.
+			h, v = 1, 1
+
+		case 3:
+			// For YCbCr images, we only support 4:4:4, 4:4:0, 4:2:2, 4:2:0,
+			// 4:1:1 or 4:1:0 chroma subsampling ratios. This implies that the
+			// (h, v) values for the Y component are either (1, 1), (1, 2),
+			// (2, 1), (2, 2), (4, 1) or (4, 2), and the Y component's values
+			// must be a multiple of the Cb and Cr component's values. We also
+			// assume that the two chroma components have the same subsampling
+			// ratio.
+			switch i {
+			case 0: // Y.
+				// We have already verified, above, that h and v are both
+				// either 1, 2 or 4, so invalid (h, v) combinations are those
+				// with v == 4.
+				if v == 4 {
+					return errUnsupportedSubsamplingRatio
+				}
+			case 1: // Cb.
+				if d.comp[0].h%h != 0 || d.comp[0].v%v != 0 {
+					return errUnsupportedSubsamplingRatio
+				}
+			case 2: // Cr.
+				if d.comp[1].h != h || d.comp[1].v != v {
+					return errUnsupportedSubsamplingRatio
+				}
+			}
+
+		case 4:
+			// For 4-component images (either CMYK or YCbCrK), we only support two
+			// hv vectors: [0x11 0x11 0x11 0x11] and [0x22 0x11 0x11 0x22].
+			// Theoretically, 4-component JPEG images could mix and match hv values
+			// but in practice, those two combinations are the only ones in use,
+			// and it simplifies the applyBlack code below if we can assume that:
+			//	- for CMYK, the C and K channels have full samples, and if the M
+			//	  and Y channels subsample, they subsample both horizontally and
+			//	  vertically.
+			//	- for YCbCrK, the Y and K channels have full samples.
+			switch i {
+			case 0:
+				if hv != 0x11 && hv != 0x22 {
+					return errUnsupportedSubsamplingRatio
+				}
+			case 1, 2:
+				if hv != 0x11 {
+					return errUnsupportedSubsamplingRatio
+				}
+			case 3:
+				if d.comp[0].h != h || d.comp[0].v != v {
+					return errUnsupportedSubsamplingRatio
+				}
+			}
+		}
+
+		d.comp[i].h = h
+		d.comp[i].v = v
+	}
+	return nil
+}
+
+// Specified in section B.2.4.1.
+func (d *decoder) processDQT(n int) error {
+loop:
+	for n > 0 {
+		n--
+		x, err := d.readByte()
+		if err != nil {
+			return err
+		}
+		tq := x & 0x0f
+		if tq > maxTq {
+			return FormatError("bad Tq value")
+		}
+		switch x >> 4 {
+		default:
+			return FormatError("bad Pq value")
+		case 0:
+			if n < blockSize {
+				break loop
+			}
+			n -= blockSize
+			if err := d.readFull(d.tmp[:blockSize]); err != nil {
+				return err
+			}
+			for i := range d.quant[tq] {
+				d.quant[tq][i] = int32(d.tmp[i])
+			}
+		case 1:
+			if n < 2*blockSize {
+				break loop
+			}
+			n -= 2 * blockSize
+			if err := d.readFull(d.tmp[:2*blockSize]); err != nil {
+				return err
+			}
+			for i := range d.quant[tq] {
+				d.quant[tq][i] = int32(d.tmp[2*i])<<8 | int32(d.tmp[2*i+1])
+			}
+		}
+	}
+	if n != 0 {
+		return FormatError("DQT has wrong length")
+	}
+	return nil
+}
+
+// Specified in section B.2.4.4.
+func (d *decoder) processDRI(n int) error {
+	if n != 2 {
+		return FormatError("DRI has wrong length")
+	}
+	if err := d.readFull(d.tmp[:2]); err != nil {
+		return err
+	}
+	d.ri = int(d.tmp[0])<<8 + int(d.tmp[1])
+	return nil
+}
+
+func (d *decoder) processApp0Marker(n int) error {
+	if n < 5 {
+		return d.ignore(n)
+	}
+	if err := d.readFull(d.tmp[:5]); err != nil {
+		return err
+	}
+	n -= 5
+
+	d.jfif = d.tmp[0] == 'J' && d.tmp[1] == 'F' && d.tmp[2] == 'I' && d.tmp[3] == 'F' && d.tmp[4] == '\x00'
+
+	if n > 0 {
+		return d.ignore(n)
+	}
+	return nil
+}
+
+func (d *decoder) processApp14Marker(n int) error {
+	if n < 12 {
+		return d.ignore(n)
+	}
+	if err := d.readFull(d.tmp[:12]); err != nil {
+		return err
+	}
+	n -= 12
+
+	if d.tmp[0] == 'A' && d.tmp[1] == 'd' && d.tmp[2] == 'o' && d.tmp[3] == 'b' && d.tmp[4] == 'e' {
+		d.adobeTransformValid = true
+		d.adobeTransform = d.tmp[11]
+	}
+
+	if n > 0 {
+		return d.ignore(n)
+	}
+	return nil
+}
+
+// decode reads a JPEG image from r and returns it as an image.Image.
+func (d *decoder) decode(r io.Reader, configOnly bool) (image.Image, error) {
+	d.r = r
+
+	// Check for the Start Of Image marker.
+	if err := d.readFull(d.tmp[:2]); err != nil {
+		return nil, err
+	}
+	if d.tmp[0] != 0xff || d.tmp[1] != soiMarker {
+		return nil, FormatError("missing SOI marker")
+	}
+
+	// Process the remaining segments until the End Of Image marker.
+	for {
+		err := d.readFull(d.tmp[:2])
+		if err != nil {
+			return nil, err
+		}
+		for d.tmp[0] != 0xff {
+			// Strictly speaking, this is a format error. However, libjpeg is
+			// liberal in what it accepts. As of version 9, next_marker in
+			// jdmarker.c treats this as a warning (JWRN_EXTRANEOUS_DATA) and
+			// continues to decode the stream. Even before next_marker sees
+			// extraneous data, jpeg_fill_bit_buffer in jdhuff.c reads as many
+			// bytes as it can, possibly past the end of a scan's data. It
+			// effectively puts back any markers that it overscanned (e.g. an
+			// "\xff\xd9" EOI marker), but it does not put back non-marker data,
+			// and thus it can silently ignore a small number of extraneous
+			// non-marker bytes before next_marker has a chance to see them (and
+			// print a warning).
+			//
+			// We are therefore also liberal in what we accept. Extraneous data
+			// is silently ignored.
+			//
+			// This is similar to, but not exactly the same as, the restart
+			// mechanism within a scan (the RST[0-7] markers).
+			//
+			// Note that extraneous 0xff bytes in e.g. SOS data are escaped as
+			// "\xff\x00", and so are detected a little further down below.
+			d.tmp[0] = d.tmp[1]
+			d.tmp[1], err = d.readByte()
+			if err != nil {
+				return nil, err
+			}
+		}
+		marker := d.tmp[1]
+		if marker == 0 {
+			// Treat "\xff\x00" as extraneous data.
+			continue
+		}
+		for marker == 0xff {
+			// Section B.1.1.2 says, "Any marker may optionally be preceded by any
+			// number of fill bytes, which are bytes assigned code X'FF'".
+			marker, err = d.readByte()
+			if err != nil {
+				return nil, err
+			}
+		}
+		if marker == eoiMarker { // End Of Image.
+			break
+		}
+		if rst0Marker <= marker && marker <= rst7Marker {
+			// Figures B.2 and B.16 of the specification suggest that restart markers should
+			// only occur between Entropy Coded Segments and not after the final ECS.
+			// However, some encoders may generate incorrect JPEGs with a final restart
+			// marker. That restart marker will be seen here instead of inside the processSOS
+			// method, and is ignored as a harmless error. Restart markers have no extra data,
+			// so we check for this before we read the 16-bit length of the segment.
+			continue
+		}
+
+		// Read the 16-bit length of the segment. The value includes the 2 bytes for the
+		// length itself, so we subtract 2 to get the number of remaining bytes.
+		if err = d.readFull(d.tmp[:2]); err != nil {
+			return nil, err
+		}
+		n := int(d.tmp[0])<<8 + int(d.tmp[1]) - 2
+		if n < 0 {
+			return nil, FormatError("short segment length")
+		}
+
+		switch marker {
+		case sof0Marker, sof1Marker, sof2Marker:
+			d.baseline = marker == sof0Marker
+			d.progressive = marker == sof2Marker
+			err = d.processSOF(n)
+			if configOnly && d.jfif {
+				return nil, err
+			}
+		case dhtMarker:
+			if configOnly {
+				err = d.ignore(n)
+			} else {
+				err = d.processDHT(n)
+			}
+		case dqtMarker:
+			if configOnly {
+				err = d.ignore(n)
+			} else {
+				err = d.processDQT(n)
+			}
+		case sosMarker:
+			if configOnly {
+				return nil, nil
+			}
+			err = d.processSOS(n)
+		case driMarker:
+			if configOnly {
+				err = d.ignore(n)
+			} else {
+				err = d.processDRI(n)
+			}
+		case app0Marker:
+			err = d.processApp0Marker(n)
+		case app14Marker:
+			err = d.processApp14Marker(n)
+		default:
+			if app0Marker <= marker && marker <= app15Marker || marker == comMarker {
+				err = d.ignore(n)
+			} else if marker < 0xc0 { // See Table B.1 "Marker code assignments".
+				err = FormatError("unknown marker")
+			} else {
+				err = UnsupportedError("unknown marker")
+			}
+		}
+		if err != nil {
+			return nil, err
+		}
+	}
+
+	if d.progressive {
+		if err := d.reconstructProgressiveImage(); err != nil {
+			return nil, err
+		}
+	}
+	if d.img1 != nil {
+		return d.img1, nil
+	}
+	if d.img3 != nil {
+		if d.blackPix != nil {
+			return d.applyBlack()
+		} else if d.isRGB() {
+			return d.convertToRGB()
+		}
+		return d.img3, nil
+	}
+	return nil, FormatError("missing SOS marker")
+}
+
+// applyBlack combines d.img3 and d.blackPix into a CMYK image. The formula
+// used depends on whether the JPEG image is stored as CMYK or YCbCrK,
+// indicated by the APP14 (Adobe) metadata.
+//
+// Adobe CMYK JPEG images are inverted, where 255 means no ink instead of full
+// ink, so we apply "v = 255 - v" at various points. Note that a double
+// inversion is a no-op, so inversions might be implicit in the code below.
+func (d *decoder) applyBlack() (image.Image, error) {
+	if !d.adobeTransformValid {
+		return nil, UnsupportedError("unknown color model: 4-component JPEG doesn't have Adobe APP14 metadata")
+	}
+
+	// If the 4-component JPEG image isn't explicitly marked as "Unknown (RGB
+	// or CMYK)" as per
+	// https://www.sno.phy.queensu.ca/~phil/exiftool/TagNames/JPEG.html#Adobe
+	// we assume that it is YCbCrK. This matches libjpeg's jdapimin.c.
+	if d.adobeTransform != adobeTransformUnknown {
+		// Convert the YCbCr part of the YCbCrK to RGB, invert the RGB to get
+		// CMY, and patch in the original K. The RGB to CMY inversion cancels
+		// out the 'Adobe inversion' described in the applyBlack doc comment
+		// above, so in practice, only the fourth channel (black) is inverted.
+		bounds := d.img3.Bounds()
+		img := image.NewRGBA(bounds)
+		imageutil.DrawYCbCr(img, bounds, d.img3, bounds.Min)
+		for iBase, y := 0, bounds.Min.Y; y < bounds.Max.Y; iBase, y = iBase+img.Stride, y+1 {
+			for i, x := iBase+3, bounds.Min.X; x < bounds.Max.X; i, x = i+4, x+1 {
+				img.Pix[i] = 255 - d.blackPix[(y-bounds.Min.Y)*d.blackStride+(x-bounds.Min.X)]
+			}
+		}
+		return &image.CMYK{
+			Pix:    img.Pix,
+			Stride: img.Stride,
+			Rect:   img.Rect,
+		}, nil
+	}
+
+	// The first three channels (cyan, magenta, yellow) of the CMYK
+	// were decoded into d.img3, but each channel was decoded into a separate
+	// []byte slice, and some channels may be subsampled. We interleave the
+	// separate channels into an image.CMYK's single []byte slice containing 4
+	// contiguous bytes per pixel.
+	bounds := d.img3.Bounds()
+	img := image.NewCMYK(bounds)
+
+	translations := [4]struct {
+		src    []byte
+		stride int
+	}{
+		{d.img3.Y, d.img3.YStride},
+		{d.img3.Cb, d.img3.CStride},
+		{d.img3.Cr, d.img3.CStride},
+		{d.blackPix, d.blackStride},
+	}
+	for t, translation := range translations {
+		subsample := d.comp[t].h != d.comp[0].h || d.comp[t].v != d.comp[0].v
+		for iBase, y := 0, bounds.Min.Y; y < bounds.Max.Y; iBase, y = iBase+img.Stride, y+1 {
+			sy := y - bounds.Min.Y
+			if subsample {
+				sy /= 2
+			}
+			for i, x := iBase+t, bounds.Min.X; x < bounds.Max.X; i, x = i+4, x+1 {
+				sx := x - bounds.Min.X
+				if subsample {
+					sx /= 2
+				}
+				img.Pix[i] = 255 - translation.src[sy*translation.stride+sx]
+			}
+		}
+	}
+	return img, nil
+}
+
+func (d *decoder) isRGB() bool {
+	if d.jfif {
+		return false
+	}
+	if d.adobeTransformValid && d.adobeTransform == adobeTransformUnknown {
+		// https://www.sno.phy.queensu.ca/~phil/exiftool/TagNames/JPEG.html#Adobe
+		// says that 0 means Unknown (and in practice RGB) and 1 means YCbCr.
+		return true
+	}
+	return d.comp[0].c == 'R' && d.comp[1].c == 'G' && d.comp[2].c == 'B'
+}
+
+func (d *decoder) convertToRGB() (image.Image, error) {
+	cScale := d.comp[0].h / d.comp[1].h
+	bounds := d.img3.Bounds()
+	img := image.NewRGBA(bounds)
+	for y := bounds.Min.Y; y < bounds.Max.Y; y++ {
+		po := img.PixOffset(bounds.Min.X, y)
+		yo := d.img3.YOffset(bounds.Min.X, y)
+		co := d.img3.COffset(bounds.Min.X, y)
+		for i, iMax := 0, bounds.Max.X-bounds.Min.X; i < iMax; i++ {
+			img.Pix[po+4*i+0] = d.img3.Y[yo+i]
+			img.Pix[po+4*i+1] = d.img3.Cb[co+i/cScale]
+			img.Pix[po+4*i+2] = d.img3.Cr[co+i/cScale]
+			img.Pix[po+4*i+3] = 255
+		}
+	}
+	return img, nil
+}
+
+// Decode reads a JPEG image from r and returns it as an image.Image.
+func Decode(r io.Reader) (image.Image, error) {
+	var d decoder
+	return d.decode(r, false)
+}
+
+// DecodeConfig returns the color model and dimensions of a JPEG image without
+// decoding the entire image.
+func DecodeConfig(r io.Reader) (image.Config, error) {
+	var d decoder
+	if _, err := d.decode(r, true); err != nil {
+		return image.Config{}, err
+	}
+	switch d.nComp {
+	case 1:
+		return image.Config{
+			ColorModel: color.GrayModel,
+			Width:      d.width,
+			Height:     d.height,
+		}, nil
+	case 3:
+		cm := color.YCbCrModel
+		if d.isRGB() {
+			cm = color.RGBAModel
+		}
+		return image.Config{
+			ColorModel: cm,
+			Width:      d.width,
+			Height:     d.height,
+		}, nil
+	case 4:
+		return image.Config{
+			ColorModel: color.CMYKModel,
+			Width:      d.width,
+			Height:     d.height,
+		}, nil
+	}
+	return image.Config{}, FormatError("missing SOF marker")
+}
+
+func init() {
+	image.RegisterFormat("jpeg", "\xff\xd8", Decode, DecodeConfig)
+}
diff --git a/src/image/jpeg/reader_test.go b/src/image/jpeg/reader_test.go
new file mode 100644
index 0000000..02a2eb6
--- /dev/null
+++ b/src/image/jpeg/reader_test.go
@@ -0,0 +1,516 @@
+// Copyright 2012 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 jpeg
+
+import (
+	"bytes"
+	"encoding/base64"
+	"fmt"
+	"image"
+	"image/color"
+	"io"
+	"math/rand"
+	"os"
+	"runtime/debug"
+	"strings"
+	"testing"
+	"time"
+)
+
+// TestDecodeProgressive tests that decoding the baseline and progressive
+// versions of the same image result in exactly the same pixel data, in YCbCr
+// space for color images, and Y space for grayscale images.
+func TestDecodeProgressive(t *testing.T) {
+	testCases := []string{
+		"../testdata/video-001",
+		"../testdata/video-001.q50.410",
+		"../testdata/video-001.q50.411",
+		"../testdata/video-001.q50.420",
+		"../testdata/video-001.q50.422",
+		"../testdata/video-001.q50.440",
+		"../testdata/video-001.q50.444",
+		"../testdata/video-005.gray.q50",
+		"../testdata/video-005.gray.q50.2x2",
+		"../testdata/video-001.separate.dc.progression",
+	}
+	for _, tc := range testCases {
+		m0, err := decodeFile(tc + ".jpeg")
+		if err != nil {
+			t.Errorf("%s: %v", tc+".jpeg", err)
+			continue
+		}
+		m1, err := decodeFile(tc + ".progressive.jpeg")
+		if err != nil {
+			t.Errorf("%s: %v", tc+".progressive.jpeg", err)
+			continue
+		}
+		if m0.Bounds() != m1.Bounds() {
+			t.Errorf("%s: bounds differ: %v and %v", tc, m0.Bounds(), m1.Bounds())
+			continue
+		}
+		// All of the video-*.jpeg files are 150x103.
+		if m0.Bounds() != image.Rect(0, 0, 150, 103) {
+			t.Errorf("%s: bad bounds: %v", tc, m0.Bounds())
+			continue
+		}
+
+		switch m0 := m0.(type) {
+		case *image.YCbCr:
+			m1 := m1.(*image.YCbCr)
+			if err := check(m0.Bounds(), m0.Y, m1.Y, m0.YStride, m1.YStride); err != nil {
+				t.Errorf("%s (Y): %v", tc, err)
+				continue
+			}
+			if err := check(m0.Bounds(), m0.Cb, m1.Cb, m0.CStride, m1.CStride); err != nil {
+				t.Errorf("%s (Cb): %v", tc, err)
+				continue
+			}
+			if err := check(m0.Bounds(), m0.Cr, m1.Cr, m0.CStride, m1.CStride); err != nil {
+				t.Errorf("%s (Cr): %v", tc, err)
+				continue
+			}
+		case *image.Gray:
+			m1 := m1.(*image.Gray)
+			if err := check(m0.Bounds(), m0.Pix, m1.Pix, m0.Stride, m1.Stride); err != nil {
+				t.Errorf("%s: %v", tc, err)
+				continue
+			}
+		default:
+			t.Errorf("%s: unexpected image type %T", tc, m0)
+			continue
+		}
+	}
+}
+
+func decodeFile(filename string) (image.Image, error) {
+	f, err := os.Open(filename)
+	if err != nil {
+		return nil, err
+	}
+	defer f.Close()
+	return Decode(f)
+}
+
+type eofReader struct {
+	data     []byte // deliver from Read without EOF
+	dataEOF  []byte // then deliver from Read with EOF on last chunk
+	lenAtEOF int
+}
+
+func (r *eofReader) Read(b []byte) (n int, err error) {
+	if len(r.data) > 0 {
+		n = copy(b, r.data)
+		r.data = r.data[n:]
+	} else {
+		n = copy(b, r.dataEOF)
+		r.dataEOF = r.dataEOF[n:]
+		if len(r.dataEOF) == 0 {
+			err = io.EOF
+			if r.lenAtEOF == -1 {
+				r.lenAtEOF = n
+			}
+		}
+	}
+	return
+}
+
+func TestDecodeEOF(t *testing.T) {
+	// Check that if reader returns final data and EOF at same time, jpeg handles it.
+	data, err := os.ReadFile("../testdata/video-001.jpeg")
+	if err != nil {
+		t.Fatal(err)
+	}
+
+	n := len(data)
+	for i := 0; i < n; {
+		r := &eofReader{data[:n-i], data[n-i:], -1}
+		_, err := Decode(r)
+		if err != nil {
+			t.Errorf("Decode with Read() = %d, EOF: %v", r.lenAtEOF, err)
+		}
+		if i == 0 {
+			i = 1
+		} else {
+			i *= 2
+		}
+	}
+}
+
+// check checks that the two pix data are equal, within the given bounds.
+func check(bounds image.Rectangle, pix0, pix1 []byte, stride0, stride1 int) error {
+	if stride0 <= 0 || stride0%8 != 0 {
+		return fmt.Errorf("bad stride %d", stride0)
+	}
+	if stride1 <= 0 || stride1%8 != 0 {
+		return fmt.Errorf("bad stride %d", stride1)
+	}
+	// Compare the two pix data, one 8x8 block at a time.
+	for y := 0; y < len(pix0)/stride0 && y < len(pix1)/stride1; y += 8 {
+		for x := 0; x < stride0 && x < stride1; x += 8 {
+			if x >= bounds.Max.X || y >= bounds.Max.Y {
+				// We don't care if the two pix data differ if the 8x8 block is
+				// entirely outside of the image's bounds. For example, this can
+				// occur with a 4:2:0 chroma subsampling and a 1x1 image. Baseline
+				// decoding works on the one 16x16 MCU as a whole; progressive
+				// decoding's first pass works on that 16x16 MCU as a whole but
+				// refinement passes only process one 8x8 block within the MCU.
+				continue
+			}
+
+			for j := 0; j < 8; j++ {
+				for i := 0; i < 8; i++ {
+					index0 := (y+j)*stride0 + (x + i)
+					index1 := (y+j)*stride1 + (x + i)
+					if pix0[index0] != pix1[index1] {
+						return fmt.Errorf("blocks at (%d, %d) differ:\n%sand\n%s", x, y,
+							pixString(pix0, stride0, x, y),
+							pixString(pix1, stride1, x, y),
+						)
+					}
+				}
+			}
+		}
+	}
+	return nil
+}
+
+func pixString(pix []byte, stride, x, y int) string {
+	s := &strings.Builder{}
+	for j := 0; j < 8; j++ {
+		fmt.Fprintf(s, "\t")
+		for i := 0; i < 8; i++ {
+			fmt.Fprintf(s, "%02x ", pix[(y+j)*stride+(x+i)])
+		}
+		fmt.Fprintf(s, "\n")
+	}
+	return s.String()
+}
+
+func TestTruncatedSOSDataDoesntPanic(t *testing.T) {
+	b, err := os.ReadFile("../testdata/video-005.gray.q50.jpeg")
+	if err != nil {
+		t.Fatal(err)
+	}
+	sosMarker := []byte{0xff, 0xda}
+	i := bytes.Index(b, sosMarker)
+	if i < 0 {
+		t.Fatal("SOS marker not found")
+	}
+	i += len(sosMarker)
+	j := i + 10
+	if j > len(b) {
+		j = len(b)
+	}
+	for ; i < j; i++ {
+		Decode(bytes.NewReader(b[:i]))
+	}
+}
+
+func TestLargeImageWithShortData(t *testing.T) {
+	// This input is an invalid JPEG image, based on the fuzzer-generated image
+	// in issue 10413. It is only 504 bytes, and shouldn't take long for Decode
+	// to return an error. The Start Of Frame marker gives the image dimensions
+	// as 8192 wide and 8192 high, so even if an unreadByteStuffedByte bug
+	// doesn't technically lead to an infinite loop, such a bug can still cause
+	// an unreasonably long loop for such a short input.
+	const input = "" +
+		"\xff\xd8\xff\xe0\x00\x10\x4a\x46\x49\x46\x00\x01\x01\x00\x00\x01" +
+		"\x00\x01\x00\x00\xff\xdb\x00\x43\x00\x10\x0b\x0c\x0e\x0c\x0a\x10" +
+		"\x0e\x89\x0e\x12\x11\x10\x13\x18\xff\xd8\xff\xe0\x00\x10\x4a\x46" +
+		"\x49\x46\x00\x01\x01\x00\x00\x01\x00\x01\x00\x00\xff\xdb\x00\x43" +
+		"\x00\x10\x0b\x0c\x0e\x0c\x0a\x10\x0e\x0d\x0e\x12\x11\x10\x13\x18" +
+		"\x28\x1a\x18\x16\x16\x18\x31\x23\x25\x1d\x28\x3a\x33\x3d\x3c\x39" +
+		"\x33\x38\x37\x40\x48\x5c\x4e\x40\x44\x57\x45\x37\x38\x50\x6d\x51" +
+		"\x57\x5f\x62\x67\x68\x67\x3e\x4d\x71\x79\x70\x64\x78\x5c\x65\x67" +
+		"\x63\xff\xc0\x00\x0b\x08\x20\x00\x20\x00\x01\x01\x11\x00\xff\xc4" +
+		"\x00\x1f\x00\x00\x01\x05\x01\x01\x01\x01\x01\x01\x00\x00\x00\x00" +
+		"\x00\x00\x00\x00\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0a\x0b\xff" +
+		"\xc4\x00\xb5\x10\x00\x02\x01\x03\x03\x02\x04\x03\x05\x05\x04\x04" +
+		"\x00\x00\x01\x7d\x01\x02\x03\x00\x04\x11\x05\x12\x21\x31\x01\x06" +
+		"\x13\x51\x61\x07\x22\x71\x14\x32\x81\x91\xa1\x08\x23\xd8\xff\xdd" +
+		"\x42\xb1\xc1\x15\x52\xd1\xf0\x24\x33\x62\x72\x82\x09\x0a\x16\x17" +
+		"\x18\x19\x1a\x25\x26\x27\x28\x29\x2a\x34\x35\x36\x37\x38\x39\x3a" +
+		"\x43\x44\x45\x46\x47\x48\x49\x4a\x53\x54\x55\x56\x57\x58\x59\x5a" +
+		"\x00\x63\x64\x65\x66\x67\x68\x69\x6a\x73\x74\x75\x76\x77\x78\x79" +
+		"\x7a\x83\x84\x85\x86\x87\x88\x89\x8a\x92\x93\x94\x95\x96\x97\x98" +
+		"\x99\x9a\xa2\xa3\xa4\xa5\xa6\xa7\xa8\xa9\xaa\xb2\xb3\xb4\xb5\xb6" +
+		"\xb7\xb8\xb9\xba\xc2\xc3\xc4\xc5\xc6\xc7\xff\xd8\xff\xe0\x00\x10" +
+		"\x4a\x46\x49\x46\x00\x01\x01\x00\x00\x01\x00\x01\x00\x00\xff\xdb" +
+		"\x00\x43\x00\x10\x0b\x0c\x0e\x0c\x0a\x10\x0e\x0d\x0e\x12\x11\x10" +
+		"\x13\x18\x28\x1a\x18\x16\x16\x18\x31\x23\x25\x1d\xc8\xc9\xca\xd2" +
+		"\xd3\xd4\xd5\xd6\xd7\xd8\xd9\xda\xe1\xe2\xe3\xe4\xe5\xe6\xe7\xe8" +
+		"\xe9\xea\xf1\xf2\xf3\xf4\xf5\xf6\xf7\xf8\xf9\xfa\xff\xda\x00\x08" +
+		"\x01\x01\x00\x00\x3f\x00\xb9\xeb\x50\xb0\xdb\xc8\xa8\xe4\x63\x80" +
+		"\xdd\x31\xd6\x9d\xbb\xf2\xc5\x42\x1f\x6c\x6f\xf4\x34\xdd\x3c\xfc" +
+		"\xac\xe7\x3d\x80\xa9\xcc\x87\x34\xb3\x37\xfa\x2b\x9f\x6a\xad\x63" +
+		"\x20\x36\x9f\x78\x64\x75\xe6\xab\x7d\xb2\xde\x29\x70\xd3\x20\x27" +
+		"\xde\xaf\xa4\xf0\xca\x9f\x24\xa8\xdf\x46\xa8\x24\x84\x96\xe3\x77" +
+		"\xf9\x2e\xe0\x0a\x62\x7f\xdf\xd9"
+
+	timer := time.AfterFunc(30*time.Second, func() {
+		debug.SetTraceback("all")
+		panic("TestLargeImageWithShortData stuck in Decode")
+	})
+	defer timer.Stop()
+
+	_, err := Decode(strings.NewReader(input))
+	if err == nil {
+		t.Fatalf("got nil error, want non-nil")
+	}
+}
+
+func TestPaddedRSTMarker(t *testing.T) {
+	// This test image comes from golang.org/issue/28717
+	const base64EncodedImage = `
+/9j/4AAhQVZJMQABAQEAeAB4AAAAAAAAAAAAAAAAAAAAAAAAAP/bAEMABAIDAwMCBAMDAwQEBAQGCgYG
+BQUGDAgJBwoODA8PDgwODxASFxMQERURDQ4UGhQVFxgZGhkPExweHBkeFxkZGP/bAEMBBAQEBgUGCwYG
+CxgQDhAYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGBgYGP/EAaIA
+AAEFAQEBAQEBAAAAAAAAAAABAgMEBQYHCAkKCxAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1Fh
+ByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNk
+ZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT
+1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6AQADAQEBAQEBAQEBAAAAAAAAAQIDBAUGBwgJCgsRAAIB
+AgQEAwQHBQQEAAECdwABAgMRBAUhMQYSQVEHYXETIjKBCBRCkaGxwQkjM1LwFWJy0QoWJDThJfEXGBka
+JicoKSo1Njc4OTpDREVGR0hJSlNUVVZXWFlaY2RlZmdoaWpzdHV2d3h5eoKDhIWGh4iJipKTlJWWl5iZ
+mqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uLj5OXm5+jp6vLz9PX29/j5+v/dAAQA
+Cv/gAAQAAP/AABEIALABQAMBIQACEQEDEQH/2gAMAwEAAhEDEQA/APnCFTk5BPPGKliAB718W7H2j3Ip
+VUuwJxzTfKXacde9VBhYRUBAyO3pTmUAbSMU5WGmybywzHGAMdelPVFC+n1qXZCuyaJADxjj2qzbBMAP
+xz1rKaVib6ltLcFvlIx2pLy0dwuAMMBnH1rFON9RNsszAZPFEYHldPzrOy3KewmBk9qUABugxjtTVmiW
+xWRcjp+VJtHXgVL3K6AgBDdM9eRTNzAZViOe1VyxaJavuf/Q8aW4mUcSGpo764AyHz+FfnnJBvVH1UsN
+CS1Q/wDte4Trip49ecA7g3FSqMW9zlqZandxJ4/EKADcSPqKni8QQMT865qOSUNjiqZdNbFiHWYXz84N
+WE1KNsfMKj2zirHHPDSj0JFvo2H36d9pUjg1sqykYOm0KbgY60omXPXmr9pFkco3zBnrQzjGcnrRzp9S
+bEbuOvao3fisZSXUpIYWGKGcbetTCSswsxnmACkYrtNSpJ2YNM//0fnK1BD7sDg9KmUHeOe/Svid3qfb
+SdmQ3AHmnr1pGBC5z19a0hohNiJkensM1J0yCKmY0yZR82e+BT1BxnpmpepN9SRCR0NSpweOoPWs6isr
+ijuWIZGBA/lVwzMVFY8ibuhXEfr+tOz8hIqUymhRnJGTSc5wBVRRDFPXBHJpB3qdmV0EX7vXmoyfl685
+p2dxWR//0vFsHZ9TQv3T618Bqz7PSwwn1phPXpSWrQEUhIx0NVXc7j0rSNwViCS4dWYpJj3BpBqVzGy7
+ZmHSq9kpblSpxkveQ+PX7uMf6wEDtU0fi24TAkX8jTeCjJaaHDUwFN7aFq28aL/GCMGrtr4xtHGGkA+t
+YTy+a+E82eAa2LsXiWzI5mXPHercOsW8hwJB+dcUqVSCOKVFxdmiwl7E2MOPzp4nQ9GH51jzNbmUoOIC
+TI4okOaUXoybDCevNBPHX8qIO4mf/9P52i4dix5zjp/n1qZFBCmviL6an2kt9CGcYnJznJpOBwegq4vQ
+L9xIUytSkfu/bv70p7j6EnQgjHSpFGVqXclkkaHb1+makUHgdazm7IFuSKOasrnjis2+oDm9qnIHlgd/
+es7gxqjkt1NLwH4xTTEhjkhutM3D15oGhkcnBGRTDKu3A7H1rS3cLn//1PEhJ8uM557UvmDaa/P7a3Ps
+xpZcZ6mo5WG45pdUC2K8ko5JIzWfcTqu7HPHrW9OLbKWhSluVLNz3wKrS3I3KfcV1Rg9CrpXK8l0F7io
+pLnLnJHGOldMYJGMpu5XNwuxjyRTBcAAjd1HWtfZmPORy3WAWWQDOM4PWtHRru6DFlmY88ZqKsVyXaKp
+QjOoa7axe28G/cWqhoHjO/n1WeJwSkS9c981wUcFCopPayFj8JC8VFbs6e38VldvmHHFaMHimAoCzDB6
+V5U8FKz5TlrZU4/CXYtdtnXIarEepW7jAcfnXGqEoKx5tfBzh0P/1fnqEAsc/wB6pI9owAD1618Qn3Ps
+35EE4UzHrx79aXaMcdaqAMWIADvj271IMeXg59KUmNLQkUDfjb1FSLxzg0pWJRLGAQAeMVIoA+uaxlaw
+0SoF/u1KowwwDUcwuo9wMjrUrY2ZPOKy0KY1T1NMdwG/CtBEFzMqnIPNUZ75FBJP5mtIQvYfoU21JFVs
+N271AurRE/e611xw73Yj/9b50GsQhOXHWnpq8JX7w4PWvjPq76H2fzHjVYCud9Q3GrRAZDUvq75kNbMz
+7vV0zjdjNZ82pqzMcj7tdlPDtIiVWKKct+AxwRxUbXi7VPJAIrZUdEZOsrsga8DFgelQtd98g5P6V0Qp
+GE6qIUut2cZ470kd2FjYc4Oce1bSpJ3Rzxq21GNcDZhSeg710ujKRbKzAg5rkxceWnqd+XtOo7bD9cl8
+qxLDPHasXwUvmyXU7Lgl8cegrnw2lGbZ14l3rU0bl3gMQCRgVU1y7WytUZQzMRwBXPRhzWRvVny3ZW8N
+6xPdXBikiZc5IOa6GG6nDsd5xnAyfaliqEacrGOHarx5pI//1/nuL754HWngEkYx1r4VWsfaMjk4mP8A
+OgnjPH1rRMLKwR4A2jH1FPA+TNRIa0ROvQcY4p4GF/pUskmi6+gqRACvPrWMnpca3JABjFSKCQOnFS2u
+o7E3XBOKcR8ucdKzUkDGSHGemKpXchVuP0rSmDMfUrl1J5rn9TvnVCc9OtelhoJtDekW0Yb6pId3zdRw
+RVT+0pAPvc57CvbhQVrHlTxD3P8A/9D4tbUpTH1I7cU/7ZdMnyqcdsV5vsErXPbWJbHLdXzYQDBY8c02
+6udQjyGVuD1FHsqfMridepZ6ED3s4IDqeD0I68VEt7J5hy3GO9aKkuhPt2BumeRjnv3pJLlgwBYE8ZqH
+T2GqujYLcuWYbhj0zTHm5B/vcGtVBLYzdRtEcUueoGB3FOjmBjcBhx2NNx1IhO+uwtqd93EgA5YcV32n
+IqwrnAz2rzMx+FI9nKldy+RmeMpfLs8DGTxTvBUKw6Csjry2WPbrXHB2wzfdnbUu8SvJF4xh1LDAJNU9
+UtVmDs3IiGB6CuenNx1R0yjfRGd4aRTqJdFG1ARXRgANg4/yK0xbvJehnhlaL9T/AP/R+e4x8xx609F+
+YZ718L6n2ju2RzqTKcYpQMjsc1pHTQWlgjUjGVH0qbkr0BqJSKRMi+uBx3p8a5HYVD8yb32JY15FSKpx
+nisp6RuNbj1BzUyrnkmo6FEqrz7U/advHOazvcRHIuSazNXDpbSSJjeqErnpmrp6CueeXusahO5zKi8c
+7VrPuGklUiSQtkd6+po0YQs0edVrTd1cqeSoJOB0xUCxpnouAecV33e558rbH//S+KdmFHTk1btywUqc
+YNcEnfc9SGl7DyAJ1AHIParx+YZ4/Guea2OmC1dhGjQn5kXNQtbQFiWiTlfSoTa2G0nuU5bG2aQ7V2jP
+JU+1RSaXC2GjuApyOorX2slYz9lF3sV/7MmViFljaoJLG6SQbkyDXQqiZzyg1Yg8i4jBJjIBpsaPyXXB
+Psea1TTMJJqysaeh2u/UUfP3QCBXdQJtTpivFzKV3FH0WURtCT8zmPHcrhkRSOWro7O28rRYIgOwGB3r
+mnph4+bOxNvEy8kWFi+ULxwKzNRkMemSPj/WMT+FckNWdLv0KPhCMmGSZl6k1ulC3zY5x2+la4r+IyKH
+wH//0/nuIHB9c9KevUAHk18La60PtHuRy/64+lOGBniqXcOlhUIxwB+NSrynSpndFImQc4A7d6lQccdR
+WcyUyWMccDnPSpVAwM1nLYaeo5BjrUyjFTugJIwd2Kfgkc59qyTs7Axkigqao3qBkYdiMVpTeugHmF7b
+hbhl2/dJB/OofJBTHfp9K+ppTbimeZOK5issYG5W7VWdBnHXB65rrjJs5JLof//U+LYtu7leM+lSpxIR
+7VwO90eskrNkqZLo3PXHoausxI4wa557o2p6JitnOCoqvI3zEkdF5qIrUuW2ogO1iWHeibazIQncHA+l
+DT0aCyaaZGNm8kA9PSl2qy9SB78Veq1ZCs9BkOGUrj86zdQGbllVMAe1bQdpGE1eBo+FoCbgtxkY966+
+E4hOeo5rycxleR72VwcaRx+t/wCmeJ7WAdDIOPpzXbSpt/dkcRr+tZ4j3aVNLzNKOteo/QjuiY7Jm6Ej
+ANYnitvL05YRxwOf8/WuXDK8l6nVUlaLZb8NQeXpijgZB/M1oIpyPzx74pV2nUbHTVoJH//V+fYhgnnv
+Txy4GBXwse59m7kMyEzkj8qkQfKatdgewIo7nIqdQAnXms52RSehMoHHPapY1wMAgVnKzFtqSwjjg4qR
+VJHXIzWc2rDiPVeeD+FSDqKh2sBKo54p+Pl61jbQG9RrDIPNU7teT6VpCztYDzfxGskWtXESdN5Pp15q
+gN2GZpB0r6ig17OL8jz535miCSPdnaxHHpVV48D7xIB7iu2LOOS7H//W+MCoeIDcc5p4VhIMDkDvXnpq
++p6zu0SwZZlVm6HJFWyRg89MdawmkmrG1NtpiMcY5OevFQ7AWOT0FSkkU9UKUPmEh8jt+VMdGLDLYAIz
+xUtrQfLo7Mj2SHjePrSspxgEk1rdGSTGJjymLEZArOjAd5GLHk9DW0NGznqa8qOj8IRHBbrnnmugu08u
+1ZiSMn868LGz/eH1GAVqKOW8LR/bfG4c8rCCx46HpXZspk88jHzMf04pY7eEfIjDO7nLz/yKmqh/sjwR
+EFwAemcVhamkmpTRxKpyCN2RWeFsveZ01FpbubsEaRWyqhAxnH5YpxIx8rf/AFuK5W3Jts2Wisf/1/n5
+SSxHOM+lP7jGa+EVz7R2IpATN1IIpwB55NaJ2FuhYzx3PvU69OQaio7sEiZOvfpU0YwmMVnJ26DRJH2G
+DipUyR361jN6FIeq8/0qUdBxWbkCRIg/D6U8j5e9Rza3ExrA8nmqt0Dkmri9BnnfjlSmvuwGQyhulYkr
+yL86DANfTYRp0o3PPr3UnYbBOWU4zz7VHIGIJVjkGu1x5Tl5ro//0Pi3fhgMHJPXFWCeQwLe9ec+jPXj
+1JIM7gw44qy+WPUjkcVjPdGsNmgdsNkjJFQ7mMhAB5FRHuXJ6WRIw+VwCc9KbtPy5JyCKgdmxhBDNj8s
+U1Cyr0J/rWultDOzTuMuSiWjlT97r7HFZ1nkk4bIPXiuiD3uc00rqzOy8Lw+XBuJPQcGrXiGYJYMwJHB
+xXz1d89U+sw8eWmkuxi/DGP/AEm+vHycYUE/jXXu6w2vzdcfiaMw1qpLsjnwi/dt+bKCn5nlw2W5Gacw
+GD2wB+dcq2O4AnyADoM80QLukUsp9f0qb6XHuf/R+fkOWNSfwjnivhUz7Nq2hFJ/retPA4PWrWuwul2L
+FjAA6VMMFeTms5PUpbEw6/hxUyfd4as5PUETRds9KlA+Xk96ym9FcaHJgGpOv1rPUpIkXg8mnNnGaz5r
+aCaEPeqtx1OT0rSL2sC1OG+ISquoQuT1UjP4/wD16wGEZUYGPevo8E26UThrpc7G7ICDzg+1ROmF+91O
+K7VKVrHNyxWx/9L4uuVAcH371JvKqScFPU1597pHqtWbZNZnc+QxI9atv8p4z9ayqPVI1pr3WyPBLDGf
+qKYnExyeQKlFPQXH7zgdetSk8rgEnis29i0lqxijLEjjt1poU7iVHHpVX3uRbsZl+2IvLX+I56U/TUUA
+KxGSfSulu0XY5oq80drpcZSzHvjpWd47fy9O2g8kjpXz0ZJ1kvM+skrUnbsSfDm1C+HlfJ/euXIPRq3b
+lleQRYBCg5HrSxk+au/IxwkbUokRw0u0cBFyR70wEEbm6sc/gK5YuyZ1PzFVgVG4ZIzmpbTaJMt07+3F
+Q9i7n//T+foic55yTipRkYBBxXwaPtXuRyZEg4pWII4qk7C6BFwAf51MhG31+lZ1Frca2LKHn8qlDY6L
+UNgl1Jo+2akQ9BWVR9xpDgffrUq+wrO7tdDsSIPUUpPHvUK1xMM8HA61WmwWOB+NXENjiPiMhE1sw9WH
+8q5vqRnjivosD/BXzOKv8QkKgZBA6ZpV27MkDOa7ObsYI//U+MdVGxlK9zninqd1sCQM45rzYaxR68vj
+YtgT5h6jvV6Q5X0+lZ1n7yLofCxhOenfFMTI3cdRzWV9DWw5ARI3qPSnMSqKCOSRUy6FRurjFLjPp9KA
+xx06n1q1qjO70Me6YtcOOcKcH1q9oqF5l75Oea6KtowOagnKol5neWcSJaocdgRzXGfEm53ERKfvEV89
+gvfxCPqcXLlw8/Q6fwkph0aCEg4VB/KrsDbmcnA4PWoxFnVkxUVaml5IgR9sMj4+ZzSTuTjcOB0/CsUt
+zo0VrhCQiF2GcAn/AD+dWLRlYZPQ8cVEk7aF+p//1fAIuvfOakxnr+NfBJ2SsfaN6jJRiUA9RSheCMfn
+VXEtgjUk/wBTVgfdBwOfSs5stbE6g7unYVLGpwAazYvUmjHHanqDx061lLazKuh6DHBFSID27VEthkin
+5cUHOPxqLvqJiEYziq8/FaQ8hHH/ABEVvIhYYyHNcsrSZG5RyOtfQYC3slfzOPEX5tAA+amHIjO31ruu
+rHNa7P8A/9b411QMIwSDnNR2xYQNkjnnkV5sLctj15JqRLZjBzweeSKuycHkD8qyq6tF0rqLI2OTnK5p
+sGWbBHQd6zWxo3qSdXLYxTpPvLnvjrWdr2LvuNYYUnj6Uxyu4/KMrVx6ky6GOGLSOwXIYmr9n58UQeFg
+svbcCRXTVty2ZzYZOVRcpvDW721tv9LsBIpAHmQNn9K4zxPqSX2sK6hljDDO6uHBYWKre0i9PxPTzDFS
+VDkmtX9x2mm65YG0REnTccAc1rx3EJgbbIpyvUGuDE0JxleSO+hXhUj7rGK/7uNcj1P6UmSU+fHPJrlS
+5bnXe9mA/wBSQDxzkfh/+qp7YbIipbaOufwrN7WLR//X8DTO7I9e9SJnAIHNfBPY+ze4yQEPnGacoHof
+amvIELCDnnqanXoBzUSaRSRMo5J744qRBxzUSelwRPEMingdsVk7pFLckA5p44OTUXuh2HqOKRskE4qV
+sKQADkYqvcj8aqD0shWOV8fqDYqc4w4rj5D365Fe7lz/AHdjlxC11CNhyRnp0pqkYOc9a9G/c5bLof/Q
++OdUH7sDnOeKrREgEN6V5tN+6exUT5rlm1IC9ec9qssBg5xzjNYzdpJmkNtRhJzgflSRDqD6flU7FkqZ
+55+holblfUVk90aLYYxJH3agvGCQPkjJ6DFaw7GMtrmdYpkcg1q6fsEwB6Ct8Q73M8DpKLZruu+3ZM43
+ADNctfeHt120cEjO3U7scmuDCVvZNux7GLw31iNmypLoN5byByrYzzxVfzdTtXcxzSAD3r0oYinW0Z41
+TB1cPrBlyy8S6nF9+QOMba1LbxepwsyFcYztrnrZdCabgb0MzlCyqGrZ+IbK4U7bgDg8E4rXtL+CQBfM
+VgTn9BXi1sJOno0e5SxVOqrxZ//R8DiGScip0HOAOa+CfkfaDJV+fpRjBIxT6iHRjnIHFTAYXBIqZ7lK
+1rEygZxUsYI68Gsm7gkSQ59PwqZRxx681lK/LYaHpjI7U8ge1ZK3Qdxw5XpijAOelGoNCDPOOaguOf8A
+CnFCOa8cJu0uTI6EHH41xLou3gYBr3culaLXmc2JV7DIx8pXHFNAxknGc16b0OM//9L44v8A/V5wMiqy
+/dbI6j1rzKb0PYq7li2C4HHHrVmT5Rx9ayqN8yLppKIwEj6+lPhAGfzNQ9rmqexKvcAYxTTjeM88isty
+xhI3mqWrNwqADB5ranozCo9GQaeBjByR1q2GIzjj2rasryMaErRVieK/kjwHbIGOBVrSbtZNRmlfAB+U
+ZP0/wrinQsnKPY9eljVJqM9zYgMTMAjoRnoRUV3ptrMrl4UY+o4rzIzcHc9Nx5kZc/hu1lUlBszWPqHh
+Z1lYxHJHQZr0cPmElpI87E5dCptoZl1o9zADuVhwfaqqz39tIPLeVMDoPpXqUqsKyPHr0KuGfun/0/Bo
+Rzlh+lToOBXwGltD7RjJQcjtSqODn9RVLyEEY96nVflzis6l3YpaEigdscCpl5XqeKiaVgJIhn/61Sgc
+AdOaylqkyluPXv7U8gZx/Koih9R+ML1ppHynIoRLEAqKYZ6njFVHRXBnP+Lk3aXN0yFz+tcM6j15r2cv
+ejOfEdLjMDb94ZHaoyvy9O/Y16mqOQ//1PjfUcBcnjJ9KqggRldo579682krxPWqO0i1AM7c59OlWJMA
+/lWVTdGtJaNiLg89afCMkjOOO1ZtGu5IpCs+emf6Uxug6Y45rO3UruiFmKydvrVK8ZZLogAn+ldNOPU5
+pvox1pH+7JzhumKlVfkOckmnNq5MI2SGScE5HTBqWDcBkMvJ6CpuraFWd+w3ULqW3gMsb4fOOO1VNP8A
+FGoxsyy7JFA6kYNKGDp1Y+8aSzCrh5e7qjesPE9s8eJh5eRnIrShvbG7AMUqZ+teVVwU6TbWqPaw+OpV
+0lezJWgV4yAVYc8GqV5pFpMMyW+O+V+lc1KtKD906q1KMlaR/9XwqPg9e9SqOQQetfnyfc+0Y2UfOKFy
+B25q1K+grCx564FTL93rilMaRKvXtz7U9OnfrWbtbQdyaLoOtSr/AFrOSstCl2HLnPNP5z1qNEA9s4FM
+7daSVhMAcZzTJlOD1prYDC8VoTpdzjP+rbH5V5w/nDkj9a9rK2mnc5sVeysQxmYMd4wDUkLM3LZFeu0u
+hxan/9b44v8AG35cnnNVArKM+tefS+HU9WoveuizCQCBuOT61ZcnGcgmsZrVG1N3TGrk4x3qWBQGPUHF
+RJspa7kjAFmwOp6Y9qSXPyDtkVgnsa9yFxhWYrj2rJUtJcFgB75rqpSvc5K0bWSRetV2Rg98c8d6UZCs
+M8AZPNTLWTuVFWiiIgjcBnDY5NToCsZAyDnqKUrtJDitWUfETbYI1Ixz7ZNYkXyzMAMevHSuvDfCceK3
+NKGIG0Ppj73pVdi0F2MMQDj7pog7uSYTXKotGxYajeQwA/aHOAeG5zV6HxJLCdk8KuvTK8dvSvNq4OFV
+6aM9ajmM6SSlqj//1/DIOpHGaljG4cYH4V+e62PtWEi45pvzcjIwKuLsS9RYwM9s1Jj5Md88molqUnYl
+Xlue1SryKlpodiWLGB6GpBWU9FcEPQd6eBgjNRYpCn7opDyM0Ru9RMQdxxzSODVR12F6mTryFrWVcfeU
+j9K82ctvxxXrZY7XRjiF7qEjUHgkYFQuuCeB19K9dO5xSSsf/9D44uwSmM556VXkXMJHGenFeZTvoexN
+Ilg4bBOTjrip5uQcY4x2qJ3ckVDSLQ6MdCPpT4uGIHGOaze1jRJbkjFgDgYpjEgqxI7dqyWhpuQ3rkQS
+HjpWfaR9/wCf1rqp6RZy1Peki8MGMfNx7d6jc7Yyc4ye1JptjvorjBneDw2QMVPEHBIOFwR3qbXF5ox/
+EDl71U4ODiqDrjJABIrsoaRSRx19ZNmrbD/RgScAjsag1CF9pcBTxWMXaTNZRbirDrF3CbSQTg9f5U+Z
+QJxwFBA/lRZc1yrvlsf/0fD4fvE1MnYk/nX5/wBD7R7iP1yMdaQLzRF9LCa7hGvOf0qXHy9qUtxkgXnt
+UqDk0pK7sCfYkUcYA6VKnIFYzWha7jhT1Hes42TsC1FIyOlNYEU4oT8hB3pGzt68U1tYVzP1IZQ56Yrz
+e6UC5lUrwGOB+Nepl3xNGddXiV0Y8jGPw6U1slSCec9a9m6OCx//0vji/wAiMY45qux3Rbh/KvNprRHr
+1N2iW3z5gHY81ZnA3AnipqfErF01aDuPhwXAPHSpME7jyOK53pKxt9kT5vmB4IPHFMkwcEnkkVPkPpqU
+75zt8sDBJ+tNtVCt8x59q6rWi7HKneSuTxAhCAQOcc0yQEA4br2FRfVqxfK7LUYw4weCu3n1qbdiNmOM
+EgY9KV9CXo2YF+VMyNx8zE5qMgsrhccdq7obK5xT3aRdtTi2VmPGMdc+lWJgr5XAGMfrWElq7G8XsmUx
+C8ZJjOQQcgdakVXZxv4PbjvijmTJlBr7z//T8RjIGPc1KhXqPXivzzdH2bWoSD5sg4INJwRzz9aEw2BM
+dM5qTjb60SdykSA88enepEOOKU31EkTREdqkH54NZT1Ra0Hd/SnHip0Q/IVjwPWmsePehPqiWC9T1pHH
+B5oTaFYoaiOOn415tqm4X842jiRh+tenlvxMzrfCV05JLDHTgD3ppYFCB6/Svb1aOHbc/wD/1PjnUv8A
+Vct3qq53IMeleZDZM9ipu0SRYDAkkewqxI24+w6VM90EL2aJYCN5AbgY/lU0R+9z+FZS0N076DZSNzc5
+yc0xyoUfh25qG9mU7XZnXLZuTgnaKsW4zjqRnriuiStHQ5Y/GOxgEVFuwpBA9+MVKsytUOjILcZ7daNQ
+ZEt3O7k8fjU3ldDbVmc/cPm4QYA/rS8gMD2zz6139EcF7tli2lVIAhVsYPap45lZiFG08DnrWLT1Ztde
+6kidMGMZHGD3pkqjcOcZ447cVjF6o0nZxbP/1fElB7HOD1qVQOOtfnSVj7VsRwNx60Ywp4PFOCVwvoIg
+HB9e1SAAjGcUpa7AiQDkf0qRRxgZGaUrWuwJYwOByalXGB6A1nOzRS3FAI70/GSOtTtuMVunemHGKUe1
+iRRgk0hwKqPkDZS1DkYFed+JEEWs3Q5GHz9M816GX6VPkZ1V7pQibOck8jrUbgDIA4r3epwJ6an/1vjf
+U9vl43YyelUwR5RGTxx0rzaauj1qrSlckiA3D5j071aYjccnHSlJXaHDZsntzkY3Z96lQjnnt9K55I3T
+GsRyDnI9ajfKKOemDxio7Iq+7M3JeVieSecYq7AMKF3Fj6V0zVlaxzQfvXGDpy2McdajjydwIO2h2Hd6
+ND41G8nAGMdD1qLV2zAFORl+R1qEtVcTdkzBl5vSCMgNwKlVG3MFJ4967r6I4bO7JY1BjJJ7Y6dDVqyQ
+CYkqTjvWEn7rN1unYtYQLntg9+tV2+aQBlbPTHrxWNPzNqux/9fxWMDP41IoUYr87WzPs3uEgBJ9uKAB
+tINCSbDoJGAMcdKkAG04/KiSXQpPQkXGcY7VINp56VLBbksYBxUgAx0rOadrDT7C45pxA4wKVkUhX7ZF
+NODzRG1iGAA5NIwGeB1FNLuBTvEBI7V5/wCNFP8Awkl1tIABUdP9kV35cl7VehFbWBlxbNuD8x9aY20Y
+xjOa96xwWstD/9D451XAhBOOTxVRlURdOa8ymrI9epu2SW23cMgDHb1qw4BcrgY9KJ25kOn8LLFsArYG
+McVIu3GdvOOTWLV2axeliOYgM+FzmoL1isAwME4GKVloDvqVoI9zEEAMMGrbIAmFUdcZrSbS0MqcWQBg
+0WMg4PalhCFCdpHWiXUFZtaDtqhywBwMY9Ko6sR5oXPTn2p09WianVIyEYebwMEnPNTL5eXBA4Hautp6
+WONSWt0TQMhUgjKj1q5bhQWKMDwBxWE20mbws3ElXy9pUEDqMd6hfaHUgE7B2HtWUO7NKiVrI//R8Vj4
+br37VMuT9TX55HY+ze+gMM5PSkzwR1AoWu4IRRg9akPC8HpSlfRFLXYeCN3TkipV60ulmBJGeR71IvTP
+Ws5K6Ghwp3FQkMH6elNJ461didOoq9SaRjxSjqwuVbgZkH1rzrxWwk8QXj9hIR+XFejl1nVt5EV/gM2I
+qCTjPTH50x+eRxya961tzzt9j//S+O9WJ+zjgDmqJYGDJYZNedSSaPXqu0tew+L7wPTvmrBO5x7Y5qZf
+EnYUPhZdhUZyO3pTgcOxGOR0rnbu7M6fhK9yVEvPy59O9QX75VBjnIoS+ETa95DLbdu6nnp7VYGSMk59
+vxrWemxlTZCvMe4BTg/nSRKHVsjGfpSe90C1VmTsG2H2x3rF1J1eVxgcHr+FVQWtyK7srFGAgOrHqDnp
+7VPGNzO/O3kkCut7pnGno0idCCNyt6Crca5ZiDjpgVzXsnc6X71rDyBgcbtucn8KhjBMpOfrz7VlTdtT
+Sqm7JH//2Q==
+`
+
+	data, err := base64.StdEncoding.DecodeString(base64EncodedImage)
+	if err != nil {
+		t.Fatalf("base64 DecodeString: %v", err)
+	}
+	if _, err = Decode(bytes.NewReader(data)); err != nil {
+		t.Fatalf("Decode: %v", err)
+	}
+}
+
+func TestExtraneousData(t *testing.T) {
+	// Encode a 1x1 red image.
+	src := image.NewRGBA(image.Rect(0, 0, 1, 1))
+	src.Set(0, 0, color.RGBA{0xff, 0x00, 0x00, 0xff})
+	buf := new(bytes.Buffer)
+	if err := Encode(buf, src, nil); err != nil {
+		t.Fatalf("encode: %v", err)
+	}
+	enc := buf.String()
+	// Sanity check that the encoded JPEG is long enough, that it ends in a
+	// "\xff\xd9" EOI marker, and that it contains a "\xff\xda" SOS marker
+	// somewhere in the final 64 bytes.
+	if len(enc) < 64 {
+		t.Fatalf("encoded JPEG is too short: %d bytes", len(enc))
+	}
+	if got, want := enc[len(enc)-2:], "\xff\xd9"; got != want {
+		t.Fatalf("encoded JPEG ends with %q, want %q", got, want)
+	}
+	if s := enc[len(enc)-64:]; !strings.Contains(s, "\xff\xda") {
+		t.Fatalf("encoded JPEG does not contain a SOS marker (ff da) near the end: % x", s)
+	}
+	// Test that adding some random junk between the SOS marker and the
+	// EOI marker does not affect the decoding.
+	rnd := rand.New(rand.NewSource(1))
+	for i, nerr := 0, 0; i < 1000 && nerr < 10; i++ {
+		buf.Reset()
+		// Write all but the trailing "\xff\xd9" EOI marker.
+		buf.WriteString(enc[:len(enc)-2])
+		// Write some random extraneous data.
+		for n := rnd.Intn(10); n > 0; n-- {
+			if x := byte(rnd.Intn(256)); x != 0xff {
+				buf.WriteByte(x)
+			} else {
+				// The JPEG format escapes a SOS 0xff data byte as "\xff\x00".
+				buf.WriteString("\xff\x00")
+			}
+		}
+		// Write the "\xff\xd9" EOI marker.
+		buf.WriteString("\xff\xd9")
+
+		// Check that we can still decode the resultant image.
+		got, err := Decode(buf)
+		if err != nil {
+			t.Errorf("could not decode image #%d: %v", i, err)
+			nerr++
+			continue
+		}
+		if got.Bounds() != src.Bounds() {
+			t.Errorf("image #%d, bounds differ: %v and %v", i, got.Bounds(), src.Bounds())
+			nerr++
+			continue
+		}
+		if averageDelta(got, src) > 2<<8 {
+			t.Errorf("image #%d changed too much after a round trip", i)
+			nerr++
+			continue
+		}
+	}
+}
+
+func benchmarkDecode(b *testing.B, filename string) {
+	data, err := os.ReadFile(filename)
+	if err != nil {
+		b.Fatal(err)
+	}
+	cfg, err := DecodeConfig(bytes.NewReader(data))
+	if err != nil {
+		b.Fatal(err)
+	}
+	b.SetBytes(int64(cfg.Width * cfg.Height * 4))
+	b.ReportAllocs()
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		Decode(bytes.NewReader(data))
+	}
+}
+
+func BenchmarkDecodeBaseline(b *testing.B) {
+	benchmarkDecode(b, "../testdata/video-001.jpeg")
+}
+
+func BenchmarkDecodeProgressive(b *testing.B) {
+	benchmarkDecode(b, "../testdata/video-001.progressive.jpeg")
+}
diff --git a/src/image/jpeg/scan.go b/src/image/jpeg/scan.go
new file mode 100644
index 0000000..94f3d3a
--- /dev/null
+++ b/src/image/jpeg/scan.go
@@ -0,0 +1,523 @@
+// Copyright 2012 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 jpeg
+
+import (
+	"image"
+)
+
+// makeImg allocates and initializes the destination image.
+func (d *decoder) makeImg(mxx, myy int) {
+	if d.nComp == 1 {
+		m := image.NewGray(image.Rect(0, 0, 8*mxx, 8*myy))
+		d.img1 = m.SubImage(image.Rect(0, 0, d.width, d.height)).(*image.Gray)
+		return
+	}
+
+	h0 := d.comp[0].h
+	v0 := d.comp[0].v
+	hRatio := h0 / d.comp[1].h
+	vRatio := v0 / d.comp[1].v
+	var subsampleRatio image.YCbCrSubsampleRatio
+	switch hRatio<<4 | vRatio {
+	case 0x11:
+		subsampleRatio = image.YCbCrSubsampleRatio444
+	case 0x12:
+		subsampleRatio = image.YCbCrSubsampleRatio440
+	case 0x21:
+		subsampleRatio = image.YCbCrSubsampleRatio422
+	case 0x22:
+		subsampleRatio = image.YCbCrSubsampleRatio420
+	case 0x41:
+		subsampleRatio = image.YCbCrSubsampleRatio411
+	case 0x42:
+		subsampleRatio = image.YCbCrSubsampleRatio410
+	default:
+		panic("unreachable")
+	}
+	m := image.NewYCbCr(image.Rect(0, 0, 8*h0*mxx, 8*v0*myy), subsampleRatio)
+	d.img3 = m.SubImage(image.Rect(0, 0, d.width, d.height)).(*image.YCbCr)
+
+	if d.nComp == 4 {
+		h3, v3 := d.comp[3].h, d.comp[3].v
+		d.blackPix = make([]byte, 8*h3*mxx*8*v3*myy)
+		d.blackStride = 8 * h3 * mxx
+	}
+}
+
+// Specified in section B.2.3.
+func (d *decoder) processSOS(n int) error {
+	if d.nComp == 0 {
+		return FormatError("missing SOF marker")
+	}
+	if n < 6 || 4+2*d.nComp < n || n%2 != 0 {
+		return FormatError("SOS has wrong length")
+	}
+	if err := d.readFull(d.tmp[:n]); err != nil {
+		return err
+	}
+	nComp := int(d.tmp[0])
+	if n != 4+2*nComp {
+		return FormatError("SOS length inconsistent with number of components")
+	}
+	var scan [maxComponents]struct {
+		compIndex uint8
+		td        uint8 // DC table selector.
+		ta        uint8 // AC table selector.
+	}
+	totalHV := 0
+	for i := 0; i < nComp; i++ {
+		cs := d.tmp[1+2*i] // Component selector.
+		compIndex := -1
+		for j, comp := range d.comp[:d.nComp] {
+			if cs == comp.c {
+				compIndex = j
+			}
+		}
+		if compIndex < 0 {
+			return FormatError("unknown component selector")
+		}
+		scan[i].compIndex = uint8(compIndex)
+		// Section B.2.3 states that "the value of Cs_j shall be different from
+		// the values of Cs_1 through Cs_(j-1)". Since we have previously
+		// verified that a frame's component identifiers (C_i values in section
+		// B.2.2) are unique, it suffices to check that the implicit indexes
+		// into d.comp are unique.
+		for j := 0; j < i; j++ {
+			if scan[i].compIndex == scan[j].compIndex {
+				return FormatError("repeated component selector")
+			}
+		}
+		totalHV += d.comp[compIndex].h * d.comp[compIndex].v
+
+		// The baseline t <= 1 restriction is specified in table B.3.
+		scan[i].td = d.tmp[2+2*i] >> 4
+		if t := scan[i].td; t > maxTh || (d.baseline && t > 1) {
+			return FormatError("bad Td value")
+		}
+		scan[i].ta = d.tmp[2+2*i] & 0x0f
+		if t := scan[i].ta; t > maxTh || (d.baseline && t > 1) {
+			return FormatError("bad Ta value")
+		}
+	}
+	// Section B.2.3 states that if there is more than one component then the
+	// total H*V values in a scan must be <= 10.
+	if d.nComp > 1 && totalHV > 10 {
+		return FormatError("total sampling factors too large")
+	}
+
+	// zigStart and zigEnd are the spectral selection bounds.
+	// ah and al are the successive approximation high and low values.
+	// The spec calls these values Ss, Se, Ah and Al.
+	//
+	// For progressive JPEGs, these are the two more-or-less independent
+	// aspects of progression. Spectral selection progression is when not
+	// all of a block's 64 DCT coefficients are transmitted in one pass.
+	// For example, three passes could transmit coefficient 0 (the DC
+	// component), coefficients 1-5, and coefficients 6-63, in zig-zag
+	// order. Successive approximation is when not all of the bits of a
+	// band of coefficients are transmitted in one pass. For example,
+	// three passes could transmit the 6 most significant bits, followed
+	// by the second-least significant bit, followed by the least
+	// significant bit.
+	//
+	// For sequential JPEGs, these parameters are hard-coded to 0/63/0/0, as
+	// per table B.3.
+	zigStart, zigEnd, ah, al := int32(0), int32(blockSize-1), uint32(0), uint32(0)
+	if d.progressive {
+		zigStart = int32(d.tmp[1+2*nComp])
+		zigEnd = int32(d.tmp[2+2*nComp])
+		ah = uint32(d.tmp[3+2*nComp] >> 4)
+		al = uint32(d.tmp[3+2*nComp] & 0x0f)
+		if (zigStart == 0 && zigEnd != 0) || zigStart > zigEnd || blockSize <= zigEnd {
+			return FormatError("bad spectral selection bounds")
+		}
+		if zigStart != 0 && nComp != 1 {
+			return FormatError("progressive AC coefficients for more than one component")
+		}
+		if ah != 0 && ah != al+1 {
+			return FormatError("bad successive approximation values")
+		}
+	}
+
+	// mxx and myy are the number of MCUs (Minimum Coded Units) in the image.
+	h0, v0 := d.comp[0].h, d.comp[0].v // The h and v values from the Y components.
+	mxx := (d.width + 8*h0 - 1) / (8 * h0)
+	myy := (d.height + 8*v0 - 1) / (8 * v0)
+	if d.img1 == nil && d.img3 == nil {
+		d.makeImg(mxx, myy)
+	}
+	if d.progressive {
+		for i := 0; i < nComp; i++ {
+			compIndex := scan[i].compIndex
+			if d.progCoeffs[compIndex] == nil {
+				d.progCoeffs[compIndex] = make([]block, mxx*myy*d.comp[compIndex].h*d.comp[compIndex].v)
+			}
+		}
+	}
+
+	d.bits = bits{}
+	mcu, expectedRST := 0, uint8(rst0Marker)
+	var (
+		// b is the decoded coefficients, in natural (not zig-zag) order.
+		b  block
+		dc [maxComponents]int32
+		// bx and by are the location of the current block, in units of 8x8
+		// blocks: the third block in the first row has (bx, by) = (2, 0).
+		bx, by     int
+		blockCount int
+	)
+	for my := 0; my < myy; my++ {
+		for mx := 0; mx < mxx; mx++ {
+			for i := 0; i < nComp; i++ {
+				compIndex := scan[i].compIndex
+				hi := d.comp[compIndex].h
+				vi := d.comp[compIndex].v
+				for j := 0; j < hi*vi; j++ {
+					// The blocks are traversed one MCU at a time. For 4:2:0 chroma
+					// subsampling, there are four Y 8x8 blocks in every 16x16 MCU.
+					//
+					// For a sequential 32x16 pixel image, the Y blocks visiting order is:
+					//	0 1 4 5
+					//	2 3 6 7
+					//
+					// For progressive images, the interleaved scans (those with nComp > 1)
+					// are traversed as above, but non-interleaved scans are traversed left
+					// to right, top to bottom:
+					//	0 1 2 3
+					//	4 5 6 7
+					// Only DC scans (zigStart == 0) can be interleaved. AC scans must have
+					// only one component.
+					//
+					// To further complicate matters, for non-interleaved scans, there is no
+					// data for any blocks that are inside the image at the MCU level but
+					// outside the image at the pixel level. For example, a 24x16 pixel 4:2:0
+					// progressive image consists of two 16x16 MCUs. The interleaved scans
+					// will process 8 Y blocks:
+					//	0 1 4 5
+					//	2 3 6 7
+					// The non-interleaved scans will process only 6 Y blocks:
+					//	0 1 2
+					//	3 4 5
+					if nComp != 1 {
+						bx = hi*mx + j%hi
+						by = vi*my + j/hi
+					} else {
+						q := mxx * hi
+						bx = blockCount % q
+						by = blockCount / q
+						blockCount++
+						if bx*8 >= d.width || by*8 >= d.height {
+							continue
+						}
+					}
+
+					// Load the previous partially decoded coefficients, if applicable.
+					if d.progressive {
+						b = d.progCoeffs[compIndex][by*mxx*hi+bx]
+					} else {
+						b = block{}
+					}
+
+					if ah != 0 {
+						if err := d.refine(&b, &d.huff[acTable][scan[i].ta], zigStart, zigEnd, 1<<al); err != nil {
+							return err
+						}
+					} else {
+						zig := zigStart
+						if zig == 0 {
+							zig++
+							// Decode the DC coefficient, as specified in section F.2.2.1.
+							value, err := d.decodeHuffman(&d.huff[dcTable][scan[i].td])
+							if err != nil {
+								return err
+							}
+							if value > 16 {
+								return UnsupportedError("excessive DC component")
+							}
+							dcDelta, err := d.receiveExtend(value)
+							if err != nil {
+								return err
+							}
+							dc[compIndex] += dcDelta
+							b[0] = dc[compIndex] << al
+						}
+
+						if zig <= zigEnd && d.eobRun > 0 {
+							d.eobRun--
+						} else {
+							// Decode the AC coefficients, as specified in section F.2.2.2.
+							huff := &d.huff[acTable][scan[i].ta]
+							for ; zig <= zigEnd; zig++ {
+								value, err := d.decodeHuffman(huff)
+								if err != nil {
+									return err
+								}
+								val0 := value >> 4
+								val1 := value & 0x0f
+								if val1 != 0 {
+									zig += int32(val0)
+									if zig > zigEnd {
+										break
+									}
+									ac, err := d.receiveExtend(val1)
+									if err != nil {
+										return err
+									}
+									b[unzig[zig]] = ac << al
+								} else {
+									if val0 != 0x0f {
+										d.eobRun = uint16(1 << val0)
+										if val0 != 0 {
+											bits, err := d.decodeBits(int32(val0))
+											if err != nil {
+												return err
+											}
+											d.eobRun |= uint16(bits)
+										}
+										d.eobRun--
+										break
+									}
+									zig += 0x0f
+								}
+							}
+						}
+					}
+
+					if d.progressive {
+						// Save the coefficients.
+						d.progCoeffs[compIndex][by*mxx*hi+bx] = b
+						// At this point, we could call reconstructBlock to dequantize and perform the
+						// inverse DCT, to save early stages of a progressive image to the *image.YCbCr
+						// buffers (the whole point of progressive encoding), but in Go, the jpeg.Decode
+						// function does not return until the entire image is decoded, so we "continue"
+						// here to avoid wasted computation. Instead, reconstructBlock is called on each
+						// accumulated block by the reconstructProgressiveImage method after all of the
+						// SOS markers are processed.
+						continue
+					}
+					if err := d.reconstructBlock(&b, bx, by, int(compIndex)); err != nil {
+						return err
+					}
+				} // for j
+			} // for i
+			mcu++
+			if d.ri > 0 && mcu%d.ri == 0 && mcu < mxx*myy {
+				// A more sophisticated decoder could use RST[0-7] markers to resynchronize from corrupt input,
+				// but this one assumes well-formed input, and hence the restart marker follows immediately.
+				if err := d.readFull(d.tmp[:2]); err != nil {
+					return err
+				}
+
+				// Section F.1.2.3 says that "Byte alignment of markers is
+				// achieved by padding incomplete bytes with 1-bits. If padding
+				// with 1-bits creates a X’FF’ value, a zero byte is stuffed
+				// before adding the marker."
+				//
+				// Seeing "\xff\x00" here is not spec compliant, as we are not
+				// expecting an *incomplete* byte (that needed padding). Still,
+				// some real world encoders (see golang.org/issue/28717) insert
+				// it, so we accept it and re-try the 2 byte read.
+				//
+				// libjpeg issues a warning (but not an error) for this:
+				// https://github.com/LuaDist/libjpeg/blob/6c0fcb8ddee365e7abc4d332662b06900612e923/jdmarker.c#L1041-L1046
+				if d.tmp[0] == 0xff && d.tmp[1] == 0x00 {
+					if err := d.readFull(d.tmp[:2]); err != nil {
+						return err
+					}
+				}
+
+				if d.tmp[0] != 0xff || d.tmp[1] != expectedRST {
+					return FormatError("bad RST marker")
+				}
+				expectedRST++
+				if expectedRST == rst7Marker+1 {
+					expectedRST = rst0Marker
+				}
+				// Reset the Huffman decoder.
+				d.bits = bits{}
+				// Reset the DC components, as per section F.2.1.3.1.
+				dc = [maxComponents]int32{}
+				// Reset the progressive decoder state, as per section G.1.2.2.
+				d.eobRun = 0
+			}
+		} // for mx
+	} // for my
+
+	return nil
+}
+
+// refine decodes a successive approximation refinement block, as specified in
+// section G.1.2.
+func (d *decoder) refine(b *block, h *huffman, zigStart, zigEnd, delta int32) error {
+	// Refining a DC component is trivial.
+	if zigStart == 0 {
+		if zigEnd != 0 {
+			panic("unreachable")
+		}
+		bit, err := d.decodeBit()
+		if err != nil {
+			return err
+		}
+		if bit {
+			b[0] |= delta
+		}
+		return nil
+	}
+
+	// Refining AC components is more complicated; see sections G.1.2.2 and G.1.2.3.
+	zig := zigStart
+	if d.eobRun == 0 {
+	loop:
+		for ; zig <= zigEnd; zig++ {
+			z := int32(0)
+			value, err := d.decodeHuffman(h)
+			if err != nil {
+				return err
+			}
+			val0 := value >> 4
+			val1 := value & 0x0f
+
+			switch val1 {
+			case 0:
+				if val0 != 0x0f {
+					d.eobRun = uint16(1 << val0)
+					if val0 != 0 {
+						bits, err := d.decodeBits(int32(val0))
+						if err != nil {
+							return err
+						}
+						d.eobRun |= uint16(bits)
+					}
+					break loop
+				}
+			case 1:
+				z = delta
+				bit, err := d.decodeBit()
+				if err != nil {
+					return err
+				}
+				if !bit {
+					z = -z
+				}
+			default:
+				return FormatError("unexpected Huffman code")
+			}
+
+			zig, err = d.refineNonZeroes(b, zig, zigEnd, int32(val0), delta)
+			if err != nil {
+				return err
+			}
+			if zig > zigEnd {
+				return FormatError("too many coefficients")
+			}
+			if z != 0 {
+				b[unzig[zig]] = z
+			}
+		}
+	}
+	if d.eobRun > 0 {
+		d.eobRun--
+		if _, err := d.refineNonZeroes(b, zig, zigEnd, -1, delta); err != nil {
+			return err
+		}
+	}
+	return nil
+}
+
+// refineNonZeroes refines non-zero entries of b in zig-zag order. If nz >= 0,
+// the first nz zero entries are skipped over.
+func (d *decoder) refineNonZeroes(b *block, zig, zigEnd, nz, delta int32) (int32, error) {
+	for ; zig <= zigEnd; zig++ {
+		u := unzig[zig]
+		if b[u] == 0 {
+			if nz == 0 {
+				break
+			}
+			nz--
+			continue
+		}
+		bit, err := d.decodeBit()
+		if err != nil {
+			return 0, err
+		}
+		if !bit {
+			continue
+		}
+		if b[u] >= 0 {
+			b[u] += delta
+		} else {
+			b[u] -= delta
+		}
+	}
+	return zig, nil
+}
+
+func (d *decoder) reconstructProgressiveImage() error {
+	// The h0, mxx, by and bx variables have the same meaning as in the
+	// processSOS method.
+	h0 := d.comp[0].h
+	mxx := (d.width + 8*h0 - 1) / (8 * h0)
+	for i := 0; i < d.nComp; i++ {
+		if d.progCoeffs[i] == nil {
+			continue
+		}
+		v := 8 * d.comp[0].v / d.comp[i].v
+		h := 8 * d.comp[0].h / d.comp[i].h
+		stride := mxx * d.comp[i].h
+		for by := 0; by*v < d.height; by++ {
+			for bx := 0; bx*h < d.width; bx++ {
+				if err := d.reconstructBlock(&d.progCoeffs[i][by*stride+bx], bx, by, i); err != nil {
+					return err
+				}
+			}
+		}
+	}
+	return nil
+}
+
+// reconstructBlock dequantizes, performs the inverse DCT and stores the block
+// to the image.
+func (d *decoder) reconstructBlock(b *block, bx, by, compIndex int) error {
+	qt := &d.quant[d.comp[compIndex].tq]
+	for zig := 0; zig < blockSize; zig++ {
+		b[unzig[zig]] *= qt[zig]
+	}
+	idct(b)
+	dst, stride := []byte(nil), 0
+	if d.nComp == 1 {
+		dst, stride = d.img1.Pix[8*(by*d.img1.Stride+bx):], d.img1.Stride
+	} else {
+		switch compIndex {
+		case 0:
+			dst, stride = d.img3.Y[8*(by*d.img3.YStride+bx):], d.img3.YStride
+		case 1:
+			dst, stride = d.img3.Cb[8*(by*d.img3.CStride+bx):], d.img3.CStride
+		case 2:
+			dst, stride = d.img3.Cr[8*(by*d.img3.CStride+bx):], d.img3.CStride
+		case 3:
+			dst, stride = d.blackPix[8*(by*d.blackStride+bx):], d.blackStride
+		default:
+			return UnsupportedError("too many components")
+		}
+	}
+	// Level shift by +128, clip to [0, 255], and write to dst.
+	for y := 0; y < 8; y++ {
+		y8 := y * 8
+		yStride := y * stride
+		for x := 0; x < 8; x++ {
+			c := b[y8+x]
+			if c < -128 {
+				c = 0
+			} else if c > 127 {
+				c = 255
+			} else {
+				c += 128
+			}
+			dst[yStride+x] = uint8(c)
+		}
+	}
+	return nil
+}
diff --git a/src/image/jpeg/writer.go b/src/image/jpeg/writer.go
new file mode 100644
index 0000000..0027f78
--- /dev/null
+++ b/src/image/jpeg/writer.go
@@ -0,0 +1,641 @@
+// Copyright 2011 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 jpeg
+
+import (
+	"bufio"
+	"errors"
+	"image"
+	"image/color"
+	"io"
+)
+
+// min returns the minimum of two integers.
+func min(x, y int) int {
+	if x < y {
+		return x
+	}
+	return y
+}
+
+// div returns a/b rounded to the nearest integer, instead of rounded to zero.
+func div(a, b int32) int32 {
+	if a >= 0 {
+		return (a + (b >> 1)) / b
+	}
+	return -((-a + (b >> 1)) / b)
+}
+
+// bitCount counts the number of bits needed to hold an integer.
+var bitCount = [256]byte{
+	0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,
+	5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
+	6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
+	6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
+	7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+	7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+	7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+	7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+	8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+	8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+	8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+	8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+	8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+	8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+	8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+	8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
+}
+
+type quantIndex int
+
+const (
+	quantIndexLuminance quantIndex = iota
+	quantIndexChrominance
+	nQuantIndex
+)
+
+// unscaledQuant are the unscaled quantization tables in zig-zag order. Each
+// encoder copies and scales the tables according to its quality parameter.
+// The values are derived from section K.1 after converting from natural to
+// zig-zag order.
+var unscaledQuant = [nQuantIndex][blockSize]byte{
+	// Luminance.
+	{
+		16, 11, 12, 14, 12, 10, 16, 14,
+		13, 14, 18, 17, 16, 19, 24, 40,
+		26, 24, 22, 22, 24, 49, 35, 37,
+		29, 40, 58, 51, 61, 60, 57, 51,
+		56, 55, 64, 72, 92, 78, 64, 68,
+		87, 69, 55, 56, 80, 109, 81, 87,
+		95, 98, 103, 104, 103, 62, 77, 113,
+		121, 112, 100, 120, 92, 101, 103, 99,
+	},
+	// Chrominance.
+	{
+		17, 18, 18, 24, 21, 24, 47, 26,
+		26, 47, 99, 66, 56, 66, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+	},
+}
+
+type huffIndex int
+
+const (
+	huffIndexLuminanceDC huffIndex = iota
+	huffIndexLuminanceAC
+	huffIndexChrominanceDC
+	huffIndexChrominanceAC
+	nHuffIndex
+)
+
+// huffmanSpec specifies a Huffman encoding.
+type huffmanSpec struct {
+	// count[i] is the number of codes of length i bits.
+	count [16]byte
+	// value[i] is the decoded value of the i'th codeword.
+	value []byte
+}
+
+// theHuffmanSpec is the Huffman encoding specifications.
+// This encoder uses the same Huffman encoding for all images.
+var theHuffmanSpec = [nHuffIndex]huffmanSpec{
+	// Luminance DC.
+	{
+		[16]byte{0, 1, 5, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0},
+		[]byte{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11},
+	},
+	// Luminance AC.
+	{
+		[16]byte{0, 2, 1, 3, 3, 2, 4, 3, 5, 5, 4, 4, 0, 0, 1, 125},
+		[]byte{
+			0x01, 0x02, 0x03, 0x00, 0x04, 0x11, 0x05, 0x12,
+			0x21, 0x31, 0x41, 0x06, 0x13, 0x51, 0x61, 0x07,
+			0x22, 0x71, 0x14, 0x32, 0x81, 0x91, 0xa1, 0x08,
+			0x23, 0x42, 0xb1, 0xc1, 0x15, 0x52, 0xd1, 0xf0,
+			0x24, 0x33, 0x62, 0x72, 0x82, 0x09, 0x0a, 0x16,
+			0x17, 0x18, 0x19, 0x1a, 0x25, 0x26, 0x27, 0x28,
+			0x29, 0x2a, 0x34, 0x35, 0x36, 0x37, 0x38, 0x39,
+			0x3a, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48, 0x49,
+			0x4a, 0x53, 0x54, 0x55, 0x56, 0x57, 0x58, 0x59,
+			0x5a, 0x63, 0x64, 0x65, 0x66, 0x67, 0x68, 0x69,
+			0x6a, 0x73, 0x74, 0x75, 0x76, 0x77, 0x78, 0x79,
+			0x7a, 0x83, 0x84, 0x85, 0x86, 0x87, 0x88, 0x89,
+			0x8a, 0x92, 0x93, 0x94, 0x95, 0x96, 0x97, 0x98,
+			0x99, 0x9a, 0xa2, 0xa3, 0xa4, 0xa5, 0xa6, 0xa7,
+			0xa8, 0xa9, 0xaa, 0xb2, 0xb3, 0xb4, 0xb5, 0xb6,
+			0xb7, 0xb8, 0xb9, 0xba, 0xc2, 0xc3, 0xc4, 0xc5,
+			0xc6, 0xc7, 0xc8, 0xc9, 0xca, 0xd2, 0xd3, 0xd4,
+			0xd5, 0xd6, 0xd7, 0xd8, 0xd9, 0xda, 0xe1, 0xe2,
+			0xe3, 0xe4, 0xe5, 0xe6, 0xe7, 0xe8, 0xe9, 0xea,
+			0xf1, 0xf2, 0xf3, 0xf4, 0xf5, 0xf6, 0xf7, 0xf8,
+			0xf9, 0xfa,
+		},
+	},
+	// Chrominance DC.
+	{
+		[16]byte{0, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0},
+		[]byte{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11},
+	},
+	// Chrominance AC.
+	{
+		[16]byte{0, 2, 1, 2, 4, 4, 3, 4, 7, 5, 4, 4, 0, 1, 2, 119},
+		[]byte{
+			0x00, 0x01, 0x02, 0x03, 0x11, 0x04, 0x05, 0x21,
+			0x31, 0x06, 0x12, 0x41, 0x51, 0x07, 0x61, 0x71,
+			0x13, 0x22, 0x32, 0x81, 0x08, 0x14, 0x42, 0x91,
+			0xa1, 0xb1, 0xc1, 0x09, 0x23, 0x33, 0x52, 0xf0,
+			0x15, 0x62, 0x72, 0xd1, 0x0a, 0x16, 0x24, 0x34,
+			0xe1, 0x25, 0xf1, 0x17, 0x18, 0x19, 0x1a, 0x26,
+			0x27, 0x28, 0x29, 0x2a, 0x35, 0x36, 0x37, 0x38,
+			0x39, 0x3a, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48,
+			0x49, 0x4a, 0x53, 0x54, 0x55, 0x56, 0x57, 0x58,
+			0x59, 0x5a, 0x63, 0x64, 0x65, 0x66, 0x67, 0x68,
+			0x69, 0x6a, 0x73, 0x74, 0x75, 0x76, 0x77, 0x78,
+			0x79, 0x7a, 0x82, 0x83, 0x84, 0x85, 0x86, 0x87,
+			0x88, 0x89, 0x8a, 0x92, 0x93, 0x94, 0x95, 0x96,
+			0x97, 0x98, 0x99, 0x9a, 0xa2, 0xa3, 0xa4, 0xa5,
+			0xa6, 0xa7, 0xa8, 0xa9, 0xaa, 0xb2, 0xb3, 0xb4,
+			0xb5, 0xb6, 0xb7, 0xb8, 0xb9, 0xba, 0xc2, 0xc3,
+			0xc4, 0xc5, 0xc6, 0xc7, 0xc8, 0xc9, 0xca, 0xd2,
+			0xd3, 0xd4, 0xd5, 0xd6, 0xd7, 0xd8, 0xd9, 0xda,
+			0xe2, 0xe3, 0xe4, 0xe5, 0xe6, 0xe7, 0xe8, 0xe9,
+			0xea, 0xf2, 0xf3, 0xf4, 0xf5, 0xf6, 0xf7, 0xf8,
+			0xf9, 0xfa,
+		},
+	},
+}
+
+// huffmanLUT is a compiled look-up table representation of a huffmanSpec.
+// Each value maps to a uint32 of which the 8 most significant bits hold the
+// codeword size in bits and the 24 least significant bits hold the codeword.
+// The maximum codeword size is 16 bits.
+type huffmanLUT []uint32
+
+func (h *huffmanLUT) init(s huffmanSpec) {
+	maxValue := 0
+	for _, v := range s.value {
+		if int(v) > maxValue {
+			maxValue = int(v)
+		}
+	}
+	*h = make([]uint32, maxValue+1)
+	code, k := uint32(0), 0
+	for i := 0; i < len(s.count); i++ {
+		nBits := uint32(i+1) << 24
+		for j := uint8(0); j < s.count[i]; j++ {
+			(*h)[s.value[k]] = nBits | code
+			code++
+			k++
+		}
+		code <<= 1
+	}
+}
+
+// theHuffmanLUT are compiled representations of theHuffmanSpec.
+var theHuffmanLUT [4]huffmanLUT
+
+func init() {
+	for i, s := range theHuffmanSpec {
+		theHuffmanLUT[i].init(s)
+	}
+}
+
+// writer is a buffered writer.
+type writer interface {
+	Flush() error
+	io.Writer
+	io.ByteWriter
+}
+
+// encoder encodes an image to the JPEG format.
+type encoder struct {
+	// w is the writer to write to. err is the first error encountered during
+	// writing. All attempted writes after the first error become no-ops.
+	w   writer
+	err error
+	// buf is a scratch buffer.
+	buf [16]byte
+	// bits and nBits are accumulated bits to write to w.
+	bits, nBits uint32
+	// quant is the scaled quantization tables, in zig-zag order.
+	quant [nQuantIndex][blockSize]byte
+}
+
+func (e *encoder) flush() {
+	if e.err != nil {
+		return
+	}
+	e.err = e.w.Flush()
+}
+
+func (e *encoder) write(p []byte) {
+	if e.err != nil {
+		return
+	}
+	_, e.err = e.w.Write(p)
+}
+
+func (e *encoder) writeByte(b byte) {
+	if e.err != nil {
+		return
+	}
+	e.err = e.w.WriteByte(b)
+}
+
+// emit emits the least significant nBits bits of bits to the bit-stream.
+// The precondition is bits < 1<<nBits && nBits <= 16.
+func (e *encoder) emit(bits, nBits uint32) {
+	nBits += e.nBits
+	bits <<= 32 - nBits
+	bits |= e.bits
+	for nBits >= 8 {
+		b := uint8(bits >> 24)
+		e.writeByte(b)
+		if b == 0xff {
+			e.writeByte(0x00)
+		}
+		bits <<= 8
+		nBits -= 8
+	}
+	e.bits, e.nBits = bits, nBits
+}
+
+// emitHuff emits the given value with the given Huffman encoder.
+func (e *encoder) emitHuff(h huffIndex, value int32) {
+	x := theHuffmanLUT[h][value]
+	e.emit(x&(1<<24-1), x>>24)
+}
+
+// emitHuffRLE emits a run of runLength copies of value encoded with the given
+// Huffman encoder.
+func (e *encoder) emitHuffRLE(h huffIndex, runLength, value int32) {
+	a, b := value, value
+	if a < 0 {
+		a, b = -value, value-1
+	}
+	var nBits uint32
+	if a < 0x100 {
+		nBits = uint32(bitCount[a])
+	} else {
+		nBits = 8 + uint32(bitCount[a>>8])
+	}
+	e.emitHuff(h, runLength<<4|int32(nBits))
+	if nBits > 0 {
+		e.emit(uint32(b)&(1<<nBits-1), nBits)
+	}
+}
+
+// writeMarkerHeader writes the header for a marker with the given length.
+func (e *encoder) writeMarkerHeader(marker uint8, markerlen int) {
+	e.buf[0] = 0xff
+	e.buf[1] = marker
+	e.buf[2] = uint8(markerlen >> 8)
+	e.buf[3] = uint8(markerlen & 0xff)
+	e.write(e.buf[:4])
+}
+
+// writeDQT writes the Define Quantization Table marker.
+func (e *encoder) writeDQT() {
+	const markerlen = 2 + int(nQuantIndex)*(1+blockSize)
+	e.writeMarkerHeader(dqtMarker, markerlen)
+	for i := range e.quant {
+		e.writeByte(uint8(i))
+		e.write(e.quant[i][:])
+	}
+}
+
+// writeSOF0 writes the Start Of Frame (Baseline Sequential) marker.
+func (e *encoder) writeSOF0(size image.Point, nComponent int) {
+	markerlen := 8 + 3*nComponent
+	e.writeMarkerHeader(sof0Marker, markerlen)
+	e.buf[0] = 8 // 8-bit color.
+	e.buf[1] = uint8(size.Y >> 8)
+	e.buf[2] = uint8(size.Y & 0xff)
+	e.buf[3] = uint8(size.X >> 8)
+	e.buf[4] = uint8(size.X & 0xff)
+	e.buf[5] = uint8(nComponent)
+	if nComponent == 1 {
+		e.buf[6] = 1
+		// No subsampling for grayscale image.
+		e.buf[7] = 0x11
+		e.buf[8] = 0x00
+	} else {
+		for i := 0; i < nComponent; i++ {
+			e.buf[3*i+6] = uint8(i + 1)
+			// We use 4:2:0 chroma subsampling.
+			e.buf[3*i+7] = "\x22\x11\x11"[i]
+			e.buf[3*i+8] = "\x00\x01\x01"[i]
+		}
+	}
+	e.write(e.buf[:3*(nComponent-1)+9])
+}
+
+// writeDHT writes the Define Huffman Table marker.
+func (e *encoder) writeDHT(nComponent int) {
+	markerlen := 2
+	specs := theHuffmanSpec[:]
+	if nComponent == 1 {
+		// Drop the Chrominance tables.
+		specs = specs[:2]
+	}
+	for _, s := range specs {
+		markerlen += 1 + 16 + len(s.value)
+	}
+	e.writeMarkerHeader(dhtMarker, markerlen)
+	for i, s := range specs {
+		e.writeByte("\x00\x10\x01\x11"[i])
+		e.write(s.count[:])
+		e.write(s.value)
+	}
+}
+
+// writeBlock writes a block of pixel data using the given quantization table,
+// returning the post-quantized DC value of the DCT-transformed block. b is in
+// natural (not zig-zag) order.
+func (e *encoder) writeBlock(b *block, q quantIndex, prevDC int32) int32 {
+	fdct(b)
+	// Emit the DC delta.
+	dc := div(b[0], 8*int32(e.quant[q][0]))
+	e.emitHuffRLE(huffIndex(2*q+0), 0, dc-prevDC)
+	// Emit the AC components.
+	h, runLength := huffIndex(2*q+1), int32(0)
+	for zig := 1; zig < blockSize; zig++ {
+		ac := div(b[unzig[zig]], 8*int32(e.quant[q][zig]))
+		if ac == 0 {
+			runLength++
+		} else {
+			for runLength > 15 {
+				e.emitHuff(h, 0xf0)
+				runLength -= 16
+			}
+			e.emitHuffRLE(h, runLength, ac)
+			runLength = 0
+		}
+	}
+	if runLength > 0 {
+		e.emitHuff(h, 0x00)
+	}
+	return dc
+}
+
+// toYCbCr converts the 8x8 region of m whose top-left corner is p to its
+// YCbCr values.
+func toYCbCr(m image.Image, p image.Point, yBlock, cbBlock, crBlock *block) {
+	b := m.Bounds()
+	xmax := b.Max.X - 1
+	ymax := b.Max.Y - 1
+	for j := 0; j < 8; j++ {
+		for i := 0; i < 8; i++ {
+			r, g, b, _ := m.At(min(p.X+i, xmax), min(p.Y+j, ymax)).RGBA()
+			yy, cb, cr := color.RGBToYCbCr(uint8(r>>8), uint8(g>>8), uint8(b>>8))
+			yBlock[8*j+i] = int32(yy)
+			cbBlock[8*j+i] = int32(cb)
+			crBlock[8*j+i] = int32(cr)
+		}
+	}
+}
+
+// grayToY stores the 8x8 region of m whose top-left corner is p in yBlock.
+func grayToY(m *image.Gray, p image.Point, yBlock *block) {
+	b := m.Bounds()
+	xmax := b.Max.X - 1
+	ymax := b.Max.Y - 1
+	pix := m.Pix
+	for j := 0; j < 8; j++ {
+		for i := 0; i < 8; i++ {
+			idx := m.PixOffset(min(p.X+i, xmax), min(p.Y+j, ymax))
+			yBlock[8*j+i] = int32(pix[idx])
+		}
+	}
+}
+
+// rgbaToYCbCr is a specialized version of toYCbCr for image.RGBA images.
+func rgbaToYCbCr(m *image.RGBA, p image.Point, yBlock, cbBlock, crBlock *block) {
+	b := m.Bounds()
+	xmax := b.Max.X - 1
+	ymax := b.Max.Y - 1
+	for j := 0; j < 8; j++ {
+		sj := p.Y + j
+		if sj > ymax {
+			sj = ymax
+		}
+		offset := (sj-b.Min.Y)*m.Stride - b.Min.X*4
+		for i := 0; i < 8; i++ {
+			sx := p.X + i
+			if sx > xmax {
+				sx = xmax
+			}
+			pix := m.Pix[offset+sx*4:]
+			yy, cb, cr := color.RGBToYCbCr(pix[0], pix[1], pix[2])
+			yBlock[8*j+i] = int32(yy)
+			cbBlock[8*j+i] = int32(cb)
+			crBlock[8*j+i] = int32(cr)
+		}
+	}
+}
+
+// yCbCrToYCbCr is a specialized version of toYCbCr for image.YCbCr images.
+func yCbCrToYCbCr(m *image.YCbCr, p image.Point, yBlock, cbBlock, crBlock *block) {
+	b := m.Bounds()
+	xmax := b.Max.X - 1
+	ymax := b.Max.Y - 1
+	for j := 0; j < 8; j++ {
+		sy := p.Y + j
+		if sy > ymax {
+			sy = ymax
+		}
+		for i := 0; i < 8; i++ {
+			sx := p.X + i
+			if sx > xmax {
+				sx = xmax
+			}
+			yi := m.YOffset(sx, sy)
+			ci := m.COffset(sx, sy)
+			yBlock[8*j+i] = int32(m.Y[yi])
+			cbBlock[8*j+i] = int32(m.Cb[ci])
+			crBlock[8*j+i] = int32(m.Cr[ci])
+		}
+	}
+}
+
+// scale scales the 16x16 region represented by the 4 src blocks to the 8x8
+// dst block.
+func scale(dst *block, src *[4]block) {
+	for i := 0; i < 4; i++ {
+		dstOff := (i&2)<<4 | (i&1)<<2
+		for y := 0; y < 4; y++ {
+			for x := 0; x < 4; x++ {
+				j := 16*y + 2*x
+				sum := src[i][j] + src[i][j+1] + src[i][j+8] + src[i][j+9]
+				dst[8*y+x+dstOff] = (sum + 2) >> 2
+			}
+		}
+	}
+}
+
+// sosHeaderY is the SOS marker "\xff\xda" followed by 8 bytes:
+//   - the marker length "\x00\x08",
+//   - the number of components "\x01",
+//   - component 1 uses DC table 0 and AC table 0 "\x01\x00",
+//   - the bytes "\x00\x3f\x00". Section B.2.3 of the spec says that for
+//     sequential DCTs, those bytes (8-bit Ss, 8-bit Se, 4-bit Ah, 4-bit Al)
+//     should be 0x00, 0x3f, 0x00<<4 | 0x00.
+var sosHeaderY = []byte{
+	0xff, 0xda, 0x00, 0x08, 0x01, 0x01, 0x00, 0x00, 0x3f, 0x00,
+}
+
+// sosHeaderYCbCr is the SOS marker "\xff\xda" followed by 12 bytes:
+//   - the marker length "\x00\x0c",
+//   - the number of components "\x03",
+//   - component 1 uses DC table 0 and AC table 0 "\x01\x00",
+//   - component 2 uses DC table 1 and AC table 1 "\x02\x11",
+//   - component 3 uses DC table 1 and AC table 1 "\x03\x11",
+//   - the bytes "\x00\x3f\x00". Section B.2.3 of the spec says that for
+//     sequential DCTs, those bytes (8-bit Ss, 8-bit Se, 4-bit Ah, 4-bit Al)
+//     should be 0x00, 0x3f, 0x00<<4 | 0x00.
+var sosHeaderYCbCr = []byte{
+	0xff, 0xda, 0x00, 0x0c, 0x03, 0x01, 0x00, 0x02,
+	0x11, 0x03, 0x11, 0x00, 0x3f, 0x00,
+}
+
+// writeSOS writes the StartOfScan marker.
+func (e *encoder) writeSOS(m image.Image) {
+	switch m.(type) {
+	case *image.Gray:
+		e.write(sosHeaderY)
+	default:
+		e.write(sosHeaderYCbCr)
+	}
+	var (
+		// Scratch buffers to hold the YCbCr values.
+		// The blocks are in natural (not zig-zag) order.
+		b      block
+		cb, cr [4]block
+		// DC components are delta-encoded.
+		prevDCY, prevDCCb, prevDCCr int32
+	)
+	bounds := m.Bounds()
+	switch m := m.(type) {
+	// TODO(wathiede): switch on m.ColorModel() instead of type.
+	case *image.Gray:
+		for y := bounds.Min.Y; y < bounds.Max.Y; y += 8 {
+			for x := bounds.Min.X; x < bounds.Max.X; x += 8 {
+				p := image.Pt(x, y)
+				grayToY(m, p, &b)
+				prevDCY = e.writeBlock(&b, 0, prevDCY)
+			}
+		}
+	default:
+		rgba, _ := m.(*image.RGBA)
+		ycbcr, _ := m.(*image.YCbCr)
+		for y := bounds.Min.Y; y < bounds.Max.Y; y += 16 {
+			for x := bounds.Min.X; x < bounds.Max.X; x += 16 {
+				for i := 0; i < 4; i++ {
+					xOff := (i & 1) * 8
+					yOff := (i & 2) * 4
+					p := image.Pt(x+xOff, y+yOff)
+					if rgba != nil {
+						rgbaToYCbCr(rgba, p, &b, &cb[i], &cr[i])
+					} else if ycbcr != nil {
+						yCbCrToYCbCr(ycbcr, p, &b, &cb[i], &cr[i])
+					} else {
+						toYCbCr(m, p, &b, &cb[i], &cr[i])
+					}
+					prevDCY = e.writeBlock(&b, 0, prevDCY)
+				}
+				scale(&b, &cb)
+				prevDCCb = e.writeBlock(&b, 1, prevDCCb)
+				scale(&b, &cr)
+				prevDCCr = e.writeBlock(&b, 1, prevDCCr)
+			}
+		}
+	}
+	// Pad the last byte with 1's.
+	e.emit(0x7f, 7)
+}
+
+// DefaultQuality is the default quality encoding parameter.
+const DefaultQuality = 75
+
+// Options are the encoding parameters.
+// Quality ranges from 1 to 100 inclusive, higher is better.
+type Options struct {
+	Quality int
+}
+
+// Encode writes the Image m to w in JPEG 4:2:0 baseline format with the given
+// options. Default parameters are used if a nil *Options is passed.
+func Encode(w io.Writer, m image.Image, o *Options) error {
+	b := m.Bounds()
+	if b.Dx() >= 1<<16 || b.Dy() >= 1<<16 {
+		return errors.New("jpeg: image is too large to encode")
+	}
+	var e encoder
+	if ww, ok := w.(writer); ok {
+		e.w = ww
+	} else {
+		e.w = bufio.NewWriter(w)
+	}
+	// Clip quality to [1, 100].
+	quality := DefaultQuality
+	if o != nil {
+		quality = o.Quality
+		if quality < 1 {
+			quality = 1
+		} else if quality > 100 {
+			quality = 100
+		}
+	}
+	// Convert from a quality rating to a scaling factor.
+	var scale int
+	if quality < 50 {
+		scale = 5000 / quality
+	} else {
+		scale = 200 - quality*2
+	}
+	// Initialize the quantization tables.
+	for i := range e.quant {
+		for j := range e.quant[i] {
+			x := int(unscaledQuant[i][j])
+			x = (x*scale + 50) / 100
+			if x < 1 {
+				x = 1
+			} else if x > 255 {
+				x = 255
+			}
+			e.quant[i][j] = uint8(x)
+		}
+	}
+	// Compute number of components based on input image type.
+	nComponent := 3
+	switch m.(type) {
+	// TODO(wathiede): switch on m.ColorModel() instead of type.
+	case *image.Gray:
+		nComponent = 1
+	}
+	// Write the Start Of Image marker.
+	e.buf[0] = 0xff
+	e.buf[1] = 0xd8
+	e.write(e.buf[:2])
+	// Write the quantization tables.
+	e.writeDQT()
+	// Write the image dimensions.
+	e.writeSOF0(b.Size(), nComponent)
+	// Write the Huffman tables.
+	e.writeDHT(nComponent)
+	// Write the image data.
+	e.writeSOS(m)
+	// Write the End Of Image marker.
+	e.buf[0] = 0xff
+	e.buf[1] = 0xd9
+	e.write(e.buf[:2])
+	e.flush()
+	return e.err
+}
diff --git a/src/image/jpeg/writer_test.go b/src/image/jpeg/writer_test.go
new file mode 100644
index 0000000..3414770
--- /dev/null
+++ b/src/image/jpeg/writer_test.go
@@ -0,0 +1,289 @@
+// Copyright 2011 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 jpeg
+
+import (
+	"bytes"
+	"fmt"
+	"image"
+	"image/color"
+	"image/png"
+	"io"
+	"math/rand"
+	"os"
+	"strings"
+	"testing"
+)
+
+// zigzag maps from the natural ordering to the zig-zag ordering. For example,
+// zigzag[0*8 + 3] is the zig-zag sequence number of the element in the fourth
+// column and first row.
+var zigzag = [blockSize]int{
+	0, 1, 5, 6, 14, 15, 27, 28,
+	2, 4, 7, 13, 16, 26, 29, 42,
+	3, 8, 12, 17, 25, 30, 41, 43,
+	9, 11, 18, 24, 31, 40, 44, 53,
+	10, 19, 23, 32, 39, 45, 52, 54,
+	20, 22, 33, 38, 46, 51, 55, 60,
+	21, 34, 37, 47, 50, 56, 59, 61,
+	35, 36, 48, 49, 57, 58, 62, 63,
+}
+
+func TestZigUnzig(t *testing.T) {
+	for i := 0; i < blockSize; i++ {
+		if unzig[zigzag[i]] != i {
+			t.Errorf("unzig[zigzag[%d]] == %d", i, unzig[zigzag[i]])
+		}
+		if zigzag[unzig[i]] != i {
+			t.Errorf("zigzag[unzig[%d]] == %d", i, zigzag[unzig[i]])
+		}
+	}
+}
+
+// unscaledQuantInNaturalOrder are the unscaled quantization tables in
+// natural (not zig-zag) order, as specified in section K.1.
+var unscaledQuantInNaturalOrder = [nQuantIndex][blockSize]byte{
+	// Luminance.
+	{
+		16, 11, 10, 16, 24, 40, 51, 61,
+		12, 12, 14, 19, 26, 58, 60, 55,
+		14, 13, 16, 24, 40, 57, 69, 56,
+		14, 17, 22, 29, 51, 87, 80, 62,
+		18, 22, 37, 56, 68, 109, 103, 77,
+		24, 35, 55, 64, 81, 104, 113, 92,
+		49, 64, 78, 87, 103, 121, 120, 101,
+		72, 92, 95, 98, 112, 100, 103, 99,
+	},
+	// Chrominance.
+	{
+		17, 18, 24, 47, 99, 99, 99, 99,
+		18, 21, 26, 66, 99, 99, 99, 99,
+		24, 26, 56, 99, 99, 99, 99, 99,
+		47, 66, 99, 99, 99, 99, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+		99, 99, 99, 99, 99, 99, 99, 99,
+	},
+}
+
+func TestUnscaledQuant(t *testing.T) {
+	bad := false
+	for i := quantIndex(0); i < nQuantIndex; i++ {
+		for zig := 0; zig < blockSize; zig++ {
+			got := unscaledQuant[i][zig]
+			want := unscaledQuantInNaturalOrder[i][unzig[zig]]
+			if got != want {
+				t.Errorf("i=%d, zig=%d: got %d, want %d", i, zig, got, want)
+				bad = true
+			}
+		}
+	}
+	if bad {
+		names := [nQuantIndex]string{"Luminance", "Chrominance"}
+		buf := &strings.Builder{}
+		for i, name := range names {
+			fmt.Fprintf(buf, "// %s.\n{\n", name)
+			for zig := 0; zig < blockSize; zig++ {
+				fmt.Fprintf(buf, "%d, ", unscaledQuantInNaturalOrder[i][unzig[zig]])
+				if zig%8 == 7 {
+					buf.WriteString("\n")
+				}
+			}
+			buf.WriteString("},\n")
+		}
+		t.Logf("expected unscaledQuant values:\n%s", buf.String())
+	}
+}
+
+var testCase = []struct {
+	filename  string
+	quality   int
+	tolerance int64
+}{
+	{"../testdata/video-001.png", 1, 24 << 8},
+	{"../testdata/video-001.png", 20, 12 << 8},
+	{"../testdata/video-001.png", 60, 8 << 8},
+	{"../testdata/video-001.png", 80, 6 << 8},
+	{"../testdata/video-001.png", 90, 4 << 8},
+	{"../testdata/video-001.png", 100, 2 << 8},
+}
+
+func delta(u0, u1 uint32) int64 {
+	d := int64(u0) - int64(u1)
+	if d < 0 {
+		return -d
+	}
+	return d
+}
+
+func readPng(filename string) (image.Image, error) {
+	f, err := os.Open(filename)
+	if err != nil {
+		return nil, err
+	}
+	defer f.Close()
+	return png.Decode(f)
+}
+
+func TestWriter(t *testing.T) {
+	for _, tc := range testCase {
+		// Read the image.
+		m0, err := readPng(tc.filename)
+		if err != nil {
+			t.Error(tc.filename, err)
+			continue
+		}
+		// Encode that image as JPEG.
+		var buf bytes.Buffer
+		err = Encode(&buf, m0, &Options{Quality: tc.quality})
+		if err != nil {
+			t.Error(tc.filename, err)
+			continue
+		}
+		// Decode that JPEG.
+		m1, err := Decode(&buf)
+		if err != nil {
+			t.Error(tc.filename, err)
+			continue
+		}
+		if m0.Bounds() != m1.Bounds() {
+			t.Errorf("%s, bounds differ: %v and %v", tc.filename, m0.Bounds(), m1.Bounds())
+			continue
+		}
+		// Compare the average delta to the tolerance level.
+		if averageDelta(m0, m1) > tc.tolerance {
+			t.Errorf("%s, quality=%d: average delta is too high", tc.filename, tc.quality)
+			continue
+		}
+	}
+}
+
+// TestWriteGrayscale tests that a grayscale images survives a round-trip
+// through encode/decode cycle.
+func TestWriteGrayscale(t *testing.T) {
+	m0 := image.NewGray(image.Rect(0, 0, 32, 32))
+	for i := range m0.Pix {
+		m0.Pix[i] = uint8(i)
+	}
+	var buf bytes.Buffer
+	if err := Encode(&buf, m0, nil); err != nil {
+		t.Fatal(err)
+	}
+	m1, err := Decode(&buf)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if m0.Bounds() != m1.Bounds() {
+		t.Fatalf("bounds differ: %v and %v", m0.Bounds(), m1.Bounds())
+	}
+	if _, ok := m1.(*image.Gray); !ok {
+		t.Errorf("got %T, want *image.Gray", m1)
+	}
+	// Compare the average delta to the tolerance level.
+	want := int64(2 << 8)
+	if got := averageDelta(m0, m1); got > want {
+		t.Errorf("average delta too high; got %d, want <= %d", got, want)
+	}
+}
+
+// averageDelta returns the average delta in RGB space. The two images must
+// have the same bounds.
+func averageDelta(m0, m1 image.Image) int64 {
+	b := m0.Bounds()
+	var sum, n int64
+	for y := b.Min.Y; y < b.Max.Y; y++ {
+		for x := b.Min.X; x < b.Max.X; x++ {
+			c0 := m0.At(x, y)
+			c1 := m1.At(x, y)
+			r0, g0, b0, _ := c0.RGBA()
+			r1, g1, b1, _ := c1.RGBA()
+			sum += delta(r0, r1)
+			sum += delta(g0, g1)
+			sum += delta(b0, b1)
+			n += 3
+		}
+	}
+	return sum / n
+}
+
+func TestEncodeYCbCr(t *testing.T) {
+	bo := image.Rect(0, 0, 640, 480)
+	imgRGBA := image.NewRGBA(bo)
+	// Must use 444 subsampling to avoid lossy RGBA to YCbCr conversion.
+	imgYCbCr := image.NewYCbCr(bo, image.YCbCrSubsampleRatio444)
+	rnd := rand.New(rand.NewSource(123))
+	// Create identical rgba and ycbcr images.
+	for y := bo.Min.Y; y < bo.Max.Y; y++ {
+		for x := bo.Min.X; x < bo.Max.X; x++ {
+			col := color.RGBA{
+				uint8(rnd.Intn(256)),
+				uint8(rnd.Intn(256)),
+				uint8(rnd.Intn(256)),
+				255,
+			}
+			imgRGBA.SetRGBA(x, y, col)
+			yo := imgYCbCr.YOffset(x, y)
+			co := imgYCbCr.COffset(x, y)
+			cy, ccr, ccb := color.RGBToYCbCr(col.R, col.G, col.B)
+			imgYCbCr.Y[yo] = cy
+			imgYCbCr.Cb[co] = ccr
+			imgYCbCr.Cr[co] = ccb
+		}
+	}
+
+	// Now check that both images are identical after an encode.
+	var bufRGBA, bufYCbCr bytes.Buffer
+	Encode(&bufRGBA, imgRGBA, nil)
+	Encode(&bufYCbCr, imgYCbCr, nil)
+	if !bytes.Equal(bufRGBA.Bytes(), bufYCbCr.Bytes()) {
+		t.Errorf("RGBA and YCbCr encoded bytes differ")
+	}
+}
+
+func BenchmarkEncodeRGBA(b *testing.B) {
+	img := image.NewRGBA(image.Rect(0, 0, 640, 480))
+	bo := img.Bounds()
+	rnd := rand.New(rand.NewSource(123))
+	for y := bo.Min.Y; y < bo.Max.Y; y++ {
+		for x := bo.Min.X; x < bo.Max.X; x++ {
+			img.SetRGBA(x, y, color.RGBA{
+				uint8(rnd.Intn(256)),
+				uint8(rnd.Intn(256)),
+				uint8(rnd.Intn(256)),
+				255,
+			})
+		}
+	}
+	b.SetBytes(640 * 480 * 4)
+	b.ReportAllocs()
+	b.ResetTimer()
+	options := &Options{Quality: 90}
+	for i := 0; i < b.N; i++ {
+		Encode(io.Discard, img, options)
+	}
+}
+
+func BenchmarkEncodeYCbCr(b *testing.B) {
+	img := image.NewYCbCr(image.Rect(0, 0, 640, 480), image.YCbCrSubsampleRatio420)
+	bo := img.Bounds()
+	rnd := rand.New(rand.NewSource(123))
+	for y := bo.Min.Y; y < bo.Max.Y; y++ {
+		for x := bo.Min.X; x < bo.Max.X; x++ {
+			cy := img.YOffset(x, y)
+			ci := img.COffset(x, y)
+			img.Y[cy] = uint8(rnd.Intn(256))
+			img.Cb[ci] = uint8(rnd.Intn(256))
+			img.Cr[ci] = uint8(rnd.Intn(256))
+		}
+	}
+	b.SetBytes(640 * 480 * 3)
+	b.ReportAllocs()
+	b.ResetTimer()
+	options := &Options{Quality: 90}
+	for i := 0; i < b.N; i++ {
+		Encode(io.Discard, img, options)
+	}
+}