21、小波域图像压缩技术详解

pytorchlight8

于 2025-11-17 16:29:21 发布

阅读量23

点赞数

CC 4.0 BY-SA版权

分类专栏： MATLAB图像视频压缩精讲文章标签：小波域图像压缩离散小波变换 DWT系数量化

本文链接：https://blog.youkuaiyun.com/pytorchlight8/article/details/155179398

MATLAB图像视频压缩精讲专栏收录该内容

28 篇文章 ¥499.90

订阅专栏¥69.90

会员秒杀 ¥9.9 重磅福利

超级会员免费看

小波域图像压缩技术详解

1. 小波域图像压缩基础

1.1 离散小波变换（DWT）系数量化

在小波域进行图像压缩时，首先需要对 L 级 DWT 系数进行均匀量化。量化过程使用从最优比特分配规则计算出的量化步长。具体操作步骤如下：
1. 输入参数 ：
- C ：通过 wavedec2 函数获得的 DWT 系数向量。
- S ：同样由 wavedec2 函数得到的记录变量。
- Qsteps ：对应于最优比特或整数比特分配的量化步长。
2. 计算 DWT 级别数 ：

L = size(S,1) - 2; % 计算 DWT 级别数

提取 DWT 系数 ：

appCoefLength = S(1,1)*S(1,2); % 计算近似系数的长度
for k = 1:L
    k1 = (k-1)*3+1;
    k2 = k1 +2;
    detCoefLength(k1:k2) = S(k+1,1)*S(k+1,2);
end
temp = C(1:appCoefLength);

量化近似系数 ：

temp = floor(temp/Qsteps(3*L+1)+0.5)*Qsteps(3*L+1);
y(1:appCoefLength) = temp;
startLoc = appCoefLength; % 近似系数向量的长度

量化细节系数 ：

for k = 1:3*L
    temp = C(startLoc+1:startLoc+detCoefLength(k));
    temp = floor(temp/Qsteps(k)+0.5)*Qsteps(k);
    y(startLoc+1:startLoc+detCoefLength(k)) = temp;
    startLoc = startLoc + detCoefLength(k);
end

1.2 量化 DWT 系数的熵编码

小波基图像压缩的最后一步是熵编码。由于量化后大多数细节系数为零，因此可以使用游程编码将系数转换为符号，以便存储或传输。可以使用任何熵编码方法，如 Huffman、算术或 Golomb - Rice 编码来对量化后的 DWT 细节系数进行编码。近似系数可以先使用预测或 DCT 方法进行压缩，然后再进行 Huffman 编码。

2. 零树编码

2.1 零树的概念

“零树”指的是由零系数构成的树。在二维 DWT 的倍频带中，每一级的 LL 带会以四叉树的方式划分为四个属于下一级的子带。因此，除 LL 带外，给定级别和方向的每个子带中的一个系数，在相同方向的下一级中都有一组对应的 4 个系数（实际为 2×2）。在 DWT 中，如果 HL、LH 或 HH 带中最高级别的某个系数为零或不显著，那么其在较低级别中对应的所有系数很可能也为零或不显著。零树编码利用了这一特性来实现高压缩率。

2.2 嵌入式零树小波（EZW）编码

EZW 是一种非常巧妙的编码方法，它将系数位置信息的编码与系数值的逐次逼近相结合。“嵌入式”意味着编码（或解码）可以在任何所需的比特率下终止，因为粗系数值及其细化信息都嵌入在编码比特流中。

EZW 编码包括两个过程：一个用于编码系数位置信息，另一个用于编码近似系数值。这两个过程会持续进行，直到达到所需的比特预算。在 EZW 中，需要扫描所有 DWT 系数，将其转换为符号进行存储或传输。扫描方法通常是从最高级别（父系数）到最低级别（远亲系数）进行。显著系数值使用位平面编码。

2.3 EZW 编码示例

以一个 8×8 的摄像师图像块为例，在进行了 128 的电平偏移后，其像素值如下表所示：
| 58 | 63 | 59 | 60 | 59 | 53 | 27 | 16 |
|----|----|----|----|----|----|----|----|
| 51 | 57 | 57 | 52 | 58 | 62 | 50 | 26 |
| 36 | 48 | 46 | 34 | 43 | 49 | 51 | 41 |
| 29 | 43 | 45 | 31 | 40 | 48 | 51 | 56 |
| 51 | 50 | 54 | 47 | 47 | 46 | 45 | 46 |
| 48 | 46 | 49 | 48 | 44 | 42 | 42 | 36 |
| 26 | 20 | 25 | 27 | 21 | 18 | 15 | 12 |
| 4 | 6 | 7 | 4 | 4 | 1 | -10 | -17 |

其三级二维 DWT 系数如下表所示：
| 300 | 69 | 31 | -12 | 7 | 5 | -5 | -18 |
|----|----|----|----|----|----|----|----|
| 29 | -11 | 71 | 82 | 4 | 1 | 3 | -8 |
| -1 | 25 | -2 | 38 | 4 | 1 | 2 | 5 |
| -4 | 12 | 1 | -13 | 17 | 19 | 17 | 27 |
| -7 | 2 | 1 | 19 | 0 | -3 | 5 | -7 |
| -14 | 14 | -8 | 3 | 0 | -2 | 1 | 7 |
| 2 | 3 | 1 | -1 | -1 | 2 | -1 | -3 |
| 2 | -1 | 3 | 4 | 4 | -3 | -1 | -3 |

EZW 编码步骤如下：
1. 确定初始阈值 ：由于最大绝对系数值为 300，初始阈值 T0 = 2⌊log2 300⌋ = 256 。
2. 第一遍扫描 ：
- 从 LL3 开始扫描系数，300 大于 T0 ，被分配符号 POS ，量化区间为 (256,512]，中点为 384，因此 300 被量化为 384。
- HL3 系数值 69 小于 T0 ，且其后代系数也不显著，被分配符号 ZTR ，量化为 0。
- 同理，LH3 和 HH3 系数及其后代系数也不显著，被分配符号 ZTR ，量化为 0。
- 第一遍扫描后的主导列表如下表所示：
| 位置 | 系数值 | 符号 | 量化值 |
|----|----|----|----|
| LL3 | 300 | POS | 384 |
| HL3 | 69 | ZTR | 0 |
| LH3 | 20 | ZTR | 0 |
| HH3 | -11 | ZTR | 0 |
- 第一遍扫描后的从属列表如下表所示：
| 系数值 | 代码 | 量化值 |
|----|----|----|
| 300 | 0 | 320 |
3. 后续扫描 ：在后续的扫描中，阈值会逐次减半，重复上述过程，直到达到所需的比特预算或阈值变为 1。

2.4 零树编码的优势

与基于 DCT 的编码器不同，基于小波的编码器需要扫描更大的块来挑选用于熵编码的系数，因此需要花费大量比特来编码系数位置信息。而 EZW 编码技术通过利用零树的特性，减少了这方面的开销。

3. JPEG2000 图像压缩标准

3.1 JPEG2000 概述

JPEG2000 是 2000 年由 ISO 标准化的一种基于 DWT 的图像压缩算法。它与 JPEG 类似，规定了解码过程、码流语法、文件格式，并提供了编码过程的指导和实际实现方法。JPEG2000 标准同时支持无损和有损压缩算法。

3.2 JPEG2000 编码流程

JPEG2000 的整体编码流程如下：

graph LR
    A[输入图像] --> B[DC 电平偏移]
    B --> C[正向分量变换]
    C --> D[分块]
    D --> E[二维 DWT]
    E --> F[标量量化]
    F --> G[熵编码]
    G --> H[压缩比特流]
    H --> I[熵解码]
    I --> J[标量反量化]
    J --> K[二维 IDWT]
    K --> L[分块转光栅扫描]
    L --> M[反向分量变换]
    M --> N[添加 DC 电平]
    N --> O[解压缩图像]

3.3 JPEG2000 编码步骤详解

3.3.1 DC 电平偏移

输入图像分量通过减去一个固定整数进行 DC 电平偏移，公式为：
[I’_k [m, n] = I_k [m, n] - 2^{size_k - 1}]
其中， size_k 是第 k 个输入分量的最大可能值。在反向分量变换后，需要将 DC 电平加回到每个像素上：
[I_k [m, n] = I’_k [m, n] + 2^{size_k - 1}]

3.3.2 正向分量变换（FCT）

为了实现更高的压缩率，将电平偏移后的图像分量（如红、绿、蓝）转换到另一组坐标系中。对于不可逆（有损）变换，FCT 表示为：
[\begin{bmatrix} Y_0 \ Y_1 \ Y_2 \end{bmatrix} = \begin{bmatrix} 0.299 & 0.587 & 0.144 \ -0.16875 & -0.33126 & 0.5 \ 0.5 & -0.41869 & -0.08131 \end{bmatrix} \begin{bmatrix} I_0 \ I_1 \ I_2 \end{bmatrix}]
对于可逆（无损）变换，FCT 由以下方程描述：
[\begin{bmatrix} Y_0 \ Y_1 \ Y_2 \end{bmatrix} = \begin{bmatrix} \frac{I_0 + 2I_1 + I_2}{4} \ I_2 - I_1 \ I_0 - I_1 \end{bmatrix}]

3.3.3 反向分量变换（ICT）

不可逆和可逆情况下的 ICT 分别由以下方程表示：
[\begin{bmatrix} I_0 \ I_1 \ I_2 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 1.402 \ 1 & -0.34413 & -0.71414 \ 1 & 1.772 & 0 \end{bmatrix} \begin{bmatrix} Y_0 \ Y_1 \ Y_2 \end{bmatrix}]
[\begin{bmatrix} I_1 \ I_0 \ I_2 \end{bmatrix} = \begin{bmatrix} Y_0 - \frac{Y_1 + Y_2}{4} \ Y_2 + I_1 \ Y_1 + I_1 \end{bmatrix}]

3.3.4 图像分量分块

将图像分量划分为块，这些块类似于 JPEG 中的宏块。块是图像分量的矩形子块，按光栅扫描顺序编号。所有块的大小相同。在分块之前，色度分量会进行 4:2:2 或 4:2:0 的子采样。

3.3.5 块分量的二维 DWT

分块后，每个分量的每个块使用 Daubechies 的 9/7 FIR 滤波器组进行 L 级二维 DWT 分解（不可逆变换）。Daubechies 的 9/7 FIR 滤波器组的分析和合成滤波器的脉冲响应分别如下表所示：
| n | LPF（分析） | HPF（分析） | LPF（合成） | HPF（合成） |
|----|----|----|----|----|
| 0 | 0.6029490182 | 1.1150870525 | 1.1150870525 | 0.6029490182 |
| ±1 | 0.2668641184 | -0.5912717631 | 0.5912717631 | -0.2668641184 |
| ±2 | -0.0782232665 | -0.0575435263 | -0.0575435263 | -0.0782232665 |
| ±3 | -0.0168641184 | 0.0912717631 | -0.0912717631 | 0.0168641184 |
| ±4 | 0.0267487574 | | | 0.0267487574 |

对于可逆压缩，使用 Le Gall 的 5/3 滤波器。Le Gall 的 5/3 滤波器的脉冲响应如下表所示：
| 采样索引 | 低通（分析） | 高通（分析） | 低通（合成） | 高通（合成） |
|----|----|----|----|----|
| 0 | 6/8 | 1 | 1 | 6/8 |
| ±1 | 2/8 | -1/2 | 1/2 | -2/8 |
| ±2 | -1/8 | -1/8 | -1/8 | -1/8 |

3.3.6 标量量化

在 JPEG2000 标准中，子带“b”的均匀量化器的步长 Q(b) 计算公式为：
[Q(b) = 2^{R_b + g_b - \epsilon_b} \left(1 + \frac{\mu_b}{2^{11}}\right)]
其中， R_b 是子带“b”的标称动态范围， g_b 是子带的增益， \epsilon_b 和 \mu_b 是系数的指数和尾数，它们可以在比特流中为每个子带显式信号，或者仅为 LL 带隐式信号。增益在低通方向为 1，在高通方向为 2。因此，每个级别的 HL 和 LH 带的增益为 2，HH 带的增益为 4。使用以下公式对给定子带中的系数进行量化：
[X_{qb}(u, v) = sign(X_b(u, v)) \left\lfloor \frac{|X_b(u, v)|}{Q(b)} \right\rfloor]

3.3.7 频率加权

为了在不牺牲视觉质量的前提下提高压缩效率，JPEG2000 建议对均匀量化器的量化步长进行频率加权。这些频率加权应用于亮度分量的不同子带，基于观看距离。JPEG2000 推荐的频率加权值如下表所示：
| 级别 | 1000 像素观看距离 | | | 2000 像素观看距离 | | | 4000 像素观看距离 | | |
|----|----|----|----|----|----|----|----|----|----|
| | HL | LH | HH | HL | LH | HH | HL | LH | HH |
| 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| 2 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0.731668 |
| 3 | 1 | 1 | 1 | 1 | 1 | 0.727203 | 0.564344 | 0.564344 | 0.285968 |
| 4 | 1 | 1 | 0.727172 | 0.560841 | 0.560841 | 0.284193 | 0.179609 | 0.179609 | 0.043903 |
| 5 | 0.560805 | 0.560805 | 0.284173 | 0.178494 | 0.178494 | 0.043631 | 0.014774 | 0.014774 | 0.000573 |

3.4 JPEG2000 压缩示例

以一个 8 bpp 的 RGB 图像为例，使用 JPEG2000 压缩方案进行有损压缩。具体步骤如下：
1. 电平偏移和分量变换 ：使用公式对图像分量进行电平偏移，并通过不可逆的正向分量变换将其转换到新的坐标系中。将色度分量转换为 4:2:2 采样格式。
2. 二维 DWT ：使用 Daubechies 的 9/7 滤波器组对分量进行三级二维 DWT。
3. 确定量化步长 ：使用公式计算所有子带的均匀量化器的量化步长。
4. 量化和频率加权 ：对每个子带中的每个系数进行量化，先将量化步长乘以适当的频率加权值，然后使用量化公式进行量化。
5. 反量化和逆变换 ：将量化后的 DWT 系数乘以修改后的量化步长进行反量化，然后使用合成滤波器组进行二维逆小波变换。将逆小波变换后的色度分量上采样到全分辨率，最后使用逆分量变换和电平偏移恢复 RGB 图像。

以下是实现该示例的 MATLAB 代码：
```matlab
% Example8_3.m
% Implements JPEG2000 wavelet transform-based compression.
% Implements only the scalar quantization normative of the
% standard. Implements up to NL levels of 2D DWT using
% Daubechies’ 9/7 filter bank.
% Implements 4:2:0 and 4:2:2 sampling formats specified
% by the user.
% Uses the frequency weightings recommended by JPEG2000 standard
% in the uniform scalar quantization of the DWT coefficients.
%
clear
NL = 3; % number of levels of 2D DWT
SamplingFormat = ‘4:2:2’;
%A = imread(‘cameraman.tif’);
A = imread(‘lighthouse.ras’);
% make sure that the image size is divisible by NL
[X,Y,Z] = size(A);
if mod(X,2^NL) ~=0
Height = floor(X/(2^NL)) (2^NL);
else
Height = X;
end
if mod(Y,2^NL) ~=0
Width = floor(Y/(2^NL)) (2^NL);
else
Width = Y;
end
Depth = Z;
clear X Y Z
% Do forward Irreversible Component Transformation (ICT)
% In JPEG2000 jargon, y0 corresponds to the Luma Y,
% y1 to Cb and y2 to Cr.
if Depth == 1
y0 = double(A(1:Height,1:Width)) - 128; % DC level shift
else
A1 = double(A(1:Height,1:Width,:))-128; % DC level shift
y0 = 0.299 A1(:,:,1) + 0.587 A1(:,:,2) + 0.144 A1(:,:,3);
y1 = -0.16875 A1(:,:,1) - 0.33126 A1(:,:,2) + 0.5 A1(:,:,3);
y2 = 0.5 A1(:,:,1) - 0.41869 A1(:,:,2) - 0.08131*A1(:,:,3);
clear A1
switch SamplingFormat
case ‘4:2:0’
y1 = imresize(y1,[Height/2 Width/2],’cubic’);
y2 = imresize(y2,[Height/2 Width/2],’cubic’);
case ‘4:2:2’
y1 = imresize(y1,[Height Width/2],’cubic’);
y2 = imresize(y2,[Height Width/2],’cubic’);
end
end
% Daubechies 9/7 Filters
% Note: The Daubechies 9/7 filter bank has 9 taps for its analysis
% Lowpass and synthesis highpass FIR filters;
% 7 taps for its analysis highpass and synthesis lowpass
% FIR filters. However, in order to obtain subbands that are
% exact multiples of 2^-n, n being a positive integer, the
% number of taps in the filters must be even. That is why you
% will notice the filter lengths to be 10 instead of 9 or 7.
% The extra tap values are zero.
%
% LO_D = analysis lowpass filter
% LO_R = synthesis lowpass filter
% HI_D = analysis highpass filter
% HI_R = synthesis highpass filter
LO_D = [0 0.0267487574 -0.0168641184 -0.0782232665 0.2668641184…
0.6029490182 0.2668641184 -0.0782232665 -0.0168641184…

0.0267487574];
HI_D = [0 0.0912717631 -0.0575435263 -0.5912717631…
1.1150870525 -0.5912717631 -0.0575435263 0.0912717631 0 0];
%
LO_R = [0 -0.0912717631 -0.0575435263 0.5912717631…
1.1150870525 0.5912717631 -0.0575435263 -0.0912717631 0 0];
HI_R = [0 0.0267487574 0.0168641184 -0.0782232665 -0.2668641184…
0.6029490182 -0.2668641184 -0.0782232665 0.0168641184…
0.0267487574];
%
% “dwt_y0” cell array has NL x 4 cells.
% In each cell at level NL, the 1st component is the LL
% coefficients,
% 2nd component is HL coefficients, 3rd component LH and
% 4th component HH coefficients.
%
% Do an NL-level 2D DWT of the Y component
dwt_y0 = cell(NL,4);
for n = 1:NL
if n == 1
[dwt_y0{n,1},dwt_y0{n,2},dwt_y0{n,3},dwt_y0{n,4}] =…
dwt2(y0,LO_D,HI_D,’mode’,’per’);
else
[dwt_y0{n,1},dwt_y0{n,2},dwt_y0{n,3},dwt_y0{n,4}] =…
dwt2(dwt_y0{n-1,1},LO_D,HI_D,’mode’,’per’);
end
end
% Quantize the Y DWT coefficients
[dwt_y0,y0Bits] = jpg2K_Quantize(dwt_y0,NL,1);
%
% Do an NL-level 2D IDWT of the Y component
% The cell array idwt_img has NL cells with each cell
% corresponding to
% the reconstructed LL image at that level.
idwt_y0 = cell(NL,1);
for n = NL:-1:1
if n == 3
idwt_y0{n} = idwt2(dwt_y0{n,1},dwt_y0{n,2},dwt_y0{n,3},…
dwt_y0{n,4},LO_R,HI_R,’mode’,’per’);
else
idwt_y0{n} = idwt2(idwt_y0{n+1},dwt_y0{n,2},dwt_y0{n,3},…
dwt_y0{n,4},LO_R,HI_R,’mode’,’per’);
end
end
SNRy0 = 20 log10(std2(y0)/std2(y0-idwt_y0{1}));
if Depth > 1
% Do an NL-level 2D DWT of the Cb and Cr components
dwt_y1 = cell(NL,4);
dwt_y2 = cell(NL,4);
for n = 1:NL
if n == 1
[dwt_y1{n,1},dwt_y1{n,2},dwt_y1{n,3},dwt_y1{n,4}] =…
dwt2(y1,LO_D,HI_D,’mode’,’per’);
[dwt_y2{n,1},dwt_y2{n,2},dwt_y2{n,3},dwt_y2{n,4}] =…
dwt2(y2,LO_D,HI_D,’mode’,’per’);
else
[dwt_y1{n,1},dwt_y1{n,2},dwt_y1{n,3},dwt_y1{n,4}] =…
dwt2(dwt_y1{n-1,1},LO_D,HI_D,’mode’,’per’);
[dwt_y2{n,1},dwt_y2{n,2},dwt_y2{n,3},dwt_y2{n,4}] =…
dwt2(dwt_y2{n-1,1},LO_D,HI_D,’mode’,’per’);
end
end
%
% Quantize the Cb and Cr DWT coefficients
[dwt_y1,y1Bits] = jpg2K_Quantize(dwt_y1,NL,2);
[dwt_y2,y2Bits] = jpg2K_Quantize(dwt_y2,NL,3);
%
% % Do an NL-level 2D IDWT of the Cb and Cr components
idwt_y1 = cell(NL,1);
idwt_y2 = cell(NL,1);
%
for n = NL:-1:1
if n == 3
idwt_y1{n} = idwt2(dwt_y1{n,1},dwt_y1{n,2},dwt_y1{n,3},…
dwt_y1{n,4},LO_R,HI_R,’mode’,’per’);
idwt_y2{n} = idwt2(dwt_y2{n,1},dwt_y2{n,2},dwt_y2{n,3},…
dwt_y2{n,4},LO_R,HI_R,’mode’,’per’);
else
idwt_y1{n} = idwt2(idwt_y1{n+1},dwt_y1{n,2},dwt_y1{n,3},…
dwt_y1{n,4},LO_R,HI_R,’mode’,’per’);
idwt_y2{n} = idwt2(idwt_y2{n+1},dwt_y2{n,2},dwt_y2{n,3},…
dwt_y2{n,4},LO_R,HI_R,’mode’,’per’);
end
end
% Compute the SNR of Cb and Cr components due to quantization
SNRy1 = 20 log10(std2(y1)/std2(y1-idwt_y1{1}));
SNRy2 = 20 log10(std2(y2)/std2(y2-idwt_y2{1}));
% Upsample Cb & Cr components to original full size
y1 = imresize(idwt_y1{1},[Height Width],’cubic’);
y2 = imresize(idwt_y2{1},[Height Width],’cubic’);
end
% Do inverse ICT & level shift
if Depth == 1
Ahat = idwt_y0{1} + 128;
figure,imshow(uint8(Ahat))
sprintf(‘Average bit rate = %5.2f bpp\n’,y0Bits/(Height Width))
sprintf(‘SNR(Y) = %4.2f dB\n’,SNRy0)
else
Ahat = zeros(Height,Width,Depth);
Ahat(:,:,1) = idwt_y0{1} + 1.402 y2 + 128;
Ahat(:,:,2) = idwt_y0{1} - 0.34413 y1 - 0.71414 y2 + 128;
Ahat(:,:,3) = idwt_y0{1} + 1.772 y1 + 128;
figure,imshow(uint8(Ahat))
sprintf(‘Average bit rate = %5.2f bpp\n’,…
(y0Bits+y1Bits+y2Bits)/(Height Width))
SNRy0 = 20 log10(std2(y0)/std2(y0-idwt_y0{1}));
sprintf(‘SNR(Y) %4.2fdB\tSNR(Cb) = %4.2f dB\tSNR(Cr) = %4.2fdB\n’,…
SNRy0,SNRy1,SNRy2)
end

function [y,TotalBits] = jpg2K_Quantize(x,NL,CompNo)
% [y,TotalBits] = jpg2K_Quantize(x,NL,CompNo)
% Quantizes using JPEG2000 scalar quantization rule
% Input:
%
% x = cell array containing the 2D DWT coefficients
%
% NL = number of levels of octave-band 2D DWT
%
% CompNo = Component number, which can be 1, 2, or 3
%
% 1 for Y, 2 for y1 (Cb) and 3 for y2 (Cr)
% Output:
%
% y = cell array of the same size as x containing
%
% quantized DWT coefficients
%
% TotalBits = Total bits used up by the coefficients
%
% The input cell array x{1,2:4} contains the 1HL, 1LH, and 1HH
% DWT coefficients at level 1, etc., and x{NL,1:4}
% contains the LL, HL, LH, and HH coefficients at level NL.
%
% The quantization step sizes are usually specified in the
% headers of the bit stream. However, here, we simply
% assign the number of bits for the uniform quantizers
% and compute the appropriate quantization steps.
%
y = cell(NL,4);% output cell array
maxVal = zeros(4,1);% maximum coefficient value
%minVal = zeros(4,1);
Qstep = zeros(NL,4);% array to store the quantization steps
Qbits = zeros(NL,4);% array to store the quantizer bits
switch CompNo
case 1
Qstep(1,2:3)=(2^(8-0)) (1+8/(2^11));
Qstep(1,4)=(2^(8-0)) (1+8/(2^11));
Qstep(2,2:3)=(2^(8-3)) (1+8/(2^11));
Qstep(2,4)=(2^(8-2)) (1+8/(2^11));
Qstep(2,4) = Qstep(2,4)/0.731668; % uses CSF
Qstep(3,2:3)=(2^(8-5)) (1+8/(2^11));
Qstep(3,2:3) = Qstep(3,2:3)/0.564344; % uses CSF
Qstep(3,4)=(2^(8-5)) (1+8/(2^11));
Qstep(3,4)= Qstep(3,4)/0.285968; % uses CSF
Qstep(3,1)=(2^(8-5)) (1+8/(2^11));
case {2,3}
Qstep(1,2:3)=(2^(8-0)) (1+8/(2^11));
Qstep(1,4)=(2^(8-1)) (1+8/(2^11));
Qstep(2,2:3)=(2^(8-3)) (1+8/(2^11));
Qstep(2,4)=(2^(8-2)) (1+8/(2^11));
Qstep(3,2:3)=(2^(8-4)) (1+8/(2^11));
Qstep(3,4)=(2^(8-3)) (1+8/(2^11));
Qstep(3,1)=(2^(8-5)) (1+8/(2^11));
end
TotalBits = 0;
for n = 1:NL
if n < NL
m1 = 2;
else
m1 = 1;
end
for m = m1:4
maxVal(m) = max(x{n,m}(:));
t = round(log2(maxVal(m)/Qstep(n,m)));
if t < 0
Qbits(n,m)= 0;
else
Qbits(n,m)= t;
end
TotalBits = TotalBits+Qbits(n,m) size(x{n,m},1)
size(x{n,m},2);
end
end
%
for n = 1:NL
if n < NL
m1 = 2;
else
m1 = 1;
end
for m = m1:4
s = sign(x{n,m});
q2 = Qstep(n,m)/2;
y{n,m} = s . round(abs(x{n,m})/q2) q2;
end
end
end

3.5 JPEG2000 压缩效果

使用 JPEG2000 压缩方案对鸟类图像进行不可逆压缩，示例结果如下：
- 平均比特率为 0.64 bpp，如果使用熵编码，比特率可能会更低。
- 亮度分量的 SNR 为 19.58 dB，色度分量的 SNR 分别为 19.78 dB 和 20.25 dB。
- 视觉质量在低比特率下非常好，没有像基于 DCT 的 JPEG 标准那样的块效应。

不同图像使用 JPEG2000 压缩方案的比特率和 SNR 如下表所示：
| 图像 | 比特率 (bpp) | Y (dB) | Cb (dB) | Cr (dB) |
| ---- | ---- | ---- | ---- | ---- |
| 摄像师 | 0.48 | 15.95 | NA | NA |
| 鸟类 | 0.64 | 19.58 | 19.78 | 20.25 |
| 航拍 | 0.52 | 12.12 | NA | NA |
| 飞机 | 0.53 | 22.21 | 14.30 | 9.62 |
| 增田 | 0.51 | 24.08 | 13.95 | 17.09 |
| 辣椒 | 0.46 | 20.84 | 19.65 | 19.00 |
| 戴帽女士 | 0.58 | 23.29 | 14.63 | 17.88 |
| 游艇 | 0.52 | 18.10 | 14.42 | 14.68 |
| 灯塔 | 0.59 | 13.18 | 14.68 | 15.35 |

从表中可以看出，整体平均比特率约为 0.5 bpp，所有重建图像都具有良好的视觉质量，没有块效应。

4. 总结

4.1 小波域图像压缩技术对比

技术	优势	劣势	适用场景
离散小波变换（DWT）系数量化	利用 DWT 特性，结合最优比特分配规则量化系数，实现图像初步压缩	量化过程可能损失部分细节信息	对图像质量要求不是极高，追求一定压缩率的场景
零树编码（EZW）	利用零树特性，结合系数位置信息编码和逐次逼近，实现高压缩率	扫描系数和编码过程相对复杂，计算量较大	对压缩率要求较高，对计算资源有一定承受能力的场景
JPEG2000	支持无损和有损压缩，采用多种技术结合，压缩效果好，视觉质量高，无块效应	编码流程复杂，实现难度较大	对图像质量和压缩率都有较高要求的场景，如专业图像存储和传输

4.2 技术选择建议

如果需要快速简单的压缩，对图像质量要求一般，可优先考虑离散小波变换（DWT）系数量化。
若追求高压缩率，且计算资源允许，零树编码（EZW）是一个不错的选择。
对于对图像质量和压缩率都有严格要求的专业应用，如医学图像、卫星图像等，JPEG2000 标准是最佳方案。

4.3 未来发展趋势

随着图像数据量的不断增加和对图像质量要求的提高，小波域图像压缩技术将不断发展。未来可能会在以下方面取得进展：
- 算法优化 ：进一步优化现有算法，提高压缩效率和图像质量，减少计算复杂度。
- 与其他技术融合 ：结合深度学习等新兴技术，实现更智能、更高效的图像压缩。
- 应用拓展 ：在更多领域得到应用，如虚拟现实、高清视频等。

通过对小波域图像压缩技术的深入了解和合理应用，可以在不同场景下实现高效的图像压缩，满足各种需求。