Python实现自适应多维Haar小波零系数消除的无损图像压缩算法

本文将深入探讨一种高效的无损图像压缩算法。该算法基于自适应多维Haar小波变换，通过对图像进行多分辨率分析，识别并消除零值或接近零值的小波系数，从而实现数据的高比例压缩，同时保证原始信息的完整重构。

算法核心思想

该算法的核心是小波变换的多分辨率特性。它首先将图像分解为不同尺度的系数表示（尺度系数和细节系数）。其中，细节系数（高频分量）通常包含大量零值或接近零值。算法的精髓在于自适应地检测这些零系数区域，并通过对离散化网格进行“粗化”（Coarsening）来合并冗余节点，从而在数学上等价地减少表示图像所需的数据总量，整个过程完全可逆，确保了压缩的无损性。

算法流程概览

以下是算法的主要步骤流程图，清晰地展示了从输入到压缩验证的完整过程：

开始
├─ 读取输入图像
│   ├─ 转换为灰度图
│   └─ 归一化像素值到[0,1]
├─ 构建离散化网格
│   ├─ 确定基础分辨率级别
│   ├─ 创建均匀细化描述符
│   └─ 建立Morton顺序线性化
├─ 数据重排序
│   └─ 将图像像素按Morton顺序重新排列
├─ 小波变换
│   ├─ 递归应用Hadamard-Haar变换
│   ├─ 生成层次化系数
│   └─ 分离尺度系数和细节系数
├─ 自适应压缩循环
│   ├─ 遍历所有节点
│   │   ├─ 检测全零细节系数的节点
│   │   ├─ 检测部分维度为零系数的节点
│   │   └─ 计划粗化操作
│   ├─ 应用粗化操作
│   │   ├─ 更新离散化网格
│   │   ├─ 重新映射系数
│   │   └─ 删除零系数维度
│   └─ 重复直到无可压缩节点
├─ 逆向变换验证
│   ├─ 从小波系数重构尺度系数
│   ├─ 重新采样生成图像
│   └─ 验证与原始图像的差异
└─ 输出压缩结果
    ├─ 显示压缩前后图像
    ├─ 计算压缩比
    └─ 验证无损性

详细算法步骤

下面，我们拆解流程中的每一个关键步骤。

1. 图像预处理
   读取输入图像（支持TIFF、PNG等格式），并将其转换为灰度图。随后将像素值归一化到[0, 1]区间，为后续的数值计算做准备。

2. 构建离散化网格与数据重排
   根据图像分辨率确定基础网格，并使用Morton顺序（Z-order） 对网格单元进行编码。这种编码能保持空间局部性，提升后续操作的效率。接着，将图像像素从传统的行优先顺序重排为Morton顺序，使空间相邻的像素在数据数组中也相邻。

3. 执行正向小波变换
   采用递归的Hadamard-Haar小波变换，自底向上地将节点系数转换为层次化的小波系数。每个父节点会生成一个表征低频信息的尺度系数和多个表征高频细节的细节系数。这是实现零系数消除算法的前提。

4. 零系数检测与自适应压缩
   这是算法的核心压缩阶段。遍历网格中的所有节点：

   • 完全粗化：若一个节点的所有细节系数均为零，则计划将其子节点完全合并。
   • 部分粗化：若一个节点仅在部分维度的细节系数为零，则计划在这些特定维度上进行合并。
     然后，应用计划好的粗化操作：更新网格结构、合并系数、删除冗余的零系数维度。此过程循环迭代，直至没有可压缩的节点为止。

5. 逆向变换与结果验证
   完成压缩后，执行逆向小波变换，从压缩后的小波系数重构出尺度系数，并上采样得到完整图像。通过计算重构图像与原始图像的像素级差异，验证误差在机器精度范围内（如小于1e-8），从而严格证明压缩的无损性。

Python代码实现关键模块

以下是使用Python实现上述核心步骤的关键代码片段，依赖于numpy、PIL及专门的dyada库进行离散化网格操作。

# 导入所需库
import numpy as np
from PIL import Image
import dyada
import dyada.discretization
import dyada.linearization

def hadamard_haar_1d() -> np.ndarray:
    """生成一维Hadamard-Haar变换矩阵 (2x2)"""
    coeff = np.ones([2, 2], dtype=np.int8)
    coeff[1, 1] = -1
    return coeff

@lru_cache(maxsize=7)
def coefficient_matrix(dimensionality: int, one_d_transform) -> np.ndarray:
    """递归生成d维Hadamard-Haar变换矩阵（Kronecker积）"""
    if dimensionality == 0:
        return np.array([[1]], dtype=np.int8)
    elif dimensionality == 1:
        return one_d_transform()
    return np.kron(one_d_transform(), coefficient_matrix(dimensionality - 1, one_d_transform))

def hierarchize(nodal_coefficients: Sequence, num_refined_dimensions: int) -> np.ndarray:
    """正向变换：将节点系数转换为层次系数"""
    # 计算层次化矩阵
    hierarchization_mat = 0.5**num_refined_dimensions * coefficient_matrix(
        num_refined_dimensions, hadamard_haar_1d
    )
    return np.matmul(hierarchization_mat, nodal_coefficients, dtype=np.float32)

def read_img_file(filename: str) -> np.ndarray:
    """读取图像并返回归一化灰度数组"""
    img = Image.open(filename).convert("L")  # 转灰度
    img_array = np.array(img, dtype=np.float32) / 255.0  # 归一化
    return img_array

def transform_to_all_wavelet_coefficients(discretization, nodal_coefficients):
    """将整个网格的节点系数转换为小波系数表示"""
    coefficients = [None for _ in range(len(discretization.descriptor))]
    # ... (反向遍历描述符，递归调用hierarchize)
    # 详细实现涉及对树结构的遍历和系数分配
    return coefficients

# ---- 主程序逻辑概览 ----
if __name__ == "__main__":
    # 1. 读取并预处理图像
    data = read_img_file("input.png")

    # 2. 构建Morton顺序离散化网格
    input_shape = data.shape
    base_resolution_level = [5, 5]  # 基础分辨率级别
    uniform_descriptor = dyada.descriptor.RefinementDescriptor(len(input_shape), base_resolution_level)
    discretization = dyada.discretization.Discretization(
        linearization=dyada.linearization.MortonOrderLinearization(),
        descriptor=uniform_descriptor,
    )

    # 3. 将图像数据重排为Morton顺序
    ordered_coefficients = np.zeros((len(discretization)), dtype=np.float32)
    # ... (将data中的像素按Morton索引填入ordered_coefficients)

    # 4. 执行正向小波变换，得到层次系数
    wavelet_coeffs = transform_to_all_wavelet_coefficients(discretization, ordered_coefficients)

    # 5. 自适应零系数消除循环 (核心压缩步骤)
    while True:
        # 计划粗化操作 (检测全零/部分零系数节点)
        planner = dyada.refinement.PlannedAdaptiveRefinement(discretization)
        # ... (遍历节点，根据系数为零的条件调用 planner.plan_coarsening)

        if no_more_coarsening_possible:
            break

        # 应用粗化：更新网格和系数映射
        new_discretization, mapping = planner.apply_refinements(track_mapping="patches")
        # ... (根据映射关系合并、删除系数，更新wavelet_coeffs)
        discretization = new_discretization

    # 6. 逆向变换验证
    # ... (从压缩后的wavelet_coeffs重构图像)
    reconstructed_image = get_resampled_image(discretization, scaling_coeffs)

    # 7. 验证无损性
    error = np.max(np.abs(original_image - reconstructed_image))
    print(f"最大重构误差: {error}")
    assert error < 1e-8, "压缩不是无损的！"

在Python中实现此类算法时，需特别注意数值精度和递归函数的效率管理。

算法特点与应用前景

该自适应算法具有以下显著优势：

• 严格无损：基于可逆的小波变换和系数粗化，确保原始数据100%恢复。
• 高压缩比：针对平滑或具有大量均匀区域的图像（如医学图像、卫星图像），能有效消除冗余，获得比传统无损编码（如PNG）更高的压缩率。
• 维度可扩展：算法设计基于多维小波，可扩展至3D体数据或更高维数据的压缩。

在人工智能与计算机视觉领域，无损压缩对于需要保留原始所有信息的任务至关重要，例如遥感影像分析、医学影像存档和科学计算数据的存储与交换。本算法提供了一种结合多分辨率分析和自适应稀疏化的高效解决方案。