简介这份资源面向数字图像处理、计算机视觉与医学成像方向的学习者和研究者提供二维图像压缩感知重构算法的完整MATLAB实现帮助读者从代码层面理解CS理论如何用远少于像素数的测量值重建图像。压缩包共26个文件以18个m脚本为核心覆盖测量矩阵生成、图像采样、ISTA/FISTA等重构算法与性能评估流程另含7个fig结果图和1个bmp测试图像整体约1.59MB结构紧凑便于逐模块调试。已有4295人学习下载说明其在教学与实验场景中具备一定参考价值。读者可据此对比高斯、伯努利等不同测量矩阵的重构质量调整松弛因子与阈值参数并借助MSE、PSNR指标和可视化结果验证算法效果也可作为向视频压缩、医学影像、遥感分析及深度学习结合方向拓展的实验起点。1. 二维图像压缩感知重构为什么它值得你花一晚上跑通一张 512×512 的灰度图按奈奎斯特采样要存 26 万个像素点可如果它本身是稀疏的——比如医学 CT、遥感切片、工业缺陷图——那大部分采样其实是在浪费存储和带宽。压缩感知Compressive Sensing, CS的核心主张就是只要信号在某个变换域稀疏就能用远少于奈奎斯特率的观测值把它恢复出来。二维图像压缩感知重构算法干的就是这件事把一张图先做稀疏变换、再降维观测、最后用重构算法从少量观测里把原图“猜”回来。这套东西适合谁做 MATLAB 图像处理大作业的学生、想给嵌入式成像链路省带宽的工程师、以及需要验证某种稀疏先验是否成立的研究者。它不挑硬件一台装了 MATLAB 的普通笔记本就能跑但它挑参数测量率、稀疏基、迭代步长任何一个设错重构结果就是一团噪声。下面我把从稀疏表示到重构代码的完整路径拆开讲内含可直接抄的 MATLAB 代码也会把那些让我翻过车的参数边界说清楚。2. 稀疏表示与观测矩阵重构能不能成这一步就定了压缩感知的数学骨架是 y Φx ΦΨs其中 x 是原图拉成列向量Ψ 是稀疏基s 是稀疏系数Φ 是观测矩阵y 是观测值。重构就是从已知的 y 和 ΦΨ 反推 s再还原 x。整条链路能不能成取决于两件事Ψ 能不能让图像足够稀疏以及 Φ 和 Ψ 的乘积是否满足 RIP有限等距性质。这一章先把这两个前提落到可操作的 MATLAB 代码上。2.1 选稀疏基DCT 还是小波别凭感觉自然图像在小波域通常比 DCT 域更稀疏因为小波能同时抓住平滑区域和边缘。但小波基的构造比 DCT 麻烦如果你只是做课程作业或快速验证DCT 是更稳的起点。下面这段代码把图像转成列向量并构造 DCT 稀疏基% 读取图像并归一化 img imread(lena_gray.png); if size(img,3) 3 img rgb2gray(img); % 彩色图先转灰度 end img im2double(img); N numel(img); % 像素总数如 512*512262144 % 构造一维 DCT 稀疏基对列向量作用 Psi dctmtx(N); % N×N 正交 DCT 矩阵 s_true Psi * img(:); % 真实稀疏系数仅用于对比dctmtx(N)生成的是正交 DCT 矩阵Psi * img(:)把图像投影到 DCT 域。注意这里对整幅图做一维 DCT等价于二维可分离 DCT 的列向量形式好处是后面写观测模型时维度统一。参数上N会直接决定矩阵大小512×512 图对应 262144×262144 的Psi内存吃紧时不要显式构造改用dct函数按块处理。2.2 观测矩阵高斯随机为什么是默认选择观测矩阵 Φ 的尺寸是 M×NM 是观测数M/N 就是测量率。常见选择有高斯随机矩阵、伯努利矩阵、部分傅里叶矩阵。高斯随机矩阵和绝大多数稀疏基都不相关RIP 满足概率高是验证阶段最省心的选择M round(0.5 * N); % 测量率 50%先从这个值试 rng(42); % 固定随机种子保证可复现 Phi randn(M, N) / sqrt(M); % 高斯随机观测矩阵按 1/sqrt(M) 归一化 % 合成观测值 y Phi * img(:); % 注意这里直接对像素观测等价于 Phi*Psi*s A Phi * Psi; % 感知矩阵重构算法真正要用的randn(M,N)/sqrt(M)的归一化很关键它让观测值的能量量级和原图接近否则后续迭代的步长要重新调。A Phi * Psi是感知矩阵所有重构算法都围绕它展开。测量率 0.5 是入门值实际能压到多少取决于图像稀疏度我后面会给出测量率扫描的验证方法。提示固定rng种子不是可选项。压缩感知重构对观测矩阵敏感不固定种子你这次跑出 28dB、下次跑出 22dB根本没法判断是算法问题还是随机波动。2.3 从观测到重构先明确你要解的是什么重构的本质是求解 s argmin ||s||_1 s.t. y A s或者带噪声时改成 argmin ||s||_1 s.t. ||y - A s||_2 ≤ ε。这是一个凸优化问题但维度高直接调linprog会慢到没法用。实际工程里用的是迭代软阈值类算法比如 ISTA、FISTA、TwIST它们每步只做矩阵乘法和软阈值MATLAB 里几十行就能实现。下一章我会把 FISTA 的完整代码写出来并解释步长和阈值参数怎么定。3. 用 FISTA 在 MATLAB 里跑通重构完整代码与参数标定FISTA快速迭代收缩阈值算法是 ISTA 的加速版核心改进是在软阈值步骤前加了一个动量项收敛速度从 O(1/k) 提升到 O(1/k²)。对于 512×512 的图像ISTA 可能要上千次迭代FISTA 通常两三百次就能到可接受的重构质量。这一章给出可直接运行的完整函数并逐段说明参数含义。3.1 FISTA 主循环每一步在做什么function [x_rec, s_rec, err_hist] fista_cs(y, A, N, lambda, max_iter, tol) % FISTA 求解 min 0.5*||y - A*s||^2 lambda*||s||_1 % 输入y 观测值(M×1)A 感知矩阵(M×N)N 信号长度 % lambda 正则化参数max_iter 最大迭代tol 收敛阈值 % 输出x_rec 重构图像s_rec 稀疏系数err_hist 每步残差 L norm(A)^2; % Lipschitz 常数用谱范数平方估计 step 1 / L; % 步长必须 1/L否则发散 s zeros(N, 1); % 稀疏系数初值 z s; % 动量变量 t 1; err_hist zeros(max_iter, 1); for k 1:max_iter grad A * (A * z - y); % 梯度A(Az - y) s_new soft_threshold(z - step * grad, lambda * step); t_new (1 sqrt(1 4 * t^2)) / 2; z s_new ((t - 1) / t_new) * (s_new - s); err_hist(k) norm(y - A * s_new); if norm(s_new - s) / max(norm(s), eps) tol s s_new; break; end s s_new; t t_new; end s_rec s; x_rec Psi_reconstruct(s); % 用稀疏基反变换回图像域 end function out soft_threshold(in, th) out sign(in) .* max(abs(in) - th, 0); endL norm(A)^2估计 Lipschitz 常数步长取1/L是保证收敛的硬约束设大了会直接发散成 NaN。soft_threshold是软阈值算子把绝对值小于阈值的系数压到零这是 L1 正则产生稀疏性的机制。动量项((t-1)/t_new)*(s_new - s)就是 FISTA 比 ISTA 快的原因它让迭代沿着历史方向加速。err_hist记录残差用来判断是否收敛。3.2 正则化参数 lambda 怎么定别用默认值lambda控制稀疏性和保真度的权衡。太大重构过度平滑、丢细节太小噪声被当成信号重构出来。理论上的上界是lambda_max norm(A * y, inf)超过这个值解会全零。实用做法是从0.01 * lambda_max开始扫lambda_max norm(A * y, inf); lambda_list [0.001, 0.005, 0.01, 0.05, 0.1] * lambda_max; psnr_list zeros(size(lambda_list)); for i 1:length(lambda_list) [x_rec, ~, ~] fista_cs(y, A, N, lambda_list(i), 300, 1e-5); x_rec reshape(x_rec, size(img)); psnr_list(i) psnr(x_rec, img); fprintf(lambda%.4f*max, PSNR%.2f dB\n, ... lambda_list(i)/lambda_max, psnr_list(i)); end跑完这个循环你会看到 PSNR 随 lambda 先升后降峰值对应的就是当前测量率下的较优值。我一般会把这个扫描固定成脚本换图像、换测量率都重跑一遍因为最优 lambda 和这两者都相关。别偷懒用固定值这是最容易翻车的地方之一。3.3 测量率扫描找到你的应用能压到多少测量率 M/N 决定了压缩程度也决定了重构质量上限。做方案评估时至少要扫 0.1 到 0.5 这几个点rates [0.1, 0.2, 0.3, 0.4, 0.5]; for r rates M round(r * N); Phi randn(M, N) / sqrt(M); A Phi * Psi; y Phi * img(:); lambda 0.01 * norm(A * y, inf); [x_rec, ~, ~] fista_cs(y, A, N, lambda, 300, 1e-5); x_rec reshape(x_rec, size(img)); fprintf(rate%.1f, PSNR%.2f dB\n, r, psnr(x_rec, img)); end经验上自然灰度图在 DCT 基下测量率 0.3 以上能到 30dB 左右0.2 以下会明显糊。如果你的应用要求 35dB 以上要么提高测量率要么换更稀疏的基小波、学习字典要么上更好的重构算法TV 正则、BM3D-CS。这个扫描结果就是你判断方案可行性的依据。4. 重构质量上不去先排查这五类问题压缩感知重构的调试比普通图像处理烦因为失败模式多、现象相似。下面五条是我实际踩过的坑按「现象 → 原因 → 解决」写遇到问题按顺序排查。4.1 重构结果是纯噪声或 NaN现象x_rec全是噪点或者迭代几步后出现 NaN。原因步长超过 1/L或者A的归一化不对。解决先打印norm(A)^2确认step 1/L确实小于它再检查Phi有没有除以sqrt(M)漏掉这步会让A的范数偏大步长相对就超了。我习惯在函数开头加一句assert(step * norm(A)^2 1.01)直接拦住这类错误。4.2 PSNR 卡在 20dB 上不去现象不管怎么调 lambdaPSNR 就在 20dB 附近。原因测量率太低或者稀疏基和图像不匹配。解决先把测量率提到 0.5 看能不能上去能上去说明是测量率问题还上不去就换小波基试。DCT 对纹理丰富的图像稀疏性差换成wavedec2做小波变换往往能涨 3 到 5dB。4.3 迭代不收敛残差震荡现象err_hist上下跳不单调下降。原因动量项系数计算错误或者tol设得太小导致在数值噪声里打转。解决检查t_new (1sqrt(14*t^2))/2有没有写错把tol从 1e-5 放宽到 1e-4通常就稳了。FISTA 的残差不保证单调轻微震荡正常但大幅震荡一定是参数问题。4.4 内存爆掉矩阵构造失败现象512×512 图跑dctmtx(N)时报内存不足。原因N262144Psi是 N×N 稠密矩阵约 512GB。解决不要显式构造Psi改用函数句柄按块做 DCT或者把图像分块如 32×32 块分别做压缩感知块间用重叠保留去块效应。分块是工程上更实际的做法。4.5 换一张图结果差很多现象在 Lena 上 32dB换到纹理图就掉到 24dB。原因压缩感知的重构质量依赖图像稀疏度纹理图在固定基下不够稀疏。解决要么针对图像类型选基遥感图用轮廓波、医学图用学习字典要么接受测量率要提高。没有一种基能对所有图都好这是压缩感知的固有边界不是代码 bug。5. 进阶把重构从「能跑」推到「能用」的三个技巧前面跑通的是标准 FISTA DCT 基线实际项目里往往还要再推一把。第一个技巧是换 TV 正则它比 L1 更擅长保边缘代价是每步要解一个子问题可以用 Chambolle 投影近似MATLAB 里多写二十行PSNR 通常能再涨 1 到 2dB。第二个技巧是 warm start先用低测量率的结果做初值再逐步提高测量率精修比一步到位收敛快。第三个是块稀疏把图像分成 16×16 块每块独立重构后拼回来内存友好且天然并行适合大图。验证重构是否真的可用别只看 PSNR。我习惯同时看三样PSNR 看整体、SSIM 看结构保真、以及把残差图abs(x_rec - img)放大看残差如果集中在边缘说明算法在丢细节如果均匀分布说明是噪声底。下面这个小脚本一次输出三个指标x_rec reshape(x_rec, size(img)); fprintf(PSNR%.2f dB, SSIM%.4f\n, psnr(x_rec, img), ssim(x_rec, img)); figure; subplot(1,3,1); imshow(img); title(原图); subplot(1,3,2); imshow(x_rec); title(重构); subplot(1,3,3); imshow(abs(x_rec - img), []); title(残差);最后说个习惯每次改参数前先把当前最优的lambda、测量率、PSNR 记到一个表格里跑完一轮再对比。压缩感知的参数空间是三维的测量率、lambda、迭代数凭记忆调参一定会乱。我早期不信邪来回改了十几轮才发现最优解就在第一次扫描的附近白折腾一晚上。希望帮到你。本文还有配套的精品资源点击获取