压缩感知重构算法之IRLS算法python实现

压缩感知重构算法之OMP算法python实现

压缩感知重构算法之CoSaMP算法python实现

压缩感知重构算法之SP算法python实现

压缩感知重构算法之IHT算法python实现

压缩感知重构算法之OLS算法python实现

压缩感知重构算法之IRLS算法python实现

IRLS（iteratively reweighted least squares）算法

（本文给出的代码未进行优化，只是为了说明算法流程 ，所以运行速度不是很快）

IRLS（iteratively reweighted least squares）算法是压缩感知重建算法当中的一个基本算法。主要是为了解决

minu||u||pp, subject to Φu=b

本文采用的代码是加入权重之后的

minu∑i=1Nwiu2i, subject to Φu=b

上式中的权重 wi 是根据前面一次 ui−1 计算得到的，具体的计算公式为：

wi=|u(n−1)i|p−2

这样上面的最优化问题可以求解得到：

u(n)=QnΦT(ΦQnΦT)−1b

其中 Qn 是一个对角矩阵，具体值从 1/wi=|u(n−1)i|2−p 得到。详细具体的解释请看参考文献1。

代码

要利用python实现，电脑必须安装以下程序

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

python代码

#coding: utf-8
'''
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# DCT基作为稀疏基，重建算法为OMP算法 ，图像按列进行处理
# email:[email protected], 
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
'''
#导入集成库
import math

# 导入所需的第三方库文件
import  numpy as np    #对应numpy包
from PIL import Image  #对应pillow包


#读取图像，并变成numpy类型的 array
im = np.array(Image.open('lena.bmp'))
#print (im.shape, im.dtype)uint8

#生成高斯随机测量矩阵
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)

#IRLS算法函数
def cs_irls(y,T_Mat):   
    L=math.floor((y.shape[])/)
    hat_x_tp=np.dot(T_Mat.T ,y)
    epsilong=
    p= # solution for l-norm p
    times=
    while (epsilong>) and (times<L):  #迭代次数
        weight=(hat_x_tp**+epsilong)**(p/-)
        Q_Mat=np.diag(/weight)
        #hat_x=Q_Mat*T_Mat'*inv(T_Mat*Q_Mat*T_Mat')*y
        temp=np.dot(np.dot(T_Mat,Q_Mat),T_Mat.T)
        temp=np.dot(np.dot(Q_Mat,T_Mat.T),np.linalg.inv(temp))
        hat_x=np.dot(temp,y)        
        if(np.linalg.norm(hat_x-hat_x_tp,) < np.sqrt(epsilong)/):
            epsilong = epsilong/
        hat_x_tp=hat_x
        times=times+
    return hat_x



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

#显示重建后的图片
image2=Image.fromarray(img_rec)
image2.show()


           

matlab代码

%matlab版本用的R2010b
function Demo_CS_IRLS()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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: R. Chartrand and W. Yin, “Iteratively Reweighted Algorithms
% for Compressed Sensing,” .
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------ 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_irls(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(,,),imagesc(img_rec_1d),title(strcat('1d rec img ',num2str(psnr),'dB'))

%****************************************
function hat_x=cs_irls(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)/

hat_x_tp=T_Mat'*y; 
epsilong=1;
p=1;                             % solution for l-norm p
times=1;
while (epsilong>10e-9) && (times<length(y)/4)
    weight=(hat_x_tp.^2+epsilong).^(p/2-1); 
    Q_Mat=diag(1./weight,0); 
    hat_x=Q_Mat*T_Mat'*inv(T_Mat*Q_Mat*T_Mat')*y;
        if(norm(hat_x-hat_x_tp,2) < sqrt(epsilong)/100)
        epsilong=epsilong/10;
    end
        hat_x_tp=hat_x;
    times=times+1;
end
           

参考文献

1、R. Chartrand and W. Yin, “Iteratively Reweighted Algorithms for Compressed Sensing,” 2008.

欢迎python爱好者加入：学习交流群 667279387