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/others/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 explained every detail clearly, leading to several pitfalls during programming. Therefore, I plan to write an article that covers the details, incorporating Python code as much as possible. 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, and how the quantization and Huffman coding mentioned above are specifically implemented, please see the introduction in 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 subsampling 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, with one Y component corresponding to one Cb and Cr component (4:4:4). In contrast, 4:2:2 corresponds to two Y components for one Cb and Cr component, and 4:2:0 corresponds to four Y components for 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 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 conversion:
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, and the data from each block is finally 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, 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 upper left corner of the quantization table are smaller, while those in the lower right corner are larger. 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 quantization process code is as follows:
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 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 must also perform zigzag scanning on the quantization table accordingly. It must be stored in this form to allow 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 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 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 component are often large, and the DC components of adjacent 8*8 blocks are usually very close to each other. Therefore, using differential encoding can save space to a greater extent. Differential encoding means storing the difference between the current block's DC component 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 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 Components)#
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 an end marker; this bug took me a long time to find, see line 23 of the code below). After run-length encoding, the above list yields: (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, 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. This is because in the encoding of section 8, 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 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. Observing the pattern, it is not difficult to 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 stored binary data stream is 010110. For the data 6, we need to use the normalized 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, 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. 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, yielding a 4-bit representation of 0111. However, unlike the DC component, we need to perform normalized 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 that the normalized Huffman encoding for 0x04 results in 1011, then (0, -8) is stored as 10110111. We then perform the same operations for (0, -6) and so on, writing the resulting data streams into the file sequentially.
To give another example, (6, 1) has 1 stored as 1, which is 1 bit long. Therefore, for 0x61, we perform normalized Huffman encoding, assuming it results in 1111011, then (6, 1) is stored as 11110111. (15, 0) is stored as the normalized Huffman encoding value of 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's 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, it is necessary to check if the data being written is 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
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: normalized 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)
# Normalized 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
Normalized Huffman Encoding#
In this article, four normalized Huffman encoding tables are introduced, used for luminance DC components, chrominance DC components, luminance AC components, and chrominance AC components.
# Luminance DC normalized 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
As seen in the code above, std_huffman_DC0, etc., 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.
Normalized Huffman encoding has its own set of rules:
- The encoding for the first number of the minimum encoding length is 0;
- The encodings of the same length are continuous;
- The first number of the next encoding length (assuming j) depends on the last number of the previous encoding length (assuming i), i.e.,
a=(b+1)<<(j-i).
According to rule 1, 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 normalized Huffman encoding table to Huffman dictionary
# data: defined normalized 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 correct normalized 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 normalized 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 index of the dictionary in the write_num function.