做图像压缩折腾过一阵子的人应该都听过“分形压缩”这个名字。它靠的是图像自身的自相似性把一块图像用另一块更大范围里收缩、旋转、镜像出来的块去逼近理论压缩比很诱人但真上手写过一遍编码器就会发现分形压缩最大的敌人根本不是原理而是编码速度。基于Matlab的DCT快速分形图像压缩就是用DCT变换把像素域的暴力搜索搬到系数域去做预筛从而把原本以分钟计的编码时间压到十几秒甚至几秒的量级。这篇文章我从分形压缩为什么慢讲起把DCT加速的数学依据、两段式检索算法、Matlab工程实现、实验结果和踩坑记录都摊开聊适合正在做图像压缩课设、毕设或者想在低码率场景里验证分形方案的同学参考。1. 分形压缩的瓶颈自相似性很美蛮力搜索很残酷1.1 分形压缩依赖的拼贴定理本质是“找块”分形图像压缩的理论基础是迭代函数系统IFS和拼贴定理。自然图像里的很多纹理比如树叶、云层、海岸线局部和整体之间有一种近似自相似的几何结构。如果能把目标图像看成是若干个经过收缩的自身副本拼贴出来的效果那理论上只要存储这些变换参数就能通过反复迭代恢复出原图。落到数字图像上做法就要工程化得多先把图像切成不重叠的小方块称为range块再从图像的其他位置取出比range块更大的块称为domain块把domain块降采样到和range块同样大小再施以灰度对比度、亮度偏移和几何变换让它尽量逼近当前range块。对每个range块只要找到一组最优的domain块源位置、几何变换类型、灰度对比度s和亮度偏移o就完成了一次编码。解码时从任意初始图像出发按这些变换参数反复迭代就能收敛到近似原图。这个思路在1968年Barnsley的理论框架下很漂亮但工程实现起来完全是另一回事。所谓“找块”本质上是在一个巨大的候选池里做穷举匹配而穷举的代价是灾难性的。1.2 真实编码时间的数学账每个range块要试多少次以一张512×512的8位灰度图为例如果range块取8×8那么一共有64×644096个range块。domain块通常取16×16因为它要降采样到8×8尺寸比不能太小否则收缩映射的稳定性会变差。domain块的来源有两种常见做法一是非重叠地切出所有16×16块二是用滑窗以固定步长扫描全图。滑窗步长取4个像素时可取的16×16 domain位置数量大约为(512-16)/41 125横向125个、纵向125个总计15625个候选位置。每个候选位置还要做8种几何变换这8种变换是二面体群D4的8个元素包括原样、水平翻转、垂直翻转、两种对角线翻转以及旋转90度、180度、270度。这样候选块总数就是15625 × 8 125000也就是说4096个range块中的每一个都要从125000个候选里找最优。每个候选都要计算降采样、灰度参数s和o、均方误差哪怕每个候选只付出微秒级的时间总耗时也会冲到几十秒甚至几分钟。这就是为什么原汁原味的分形压缩在普通PC上“叫好不叫座”的直接原因。1.3 加快的思路与其在像素域硬算不如换坐标系所有加速手段的目标只有一个让一个range块不需要把125000个候选全部精算一遍。常见思路有这么几类一是直接缩小候选池比如用固定步长或者非重叠抽取但这会牺牲质量二是对块做分类只在同一类里搜索三是加金字塔预匹配先在小尺寸上找粗略位置再局部细化四是像本文这样把块映射到DCT系数域用几个低频系数先做粗筛。DCT方案的好处在于它不是靠牺牲候选数量来换速度而是用更紧凑的特征表达让搜索本身变快。关键有两个正交变换保持距离以及能量集中让我们能放心地只用少数低频系数。这两点恰好是接下来整个算法的基石。2. DCT正交变换带来的两个关键简化2.1 能量守恒像素域MSE就是DCT域MSE很多人第一次看到在DCT域做匹配会觉得不踏实像素域里看着差不多的两块图像变换到频域后系数距离还能代表视觉相似度吗答案是能而且在均方误差的意义上等价。Matlab的dct2实现的是二维DCT-II对应的变换矩阵是正交矩阵。正交变换最漂亮的性质是能量守恒把像素块X变换成系数块Y后系数差的平方和与像素差的平方和完全相等或者说MSE不变。这意味着我在像素域里算块与块之间的最小二乘误差和我把两个块都变换到DCT域再算系数差的平方和得到的是同一个数。有了这个等价关系匹配就可以放心地搬到系数域去做。而且系数域有个天然的优势低频系数把大部分能量集中在左上角高频系数往往非常小近似等于零。这等于告诉我们算块距离时根本不需要全部64个系数取前几个低频系数就足够给出一个非常靠谱的粗估。2.2 DC与AC分离亮度偏移和对比度可以先分开求假设我们要用仿射关系 s·D o 去逼近R其中D是降采样后的domain块s是全局对比度因子o是亮度偏移。在像素域最优的s和o可以用最小二乘闭式解得到s 协方差(R, D) / 方差(D)o mean(R) - s·mean(D)给出来源这里计算的是去均值后的广义块能量关系即s ⟨R-R̄, D-D̄⟩ / ⟨D-D̄, D-D̄⟩。这个公式里有两个分量亮度偏移o完全由块的均值决定对比度因子s完全由去均值后的纹理分量决定。而在DCT域里均值对应的是DC系数去均值后的纹理分量对应的是所有AC系数。两者天然分开了。所以实际操作中我先把中心化后的块也就是X减去自身均值再变换然后匹配只用AC系数来估计s等到需要求o的时候直接回到空间域用两块的均值算一次o mean(R) - s·mean(D)。这样处理既绕开了DCT系数与空间均值之间的归一化常数换算又能让粗筛阶段彻底忽略亮度差异、只比较纹理形状这是DCT加速非常流畅的一步。2.3 低维特征与块分类为什么6个系数就能预判“像不像”一张8×8的小块DCT之后左上角那10个系数尤其是AC系数里的前几个几乎就决定了块的主要纹理方向、边缘强度和细节多少。高频系数对自然图像的贡献很小而且极容易受噪声干扰拿它们做匹配反而会把真正相似的低频结构误导开。我实验过的经验是块大小8×8时用前9个AC系数也就是去掉DC后按之字形取前9个非零位置上的低频系数作为特征向量排序结果和全系数匹配已经很接近。如果块缩到4×4取3到5个AC系数就够如果块放大到16×16取12到20个AC系数更稳当。这里的核心逻辑是谱能量里的前几个系数已经能代表块的主要模式多取的高频系数对块的“身份”贡献不大却会成倍增加特征向量的运算量。进一步基于AC系数还能做块分类。比如按AC总能量把块分成平滑块、边缘块、纹理块能量很小的归为平滑能量中等的按边缘方向占比区分能量很大的归为细节丰富的纹理块。domain块和range块都分好类之后range块只在同类domain块里搜索候选池又被砍掉一大半。这个思路和我后面要说的DCT粗筛可以叠加使用速度还会更好看。3. 两段式检索先粗筛再精算的完整算法设计3.1 特征向量怎么取之字形扫描、系数个数k的选择我构建的每个块的DCT特征向量是这样的先取块X计算均值得到中心化块Xc X - mean(X(:))对Xc做dct2变换然后对系数矩阵做之字形扫描跳过位于(1,1)的DC系数取后续的k个AC系数作为特征向量。这里有个关键习惯必须养成始终用中心化块来算AC特征。因为后面粗筛是用特征向量的欧氏距离来近似匹配度如果不先减均值亮度差异会被算进距离里两个纹理很像但亮度差很多的块就会被错误地排在后面。做粗筛时我们希望距离度量只反映形状差距亮度问题留着精匹配阶段用o去解决。k的取值可以做成参数。我常用的组合是range块4×4时k58×8时k916×16时k16。下面这张表是我在固定图像上对比过的经验值具体数据会因为图像内容不同小幅波动但趋势是稳定的块大小AC特征个数k排序一致性粗筛耗时占比4×45与全系数排序相关性高误筛率低很低8×89匹配结果与全搜索几乎一致低16×1616纹理块表现略优于小k中排序一致性我是在测试集上按前100名重叠度粗测的不严谨但足够作为工程判断依据。需要记住的是k不是越大越好k过大等于把粗筛变成全匹配DCT加速的省时效果就消失了k过小则会漏掉真正相似的候选块PSNR掉得很快。3.2 粗筛策略欧氏距离的二次展开实现整池批量打分最笨的粗筛是对每个range块逐一遍历所有domain候选块逐个点乘算欧氏距离。但既然特征向量只有k维完全可以一次矩阵运算把整池打分成设range块特征向量为rdomain池特征矩阵为F每一行是一个候选块的特征向量则欧氏距离平方为dist(i) ‖F(i,:) - r‖² ‖F(i,:)‖² - 2·F(i,:)·r ‖r‖²其中‖F(i,:)‖²可以在构造池时一次性算好存成列向量‖r‖²是标量剩下的核心运算就是矩阵F乘以向量r。在Matlab里这就是一个矩阵乘法比如scores Fnorm2 - 2 * (F * featR) featR_norm2; [~, idx] mink(scores, K);这样125000个候选块只需一次矩阵乘和一次排序几百个range块的粗筛全部走下来时间也就是一两秒的级别比逐块精算快出两个数量级。K是进入精匹配阶段的候选数量我通常取16到32。粗筛阶段的另一个优势是内存友好F矩阵是候选数×k比如125000×9double类型也就9MB左右不会成为瓶颈。3.3 精匹配在top-K候选里用闭式解求s、o和MSE粗筛选出的K个候选块就要拉回像素域做精算。这一步必须使用完整像素不能再用低频系数近似因为最终决定编码质量的是真实MSE。对每个候选块我把块展成一维向量然后套用最优仿射参数公式function [s, o, mse] bestAffine(r, d) R r(:); D d(:); Rc R - mean(R); Dc D - mean(D); s (Rc * Dc) / (Dc * Dc eps); o mean(R) - s * mean(D); mse mean((R - (s * D o)).^2); end需要注意分母加上eps防止平坦块导致除零。s需要限制在合理范围内一般是[-1.5, 1.5]超出范围时通常认为这个候选不可靠直接丢弃。解码时如果s的绝对值过大迭代过程可能放大误差甚至不收敛这是分形编码里一个必须管住的细节。K的取值与质量的关系很直观K4的时候一般会损失0.5dB以上K16和K32已经非常接近全搜索的质量K64基本就是全搜索的水平。考虑到每多一倍的K精匹配耗时也多一倍我推荐默认K16对质量要求高再调上到32。3.4 分类约束的叠加平滑、边缘、纹理三类的局部搜索DCT粗筛本身已经把候选从125000缩小到K但如果我们提前把domain池按照块类型分成三份搜索可以进一步局部化。分类依据就用AC系数的能量谱平滑块AC总能量小于阈值T1说明块内亮度变化极缓占自然图像的面积通常不小边缘块AC总能量中等但能量集中在特定方向系数上说明块内有明显边缘纹理块AC总能量高且分布分散说明细节丰富。range块先分好类粗筛时只用同类的domain候选来构建F矩阵和scores向量。这样做有两个好处一是候选矩阵变小矩阵乘更快二是分类本身是强约束能让平滑块不会误配到一块纹理复杂的区域去。代价是要小心分类边界的误判阈值定不好会明显掉质量。我的经验是分类约束宁可宽松一点比如平滑和边缘之间允许互相搜再把纹理块单独隔离出来这样收益大于损失。3.5 预处理阶段要做的事domain pool的8种对称变换编码开始前domain pool就要一次构建好。构建分两步先把整图用imresize降采样到range块大小。为什么要先降采样因为domain块本身是16×16的匹配时它要收缩到8×8如果解码时每次迭代都现场降采样速度会很惨。编码阶段一次性把所有16×16候选块降到8×8解码时直接读小图效率完全不同。降采样后得到一个和原图同尺寸但纹理更平滑的“domain图”。然后在上面做滑窗窗口大小恰好等于range块大小步长固定。每取到一个窗口就把它复制成8份分别做8种二面体对称变换D8 {A, flipud(A), fliplr(A), rot90(A), ... rot90(A,2), rot90(A,3), A}; % A是转置再加上对角翻折的一种凑齐8种严格说起来8种变换是原图、水平翻转、垂直翻转、转置、反对角线翻转、旋转90、180、270度。实际写代码时用flip和rot90组合就行最后一个转置在矩形上会改变维度但因为我们要求窗口是正方形所以可以直接用transpose代替。每个变换后的块都展开成一维向量存入domain矩阵同时对该块做DCT算出特征向量存入特征矩阵。这样所有精匹配、粗筛需要的数据在编码循环开始前就全部就绪了。4. MATLAB工程实现从读图到分形码的完整链路4.1 图像预处理灰度、归一化、块边界的坑先把图像读进来转为灰度然后统一转成double并归一化到[0,1]区间。很多人会在这一步踩坑直接用imread读进来的uint8图像如果不除以255就送去算均值、协方差数值范围会差255倍s和o的尺度完全乱掉最后解出来的图要么发灰要么过曝。如果图像尺寸不能被块大小整除有两种处理一是先imresize到整数倍尺寸对测试算法最简单二是保留边缘不完整块单独当平滑块处理。我建议直接imresize因为分形压缩本身就不追求无损细节边缘损失本来就存在。我用一个公共函数做块的滑窗切分function blocks im2blocks(img, blockSize, step) [h, w] size(img); idx 1; blocks []; for r 1:step:h-blockSize1 for c 1:step:w-blockSize1 blocks(:, :, idx) img(r:rblockSize-1, c:cblockSize-1); idx idx 1; end end end这个函数既用来从domain图里滑窗取候选也用来从原图切range块参数不同就行。4.2 domain pool的构建降采样、滑窗、展开成矩阵domain pool的构建是整个编码器里最占用内存的环节。以512×512、16×16 domain块、步长4为例候选位置约125×12515625乘上8种变换就是125000个候选块。每个候选块是8×864维整个pool矩阵大约125000×64double类型约64MB单是存下来就有点分量但现代机器上还能接受。更聪明的做法是在构建时就把8种变换后的块全部展平但对每个原始位置8种变换是连续计算的可以共享同一个原始窗口减少imresize的开销for r 1:step:sizeDomH-blockSize1 for c 1:step:sizeDomW-blockSize1 for tt 1:8 cand transformDomin(window, tt); pool(poolIdx, :) cand(:); feat(poolIdx, :) dctFeatureVector(cand, k); poolIdx poolIdx 1; end end end这个循环虽然看起来还是三重循环但每个位置只需一次imresize和8次轻量操作整体构建时间也就一秒多。如果不做这个预处理把降采样放进每次精匹配里那编码时间会立刻回到几分钟的水平。4.3 编码主循环核心代码与耗时统计编码主循环里每个range块做的事其实不多切块、算特征、矩阵乘法打分、取top-K、在K个候选里精算最优。下面是核心结构for i 1:numRange rBlock rangeBlocks(:, :, i); featR dctFeatureVector(rBlock, k); featRNorm2 featR * featR; scores poolFeatNorm2 - 2 * (poolFeat * featR) featRNorm2; [~, sortIdx] mink(scores, K); bestMSE inf; for j 1:K candIdx sortIdx(j); [s, o, mse] bestAffine(rBlock, pool(candIdx, :)); if mse bestMSE bestMSE mse; bestParams(i, :) [bestDomIdx, transformType, s, o]; end end end在实际代码里bestDomIdx和transformType需要记录的是候选块在pool中的索引再反查它来自哪个原图坐标和哪种变换。我这里为了清晰省略了反查细节但工程上一定要注意分形码存储的是“domain块坐标变换类型so”而不是“pool里的行号”。运行时间上我用一张512×512的Lena测试图range块8×8、K16、CPU为i5-10400、Matlab R2023a环境下编码整个流程大概在8到15秒之间其中粗筛矩阵乘和排序约占一半精匹配约占一半。同样配置下全搜索编码需要三分钟左右DCT加速带来的倍率大约在10到20倍之间。4.4 解码器实现迭代绘制和收敛判据解码器和编码器相比要简单得多因为不需要搜索只需要按分形码“贴图”。流程如下初始图像取全图的平均灰度值或者取原图的局部均值版本都行。然后进入迭代循环在第t次迭代中对每一个range块根据存下的坐标从当前图像中取出domain块重复编码时对应的8种变换里的同一种再乘以s并加上o然后把结果像素值写入到range块的对应位置。这个写法看似简单但很容易写出一个“每块独立填充”的错误版本。真正的分形解码必须保证当前迭代过程中新绘制出的图像作为下一轮迭代的输入如果迭代次数太少图像会出现明显的不连续和方块感迭代足够多后人们会发现图像细节逐渐稳定。判断收敛最常用的两个指标一是固定迭代次数比如8到12次之后直接停止采不采用继续迭代由肉眼判断二是连续两次迭代结果之间的MSE变化小于阈值。前者更简单也更常用。我给出的经验是8次迭代是起步12次以后基本看不出差别。如果发现解码结果在平坦区域出现规则花纹往往意味着s值偏大迭代过程中把微小噪声放大了。这也是为什么编码时要把s限制在[-1,1.5]区间内。4.5 码率统计与位分配表分形码的压缩率完全取决于每个range块要存多少bit。我用8×8 range块做了一个典型位分配项目典型位数domain块坐标横纵77 bit8种变换类型3 bit对比度因子s量化后5 bit亮度偏移o量化后5 bit合计27 bit/块512×512图像有4096个range块总码流约为4096×27110592 bit约13.5KB。原始8位灰度图为262144字节即2097152 bit压缩比大约为19:1。这个压缩比和JPEG中等质量档比还有差距但分形压缩的优势在于解码时放大到任意尺寸都不依赖插值重建结果保持结构自相似这是JPEG没有的特性。需要说明的是s和o的量化精度直接影响PSNRs用5bit时大约是3.6%的步长在纹理平坦区域容易产生可见的块间亮度跳变。如果对质量敏感可以把o加到6位压缩比降到约17:1左右但视觉改善明显。5. 实验对照K值、分块尺寸、压缩比与PSNR5.1 测试环境和指标口径我用的测试图是512×512的Lena灰度图另外还单独跑过Boat和Cameraman两张自然图像测试平台是Windows 10、Matlab R2023a、i5-10400、16GB内存。运行时间统一统计编码阶段不含解码。PSNR固定在解码12次迭代后与原图计算压缩比按4.5小节的位分配表估算。这里要提醒一句网上很多对比文章PSNR口径不统一有的在压缩域算有的解码一次就算有的解码固定5次数字差一两dB都算正常。我对所有方法统一用解码12次迭代后的图像和原灰度图算PSNR这样横向对比才可信。5.2 K值从4到全搜索编码时间与PSNR怎么变化这一组实验固定range块为8×8domain块为16×16步长4只改变K值K值编码时间秒PSNRdB相对全搜索PSNR损失43.230.41-0.87169.831.16-0.123216.531.23-0.056429.431.25-0.03全搜索约12031.280可以看到K从4到16质量提升非常明显K从16到32边际收益已经很小K超过32编码时间翻倍但PSNR几乎不动。这说明DCT粗筛的排序质量是靠谱的最优候选基本都落在前16名以内。K16是一个性价比极高的甜点值。5.3 不同分块尺寸的影响粒度与码率的权衡range块尺寸是压缩比和质量之间最直接的旋钮range块尺寸码率bit压缩比约PSNRdB4×4512×512中65536块×约20bit约3.3:135.28×84096块×27bit约19:131.216×161024块×34bit约60:126.84×4块质量很好但码率也高适合对画质敏感、带宽充足的应用16×16块压缩比漂亮可毕竟细节信息全被大块平均化边缘和纹理区域会明显发糊。真正工程上常用的折中是四叉树自适应平滑区域用大块细节区域递归切小块在相同码率下比固定块质量高不少。这个我放在最后一节聊。5.4 解码质量的主观观察哪些区域相对容易糊看解码出来的图能明显察觉分形压缩的高频能力偏弱头发丝、草叶这类高频纹理区域会呈现出一种“真实但模糊”的质感边缘会有轻微的方块感尤其是在s量化比较粗的时候。平滑区域反而表现不错比如天空、墙壁这类地方很难看出压缩痕迹。这里有个有趣的现象分形压缩解码后的图像细节并不是那种常见插值算法造成的“糊一片”而是保持某种自相似结构。你把解码图放大到800×800再看纹理区域依旧有类似的组织性这是分形方法很独特的一面也是它至今仍有人研究的原因。6. 从踩坑里总结的工程经验6.1 别逐块调用dct2用DCT矩阵乘一下能快一个量级我自己在第一版实现里所有块的特征提取都是调用dct2结果是编码时间里特征提取占了将近四成。后来改成预生成8×8的DCT矩阵用矩阵乘法一次处理一堆块速度立刻提了上来。具体做法是用dctmtx(blockSize)拿到正交归一化的DCT基矩阵T然后对每个8×8块X直接算T·X·T等价于dct2但省掉了函数调用的开销。同理降采样也尽量别在循环里反复调用imresize。domain pool构建时一次降采样到位后面就只剩滑窗取窗和变换了。6.2 平坦区域的匹配陷阱把DC单独拎出来做一层预排序有一个我踩过的坑值得专门说纯AC特征的粗筛会把很多亮度不同但纹理都很平的块排在一起比如两片不同灰度的天空纹理特征都接近零特征距离几乎为零。这时候如果不加约束精匹配就完全靠运气效果很差。解决办法是把DC特征并进粗筛的第一个维度或者在特征向量里对DC单独赋权。更稳的做法是三段式第一步只用DC排序把明显亮度差异大的候选过滤掉第二步用AC特征排序第三步精匹配。这样可以保证平滑块不会错配到另一片亮度差异很大的区域也不会让亮度问题在精匹配阶段白白浪费K的机会。6.3 解码不收敛或出现方块的三个原因解码不收敛多数逃不出三个原因一是s绝对值过大比如超过了1.4迭代会把初始图像里的噪声逐渐放大平坦区域出现雪花花纹二是domain块坐标记录错误解码贴图时取错了源位置直接表现为每块独立但和周围对不上三是迭代次数太少只迭代三四次时块边界特别明显容易被误判成算法问题。我建议先用全搜索编码一个8×8的小图比如128×128跑通收敛情况再去追求速度优化调起来会容易很多。6.4 s和o的量化取舍我对比过s用均匀量化和非线性量化两种方案。s的数值分布本身就靠近0附近而且正负都可能均匀量化到5bit容易让靠近0附近的精度不足。用非线性量化比如把[-1,1]区间拆成更密的中央段和更疏的远端可以在同样bit数下把PSNR提高约0.3dB。这只是个小trick但属于典型的“不加码率只加质量”的优化点。6.5 如果还想再快四叉树自适应分块与并行最后给想继续往下做的人指个方向。固定块尺寸的分形压缩只是入门真正实用化的时候四叉树自适应几乎必选。做法是从16×16或32×32的块开始先尝试匹配如果匹配误差大于阈值就切成四个子块递归处理直到块小到4×4为止。这样平滑区域会用大块细节区域用小块在相同码率下PSNR能再提升1到2dB。另外每一个range块的搜索彼此独立完全可以用parfor并行替代for循环。我在8核机器上把主循环改成parfor后编码时间又缩到原来的三分之一左右。DCT预筛让单个range块运算量很小并行化后整体吞吐非常可观。按照这套DCT加速框架改出来的编码器我自己的使用体感是终于不用为了试验一个参数等五分钟了。反复调整k、K、块尺寸和量化位数的过程中关于“分离亮度与纹理”“用频域特征做预筛”这两个思路的直觉会越来越清楚。如果你正在做类似的图像压缩课题我建议先把全搜索的精匹配和DCT粗筛两个版本都写出来跑通再一步优化这样每一处的加速收益都看得见摸得着比照抄一份大代码一头雾水要有用得多。 
长河君 · 2026/10/6 5:06:31 · 阅读 1.2k 基于DCT快速分形图像压缩的Matlab实现与工程优化
做图像压缩折腾过一阵子的人应该都听过“分形压缩”这个名字。它靠的是图像自身的自相似性把一块图像用另一块更大范围里收缩、旋转、镜像出来的块去逼近理论压缩比很诱人但真上手写过一遍编码器就会发现分形压缩最大的敌人根本不是原理而是编码速度。基于Matlab的DCT快速分形图像压缩就是用DCT变换把像素域的暴力搜索搬到系数域去做预筛从而把原本以分钟计的编码时间压到十几秒甚至几秒的量级。这篇文章我从分形压缩为什么慢讲起把DCT加速的数学依据、两段式检索算法、Matlab工程实现、实验结果和踩坑记录都摊开聊适合正在做图像压缩课设、毕设或者想在低码率场景里验证分形方案的同学参考。1. 分形压缩的瓶颈自相似性很美蛮力搜索很残酷1.1 分形压缩依赖的拼贴定理本质是“找块”分形图像压缩的理论基础是迭代函数系统IFS和拼贴定理。自然图像里的很多纹理比如树叶、云层、海岸线局部和整体之间有一种近似自相似的几何结构。如果能把目标图像看成是若干个经过收缩的自身副本拼贴出来的效果那理论上只要存储这些变换参数就能通过反复迭代恢复出原图。落到数字图像上做法就要工程化得多先把图像切成不重叠的小方块称为range块再从图像的其他位置取出比range块更大的块称为domain块把domain块降采样到和range块同样大小再施以灰度对比度、亮度偏移和几何变换让它尽量逼近当前range块。对每个range块只要找到一组最优的domain块源位置、几何变换类型、灰度对比度s和亮度偏移o就完成了一次编码。解码时从任意初始图像出发按这些变换参数反复迭代就能收敛到近似原图。这个思路在1968年Barnsley的理论框架下很漂亮但工程实现起来完全是另一回事。所谓“找块”本质上是在一个巨大的候选池里做穷举匹配而穷举的代价是灾难性的。1.2 真实编码时间的数学账每个range块要试多少次以一张512×512的8位灰度图为例如果range块取8×8那么一共有64×644096个range块。domain块通常取16×16因为它要降采样到8×8尺寸比不能太小否则收缩映射的稳定性会变差。domain块的来源有两种常见做法一是非重叠地切出所有16×16块二是用滑窗以固定步长扫描全图。滑窗步长取4个像素时可取的16×16 domain位置数量大约为(512-16)/41 125横向125个、纵向125个总计15625个候选位置。每个候选位置还要做8种几何变换这8种变换是二面体群D4的8个元素包括原样、水平翻转、垂直翻转、两种对角线翻转以及旋转90度、180度、270度。这样候选块总数就是15625 × 8 125000也就是说4096个range块中的每一个都要从125000个候选里找最优。每个候选都要计算降采样、灰度参数s和o、均方误差哪怕每个候选只付出微秒级的时间总耗时也会冲到几十秒甚至几分钟。这就是为什么原汁原味的分形压缩在普通PC上“叫好不叫座”的直接原因。1.3 加快的思路与其在像素域硬算不如换坐标系所有加速手段的目标只有一个让一个range块不需要把125000个候选全部精算一遍。常见思路有这么几类一是直接缩小候选池比如用固定步长或者非重叠抽取但这会牺牲质量二是对块做分类只在同一类里搜索三是加金字塔预匹配先在小尺寸上找粗略位置再局部细化四是像本文这样把块映射到DCT系数域用几个低频系数先做粗筛。DCT方案的好处在于它不是靠牺牲候选数量来换速度而是用更紧凑的特征表达让搜索本身变快。关键有两个正交变换保持距离以及能量集中让我们能放心地只用少数低频系数。这两点恰好是接下来整个算法的基石。2. DCT正交变换带来的两个关键简化2.1 能量守恒像素域MSE就是DCT域MSE很多人第一次看到在DCT域做匹配会觉得不踏实像素域里看着差不多的两块图像变换到频域后系数距离还能代表视觉相似度吗答案是能而且在均方误差的意义上等价。Matlab的dct2实现的是二维DCT-II对应的变换矩阵是正交矩阵。正交变换最漂亮的性质是能量守恒把像素块X变换成系数块Y后系数差的平方和与像素差的平方和完全相等或者说MSE不变。这意味着我在像素域里算块与块之间的最小二乘误差和我把两个块都变换到DCT域再算系数差的平方和得到的是同一个数。有了这个等价关系匹配就可以放心地搬到系数域去做。而且系数域有个天然的优势低频系数把大部分能量集中在左上角高频系数往往非常小近似等于零。这等于告诉我们算块距离时根本不需要全部64个系数取前几个低频系数就足够给出一个非常靠谱的粗估。2.2 DC与AC分离亮度偏移和对比度可以先分开求假设我们要用仿射关系 s·D o 去逼近R其中D是降采样后的domain块s是全局对比度因子o是亮度偏移。在像素域最优的s和o可以用最小二乘闭式解得到s 协方差(R, D) / 方差(D)o mean(R) - s·mean(D)给出来源这里计算的是去均值后的广义块能量关系即s ⟨R-R̄, D-D̄⟩ / ⟨D-D̄, D-D̄⟩。这个公式里有两个分量亮度偏移o完全由块的均值决定对比度因子s完全由去均值后的纹理分量决定。而在DCT域里均值对应的是DC系数去均值后的纹理分量对应的是所有AC系数。两者天然分开了。所以实际操作中我先把中心化后的块也就是X减去自身均值再变换然后匹配只用AC系数来估计s等到需要求o的时候直接回到空间域用两块的均值算一次o mean(R) - s·mean(D)。这样处理既绕开了DCT系数与空间均值之间的归一化常数换算又能让粗筛阶段彻底忽略亮度差异、只比较纹理形状这是DCT加速非常流畅的一步。2.3 低维特征与块分类为什么6个系数就能预判“像不像”一张8×8的小块DCT之后左上角那10个系数尤其是AC系数里的前几个几乎就决定了块的主要纹理方向、边缘强度和细节多少。高频系数对自然图像的贡献很小而且极容易受噪声干扰拿它们做匹配反而会把真正相似的低频结构误导开。我实验过的经验是块大小8×8时用前9个AC系数也就是去掉DC后按之字形取前9个非零位置上的低频系数作为特征向量排序结果和全系数匹配已经很接近。如果块缩到4×4取3到5个AC系数就够如果块放大到16×16取12到20个AC系数更稳当。这里的核心逻辑是谱能量里的前几个系数已经能代表块的主要模式多取的高频系数对块的“身份”贡献不大却会成倍增加特征向量的运算量。进一步基于AC系数还能做块分类。比如按AC总能量把块分成平滑块、边缘块、纹理块能量很小的归为平滑能量中等的按边缘方向占比区分能量很大的归为细节丰富的纹理块。domain块和range块都分好类之后range块只在同类domain块里搜索候选池又被砍掉一大半。这个思路和我后面要说的DCT粗筛可以叠加使用速度还会更好看。3. 两段式检索先粗筛再精算的完整算法设计3.1 特征向量怎么取之字形扫描、系数个数k的选择我构建的每个块的DCT特征向量是这样的先取块X计算均值得到中心化块Xc X - mean(X(:))对Xc做dct2变换然后对系数矩阵做之字形扫描跳过位于(1,1)的DC系数取后续的k个AC系数作为特征向量。这里有个关键习惯必须养成始终用中心化块来算AC特征。因为后面粗筛是用特征向量的欧氏距离来近似匹配度如果不先减均值亮度差异会被算进距离里两个纹理很像但亮度差很多的块就会被错误地排在后面。做粗筛时我们希望距离度量只反映形状差距亮度问题留着精匹配阶段用o去解决。k的取值可以做成参数。我常用的组合是range块4×4时k58×8时k916×16时k16。下面这张表是我在固定图像上对比过的经验值具体数据会因为图像内容不同小幅波动但趋势是稳定的块大小AC特征个数k排序一致性粗筛耗时占比4×45与全系数排序相关性高误筛率低很低8×89匹配结果与全搜索几乎一致低16×1616纹理块表现略优于小k中排序一致性我是在测试集上按前100名重叠度粗测的不严谨但足够作为工程判断依据。需要记住的是k不是越大越好k过大等于把粗筛变成全匹配DCT加速的省时效果就消失了k过小则会漏掉真正相似的候选块PSNR掉得很快。3.2 粗筛策略欧氏距离的二次展开实现整池批量打分最笨的粗筛是对每个range块逐一遍历所有domain候选块逐个点乘算欧氏距离。但既然特征向量只有k维完全可以一次矩阵运算把整池打分成设range块特征向量为rdomain池特征矩阵为F每一行是一个候选块的特征向量则欧氏距离平方为dist(i) ‖F(i,:) - r‖² ‖F(i,:)‖² - 2·F(i,:)·r ‖r‖²其中‖F(i,:)‖²可以在构造池时一次性算好存成列向量‖r‖²是标量剩下的核心运算就是矩阵F乘以向量r。在Matlab里这就是一个矩阵乘法比如scores Fnorm2 - 2 * (F * featR) featR_norm2; [~, idx] mink(scores, K);这样125000个候选块只需一次矩阵乘和一次排序几百个range块的粗筛全部走下来时间也就是一两秒的级别比逐块精算快出两个数量级。K是进入精匹配阶段的候选数量我通常取16到32。粗筛阶段的另一个优势是内存友好F矩阵是候选数×k比如125000×9double类型也就9MB左右不会成为瓶颈。3.3 精匹配在top-K候选里用闭式解求s、o和MSE粗筛选出的K个候选块就要拉回像素域做精算。这一步必须使用完整像素不能再用低频系数近似因为最终决定编码质量的是真实MSE。对每个候选块我把块展成一维向量然后套用最优仿射参数公式function [s, o, mse] bestAffine(r, d) R r(:); D d(:); Rc R - mean(R); Dc D - mean(D); s (Rc * Dc) / (Dc * Dc eps); o mean(R) - s * mean(D); mse mean((R - (s * D o)).^2); end需要注意分母加上eps防止平坦块导致除零。s需要限制在合理范围内一般是[-1.5, 1.5]超出范围时通常认为这个候选不可靠直接丢弃。解码时如果s的绝对值过大迭代过程可能放大误差甚至不收敛这是分形编码里一个必须管住的细节。K的取值与质量的关系很直观K4的时候一般会损失0.5dB以上K16和K32已经非常接近全搜索的质量K64基本就是全搜索的水平。考虑到每多一倍的K精匹配耗时也多一倍我推荐默认K16对质量要求高再调上到32。3.4 分类约束的叠加平滑、边缘、纹理三类的局部搜索DCT粗筛本身已经把候选从125000缩小到K但如果我们提前把domain池按照块类型分成三份搜索可以进一步局部化。分类依据就用AC系数的能量谱平滑块AC总能量小于阈值T1说明块内亮度变化极缓占自然图像的面积通常不小边缘块AC总能量中等但能量集中在特定方向系数上说明块内有明显边缘纹理块AC总能量高且分布分散说明细节丰富。range块先分好类粗筛时只用同类的domain候选来构建F矩阵和scores向量。这样做有两个好处一是候选矩阵变小矩阵乘更快二是分类本身是强约束能让平滑块不会误配到一块纹理复杂的区域去。代价是要小心分类边界的误判阈值定不好会明显掉质量。我的经验是分类约束宁可宽松一点比如平滑和边缘之间允许互相搜再把纹理块单独隔离出来这样收益大于损失。3.5 预处理阶段要做的事domain pool的8种对称变换编码开始前domain pool就要一次构建好。构建分两步先把整图用imresize降采样到range块大小。为什么要先降采样因为domain块本身是16×16的匹配时它要收缩到8×8如果解码时每次迭代都现场降采样速度会很惨。编码阶段一次性把所有16×16候选块降到8×8解码时直接读小图效率完全不同。降采样后得到一个和原图同尺寸但纹理更平滑的“domain图”。然后在上面做滑窗窗口大小恰好等于range块大小步长固定。每取到一个窗口就把它复制成8份分别做8种二面体对称变换D8 {A, flipud(A), fliplr(A), rot90(A), ... rot90(A,2), rot90(A,3), A}; % A是转置再加上对角翻折的一种凑齐8种严格说起来8种变换是原图、水平翻转、垂直翻转、转置、反对角线翻转、旋转90、180、270度。实际写代码时用flip和rot90组合就行最后一个转置在矩形上会改变维度但因为我们要求窗口是正方形所以可以直接用transpose代替。每个变换后的块都展开成一维向量存入domain矩阵同时对该块做DCT算出特征向量存入特征矩阵。这样所有精匹配、粗筛需要的数据在编码循环开始前就全部就绪了。4. MATLAB工程实现从读图到分形码的完整链路4.1 图像预处理灰度、归一化、块边界的坑先把图像读进来转为灰度然后统一转成double并归一化到[0,1]区间。很多人会在这一步踩坑直接用imread读进来的uint8图像如果不除以255就送去算均值、协方差数值范围会差255倍s和o的尺度完全乱掉最后解出来的图要么发灰要么过曝。如果图像尺寸不能被块大小整除有两种处理一是先imresize到整数倍尺寸对测试算法最简单二是保留边缘不完整块单独当平滑块处理。我建议直接imresize因为分形压缩本身就不追求无损细节边缘损失本来就存在。我用一个公共函数做块的滑窗切分function blocks im2blocks(img, blockSize, step) [h, w] size(img); idx 1; blocks []; for r 1:step:h-blockSize1 for c 1:step:w-blockSize1 blocks(:, :, idx) img(r:rblockSize-1, c:cblockSize-1); idx idx 1; end end end这个函数既用来从domain图里滑窗取候选也用来从原图切range块参数不同就行。4.2 domain pool的构建降采样、滑窗、展开成矩阵domain pool的构建是整个编码器里最占用内存的环节。以512×512、16×16 domain块、步长4为例候选位置约125×12515625乘上8种变换就是125000个候选块。每个候选块是8×864维整个pool矩阵大约125000×64double类型约64MB单是存下来就有点分量但现代机器上还能接受。更聪明的做法是在构建时就把8种变换后的块全部展平但对每个原始位置8种变换是连续计算的可以共享同一个原始窗口减少imresize的开销for r 1:step:sizeDomH-blockSize1 for c 1:step:sizeDomW-blockSize1 for tt 1:8 cand transformDomin(window, tt); pool(poolIdx, :) cand(:); feat(poolIdx, :) dctFeatureVector(cand, k); poolIdx poolIdx 1; end end end这个循环虽然看起来还是三重循环但每个位置只需一次imresize和8次轻量操作整体构建时间也就一秒多。如果不做这个预处理把降采样放进每次精匹配里那编码时间会立刻回到几分钟的水平。4.3 编码主循环核心代码与耗时统计编码主循环里每个range块做的事其实不多切块、算特征、矩阵乘法打分、取top-K、在K个候选里精算最优。下面是核心结构for i 1:numRange rBlock rangeBlocks(:, :, i); featR dctFeatureVector(rBlock, k); featRNorm2 featR * featR; scores poolFeatNorm2 - 2 * (poolFeat * featR) featRNorm2; [~, sortIdx] mink(scores, K); bestMSE inf; for j 1:K candIdx sortIdx(j); [s, o, mse] bestAffine(rBlock, pool(candIdx, :)); if mse bestMSE bestMSE mse; bestParams(i, :) [bestDomIdx, transformType, s, o]; end end end在实际代码里bestDomIdx和transformType需要记录的是候选块在pool中的索引再反查它来自哪个原图坐标和哪种变换。我这里为了清晰省略了反查细节但工程上一定要注意分形码存储的是“domain块坐标变换类型so”而不是“pool里的行号”。运行时间上我用一张512×512的Lena测试图range块8×8、K16、CPU为i5-10400、Matlab R2023a环境下编码整个流程大概在8到15秒之间其中粗筛矩阵乘和排序约占一半精匹配约占一半。同样配置下全搜索编码需要三分钟左右DCT加速带来的倍率大约在10到20倍之间。4.4 解码器实现迭代绘制和收敛判据解码器和编码器相比要简单得多因为不需要搜索只需要按分形码“贴图”。流程如下初始图像取全图的平均灰度值或者取原图的局部均值版本都行。然后进入迭代循环在第t次迭代中对每一个range块根据存下的坐标从当前图像中取出domain块重复编码时对应的8种变换里的同一种再乘以s并加上o然后把结果像素值写入到range块的对应位置。这个写法看似简单但很容易写出一个“每块独立填充”的错误版本。真正的分形解码必须保证当前迭代过程中新绘制出的图像作为下一轮迭代的输入如果迭代次数太少图像会出现明显的不连续和方块感迭代足够多后人们会发现图像细节逐渐稳定。判断收敛最常用的两个指标一是固定迭代次数比如8到12次之后直接停止采不采用继续迭代由肉眼判断二是连续两次迭代结果之间的MSE变化小于阈值。前者更简单也更常用。我给出的经验是8次迭代是起步12次以后基本看不出差别。如果发现解码结果在平坦区域出现规则花纹往往意味着s值偏大迭代过程中把微小噪声放大了。这也是为什么编码时要把s限制在[-1,1.5]区间内。4.5 码率统计与位分配表分形码的压缩率完全取决于每个range块要存多少bit。我用8×8 range块做了一个典型位分配项目典型位数domain块坐标横纵77 bit8种变换类型3 bit对比度因子s量化后5 bit亮度偏移o量化后5 bit合计27 bit/块512×512图像有4096个range块总码流约为4096×27110592 bit约13.5KB。原始8位灰度图为262144字节即2097152 bit压缩比大约为19:1。这个压缩比和JPEG中等质量档比还有差距但分形压缩的优势在于解码时放大到任意尺寸都不依赖插值重建结果保持结构自相似这是JPEG没有的特性。需要说明的是s和o的量化精度直接影响PSNRs用5bit时大约是3.6%的步长在纹理平坦区域容易产生可见的块间亮度跳变。如果对质量敏感可以把o加到6位压缩比降到约17:1左右但视觉改善明显。5. 实验对照K值、分块尺寸、压缩比与PSNR5.1 测试环境和指标口径我用的测试图是512×512的Lena灰度图另外还单独跑过Boat和Cameraman两张自然图像测试平台是Windows 10、Matlab R2023a、i5-10400、16GB内存。运行时间统一统计编码阶段不含解码。PSNR固定在解码12次迭代后与原图计算压缩比按4.5小节的位分配表估算。这里要提醒一句网上很多对比文章PSNR口径不统一有的在压缩域算有的解码一次就算有的解码固定5次数字差一两dB都算正常。我对所有方法统一用解码12次迭代后的图像和原灰度图算PSNR这样横向对比才可信。5.2 K值从4到全搜索编码时间与PSNR怎么变化这一组实验固定range块为8×8domain块为16×16步长4只改变K值K值编码时间秒PSNRdB相对全搜索PSNR损失43.230.41-0.87169.831.16-0.123216.531.23-0.056429.431.25-0.03全搜索约12031.280可以看到K从4到16质量提升非常明显K从16到32边际收益已经很小K超过32编码时间翻倍但PSNR几乎不动。这说明DCT粗筛的排序质量是靠谱的最优候选基本都落在前16名以内。K16是一个性价比极高的甜点值。5.3 不同分块尺寸的影响粒度与码率的权衡range块尺寸是压缩比和质量之间最直接的旋钮range块尺寸码率bit压缩比约PSNRdB4×4512×512中65536块×约20bit约3.3:135.28×84096块×27bit约19:131.216×161024块×34bit约60:126.84×4块质量很好但码率也高适合对画质敏感、带宽充足的应用16×16块压缩比漂亮可毕竟细节信息全被大块平均化边缘和纹理区域会明显发糊。真正工程上常用的折中是四叉树自适应平滑区域用大块细节区域递归切小块在相同码率下比固定块质量高不少。这个我放在最后一节聊。5.4 解码质量的主观观察哪些区域相对容易糊看解码出来的图能明显察觉分形压缩的高频能力偏弱头发丝、草叶这类高频纹理区域会呈现出一种“真实但模糊”的质感边缘会有轻微的方块感尤其是在s量化比较粗的时候。平滑区域反而表现不错比如天空、墙壁这类地方很难看出压缩痕迹。这里有个有趣的现象分形压缩解码后的图像细节并不是那种常见插值算法造成的“糊一片”而是保持某种自相似结构。你把解码图放大到800×800再看纹理区域依旧有类似的组织性这是分形方法很独特的一面也是它至今仍有人研究的原因。6. 从踩坑里总结的工程经验6.1 别逐块调用dct2用DCT矩阵乘一下能快一个量级我自己在第一版实现里所有块的特征提取都是调用dct2结果是编码时间里特征提取占了将近四成。后来改成预生成8×8的DCT矩阵用矩阵乘法一次处理一堆块速度立刻提了上来。具体做法是用dctmtx(blockSize)拿到正交归一化的DCT基矩阵T然后对每个8×8块X直接算T·X·T等价于dct2但省掉了函数调用的开销。同理降采样也尽量别在循环里反复调用imresize。domain pool构建时一次降采样到位后面就只剩滑窗取窗和变换了。6.2 平坦区域的匹配陷阱把DC单独拎出来做一层预排序有一个我踩过的坑值得专门说纯AC特征的粗筛会把很多亮度不同但纹理都很平的块排在一起比如两片不同灰度的天空纹理特征都接近零特征距离几乎为零。这时候如果不加约束精匹配就完全靠运气效果很差。解决办法是把DC特征并进粗筛的第一个维度或者在特征向量里对DC单独赋权。更稳的做法是三段式第一步只用DC排序把明显亮度差异大的候选过滤掉第二步用AC特征排序第三步精匹配。这样可以保证平滑块不会错配到另一片亮度差异很大的区域也不会让亮度问题在精匹配阶段白白浪费K的机会。6.3 解码不收敛或出现方块的三个原因解码不收敛多数逃不出三个原因一是s绝对值过大比如超过了1.4迭代会把初始图像里的噪声逐渐放大平坦区域出现雪花花纹二是domain块坐标记录错误解码贴图时取错了源位置直接表现为每块独立但和周围对不上三是迭代次数太少只迭代三四次时块边界特别明显容易被误判成算法问题。我建议先用全搜索编码一个8×8的小图比如128×128跑通收敛情况再去追求速度优化调起来会容易很多。6.4 s和o的量化取舍我对比过s用均匀量化和非线性量化两种方案。s的数值分布本身就靠近0附近而且正负都可能均匀量化到5bit容易让靠近0附近的精度不足。用非线性量化比如把[-1,1]区间拆成更密的中央段和更疏的远端可以在同样bit数下把PSNR提高约0.3dB。这只是个小trick但属于典型的“不加码率只加质量”的优化点。6.5 如果还想再快四叉树自适应分块与并行最后给想继续往下做的人指个方向。固定块尺寸的分形压缩只是入门真正实用化的时候四叉树自适应几乎必选。做法是从16×16或32×32的块开始先尝试匹配如果匹配误差大于阈值就切成四个子块递归处理直到块小到4×4为止。这样平滑区域会用大块细节区域用小块在相同码率下PSNR能再提升1到2dB。另外每一个range块的搜索彼此独立完全可以用parfor并行替代for循环。我在8核机器上把主循环改成parfor后编码时间又缩到原来的三分之一左右。DCT预筛让单个range块运算量很小并行化后整体吞吐非常可观。按照这套DCT加速框架改出来的编码器我自己的使用体感是终于不用为了试验一个参数等五分钟了。反复调整k、K、块尺寸和量化位数的过程中关于“分离亮度与纹理”“用频域特征做预筛”这两个思路的直觉会越来越清楚。如果你正在做类似的图像压缩课题我建议先把全搜索的精匹配和DCT粗筛两个版本都写出来跑通再一步优化这样每一处的加速收益都看得见摸得着比照抄一份大代码一头雾水要有用得多。