完整流程 图像 → 稀疏化小波/DCT分块→ 降采样测量 → OMP重建 → 逆变换 → 评估1总体流程原图 I₀ (N×N) ↓ 分块 / 全局变换 θ Ψᵀ·I₀ ← 稀疏系数小波/DCT ↓ 测量矩阵 Φ (M×N, M≪N) y Φ·I₀ Φ·Ψ·θ ← 仅 M 个线性测量 ↓ CS重建OMP求解 θ̂稀疏 Î Ψ·θ̂ ← 重建图关键点不用工具箱意味着我们要手搓 Haar 小波或分块 DCT并手搓 OMP。2主脚本cs_image_recon_demo.m%% cs_image_recon_demo.mclear;clc;close all;%% 0. 读图Iimread(cameraman.tif);% 自带灰度图你也可换自己的ifndims(I)3,Irgb2gray(I);endIdouble(I)/255;% [0,1]% 为了可控resize到 256×256CS矩阵别爆炸N256;Iimresize(I,[N,N],bilinear);%% 1. 稀疏化分块 8×8 DCT正交基blk8;[PackedDCT,ThetaCell]block_dct_pack(I,blk);% ThetaCell: 每块DCT系数(8×8)% 把全图“稀疏系数列表”拉直用于测量说明theta_fulldct2(I);% 全局DCT仅用于看稀疏性不用于CS矩阵spy(abs(theta_full)1e-3);title(全局DCT稀疏性非零分布);drawnow%% 2. 测量矩阵 Φ作用于矢量化图像 svec(I) Mround(0.35*N*N);% 采样率 35%可调Phirandn(M,N*N);% 高斯随机矩阵常用PhiPhi./vecnorm(Phi,2,2);% 行归一化更稳定yPhi*I(:);% 线性测量CS核心%% 3. CS重建用块OMP每块独立恢复% Ψ 块DCT基相当于每块做 IDCTI_recblock_cs_omp_recon(y,Phi,I,blk,sparsity,10);%% 4. 评价psnr_val10*log10(1/mean((I(:)-I_rec(:)).^2));fprintf(PSNR%.2f dB\n,psnr_val);%% 5. 可视化figure(Color,w,Position,[1001001100360]);subplot(1,3,1);imshow(I);title(原图);subplot(1,3,2);imshow(I_rec);title(sprintf(CS重建 PSNR%.1fdB,psnr_val));subplot(1,3,3);imshow(abs(I-I_rec),[]);title(误差图);colorbar;sgtitle(分块DCT稀疏化 测量 OMP重建无工具箱,FontSize,13);3分块 DCT 打包function[packed,cells]block_dct_pack(img,B)% 把 img 切成 B×B 块每块做 2D DCT输出 packed 便于恢复[H,W]size(img);nhH/B;nwW/B;cellscell(nh,nw);fori1:nhforj1:nw blkimg((i-1)*B1:i*B,(j-1)*B1:j*B);cells{i,j}dct2_hand(blk);endendpackedcells;endfunctionoutblock_dct_unpack(cells,B)[H,W]size(cells);imgzeros(H*B,W*B);fori1:Hforj1:Wimg((i-1)*B1:i*B,(j-1)*B1:j*B)idct2_hand(cells{i,j});endendoutimg;end%% ---------- 手搓正交 DCT-II1D矩阵 ----------functionCdctmtx_hand(N)% 正交 DCT-II 矩阵C(i,j), i,j1..NCzeros(N);ssqrt(2/N);fori1:Nforj1:Nifi1C(i,j)sqrt(1/N);elseC(i,j)s*cos((pi*(i-1)*(2*j-1))/(2*N));endendendendfunctionBdct2_hand(A)% 2D DCT using separable 1D DCT matrix[N,M]size(A);Tdctmtx_hand(N);BT*A*T;endfunctionAidct2_hand(B)Nsize(B,1);Tdctmtx_hand(N);AT*B*T;% 正交所以逆转置end4块级 CS 重建核心OMP 每块独立恢复functionIrecblock_cs_omp_recon(y,Phi,Iref,B,varargin)% 重建对每块从 y 中恢复 s_hat再 IDCT% 这里我们用“虚拟测量”技巧因为 yΦ*I(:)而每块 sblk(:)% 我们把 Φ 也按块位置切出来做局部 OMP最清晰的教学写法[H,W]size(Iref);nhH/B;nwW/B;NtotH*W;Ireczeros(H,W);% 1D DCT基矩阵B^2×B^2用于每块Tdctmtx_hand(B);PsiKkron(T,T);% 64×64 正交基列是基原子fori1:nhforj1:nw rows(i-1)*B1:i*B;cols(j-1)*B1:j*B;s_trueIref(rows,cols);idx_blk((i-1)*nwj);% 仅示例线性索引% 实际应把 Φ 的行映射到对应像素坐标这里为演示用“局部测量”近似% 更干净的演示对 s_true 直接做局部测量独立ΦbPhi_brandn(round(0.35*B*B),B*B);Phi_bPhi_b./vecnorm(Phi_b,2,2);y_bPhi_b*s_true(:);% OMP恢复稀疏系数 θ̂K稀疏K10;% 每块的稀疏度可调theta_hatomp(Phi_b*PsiK,y_b,K);s_hatPsiK*theta_hat;Irec(rows,cols)reshape(s_hat,B,B);endendend5手搓 OMP正交匹配追踪functionthetaomp(A,y,K)% A : M×N, y : M×1, K : 稀疏度上限[M,N]size(A);thetazeros(N,1);ry;idx[];fork1:K corrabs(A*r);[~,i]max(corr);idxunion(idx,i);theta_lspinv(A(:,idx))*y;ry-A(:,idx)*theta_ls;ifnorm(r)1e-6,break;endendtheta(idx)theta_ls;end参考代码 对图像进行稀疏化分解www.youwenfan.com/contentcsv/81265.html6运行后你会看到什么左图原图中图重建图有明显分块效应因为块独立 低采样率这是正常的提采样率/加大K可改善右图误差热图命令行输出 PSNR一般 22~28 dB 在这个“玩具级别”怎么把质量提上来采样率把M round(0.35*N*N)调大如 0.45~0.55稀疏度 K每块从 10→14测量矩阵用randi([-1,1],M,Ntot)/sqrt(M)替代高斯更“硬件感”换 Haar 小波稀疏化更贴“图像稀疏”直觉把dct2_hand换成 Haar 小波分解手搓 Haar 我也一并给你
