# 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/vertical sampling factor: 1  Quantization table ID used: 0
    fp.write(pack(">B", 2))                # Color component ID: 2
    fp.write(pack(">H", 0x1101))        # Horizontal/vertical sampling factor: 1  Quantization table ID used: 1
    fp.write(pack(">B", 3))                # Color component ID: 3
    fp.write(pack(">H", 0x1101))        # Horizontal/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 used
    fp.write(pack(">H", 0x0211))        # Color component 2 DC, AC Huffman table ID used
    fp.write(pack(">H", 0x0311))        # Color component 3 DC, AC Huffman table ID used
    fp.write(pack(">B", 0x00))
    fp.write(pack(">B", 0x3f))
    fp.write(pack(">B", 0x00))            # Fixed value
    fp.close()

# Convert image size, 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)

YCrCb = cv2.cvtColor(BGR, cv2.COLOR_BGR2YCrCb)
npdata = np.array(YCrCb, np.int16)

for i in range(0, h, 8):
    for j in range(0, w, 8):
        ...

now_block = npdata[i:i+8, j:j+8, 0] - 128        # Extract an 8*8 block and subtract 128 Y component
now_block = npdata[i:i+8, j:j+8, 2] - 128        # Extract an 8*8 block and subtract 128 Cb component
now_block = npdata[i:i+8, j:j+8, 1] - 128        # Extract an 8*8 block and subtract 128 Cr component
now_block_dct = cv2.dct(np.float32(now_block))    # DCT transform

# 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
)

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)

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

last_block_ydc = 0
last_block_cbdc = 0
last_block_crdc = 0

now_block_dc = now_block_zz[0] - last_block_ydc # DC component recorded in differential form
last_block_ydc = now_block_zz[0]                    # Record 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]

now_block_ac = RLE(now_block_zz[1:])

# Run-length encoding of zeros
# AClist: AC data to encode
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 at the end as end marker
    if(AClist[-1] == 0):
        re.extend([0, 0])
    return np.array(re, np.int16)

             Value              Bit length         Actual saved 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                  ...

# Special binary encoding format
# num: number to encode
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

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):                            # Data too long will cause memory overflow
        num = int(s[0:8], 2)                    # Slower execution speed
        fp.write(pack(">B", num))
        if(num == 0xff):                        # To prevent marker conflict
            fp.write(pack(">B", 0))                # When 0xff appears in data, need to 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 means DC)
# num: data to write
# 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 count of zeros (AC) and the bit length of num
    s += tbl[tnum].str_code
    # Special format data storage for num
    s += tobin(num)
    return s

# 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
)
...
# Calculate 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

# Class for recording Huffman dictionary
# symbol: original data
# code: corresponding encoded data
# n_bit: number of binary bits in encoding
# str_code: binary data of 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 basis"""
    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 canonical Huffman encoding table to Huffman dictionary
# data: defined canonical Huffman encoding table
def DHT2tbl(data):
    numbers = data[0:16]                # Count of codes for 1~16 bit lengths
    symbols = data[16:len(data)]        # Original data
    if(sum(numbers) != len(symbols)):    # Check if it's a valid canonical Huffman encoding table
        print("Wrong DHT!")
        exit()
    code = 0
    SC = []                                # List for recording dictionary
    for n_bit in range(1, 17):
        # Calculate dictionary according to 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)