我要提问
ARTICLE DETAIL

资讯详情

前沿编程新知与开发实战干货的深度解读。

二维图像压缩感知重构MATLAB实战:OMP、CoSaMP与IRLS算法实现与调优

二维图像压缩感知重构MATLAB实战:OMP、CoSaMP与IRLS算法实现与调优 简介这份资源面向数字图像处理、计算机视觉与医学成像方向的学习者和研究者提供二维图像压缩感知重构算法的完整MATLAB实现帮助读者从代码层面理解CS理论如何用远少于像素数的测量值重建图像。压缩包共26个文件以18个m脚本为核心覆盖测量矩阵生成、图像采样、ISTA/FISTA等重构算法与性能评估流程另含7个fig结果图和1个bmp测试图像整体约1.59MB结构紧凑便于逐文件调试。已有4295人学习下载说明其在教学与实验场景中具备一定参考价值。读者可据此对比高斯、伯努利等不同测量矩阵的重构质量调整松弛因子与阈值参数并借助MSE、PSNR指标和可视化结果验证算法效果也可作为向视频压缩、医学影像、遥感分析及深度学习结合方向拓展的实验起点。1. 二维图像压缩感知重构一份能直接跑通的 MATLAB 代码包做图像处理的人迟早会撞上压缩感知这个坎。采样率明明压到了 30% 甚至更低重构出来的图却要跟原图肉眼难辨这中间的落差就是重构算法要填的坑。这份资源给的不是论文里的伪代码而是一套完整的 MATLAB 程序代码专门针对二维图像做压缩感知重构。它把测量矩阵生成、稀疏基构造、重构求解、PSNR 对比这几步全串起来了拿到手改个图片路径就能出结果。适合正在做图像处理大作业、想验证 OMP 或 CoSaMP 这类算法实际效果的人也适合需要快速搭一个重构 baseline 再往上改的研究场景。MATLAB 版本不用追新R2016b 之后的都能跑核心函数没用到太偏的工具箱。2. 压缩感知重构的数学骨架从测量矩阵到稀疏基的选型逻辑2.1 为什么二维图像不能直接套一维压缩感知一维压缩感知的模型是 y Φxx 是长度 N 的列向量Φ 是 M×N 的测量矩阵M 远小于 N。但二维图像是 N×N 的矩阵如果直接按列拉伸成 N²×1 的向量测量矩阵的尺寸会膨胀到 M×N²内存和计算量都吃不消。常见做法是分块处理把图像切成 32×32 或 64×64 的小块每块单独做压缩感知最后拼回去。另一种做法是用可分离的测量矩阵行方向和列方向各乘一个矩阵把二维运算拆成两次一维运算。这份代码走的是分块路线块大小在参数区可以改默认给的是 32×32兼顾了重构质量和速度。分块之后有个细节容易被忽略块与块之间的边界效应。如果每块独立重构再硬拼块边缘会出现明显的方格感。代码里用了一个简单的重叠分块策略块间重叠 4 个像素重构时对重叠区域做加权平均。这个处理不复杂但有没有它视觉效果差很多。2.2 测量矩阵的选取高斯矩阵还是伯努利矩阵测量矩阵要满足 RIP 条件实际中没法直接验证所以选的都是大概率满足 RIP 的随机矩阵。高斯矩阵是最常用的每个元素服从 N(0, 1/M)跟绝大多数稀疏基都不相关。伯努利矩阵元素取 ±1/√M硬件实现时更友好但重构质量在高压缩比下会略逊于高斯矩阵。代码里默认用的是高斯矩阵生成方式如下% 生成 M x N 的高斯测量矩阵M 为测量数N 为信号长度 function Phi generate_gaussian_matrix(M, N) % 每个元素服从 N(0, 1/M)保证列向量近似正交 Phi randn(M, N) / sqrt(M); % 对每一列做归一化避免数值量级差异过大 Phi Phi ./ repmat(sqrt(sum(Phi.^2, 1)), M, 1); end这里做了列归一化原因是后续重构算法对矩阵列向量的能量敏感。如果不归一化某些列的范数偏大会导致求解时权重失衡重构出来的图像局部过亮或过暗。归一化之后每列能量为 1求解稳定性明显提升。参数 M 就是测量数M/N 就是采样率代码里通过采样率反算 M比如采样率 0.3、块大小 32×32 时N1024M307。2.3 稀疏基的选择DCT 还是小波自然图像在 DCT 基下是近似稀疏的在小波基下更稀疏。DCT 基构造简单就是余弦函数采样计算快小波基需要滤波器组多一层分解但压缩感知重构质量通常更好。代码里两种基都给了通过一个 flag 切换。DCT 基的构造代码% 构造 N x N 的 DCT 正交基 function Psi generate_dct_basis(N) Psi zeros(N, N); for k 0:N-1 for n 0:N-1 if k 0 Psi(k1, n1) 1 / sqrt(N); else Psi(k1, n1) sqrt(2/N) * cos(pi * k * (2*n1) / (2*N)); end end end end这个双重循环在 N1024 时会比较慢实际跑的时候可以用 dctmtx 函数替代MATLAB 自带一行就能生成。代码里保留手写版本是为了让读者看清 DCT 基的每个元素怎么来的实际使用时建议换成 dctmtx。小波基用的是 sym4 小波通过 wavedec2 做两层分解再把系数拉成列向量。注意小波分解后的系数长度跟原信号长度一致但排列顺序会影响稀疏性代码里做了按能量排序的处理把大系数排在前面这样重构时阈值收缩更有效。3. 重构算法实战OMP、CoSaMP 与 IRLS 的代码实现与参数调优3.1 OMP 算法的 MATLAB 实现与停止条件设置OMP 是正交匹配追踪思路是每次迭代选一个跟残差最相关的原子然后对已选原子做最小二乘更新。它的优点是实现简单、速度快缺点是需要知道稀疏度 KK 估不准就重构不好。代码里的 OMP 核心循环function x_rec omp_reconstruct(y, Phi, K) % y: 测量向量 M x 1 % Phi: 测量矩阵 M x N % K: 稀疏度 [M, N] size(Phi); residual y; support_set []; x_rec zeros(N, 1); for iter 1:K % 计算残差与各原子的内积找最大相关原子 correlations abs(Phi * residual); [~, idx] max(correlations); support_set [support_set, idx]; % 在已选原子集上做最小二乘 Phi_sub Phi(:, support_set); x_sub pinv(Phi_sub) * y; % 更新残差 residual y - Phi_sub * x_sub; % 残差能量低于阈值时提前退出 if norm(residual) 1e-6 break; end end x_rec(support_set) x_sub; end停止条件有两个迭代次数达到 K或者残差能量低于 1e-6。实际跑的时候如果 K 设得偏大OMP 会把噪声也选进来重构图像出现颗粒感。K 设得偏小细节丢失。经验值是 K 取 N 的 10% 到 20%具体看图像纹理复杂度。代码里 K 是通过稀疏度估计函数自动算的用 DCT 系数排序后取累计能量 95% 的位置作为 K比固定值靠谱。3.2 CoSaMP 算法的回溯机制与参数影响CoSaMP 跟 OMP 的区别在于每次迭代选 2K 个原子然后回溯剪枝只保留 K 个。这个回溯机制让它在稀疏度估计不准时更鲁棒代价是每次迭代计算量翻倍。代码里 CoSaMP 的关键步骤function x_rec cosamp_reconstruct(y, Phi, K) [M, N] size(Phi); residual y; support_set []; x_rec zeros(N, 1); max_iter 50; for iter 1:max_iter % 步骤1形成代理信号找 2K 个最大相关原子 proxy Phi * residual; [~, sort_idx] sort(abs(proxy), descend); new_support sort_idx(1:min(2*K, N)); % 步骤2合并支撑集 support_set union(support_set, new_support); % 步骤3在合并支撑集上做最小二乘 Phi_sub Phi(:, support_set); x_sub pinv(Phi_sub) * y; % 步骤4剪枝只保留 K 个最大系数 [~, keep_idx] sort(abs(x_sub), descend); keep_idx keep_idx(1:min(K, length(keep_idx))); support_set support_set(keep_idx); x_sub x_sub(keep_idx); % 步骤5更新残差 residual y - Phi(:, support_set) * x_sub; if norm(residual) 1e-6 break; end end x_rec(support_set) x_sub; endmax_iter 默认 50实际跑的时候大多数图像在 20 次以内就收敛了。如果跑到 50 次还没收敛通常是测量矩阵条件数太差或者稀疏基选得不合适。回溯剪枝那一步是 CoSaMP 的核心它保证了支撑集大小始终不超过 K避免 OMP 那种选错原子就一路错到底的问题。代价是每次迭代要做两次排序N 大时排序开销不可忽略。3.3 IRLS 算法用迭代重加权最小二乘逼近 L0 范数IRLS 不走贪婪路线它用一系列加权 L2 问题去逼近 L0 范数。每次迭代根据当前解的幅度调整权重幅度小的系数权重加大逼着它们往零走。function x_rec irls_reconstruct(y, Phi, max_iter, epsilon) [M, N] size(Phi); x_rec Phi * y; % 初始解用匹配滤波 weights ones(N, 1); for iter 1:max_iter % 构造加权矩阵 W diag(weights); % 求解加权最小二乘min ||Wx||_2 s.t. Phi*x y % 用拉格朗日乘子法转为无约束问题 x_rec W \ (Phi * ((Phi / W) * Phi \ y)); % 更新权重epsilon 防止除零 weights 1 ./ (abs(x_rec) epsilon); % 权重变化很小时退出 if norm(weights - 1./(abs(x_rec)epsilon)) 1e-4 break; end end endepsilon 是个关键参数默认给 1e-3。太大则稀疏性不够太小则数值不稳定。IRLS 的优点是重构精度通常比 OMP 高缺点是要解线性方程组每次迭代计算量比 OMP 大一个量级。图像块小的时候32×32还能接受块大到 128×128 就明显慢了。代码里给 IRLS 设了最大迭代 30 次实际跑下来 15 次左右权重就基本稳定了。4. 避坑与排查重构图像质量不达标时先查这五处4.1 重构图像全黑或全白现象跑完代码imshow 出来的图要么全黑要么全白PSNR 低到个位数。原因测量矩阵没有归一化或者稀疏基的逆变换用错了。全黑通常是重构系数全为零全白是系数爆炸。检查 generate_gaussian_matrix 里有没有做列归一化以及 DCT 基的转置有没有用对。解决在重构完成后加一句 x_rec x_rec / max(abs(x_rec(:)))先归一化再显示。如果还是不对把测量矩阵换成 dctmtx 的前 M 行这个是确定性的能排除随机矩阵的问题。4.2 PSNR 卡在 20dB 上不去现象采样率提到 0.5 了PSNR 还是 20 出头跟论文里的 35dB 差很远。原因稀疏度 K 设得不对或者分块大小太小。32×32 的块在采样率 0.5 时M512对纹理复杂的块来说测量数还是不够。另外 OMP 的 K 如果按固定值给遇到平滑块会过拟合遇到纹理块会欠拟合。解决把块大小改成 64×64采样率不变的情况下每块测量数翻倍。K 改成自适应用 DCT 系数累计能量 95% 来定。如果还不行换 CoSaMP它对 K 的估计误差容忍度更高。4.3 重构时间随图像尺寸指数增长现象512×512 的图跑了十分钟还没出结果内存占用飙升。原因没有分块直接把整图拉成 262144×1 的向量测量矩阵尺寸 131072×262144内存直接爆。或者分块了但每块都重新生成测量矩阵重复计算。解决确认代码走的是分块路径块大小不超过 64×64。测量矩阵生成一次就存下来所有块共用。如果内存还是紧张把测量矩阵改成稀疏存储MATLAB 的 sparse 函数能省不少内存。4.4 块边界出现方格状伪影现象重构图像放大看块与块之间有明显的亮度跳变像马赛克。原因分块时没有重叠每块独立重构后硬拼块边缘的系数不连续。解决启用重叠分块重叠像素数设为块大小的 1/8 到 1/4。重构后对重叠区域做加权平均权重用余弦窗中间高两边低。代码里 overlap 参数默认是 4块大小 32 时相当于 12.5% 重叠够用了。4.5 换图片后重构完全失效现象用代码自带的测试图跑得好好的换成自己的图就重构不出来。原因图片尺寸不是块大小的整数倍分块时边界处理出错。或者图片是彩色三通道代码只处理了灰度。解决读图后先做尺寸裁剪把长宽都裁成块大小的整数倍。彩色图先用 rgb2gray 转灰度或者对三个通道分别做压缩感知再合并。代码里加一句 size_check 函数自动裁剪并给出提示。5. 从重构到评估PSNR、SSIM 与采样率曲线的自动化脚本跑通单张图的重构只是第一步真正做实验的时候需要批量跑不同采样率、不同算法然后画对比曲线。这份代码里带了一个评估脚本把 PSNR 和 SSIM 的计算封装好了直接调用就行。PSNR 的计算function psnr_val calc_psnr(original, reconstructed) % 确保数据类型一致 original double(original); reconstructed double(reconstructed); mse mean((original(:) - reconstructed(:)).^2); if mse 0 psnr_val Inf; return; end max_val max(original(:)); psnr_val 10 * log10(max_val^2 / mse); end注意 max_val 取的是原图最大值不是固定的 255。如果图像归一化到 [0,1] 了max_val 就是 1PSNR 数值会跟 [0,255] 范围差 20dB 左右对比论文数据时要统一量纲。SSIM 的计算用 MATLAB 自带的 ssim 函数但要注意它默认是彩色图模式灰度图要加参数ssim_val ssim(reconstructed, original, DynamicRange, max(original(:)));批量实验的脚本结构是三层循环外层采样率从 0.1 到 0.5 步进 0.1中层算法遍历 OMP、CoSaMP、IRLS内层遍历测试图片。每次跑完把 PSNR 和 SSIM 存到矩阵里最后用 plot 画曲线。我一般会加一个计时器记录每种算法的平均重构时间。因为实际选算法的时候不只看质量还要看能不能在可接受的时间内跑完。OMP 在采样率 0.3 时单块重构大概 0.05 秒CoSaMP 是 0.12 秒IRLS 要 0.8 秒。整图 512×512 分 256 个块OMP 总共 13 秒左右IRLS 要三分钟以上。如果做实时性要求高的场景IRLS 基本不用考虑。还有一个技巧是固定随机种子。测量矩阵是随机生成的不同种子跑出来的 PSNR 会有 0.5dB 左右的波动。做对比实验时把 rng(42) 写在循环外面保证所有算法用的是同一个测量矩阵这样对比才公平。这个细节论文里经常不写但实际复现时不知道就会觉得结果对不上。从那以后我每次跑压缩感知实验都先把随机种子固定、量纲统一、块大小和重叠参数确认一遍这三步走完再开始跑数据。希望帮到你。本文还有配套的精品资源点击获取
返回列表