lena.rar_lena_site:www.pudn.com_压缩感知 图像_图像压缩感知_稀疏表示

clc; clear all; tic; im=imread('lena256.jpg'); %im=rgb2gray(im); %读图像，转换为灰度图 [m n]=size(im); im2=zeros(m,n); bm=16; bn=16; %分块大小 p=0.9; %采样率 d=bm*bn; %信号长度 N=round(d*p); %测量值个数 bx=ceil(m/bm); by=ceil(n/bn); %分块数 imc=zeros(bx*bm,by*bn); %这里使用相乘，而不使用m和n，主要是有可能这个分快法，没法刚好分完整 for i=1:bx*bm %分块没有覆盖到的地方补0 for j=1:by*bn if i<=m&&j<=n imc(i,j)=im(i,j); else imc(i,j)=0; end%if end%for end%for for x=1:bx; %对每个小块处理 for y=1:by; b=imc((1+(x-1)*bm):(x*bm),(1+(y-1)*bn):(y*bn)); %从每一个小块的第一列算起，也就是分块处理 xin=reshape(b,d,1); %将所要处理的方块矩阵，改变成列向量，即将其变为了一维的信号 xin=double(xin); %输入信号 %下面四句实现了稀疏变换 A=dctmtx(d); %DCT矩阵，这个只是得到了DCT变换矩阵，而Y=TXT'则是实现了DCT变换 Phi=randn(N,d); %高斯随机矩阵（以随机矩阵作为观测矩阵） T=Phi*inv(A); %传感矩阵(随机矩阵和逆矩阵的乘积) s=Phi*xin; %测量值 %下面利用OMP算法进行图像重构 %OMP算法 hat_y=zeros(1,d); % 待重构的谱域(变换域)向量 Aug_t=[]; % 增量矩阵(初始值为空矩阵) aug_t=[]; pos_array=[]; r_n=s; times=1; while norm(r_n)>0.01 for col=1:d; product(col)=abs(T(:,col)'*r_n); % 恢复矩阵的列向量和残差的投影系数(内积值) end%for [val,pos]=max(product); [val,pos]=max(product); % 最大投影系数对应的位置 Aug_t=[Aug_t,T(:,pos)]; % 矩阵扩充 T(:,pos)=zeros(N,1); % 选中的列置零（实质上应该去掉，为了简单我把它置零） aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s; % 最小二乘,使残差最小 r_n=s-Aug_t*aug_y; % 残差 pos_array(times)=pos; % 纪录最大投影系数的位置 times=times+1; end%while hat_y(pos_array)=aug_y; % 重构的向量 hat_x=hat_y; hat_x2=inv(A)*hat_x'; %复原信号 b2=reshape(hat_x2,bm,bn); %将计算好的分块组合 imc2((1+(x-1)*bm):(x*bm),(1+(y-1)*bn):(y*bn))=b2; end%for end%for imc2=uint8(imc2); im2=imc2(1:m,1:n); mse=sum(sum(abs(im-im2).^2))/(m*n); %计算误差 toc; %计算时间 figure(1); %显示图像 imshow(im2); title(strcat('采样率=',num2str(p),'分块=',... num2str(bm),'×',num2str(bn),'计算时间=',... num2str(round(toc)),'s MSE=',num2str(round(mse*10)/10))); colormap gray;