算法遇到了问题,被卡了挺久。后来发现是个小问题,就是忘记将最终数据转换为uint8,导致图像老是无法正常显示,但是检查代码的整个逻辑过程又没有错误。这是因为在图像中,数据是以两个字节保存的。加上今天的去兼职,时间不太够,代码也是在地铁上给调试成功的。
原理我就不多说了。就是我发现压缩再复原,图像会出现很多噪点,而且压缩率越大,恢复出来的图像噪点越多。
具体原因还不是很清楚,我猜大致以下原因:
1.不可避免会产生噪点,也就是无法做到百分百完美复原。我看参考论文好像也没办法完美复原。
2.恢复矩阵或者压缩矩阵选择得不好。这个有待商榷,不过看了几篇论文,选择的矩阵都不同,不知这个会不会有一定影响。
下图,为采样率为0.75时,恢复图(左)和原图(右)的对比。
代码参考了别人写的代码,原文链接实在找不到了。代码如下:
'''
基于图像的压缩感知代码的实现(正交匹配跟踪算法)
参考文献: 任晓馨. 压缩感知贪婪匹配追踪类重建算法研究[D].
作者: kyda
时间: 2019.7.22
'''
import cv2
import numpy as np
import numpy.linalg as lg
import math
def get_matri_new_plus(A):
'''
求A矩阵的伪逆
:param A: 输入矩阵
:return A_new: 求解的伪逆矩阵
'''
A_new = A.T.dot(A)
if A_new.ndim < 2:
A_new = 1 / A_new
A_new = A_new * A.T
else:
A_new = lg.inv(A_new)
A_new = A_new.dot(A.T)
return A_new
def get_A_ba(A):
'''
求解A_ba,也就是矩阵列向量单位化
:param A: 输入矩阵
:return A_new: 返回列向量单位后的矩阵
'''
A_new = A.T
for row in range(A_new.shape[0]):
norm2 = lg.norm(A_new[row], ord=2) # 求向量的2范数
A_new[row] /= norm2 # 单位化
return A_new
def find_pos_VecInMatrix(A, v):
'''
求解某个向量在矩阵中的相应行位置
:param A: 输入矩阵
:param v: 输入向量
:return pos: 位置
'''
A_rownum = A.shape[0] # 矩阵行的数量
for pos in range(A_rownum):
if (v == A[pos]).all(): # 判断向量v和矩阵中的某一行是否相等
return pos
N = 256 # 向量长度
sampleRate = 0.75 # 采样率,也就是256长度的向量只采样64个
im = cv2.imread('lena.png') # 读取图像
im = cv2.resize(im, (N, N)) # 将图像修改成256X256的大小,这么做主要是方便计算
Phi=np.random.randn(int(sampleRate*N), N) # 生成高斯随机采样矩阵
s = Phi.dot(im) # 获得线性测量,s其实就是我们要采样的数据
#生成稀疏基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)
P = Phi.dot(mat_dct_1d.T) # 恢复矩阵(测量矩阵*正交反变换矩阵)
def cs_omp(s, mat):
'''
压缩感知主体算法
:param s: 待重构向量
:param mat: 重构矩阵
:param N: 重构向量的长度
:return: 重构完成的向量
'''
global N, K
m = math.floor(3*(s.shape[0])/4) # 算法迭代次数(m>=K)
hat_y = np.zeros(N) # 待重构的谱域(变换域)向量
Aug_t = None # 增量矩阵(初始值为空矩阵)
r_n = s # 残差值
T_tran = mat.T
A_ba = T_tran.copy() # 初始化A矩阵
A_ba1 = A_ba # 复制
pos_selected_v = []
for times in range(m):
w = A_ba.dot(r_n) # 求各列向量和残差的投影系数
pos = np.argmax(abs(w)) # 最大投影系数所对应的位置
b_k = A_ba[pos]
A_ba = np.delete(A_ba, pos, axis=0) # 选中的列删除
pos1 = find_pos_VecInMatrix(A_ba1, b_k)
pos_selected_v.append(pos1) # 纪录最大投影系数的位置
if Aug_t is None: # 将最大贡献向量(未单位化的)加入Aug_t矩阵
Aug_t = T_tran[pos1]
else:
Aug_t = np.c_[Aug_t, T_tran[pos1]]
A_new_plus = get_matri_new_plus(Aug_t)
aug_y = A_new_plus.dot(s) # 解最小二乘解
r_n = s - Aug_t.dot(aug_y) # 求残差
for item in range(m): # 重构的谱域向量
pos = pos_selected_v[item]
x_r = aug_y[item]
try:
hat_y[pos] = x_r
except np.ComplexWarning:
pass
hat_x = mat_dct_1d.T.dot(hat_y) # 做逆傅里叶变换重构得到时域信号
return hat_x
restruct_im = np.zeros([256, 256, 3]) # 重构的图像
for x in range(256):
for y in range(3):
restruct_s = cs_omp(s[:, x].T[y], P) # 重构向量
restruct_im[x].T[y] = restruct_s # 将重构的向量写入图像
restruct_im = np.array(restruct_im, dtype=np.uint8) # 这里一定要将矩阵的数据形式转换为uint8,不然图片无法正常显示
cv2.imshow('im_re', restruct_im)
cv2.imshow('im_ss', im)
cv2.waitKey(0)
欢迎关注我的公众号,可以交流问题: