用Python实战DCT与DWT：图像压缩与去噪的高效实现

当你第一次看到DCT（离散余弦变换）和DWT（离散小波变换）这两个术语时，是否被复杂的数学公式吓退？其实，这些看似高深的技术早已渗透到我们日常使用的JPEG图像、视频编码甚至医疗影像处理中。本文不会让你死记硬背任何公式，而是通过Python代码带你亲手实现图像压缩和去噪的完整流程。我们将使用OpenCV和PyWavelets这两个强大的库，从环境搭建到效果对比，一步步拆解这两个变换的实战应用。

1. 环境准备与基础概念速成

在开始编码前，我们需要快速理解DCT和DWT的核心思想。DCT就像把图像分解成不同频率的余弦波组合，而DWT则像用不同大小的"显微镜"观察图像细节。这两种变换都能将图像从空间域转换到频域，让我们能够有针对性地处理不同频率的成分。

1.1 安装必要的Python库

打开你的终端或命令提示符，运行以下命令安装所需库：

pip install opencv-python numpy matplotlib pywavelets scikit-image

验证安装是否成功：

import cv2
import pywt
print(f"OpenCV版本: {cv2.__version__}")
print(f"PyWavelets版本: {pywt.__version__}")

1.2 准备测试图像

我们将使用一张包含丰富细节和平滑区域的测试图像。你可以使用自己的图片，或者从skimage库加载示例图像：

from skimage import data
import matplotlib.pyplot as plt

# 加载示例图像
image = data.camera()
plt.imshow(image, cmap='gray')
plt.title("原始图像")
plt.show()

提示：选择测试图像时，最好包含清晰的边缘、纹理区域和平滑区域，这样能更直观地观察变换效果。

2. DCT实战：JPEG压缩的核心技术

DCT是JPEG图像压缩的核心算法，它通过保留重要的低频信息，舍弃对人眼不敏感的高频细节来实现压缩。让我们一步步实现这个过程。

2.1 分块DCT变换与可视化

JPEG标准使用8×8的分块DCT变换。以下是实现代码：

import numpy as np

def apply_block_dct(image, block_size=8):
    # 将图像转换为float32类型
    img_float = np.float32(image) / 255.0
    # 获取图像尺寸
    h, w = img_float.shape
    # 初始化DCT系数矩阵
    dct_blocks = np.zeros_like(img_float)
    
    # 对每个8x8块应用DCT
    for i in range(0, h, block_size):
        for j in range(0, w, block_size):
            block = img_float[i:i+block_size, j:j+block_size]
            dct_block = cv2.dct(block)
            dct_blocks[i:i+block_size, j:j+block_size] = dct_block
    
    return dct_blocks

# 应用DCT变换
dct_coeff = apply_block_dct(image)

可视化DCT系数矩阵：

plt.figure(figsize=(12,6))
plt.subplot(121), plt.imshow(image, cmap='gray'), plt.title("原始图像")
plt.subplot(122), plt.imshow(np.log1p(np.abs(dct_coeff)), cmap='jet'), plt.title("DCT系数矩阵(对数尺度)")
plt.colorbar()
plt.show()

你会注意到DCT系数矩阵的左上角区域值较大，这对应图像的低频成分（整体亮度和大致轮廓），而右下角的高频系数值较小，对应图像的细节和边缘。

2.2 实现简单的JPEG压缩模拟

JPEG压缩的关键步骤是量化 - 即舍弃小的DCT系数。我们来实现一个简单的量化过程：

def quantize_dct(dct_coeff, keep_fraction=0.25):
    # 创建量化掩码
    mask = np.zeros_like(dct_coeff)
    block_size = 8
    h, w = dct_coeff.shape
    
    # 对每个8x8块应用相同的量化模式
    for i in range(0, h, block_size):
        for j in range(0, w, block_size):
            # 保留左上角的部分系数
            for x in range(block_size):
                for y in range(block_size):
                    if x + y < block_size * keep_fraction:
                        mask[i+x, j+y] = 1
    return dct_coeff * mask

# 只保留25%的DCT系数
quantized_dct = quantize_dct(dct_coeff, 0.25)

然后进行逆DCT变换重建图像：

def apply_block_idct(dct_blocks, block_size=8):
    h, w = dct_blocks.shape
    reconstructed = np.zeros_like(dct_blocks, dtype=np.float32)
    
    for i in range(0, h, block_size):
        for j in range(0, w, block_size):
            block = dct_blocks[i:i+block_size, j:j+block_size]
            idct_block = cv2.idct(block)
            reconstructed[i:i+block_size, j:j+block_size] = idct_block
    
    return np.clip(reconstructed * 255, 0, 255).astype(np.uint8)

# 重建图像
reconstructed_img = apply_block_idct(quantized_dct)

比较原始图像和重建图像：

plt.figure(figsize=(12,6))
plt.subplot(121), plt.imshow(image, cmap='gray'), plt.title("原始图像")
plt.subplot(122), plt.imshow(reconstructed_img, cmap='gray'), plt.title(f"压缩后图像(保留25%系数)")
plt.show()

2.3 DCT参数调优实战

DCT压缩的效果很大程度上取决于保留多少系数以及如何选择这些系数。让我们尝试不同的量化策略：

def zigzag_mask(block_size, keep_count):
    """生成Z字形扫描的掩码"""
    mask = np.zeros((block_size, block_size), dtype=bool)
    rows, cols = 0, 0
    direction = 1  # 1表示向上移动，-1表示向下移动
    
    for i in range(keep_count):
        mask[rows, cols] = True
        if direction == 1:
            if cols == block_size - 1:
                rows += 1
                direction = -1
            elif rows == 0:
                cols += 1
                direction = -1
            else:
                rows -= 1
                cols += 1
        else:
            if rows == block_size - 1:
                cols += 1
                direction = 1
            elif cols == 0:
                rows += 1
                direction = 1
            else:
                rows += 1
                cols -= 1
    return mask

def advanced_quantize(dct_coeff, keep_fraction=0.25):
    block_size = 8
    h, w = dct_coeff.shape
    mask = np.zeros_like(dct_coeff, dtype=bool)
    keep_per_block = int(block_size**2 * keep_fraction)
    
    for i in range(0, h, block_size):
        for j in range(0, w, block_size):
            block_mask = zigzag_mask(block_size, keep_per_block)
            mask[i:i+block_size, j:j+block_size] = block_mask
    
    return dct_coeff * mask

# 使用Z字形扫描保留系数
advanced_quantized = advanced_quantize(dct_coeff, 0.25)
advanced_reconstructed = apply_block_idct(advanced_quantized)

比较两种量化方法的效果：

plt.figure(figsize=(15,5))
plt.subplot(131), plt.imshow(image, cmap='gray'), plt.title("原始图像")
plt.subplot(132), plt.imshow(reconstructed_img, cmap='gray'), plt.title("简单量化(保留左上角)")
plt.subplot(133), plt.imshow(advanced_reconstructed, cmap='gray'), plt.title("Z字形量化")
plt.show()

你会发现Z字形扫描的量化方法在相同压缩率下能保留更多视觉上重要的信息，这正是实际JPEG压缩采用的策略。

3. DWT实战：小波变换在图像去噪中的应用

与DCT不同，DWT（离散小波变换）能同时提供频率和位置信息，使其在图像去噪、边缘检测等应用中表现优异。让我们探索如何使用DWT进行图像去噪。

3.1 多级小波分解与重构

我们将使用PyWavelets库实现二维小波变换：

def dwt_decomposition(image, wavelet='haar', level=2):
    # 转换为浮点型
    img_float = np.float32(image) / 255.0
    # 进行小波分解
    coeffs = pywt.wavedec2(img_float, wavelet, level=level)
    return coeffs

def dwt_reconstruction(coeffs, wavelet='haar'):
    # 小波重构
    reconstructed = pywt.waverec2(coeffs, wavelet)
    # 裁剪到[0,1]范围并转换回uint8
    return np.clip(reconstructed * 255, 0, 255).astype(np.uint8)

# 2级小波分解
coeffs = dwt_decomposition(image, 'db1', 2)

可视化小波系数：

def plot_dwt_coeffs(coeffs):
    titles = ['Approximation', 'Horizontal', 'Vertical', 'Diagonal']
    fig, axes = plt.subplots(1, 4, figsize=(16,4))
    for i, (ax, title) in enumerate(zip(axes, titles)):
        if i == 0:
            # 近似分量
            ax.imshow(coeffs[0], cmap='gray')
        else:
            # 细节分量
            ax.imshow(np.abs(coeffs[1][i-1]), cmap='jet')
        ax.set_title(title)
        ax.axis('off')
    plt.show()

plot_dwt_coeffs(coeffs)

你会看到四个子图：近似分量（低频信息）、水平细节、垂直细节和对角细节。噪声通常表现在高频细节分量中。

3.2 小波阈值去噪实现

小波去噪的核心思想是对细节系数进行阈值处理：

def wavelet_denoise(image, wavelet='db1', level=2, threshold=0.1, mode='soft'):
    # 小波分解
    coeffs = dwt_decomposition(image, wavelet, level)
    # 阈值处理细节系数
    new_coeffs = [coeffs[0]]
    for i in range(1, len(coeffs)):
        new_detail = []
        for detail in coeffs[i]:
            if mode == 'soft':
                # 软阈值
                new_detail.append(pywt.threshold(detail, threshold * np.max(detail)))
            else:
                # 硬阈值
                new_detail.append(detail * (np.abs(detail) > (threshold * np.max(detail))))
        new_coeffs.append(tuple(new_detail))
    # 重构图像
    return dwt_reconstruction(new_coeffs, wavelet)

# 添加一些高斯噪声
noisy_image = image + np.random.normal(0, 25, image.shape).astype(np.uint8)
# 去噪
denoised_image = wavelet_denoise(noisy_image, 'sym4', 3, 0.15)

比较去噪效果：

plt.figure(figsize=(15,5))
plt.subplot(131), plt.imshow(image, cmap='gray'), plt.title("原始图像")
plt.subplot(132), plt.imshow(noisy_image, cmap='gray'), plt.title("带噪声图像")
plt.subplot(133), plt.imshow(denoised_image, cmap='gray'), plt.title("小波去噪后")
plt.show()

3.3 小波基选择与参数优化

不同的小波基函数会产生不同的去噪效果。让我们比较几种常用小波基：

wavelets = ['haar', 'db2', 'sym4', 'coif2']
results = {}
for w in wavelets:
    results[w] = wavelet_denoise(noisy_image, w, 3, 0.15)

plt.figure(figsize=(15,10))
for i, (name, img) in enumerate(results.items()):
    plt.subplot(2,2,i+1)
    plt.imshow(img, cmap='gray')
    plt.title(f"{name}小波去噪效果")
plt.tight_layout()
plt.show()

选择合适的小波基需要考虑以下因素：

• 对称性：对称小波(sym, coif)能减少边界失真
• 支撑长度：较长的小波能更好捕捉平滑区域特征
• 正则性：影响重建图像的平滑度

注意：在实际应用中，通常需要通过交叉验证选择最优的小波基和阈值参数。对于自然图像，sym4和coif2通常是不错的选择。

4. DCT与DWT的综合对比与应用选择

现在我们已经实践了两种变换，让我们从几个关键维度进行对比：

4.1 技术特性对比

特性	DCT	DWT
变换类型	块变换(通常8×8)	全局多分辨率变换
频率定位	仅频率信息	频率+空间位置信息
压缩效率	高(JPEG标准)	更高(JPEG2000标准)
计算复杂度	较低	较高
边界效应	块效应明显	边界扩展可控制
适用场景	有损压缩	压缩、去噪、特征提取
专利状态	专利已过期	部分小波基仍有专利限制

4.2 实际应用选择指南

根据项目需求选择合适的变换：

• 选择DCT当：

  • 需要兼容JPEG标准
  • 计算资源有限
  • 可接受轻度块效应
  • 主要目标是压缩而非分析

• 选择DWT当：

  • 需要更高质量的压缩(如医学影像)
  • 同时需要去噪或特征提取
  • 能承受稍高的计算成本
  • 需要多分辨率分析

4.3 混合应用案例：JPEG2000压缩

JPEG2000标准结合了DWT和先进的编码技术，实现了比传统JPEG更好的压缩效率。以下是简化的实现思路：

def jpeg2000_like_compression(image, wavelet='bior4.4', level=5, keep_fraction=0.1):
    # 小波分解
    coeffs = pywt.wavedec2(np.float32(image)/255.0, wavelet, level=level)
    # 系数排序和截断
    coeff_arr, coeff_slices = pywt.coeffs_to_array(coeffs)
    sorted_coeff = np.sort(np.abs(coeff_arr.ravel()))[::-1]
    threshold = sorted_coeff[int(len(sorted_coeff) * keep_fraction)]
    # 阈值处理
    new_coeff_arr = coeff_arr * (np.abs(coeff_arr) >= threshold)
    new_coeffs = pywt.array_to_coeffs(new_coeff_arr, coeff_slices, output_format='wavedec2')
    # 重构
    reconstructed = pywt.waverec2(new_coeffs, wavelet)
    return np.clip(reconstructed * 255, 0, 255).astype(np.uint8)

# 10%系数保留
j2k_image = jpeg2000_like_compression(image, keep_fraction=0.1)

比较JPEG和JPEG2000风格的压缩：

plt.figure(figsize=(15,5))
plt.subplot(131), plt.imshow(image, cmap='gray'), plt.title("原始图像")
plt.subplot(132), plt.imshow(reconstructed_img, cmap='gray'), plt.title("DCT压缩(25%系数)")
plt.subplot(133), plt.imshow(j2k_image, cmap='gray'), plt.title("类JPEG2000压缩(10%系数)")
plt.show()

尽管JPEG2000风格的压缩保留了更少的系数(10% vs 25%)，但视觉质量却更好，特别是边缘和纹理区域。