Initial QSfera import

This commit is contained in:
Курнат Андрей
2026-06-07 10:20:04 +03:00
commit 2315f25754
16485 changed files with 4826827 additions and 0 deletions
+523
View File
@@ -0,0 +1,523 @@
// Copyright 2025 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
// Discrete Cosine Transformation (DCT) implementations using the algorithm from
// Christoph Loeffler, Adriaan Lightenberg, and George S. Mostchytz,
// "Practical Fast 1-D DCT Algorithms with 11 Multiplications," ICASSP 1989.
// https://ieeexplore.ieee.org/document/266596
//
// Since the paper is paywalled, the rest of this comment gives a summary.
//
// A 1-dimensional forward DCT (1D FDCT) takes as input 8 values x0..x7
// and transforms them in place into the result values.
//
// The mathematical definition of the N-point 1D FDCT is:
//
// X[k] = α_k Σ_n x[n] * cos (2n+1)*k*π/2N
//
// where α₀ = √2 and α_k = 1 for k > 0.
//
// For our purposes, N=8, so the angles end up being multiples of π/16.
// The most direct implementation of this definition would require 64 multiplications.
//
// Loeffler's paper presents a more efficient computation that requires only
// 11 multiplications and works in terms of three basic operations:
//
// - A "butterfly" x0, x1 = x0+x1, x0-x1.
// The inverse is x0, x1 = (x0+x1)/2, (x0-x1)/2.
//
// - A scaling of x0 by k: x0 *= k. The inverse is scaling by 1/k.
//
// - A rotation of x0, x1 by θ, defined as:
// x0, x1 = x0 cos θ + x1 sin θ, -x0 sin θ + x1 cos θ.
// The inverse is rotation by -θ.
//
// The algorithm proceeds in four stages:
//
// Stage 1:
// - butterfly x0, x7; x1, x6; x2, x5; x3, x4.
//
// Stage 2:
// - butterfly x0, x3; x1, x2
// - rotate x4, x7 by 3π/16
// - rotate x5, x6 by π/16.
//
// Stage 3:
// - butterfly x0, x1; x4, x6; x7, x5
// - rotate x2, x3 by 6π/16 and scale by √2.
//
// Stage 4:
// - butterfly x7, x4
// - scale x5, x6 by √2.
//
// Finally, the values are permuted. The permutation can be read as either:
// - x0, x4, x2, x6, x7, x3, x5, x1 = x0, x1, x2, x3, x4, x5, x6, x7 (paper's form)
// - x0, x1, x2, x3, x4, x5, x6, x7 = x0, x7, x2, x5, x1, x6, x3, x4 (sorted by LHS)
//
// The code below uses the second form to make it easier to merge adjacent stores.
// (Note that unlike in recursive FFT implementations, the permutation here is
// not always mapping indexes to their bit reversals.)
//
// As written above, the rotation requires four multiplications, but it can be
// reduced to three by refactoring (see [dctBox] below), and the scaling in
// stage 3 can be merged into the rotation constants, so the overall cost
// of a 1D FDCT is 11 multiplies.
//
// The 1D inverse DCT (IDCT) is the 1D FDCT run backward
// with all the basic operations inverted.
// dctBox implements a 3-multiply, 3-add rotation+scaling.
// Given x0, x1, k*cos θ, and k*sin θ, dctBox returns the
// rotated and scaled coordinates.
// (It is called dctBox because the rotate+scale operation
// is drawn as a box in Figures 1 and 2 in the paper.)
func dctBox(x0, x1, kcos, ksin int32) (y0, y1 int32) {
// y0 = x0*kcos + x1*ksin
// y1 = -x0*ksin + x1*kcos
ksum := kcos * (x0 + x1)
y0 = ksum + (ksin-kcos)*x1
y1 = ksum - (kcos+ksin)*x0
return y0, y1
}
// A block is an 8x8 input to a 2D DCT (either the FDCT or IDCT).
// The input is actually only 8x8 uint8 values, and the outputs are 8x8 int16,
// but it is convenient to use int32s for intermediate storage,
// so we define only a single block type of [8*8]int32.
//
// A 2D DCT is implemented as 1D DCTs over the rows and columns.
//
// dct_test.go defines a String method for nice printing in tests.
type block [blockSize]int32
const blockSize = 8 * 8
// Note on Numerical Precision
//
// The inputs to both the FDCT and IDCT are uint8 values stored in a block,
// and the outputs are int16s in the same block, but the overall operation
// uses int32 values as fixed-point intermediate values.
// In the code comments below, the notation "QN.M" refers to a
// signed value of 1+N+M significant bits, one of which is the sign bit,
// and M of which hold fractional (sub-integer) precision.
// For example, 255 as a Q8.0 value is stored as int32(255),
// while 255 as a Q8.1 value is stored as int32(510),
// and 255.5 as a Q8.1 value is int32(511).
// The notation UQN.M refers to an unsigned value of N+M significant bits.
// See https://en.wikipedia.org/wiki/Q_(number_format) for more.
//
// In general we only need to keep about 16 significant bits, but it is more
// efficient and somewhat more precise to let unnecessary fractional bits
// accumulate and shift them away in bulk rather than after every operation.
// As such, it is important to keep track of the number of fractional bits
// in each variable at different points in the code, to avoid mistakes like
// adding numbers with different fractional precisions, as well as to keep
// track of the total number of bits, to avoid overflow. A comment like:
//
// // x[123] now Q8.2.
//
// means that x1, x2, and x3 are all Q8.2 (11-bit) values.
// Keeping extra precision bits also reduces the size of the errors introduced
// by using right shift to approximate rounded division.
// Constants needed for the implementation.
// These are all 60-bit precision fixed-point constants.
// The function c(val, b) rounds the constant to b bits.
// c is simple enough that calls to it with constant args
// are inlined and constant-propagated down to an inline constant.
// Each constant is commented with its Ivy definition (see robpike.io/ivy),
// using this scaling helper function:
//
// op fix x = floor 0.5 + x * 2**60
const (
cos1 = 1130768441178740757 // fix cos 1*pi/16
sin1 = 224923827593068887 // fix sin 1*pi/16
cos3 = 958619196450722178 // fix cos 3*pi/16
sin3 = 640528868967736374 // fix sin 3*pi/16
sqrt2 = 1630477228166597777 // fix sqrt 2
sqrt2_cos6 = 623956622067911264 // fix (sqrt 2)*cos 6*pi/16
sqrt2_sin6 = 1506364539328854985 // fix (sqrt 2)*sin 6*pi/16
sqrt2inv = 815238614083298888 // fix 1/sqrt 2
sqrt2inv_cos6 = 311978311033955632 // fix (1/sqrt 2)*cos 6*pi/16
sqrt2inv_sin6 = 753182269664427492 // fix (1/sqrt 2)*sin 6*pi/16
)
func c(x uint64, bits int) int32 {
return int32((x + (1 << (59 - bits))) >> (60 - bits))
}
// fdct implements the forward DCT.
// Inputs are UQ8.0; outputs are Q13.0.
func fdct(b *block) {
fdctCols(b)
fdctRows(b)
}
// fdctCols applies the 1D DCT to the columns of b.
// Inputs are UQ8.0 in [0,255] but interpreted as [-128,127].
// Outputs are Q10.18.
func fdctCols(b *block) {
for i := range 8 {
x0 := b[0*8+i]
x1 := b[1*8+i]
x2 := b[2*8+i]
x3 := b[3*8+i]
x4 := b[4*8+i]
x5 := b[5*8+i]
x6 := b[6*8+i]
x7 := b[7*8+i]
// x[01234567] are UQ8.0 in [0,255].
// Stage 1: four butterflies.
// In general a butterfly of QN.M inputs produces Q(N+1).M outputs.
// A butterfly of UQN.M inputs produces a UQ(N+1).M sum and a QN.M difference.
x0, x7 = x0+x7, x0-x7
x1, x6 = x1+x6, x1-x6
x2, x5 = x2+x5, x2-x5
x3, x4 = x3+x4, x3-x4
// x[0123] now UQ9.0 in [0, 510].
// x[4567] now Q8.0 in [-255,255].
// Stage 2: two boxes and two butterflies.
// A box on QN.M inputs with B-bit constants
// produces Q(N+1).(M+B) outputs.
// (The +1 is from the addition.)
x4, x7 = dctBox(x4, x7, c(cos3, 18), c(sin3, 18))
x5, x6 = dctBox(x5, x6, c(cos1, 18), c(sin1, 18))
// x[47] now Q9.18 in [-354, 354].
// x[56] now Q9.18 in [-300, 300].
x0, x3 = x0+x3, x0-x3
x1, x2 = x1+x2, x1-x2
// x[01] now UQ10.0 in [0, 1020].
// x[23] now Q9.0 in [-510, 510].
// Stage 3: one box and three butterflies.
x2, x3 = dctBox(x2, x3, c(sqrt2_cos6, 18), c(sqrt2_sin6, 18))
// x[23] now Q10.18 in [-943, 943].
x0, x1 = x0+x1, x0-x1
// x0 now UQ11.0 in [0, 2040].
// x1 now Q10.0 in [-1020, 1020].
// Store x0, x1, x2, x3 to their permuted targets.
// The original +128 in every input value
// has cancelled out except in the "DC signal" x0.
// Subtracting 128*8 here is equivalent to subtracting 128
// from every input before we started, but cheaper.
// It also converts x0 from UQ11.18 to Q10.18.
b[0*8+i] = (x0 - 128*8) << 18
b[4*8+i] = x1 << 18
b[2*8+i] = x2
b[6*8+i] = x3
x4, x6 = x4+x6, x4-x6
x7, x5 = x7+x5, x7-x5
// x[4567] now Q10.18 in [-654, 654].
// Stage 4: two √2 scalings and one butterfly.
x5 = (x5 >> 12) * c(sqrt2, 12)
x6 = (x6 >> 12) * c(sqrt2, 12)
// x[56] still Q10.18 in [-925, 925] (= 654√2).
x7, x4 = x7+x4, x7-x4
// x[47] still Q10.18 in [-925, 925] (not Q11.18!).
// This is not obvious at all! See "Note on 925" below.
// Store x4 x5 x6 x7 to their permuted targets.
b[1*8+i] = x7
b[3*8+i] = x5
b[5*8+i] = x6
b[7*8+i] = x4
}
}
// fdctRows applies the 1D DCT to the rows of b.
// Inputs are Q10.18; outputs are Q13.0.
func fdctRows(b *block) {
for i := range 8 {
x := b[8*i : 8*i+8 : 8*i+8]
x0 := x[0]
x1 := x[1]
x2 := x[2]
x3 := x[3]
x4 := x[4]
x5 := x[5]
x6 := x[6]
x7 := x[7]
// x[01234567] are Q10.18 [-1020, 1020].
// Stage 1: four butterflies.
x0, x7 = x0+x7, x0-x7
x1, x6 = x1+x6, x1-x6
x2, x5 = x2+x5, x2-x5
x3, x4 = x3+x4, x3-x4
// x[01234567] now Q11.18 in [-2040, 2040].
// Stage 2: two boxes and two butterflies.
x4, x7 = dctBox(x4>>14, x7>>14, c(cos3, 14), c(sin3, 14))
x5, x6 = dctBox(x5>>14, x6>>14, c(cos1, 14), c(sin1, 14))
// x[47] now Q12.18 in [-2830, 2830].
// x[56] now Q12.18 in [-2400, 2400].
x0, x3 = x0+x3, x0-x3
x1, x2 = x1+x2, x1-x2
// x[01234567] now Q12.18 in [-4080, 4080].
// Stage 3: one box and three butterflies.
x2, x3 = dctBox(x2>>14, x3>>14, c(sqrt2_cos6, 14), c(sqrt2_sin6, 14))
// x[23] now Q13.18 in [-7539, 7539].
x0, x1 = x0+x1, x0-x1
// x[01] now Q13.18 in [-8160, 8160].
x4, x6 = x4+x6, x4-x6
x7, x5 = x7+x5, x7-x5
// x[4567] now Q13.18 in [-5230, 5230].
// Stage 4: two √2 scalings and one butterfly.
x5 = (x5 >> 14) * c(sqrt2, 14)
x6 = (x6 >> 14) * c(sqrt2, 14)
// x[56] still Q13.18 in [-7397, 7397] (= 5230√2).
x7, x4 = x7+x4, x7-x4
// x[47] still Q13.18 in [-7395, 7395] (= 2040*3.6246).
// See "Note on 925" below.
// Cut from Q13.18 to Q13.0.
x0 = (x0 + 1<<17) >> 18
x1 = (x1 + 1<<17) >> 18
x2 = (x2 + 1<<17) >> 18
x3 = (x3 + 1<<17) >> 18
x4 = (x4 + 1<<17) >> 18
x5 = (x5 + 1<<17) >> 18
x6 = (x6 + 1<<17) >> 18
x7 = (x7 + 1<<17) >> 18
// Note: Unlike in fdctCols, saved all stores for the end
// because they are adjacent memory locations and some systems
// can use multiword stores.
x[0] = x0
x[1] = x7
x[2] = x2
x[3] = x5
x[4] = x1
x[5] = x6
x[6] = x3
x[7] = x4
}
}
// "Note on 925", deferred from above to avoid interrupting code.
//
// In fdctCols, heading into stage 2, the values x4, x5, x6, x7 are in [-255, 255].
// Let's call those specific values b4, b5, b6, b7, and trace how x[4567] evolve:
//
// Stage 2:
//
// x4 = b4*cos3 + b7*sin3
// x7 = -b4*sin3 + b7*cos3
// x5 = b5*cos1 + b6*sin1
// x6 = -b5*sin1 + b6*cos1
//
// Stage 3:
//
// x4 = x4+x6 = b4*cos3 + b7*sin3 - b5*sin1 + b6*cos1
// x6 = x4-x6 = b4*cos3 + b7*sin3 + b5*sin1 - b6*cos1
// x7 = x7+x5 = -b4*sin3 + b7*cos3 + b5*cos1 + b6*sin1
// x5 = x7-x5 = -b4*sin3 + b7*cos3 - b5*cos1 - b6*sin1
//
// Stage 4:
//
// x7 = x7+x4 = -b4*sin3 + b7*cos3 + b5*cos1 + b6*sin1 + b4*cos3 + b7*sin3 - b5*sin1 + b6*cos1
// = b4*(cos3-sin3) + b5*(cos1-sin1) + b6*(cos1+sin1) + b7*(cos3+sin3)
// < 255*(0.2759 + 0.7857 + 1.1759 + 1.3871) = 255*3.6246 < 925.
//
// x4 = x7-x4 = -b4*sin3 + b7*cos3 + b5*cos1 + b6*sin1 - b4*cos3 - b7*sin3 + b5*sin1 - b6*cos1
// = -b4*(cos3+sin3) + b5*(cos1+sin1) + b6*(sin1-cos1) + b7*(cos3-sin3)
// < same 925.
//
// The fact that x5, x6 are also at most 925 is not a coincidence: we are computing
// the same kinds of numbers for all four, just with different paths to them.
//
// In fdctRows, the same analysis applies, but the initial values are
// in [-2040, 2040] instead of [-255, 255], so the bound is 2040*3.6246 < 7395.
// idct implements the inverse DCT.
// Inputs are UQ8.0; outputs are Q10.3.
func idct(b *block) {
// A 2D IDCT is a 1D IDCT on rows followed by columns.
idctRows(b)
idctCols(b)
}
// idctRows applies the 1D IDCT to the rows of b.
// Inputs are UQ8.0; outputs are Q9.20.
func idctRows(b *block) {
for i := range 8 {
x := b[8*i : 8*i+8 : 8*i+8]
x0 := x[0]
x7 := x[1]
x2 := x[2]
x5 := x[3]
x1 := x[4]
x6 := x[5]
x3 := x[6]
x4 := x[7]
// Run FDCT backward.
// Independent operations have been reordered somewhat
// to make precision tracking easier.
//
// Note that "x0, x1 = x0+x1, x0-x1" is now a reverse butterfly
// and carries with it an implicit divide by two: the extra bit
// is added to the precision, not the value size.
// x[01234567] are UQ8.0 in [0, 255].
// Stages 4, 3, 2: x0, x1, x2, x3.
x0 <<= 17
x1 <<= 17
// x0, x1 now UQ8.17.
x0, x1 = x0+x1, x0-x1
// x0 now UQ8.18 in [0, 255].
// x1 now Q7.18 in [-127½, 127½].
// Note: (1/sqrt 2)*((cos 6*pi/16)+(sin 6*pi/16)) < 0.924, so no new high bit.
x2, x3 = dctBox(x2, x3, c(sqrt2inv_cos6, 18), -c(sqrt2inv_sin6, 18))
// x[23] now Q8.18 in [-236, 236].
x1, x2 = x1+x2, x1-x2
x0, x3 = x0+x3, x0-x3
// x[0123] now Q8.19 in [-246, 246].
// Stages 4, 3, 2: x4, x5, x6, x7.
x4 <<= 7
x7 <<= 7
// x[47] now UQ8.7
x7, x4 = x7+x4, x7-x4
// x7 now UQ8.8 in [0, 255].
// x4 now Q7.8 in [-127½, 127½].
x6 = x6 * c(sqrt2inv, 8)
x5 = x5 * c(sqrt2inv, 8)
// x[56] now UQ8.8 in [0, 181].
// Note that 1/√2 has five 0s in its binary representation after
// the 8th bit, so this multipliy is actually producing 12 bits of precision.
x7, x5 = x7+x5, x7-x5
x4, x6 = x4+x6, x4-x6
// x[4567] now Q8.9 in [-218, 218].
x4, x7 = dctBox(x4>>2, x7>>2, c(cos3, 12), -c(sin3, 12))
x5, x6 = dctBox(x5>>2, x6>>2, c(cos1, 12), -c(sin1, 12))
// x[4567] now Q9.19 in [-303, 303].
// Stage 1.
x0, x7 = x0+x7, x0-x7
x1, x6 = x1+x6, x1-x6
x2, x5 = x2+x5, x2-x5
x3, x4 = x3+x4, x3-x4
// x[01234567] now Q9.20 in [-275, 275].
// Note: we don't need all 20 bits of "precision",
// but it is faster to let idctCols shift it away as part
// of other operations rather than downshift here.
x[0] = x0
x[1] = x1
x[2] = x2
x[3] = x3
x[4] = x4
x[5] = x5
x[6] = x6
x[7] = x7
}
}
// idctCols applies the 1D IDCT to the columns of b.
// Inputs are Q9.20.
// Outputs are Q10.3. That is, the result is the IDCT*8.
func idctCols(b *block) {
for i := range 8 {
x0 := b[0*8+i]
x7 := b[1*8+i]
x2 := b[2*8+i]
x5 := b[3*8+i]
x1 := b[4*8+i]
x6 := b[5*8+i]
x3 := b[6*8+i]
x4 := b[7*8+i]
// x[012345678] are Q9.20.
// Start by adding 0.5 to x0 (the incoming DC signal).
// The butterflies will add it to all the other values,
// and then the final shifts will round properly.
x0 += 1 << 19
// Stages 4, 3, 2: x0, x1, x2, x3.
x0, x1 = (x0+x1)>>2, (x0-x1)>>2
// x[01] now Q9.19.
// Note: (1/sqrt 2)*((cos 6*pi/16)+(sin 6*pi/16)) < 1, so no new high bit.
x2, x3 = dctBox(x2>>13, x3>>13, c(sqrt2inv_cos6, 12), -c(sqrt2inv_sin6, 12))
// x[0123] now Q9.19.
x1, x2 = x1+x2, x1-x2
x0, x3 = x0+x3, x0-x3
// x[0123] now Q9.20.
// Stages 4, 3, 2: x4, x5, x6, x7.
x7, x4 = x7+x4, x7-x4
// x[47] now Q9.21.
x5 = (x5 >> 13) * c(sqrt2inv, 14)
x6 = (x6 >> 13) * c(sqrt2inv, 14)
// x[56] now Q9.21.
x7, x5 = x7+x5, x7-x5
x4, x6 = x4+x6, x4-x6
// x[4567] now Q9.22.
x4, x7 = dctBox(x4>>14, x7>>14, c(cos3, 12), -c(sin3, 12))
x5, x6 = dctBox(x5>>14, x6>>14, c(cos1, 12), -c(sin1, 12))
// x[4567] now Q10.20.
x0, x7 = x0+x7, x0-x7
x1, x6 = x1+x6, x1-x6
x2, x5 = x2+x5, x2-x5
x3, x4 = x3+x4, x3-x4
// x[01234567] now Q10.21.
x0 >>= 18
x1 >>= 18
x2 >>= 18
x3 >>= 18
x4 >>= 18
x5 >>= 18
x6 >>= 18
x7 >>= 18
// x[01234567] now Q10.3.
b[0*8+i] = x0
b[1*8+i] = x1
b[2*8+i] = x2
b[3*8+i] = x3
b[4*8+i] = x4
b[5*8+i] = x5
b[6*8+i] = x6
b[7*8+i] = x7
}
}
+245
View File
@@ -0,0 +1,245 @@
// 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.ErrUnexpectedEOF {
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.
clear(h.lut[:])
var x, code uint32
for i := range uint32(lutSize) {
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
}
+814
View File
@@ -0,0 +1,814 @@
// 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"
"io"
"github.com/kovidgoyal/go-parallel"
"github.com/kovidgoyal/imaging/nrgb"
"github.com/kovidgoyal/imaging/nrgba"
)
// 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.
expand struct{ h, v int } // subsampleRatio for this component
}
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
flex bool // uses a non-standard subsampleRatio
force_flex bool // used for testing
maxH, maxV int
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 {
return nil
}
if err == io.EOF {
err = io.ErrUnexpectedEOF
}
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 {
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 := min(d.bytes.j-d.bytes.i, n)
d.bytes.i += m
n -= m
if n == 0 {
break
}
if err := d.fill(); err != nil {
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
}
d.maxH, d.maxV = max(d.maxH, h), max(d.maxV, v)
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:
if i == 0 && v == 4 {
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
}
if d.nComp == 3 {
for i := range 3 {
if d.maxH%d.comp[i].h != 0 || d.maxV%d.comp[i].v != 0 {
return errUnsupportedSubsamplingRatio
}
}
}
for i := range d.nComp {
d.comp[i].expand.h = d.maxH / d.comp[i].h
d.comp[i].expand.v = d.maxV / d.comp[i].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)
src := nrgba.NewNRGBAScanner(d.img3)
w, h := img.Bounds().Dx(), img.Bounds().Dy()
size := w * 4
if err := parallel.Run_in_parallel_over_range(0, func(start, limit int) {
for y := start; y < limit; y++ {
i := y * img.Stride
src.Scan(0, y, w, y+1, img.Pix[i:i+size])
}
}, 0, h); err != nil {
panic(err)
}
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 := nrgb.NewNRGB(bounds)
parallel.Run_in_parallel_over_range(0, func(start, limit int) {
for y := start; y < limit; 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+3*i+0] = d.img3.Y[yo+i]
img.Pix[po+3*i+1] = d.img3.Cb[co+i/cScale]
img.Pix[po+3*i+2] = d.img3.Cr[co+i/cScale]
}
}
}, bounds.Min.Y, bounds.Max.Y)
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 = nrgb.Model
}
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")
}
+596
View File
@@ -0,0 +1,596 @@
// 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
}
subsampleRatio := image.YCbCrSubsampleRatio444
if d.comp[1].h != d.comp[2].h || d.comp[1].v != d.comp[2].v || d.maxH != d.comp[0].h || d.maxV != d.comp[0].v {
d.flex = true
} else {
if d.force_flex {
d.flex = true
} else {
hRatio := d.maxH / d.comp[1].h
vRatio := d.maxV / d.comp[1].v
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:
d.flex = true
}
}
}
m := image.NewYCbCr(image.Rect(0, 0, 8*d.maxH*mxx, 8*d.maxV*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 := range nComp {
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 := range i {
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 := range nComp {
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 := range myy {
for mx := range mxx {
for i := range nComp {
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 {
// For well-formed input, the RST[0-7] restart marker follows
// immediately. For corrupt input, call findRST to try to
// resynchronize.
if err := d.readFull(d.tmp[:2]); err != nil {
return err
} else if d.tmp[0] != 0xff || d.tmp[1] != expectedRST {
if err := d.findRST(expectedRST); err != nil {
return err
}
}
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
}
func (d *decoder) storeFlexBlock(b *block, bx, by, compIndex int) {
h, v := d.comp[compIndex].expand.h, d.comp[compIndex].expand.v
dst, stride := []byte(nil), 0
bx, by = bx*h, by*v
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
}
for y := range 8 {
y8 := y * 8
yv := y * v
for x := range 8 {
c := b[y8+x]
var val uint8
if c < -128 {
val = 0
} else if c > 127 {
val = 255
} else {
val = uint8(c + 128)
}
xh := x * h
for yy := range v {
for xx := range h {
dst[(yv+yy)*stride+xh+xx] = val
}
}
}
}
}
// 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 := range blockSize {
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 {
if d.flex {
d.storeFlexBlock(b, bx, by, compIndex)
return nil
}
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 := range 8 {
y8 := y * 8
yStride := y * stride
for x := range 8 {
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
}
// findRST advances past the next RST restart marker that matches expectedRST.
// Other than I/O errors, it is also an error if we encounter an {0xFF, M}
// two-byte marker sequence where M is not 0x00, 0xFF or the expectedRST.
//
// This is similar to libjpeg's jdmarker.c's next_marker function.
// https://github.com/libjpeg-turbo/libjpeg-turbo/blob/2dfe6c0fe9e18671105e94f7cbf044d4a1d157e6/jdmarker.c#L892-L935
//
// Precondition: d.tmp[:2] holds the next two bytes of JPEG-encoded input
// (input in the d.readFull sense).
func (d *decoder) findRST(expectedRST uint8) error {
for {
// i is the index such that, at the bottom of the loop, we read 2-i
// bytes into d.tmp[i:2], maintaining the invariant that d.tmp[:2]
// holds the next two bytes of JPEG-encoded input. It is either 0 or 1,
// so that each iteration advances by 1 or 2 bytes (or returns).
i := 0
if d.tmp[0] == 0xff {
if d.tmp[1] == expectedRST {
return nil
} else if d.tmp[1] == 0xff {
i = 1
} else if d.tmp[1] != 0x00 {
// libjpeg's jdmarker.c's jpeg_resync_to_restart does something
// fancy here, treating RST markers within two (modulo 8) of
// expectedRST differently from RST markers that are 'more
// distant'. Until we see evidence that recovering from such
// cases is frequent enough to be worth the complexity, we take
// a simpler approach for now. Any marker that's not 0x00, 0xff
// or expectedRST is a fatal FormatError.
return FormatError("bad RST marker")
}
} else if d.tmp[1] == 0xff {
d.tmp[0] = 0xff
i = 1
}
if err := d.readFull(d.tmp[i:2]); err != nil {
return err
}
}
}