壓縮感知重構演算法之OLS演算法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: [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*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
相關推薦
壓縮感知重構演算法之OLS演算法python實現
Orthogonal Least Squares (OLS)演算法流程 實驗 要利用python實現,電腦必須安裝以下程式 python (本文用的python版本為3.5.1) numpy python包(本文用
演算法之氣泡排序-python實現
大家好,歡迎收看我的文章,如果感覺我的文章能對您有所幫助,您可以點選關注我,您的支援就是我堅持創作的動力 氣泡排序演算法 比如有6個數: [22,44,33,55,66,77]從大到小排序,對相鄰的兩位進行比較 第一輪 第一次比較: 44,22,33,55
壓縮感知重構演算法之IRLS演算法python實現
IRLS(iteratively reweighted least squares)演算法 (本文給出的程式碼未進行優化,只是為了說明演算法流程 ,所以執行速度不是很快) IRLS(iteratively reweighted least squar
壓縮感知重構演算法之OMP演算法python實現
本文主要簡單介紹了利用python程式碼實現壓縮感知的過程。 壓縮感知簡介 【具體可以參考這篇文章】 假設一維訊號x長度為N,稀疏度為K。Φ 為大小M×N矩陣(M<<N)。y=Φ×x為長度M的一維測量值。壓縮感知問題就是已知測量值y和測
壓縮感知重構演算法之CoSaMP演算法python實現
演算法流程 演算法分析 python程式碼 要利用python實現,電腦必須安裝以下程式 python (本文用的python版本為3.5.1) numpy python包(本文用的版本為1.10.4) scipy python
壓縮感知重構演算法之壓縮取樣匹配追蹤(CoSaMP)
題目:壓縮感知重構演算法之壓縮取樣匹配追蹤(CoSaMP) 壓縮取樣匹配追蹤(CompressiveSampling MP)是D. Needell繼ROMP之後提出的又一個具有較大影響力的重構演算法。CoSaMP也是對OMP的一種改進,每次迭代選擇多個原子,除了原子的選擇
[轉]壓縮感知重構算法之分段正交匹配追蹤(StOMP)
參數配置 組成 jaf second red [1] figure nor 拉伸 分段正交匹配追蹤(StagewiseOMP)或者翻譯為逐步正交匹配追蹤,它是OMP另一種改進算法,每次叠代可以選擇多個原子。此算法的輸入參數中沒有信號稀疏度K,因此相比於ROMP及CoSaMP
python 資料結構與演算法之歸併演算法
def merge_sort(alist): n=len(alist) if n<=1: return alist mid=n//2 left_list =merge_sort(alist[:mid]) right_list =mer
小白之KMP演算法詳解及python實現
在看子串匹配問題的時候,書上的關於KMP的演算法的介紹總是理解不了。看了一遍程式碼總是很快的忘掉,後來決定好好分解一下KMP演算法,算是給自己加深印象。 ------------------------- 分割線-------------------------------
聚類演算法之DBSCAN演算法之二:高維資料剪枝應用NQ-DBSCAN
一、經典DBSCAN的不足 1.由於“維度災難”問題,應用高維資料效果不佳 2.執行時間在尋找每個點的最近鄰和密度計算,複雜度是O(n2)。當d>=3時,由於BCP等數學問題出現,時間複雜度會急劇上升到Ω(n的四分之三次方)。 二、DBSCAN在高維資料的改進 目前的研究有
聚類演算法之DBSCAN演算法之一:經典DBSCAN
DBSCAN是基於密度空間的聚類演算法,與KMeans演算法不同,它不需要確定聚類的數量,而是基於資料推測聚類的數目,它能夠針對任意形狀產生聚類。 1.epsilon-neighborhood epsoiln-neighborhood(簡稱e-nbhd)可理解為密度空間,表示半徑為e
WF曲速未來:區塊鏈核心演算法之Paxos演算法
Paxos演算法解決的問題是在一個可能發生訊息可能會延遲、丟失、重複的分散式系統中如何就某個值達成一致,保證不論發生以上任何異常,都不會破壞決議的一致性。 先帶你會看一下libpaxos3的程式碼: 第一步獲取和編譯LibPaxos3所需的基本步驟: 執行示例  
梅森演算法生成隨機數的Python實現
import time class Util(object): def __init__(self): self.index = 624 self.MT = [0] * 624 def inter(self,t): return (0xFFFFFFFF
資料探勘演算法之K_means演算法
轉載地址:https://blog.csdn.net/baimafujinji/article/details/50570824 聚類是將相似物件歸到同一個簇中的方法,這有點像全自動分類。簇內的物件越相似,聚類的效果越好。支援向量機、神經網路所討論的分類問題都是有監督的學習方式
連結分析演算法之 HITS演算法
分享一下我老師大神的人工智慧教程!零基礎,通俗易懂!http://blog.csdn.net/jiangjunshow 也歡迎大家轉載本篇文章。分享知識,造福人民,實現我們中華民族偉大復興!  
文字相似度bm25演算法的原理以及Python實現(jupyter notebook)
今天我們一起來學習一下自然語言處理中的bm25演算法,bm25演算法是常見的用來計算query和文章相關度的相似度的。其實這個演算法的原理很簡單,就是將需要計算的query分詞成w1,w2,…,wn,然後求出每一個詞和文章的相關度,最後將這些相關度進行累加,最終就可以的得到文字相似度計算
字串匹配演算法之KMP演算法詳情
package demo; /* 字串匹配演算法 */ public class StringKMP { //找出從第一個字元開始 子串T在主串S的第一個位置 如果沒有則返回-1 public static int index(String S, String T)
在Object-C中學習資料結構與演算法之排序演算法
筆者在學習資料結構與演算法時,嘗試著將排序演算法以動畫的形式呈現出來更加方便理解記憶,本文配合Demo 在Object-C中學習資料結構與演算法之排序演算法閱讀更佳。 目錄 選擇排序 氣泡排序 插入排序 快速排序 雙路快速排序 三路快速排序 堆排序 總結與收穫
SVM演算法 K-means的python實現
arg argument of the maximum/minimum arg max f(x): 當f(x)取最大值時,x的取值 arg min f(x):當f(x)取最小值時,x的取值 s.t.是subject to (such that)的縮
大資料探勘領域十大經典演算法之—CART演算法(附程式碼)
簡介 CART與C4.5類似,是決策樹演算法的一種。此外,常見的決策樹演算法還有ID3,這三者的不同之處在於特徵的劃分: ID3:特徵劃分基於資訊增益 C4.5:特徵劃分基於資訊增益比 CART:特徵劃分基於基尼指數 基本思想 CART假設決策樹是二叉樹,