壓縮感知重構演算法之IRLS演算法python實現
阿新 • • 發佈:2018-12-30
IRLS(iteratively reweighted least squares)演算法
(本文給出的程式碼未進行優化,只是為了說明演算法流程 ,所以執行速度不是很快)
IRLS(iteratively reweighted least squares)演算法是壓縮感知重建演算法當中的一個基本演算法。主要是為了解決
本文采用的程式碼是加入權重之後的
上式中的權重
這樣上面的最優化問題可以求解得到:
其中
程式碼
要利用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 = 0.7 #取樣率
Phi = np.random.randn(256, 256)
u, s, vh = np.linalg.svd(Phi)
Phi = u[:256*sampleRate,] #將測量矩陣正交化
#生成稀疏基DCT矩陣
mat_dct_1d=np.zeros((256,256))
v=range(256)
for k in range(0,256):
dct_1d=np.cos(np.dot(v,k*math.pi/256))
if k>0:
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[0])/4)
hat_x_tp=np.dot(T_Mat.T ,y)
epsilong=1
p=1 # solution for l-norm p
times=1
while (epsilong>10e-9) and (times<L): #迭代次數
weight=(hat_x_tp**2+epsilong)**(p/2-1)
Q_Mat=np.diag(1/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,2) < np.sqrt(epsilong)/100):
epsilong = epsilong/10
hat_x_tp=hat_x
times=times+1
return hat_x
#重建
sparse_rec_1d=np.zeros((256,256)) # 初始化稀疏係數矩陣
Theta_1d=np.dot(Phi,mat_dct_1d) #測量矩陣乘上基矩陣
for i in range(256):
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, roy@mail.ustc.edu.cn, http://home.ustc.edu.cn/~roy
% Reference: R. Chartrand and W. Yin, “Iteratively Reweighted Algorithms
% for Compressed Sensing,” 2008.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------ 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/3),width); % only keep one third of the original data
Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(height/3),1]); % normalize each column
mat_dct_1d=zeros(256,256); % 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(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(1)
subplot(2,2,1),imagesc(img),title('original image')
subplot(2,2,2),imagesc(Phi),title('measurement mat')
subplot(2,2,3),imagesc(mat_dct_1d),title('1d dct mat')
psnr = 20*log10(255/sqrt(mean((img(:)-img_rec_1d(:)).^2)))
subplot(2,2,4),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)/4
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