1. 程式人生 > >壓縮感知重構演算法之OLS演算法python實現

壓縮感知重構演算法之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假設決策樹是二叉樹,