一、前言
与本科期间导师给的课题相关，暑假期间学了一些关于DCT量化的内容。上网查阅了之后发现已经有博主使用MATLAB实现了这个功能。但是少有人用Python实现，于是我就动手写了一下。
实现起来并不复杂，主要的难点就是如何对图像进行均分快。由于二维DCT变换的性质，分割图像可以大幅提高效率。对图像进行DCT变换的主要思路就是将图像分割为多个小块，再对每一小块进行DCT变换。下面笔者先简单介绍一下DCT变换。(ps: 这篇文章在我在CSDN上已经发布过了，ID也是Icathia7，具体可以看链接：基于DCT的图像量化-Python简单实现)。
二、DCT变换
图像中的低频信号包含了图像主要信息(空间、大致轮廓、颜色等)，而高频信号包含图像的细节信息(细节轮廓，噪声)等。而人眼对细节信息，也就是高频信号并不敏感。对图像信号进行DCT变换可以将大部分能量聚集于低频信号，而减少高频信号。这方便了后续进行量化的操作。
我们常用的DCT变换的公式如下：
其中
DCT变换本质上是DFT变换在输入实偶函数下的一种特殊情况。即其数学性质与DFT相同。但DCT变换在做周期延拓之前增加了镜像延拓。这有效避免了延拓后的信号出现跳变。而时域上的跳变在频域上对应着高频分量。DCT变换有效避免了高频信号带来的干扰，相比DFT有更好的能量聚集。下图很好的说明了这一点：
此外，由于DCT变换是在实信号上进行的变换，相比较DFT的即有实数部分也有虚数部分，DCT的效率比DFT要高。因此DCT就更适合于图像、声音和视频等可能需要实时处理的信号。
因为图像属于二维信号，因此要做二维DCT变换。下面为二维DCT变换公式：
其中
不难看出二维DCT变换实际上是对图像所有行做DCT变换，再对所有列做DCT变换。直觉上来看再对所有列做变换似乎是多余的。但如果只做行或列变换的话，那么行与行(或列与列)之间的连接处就可能出现信号的跳变，从而引入高频分量，这不是我们所希望的。
因为在计算机上进行矩阵运算要比遍历乘法求和运算快得多，而二维DCT变换可以使用矩阵进行运算，下面为矩阵形式的DCT变换：
令变换矩阵A为
那么二维DCT变换可写作
完成DCT变换后就可以进行量化了，本文使用遮罩(mask)进行量化操作。
三、量化
图像经过DCT变换后，主要的低频能量都聚集于左上角，少量的高频能量于右下角。如下图所示。左上角为低频信号所在位置，右下角为高频信号所在位置。可以看到左上角颜色更亮，而右下角颜色更暗。亮色说明系数大，能量多，反之说明系数小，能量少。
这里量化的方式为保留左上角的低频信号，过滤右下角的高频信号。本文通过遮罩进行逐元素相乘来实现过滤。如下图所示。通过乘法，将遮罩中1的位置数据保留，0的部分删除。即可实现量化。或者可以通过系数的大小进行筛选，系数大于某个值就保留，否则删除，这里不多赘述，大家可以自行实现。
四、Python实现
以下为代码
# -*- coding: utf-8 -*-
"""
Created on Mon Sep 4 19:08:36 2023
@author: Icathia7
"""
import cv2
import numpy as np
import matplotlib.pyplot as plt
#%%
def DCT_quantify(data, block_size=8, set_mask=None):
# 获取图像大小
h, w = data.shape
# 获取遮罩
if type(set_mask)==np.ndarray:
mask = set_mask
else:
# 若不设置遮罩，默认为全1矩阵
mask = np.ones([block_size, block_size])
# 填充图像以便分割
if (h_to_pad:=(h % block_size)) != 0:
data = np.pad(data, ((0, block_size - h_to_pad), (0, 0)))
if (w_to_pad:=(w % block_size)) != 0:
data = np.pad(data, ((0, 0), (0, block_size - w_to_pad)))
new_h, new_w = data.shape
# 计算分割份数
v_slices_num = new_h // block_size
h_slices_num = new_w // block_size
# 按行分割
hori_data = np.vsplit(data, indices_or_sections=v_slices_num)
for i, row in enumerate(hori_data):
# 按列分割
vert_data = np.hsplit(row, indices_or_sections=h_slices_num)
first_v_block = cv2.dct(vert_data[0].astype(np.float32))
# 遮罩
first_v_block = np.multiply(first_v_block, mask)
first_v_iblock = cv2.idct(first_v_block)
dct_block_rows = first_v_block
idct_block_rows = first_v_iblock
for j, block in enumerate(vert_data[1:]):
# 进行DCT变换
single_block = cv2.dct(block.astype(np.float32))
# 遮罩
single_block = np.multiply(single_block, mask)
# 反变换
single_iblock = cv2.idct(single_block)
# 横向合并块成行
dct_block_rows = np.hstack([dct_block_rows, single_block])
idct_block_rows = np.hstack([idct_block_rows, single_iblock])
if i == 0:
dct_img = dct_block_rows
idct_img = idct_block_rows
else:
# 纵向合并行
dct_img = np.vstack([dct_img, dct_block_rows])
idct_img = np.vstack([idct_img, idct_block_rows])
return dct_img, idct_img
此函数有三个参数data, block_size, set_mask。其中data为图像数据；block_size为每一分块的大小，单位为像素，默认为8；set_mask为遮罩，默认为None，此时置遮罩为一全1矩阵，相当于未使用遮罩。需要注意的是，对于图像长宽无法被整除的情况，函数里做了一个判断，对无法整除的图像进行填充0的处理，这使得输出图像的长和宽可能会比输入图像大几个像素。
接下来进行测试：
测试代码，这里先不使用遮罩。
以下为代码
if __name__ == "__main__":
# 读取图像
img = plt.imread("output_png.png")
R = img[:, :, 0]
G = img[:, :, 1]
B = img[:, :, 2]
plt.figure(1)
plt.imshow(img)
#%%
# 自定义遮罩
mask = np.array([[1, 1, 0, 0, 0, 0, 0, 0],
[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, 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, 0]], dtype=np.float32)
# 分别对三个通道进行量化
B_dct, B_idct = DCT_quantify(B)
G_dct, G_idct = DCT_quantify(G)
R_dct, R_idct = DCT_quantify(R)
#%%
# 合并三通道图像
rebuild_img = np.stack((R_idct, G_idct, B_idct), axis=2)
plt.figure(2)
plt.imshow(rebuild_img)
测试结果：
原图：
还原后的图像：
可以看出，在不使用遮罩的情况下利用反DCT还原出来的图像与原图视觉上几乎没有差别。
接下来使用遮罩再测试一遍，代码如下：
以下为代码
# 自定义遮罩
mask = np.array([[1, 1, 0, 0, 0, 0, 0, 0],
[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, 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, 0]], dtype=np.float32)
# 分别对三个通道进行量化
B_dct, B_idct = DCT_quantify(B, set_mask=mask)
G_dct, G_idct = DCT_quantify(G, set_mask=mask)
R_dct, R_idct = DCT_quantify(R, set_mask=mask)
这里选择保留左上角3个系数，还原结果如下图：
压缩比约为: 1.7392
可以看出，在保留左上三个系数的情况下，图像大致内容都被还原出来了。要提升清晰度，我们可以增加保留的系数数量。
保留左上6个系数：
压缩比约为: 1.3274
保留左上10个系数：
压缩比约为1.1703
不难看出，我们还原出的图像清晰度随着系数的增加而增加，而压缩比随之降低。
接下来我们测试仅保留高频信息时还原出的图像：
保留右下49个数据：
保留右下54个系数：
保留右下58个系数：
可以看出，即使保留高频的系数比低频多得多，视觉上能还原的数据却少的多。 说明绝大部分能量都聚集在低频部分。
