压缩感知重构算法之OLS算法python实现
2016-04-05 21:00
423 查看
压缩感知重构算法之OMP算法python实现
压缩感知重构算法之CoSaMP算法python实现
压缩感知重构算法之SP算法python实现
压缩感知重构算法之IHT算法python实现
压缩感知重构算法之OLS算法python实现
压缩感知重构算法之IRLS算法python实现
python (本文用的python版本为3.5.1)
numpy python包(本文用的版本为1.10.4)
scipy python包(本文用的版本为0.17.0)
pillow python包(本文用的版本为3.1.1)
2、Hashemi A, Vikalo H. Sparse Linear Regression via Generalized Orthogonal Least-Squares[J]. arXiv preprint arXiv:1602.06916, 2016.
欢迎python爱好者加入:学习交流群 667279387
压缩感知重构算法之CoSaMP算法python实现
压缩感知重构算法之SP算法python实现
压缩感知重构算法之IHT算法python实现
压缩感知重构算法之OLS算法python实现
压缩感知重构算法之IRLS算法python实现
Orthogonal Least Squares (OLS)算法流程
实验
要利用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:ncuzhangben@qq.com, #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ''' #导入集成库 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*sampleRate,256) #生成稀疏基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) #OLS算法函数(未完成待修改) def cs_ols(y,D): L=math.floor(3*(y.shape[0])/4) residual=y #初始化残差 index=np.zeros((L),dtype=int) for i in range(L): index[i]= -1 result=np.zeros((256)) for j in range(L): #迭代次数 pos_temp=np.argsort(np.fabs(np.dot(np.linalg.pinv(D[:,pos]),y)))#对应步骤2 product=np.fabs(np.dot(D.T,residual)) pos=np.argmax(product) #最大投影系数对应的位置 index[j]=pos my=np.linalg.pinv(D[:,index>=0]) a=np.dot(my,y) residual=y-np.dot(D[:,index>=0],a) result[index>=0]=a 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_ols(img_cs_1d[:,i],Theta_1d) #利用OMP算法计算稀疏系数 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_OLS() %------------ 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(0.5*height),width); % only keep one third of the original data Phi = Phi./repmat(sqrt(sum(Phi.^2,1)),[floor(0.5*height),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); % height*width的0矩阵 Theta_1d=Phi*mat_dct_1d;%测量矩阵乘上基矩阵 for i=1:width column_rec=OLS(img_cs_1d(:,i),Theta_1d);%算法的目的是得到稀疏系数 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),imshow(uint8(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')); %************************************************************************% function [s, residual] = OLS(y, A, err) % Orthogonal Least Squares [1] for Sparse Signal Reconstruction % Input % A = N X d dimensional measurment matrix % y = N dimensional observation vector % m = sparsity of the underlying signal % Output % s = estimated sparse signal % r = residual % [1] T. Blumensath, M. E. Davies; "On the Difference Between Orthogonal % Matching Pursuit and Orthogonal Least Squares", manuscript 2007 if nargin < 3 err = 1e-5; end n1=length(y); m=floor(3*n1/4); s = zeros(size(A,2),1); r(:,1) = y; L = []; Psi = []; normA=(sum(A.^2,1)).^0.5; NI = 1:size(A,2); for i = 2:m+1 DR = A'*r(:,i-1); [v, I] = max(abs(DR(NI))./normA(NI)'); INI = NI(I); L = [L' INI']'; NI(I)=[]; Psi = A(:,L); x = Psi\y; yApprox = Psi*x; r(:,i) = y - yApprox; if norm(r(:,end)) < err break; end end s(L) = x; residual = r(:,end);
参考文章
1、Blumensath T, Davies M E. On the difference between orthogonal matching pursuit and orthogonal least squares[J]. 2007.2、Hashemi A, Vikalo H. Sparse Linear Regression via Generalized Orthogonal Least-Squares[J]. arXiv preprint arXiv:1602.06916, 2016.
欢迎python爱好者加入:学习交流群 667279387
相关文章推荐
- Python与Selenium初试
- Python MySQL操作
- python读取txt文件最后一行(文件大+文件小)
- Python 对新浪微博的博文元素 (Word, Screen Name)的频率分析
- Python列表和元组
- GDB 编译--with-python unusable python问题
- python get方法
- 关于Python 中的 map()函数
- Python优雅编程技巧
- 数据分类K—means 算法的python代码实现
- Python~~~关键字~~~
- python time 与datetime之间的区别与联系
- Python~迭代
- python time 与datetime之间的区别与联系
- Python描述符(descriptor)解密
- Python~切片Slice
- python version 2.7 required,which was not found in the registry
- 理解Python中的with…as…语法
- Python如何安装egg组件
- python常用内置模块,执行系统命令的模块