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;
