This article is synchronized and updated to xLog by Mix Space
For the best browsing experience, it is recommended to visit the original link
https://www.do1e.cn/posts/codec/JPEG-detail
Preface#
This blog was originally published on 2021-08-22 on CSDN, and has been copied here with some formatting issues corrected.
Recently, I have been learning how to perform JPEG encoding. I found many articles online, but few explain every detail clearly, leading to many pitfalls during programming. Therefore, I plan to write an article that covers the details as much as possible, combined with Python code. The specific program can be referenced in my open-source project on GitHub.
Of course, my introduction and code are not very complete and may even contain some errors; they can only serve as a beginner's guide, so please forgive me.
Various Markers in JPEG Files#
Many articles introduce the markers in JPEG files. I have also uploaded a document that annotates an actual image (click to download) for reference.
All markers start with 0xff (hexadecimal 255), followed by the byte count representing this block of data and the data describing this block's information, as shown in the figure below:
# Write JPEG format decoding information
# filename: output file name
# h: image height
# w: image width
def write_head(filename, h, w):
# Open file in binary write mode (overwrite)
fp = open(filename, "wb")
# SOI
fp.write(pack(">H", 0xffd8))
# APP0
fp.write(pack(">H", 0xffe0))
fp.write(pack(">H", 16)) # APP0 byte count
fp.write(pack(">L", 0x4a464946)) # JFIF
fp.write(pack(">B", 0)) # 0
fp.write(pack(">H", 0x0101)) # Version: 1.1
fp.write(pack(">B", 0x01)) # Pixel density unit: pixels/inch
fp.write(pack(">L", 0x00480048)) # XY direction pixel density
fp.write(pack(">H", 0x0000)) # No thumbnail information
# DQT_0
fp.write(pack(">H", 0xffdb))
fp.write(pack(">H", 64+3)) # Quantization table byte count
fp.write(pack(">B", 0x00)) # Quantization table precision: 8bit(0) Quantization table ID: 0
tbl = block2zz(std_luminance_quant_tbl)
for item in tbl:
fp.write(pack(">B", item)) # Quantization table 0 content
# DQT_1
fp.write(pack(">H", 0xffdb))
fp.write(pack(">H", 64+3)) # Quantization table byte count
fp.write(pack(">B", 0x01)) # Quantization table precision: 8bit(0) Quantization table ID: 1
tbl = block2zz(std_chrominance_quant_tbl)
for item in tbl:
fp.write(pack(">B", item)) # Quantization table 1 content
# SOF0
fp.write(pack(">H", 0xffc0))
fp.write(pack(">H", 17)) # Frame image information byte count
fp.write(pack(">B", 8)) # Precision: 8bit
fp.write(pack(">H", h)) # Image height
fp.write(pack(">H", w)) # Image width
fp.write(pack(">B", 3)) # Number of color components: 3(YCrCb)
fp.write(pack(">B", 1)) # Color component ID: 1
fp.write(pack(">H", 0x1100)) # Horizontal and vertical sampling factor: 1 Quantization table ID used: 0
fp.write(pack(">B", 2)) # Color component ID: 2
fp.write(pack(">H", 0x1101)) # Horizontal and vertical sampling factor: 1 Quantization table ID used: 1
fp.write(pack(">B", 3)) # Color component ID: 3
fp.write(pack(">H", 0x1101)) # Horizontal and vertical sampling factor: 1 Quantization table ID used: 1
# DHT_DC0
fp.write(pack(">H", 0xffc4))
fp.write(pack(">H", len(std_huffman_DC0)+3)) # Huffman table byte count
fp.write(pack(">B", 0x00)) # DC0
for item in std_huffman_DC0:
fp.write(pack(">B", item)) # Huffman table content
# DHT_AC0
fp.write(pack(">H", 0xffc4))
fp.write(pack(">H", len(std_huffman_AC0)+3)) # Huffman table byte count
fp.write(pack(">B", 0x10)) # AC0
for item in std_huffman_AC0:
fp.write(pack(">B", item)) # Huffman table content
# DHT_DC1
fp.write(pack(">H", 0xffc4))
fp.write(pack(">H", len(std_huffman_DC1)+3)) # Huffman table byte count
fp.write(pack(">B", 0x01)) # DC1
for item in std_huffman_DC1:
fp.write(pack(">B", item)) # Huffman table content
# DHT_AC1
fp.write(pack(">H", 0xffc4))
fp.write(pack(">H", len(std_huffman_AC1)+3)) # Huffman table byte count
fp.write(pack(">B", 0x11)) # AC1
for item in std_huffman_AC1:
fp.write(pack(">B", item)) # Huffman table content
# SOS
fp.write(pack(">H", 0xffda))
fp.write(pack(">H", 12)) # Scan start information byte count
fp.write(pack(">B", 3)) # Number of color components: 3
fp.write(pack(">H", 0x0100)) # Color component 1 DC, AC Huffman table ID
fp.write(pack(">H", 0x0211)) # Color component 2 DC, AC Huffman table ID
fp.write(pack(">H", 0x0311)) # Color component 3 DC, AC Huffman table ID
fp.write(pack(">B", 0x00))
fp.write(pack(">B", 0x3f))
fp.write(pack(">B", 0x00)) # Fixed value
fp.close()
At this point, we have only the image data part left to write, but how the image data part is encoded, as well as how the quantization and Huffman encoding mentioned above are specifically implemented, please see the introduction in the next section.
JPEG Encoding Process#
Since JPEG encoding requires the image to be divided into 8*8 blocks, the height and width of the image must be multiples of 8. Therefore, we can use image interpolation or subsampling methods to make slight changes to the image, transforming it into an image whose height and width are both multiples of 8. For an image with thousands or tens of thousands of pixels, this operation does not significantly affect the overall aspect ratio of the image.
# Resize the image, must be divisible into 8*8 blocks
if((h % 8 == 0) and (w % 8 == 0)):
nblock = int(h * w / 64)
else:
h = h // 8 * 8
w = w // 8 * 8
YCrCb = cv2.resize(YCrCb, [h, w], cv2.INTER_CUBIC)
nblock = int(h * w / 64)
Color Space Conversion#
JPEG images uniformly use the YCbCr color space because the human eye is more sensitive to brightness than to chromaticity. Therefore, we selectively increase the compression of the Cb and Cr components, which can ensure that the visual quality of the image remains unchanged while significantly reducing the image size. After transforming to the YCbCr space, we can sample the Cb and Cr color components to reduce their pixel count, achieving greater compression. Common sampling formats include 4:4:4, 4:2:2, and 4:2:0. This corresponds to the horizontal and vertical sampling factors in the SOF0 marker. For simplicity, all sampling factors in this article are set to 1, meaning no sampling is performed, and one Y component corresponds to one Cb and Cr component (4:4:4). In 4:2:2, two Y components correspond to one Cb and Cr component, while in 4:2:0, four Y components correspond to one Cb and Cr component. As shown in the figure below, each cell corresponds to a Y component, while the blue cells are the sampled pixels of the Cb and Cr components.
The formulas for color space conversion are:
The above calculations are all rounded to the nearest integer. In a 24-bit RGB BMP image, the ranges of R, G, and B components are all 0-255. Through simple mathematical relationships, we can see that the ranges of Y, Cb, and Cr components are also 0-255. In JPEG images, we usually need to subtract 128 from each component to allow for both positive and negative ranges.
In Python, we can use functions from the OpenCV library to perform color space transformations:
YCrCb = cv2.cvtColor(BGR, cv2.COLOR_BGR2YCrCb)
npdata = np.array(YCrCb, np.int16)
8*8 Block Division#
In JPEG encoding, each 8*8 block is processed in order from top to bottom and left to right. Finally, the data from each block is combined together. For each block's Y, Cb, and Cr color components, the same operations are performed in the order of Y, Cb, and Cr (the quantization tables and Huffman tables used may differ).
for i in range(0, h, 8):
for j in range(0, w, 8):
...
DCT Transformation#
DCT (Discrete Cosine Transform) converts spatial domain data into frequency domain data, allowing us to selectively reduce high-frequency component data in the frequency domain without significantly affecting the visual quality of the image. Compared to the discrete Fourier transform, the discrete cosine transform operates entirely in the real number domain, which is more advantageous for computer calculations. The formula for the discrete cosine transform is as follows:
where $C(u)=\begin{cases}\frac{1}{\sqrt{2}}&u=0\\1&u\neq0\end{cases}$. In JPEG, $M=N=8$.
Of course, we can also use functions from the OpenCV library:
now_block = npdata[i:i+8, j:j+8, 0] - 128 # Extract an 8*8 block and subtract 128 from the Y component
now_block = npdata[i:i+8, j:j+8, 2] - 128 # Extract an 8*8 block and subtract 128 from the Cb component
now_block = npdata[i:i+8, j:j+8, 1] - 128 # Extract an 8*8 block and subtract 128 from the Cr component
now_block_dct = cv2.dct(np.float32(now_block)) # DCT transformation
Quantization#
After the DCT transformation, the DC component is the first element of the 88 block, with low-frequency components concentrated in the top left corner and high-frequency components concentrated in the bottom right corner. To selectively remove high-frequency components, we can perform quantization, which essentially involves dividing each element in the 88 block by a fixed value. The elements in the upper left corner of the quantization table are smaller, while those in the lower right corner are larger. An example of a quantization table is shown below (the Y component and Cb Cr components use different quantization tables):
# Luminance quantization table
std_luminance_quant_tbl = np.array(
[
[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]
],
np.uint8
)
# Chrominance quantization table
std_chrominance_quant_tbl = np.array(
[
[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]
],
np.uint8
)
Quantization process code:
now_block_qut = quantize(now_block_dct, 0) # Y component quantization
now_block_qut = quantize(now_block_dct, 2) # Cb component quantization
now_block_qut = quantize(now_block_dct, 1) # Cr component quantization
# Quantization
# block: current 8*8 block data
# dim: dimension 0:Y 1:Cr 2:Cb
def quantize(block, dim):
if(dim == 0):
# Use luminance quantization table
qarr = std_luminance_quant_tbl
else:
# Use chrominance quantization table
qarr = std_chrominance_quant_tbl
return (block / qarr).round().astype(np.int16)
After quantization, many zeros appear in the lower right corner of the 8*8 block. To concentrate these zeros and allow run-length encoding to produce less data, we next perform zigzag scanning.
Zigzag Scanning#
Zigzag scanning refers to transforming the 8*8 block into a list of 64 items in the following order.
Ultimately, we obtain a list of length 64: (41, -8, -6, -5, 13, 11, -1, 1, 2, -2, -3, -5, 1, 1, -5, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0). The subsequent operations will use this list as an example.
It is important to note that when storing the quantization table, we also need to perform zigzag scanning on the quantization table to store it in this form so that the image viewer can decode the image correctly (I initially spent a lot of debugging time on this detail), as seen in the code for writing the identifiers at the beginning of this article.
now_block_zz = block2zz(now_block_qut) # Zigzag scanning
# Zigzag scanning
# block: current 8*8 block data
def block2zz(block):
re = np.empty(64, np.int16)
# Current position in block
pos = np.array([0, 0])
# Define four scanning directions
R = np.array([0, 1])
LD = np.array([1, -1])
D = np.array([1, 0])
RU = np.array([-1, 1])
for i in range(0, 64):
re[i] = block[pos[0], pos[1]]
if(((pos[0] == 0) or (pos[0] == 7)) and (pos[1] % 2 == 0)):
pos = pos + R
elif(((pos[1] == 0) or (pos[1] == 7)) and (pos[0] % 2 == 1)):
pos = pos + D
elif((pos[0] + pos[1]) % 2 == 0):
pos = pos + RU
else:
pos = pos + LD
return re
Differential Encoding (DC Component)#
The values of the DC components are often large, and the DC components of adjacent 8*8 blocks are often very similar. Therefore, using differential encoding can save more space. Differential encoding means storing the difference between the current block and the previous block's DC component, while the first block stores its own value. It is important to note that the Y, Cb, and Cr components are differentially encoded correspondingly, meaning each component is subtracted from its counterpart. The following section will introduce how to encode and store the DC component now_block_dc.
last_block_ydc = 0
last_block_cbdc = 0
last_block_crdc = 0
now_block_dc = now_block_zz[0] - last_block_ydc # Record the DC component in differential form
last_block_ydc = now_block_zz[0] # Record the current value
now_block_dc = now_block_zz[0] - last_block_cbdc
last_block_cbdc = now_block_zz[0]
now_block_dc = now_block_zz[0] - last_block_crdc
last_block_crdc = now_block_zz[0]
Run-Length Encoding of Zeros (AC Component)#
After zigzag scanning, many zeros are concentrated together. The AC component list is: (-8, -6, -5, 13, 11, -1, 1, 2, -2, -3, -5, 1, 1, -5, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0).
In run-length encoding of zeros, two numbers are stored each time: the second number is a non-zero number, and the first number indicates how many zeros precede this non-zero number. Finally, two zeros are added as an end identifier (especially note that when the input data does not end with zero, two zeros are not needed as an end identifier; this bug took me a long time to find, as seen in line 23 of the code). After run-length encoding, the above list becomes: (0, -8), (0, -6), (0, -5), (0, 13), (0, 11), (0, -1), (0, 1), (0, 2), (0, -2), (0, -3), (0, -5), (0, 1), (0, 1), (0, -5), (0, 1), (3, -1), (6, 1), (0, 1), (0, -1), (27, 1), (0, 0). This data length is 42, which is a slight reduction from the original 63. Of course, this is a special case; actual data will have more zeros, even all zeros, and the encoded size can be smaller.
You may have noticed that the data (27, 1) is highlighted in red because in the encoding of the eighth part, the first number is stored as a 4-bit number, so the range is 0-15. Here, it clearly exceeds that, so we need to split it into (15, 0) and (11, 1), where (15, 0) represents 16 zeros, and (11, 1) represents 11 zeros followed by a 1.
now_block_ac = RLE(now_block_zz[1:])
# Run-length encoding of zeros
# AClist: data to be encoded
def RLE(AClist: np.ndarray) -> np.ndarray:
re = []
cnt = 0
for i in range(0, 63):
if(AClist[i] == 0 and cnt != 15):
cnt += 1
else:
re.append(cnt)
re.append(AClist[i])
cnt = 0
# Remove all trailing [15 0]
while(re[-1] == 0):
re.pop()
re.pop()
if(len(re) == 0):
break
# Add two zeros as an end marker
if(AClist[-1] == 0):
re.extend([0, 0])
return np.array(re, np.int16)
Special Binary Encoding for JPEG#
After the above groundwork, this section will truly introduce how the encoded DC and AC components are written into the file as a data stream.
In JPEG encoding, the following binary encoding format is used:
Value Bit Length Actual Stored Value
0 0 None
-1,1 1 0,1
-3,-2,2,3 2 00,01,10,11
-7,-6,-5,-4,4,5,6,7 3 000,001,010,011,100,101,110,111
-15,..,-8,8,..,15 4 0000,..,0111,1000,..,1111
-31,..,-16,16,..,31 5 00000,..,01111,10000,..,11111
-63,..,-32,32,..,63 6 ...
-127,..,-64,64,..,127 7 ...
-255,..,-128,128,..,255 8 ...
-511,..,-256,256,..,511 9 ...
-1023,..,-512,512,..,1023 10 ...
-2047,..,-1024,1024,..,2047 11 ...
For a number to be stored, we need to determine the bit length and the actual binary value to be stored according to the above format. Observing the pattern, we can see that the stored value for positive numbers is their actual binary representation, and the bit length is also their actual bit length. The same applies to negative numbers, where the bit length is the same and the binary value is the bitwise negation of the number. Zero does not need to be stored.
# Special binary encoding format
# num: number to be encoded
def tobin(num):
s = ""
if(num > 0):
while(num != 0):
s += '0' if(num % 2 == 0) else '1'
num = int(num / 2)
s = s[::-1] # Reverse
elif(num < 0):
num = -num
while(num != 0):
s += '1' if(num % 2 == 0) else '0'
num = int(num / 2)
s = s[::-1]
return s
For the DC component, assuming the value after differential encoding is -41, we can determine that its bit length is 6, and the binary data stream to be stored is 010110. For the data 6, we need to use the canonical Huffman encoding to store its binary data stream. This part will be introduced in section 9. For now, let's assume that the binary data stream for 6 is stored as 100, so the DC value for a color component of this 8*8 block is stored as 100010110.
After writing the binary data stream for the DC component into the file, we then encode the AC values for this color component of the 8*8 block. After run-length encoding, the values are: (0, -8), (0, -6), (0, -5), (0, 13), (0, 11), (0, -1), (0, 1), (0, 2), (0, -2), (0, -3), (0, -5), (0, 1), (0, 1), (0, -5), (0, 1), (3, -1), (6, 1), (0, 1), (0, -1), (15, 0), (11, 1), (0, 0).
First, we store (0, -8). For the second number, we perform the same operation and obtain 4 bits and 0111. However, unlike the DC component, we need to perform canonical Huffman encoding on 0x04, where the high four bits represent the first number of (0, -8), and the fourth bit represents the bit length of the second number. Assuming the canonical Huffman encoding for 0x04 is stored as 1011, then (0, -8) is stored as 10110111. Next, we perform the same operations for (0, -6), and so on, writing the resulting data streams into the file.
To give another example, (6, 1) has 1 stored as 1, which is 1 bit long. Therefore, for 0x61, we perform canonical Huffman encoding, assuming it is 1111011, so (6, 1) is stored as 11110111. (15, 0) only stores the canonical Huffman encoding value for 0xf0.
After writing the data for one color component (assumed to be Y) of an 88 block, we then write the data for the Cb color component of this block, followed by the Cr component data. Using the same method, we write the data for each 88 block from left to right and top to bottom, and finally write the EOI marker (0xffd9) to indicate the end of the image.
Note: During the data writing process, we need to check if we are writing 0xff. To prevent marker conflicts, we need to append 0x00 after it.
s = write_num(s, -1, now_block_dc, DC0) # Write DC data according to encoding method
for l in range(0, len(now_block_ac), 2): # Write AC data
s = write_num(s, now_block_ac[l], now_block_ac[l+1], AC0)
while(len(s) >= 8): # Record data too long will cause memory overflow
num = int(s[0:8], 2) # Running speed slows down
fp.write(pack(">B", num))
if(num == 0xff): # To prevent marker conflicts
fp.write(pack(">B", 0)) # If data contains 0xff, append two 0x00 after it
s = s[8:len(s)]
# Write data according to encoding method
# s: binary data not yet written to file
# n: number of leading zeros (-1 represents DC)
# num: data to be written
# tbl: canonical Huffman encoding dictionary
def write_num(s, n, num, tbl):
bit = 0
tnum = num
while(tnum != 0):
bit += 1
tnum = int(tnum / 2)
if(n == -1): # DC
tnum = bit
if(tnum > 11):
print("Write DC data Error")
exit()
else: # AC
if((n > 15) or (bit > 11) or (((n != 0) and (n != 15)) and (bit == 0))):
print("Write AC data Error")
exit()
tnum = n * 10 + bit + (0 if(n != 15) else 1)
# Canonical Huffman encoding records the number of zeros (AC) and the bit length of num
s += tbl[tnum].str_code
# Special form of data storage for num
s += tobin(num)
return s
Canonical Huffman Encoding#
In this article, four encoding tables are introduced for canonical Huffman encoding, specifically for luminance DC components, chrominance DC components, luminance AC components, and chrominance AC components.
# Luminance DC canonical Huffman encoding table
std_huffman_DC0 = np.array(
[0, 0, 7, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
4, 5, 3, 2, 6, 1, 0, 7, 8, 9, 10, 11],
np.uint8
)
...
# Convert to Huffman dictionary
DC0 = DHT2tbl(std_huffman_DC0) # Luminance DC component
DC1 = DHT2tbl(std_huffman_DC1) # Chrominance DC component
AC0 = DHT2tbl(std_huffman_AC0) # Luminance AC component
AC1 = DHT2tbl(std_huffman_AC1) # Chrominance AC component
In the above code, std_huffman_DC0 and others are the actual values stored in the identifiers. The first 16 numbers (0, 0, 7, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0) represent how many numbers there are for each length from 1 to 16 bits, while the following 12 numbers are exactly the sum of the previous 16 numbers. What std_huffman_DC0 describes is shown in the figure below:
Now we only know the length of the encoded data for each original data, but we do not know what the actual values are.
Canonical Huffman encoding has its own set of rules:
- The encoding for the smallest length is 0;
- Encodings of the same length are continuous;
- The first number of the next encoding length (let's say j) depends on the last number of the previous encoding length (let's say i), which is
a=(b+1)<<(j-i).
Based on rule 1, we know that the encoding for 4 is 000. Based on rule 2, the encoding for 5 is 001, the encoding for 3 is 010, the encoding for 2 is 011..., and the encoding for 0 is 110. Based on rule 3, the encoding for 7 is 1110, and the encoding for 8 is 11110...
# Class to record the Huffman dictionary
# symbol: original data
# code: corresponding encoded data
# n_bit: number of bits in the encoding
# str_code: binary data of the encoding
class Sym_Code():
def __init__(self, symbol, code, n_bit):
self.symbol = symbol
self.code = code
str_code=''
mask = 1 << (n_bit - 1)
for i in range(0, n_bit):
if(mask & code):
str_code += '1'
else:
str_code += '0'
mask >>= 1
self.str_code = str_code
"""Define output format"""
def __str__(self):
return "0x{:0>2x} | {}".format(self.symbol, self.str_code)
"""Define sorting criteria"""
def __eq__(self, other):
return self.symbol == other.symbol
def __le__(self, other):
return self.symbol < other.symbol
def __gt__(self, other):
return self.symbol > other.symbol
# Convert the canonical Huffman encoding table to a Huffman dictionary
# data: defined canonical Huffman encoding table
def DHT2tbl(data):
numbers = data[0:16] # Number of encodings corresponding to lengths of 1-16 bits
symbols = data[16:len(data)] # Original data
if(sum(numbers) != len(symbols)): # Check if it is a valid canonical Huffman encoding table
print("Wrong DHT!")
exit()
code = 0
SC = [] # List to record the dictionary
for n_bit in range(1, 17):
# Calculate the dictionary according to the canonical Huffman encoding rules
for symbol in symbols[sum(numbers[0:n_bit-1]):sum(numbers[0:n_bit])]:
SC.append(Sym_Code(symbol, code, n_bit))
code += 1
code <<= 1
return sorted(SC)
The final Huffman dictionary obtained is quite long; you can check it in my GitHub project. Understanding the patterns within it can clarify how I derived the index of the dictionary in the write_num function.