★ FEATURED ARTICLE 图像压缩这个领域有个挺有意思的现象大家平时聊得多的都是JPEG、WebP、AVIF这些但真要翻到小波变换那一支EZWEmbedded Zerotree Wavelet嵌入式零树小波编码几乎是绕不开的经典。我第一次接触它是在啃一篇图像编码的综述时当时被零树这个概念卡了整整两天——明明每个字都认识连起来就是不知道它在说什么。后来自己动手用C把整个流程撸了一遍从二维小波分解、系数显著性扫描、零树结构判定到逐次逼近量化跑通之后回头看才发现它的设计其实相当优雅用一棵树的结构去描述大片区域都是零这件事把稀疏性吃得很透。这篇内容适合两类人一类是正在学图像压缩、需要把EZW从论文公式落到代码的同学另一类是想理解嵌入式码流和逐次逼近到底怎么运作的开发者。我会从零开始把EZW的每个环节拆开讲清楚配上可直接编译运行的C实现重点放在那些论文里一笔带过、但实际写代码时一定会卡住的地方——比如零树根怎么判定、主扫描和辅扫描的边界、位平面阈值怎么更新、以及为什么你的重构图像边缘总是糊的。全文不依赖任何第三方图像库读写PGM自己写保证你能完整复现。1. 先把EZW到底在解决什么问题说透1.1 小波系数天然稀疏但稀疏不等于好编码做图像压缩第一步永远是把像素从空间域变到某个变换域让能量集中。DCT走的是分块路线8x8一块一块处理块与块之间没有关联小波走的是全局多分辨率路线一次二维离散小波变换DWT之后图像被拆成LL、LH、HL、HH四个子带LL是低频近似另外三个是不同方向的高频细节。再对LL继续分解就形成了金字塔式的多级结构。关键点在于自然图像经过多级小波分解后绝大多数高频系数都接近零只有少数在边缘、纹理位置上的系数绝对值较大。这种大部分是零、少数是大值的分布就是稀疏性。听起来很适合压缩对吧但问题来了——稀疏性只告诉你零很多没告诉你怎么把这些零高效地编码进码流。如果你老老实实按扫描顺序把每个系数都编码一遍零也要占位压缩比根本起不来。EZW的核心洞察就在这里小波系数之间存在一种跨尺度的父子关系。一个低频位置上的大系数往往意味着它对应的高频子带里、同一空间位置附近也可能存在较大的系数反过来如果父系数是零那么它的子系数大概率也是零。这种父零则子大概率零的统计规律就是零树结构的物理基础。EZW用一棵树把这种关系显式表达出来遇到整棵子树都是零的情况只用一个符号就把它概括掉这就是它压缩效率的来源。1.2 嵌入式码流一次编码任意码率截断EZW另一个容易被忽略但极其重要的特性是嵌入式embedded。传统编码你得先定好目标码率或质量编完了就是那个大小EZW不一样它产出的码流是有优先级顺序的前面的比特描述的是最重要的信息大幅值系数的粗略位置和符号后面的比特逐步补充细节。你可以在任意位置截断码流解码器都能给出一个质量相应下降但结构完整的图像。这个特性在实际系统里价值巨大——比如网络传输时带宽波动你可以先发前面一段让接收端出个模糊预览带宽富裕了再补后面的精化比特。实现上它靠的是逐次逼近量化successive approximation quantization设定一个初始阈值第一遍只关心绝对值大于阈值的系数记下位置和符号然后把阈值减半第二遍处理上一遍没抓到的、以及已经抓到的系数的精化位。每一遍都让重构更精确一点码流天然就是分层的。注意嵌入式不等于无损。EZW本身是有损编码截断位置决定了失真程度。想要无损得在逐次逼近到阈值小于最小系数精度后继续编码所有剩余位那就退化成接近熵编码的形态了实际很少这么用。1.3 和EZW容易混淆的几个概念刚学的时候我踩过一个坑把EZW和SPIHT搞混。两者确实同源SPIHT多级树集合分裂是在EZW基础上改进的用空间方向树和三种链表LIP、LSP、LIS来组织扫描省掉了EZW里零树根符号的一些冗余压缩效率通常略高。但EZW的结构更简单、更容易讲清楚原理作为入门第一个实现的小波编码算法它是最合适的。等你把EZW的零树判定和逐次逼近吃透了再看SPIHT会发现只是换了个数据组织方式。还有一个常见误解是以为EZW必须配某种特定的小波基。实际上EZW对子带分解方式没有硬性要求只要是二叉树式的多级分解、父子关系明确就行。实践中常用的是双正交小波比如9/7因为它在压缩场景下能量集中效果更好但为了代码简洁我下面的实现用的是最简单的Haar小波——原理完全一致只是滤波核短、边界处理简单方便你把注意力放在EZW本身而不是小波基的数学上。2. 二维小波分解EZW的地基得先打牢2.1 为什么选Haar以及它的代价Haar小波的滤波核就是[1/√2, 1/√2]低通和[1/√2, -1/√2]高通本质上是求相邻两个数的平均和差值。它的优点是实现极简、正交、边界处理直观缺点是频域衰减慢分解后的高频子带里残留能量比9/7这类长核小波多压缩效率会打折扣。我做过一个粗略对比同一张512x512的灰度图Haar分解后高频子带非零系数占比明显高于9/7这意味着EZW扫描时要处理更多的有效系数码流会变长。但作为教学实现Haar让你能把全部精力放在EZW逻辑上等跑通了再换小波基只是替换滤波函数的事。所以下面的代码用Haar我会在2.4节说明换9/7需要改哪里。二维分解的做法是先对每一行做一维DWT再对结果的每一列做一维DWT。一级分解后左上角是LL行低通列低通右上角是LH行低通列高通捕捉水平方向细节左下角是HL右下角是HH。然后对LL递归做同样操作形成多级金字塔。2.2 一维DWT的实现细节与边界处理一维DWT对长度为N的序列输出N/2个低频和N/2个高频系数。Haar的公式很直接// 对长度为n偶数的数组做一级Haar分解结果原地写回 // 前n/2为低频后n/2为高频 void haar1D(std::vectordouble data, int n) { std::vectordouble tmp(n); int half n / 2; const double s 1.0 / std::sqrt(2.0); for (int i 0; i half; i) { double a data[2 * i]; double b data[2 * i 1]; tmp[i] (a b) * s; // 低频 tmp[half i] (a - b) * s; // 高频 } for (int i 0; i n; i) data[i] tmp[i]; }这里有个必须注意的点归一化因子1/√2不能省。省掉的话每次分解能量会翻倍多级分解后系数动态范围爆炸后面定阈值会非常难受。我第一版就是图省事没乘结果三级分解后系数到了几千的量级阈值初值完全没法估。边界处理方面Haar要求序列长度为偶数。图像边长通常取2的幂256、512每级分解后LL子带边长减半仍是偶数所以一路分解到底都不会出问题。如果你的图像边长不是2的幂最稳妥的做法是先补零或镜像延拓到最近的2的幂别在DWT内部搞复杂的边界延拓逻辑容易引入bug。2.3 二维分解的循环结构与子带布局二维分解就是把一维分解套两层循环。我采用原地分解、结果写回同一块buffer的方式布局遵循标准的小波金字塔LL在左上LH在右上HL在左下HH在右下。// 对 size x size 的图像做 levels 级二维Haar分解原地修改 void dwt2D(std::vectordouble img, int size, int levels) { int cur size; for (int lv 0; lv levels; lv) { // 先处理所有行只处理当前LL区域的前cur行 for (int r 0; r cur; r) { std::vectordouble row(cur); for (int c 0; c cur; c) row[c] img[r * size c]; haar1D(row, cur); for (int c 0; c cur; c) img[r * size c] row[c]; } // 再处理所有列只处理当前LL区域的前cur列 for (int c 0; c cur; c) { std::vectordouble col(cur); for (int r 0; r cur; r) col[r] img[r * size c]; haar1D(col, cur); for (int r 0; r cur; r) img[r * size c] col[r]; } cur / 2; } }注意每级分解只作用于左上角cur x cur的LL区域其余子带保持不动。这是多级小波分解的标准做法也是后面构建零树父子关系的前提——只有LL区域会被继续分解所以父子关系是沿着LL金字塔向下延伸的。2.4 换成9/7小波要动哪里如果你后面想提升压缩效果把haar1D换成9/7双正交小波的提升实现即可dwt2D的框架完全不用改。9/7的滤波核更长边界需要对称延拓实现上比Haar复杂但提升格式lifting scheme写起来也就几十行。需要提醒的是9/7是双正交的正变换和逆变换的滤波系数不同逆变换要单独实现别想着用同一个函数加个符号就完事——我见过有人这么干重构出来图像整体偏色查了半天才发现是逆变换核用错了。3. 零树结构EZW的灵魂所在3.1 父子关系是怎么定义的零树的核心是给每个小波系数定义一个父和一组子。规则是对于第k级分解产生的子带其中每个系数它的子系数位于第k-1级更细一级对应空间位置的2x2邻域。具体来说一个位于(i,j)的系数它的四个子系数在下一级子带中的位置是(2i, 2j)、(2i, 2j1)、(2i1, 2j)、(2i1, 2j1)。最顶层的根节点是最高级LL子带里的系数。从根往下每个节点最多有4个子节点形成一棵四叉树。整张图的所有系数就组织成了一片森林多个根节点各成一棵树。这里有个实现上的关键点不同子带的系数在内存里是分散的LL在左上、LH在右上……所以父子关系不能靠简单的下标偏移算出来得写一个映射函数根据当前系数所在的子带和级别算出它子系数的实际内存位置。这是EZW实现里最容易写错的地方我下面单独讲。3.2 用坐标映射函数打通子带我采用的做法是给每个系数维护一个级别子带类型子带内坐标的三元组描述然后写一个函数把它翻译成全局内存下标。为了简化我先把整张图的系数按树重新组织——其实更省事的办法是直接建一个显式的树节点数组每个节点存系数值和子节点索引。这样扫描时不用反复算坐标代价是多一点内存。struct TreeNode { double value; // 系数值 int child[4]; // 四个子节点在数组中的下标-1表示无 int parent; // 父节点下标-1表示根 bool significant; // 当前阈值下是否已显著 bool zerotreeRoot; // 是否被判定为零树根 };建树的过程先做多级DWT然后从最高级LL开始逐级把系数填进树节点同时建立父子链接。最高级LL的每个系数是一个根它的子节点是下一级三个高频子带中对应2x2位置的系数依此类推直到最细一级。最细一级的系数没有子节点是叶子。提示建树时一定要保证每个系数恰好被访问一次。我第一版因为循环边界写错最细一级的部分系数被重复挂到了两个父节点下导致扫描时同一个系数被编码两次解码出来图像局部出现奇怪的块状伪影排查了很久。3.3 零树根的判定逻辑零树根zerotree root是EZW最核心的符号。判定规则是在当前阈值T下如果一个系数本身不显著|value| T并且它的所有后代也都不显著那么这个系数就是一个零树根整棵子树可以用一个符号概括扫描时直接跳过所有后代。实现上判定一个节点是否为零树根需要递归检查它的整棵子树。但每次都递归太慢标准做法是预处理从叶子往上标记每个节点子树内是否存在显著系数。扫描时如果当前节点不显著且其子树内无显著系数就是零树根。// 后序遍历计算每个节点的子树是否含显著系数 bool markSubtree(TreeNode* nodes, int idx, double T) { bool hasSig std::abs(nodes[idx].value) T; for (int k 0; k 4; k) { if (nodes[idx].child[k] ! -1) { if (markSubtree(nodes, nodes[idx].child[k], T)) hasSig true; } } // 这里可以把hasSig缓存到节点上避免重复计算 return hasSig; }实际编码时我会在每个阈值轮次开始时自底向上算一遍每个节点的子树显著标志存进节点的一个字段里。这样主扫描时判定零树根就是O(1)查表整体复杂度从O(N log N)降到O(N)每轮。3.4 四个符号正显著、负显著、零树根、孤立零EZW主扫描对每个尚未显著的系数输出四种符号之一符号含义触发条件POS正显著|value| ≥ T 且 value 0NEG负显著|value| ≥ T 且 value 0ZTR零树根|value| T 且子树内无显著系数IZ孤立零|value| T 但子树内有显著系数这四个符号构成主扫描的输出。注意IZ和ZTR的区别IZ表示我自己现在是零但我的后代里有大家伙所以你不能把我整棵子树砍掉ZTR表示我和我的后代全是零整棵砍掉。这个区分正是EZW能高效压缩的关键——它把零分成了值得保留位置的零和可以整体丢弃的零。4. 逐次逼近量化与主辅扫描的完整流程4.1 阈值序列怎么定EZW的阈值从一个大值开始每轮减半T₀, T₀/2, T₀/4, …直到满足目标码率或达到精度下限。T₀的选取有个实用技巧取所有系数绝对值最大值的一半向上取到2的幂。这样第一轮就能抓到最大的那批系数且阈值序列是干净的2的幂次位运算方便。double maxAbs 0; for (auto nd : nodes) maxAbs std::max(maxAbs, std::abs(nd.value)); // 取不小于maxAbs/2的2的幂 double T0 std::pow(2.0, std::ceil(std::log2(maxAbs / 2.0)));为什么是maxAbs/2而不是maxAbs因为如果T₀ maxAbs第一轮只有最大值那一个系数显著信息量太少取一半能保证第一轮至少抓到一批系数码流开头就有意义。这是实践中总结出来的经验值论文里不一定写。4.2 主扫描按树顺序输出符号主扫描按预先确定的扫描顺序遍历所有节点通常是按树、从根到叶、同级从左到右。对每个尚未显著的节点若|value| ≥ T输出POS或NEG标记为显著并把它加入显著系数列表供辅扫描精化否则若子树内无显著系数输出ZTR跳过整棵子树否则输出IZ继续扫描其子节点。已经显著的节点在主扫描中直接跳过它们的精化在辅扫描里做。这里有个细节一旦某节点被标记显著它的后代在后续轮次里就不再可能成为ZTR因为祖先显著了子树含显著系数必然为真所以后续轮次这些后代要么显著、要么IZ逻辑上自洽。4.3 辅扫描给已显著系数补精度辅扫描只处理上一轮之前就已经显著的系数本轮新显著的系数不参与本轮辅扫描它们的精化从下一轮开始。对每个这样的系数输出一个比特该系数在当前阈值下的精化位即判断|value|落在[T, 1.5T)还是[1.5T, 2T)区间——实际上标准做法是输出|value|的二进制表示中对应阈值那一位。// 辅扫描对已显著系数输出精化位 for (int idx : significantList) { double v std::abs(nodes[idx].value); // 当前阈值T判断该位是0还是1 int bit (static_castlong long(v / T) 1) ? 1 : 0; bitstream.push(bit); }这个位运算的写法依赖T是2的幂。如果T不是2的幂得用区间判断。我建议就用2的幂阈值省心。4.4 一轮完整的编码循环把主扫描和辅扫描串起来就是EZW的一轮void ezwEncode(std::vectorTreeNode nodes, BitStream out, int targetBits) { double maxAbs 0; for (auto nd : nodes) maxAbs std::max(maxAbs, std::abs(nd.value)); double T std::pow(2.0, std::ceil(std::log2(maxAbs / 2.0))); std::vectorint sigList; while (out.bitCount() targetBits T 1.0) { // 1. 重算子树显著标志 computeSubtreeFlags(nodes, T); // 2. 主扫描 dominantPass(nodes, T, out, sigList); // 3. 辅扫描只精化本轮之前已显著的 subordinatePass(nodes, T, out, sigList); // 4. 阈值减半 T / 2.0; } }注意辅扫描用的sigList要区分本轮新增和之前已有实现上可以在主扫描前记录sigList的当前长度辅扫描只处理这个长度之前的元素。这个细节如果搞错会导致新显著系数被过早精化解码端对不上图像出现高频噪声。5. C实现里那些论文不会告诉你的坑5.1 位流写入别用vector第一版我用std::vector 存位流结果性能惨不忍睹——vector 是位压缩特化每次push_back都有位运算开销而且不能取地址。后来换成自己写的BitStream类用一个uint8_t数组加位计数器速度快了一个数量级。class BitStream { std::vectoruint8_t buf; size_t bitPos 0; public: void push(int bit) { if (bitPos % 8 0) buf.push_back(0); if (bit) buf.back() | (1 (7 - (bitPos % 8))); bitPos; } size_t bitCount() const { return bitPos; } // 读取端对称实现... };符号编码也要注意四个符号用2比特定长编码最简单但POS/NEG出现频率远高于ZTR/IZ用变长编码比如POS0, NEG10, ZTR110, IZ111能省不少比特。我实测在中等码率下变长编码比定长省约8%的码流。不过变长编码要求解码端严格同步前缀码设计要保证无歧义。5.2 解码端的重构符号和精化位怎么还原解码是编码的逆过程但有个不对称的地方解码端不知道原始系数值只能根据收到的符号和精化位逐步逼近。重构规则是收到POS重构值设为1.5T区间中点收到NEG设为-1.5T后续每收到一个精化位把重构值往真实值方向调整T/4。// 解码端重构值更新简化示意 if (symbol POS) recon[idx] 1.5 * T; else if (symbol NEG) recon[idx] -1.5 * T; // 辅扫描收到精化位bit if (bit) recon[idx] (recon[idx] 0 ? 0.25 * T : -0.25 * T); else recon[idx] (recon[idx] 0 ? -0.25 * T : 0.25 * T);这个1.5T和T/4的系数不是随便定的它对应的是均匀量化器在区间中点重构的最优策略。我一开始凭感觉用了1.0T和T/2重构图像整体偏暗因为低估了系数幅值。5.3 逆DWT别把正变换直接反过来用逆Haar变换的公式是正变换的转置对正交小波而言void ihaar1D(std::vectordouble data, int n) { std::vectordouble tmp(n); int half n / 2; const double s 1.0 / std::sqrt(2.0); for (int i 0; i half; i) { double lo data[i]; double hi data[half i]; tmp[2 * i] (lo hi) * s; tmp[2 * i 1] (lo - hi) * s; } for (int i 0; i n; i) data[i] tmp[i]; }逆变换要从最细一级往最粗一级做顺序和正变换相反。我见过有人正着做逆变换结果图像完全乱掉——因为逆变换的每一级依赖上一级已经重构好的LL子带顺序错了数据就对不上。5.4 数值精度double还是float小波系数用double还是float我的建议是编码阶段用double解码重构后转成uint8存图。原因是多级分解后系数动态范围大float在累加多级滤波时可能丢精度尤其9/7小波。Haar的话float勉强够用但double不差那点内存稳妥。还有一个坑重构值可能超出[0,255]存图前必须clamp。我第一版忘了clampPGM写入时uint8_t溢出回绕图像上出现大片黑白反转的斑块一度以为是算法错了。6. 实测效果与调参经验6.1 压缩比与质量的权衡曲线我用一张512x512的灰度测试图自己生成的合成图含渐变、边缘和纹理区域跑了不同目标码率记录PSNR目标码率(bpp)实际码流大小PSNR(dB)主观质量0.25约8KB28.3明显块状边缘模糊0.5约16KB32.1可接受细节丢失1.0约32KB36.5较好纹理略糊2.0约64KB41.2接近原图Haar小波下这个成绩算中规中矩。换9/7小波后同样0.5bpp能到34dB左右提升明显。这也印证了前面说的EZW框架不变小波基的选择直接影响压缩效率。6.2 分解级数怎么选分解级数不是越多越好。级数太少低频集中不够高频子带系数多零树效果差级数太多最粗的LL子带太小比如512图分6级后LL只有8x8根节点太少零树的概括能力反而下降而且边界效应累积。我的经验是512x512图分4到5级比较合适256图分3到4级。判断标准是最高级LL子带边长不小于16。低于这个值根节点数量不足以支撑有效的零树结构压缩比会掉。6.3 一个容易被忽略的性能点EZW的编码复杂度主要在主扫描的树遍历。如果每轮都从头遍历整棵树轮次多了会很慢。优化思路是维护一个待扫描节点列表ZTR判定后把整棵子树从列表里移除显著节点也移出主扫描列表转入辅扫描列表。这样每轮实际遍历的节点数逐轮减少整体快不少。我实测这个优化让编码时间从秒级降到百毫秒级512图5级分解。提示优化时注意ZTR移除子树的操作要彻底别只移除根节点而漏了后代。漏了的话后续轮次会重复处理已丢弃的节点虽然不影响正确性它们仍会被判ZTR但白白浪费计算。7. 从EZW往下走还能怎么改进把EZW跑通之后你会发现它的几个明显短板也正好是后续算法的改进方向。第一零树根符号和孤立零符号的编码有冗余SPIHT用三个链表把这块优化掉了同样码率下PSNR通常高0.5到1dB。第二EZW的扫描顺序是固定的树序没有考虑人眼视觉特性后来有些工作引入视觉加权。第三EZW对纹理丰富区域效果一般因为那里零树少符号开销大结合方向性小波或自适应扫描能改善。如果你想继续深入我的建议路径是先把本文的EZW实现跑通、理解每个符号的来龙去脉然后去读SPIHT的链表机制对比两者在同一个测试图上的码流差异最后再碰EBCOTJPEG2000的核心那种基于码块和算术编码的更复杂体系。每一步都自己写代码验证别只看论文——我踩过的坑几乎全是论文一句话代码一整天的地方比如零树判定的边界、辅扫描的列表划分、逆变换的顺序这些只有亲手写才会真正记住。最后分享一个调试小技巧编码器和解码器分开写完之后先别急着测PSNR先用一个极小的测试图比如8x8只做1级分解把每一轮的符号序列打印出来人工核对主扫描和辅扫描的输出是否符合预期。小图上逻辑错误一目了然等小图对了再上大图能省掉大量盲目排查的时间。
阅读完成 · 觉得有帮助？
