数字图像处理100问—38 JPEG 压缩——第二步：离散余弦变换+量化

提示：内容整理自：https://github.com/gzr2017/ImageProcessing100Wen
CV小白从0开始学数字图像处理

38 JPEG 压缩——第二步：离散余弦变换+量化

量化离散余弦变换系数并使用 离散余弦逆变换恢复。再比较变换前后图片的大小。

量化离散余弦变换系数是用于编码 JPEG 图像的技术。

量化即在对值在预定义的区间内舍入，其中floor、ceil、round等是类似的计算。

在 JPEG 图像中，根据下面所示的量化矩阵量化离散余弦变换系数。该量化矩阵取自 JPEG 软件开发联合会组织颁布的标准量化表。在量化中，将8x 8的系数除以（量化矩阵） Q 并四舍五入。之后然后再乘以 Q 。对于离散余弦逆变换，应使用所有系数。

Q = 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)), dtype=np.float32)

由于量化降低了图像的大小，因此可以看出数据量已经减少。

代码如下：

1.引入库

CV2计算机视觉库

import cv2
import numpy as np
import matplotlib.pyplot as plt

2.读入数据

img = cv2.imread("imori.jpg").astype(np.float32)
H, W, C = img.shape

3.灰度化

gray = 0.2126 * img[..., 2] + 0.7152 * img[..., 1] + 0.0722 * img[..., 0]

4.DCT

T = 8
K = 8
X = np.zeros((H, W), dtype=np.float64)
#indx = np.tile(np.arange(T), (T, 1))
#indy = np.arange(T).repeat(T).reshape(T, -1)
#dct = np.ones_like(indx, dtype=np.float32)
#dct[:, 0] /= np.sqrt(2)
#dct[0] /= np.sqrt(2)

Q = 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)), dtype=np.float32)

def w(x, y, u, v):
    cu = 1.
    cv = 1.
    if u == 0:
        cu /= np.sqrt(2)
    if v == 0:
        cv /= np.sqrt(2)
    theta = np.pi / (2 * T)
    return (( 2 * cu * cv / T) * np.cos((2*x+1)*u*theta) * np.cos((2*y+1)*v*theta))
    
for yi in range(0, H, T):
    for xi in range(0, W, T):
        for v in range(T):
            for u in range(T):
                for y in range(T):
                    for x in range(T):
                        X[v+yi, u+xi] += gray[y+yi, x+xi] * w(x,y,u,v)
        X[yi:yi+T, xi:xi+T] = np.round(X[yi:yi+T, xi:xi+T] / Q) * Q
                
                #_x = indx + xi * T
                #_y = indy + yi * T
                #_u = u + xi * T
                #_v = v + yi * T
                #X[_v, _u] = np.sum(C * gray[_y, _x] * np.cos((2*indx+1)*u*np.pi/(2*T)) * np.cos((2*indy+1)*v*np.pi/(2*T)))

5.IDCT

out = np.zeros((H, W), dtype=np.float64)

for yi in range(0, H, T):
    for xi in range(0, W, T):
        for y in range(T):
            for x in range(T):
                for v in range(K):
                    for u in range(K):
                        out[y+yi, x+xi] += X[v+yi, u+xi] * w(x,y,u,v)
                """
                _u = indx + xi * T
                _v = indy + yi * T
                _x = x + yi * T
                _y = y + xi * T
                out[_y, _x] = np.sum(C * X[_v, _u] * np.cos((2*x+1)*indx*np.pi/(2*T)) * np.cos((2*y+1)*indy*np.pi/(2*T))) * 4. / (T ** 2)
                """
out[out>255] = 255
out = np.floor(out).astype(np.uint8)

6.MSE

v_max = 255.
mse = np.sum(np.power(np.abs(gray.astype(np.float32) - out.astype(np.float32)), 2)) / (H * W)
psnr = 10 * np.log10(v_max ** 2 / mse)

print("PSNR >>", psnr)

bitrate = 1. * T * K ** 2 / (T ** 2)
print("bitrate >>", bitrate)

7.保存结果

cv2.imshow("result", out)
cv2.waitKey(0)
cv2.imwrite("out.jpg", out)

8.结果