壓縮感覺重構算法之CoSaMP算法python實作

壓縮感覺重構算法之OMP算法python實作

壓縮感覺重構算法之CoSaMP算法python實作

壓縮感覺重構算法之SP算法python實作

壓縮感覺重構算法之IHT算法python實作

壓縮感覺重構算法之OLS算法python實作

壓縮感覺重構算法之IRLS算法python實作

算法流程

算法分析

python代碼

要利用python實作，電腦必須安裝以下程式

• python （本文用的python版本為3.5.1）
• numpy python包（本文用的版本為1.10.4）
• scipy python包（本文用的版本為0.17.0）
• pillow python包（本文用的版本為3.1.1）

#coding:utf-
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# DCT基作為稀疏基，重建算法為CoSaMP算法,圖像按列進行處理
# 參考文獻: D. Deedell andJ. Tropp, “COSAMP: Iterative Signal Recovery from
#Incomplete and Inaccurate Samples,” 
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#導入內建庫
import math

# 導入所需的第三方庫檔案
import  numpy as np    #對應numpy包
from PIL import Image  #對應pillow包


#讀取圖像，并變成numpy類型的 array
im = np.array(Image.open('lena.bmp'))#圖檔大小*256

#生成高斯随機測量矩陣
sampleRate=  #采樣率
Phi=np.random.randn(*sampleRate,)
# Phi=np.random.randn(,)
# u, s, vh = np.linalg.svd(Phi)
# Phi = u[:*sampleRate,] #将測量矩陣正交化


#生成稀疏基DCT矩陣
mat_dct_1d=np.zeros((,))
v=range()
for k in range(,):  
    dct_1d=np.cos(np.dot(v,k*math.pi/))
    if k>:
        dct_1d=dct_1d-np.mean(dct_1d)
    mat_dct_1d[:,k]=dct_1d/np.linalg.norm(dct_1d)

#随機測量
img_cs_1d=np.dot(Phi,im)

#CoSaMP算法函數
def cs_CoSaMP(y,D):     
    S=math.floor(y.shape[]/)  #稀疏度    
    residual=y  #初始化殘差
    pos_last=np.array([],dtype=np.int64)
    result=np.zeros(())

    for j in range(S):  #疊代次數
        product=np.fabs(np.dot(D.T,residual))       
        pos_temp=np.argsort(product)
        pos_temp=pos_temp[::-]#反向，得到前面L個大的位置
        pos_temp=pos_temp[:*S]#對應步驟
        pos=np.union1d(pos_temp,pos_last)   

        result_temp=np.zeros(())
        result_temp[pos]=np.dot(np.linalg.pinv(D[:,pos]),y)

        pos_temp=np.argsort(np.fabs(result_temp))
        pos_temp=pos_temp[::-]#反向，得到前面L個大的位置
        result[pos_temp[:S]]=result_temp[pos_temp[:S]]
        pos_last=pos_temp
        residual=y-np.dot(D,result)

    return  result



#重建
sparse_rec_1d=np.zeros((,))   # 初始化稀疏系數矩陣    
Theta_1d=np.dot(Phi,mat_dct_1d)   #測量矩陣乘上基矩陣
for i in range():
    print('正在重建第',i,'列。。。')
    column_rec=cs_CoSaMP(img_cs_1d[:,i],Theta_1d)  #利用CoSaMP算法計算稀疏系數
    sparse_rec_1d[:,i]=column_rec;        
img_rec=np.dot(mat_dct_1d,sparse_rec_1d)          #稀疏系數乘上基矩陣

#顯示重建後的圖檔
image2=Image.fromarray(img_rec)
image2.show()

           

matlab代碼

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, [email protected], http://home.ustc.edu.cn/~roy
% Reference: D. Deedell andJ. Tropp, “COSAMP: Iterative Signal Recovery from
% Incomplete and Inaccurate Samples,” .
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------ read in the image --------------
img=imread('lena.bmp');     % testing image
img=double(img);
[height,width]=size(img);


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


mat_dct_1d=zeros(,);  % building the DCT basis (corresponding to each column)
for k=:: 
    dct_1d=cos([::]'*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(height,width);            
Theta_1d=Phi*mat_dct_1d;
for i=1:width
    column_rec=cs_cosamp(img_cs_1d(:,i),Theta_1d,height);
    sparse_rec_1d(:,i)=column_rec';           % sparse representation
end
img_rec_1d=mat_dct_1d*sparse_rec_1d;          % inverse transform


%------------ show the results --------------------
figure()
subplot(,,),imagesc(img),title('original image')
subplot(,,),imagesc(Phi),title('measurement mat')
subplot(,,),imagesc(mat_dct_1d),title('1d dct mat')
psnr = *log1(/sqrt(mean((img(:)-img_rec_1d(:)).^)));
subplot(,,),imshow(uint8(img_rec_1d));
title(strcat('PSNR=',num2str(psnr),'dB'));
disp('over')


%************************************************************************%
function hat_x=cs_cosamp(y,T_Mat,m)
% 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)/

n=length(y);                           % length of measurements
s=floor(n/);                                 % sparsity                  
r_n=y;                                 % initial residuals

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

for times=:s                          % number of iterations

    product=abs(T_Mat'*r_n);
    [val,pos]=sort(product,'descend');
    sig_pos_cr=pos(1:2*s);             % significant pos for curretn iteration

    sig_pos=union(sig_pos_cr,sig_pos_lt);

    Aug_t=T_Mat(:,sig_pos);            % current selected entries of T_Mat 

    aug_x_cr=zeros(m,1);               
    aug_x_cr(sig_pos)=(Aug_t'*Aug_t)^(-)*Aug_t'*y;  % temp recovered x (sparse)

    [val,pos]=sort(abs(aug_x_cr),'descend');

    hat_x=zeros(1,m);
    hat_x(pos(1:s))=aug_x_cr(pos(1:s));% recovered x with s sparsity  

    sig_pos_lt=pos(1:s);               % refresh the significant positions

    r_n=y-T_Mat*hat_x';
end
           

參考文獻

1、D. Deedell andJ. Tropp, “COSAMP: Iterative Signal Recovery from Incomplete and Inaccurate Samples,” 2008.

歡迎python愛好者加入：學習交流群 667279387