周报：2021-3-1~2021-3-7

本周工作：

  • 利用小波基优化了前期OMP算法的程序，2D图像重建结果得到显著提升

  • 实现了CASSI的TwIST算法，并应用到了Columbia University的CAVE多光谱数据集上，图像重建各波段光谱恢复良好，但图像分辨率较低

基于小波基的2D图像CS重建(OMP算法)

最近下载到了一些CS的算法代码，也看了一些博客，学习到了如何对2D图像在稀疏基作用下进行稀疏表示，并用OMP算法进行重建的步骤。将其对我前期的OMP算法代码进行了改进，发现重构效果也不差，噪声更小。

参考Blog：

原始在DCT基下重建的效果如下：

image-20210309095247973

在DWT基下的重建效果如下：

image-20210309095401761

其中生成DWT基的代码如下（并附上我自己改的Python版）：

Matlab code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
%  程序作者：沙威，香港大学电气电子工程学系，wsha@eee.hku.hk
%  参考文献：小波分析理论与MATLAB R2007实现，葛哲学，沙威，第20章  小波变换在矩阵方程求解中的应用（沙威、陈明生编写）.

%  构造正交小波变换矩阵，图像大小N*N，N=2^P，P是整数。

function ww=DWT(N)

[h,g]= wfilters('sym8','d');       %  分解低通和高通滤波器

% N=256;                           %  矩阵维数(大小为2的整数幂次)
L=length(h);                       %  滤波器长度
rank_max=log2(N);                  %  最大层数
rank_min=double(int8(log2(L)))+1;  %  最小层数
ww=1;   %  预处理矩阵

%  矩阵构造
for jj=rank_min:rank_max
    
    nn=2^jj;
    
    %  构造向量
    p1_0=sparse([h,zeros(1,nn-L)]);
    p2_0=sparse([g,zeros(1,nn-L)]);
    
    %  向量圆周移位
    for ii=1:nn/2
        p1(ii,:)=circshift(p1_0',2*(ii-1))';
        p2(ii,:)=circshift(p2_0',2*(ii-1))';
    end
    
    %  构造正交矩阵
    w1=[p1;p2];
    mm=2^rank_max-length(w1);
    w=sparse([w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)]);
    ww=ww*w;
    
    clear p1;clear p2;
end
Python code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
import numpy as np
import pywt

def DWT(N):
    sym8 = pywt.Wavelet('haar') # 选sym8小波
    lod,hid = sym8.dec_lo, sym8.dec_hi # 分解低通和高通
    L = len(lod) # 小波长度
    # print(L)
    rank_max = np.log2(N) # 最大层数
    rank_min = float(np.int8(np.log2(L))) + 1 # 最小层数
    res = 1 # result
    for jj in range(int(rank_min),int(rank_max)+1):
        nn = 2 ** jj
        p1_0 = np.r_[lod, np.zeros(nn-L)]
        p2_0 = np.r_[hid, np.zeros(nn-L)]
        p1 = np.zeros((nn//2, nn))
        p2 = np.zeros((nn//2, nn))
        for ii in range(nn//2):
            p1[ii, :] = np.roll(p1_0.T, 2*ii, axis=0).T
            p2[ii, :] = np.roll(p2_0.T, 2*ii, axis=0).T
        w1 = np.r_[p1,p2]
        mm = int(2**rank_max - len(w1))
        uphalf = np.c_[w1, np.zeros((len(w1),mm))]
        downhalf = np.c_[np.zeros((mm,len(w1))), np.eye(mm)]
        w = np.r_[uphalf,downhalf]
        res = np.dot(res, w)
        
    return res

重建代码也附上：

Matlab (Gauss Random Matrix) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
%  本程序实现图像LENA的压缩传感
%  程序作者：沙威，香港大学电气电子工程学系，wsha@eee.hku.hk
%  算法采用正交匹配法，参考文献 Joel A. Tropp and Anna C. Gilbert 
%  Signal Recovery From Random Measurements Via Orthogonal Matching
%  Pursuit，IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 53, NO. 12,
%  DECEMBER 2007.
%  该程序没有经过任何优化

function Wavelet_OMP

clc;clear

%  读文件
X=imread('lena256.bmp');
X=double(X);
[a,b]=size(X);

%  小波变换矩阵生成
ww=DWT(a);

%  小波变换让图像稀疏化（注意该步骤会耗费时间，但是会增大稀疏度）
X1=ww*sparse(X)*ww';
X1=full(X1);



%  随机矩阵生成
M=190;
R= GaussMtx(M,a);
    
%  测量
Y=R*X1;

%  OMP算法
X2=zeros(a,b);  %  恢复矩阵
for i=1:b  %  列循环       
    rec=omp(Y(:,i),R,a);
    X2(:,i)=rec;
end


%  原始图像
figure(1);
imshow(uint8(X));


%  变换图像   %小波变换后的图像
figure(2);
imshow(uint8(X1));


%测量值图像
figure(3);
imshow(uint8(Y));


%  压缩传感恢复的图像
figure(4);
X3=ww'*sparse(X2)*ww;  %  小波反变换
X3=full(X3);
imshow(uint8(X3));


%  误差(PSNR)
errorx=sum(sum(abs(X3-X).^2));       %  MSE误差
psnr=10*log10(255*255/(errorx/a/b))   %  PSNR
 dNMSE = nmse(uint8(X),uint8(X3))   %NMSE误差

%  OMP的函数
%  s-测量；T-观测矩阵；N-向量大小
function hat_y=omp(s,T,N)

Size=size(T);                                     %  观测矩阵大小
M=Size(1);                                        %  测量
hat_y=zeros(1,N);                                 %  待重构的谱域(变换域)向量                     
Aug_t=[];                                         %  增量矩阵(初始值为空矩阵)
r_n=s;                                            %  残差值

for times=1:M/4;                                  %  迭代次数(稀疏度是测量的1/4)
    for col=1:N;                                  %  恢复矩阵的所有列向量
        product(col)=abs(T(:,col)'*r_n);          %  恢复矩阵的列向量和残差的投影系数(内积值) 
    end
    [val,pos]=max(product);                       %  最大投影系数对应的位置
    Aug_t=[Aug_t,T(:,pos)];                       %  矩阵扩充
    T(:,pos)=zeros(M,1);                          %  选中的列置零（实质上应该去掉，为了简单我把它置零）
    aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s;           %  最小二乘,使残差最小
    r_n=s-Aug_t*aug_y;                            %  残差
    pos_array(times)=pos;                         %  纪录最大投影系数的位置
    
    if (norm(r_n)<9)                              %  残差足够小
        break;
    end
end
hat_y(pos_array)=aug_y;                           %  重构的向量
Python (Gauss Random Matrix) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
import random
import numpy as np
import time
import sys
import cv2
import matplotlib.pyplot as plt
import scipy.optimize as spopt
import scipy.fftpack as spfft
import scipy.ndimage as spimg
from scipy import linalg
from DWT import DWT

# OMP algorithm
def OMP(y,A,K):
	col = A.shape[1]
	residual = y
	ind = []
	for i in range(K):
		prod = np.fabs(np.dot(A.T,residual))
		pos = np.argmax(prod)
		ind.append(pos)
		a = np.dot(np.linalg.pinv(A[:,ind]),y)
		residual = y-np.dot(A[:,ind],a)
	
	Res = np.zeros((col,))
	Res[ind] = a
	# print(Res.shape)
	return Res

img = cv2.imread(r"E:\Graduation-Project\OMP\lena.jpg")
img = img[128:128+256, 128:128+256, 0]

# cv2.imshow("img",img)
# cv2.waitKey(0)

# radom gaussian matrix ----> undersample
sample_rate = 0.7
N = 256
Phi = np.random.randn(int(sample_rate*N),N)
# Phi = linalg.orth(Phi.T).T 
# Phi /= np.linalg.norm(Phi)

# plt.imshow(Phi)
# plt.colorbar()
# plt.show()
# exit(0)


# sym8 Wavelet ----> sparse base
Psi = DWT(N)
img = np.dot(np.dot(Psi,img),Psi.T)

# plt.imshow(Psi)
# plt.colorbar()
# plt.show()
# exit(0)


# measurement ----> undersample image
measure_mat = np.dot(Phi,img)

# K coefficient for N colomns
sparse_coe = np.zeros((N,N))

# θ = Φ ψ
# Theta = np.dot(Phi,Psi)

Theta = Phi


time_consume = []
# OMP for every colomn
for k in range(30,40,10):
	st = time.perf_counter()
	for i in range(N):
		sparse_coe[:,i] = OMP(measure_mat[:,i],Theta,k)
	en = time.perf_counter()
	img_rec = np.dot(np.dot(Psi.T,sparse_coe),Psi)
	# img_rec /= img_rec.max()
	# img_rec *= 255
	img_rec = img_rec.astype(np.uint8)
	# print(img_rec.shape)
	# print(img.shape)
	# cv2.imshow("img",img)
	cv2.imshow("rate=%.1f,k=%d"%(sample_rate,k),img_rec)
	sys.stdout.writelines(f"rate={sample_rate:.1f},k={k} runtime:{en-st}s\n")
	time_consume.append(en-st)
cv2.waitKey(0)

CASSI系统TwIST算法对高光谱数据集的模拟试验

在Github上下载了对应的原版CASSI源码后，花了两天时间将其改为Python版，无奈跑得太慢，只好屈服于MATLAB。。。

重建步骤梳理

1. 读取图像、构造编码孔径、预处理（shearing & smash）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
%%
close all;clear all;clc;

% 256 x 256 x 31 的多光谱图像
R = 256;
C = 256;
L = 31;

coded_aperture = randi([0.,1.],[R,C]); % 随机01编码孔径
coded_aperture = repmat(coded_aperture,[1,1,L]); % 补充光谱维度
coded_aperture_shift = shift(coded_aperture, 1); % shearing

% 读取CAVE数据集并整合成mat
path = '..\CAVE_dataset\glass_tiles_ms\glass_tiles_ms\glass_tiles_ms_';

for i=1:31;
    tp = imread([path num2str(i,'%02d') '.png']);
    origin_img(:,:,i) = double(tp(128:128+R-1,128:128+C-1)); % 转成double好算点积
end
origin_img = origin_img/max(origin_img(:)); % 对原图归一化
save(['orig_data_1.mat'],'origin_img');

origin_img_shift = shift(origin_img, 1); % shearing
y = sum(origin_img_shift.*coded_aperture_shift, 3); % smash
y = y/max(y(:)); % 归一化
figure;imagesc(y);colormap gray;axis image;colorbar;set(gcf,'color','w'); % 作图（模拟的detector探测结果）

其中，shift函数如下：

1
2
3
4
5
6
function y = shift(x, step)
    [R,C,L] = size(x);
    for i=1:L
        y(:,i*step:i*step+C-1,i) = x(:,:,i);
    end
end

2. 设定正则化参数和迭代次数

1
2
3
4
5
6
%%
%Regularization parameter
tau = 0.05; % L1正则化参数，越大解越稀疏

%Number of TwIST main loop iterations
maxiterations = 50; % 最大迭代次数

3. 构建CASSI系统的传递模型及其转置

1
2
3
4
%%
%Define the functional forms for the Forward and Forward^T models
A = @(x) Rfuntwist(x,coded_aperture_shift);
AT = @(x) RTfuntwist(x,coded_aperture_shift);

其中两个函数分别为正向传递的模型和它的转置：

1
2
3
4
5
6
7
8
9
10
11
12
13
function y = Rfuntwist(f,C)
    [n1,n2,m] = size(C);
    f = reshape(f,[n1,n2,m]); % reshape成数据立方体格式
    gp = f.*C; % 对每个光谱段图像编码
    y = sum(gp,3); % 沿光谱维度叠加投影到CCD上
end

function f = RTfuntwist(y,C)
    [n1,n2,m] = size(C);
    y=reshape(y,n1,n2);
    yp = repmat(y,[1,1,m]); % 沿光谱维度复制m份
    f = yp.*C; % 编码
end

4. 定义去噪函数 (放弃软阈值，采用全变分去噪)

1
2
3
4
5
6
%%Preparing to run TwIST
%Number of Chambolle iterations to perform for total variation minimization
piter = 4; % 全变分去噪迭代次数

Psi = @(x,th) mycalltoTVnew(x,th,piter); % 全变分去噪
Phi = @(x) TVnormspectralimaging(x); % 定义去噪问题的正则项为光谱图像梯度和

全变分去噪采用Chambolle projection's algorithm

5. TwIST重建

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
tolA = 1e-8; % 收敛标准
%Run TwIST
[x_twist_orig,dummy,obj_twist,...
    times_twist,dummy,mse_twist]= ...
    TwISTmod(y,A,tau,...
    'AT', AT, ...
    'Psi', Psi, ...
    'Phi',Phi, ...
    'Initialization',2,...
    'Monotone',1,...
    'StopCriterion',1,...
    'MaxIterA',maxiterations,...
    'ToleranceA',tolA,...
    'Debias',0,...
    'Verbose', 1);

6. 图像复原

1
2
3
4
5
6
%x_twist is the estimate of the data cube. Prepare it to be displayed as a colored cube.
x_twist = x_twist_orig;
x_twist = x_twist.*(x_twist>=0); % 去掉负值
x_twist = x_twist/max(x_twist(:));
x_twist = shift_back(x_twist,1);
save(['rec_result_1.mat'],'x_twist');

其中，shift_back函数如下：

1
2
3
4
5
6
function y = shift_back(x, step)
    [R,C,L] = size(x);
    for i=1:L
        y(:,:,i) = x(:,i*step:i*step+R-1,i);
    end
end

重建结果

reconstruction_res1_img

image-20210310192447973

项目记录
周报

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议 ，转载请注明出处！