压缩感知重构图像的几种方法

这个程序是基于压缩感知理论的图像重构算法，采用离散余弦变换将原始图像稀疏化，在采用高斯随机矩阵进行采样，最后使用MP算法、OMP算法、coamap算法等几种算法重构

function Demo_CS_CoSaMP（）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the DCT basis is selected as the sparse representation dictionary
% instead of seting the whole image as a vector I process the image in the
% fashion of column-by-column so as to reduce the complexity.

% Author: Chengfu Huo roy@mail.ustc.edu.cn http://home.ustc.edu.cn/~roy
% Reference: D. Deedell andJ. Tropp “COSAMP: Iterative Signal Recovery from
% Incomplete and Inaccurate Samples” 2008.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------ read in the image --------------
img=imread（‘lena.bmp‘）;     % testing image
img=double（img）;
[heightwidth]=size（img）;


%------------ form the measurement matrix and base matrix ---------------
Phi=randn（floor（height/3）width）;  % only keep one third of the original data  
Phi = Phi./repmat（sqrt（sum（Phi.^21））[floor（height/3）1]）; % normalize each column


mat_dct_1d=zeros（256256）;  % building the DCT basis （corresponding to each column）
for k=0:1:255 
    dct_1d=cos（[0:1:255]‘*k*pi/256）;
    if k>0
        dct_1d=dct_1d-mean（dct_1d）; 
    end;
    mat_dct_1d（:k+1）=dct_1d/norm（dct_1d）;
end


%--------- projection ---------
img_cs_1d=Phi*img;          % treat each column as a independent signal


%-------- recover using omp ------------
sparse_rec_1d=zeros（heightwidth）;            
Theta_1d=Phi*mat_dct_1d;
for i=1:width
    column_rec=cs_cosamp（img_cs_1d（:i）Theta_1dheight）;
    sparse_rec_1d（:i）=column_rec‘;           % sparse representation
end
img_rec_1d=mat_dct_1d*sparse_rec_1d;          % inverse transform


%------------ show the results --------------------
figure（1）
subplot（221）imagesc（img）title（‘original image‘）
subplot（222）imagesc（Phi）title（‘measurement mat‘）
subplot（223）imagesc（mat_dct_1d）title（‘1d dct mat‘）
psnr = 20*log10（255/sqrt（mean（（img（:）-img_rec_1d（:））.^2）））
subplot（224）imagesc（img_rec_1d）title（strcat（‘1d rec img ‘num2str（psnr）‘dB‘））

disp（‘over‘）


%************************************************************************%
function hat_x=cs_cosamp（yT_Matm）
% y=T_Mat*x T_Mat is n-by-m
% y - measurements
% T_Mat - combination of random matrix and sparse representation basis
% m - size of the original signal
% the sparsity is length（y）/4

n=length（y）;                           % length of measurements
s=floor（n/4）;                                 % sparsity                  
r_n=y;                                 % initial residuals

sig_pos_lt=[];                         % significant pos for last time iteration

for times=1:s                          % number of iterations
    
    product=abs（T_Mat‘*r_n）;
    [valpos]=sort（product‘descend‘）;
    sig_pos_cr=pos（1:2*s）;             % significant pos for curretn iteration
    
    sig_pos=union（sig_pos_crsig_pos_lt）;
    
    Aug_t=T_Mat（:sig_pos）;            % current selected entries 

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件     428181  2010-07-23 10:17  CS_recovery_algorithms_OMP_SP_IHT\CS Recovery Algorithms.pdf

     文件       3456  2010-07-16 14:45  CS_recovery_algorithms_OMP_SP_IHT\software\Demo_CS_CoSaMP.m

     文件       8278  2010-07-15 10:02  CS_recovery_algorithms_OMP_SP_IHT\software\Demo_CS_GBP.m

     文件       2904  2010-07-22 14:16  CS_recovery_algorithms_OMP_SP_IHT\software\Demo_CS_IHT.m

     文件       2903  2010-07-15 10:02  CS_recovery_algorithms_OMP_SP_IHT\software\Demo_CS_IRLS.m

     文件       3340  2010-07-16 14:47  CS_recovery_algorithms_OMP_SP_IHT\software\Demo_CS_OMP.m

     文件       3434  2010-07-15 10:03  CS_recovery_algorithms_OMP_SP_IHT\software\Demo_CS_SP.m

     文件      66614  2010-05-11 20:29  CS_recovery_algorithms_OMP_SP_IHT\software\lena.bmp

     目录          0  2010-07-23 10:18  CS_recovery_algorithms_OMP_SP_IHT\software

     目录          0  2010-07-23 10:18  CS_recovery_algorithms_OMP_SP_IHT

----------- ---------  ---------- -----  ----

               519110                    10