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 very few explained every detail clearly, which led 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 bear with 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 the specific implementations of quantization and Huffman coding mentioned above, will be explained 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 will not significantly change 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 chroma. 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 contrast, 4:2:2 means two Y components correspond to one Cb and Cr component, while 4:2:0 means 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 represent the sampled pixels of the Cb and Cr components.
The formulas for color space conversion are:
All calculations are rounded to the nearest integer. In a 24-bit RGB BMP image, the ranges of the R, G, and B components are all 0-255. Through simple mathematical relationships, we can see that the ranges of the Y, Cb, and Cr components are also 0-255. In JPEG images, we typically 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 (the quantization 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, making it more suitable 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 quantization table are smaller in the top left corner and larger in the bottom right corner. An example of a quantization table is as follows (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
)
The code for the quantization process is:
now_block_qut = quantize(now_block_dct, 0) # Quantization for Y component
now_block_qut = quantize(now_block_dct, 2) # Quantization for Cb component
now_block_qut = quantize(now_block_dct, 1) # Quantization for Cr component
# 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 bottom right corner of the 8*8 block. To concentrate these zeros and allow run-length encoding to produce less data, we will 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 accordingly. It must be stored in this format for the image viewer to decode the image correctly (I spent a lot of debugging time on this detail at the beginning, as seen in the code for writing identifiers at the beginning).
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 usually very similar. Therefore, using differential encoding can save space to a greater extent. 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 correspondingly differentially encoded, meaning that 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 this 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, we store two numbers 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 used as an end marker (especially note that when the input data does not end with zero, two zeros are not needed as an end marker; 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 slightly reduced from the original 63. Of course, this is a special case; actual data will have more zeros, and the encoded size can be even smaller.
You may have noticed that the data (27, 1) is highlighted in red. This is because in the encoding of section 8, the first number stored is a 4-bit number, so its range is 0-15. Here, it clearly exceeds that range, 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 preparations, this section will truly introduce how the encoded DC and AC components are written into the file as data streams.
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 required bit length and the actual binary value to be stored according to the above format. It is not difficult to observe the pattern: the stored value of a positive number is its actual binary representation, and the bit length is also its actual bit length. The corresponding negative number has the same bit length and its binary value is the bitwise negation of the positive 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 Huffman encoding to store its binary data stream, which will be introduced in section 9. For now, let's assume that the binary data stream for 6 is stored as 100, then the DC value for a color component of this 8*8 block would be stored as 100010110.
After writing the binary data stream for the DC component into the file, we will then encode the AC values for this color component of the 8*8 block. The values obtained after run-length encoding 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, obtaining a 4-bit value and 0111. However, unlike the DC component, we need to perform 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 stored. Assuming the Huffman encoding for 0x04 results in 1011, then (0, -8) would be stored as 10110111. We then perform the same operations for (0, -6), and so on, writing the resulting data streams into the file.
To illustrate further, for (6, 1), the number 1 is stored as 1, which is 1 bit long. Therefore, for 0x61, we perform Huffman encoding, assuming it results in 1111011, then (6, 1) would be stored as 11110111. (15, 0) would only store the Huffman encoding value for 0xf0.
After writing the data for one color component (assumed to be Y), we then write the data for the Cb color component of this 88 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 0xff appears in data, append two 0x00
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 for DC)
# num: data to be written
# tbl: 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)
# Record the number of leading zeros (AC) and the bit length of num
s += tbl[tnum].str_code
# Store num in special format
s += tobin(num)
return s
Standard Huffman Encoding#
In this article, four Huffman encoding tables are introduced, which are used for luminance DC components, chrominance DC components, luminance AC components, and chrominance AC components.
# Luminance DC standard 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
)
...
# Calculate the 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. The std_huffman_DC0 describes the following image:
Now we only know the length of the encoded data for each original data, but we do not know what the actual data is.
Standard Huffman encoding has its own set of rules:
- The encoding for the smallest length is 0;
- The encodings of the same length are continuous;
- The first encoding of the next length (let's say j) depends on the last encoding (let's say i) of the previous length, i.e.,
a=(b+1)<<(j-i).
According to rule 1, we know that the encoding for 4 is 000. According to rule 2, the encoding for 5 is 001, for 3 is 010, for 2 is 011..., and for 0 is 110. According to rule 3, the encoding for 7 is 1110, 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 standard Huffman encoding table to a Huffman dictionary
# data: defined standard Huffman encoding table
def DHT2tbl(data):
numbers = data[0:16] # Number of encodings corresponding to 1~16 bit lengths
symbols = data[16:len(data)] # Original data
if(sum(numbers) != len(symbols)): # Check if it is a correct standard 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 standard 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 resulting Huffman dictionary is quite long and can be viewed in my GitHub project. Understanding the patterns within it can clarify how I derived the indices in the write_num function.