OMP步骤及可运行的python代码（个人理解版）

第一步：初始化残差，以及已经使用过的列向量空间；
第二步：对列向量进行初始化，也就是归一化；（因为求贡献时，要用到内积的绝对值除以列向量的模，因此先归一化后，就可以直接求内积的绝对值，为什么是内积？内积相当于残差在该列向量方向的投影长度，也就是残差对这个列向量的贡献大小）
第三步：在未用的列向量空间中，找到对残差贡献最大的列向量
第四步：将找到的列向量加入已经用列向量空间
第五步：用上述的列向量空间对b进行求最小值的操作，从而得到x的值（未使用的列向量上的x值为0），从而得到新一轮的残差（相当于把贡献最大的方向去掉，找第二，第三，。。。，贡献的方向）
第六步：重复上述步骤，直至残差达到给定的阈值

此时我们认为，结束时，除了已使用的列向量空间对应的x不为0，其他的未使用列向量空间对应的x都为0，因此，算法执行循环几次，得到的x的稀疏度就是多少。
OMP与MP最大的不同在于图中STEP-4这步骤，MP更新的是一个X_K，而OMP是将所有已经选过的列空间对应的x都进行更新；
用最小二乘的时候就已经保证了残差和每次更新的支撑向量集（已使用的列向量空间）正交了。最小二乘的几何意义就是正交投影，在这里投影空间就是支撑向量集张成的空间（已使用的列向量空间），残差y-Ax和投影空间正交，也就是说和投影空间的每个向量正交。因此才有了正交（0）。

处理二维图像的过程跟处理一维信号类似，OMP把二维图像的每一列进行类似一维信号的操作，再将所有处理完的列进行拼接成图像的操作；为什么图片稀疏恢复效果会比较差呢？个人觉得：1，图片一定每列都具有相同的稀疏性（代码里面取K=50，相当与默认在稀疏基作用下，每一列只有50个非零元）2.稀疏集的选取，可能在DCT这种稀疏基下，不具备稀疏性； 3. 这种对图像每一列进行一维信号重构的方法对图像的结构进行了破坏，因此回恢复效果比较差；（上述观点是个人意见，如有错误，麻烦指出，谢谢）

def cs_omp(y,Phi,K):    residual=y  #初始化残差(M,N) = Phi.shapeindex=np.zeros(N,dtype=int)for i in range(N): #第i列被选中就是1，未选中就是-1index[i]= -1result=np.zeros((N,1))for j in range(K):  #迭代次数product=np.fabs(np.dot(Phi.T,residual))pos=np.argmax(product)  #最大投影系数对应的位置        index[pos]=1 #对应的位置取1my=np.linalg.pinv(Phi[:,index>=0]) #广义逆矩阵（伪逆）(A^T*A)^(-1)*A^T，最小二乘法得到的解和伪逆矩阵是等价的。         a=np.dot(my,y) #最小二乘 residual=y-np.dot(Phi[:,index>=0],a)result[index>=0]=aCandidate = np.where(index>=0) #返回所有选中的列return  result, Candidate

#生成稀疏基DCT矩阵
mat_dct_1d=np.zeros((N,N))
v=range(N)
for k in range(0,N):  dct_1d=np.cos(np.dot(v,k*math.pi/N))if k>0:dct_1d=dct_1d-np.mean(dct_1d)mat_dct_1d[:,k]=dct_1d/np.linalg.norm(dct_1d)
print(mat_dct_1d)

from PIL import Image
#读取图像，并变成numpy类型的 array
im = Image.open('img.jpg').convert('L') #模式“RGB”转换为其他不同模式#图片大小256*256，为灰度图像，每个像素用8个bit表示，0表示黑，255表示白，其他数字表示不同的灰度。
#转换公式：L = R * 299/1000 + G * 587/1000+ B * 114/1000。
plt.imshow(im,cmap='gray')   #cmap参数: 为调整显示颜色  gray为黑白色，加_r取反为白黑色
N = 256
M = 128
K = 50
im = np.array(im)
# 观测矩阵
Phi=np.random.randn(M,N)/np.sqrt(M)
# 随机测量
img_cs_1d=np.dot(Phi,im)   #（M*N）

# 重建
sparse_rec_1d=np.zeros((N,N))   # 初始化稀疏系数矩阵
Theta_1d=np.dot(Phi,mat_dct_1d)   #测量矩阵乘上稀疏基矩阵
for i in range(N):if i%32==0:print('正在重建第',i,'列。')y=np.reshape(img_cs_1d[:,i],(M,1))      #将图像的第I列进行重构column_rec, Candidate=cs_omp(y,Theta_1d,K) #利用OMP算法计算稀疏系数，column_rec为重构的稀疏系数的值，不是x的值，x=稀疏基与稀疏系数的值#做内积x_pre = np.reshape(column_rec,(N))     #常用于矩阵规格变换，将矩阵转换为特定的行和列的矩阵sparse_rec_1d[:,i]=x_pre

img_rec=np.dot(mat_dct_1d,sparse_rec_1d)          #稀疏系数乘上基矩阵，对图像进行重构
#显示重建后的图片
img_pre=Image.fromarray(img_rec)
plt.imshow(img_pre,cmap='gray')
error = np.linalg.norm(img_rec-im)/np.linalg.norm(im)     #采样率为M/N=50%,说明在余弦变换下或许不满足稀疏度假设，也就是在此稀疏基下不稀疏
error

OMP步骤及可运行的python代码（个人理解版）相关推荐

王者荣耀——bat批处理文件，自动刷金币版（脱胎于30行Python代码刷金币版），Windows双击即可运行！
参考<30行Python代码刷王者荣耀金币>:https://segmentfault.com/a/1190000012520431 1.源代码以下是源代码部分,全部复制到文本文档, 用 ...
终于从树堆里爬出来了——堆排序(基于二叉树)基本思想、步骤、复杂度及python代码，欢迎交流
欢迎关注,敬请点赞! 树堆逃生记一.动图演示二.思路分析 1. 相关概念 2. 基本思想 3. 步骤 [步骤一] 构造初始堆 [步骤二] 将堆顶元素与末尾元素进行交换,使末尾元素最大. [步骤三] ...
详解200行Python代码实现控制台版2048【总有一款坑适合你】【超详细】
跟着实验楼学习了2048的Python实现,先丢个地址 200行Python代码实现2048 我接触Python时间不长,只了解一些基本的语法和容器,在学习的过程中遇到不少问题,这里做一个记录. cu ...
【20210906】让实验室服务器运行本地python代码
从零开始配置实验室电脑的python环境 1. 电脑信息 2. 电脑环境配置 (1)Pycharm (2)anaconda (3)配置Anaconda+pycharm环境 3. 服务器环境配置小结 ...
※ 用一个代码同时运行其他python代码
※ 用一个代码执行指定python程序本文主要介绍一个简单的小知识,即利用一个代码去执行所有你所写好的代码程序. 直接开工! import os os.system("python 执行的 ...
用来表示python代码块的是什么_三分钟带你用简单的Python代码深入理解Python中的元类...
互联网的数据爆炸式的增长,而利用 Python 爬虫我们可以获取大量有价值的数据类也是对象在理解元类前,需要先掌握Python中的类.在大多数编程语言中,类就是一组描述如何生成对象的代码段.在Py ...
python调用短信宝API发送短信（附python代码易理解）
原版API如下:接口说明_马上使用更好的短信服务-短信宝官网 (smsbao.com) 直接上代码复制过去就行 import requestscontent = str("[短信宝]您的验 ...
python代码写好了怎么运行不了-python代码可以直接运行吗 Python写了代码如何运行...
先下载python,然后打开命令行,输入 python 你的代码文件名. 有python代码怎么编成可执行的exe程序? 如果可以能否帮小编做成可执行的exe程序儿女情长什么的,真的很影响小编行走江湖 ...

OMP步骤及可运行的python代码（个人理解版）

OMP步骤及可运行的python代码（个人理解版）相关推荐

最新文章

热门文章