1. 为什么卡通纹理分解能同时去噪和保住织物细节1.1 一张含噪图不是只有“信号”和“噪声”做图像去噪最让人头疼的往往不是噪声本身而是噪声和真实纹理纠缠在一起。用Matlab做图像处理时我经常收到类似的需求把一张带噪的织物照片恢复干净同时不能把纱线的纹理、纸张的颗粒感也一起抹掉。普通高斯滤波会把纹理当成噪声全变分TV滤波又容易把细节压平最后我转向了卡通纹理图像分解的思路——把图像拆成平滑的卡通分量和振荡的纹理分量再对纹理分量用块状低秩纹理表征去做恢复。这篇就聊聊我基于Matlab的实现过程。很多教科书把加性噪声描述成一个与图像内容无关的额外信号观测图 f 干净图 噪声所以只要做一个低通滤波就能把高频噪声“滤”掉。这种说法在理想频带分离的情况下成立但真实图像根本不会这么听话。显微图像中的丝状结构、布料纹理、不规则表面颗粒它们的高频成分和随机噪声的高频成分重叠在一起。低通滤波之后图像确实安静了但织物纹理、皮肤毛孔、草地质感也一起被磨平了。这不是去噪这是把图像里的“有意义的震荡”和“无意义的震荡”一起删掉。卡通纹理分解的思路不一样它把图像显式地写成卡通分量 u 和纹理分量 v 的叠加即 f ≈ u v。u 负责描述大块平滑区域和明显的边缘轮廓v 负责描述周期、震荡、重复式的细节。这样一来去噪不再是一味地抹高频而是对 u 和 v 分别施加合适的先验让噪声在两者中被重新分配。1.2 常见卡通纹理分解模型的短板经典的全变分去噪模型 ROF即最小化梯度 L1 范数min_u ||∇u||_1 λ/2 ||f - u||_2^2它的优点是能在保留边缘的同时抑制噪声。问题在于TV 模型从设计上就默认“图像内容主要是分段平滑”遇到纹理密集区域时它会把这些振荡细节当成需要消除的梯度跳变于是纹理被“压平”。传统做法里如果只输出 u 作为去噪结果织物类图像基本必糊。后来有研究者把纹理显式建模进目标比如 Meyer 的 G 范数模型用散度空间来描述纹理或者用局部傅里叶基、小波基去拟合震荡成分。这些方法比单 TV 更能保护细节但真正的难点不在“能不能把纹理分离出来”而在“分离出来的纹理 v 本身也是带噪的”。如果直接把 v 拿去重建噪声依然存在。所以在卡通纹理分解框架下真正要解决的是两个问题第一让卡通 u 尽量干净且保留主轮廓第二让纹理 v 尽量真实同时把混合在其中的随机噪声剔除掉。第二个问题正是“块状低秩纹理表征”能发力的地方。1.3 低秩先验与块状策略为什么有效低秩的核心观察是自然纹理具有很强的局部自相似性。一块规律排列的条纹在矩阵形式下任意两行都近似线性相关所以这个纹理块矩阵的秩很低一片布纹、砖墙、水波纹局部看也能用少数几个基模式叠加出来。反过来高斯白噪声是逐像素独立的对应的矩阵几乎没有低秩结构它的奇异值从第一个到最后一个都弥散地分布着。如果把整幅图直接做矩阵低秩分解通常不成立因为图像全局并不低秩。但把图像切成很多局部块每个块内的纹理往往满足低秩假设。这就是“块状”的意义局部地应用低秩先验比全局低秩更符合图像的实际结构。更进一步还可以把相似块聚成组组内所有块拉成矩阵的列再对这个组矩阵做低秩逼近这也就是 WNNM、BM3D 这类算法的底层逻辑。在本篇实现里我以“分块低秩”为基础因为它在 Matlab 里实现直观、调参简单也足够支撑卡通纹理分解中的纹理去噪。2. 块状低秩纹理表征从直觉到数学建模2.1 纹理块矩阵为什么是低秩的为了说明问题可以做一个简单实验任意取一张条纹图像的一块 16×16 区域把像素值直接当成矩阵。由于条纹沿着某个方向反复出现这个矩阵的列向量之间高度相关SVD 之后通常前两三个奇异值就占据了 99% 的能量。随机噪声块则完全不同它的奇异值衰减非常慢几乎均匀分布。在卡通纹理分解框架里纹理分量 v 对应的是图像中振荡、重复的结构。局部块内这些结构天然低秩。所以对纹理块做奇异值软阈值实际上就是在说我相信真实纹理只由少数几个主要模式构成那些零散分布在所有奇异值上的小尾巴大概率是噪声。如果改用“相似块分组”的低秩模型逻辑也类似只是把多个块向量排成矩阵用块与块之间的相关性来强化低秩性。实际效果通常比单块低秩更强但计算量和内存开销也更大。本篇先处理单块版本让流程先跑通。2.2 核范数与奇异值软阈值最自然的低秩逼近工具矩阵的低秩约束在数学上很难直接优化但它的凸松弛形式——核范数 ||X||_* Σ σ_i(X)也就是奇异值之和——计算非常优雅。要解决形如min_X 1/2 ||X - Y||F^2 τ ||X||*的优化问题时最优解是奇异值软阈值算子也叫 SVTD_τ(Y) U · max(S - τ, 0) · V^T其中 Y U S V^T 是 Y 的奇异值分解S 是奇异值矩阵τ 是阈值。直观理解把奇异值谱整体往下压 τ负的直接归零。奇异值大的主要结构保留奇异值小的噪声分量被削掉。这个算子与标量信号的软阈值去噪异曲同工只是作用在矩阵的奇异值谱上。另一个需要考虑的替代方案是硬阈值也就是只保留大于 τ 的奇异值不收缩。硬阈值在“我知道确切的秩”时更好但真实图像纹理的秩往往不确定硬阈值容易产生振铃和不连续。软阈值因为有一个连续收缩的过程在实际的交替迭代里稳定很多。2.3 整体优化目标与交替求解思路把卡通纹理分解和块状低秩纹理表征组合起来可以得到一个联合优化目标min_{u,v} 1/2 ||f - u - v||F^2 α TV(u) β Σ_i ||R_i v||*其中 R_i 表示取第 i 个图像块的操作TV(u) 是卡通分量的全变分正则后一项是纹理分量的块状核范数和f 是含噪观测图u v 是对去噪重建图的估计。数据项保证分解之后两者加起来仍然忠于观测TV 项让 u 分段平滑核范数项让 v 的每个局部块低秩。这个目标函数耦合了 u 和 v直接一起优化很麻烦所以我采用交替方向的思想固定 v更新 u此时问题退化为对残差图 g f - v 做一个标准 TV 去噪。固定 u更新 v此时 r f - u 被看作“待恢复纹理的含噪版本”对 r 的每个局部块做 SVT再重建为完整的纹理分量。重复以上两步直到残差能量稳定。严格来说第 2 步的“直接把 r 做近端算子”是近端梯度法的一种简化形式在块与块重叠的情况下还包含一个平均投影。为了收敛更稳我通常会在更新 v 时加入一个松弛系数 ηv_new (1 - η)·v_old η·SVT_based(r)。这样能有效抑制振荡后面实验部分我会详细说明。3. Matlab主流程交替迭代的骨架设计3.1 主函数结构与核心代码在 Matlab 中我把整个算法封装成一个主函数输入是灰度图 I 和结构体 opts输出是卡通分量 U、纹理分量 V 以及重建去噪图 ImRec。下面是最基本的骨架function [U, V, ImRec] cartoonTextureDecomp(I, opts) % 卡通纹理分解 块状低秩纹理去噪 % I : double型灰度图取值 [0,1] % opts 字段 % lambdaTV : 卡通分量TV正则强度 % tvIter : 内部TV迭代次数 % lambdaLR : 纹理分块低秩阈值 % patchSize: 纹理块尺寸 % step : 分块步长 % iterNum : 外层交替迭代次数 % eta : 松弛系数 (0,1] f I; U f; % 卡通分量初值 V zeros(size(f)); % 纹理分量初值 energy zeros(opts.iterNum, 1); for k 1:opts.iterNum % 1. 固定纹理V更新卡通U对 f - V 做TV去噪 U tvDenoise(f - V, opts.lambdaTV, opts.tvIter); % 2. 固定卡通U更新纹理V对 f - U 做块状低秩近端 R f - U; Vnew blockLowRankProx(R, opts.patchSize, opts.step, opts.lambdaLR); % 松弛更新提升稳定性 V (1 - opts.eta) * V opts.eta * Vnew; % 3. 记录能量与收敛情况 residue f - U - V; energy(k) 0.5 * norm(residue(:), 2)^2; if k 1 abs(energy(k) - energy(k-1)) 1e-5 break; end end ImRec U V; end实际使用时我建议把tvDenoise和blockLowRankProx单独写成函数方便逐个模块调试。主循环里的energy数组很重要它不仅能判断收敛还能帮你发现参数是否设置得不合理如果能量曲线一直剧烈上下跳说明松弛系数太大或者阈值过强。3.2 为什么用交替方向而不是一步到位有人会问既然目标函数已经写出来了为什么不用现成的凸优化求解器一次性求解原因很实际图像尺寸稍微大一点比如 512×512目标函数里的每个图像块是一个变量子空间所有块重叠在一起直接组合成一个大规模优化问题会非常笨重。而 TV 问题和低秩问题各自都有非常成熟的快速算子TV 可以用 Chambolle 对偶投影低秩可以用奇异值软阈值。交替方向的价值就是把一个复杂耦合问题拆成两个能“闭式求解”的子问题每次迭代都只调用成熟算子而不是去硬解一个全局大矩阵。这种拆解在工程上还有一个好处便于监控。第 1 次迭代后你能直接看到卡通分量和纹理分量分别长什么样如果纹理分量出现大块残渣说明 TV 的正则太弱卡通把噪声漏下去了如果卡通分量里还残留条纹说明 TV 太强把纹理过早地压进了卡通。这种“中间过程可视”的调试体验是直接丢给优化器无法获得的。3.3 收敛判断与中间结果监视除了用残差能量判断收敛我强烈建议在迭代过程中用 subplot 实时显示 U、V、residue 三张图。第一次跑通时你大概率会遇到两种现象残差图里出现明显的“污染区”或者纹理分量快速变成全图灰蒙蒙一片。这些都需要通过观察中间结果来定位问题。我的监视代码很简单if mod(k, 10) 1 subplot(1,3,1); imshow(U, []); title(卡通 U); subplot(1,3,2); imshow(V, []); title(纹理 V); subplot(1,3,3); imshow(f - U - V, []); title(残差); drawnow; end如果残差图接近零均值白噪声说明分解已经比较充分如果残差里还有明显的结构说明当前参数下卡通和纹理之间没有分配好。一般迭代 30 到 50 次就足够稳定了不需要追求理论上的完全收敛。4. 两个核心函数的Matlab实现细节4.1 卡通分量用全变分压住噪声但保留边缘TV 求解我使用的是 Chambolle 对偶投影法。它不是最快的算法但代码量小对这个小项目足够用。这个算法的关键点是引入一个对偶变量 p (px, py)迭代更新 p 后卡通分量用 u g - λ·div(p) 恢复。function u tvDenoise(g, lambda, numIter) % Chambolle对偶投影法求解 ROF TV去噪 [rows, cols] size(g); px zeros(rows, cols); py zeros(rows, cols); tau 0.25; % 稳定性步长一般不超过0.25 for k 1:numIter % 计算 div(p) divp circshift(px, [0 1]) - px circshift(py, [1 0]) - py; % 对偶变量目标中的梯度项 term divp - g / lambda; gx circshift(term, [0 -1]) - term; gy circshift(term, [-1 0]) - term; % 梯度上升更新 pxn px tau * gx; pyn py tau * gy; % 投影到单位圆保证 |p| 1 nrm max(1, sqrt(pxn.^2 pyn.^2)); px pxn ./ nrm; py pyn ./ nrm; end divp circshift(px, [0 1]) - px circshift(py, [1 0]) - py; u g - lambda * divp; end这里特别提醒circshift的位移方向决定了梯度算子是前向差分还是后向差分不能随手乱改。term计算如果发现除以 lambda 后数值范围不稳可以把图像先归一化到 [0,1]lambda 取 0.010.1 量级迭代 3060 次就够。TV 函数内部不需要追求完全收敛因为外层交替迭代还会反复调用它每次到“差不多”的状态即可。4.2 纹理分量分块SVT低秩更新纹理分量的核心是 blockLowRankProx 函数。它对输入 R当前残差做重叠分块每个块做一次 SVD 和奇异值软阈值再放回原位置。代码实现时要注意Matlab 里[U, S, Vt] svd(blk, econ)返回的第三个变量是 V^T重建时要写成U * S * Vt这个细节容易写错。function Vout blockLowRankProx(R, patchSize, step, tau) [rows, cols] size(R); acc zeros(rows, cols); % 累加器 cnt zeros(rows, cols); % 计数权重 for i 1:step:rows-patchSize1 for j 1:step:cols-patchSize1 blk R(i:ipatchSize-1, j:jpatchSize-1); [U, S, Vt] svd(blk, econ); S max(S - tau, 0); % 奇异值软阈值 newBlk U * S * Vt; % 重建低秩块 acc(i:ipatchSize-1, j:jpatchSize-1) ... acc(i:ipatchSize-1, j:jpatchSize-1) newBlk; cnt(i:ipatchSize-1, j:jpatchSize-1) ... cnt(i:ipatchSize-1, j:jpatchSize-1) 1; end end Vout acc ./ max(cnt, 1); end这个函数有两点值得展开。第一阈值 tau 对应理论部分的 τ它直接控制纹理去噪强度。tau 太大弱纹理会被当作噪声削掉tau 太小噪声残留在纹理里。第二因为每个像素会被多个重叠块覆盖所以不是简单“覆盖回去”而是累加后取平均否则会出现明显的块状网格伪影。max(cnt,1)是为了避免边界处计数为零时出现 NaN。4.3 重叠块步长与块尺寸的实际选法步长 step 决定了块与块之间的重叠程度。step 1 时每个像素都被大量块覆盖结果最平滑但计算量会暴涨step patchSize 时块完全不重叠速度最快但重建图像容易出现拼缝。我一般取 patchSize/4比如 8 的块用 step216 的块用 step4。这样既保证了块间连续性又不至于让计算量翻太多倍。块尺寸 patchSize 的选取需要参考纹理尺度。条纹很细的布纹8×8 或 10×10 就够用纹理周期比较大的砖墙、波纹建议用 16×16 甚至 20×20。理论上块至少得覆盖一个完整纹理周期低秩结构才明显。如果块太小块内基本就是常数SVT 对噪声的抑制有限如果块太大局部低秩假设被破坏blockLowRankProx 会把不相关的结构强行“压秩”导致纹理失真。5. 参数怎么调一张表看清楚影响方向5.1 先搭一个合成噪声实验调参之前建议先构造一个可复现的合成噪声实验。我常用下面的代码生成测试图orig im2double(imread(barbara.png)); rng(2024); noisy orig 0.05 * randn(size(orig));然后用 PSNR 和 SSIM 来评价输出ImRec U V相对orig的质量。为什么要额外保留一个干净原图因为去噪质量不能靠肉眼“觉得”必须有一个量化基线。真实噪音场景里没有原图但在开发算法阶段量化评价能快速暴露参数问题。5.2 参数影响表我把这个实现里最容易影响结果的参数整理成下面这张表方便大家对照调整参数数值增大数值减小参考范围patchSize更适合大尺度纹理低秩更明显保留更小尺度细节但去噪能力变弱816step计算更快块间连续性变差重叠更多结果更平滑耗时上升24lambdaTV卡通更平滑纹理更容易进入v卡通保留细节但可能留噪0.010.1lambdaLR纹理更干净弱纹理可能被削掉弱纹理保留更好但噪声也会留下0.010.2tvIter每次TV子问题解得更透耗时上升子问题欠收敛卡通不够干净3060iterNum内外迭代更充分超过后会过平滑欠迭代卡通纹理分离不彻底3050eta接近1时收敛快但容易振荡接近0时稳定但收敛慢0.50.9以我自己实验为例对 Barbara 测试图加 σ0.05 的高斯噪声经过简单网格搜索选择 patchSize8、step2、lambdaTV0.02、lambdaLR0.08、eta0.8得到的结果比直接高斯滤波高大约 2dB比纯 TV 去噪高约 0.8dB。最让我满意的是布纹和裤子纹理的视觉保留程度明显好于 TV这是 PSNR 数据之外更容易感知的差异。5.3 我习惯的调参顺序很多新手拿到算法就一头扎进参数网格搜索我觉得效率不高。我自己的顺序是先调 lambdaTV让卡通分量看起来“干净但不失边缘”再调 lambdaLR让纹理分量不再明显含噪同时不把主要条纹洗掉接着调 patchSize观察纹理细节的还原情况最后调 step 和 eta 来平衡平滑度和计算速度。有个反直觉的经验是lambdaLR 并不是越大越好。我第一次调参时把 lambdaLR 调到 0.2纹理分量确实非常平滑但原本清晰的布纹也被磨成了“扁平的不规则块”。原因是弱纹理的奇异值本身就不大阈值过大时会和噪声一起被削掉。正确的做法是保持一个适度的阈值把噪声的“长尾奇异值”削掉而不是把整个谱都压下去。6. 我在Matlab实现中踩过的几个坑6.1 SVD之前忘记归一化导致阈值完全失效第一次跑通时我用imread直接读图没有转成 double也没有归一化然后沿用 [0,1] 尺度下调好的 tau0.05。结果纹理分量基本没有变化去噪效果约等于零。排查了很久才发现uint8 图像的像素范围是 0255纹理块的奇异值量级在 1000 左右阈值 0.05 对它来说小到可以忽略。修复很简单进入算法前统一执行I im2double(I)所有后续参数都在 [0,1] 量级下校准。这类问题隐蔽在结果上不仔细看奇异值谱很难发现。建议在写函数时第一行就检查输入类型if ~isa(I, double), I im2double(I); end6.2 边界填充不当边缘出现一圈暗纹分块提取时如果不做边界处理循环到i rows-patchSize1就停意味着图像最后几行和最后几列像素永远参与不到纹理分量的低秩更新。结果就是后来我看到的去噪重建图右边和底边明显比内部暗像加了边框。其实这是处理边界时的索引截断问题。修复方法是先对 R 做padarray用replicate模式扩展边缘再分块处理处理完裁剪回原始尺寸。这样边缘像素也有对应的块参与重建暗纹基本消失。这里不建议用symmetric或者补零因为复制边缘在实际图像上更自然。padSize floor(patchSize/2); Rp padarray(R, [padSize, padSize], replicate); % 分块循环处理Rp最后裁剪 Vout Vout(padSize1:rowspadSize, ...)6.3 重叠块不累加权重直接覆盖导致“马赛克”这是我踩过最直观的坑。起初我没写acc和cnt而是直接把newBlk写回V的对应位置。结果相邻块各自为政重建出来的纹理分量充满了拼接痕迹像低分辨率马赛克。原因是重叠区域的每个像素被多次赋值最后一次覆盖把前面所有块的贡献都丢掉了。正确的做法就是用累加器求和最后除以每个像素被覆盖的次数。这本质上是对同一像素的多个块估计做一个平均既能消除拼接缝又能降低单块估计的方差。代码就两行为acc和cnt但缺了它整个结果质量会下降一个档次。6.4 外循环能量曲线抖动越迭代越差有几次调参时我发现能量曲线在前 10 次迭代快速下降之后突然开始上下抖动PSNR 不升反降。检查后发现eta我设成了 1相当于每次完全用最新的低秩纹理替换旧纹理。在强噪声场景下这种方式很容易让卡通和纹理之间产生“抢能量”振荡卡通分量把纹理抢走一部分下一次纹理低秩更新又把纹理抢回来来回拉扯。解决方法是把eta降到 0.7 左右每次只更新一部分。这相当于给交替迭代装了个阻尼器虽然收敛稍慢但稳定很多。如果还抖就继续降到 0.5。这个技巧在我处理很多类似分裂算法时都非常管用。6.5 分块循环太慢怎么提速Matlab 的双重 for 循环在小图上看不出问题一旦图变成 1024×1024blockLowRankProx 可能是整段程序最耗时的部分。我先用profile定位到它然后做了三件事预处理分配acc和cnt避免循环内动态扩容把svd的输入输出从[U,S,V]改成[U,S,Vt]避免多余的共轭转置如果电脑有多核把内层行列循环改成parfor。但如果使用parfor要注意循环里acc和cnt的累加不是并行的安全操作。更稳妥的优化方案是用im2col一次性取得所有块把循环转成矩阵运算再统一做批量 SVD。不过批量 SVD 在 Matlab 里仍需循环实际提升有限。对我个人项目来说parfor结合合适步长已经够用毕竟核心目的不是写一个生产级高速库而是把算法逻辑验证清楚。最后分享一点个人体会这类方法的调参就像在“卡通”和“纹理”之间找一个平衡点噪声水平、图像内容、纹理尺度都会影响最终参数。我在实际项目中逐渐放弃了对所有图像用一组固定参数的想法而是针对每一类图像做一次小范围参数搜索往往比追求一个“万能参数组合”省时间得多。如果能搭配一个简单的 UI让使用者实时拖动 lambdaTV 和 lambdaLR 两个滑动条观察效果这个 Matlab 实现的可玩性会更高也更接近一个真正的图像处理工具。 Matlab实现卡通纹理分解与块状低秩图像去噪
1. 为什么卡通纹理分解能同时去噪和保住织物细节1.1 一张含噪图不是只有“信号”和“噪声”做图像去噪最让人头疼的往往不是噪声本身而是噪声和真实纹理纠缠在一起。用Matlab做图像处理时我经常收到类似的需求把一张带噪的织物照片恢复干净同时不能把纱线的纹理、纸张的颗粒感也一起抹掉。普通高斯滤波会把纹理当成噪声全变分TV滤波又容易把细节压平最后我转向了卡通纹理图像分解的思路——把图像拆成平滑的卡通分量和振荡的纹理分量再对纹理分量用块状低秩纹理表征去做恢复。这篇就聊聊我基于Matlab的实现过程。很多教科书把加性噪声描述成一个与图像内容无关的额外信号观测图 f 干净图 噪声所以只要做一个低通滤波就能把高频噪声“滤”掉。这种说法在理想频带分离的情况下成立但真实图像根本不会这么听话。显微图像中的丝状结构、布料纹理、不规则表面颗粒它们的高频成分和随机噪声的高频成分重叠在一起。低通滤波之后图像确实安静了但织物纹理、皮肤毛孔、草地质感也一起被磨平了。这不是去噪这是把图像里的“有意义的震荡”和“无意义的震荡”一起删掉。卡通纹理分解的思路不一样它把图像显式地写成卡通分量 u 和纹理分量 v 的叠加即 f ≈ u v。u 负责描述大块平滑区域和明显的边缘轮廓v 负责描述周期、震荡、重复式的细节。这样一来去噪不再是一味地抹高频而是对 u 和 v 分别施加合适的先验让噪声在两者中被重新分配。1.2 常见卡通纹理分解模型的短板经典的全变分去噪模型 ROF即最小化梯度 L1 范数min_u ||∇u||_1 λ/2 ||f - u||_2^2它的优点是能在保留边缘的同时抑制噪声。问题在于TV 模型从设计上就默认“图像内容主要是分段平滑”遇到纹理密集区域时它会把这些振荡细节当成需要消除的梯度跳变于是纹理被“压平”。传统做法里如果只输出 u 作为去噪结果织物类图像基本必糊。后来有研究者把纹理显式建模进目标比如 Meyer 的 G 范数模型用散度空间来描述纹理或者用局部傅里叶基、小波基去拟合震荡成分。这些方法比单 TV 更能保护细节但真正的难点不在“能不能把纹理分离出来”而在“分离出来的纹理 v 本身也是带噪的”。如果直接把 v 拿去重建噪声依然存在。所以在卡通纹理分解框架下真正要解决的是两个问题第一让卡通 u 尽量干净且保留主轮廓第二让纹理 v 尽量真实同时把混合在其中的随机噪声剔除掉。第二个问题正是“块状低秩纹理表征”能发力的地方。1.3 低秩先验与块状策略为什么有效低秩的核心观察是自然纹理具有很强的局部自相似性。一块规律排列的条纹在矩阵形式下任意两行都近似线性相关所以这个纹理块矩阵的秩很低一片布纹、砖墙、水波纹局部看也能用少数几个基模式叠加出来。反过来高斯白噪声是逐像素独立的对应的矩阵几乎没有低秩结构它的奇异值从第一个到最后一个都弥散地分布着。如果把整幅图直接做矩阵低秩分解通常不成立因为图像全局并不低秩。但把图像切成很多局部块每个块内的纹理往往满足低秩假设。这就是“块状”的意义局部地应用低秩先验比全局低秩更符合图像的实际结构。更进一步还可以把相似块聚成组组内所有块拉成矩阵的列再对这个组矩阵做低秩逼近这也就是 WNNM、BM3D 这类算法的底层逻辑。在本篇实现里我以“分块低秩”为基础因为它在 Matlab 里实现直观、调参简单也足够支撑卡通纹理分解中的纹理去噪。2. 块状低秩纹理表征从直觉到数学建模2.1 纹理块矩阵为什么是低秩的为了说明问题可以做一个简单实验任意取一张条纹图像的一块 16×16 区域把像素值直接当成矩阵。由于条纹沿着某个方向反复出现这个矩阵的列向量之间高度相关SVD 之后通常前两三个奇异值就占据了 99% 的能量。随机噪声块则完全不同它的奇异值衰减非常慢几乎均匀分布。在卡通纹理分解框架里纹理分量 v 对应的是图像中振荡、重复的结构。局部块内这些结构天然低秩。所以对纹理块做奇异值软阈值实际上就是在说我相信真实纹理只由少数几个主要模式构成那些零散分布在所有奇异值上的小尾巴大概率是噪声。如果改用“相似块分组”的低秩模型逻辑也类似只是把多个块向量排成矩阵用块与块之间的相关性来强化低秩性。实际效果通常比单块低秩更强但计算量和内存开销也更大。本篇先处理单块版本让流程先跑通。2.2 核范数与奇异值软阈值最自然的低秩逼近工具矩阵的低秩约束在数学上很难直接优化但它的凸松弛形式——核范数 ||X||_* Σ σ_i(X)也就是奇异值之和——计算非常优雅。要解决形如min_X 1/2 ||X - Y||F^2 τ ||X||*的优化问题时最优解是奇异值软阈值算子也叫 SVTD_τ(Y) U · max(S - τ, 0) · V^T其中 Y U S V^T 是 Y 的奇异值分解S 是奇异值矩阵τ 是阈值。直观理解把奇异值谱整体往下压 τ负的直接归零。奇异值大的主要结构保留奇异值小的噪声分量被削掉。这个算子与标量信号的软阈值去噪异曲同工只是作用在矩阵的奇异值谱上。另一个需要考虑的替代方案是硬阈值也就是只保留大于 τ 的奇异值不收缩。硬阈值在“我知道确切的秩”时更好但真实图像纹理的秩往往不确定硬阈值容易产生振铃和不连续。软阈值因为有一个连续收缩的过程在实际的交替迭代里稳定很多。2.3 整体优化目标与交替求解思路把卡通纹理分解和块状低秩纹理表征组合起来可以得到一个联合优化目标min_{u,v} 1/2 ||f - u - v||F^2 α TV(u) β Σ_i ||R_i v||*其中 R_i 表示取第 i 个图像块的操作TV(u) 是卡通分量的全变分正则后一项是纹理分量的块状核范数和f 是含噪观测图u v 是对去噪重建图的估计。数据项保证分解之后两者加起来仍然忠于观测TV 项让 u 分段平滑核范数项让 v 的每个局部块低秩。这个目标函数耦合了 u 和 v直接一起优化很麻烦所以我采用交替方向的思想固定 v更新 u此时问题退化为对残差图 g f - v 做一个标准 TV 去噪。固定 u更新 v此时 r f - u 被看作“待恢复纹理的含噪版本”对 r 的每个局部块做 SVT再重建为完整的纹理分量。重复以上两步直到残差能量稳定。严格来说第 2 步的“直接把 r 做近端算子”是近端梯度法的一种简化形式在块与块重叠的情况下还包含一个平均投影。为了收敛更稳我通常会在更新 v 时加入一个松弛系数 ηv_new (1 - η)·v_old η·SVT_based(r)。这样能有效抑制振荡后面实验部分我会详细说明。3. Matlab主流程交替迭代的骨架设计3.1 主函数结构与核心代码在 Matlab 中我把整个算法封装成一个主函数输入是灰度图 I 和结构体 opts输出是卡通分量 U、纹理分量 V 以及重建去噪图 ImRec。下面是最基本的骨架function [U, V, ImRec] cartoonTextureDecomp(I, opts) % 卡通纹理分解 块状低秩纹理去噪 % I : double型灰度图取值 [0,1] % opts 字段 % lambdaTV : 卡通分量TV正则强度 % tvIter : 内部TV迭代次数 % lambdaLR : 纹理分块低秩阈值 % patchSize: 纹理块尺寸 % step : 分块步长 % iterNum : 外层交替迭代次数 % eta : 松弛系数 (0,1] f I; U f; % 卡通分量初值 V zeros(size(f)); % 纹理分量初值 energy zeros(opts.iterNum, 1); for k 1:opts.iterNum % 1. 固定纹理V更新卡通U对 f - V 做TV去噪 U tvDenoise(f - V, opts.lambdaTV, opts.tvIter); % 2. 固定卡通U更新纹理V对 f - U 做块状低秩近端 R f - U; Vnew blockLowRankProx(R, opts.patchSize, opts.step, opts.lambdaLR); % 松弛更新提升稳定性 V (1 - opts.eta) * V opts.eta * Vnew; % 3. 记录能量与收敛情况 residue f - U - V; energy(k) 0.5 * norm(residue(:), 2)^2; if k 1 abs(energy(k) - energy(k-1)) 1e-5 break; end end ImRec U V; end实际使用时我建议把tvDenoise和blockLowRankProx单独写成函数方便逐个模块调试。主循环里的energy数组很重要它不仅能判断收敛还能帮你发现参数是否设置得不合理如果能量曲线一直剧烈上下跳说明松弛系数太大或者阈值过强。3.2 为什么用交替方向而不是一步到位有人会问既然目标函数已经写出来了为什么不用现成的凸优化求解器一次性求解原因很实际图像尺寸稍微大一点比如 512×512目标函数里的每个图像块是一个变量子空间所有块重叠在一起直接组合成一个大规模优化问题会非常笨重。而 TV 问题和低秩问题各自都有非常成熟的快速算子TV 可以用 Chambolle 对偶投影低秩可以用奇异值软阈值。交替方向的价值就是把一个复杂耦合问题拆成两个能“闭式求解”的子问题每次迭代都只调用成熟算子而不是去硬解一个全局大矩阵。这种拆解在工程上还有一个好处便于监控。第 1 次迭代后你能直接看到卡通分量和纹理分量分别长什么样如果纹理分量出现大块残渣说明 TV 的正则太弱卡通把噪声漏下去了如果卡通分量里还残留条纹说明 TV 太强把纹理过早地压进了卡通。这种“中间过程可视”的调试体验是直接丢给优化器无法获得的。3.3 收敛判断与中间结果监视除了用残差能量判断收敛我强烈建议在迭代过程中用 subplot 实时显示 U、V、residue 三张图。第一次跑通时你大概率会遇到两种现象残差图里出现明显的“污染区”或者纹理分量快速变成全图灰蒙蒙一片。这些都需要通过观察中间结果来定位问题。我的监视代码很简单if mod(k, 10) 1 subplot(1,3,1); imshow(U, []); title(卡通 U); subplot(1,3,2); imshow(V, []); title(纹理 V); subplot(1,3,3); imshow(f - U - V, []); title(残差); drawnow; end如果残差图接近零均值白噪声说明分解已经比较充分如果残差里还有明显的结构说明当前参数下卡通和纹理之间没有分配好。一般迭代 30 到 50 次就足够稳定了不需要追求理论上的完全收敛。4. 两个核心函数的Matlab实现细节4.1 卡通分量用全变分压住噪声但保留边缘TV 求解我使用的是 Chambolle 对偶投影法。它不是最快的算法但代码量小对这个小项目足够用。这个算法的关键点是引入一个对偶变量 p (px, py)迭代更新 p 后卡通分量用 u g - λ·div(p) 恢复。function u tvDenoise(g, lambda, numIter) % Chambolle对偶投影法求解 ROF TV去噪 [rows, cols] size(g); px zeros(rows, cols); py zeros(rows, cols); tau 0.25; % 稳定性步长一般不超过0.25 for k 1:numIter % 计算 div(p) divp circshift(px, [0 1]) - px circshift(py, [1 0]) - py; % 对偶变量目标中的梯度项 term divp - g / lambda; gx circshift(term, [0 -1]) - term; gy circshift(term, [-1 0]) - term; % 梯度上升更新 pxn px tau * gx; pyn py tau * gy; % 投影到单位圆保证 |p| 1 nrm max(1, sqrt(pxn.^2 pyn.^2)); px pxn ./ nrm; py pyn ./ nrm; end divp circshift(px, [0 1]) - px circshift(py, [1 0]) - py; u g - lambda * divp; end这里特别提醒circshift的位移方向决定了梯度算子是前向差分还是后向差分不能随手乱改。term计算如果发现除以 lambda 后数值范围不稳可以把图像先归一化到 [0,1]lambda 取 0.010.1 量级迭代 3060 次就够。TV 函数内部不需要追求完全收敛因为外层交替迭代还会反复调用它每次到“差不多”的状态即可。4.2 纹理分量分块SVT低秩更新纹理分量的核心是 blockLowRankProx 函数。它对输入 R当前残差做重叠分块每个块做一次 SVD 和奇异值软阈值再放回原位置。代码实现时要注意Matlab 里[U, S, Vt] svd(blk, econ)返回的第三个变量是 V^T重建时要写成U * S * Vt这个细节容易写错。function Vout blockLowRankProx(R, patchSize, step, tau) [rows, cols] size(R); acc zeros(rows, cols); % 累加器 cnt zeros(rows, cols); % 计数权重 for i 1:step:rows-patchSize1 for j 1:step:cols-patchSize1 blk R(i:ipatchSize-1, j:jpatchSize-1); [U, S, Vt] svd(blk, econ); S max(S - tau, 0); % 奇异值软阈值 newBlk U * S * Vt; % 重建低秩块 acc(i:ipatchSize-1, j:jpatchSize-1) ... acc(i:ipatchSize-1, j:jpatchSize-1) newBlk; cnt(i:ipatchSize-1, j:jpatchSize-1) ... cnt(i:ipatchSize-1, j:jpatchSize-1) 1; end end Vout acc ./ max(cnt, 1); end这个函数有两点值得展开。第一阈值 tau 对应理论部分的 τ它直接控制纹理去噪强度。tau 太大弱纹理会被当作噪声削掉tau 太小噪声残留在纹理里。第二因为每个像素会被多个重叠块覆盖所以不是简单“覆盖回去”而是累加后取平均否则会出现明显的块状网格伪影。max(cnt,1)是为了避免边界处计数为零时出现 NaN。4.3 重叠块步长与块尺寸的实际选法步长 step 决定了块与块之间的重叠程度。step 1 时每个像素都被大量块覆盖结果最平滑但计算量会暴涨step patchSize 时块完全不重叠速度最快但重建图像容易出现拼缝。我一般取 patchSize/4比如 8 的块用 step216 的块用 step4。这样既保证了块间连续性又不至于让计算量翻太多倍。块尺寸 patchSize 的选取需要参考纹理尺度。条纹很细的布纹8×8 或 10×10 就够用纹理周期比较大的砖墙、波纹建议用 16×16 甚至 20×20。理论上块至少得覆盖一个完整纹理周期低秩结构才明显。如果块太小块内基本就是常数SVT 对噪声的抑制有限如果块太大局部低秩假设被破坏blockLowRankProx 会把不相关的结构强行“压秩”导致纹理失真。5. 参数怎么调一张表看清楚影响方向5.1 先搭一个合成噪声实验调参之前建议先构造一个可复现的合成噪声实验。我常用下面的代码生成测试图orig im2double(imread(barbara.png)); rng(2024); noisy orig 0.05 * randn(size(orig));然后用 PSNR 和 SSIM 来评价输出ImRec U V相对orig的质量。为什么要额外保留一个干净原图因为去噪质量不能靠肉眼“觉得”必须有一个量化基线。真实噪音场景里没有原图但在开发算法阶段量化评价能快速暴露参数问题。5.2 参数影响表我把这个实现里最容易影响结果的参数整理成下面这张表方便大家对照调整参数数值增大数值减小参考范围patchSize更适合大尺度纹理低秩更明显保留更小尺度细节但去噪能力变弱816step计算更快块间连续性变差重叠更多结果更平滑耗时上升24lambdaTV卡通更平滑纹理更容易进入v卡通保留细节但可能留噪0.010.1lambdaLR纹理更干净弱纹理可能被削掉弱纹理保留更好但噪声也会留下0.010.2tvIter每次TV子问题解得更透耗时上升子问题欠收敛卡通不够干净3060iterNum内外迭代更充分超过后会过平滑欠迭代卡通纹理分离不彻底3050eta接近1时收敛快但容易振荡接近0时稳定但收敛慢0.50.9以我自己实验为例对 Barbara 测试图加 σ0.05 的高斯噪声经过简单网格搜索选择 patchSize8、step2、lambdaTV0.02、lambdaLR0.08、eta0.8得到的结果比直接高斯滤波高大约 2dB比纯 TV 去噪高约 0.8dB。最让我满意的是布纹和裤子纹理的视觉保留程度明显好于 TV这是 PSNR 数据之外更容易感知的差异。5.3 我习惯的调参顺序很多新手拿到算法就一头扎进参数网格搜索我觉得效率不高。我自己的顺序是先调 lambdaTV让卡通分量看起来“干净但不失边缘”再调 lambdaLR让纹理分量不再明显含噪同时不把主要条纹洗掉接着调 patchSize观察纹理细节的还原情况最后调 step 和 eta 来平衡平滑度和计算速度。有个反直觉的经验是lambdaLR 并不是越大越好。我第一次调参时把 lambdaLR 调到 0.2纹理分量确实非常平滑但原本清晰的布纹也被磨成了“扁平的不规则块”。原因是弱纹理的奇异值本身就不大阈值过大时会和噪声一起被削掉。正确的做法是保持一个适度的阈值把噪声的“长尾奇异值”削掉而不是把整个谱都压下去。6. 我在Matlab实现中踩过的几个坑6.1 SVD之前忘记归一化导致阈值完全失效第一次跑通时我用imread直接读图没有转成 double也没有归一化然后沿用 [0,1] 尺度下调好的 tau0.05。结果纹理分量基本没有变化去噪效果约等于零。排查了很久才发现uint8 图像的像素范围是 0255纹理块的奇异值量级在 1000 左右阈值 0.05 对它来说小到可以忽略。修复很简单进入算法前统一执行I im2double(I)所有后续参数都在 [0,1] 量级下校准。这类问题隐蔽在结果上不仔细看奇异值谱很难发现。建议在写函数时第一行就检查输入类型if ~isa(I, double), I im2double(I); end6.2 边界填充不当边缘出现一圈暗纹分块提取时如果不做边界处理循环到i rows-patchSize1就停意味着图像最后几行和最后几列像素永远参与不到纹理分量的低秩更新。结果就是后来我看到的去噪重建图右边和底边明显比内部暗像加了边框。其实这是处理边界时的索引截断问题。修复方法是先对 R 做padarray用replicate模式扩展边缘再分块处理处理完裁剪回原始尺寸。这样边缘像素也有对应的块参与重建暗纹基本消失。这里不建议用symmetric或者补零因为复制边缘在实际图像上更自然。padSize floor(patchSize/2); Rp padarray(R, [padSize, padSize], replicate); % 分块循环处理Rp最后裁剪 Vout Vout(padSize1:rowspadSize, ...)6.3 重叠块不累加权重直接覆盖导致“马赛克”这是我踩过最直观的坑。起初我没写acc和cnt而是直接把newBlk写回V的对应位置。结果相邻块各自为政重建出来的纹理分量充满了拼接痕迹像低分辨率马赛克。原因是重叠区域的每个像素被多次赋值最后一次覆盖把前面所有块的贡献都丢掉了。正确的做法就是用累加器求和最后除以每个像素被覆盖的次数。这本质上是对同一像素的多个块估计做一个平均既能消除拼接缝又能降低单块估计的方差。代码就两行为acc和cnt但缺了它整个结果质量会下降一个档次。6.4 外循环能量曲线抖动越迭代越差有几次调参时我发现能量曲线在前 10 次迭代快速下降之后突然开始上下抖动PSNR 不升反降。检查后发现eta我设成了 1相当于每次完全用最新的低秩纹理替换旧纹理。在强噪声场景下这种方式很容易让卡通和纹理之间产生“抢能量”振荡卡通分量把纹理抢走一部分下一次纹理低秩更新又把纹理抢回来来回拉扯。解决方法是把eta降到 0.7 左右每次只更新一部分。这相当于给交替迭代装了个阻尼器虽然收敛稍慢但稳定很多。如果还抖就继续降到 0.5。这个技巧在我处理很多类似分裂算法时都非常管用。6.5 分块循环太慢怎么提速Matlab 的双重 for 循环在小图上看不出问题一旦图变成 1024×1024blockLowRankProx 可能是整段程序最耗时的部分。我先用profile定位到它然后做了三件事预处理分配acc和cnt避免循环内动态扩容把svd的输入输出从[U,S,V]改成[U,S,Vt]避免多余的共轭转置如果电脑有多核把内层行列循环改成parfor。但如果使用parfor要注意循环里acc和cnt的累加不是并行的安全操作。更稳妥的优化方案是用im2col一次性取得所有块把循环转成矩阵运算再统一做批量 SVD。不过批量 SVD 在 Matlab 里仍需循环实际提升有限。对我个人项目来说parfor结合合适步长已经够用毕竟核心目的不是写一个生产级高速库而是把算法逻辑验证清楚。最后分享一点个人体会这类方法的调参就像在“卡通”和“纹理”之间找一个平衡点噪声水平、图像内容、纹理尺度都会影响最终参数。我在实际项目中逐渐放弃了对所有图像用一组固定参数的想法而是针对每一类图像做一次小范围参数搜索往往比追求一个“万能参数组合”省时间得多。如果能搭配一个简单的 UI让使用者实时拖动 lambdaTV 和 lambdaLR 两个滑动条观察效果这个 Matlab 实现的可玩性会更高也更接近一个真正的图像处理工具。