压缩感知重构算法之CoSaMP算法python实现
来源:互联网 发布:网络按地域分类 编辑:程序博客网 时间:2024/06/05 04:17
算法流程
算法分析
python代码
要利用python实现,电脑必须安装以下程序
- python (本文用的python版本为3.5.1)
- numpy python包(本文用的版本为1.10.4)
- scipy python包(本文用的版本为0.17.0)
- pillow python包(本文用的版本为3.1.1)
#coding:utf-8#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%# DCT基作为稀疏基,重建算法为CoSaMP算法,图像按列进行处理# 参考文献: D. Deedell andJ. Tropp, “COSAMP: Iterative Signal Recovery from#Incomplete and Inaccurate Samples,” 2008.#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#导入集成库import math# 导入所需的第三方库文件import numpy as np #对应numpy包from PIL import Image #对应pillow包#读取图像,并变成numpy类型的 arrayim = np.array(Image.open('lena.bmp'))#图片大小256*256#生成高斯随机测量矩阵sampleRate=0.5 #采样率Phi=np.random.randn(256*sampleRate,256)# 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)#CoSaMP算法函数def cs_CoSaMP(y,D): S=math.floor(y.shape[0]/4) #稀疏度 residual=y #初始化残差 pos_last=np.array([],dtype=np.int64) result=np.zeros((256)) for j in range(S): #迭代次数 product=np.fabs(np.dot(D.T,residual)) pos_temp=np.argsort(product) pos_temp=pos_temp[::-1]#反向,得到前面L个大的位置 pos_temp=pos_temp[0:2*S]#对应步骤3 pos=np.union1d(pos_temp,pos_last) result_temp=np.zeros((256)) result_temp[pos]=np.dot(np.linalg.pinv(D[:,pos]),y) pos_temp=np.argsort(np.fabs(result_temp)) pos_temp=pos_temp[::-1]#反向,得到前面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((256,256)) # 初始化稀疏系数矩阵 Theta_1d=np.dot(Phi,mat_dct_1d) #测量矩阵乘上基矩阵for i in range(256): 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, 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 imageimg=double(img);[height,width]=size(img);%------------ form the measurement matrix and base matrix ---------------Phi=randn(floor(height/2),width); % only keep one third of the original data Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(height/2),1]); % normalize each columnmat_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_cosamp(img_cs_1d(:,i),Theta_1d,height); sparse_rec_1d(:,i)=column_rec'; % sparse representationendimg_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),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)/4n=length(y); % length of measurementss=floor(n/4); % sparsity r_n=y; % initial residualssig_pos_lt=[]; % significant pos for last time iterationfor times=1: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)^(-1)*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.
0 0
- 压缩感知重构算法之CoSaMP算法python实现
- 压缩感知重构算法之压缩感知匹配追踪(CoSaMp)
- 压缩感知重构算法之压缩采样匹配追踪(CoSaMP)
- 压缩感知重构算法之压缩采样匹配追踪(CoSaMP)
- 压缩感知重构算法之压缩采样匹配追踪(CoSaMP)
- 压缩感知重构算法之OMP算法python实现
- 压缩感知重构算法之SP算法python实现
- 压缩感知重构算法之IHT算法python实现
- 压缩感知重构算法之OLS算法python实现
- 压缩感知重构算法之IRLS算法python实现
- 压缩感知重构算法之正交匹配追踪(OMP)
- 压缩感知重构算法之正交匹配追踪(OMP)
- 压缩感知重构算法之迭代硬阈值(IHT)
- 压缩感知重构算法之正交匹配追踪(OMP)
- 压缩感知重构算法之迭代软阈值(IST)
- 压缩感知重构算法之正交匹配追踪(OMP)
- CS重构之CoSaMP
- 压缩感知重构算法——SP算法
- C++实验2-模拟ATM
- 再谈数据库优化(database tuning)的真谛和误区
- Android--Intent在活动间传递数据
- hdu2825(AC自动机+状态压缩dp)
- 【开源框架】Android之史上最全最简单最有用的第三方开源库收集整理,有助于快速开发
- 压缩感知重构算法之CoSaMP算法python实现
- Android WebView中Header与Cookie应用场景深入剖析
- 如何取出Map的key和value值
- Android Studio 报错:Default activity not found
- 用户不在sudoers文件中的解决方法
- 【leetcode】Array——Find Peak Element(162)
- 【STM32】HardFault_Handler问题解决的一种办法
- 利用委托加观察者模式实现老板状态变化通知C++
- IOS学习 Block语法的使用:声明、创建、引用局部变量,使用block从后往前传值