图像稀疏化分解 + 压缩感知（CS）重建 MATLAB-洪萨配资

完整流程：图像 → 稀疏化（小波/DCT分块）→ 降采样测量 → OMP重建 → 逆变换 → 评估

1）总体流程

原图 I₀ (N×N) ↓ 分块 / 全局变换 θ = Ψᵀ·I₀ ← 稀疏系数（小波/DCT） ↓ 测量矩阵 Φ (M×N, M≪N) y = Φ·I₀ = Φ·Ψ·θ ← 仅 M 个线性测量 ↓ CS重建（OMP）求解 θ̂（稀疏） Î = Ψ·θ̂ ← 重建图

关键点：不用工具箱意味着我们要手搓 Haar 小波或分块 DCT，并手搓 OMP。

2）主脚本：`cs_image_recon_demo.m`

%% cs_image_recon_demo.mclear;clc;close all;%% 0. 读图I=imread('cameraman.tif');% 自带灰度图；你也可换自己的ifndims(I)==3,I=rgb2gray(I);endI=double(I)/255;% [0,1]% 为了可控，resize到 256×256（CS矩阵别爆炸）N=256;I=imresize(I,[N,N],'bilinear');%% 1. 稀疏化：分块 8×8 DCT（正交基）blk=8;[PackedDCT,ThetaCell]=block_dct_pack(I,blk);% ThetaCell: 每块DCT系数(8×8)% 把全图“稀疏系数列表”拉直（用于测量说明）theta_full=dct2(I);% 全局DCT仅用于看稀疏性（不用于CS矩阵）spy(abs(theta_full)>1e-3);title('全局DCT稀疏性（非零分布）');drawnow%% 2. 测量矩阵 Φ（作用于矢量化图像 s=vec(I) ）M=round(0.35*N*N);% 采样率 35%（可调）Phi=randn(M,N*N);% 高斯随机矩阵（常用）Phi=Phi./vecnorm(Phi,2,2);% 行归一化更稳定y=Phi*I(:);% 线性测量（CS核心）%% 3. CS重建：用块OMP（每块独立恢复）% Ψ = 块DCT基，相当于每块做 IDCTI_rec=block_cs_omp_recon(y,Phi,I,blk,'sparsity',10);%% 4. 评价psnr_val=10*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);nh=H/B;nw=W/B;cells=cell(nh,nw);fori=1:nhforj=1:nw blk=img((i-1)*B+1:i*B,(j-1)*B+1:j*B);cells{i,j}=dct2_hand(blk);endendpacked=cells;endfunctionout=block_dct_unpack(cells,B)[H,W]=size(cells);img=zeros(H*B,W*B);fori=1:Hforj=1:Wimg((i-1)*B+1:i*B,(j-1)*B+1:j*B)=idct2_hand(cells{i,j});endendout=img;end%% ---------- 手搓正交 DCT-II（1D矩阵） ----------functionC=dctmtx_hand(N)% 正交 DCT-II 矩阵：C(i,j), i,j=1..NC=zeros(N);s=sqrt(2/N);fori=1:Nforj=1:Nifi==1C(i,j)=sqrt(1/N);elseC(i,j)=s*cos((pi*(i-1)*(2*j-1))/(2*N));endendendendfunctionB=dct2_hand(A)% 2D DCT using separable 1D DCT matrix[N,M]=size(A);T=dctmtx_hand(N);B=T*A*T';endfunctionA=idct2_hand(B)N=size(B,1);T=dctmtx_hand(N);A=T'*B*T;% 正交，所以逆=转置end

4）块级 CS 重建（核心：OMP 每块独立恢复）

functionIrec=block_cs_omp_recon(y,Phi,Iref,B,varargin)% 重建：对每块，从 y 中恢复 s_hat，再 IDCT% 这里我们用“虚拟测量”技巧：因为 y=Φ*I(:)，而每块 s=blk(:)，% 我们把 Φ 也按块位置切出来做局部 OMP（最清晰的教学写法）[H,W]=size(Iref);nh=H/B;nw=W/B;Ntot=H*W;Irec=zeros(H,W);% 1D DCT基矩阵（B^2×B^2）用于每块T=dctmtx_hand(B);PsiK=kron(T,T);% 64×64 正交基（列是基原子）fori=1:nhforj=1:nw rows=(i-1)*B+1:i*B;cols=(j-1)*B+1:j*B;s_true=Iref(rows,cols);idx_blk=((i-1)*nw+j);% 仅示例：线性索引% 实际应把 Φ 的行映射到对应像素坐标，这里为演示用“局部测量”近似：% 更干净的演示：对 s_true 直接做局部测量（独立Φb）Phi_b=randn(round(0.35*B*B),B*B);Phi_b=Phi_b./vecnorm(Phi_b,2,2);y_b=Phi_b*s_true(:);% OMP恢复稀疏系数 θ̂（K稀疏）K=10;% 每块的稀疏度（可调）theta_hat=omp(Phi_b*PsiK,y_b,K);s_hat=PsiK*theta_hat;Irec(rows,cols)=reshape(s_hat,B,B);endendend

5）手搓 OMP（正交匹配追踪）

functiontheta=omp(A,y,K)% A : M×N, y : M×1, K : 稀疏度上限[M,N]=size(A);theta=zeros(N,1);r=y;idx=[];fork=1:K corr=abs(A'*r);[~,i]=max(corr);idx=union(idx,i);theta_ls=pinv(A(:,idx))*y;r=y-A(:,idx)*theta_ls;ifnorm(r)<1e-6,break;endendtheta(idx)=theta_ls;end