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 explain every detail clearly, which led to many pitfalls while 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 only have the image data part left to write, but how the image data part is encoded, as well as how the quantization and Huffman coding mentioned above are specifically implemented, please see the next section.
JPEG Encoding Process#
Since the JPEG encoding process requires the image to be divided into 8*8 blocks, this requires that both the height and width of the image be multiples of 8. Therefore, we can use image interpolation or sampling methods to make slight changes to the image, transforming it into an image with both height and width as 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 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 size of the image. After transforming to the YCbCr space, we can sample the Cb and Cr color components to reduce their number of points, thus 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, with one Y component corresponding to one Cb and Cr component (4:4:4). In 4:2:2, two Y components correspond to one Cb and Cr component, and 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 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 for subsequent data processing. 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 for processing, allowing us to selectively reduce the data of high-frequency components 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 favorable 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 Y component
now_block = npdata[i:i+8, j:j+8, 2] - 128 # Extract an 8*8 block and subtract 128 from Cb component
now_block = npdata[i:i+8, j:j+8, 1] - 128 # Extract an 8*8 block and subtract 128 from 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 set of quantization tables 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
)
The code for the quantization process:
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 will perform zigzag scanning next.
Zigzag Scanning#
Zigzag scanning refers to converting 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 must also perform zigzag scanning on the quantization table to store it in this format so that the image viewer can decode the image correctly (I spent a lot of debugging time on this detail at the beginning), as seen in the code that writes 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 the 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 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 corresponding previous value. In the following section, I 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 # Differentially record the DC component
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 added as an end marker (especially note that when the input data does not end with zero, two zeros are not needed as end markers; this bug took me a long time to find, see line 23 of the code below). 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 has a length of 42, which is a slight reduction compared to 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 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: AC 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 in 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, there are the following binary encoding formats:
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 for positive numbers is their actual binary representation, and the bit length is also their actual bit length. The corresponding negative numbers have the same bit length and their 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, suppose the value after differential encoding is -41. Following the above operations, we can find that its bit length is 6, and the stored binary data stream 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 stored binary data stream for 6 is 100, then 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. 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 to obtain a 4-bit representation of 0111, but 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 that the Huffman encoding for 0x04 results in 1011, then (0, -8) is stored as 10110111. We continue this process for (0, -6) and so on, writing the resulting data streams into the file.
To give another example, for (6, 1), where 1 is stored as 1 (1 bit), we perform Huffman encoding on 0x61, assuming it results in 1111011, then (6, 1) is stored as 11110111. (15, 0) only stores the Huffman encoding value for 0xf0.
After writing the data for one color component (let's assume Y) of an 88 block, we then write the data for the Cb color component, followed by the Cr component. 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 afterwards.
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 preceding zeros (-1 represents 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 zeros for Huffman encoding (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 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, etc., are the actual values stored in the identifiers, which can be seen in the code in the identifier introduction. 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 of 1-16 bits after encoding, followed by 12 numbers that exactly match 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 encoding of the next length (let's say j) depends on the last encoding (let's say i) of the previous length, meaning
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, the encoding for 3 is 010, the encoding for 2 is 011... and the encoding for 0 is 110. According to rule 3, the encoding for 7 is 1110, 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 standard Huffman encoding table to Huffman dictionary
# data: defined standard 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 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 for the dictionary in the write_num function.