Image Super-Resolution via Sparse Representation——基于稀疏表示的超分辨率重建

　　经典超分辨率重建论文，基于稀疏表示。下面首先介绍稀疏表示，然后介绍论文的基本思想和算法优化过程，最后使用python进行实验。

稀疏表示

　　稀疏表示是指，使用过完备字典中少量向量的线性组合来表示某个元素。过完备字典是一个列数大于行数的行满秩矩阵，也就是说，它的列向量有无数种线性组合来表达列向量空间中的任意点。由于它的列数通常远大于行数，可以使用占比很小的列向量来表示特定的向量，我们称这种表示为稀疏表示。

　　那么如何获得这个字典呢？它在特定的任务下有特定的取值。和炼丹类似，我们先要用大量数据来训练这个矩阵，让它提取出能稀疏表示这些数据的特征，进而拥有稀疏表示其它相似数据的能力。

　　训练过完备字典的过程称为稀疏编码。设训练数据集为矩阵$X=(x_1,x_2,...,x_n)in R^{m imes n}$，待训练矩阵为$Ain R^{m imes K}$，矩阵对每一数据的表示权重为$alpha = (alpha_1,alpha_2,...,alpha_n)in R^{K imes n}$。进行如下优化：

egin{align} minlimits_{A,alpha}|alpha|_0;;;s.t. ; |Aalpha - X|_2^2le epsilon end{align}

　　容易理解，约束重建向量与原始向量差异的同时，优化表示权重$alpha$的稀疏性。通过优化获得所需的过完备字典后，我们就可以用它来稀疏表示新的数据了。对于某个向量$y$，我们可以进行类似的优化来获得它在字典中的稀疏表示：

egin{align} minlimits_{alpha}|alpha|_0;;;s.t. ; |Aalpha - y|_2^2le epsilon end{align}

　　因为零范数的离散性，以上优化是NP-难问题，无法用正常的优化算法解决。所以通常会对上式进行调整，以方便优化。而在线性约束下，对一范数的优化有稀疏性（点击链接理解），因此可以转换为对一范数的优化。然后根据拉格朗日乘子法（据说如此，但我觉得这个转换并不等价），不等式约束可以移入优化中：

egin{align} minlimits_{alpha} lambda |alpha|_1 + |Aalpha - y|_2^2 end{align}

　　同样，$(1)$式也可以进行类似的转换。

　　以上就是稀疏表示的流程，看到这个提取特征然后重建的过程，我们会联想到主成分分析（PCA）。PCA能使我们方便地找到一组“完备”基向量，但是这里我们要做的是找到一组“过完备”的基向量来表示输入向量。过完备基的好处是它们能更有效地找出隐含在输入数据内部的结构与模式。与PCA不同，对于过完备基来说，系数$alpha$不再由输入向量$y$单独确定。因此，在稀疏编码算法中，我们另加一个评判标准“稀疏性”来解决因过完备而导致的退化（degeneracy）问题。

　　上面这段是百度百科原话。我觉得，把过完备字典与神经网络进行对比，可以把这个待训练的很“宽”的矩阵看做参数量很大的网络。我们知道参数量大而训练数据不充足的时候模型很容易过拟合，为了防止过拟合就要加上正则项，以使参数能专注于学习更有共性的特征。我们可以把上面的稀疏性看做正则化来理解，使字典的列向量能表达一些更有“特点”的信息。

论文原理及实现流程

基本思想

　　在训练阶段，论文同时对LR训练集$Y= (y_1,y_2,...,y_n)$和对应的HR训练集$X = (x_1,x_2,...,x_n)$分别训练两个过完备字典$D_l,D_h$，使得LR数据$y_i$和它对应的HR数据$x_i$能以相同的稀疏编码$alpha_i$分别被$D_l$和$D_h$表示。也就是

$left{ egin{aligned} &D_lalpha_i approx y_i \ &D_halpha_i approx x_i end{aligned} ight.$

　　在测试阶段，我们已经有了训练好的$D_l$和相对应的$D_h$。对于测试图像$y_t$，首先通过优化获得$y_t$在$D_l$中的稀疏表示$alpha_t$，此时有$D_lalpha_t approx y_t$。然后用这个表示通过$D_h$映射出对应的SR图像，即$hat{x}_t=D_halpha_t$。

训练过程

　　训练过程就是训练上述的过完备字典对。因为性能的因素，我们不可能直接对整张图进行稀疏编码，论文是将图像分为方形的区块（patch）进行编码的。因此，用于训练的成对数据不是整张的LR-HR图像对，而是所有图像分割开来的区块对。现在把LR训练集的所有区块表示为$Y=(y_1,y_2,...,y_n)in R^{M imes n}$，相应的HR训练集区块表示为$X = (x_1,x_2,...,x_n)in R^{N imes n}$。如果放大倍数为$t$倍，则有$N=t^2M$。由于论文中是先将LR图像通过Bicubic放大到与HR图像一致，则$N=M$。

　　优化式很直观：

egin{align} minlimits_{D_l,D_h,alpha}frac{1}{N}|X - D_halpha|_2^2+frac{1}{M}|Y - D_lalpha|_2^2 + lambda |alpha|_1 end{align}

　　其中$D_lin R^{M imes K},D_hin R^{N imes K}$，分别表示待训练的LR和HR字典，$K$表示字典的原子数；$alphain R^{K imes n}$为表示矩阵；$lambda$为平衡稀疏性和重建一致性的系数。一二两项惩罚使用相同表示的重建差异，第三项用来优化表示的稀疏性。

　　论文说同时优化$D_c$和$alpha$非凸，但是固定其中一个变量，然后对另一个变量的优化则是凸的。因此可以将凸优化交替进行，最终达到一个局部最优点。既然只是达到局部最优，我们也可以选择使用梯度下降。

推理过程

　　在获得$D_l$和$D_h$后，就可以用它们对LR图像进行重建了。论文采用扫描的方式，一个区块一个区块从上到下、从左到右对测试图像进行超分辨率重建。另外，为了使相邻区块之间能相互匹配，防止颜色上的冲突，前后两个区块之间设有重叠部分，在重建一个区块时，重叠部分要和上一个区块一致。具体优化方式如下。

　　首先将测试图像$y$按顺序划分为$m$个区块$(p_1,p_2,...,p_m)$，再设区块在$D_l$中的表示为$(alpha_1,alpha_2,...,alpha_m)$。然后按顺序对所有区块执行优化，对于第$i$个区块，优化式如下：

egin{align} minlimits_{alpha_i} lambda|alpha_i|_1 + |FD_lalpha_i - Fp_i|_2^2+ |PD_halpha_i - w|_2^2 end{align}

　　其中$w$表示已重建区块和当前区块的重叠部分，$P$表示选择当前区块与之前区块的重叠部分的操作。至于$F$，论文说是一个线性提取器，使图像中感知效果更好。论文实验时用的是一阶、二阶导数滤波器，但没说清楚具体如何操作，我实验就没有用。

　　式子意义很明显，第一项保证表示的稀疏性，第二项优化原始LR图像的重建一致性，第三项优化SR图像相邻区块重叠部分的一致性。获得SR图像所有区块的稀疏表示后，左乘$D_h$，然后将区块拼接起来，就是最终的SR图像了。

　　除了以上步骤以外，论文还多了一个所谓全局重建约束。用什么反向投影法，通过迭代让SR图像退化后能和原始图像更相似。由于说的很不清楚，这里就不加了，而且我觉得这不是这篇论文的主要内容。

　　另外，论文在执行推理过程之前，先将原图减去自身元素的均值，以使模型能更专注于纹理的重建，在重建完的SR图像上再加回这个均值。但是这个策略只在推理章节提了一句，在训练$D_l,D_h$时是否使用标准化的图像并没有说明。

实验与分析

　　实验使用LSUN数据集中的bedroom作为训练集，从中选取1024张长宽都大于256像素的图片，居中裁剪至256x256，获得HR训练集。然后对HR使用Bicubic缩小4倍至64x64，再放大回256x256，获得LR训练集（也就是四倍超分辨率）。将LR/HR区块划分为16x16大小，则每张图片都可被划分为16x16个区块。由于算力有限，定义$D_l$和$D_h$的原子（列向量）数为2560。图像通道的选择与论文中一致，先将图像色彩空间转换到YCrCb，然后单取亮度Y通道用于训练与推理，因此$D_l,D_h$原子维度为$16 imes 16 imes 1=256$。另外，对于每个图像区块中像素的亮度值，都会减去区块的亮度平均值，使得模型能专注于纹理的训练与生成，这一点也是与论文一致的。

训练

　　综上，对于$(4)$式中的各个参数，有：

egin{equation} left{ egin{aligned} &X in R^{(16 imes 16) imes (1024 imes 16^2)}\ &Y in R^{(16 imes 16) imes (1024 imes 16^2)}\ &D_h in R^{(16 imes 16) imes 2560}\ &D_l in R^{(16 imes 16) imes 2560}\ &alpha in R^{2560 imes (1024 imes 16^2)}\ end{aligned} ight. end{equation}

　　另设$(4)$式$lambda=0.1$，使用RMSProp对$(4)$式进行优化，根据以上所列参数，Pytorch代码如下：

  1 #%%
  2 import torch,os,cv2
  3 import numpy as np
  4 import matplotlib.pyplot as plt
  5 from torch import optim,cuda 
  6 
  7 #读取图像并将LR放大、获得亮度单通道
  8 LR_path = r'E:DataSetsSRTestLR'
  9 HR_path = r'E:DataSetsSRTestHR'
 10 LR_imgs_YUV = np.zeros([1024,256,256,3],dtype=np.uint8)
 11 HR_imgs_YUV = np.zeros([1024,256,256,3],dtype=np.uint8)
 12 
 13 for i, j in zip(os.listdir(LR_path),range(1024)):
 14   img_path = os.path.join(LR_path,i)
 15   img = cv2.resize(plt.imread(img_path),(256,256),interpolation=cv2.INTER_CUBIC)
 16   LR_imgs_YUV[j] =  cv2.cvtColor(img, cv2.COLOR_RGB2YCrCb)
 17 for i, j in zip(os.listdir(HR_path),range(1024)):
 18   img_path = os.path.join(HR_path,i)
 19   HR_imgs_YUV[j] = cv2.cvtColor(plt.imread(img_path), cv2.COLOR_RGB2YCrCb)
 20 LR_imgs_YUV = torch.tensor(LR_imgs_YUV,dtype=torch.float32,device='cuda')
 21 HR_imgs_YUV = torch.tensor(HR_imgs_YUV,dtype=torch.float32,device='cuda') 
 22  
 23 #%%
 24 #减去区块平均值，获得预处理训练数据，定义各个变量 
 25 def imgs2patches(imgs, patch_size):
 26   #将图像集转换为区块集
 27   imgs_n = len(imgs)
 28   patch_n = int(imgs.shape[1]/patch_size)
 29   patches = torch.zeros([imgs_n*patch_n**2, patch_size*patch_size],device='cuda') 
 30   for i in range(patch_n): 
 31     for j in range(patch_n):
 32       t = imgs[:,i*patch_size:(i+1)*patch_size,j*patch_size:(j+1)*patch_size]
 33       t = torch.reshape(t,[imgs_n,-1]) 
 34       now = i * patch_n + j
 35       patches[imgs_n*now:imgs_n*(now+1),:] = t
 36   return patches.T
 37 
 38 LR_train = imgs2patches(LR_imgs_YUV[:,:,:,0], 16)#将训练集转换为矩阵
 39 HR_train = imgs2patches(HR_imgs_YUV[:,:,:,0], 16)    
 40  
 41 atom_n = 2560
 42 
 43 X = HR_train - HR_train.mean(0).reshape([1,-1]) #各个区块减去平均值
 44 Y = LR_train - LR_train.mean(0).reshape([1,-1])
 45 Dh = torch.normal(0,1,[16*16,atom_n],device='cuda') 
 46 Dl = torch.normal(0,1,[16*16,atom_n],device='cuda') 
 47  
 48 alpha = torch.normal(0,1,[atom_n,1024*16*16],device='cuda') 
 49 Dh.requires_grad_(True)
 50 Dl.requires_grad_(True)
 51 alpha.requires_grad_(True)
 52 opt = optim.RMSprop([Dh,Dl,alpha]) 
 53 #%%
 54 #训练模型
 55 from torch.utils.tensorboard import SummaryWriter 
 56 
 57 writer = SummaryWriter('TrainLogs/')
 58 def iter_one_epoch(lamb=0.01):
 59   L1 = torch.mean((X - torch.matmul(Dh,alpha))**2)
 60   L2 = torch.mean((Y - torch.matmul(Dl,alpha))**2)
 61   L3 = torch.mean(torch.abs(alpha))
 62   Loss = L1 + L2 + lamb * L3
 63   opt.zero_grad()
 64   Loss.backward()
 65   opt.step()
 66   return L1,L2,L3,Loss
 67  
 68 for i in range(1, 1500): 
 69   L1,L2,L3,Loss = iter_one_epoch(lamb=1)
 70   print(i,Loss.cpu().detach().numpy())  
 71   writer.add_scalar('L1', L1, int(i))
 72   writer.add_scalar('L2', L2, int(i))
 73   writer.add_scalar('L3', L3, int(i))
 74   writer.add_scalar('Loss', Loss, int(i))
 75   if i % 200 == 0 and i < 700: 
 76     for i in opt.param_groups: 
 77       i['lr'] = i['lr']*0.5  
 78 #%%
 79 print("保存字典")
 80 torch.save(Dl,'dictionaries/Dic_LR') #保存两个字典 
 81 torch.save(Dh,'dictionaries/Dic_HR') 
 82 #%%
 83 #用Dh重建HR图像验证训练结果
 84 def get_recon_LR_HR(n): 
 85   print(Dh.shape,Dl.shape)
 86  
 87   recon_LR = torch.matmul(Dl, alpha)
 88   recon_HR = torch.matmul(Dh, alpha) 
 89   LR = torch.zeros([1024,256,256],device='cuda')
 90   HR = torch.zeros([1024,256,256],device='cuda') 
 91   
 92   for i in range(n): 
 93     print(i)
 94     for j in range(16):
 95       for k in range(16): 
 96         LR[i,16*j:16*(j+1),16*k:16*(k+1)] = recon_LR[:,i+(j*16+k)*1024].reshape([16,16])
 97         HR[i,16*j:16*(j+1),16*k:16*(k+1)] = recon_HR[:,i+(j*16+k)*1024].reshape([16,16])
 98   return LR,HR 
 99 lr,hr = get_recon_LR_HR(100) 
100 
101 n = 10
102 fig = plt.figure(figsize=(15,15))
103 ax1,ax2,ax3,ax4 = fig.add_subplot(221),fig.add_subplot(222),fig.add_subplot(223),fig.add_subplot(224)
104 ax1.imshow(LR_imgs_YUV[n].cpu()[:,:,0]/255,cmap='gray')
105 ax2.imshow(lr[n].cpu().detach(),cmap='gray')
106 ax3.imshow(HR_imgs_YUV[n].cpu()[:,:,0]/255,cmap='gray')
107 ax4.imshow(hr[n].cpu().detach(),cmap='gray')
108 ax1.set_title('LR image',fontsize=20)
109 ax2.set_title('Reconstructed LR image',fontsize=20)
110 ax3.set_title('HR image',fontsize=20)
111 ax4.set_title('Reconstructed HR image',fontsize=20)
112 plt.show()

　　我另外对比过Adam和原始GD，迭代速度都没有RMSProp快，而且loss在稳定后是最小的。算法总共迭代了1500次，使用的是RMSProp默认的学习率，但每200次迭代都会下调为原来的一半。整个迭代在3090下用时10分钟。

　　以下是训练集的LR和它对应的HR图像的重建效果（只显示亮度通道），因为重建后没有加回均值，所以会有一块一块的伪影：

推理

顺序优化区块

　　论文的推理策略是按顺序重建SR图像的各个区块，同时约束相邻区块的重叠部分的相似性。定义LR图像放大前相邻区块之间的重叠为1像素宽，则放大后以及HR图像相邻区块之间有4像素宽的重叠，并且图像能被分成21x21个区块。则$(5)$式各个变量的规模如下

egin{equation} left{ egin{aligned} &alpha_{ij}in R^{2560 imes 1},;;i,j=1,2,...,21\ &p_{ij}in R^{(16 imes 16) imes 1},;;i,j=1,2,...,21\ &D_lin R^{(16 imes 16) imes2560}\ &D_hin R^{(16 imes 16) imes 2560}\ end{aligned} ight. end{equation}

　　优化代码如下：

 1 #%%
 2 import torch,cv2
 3 import matplotlib.pyplot as plt
 4 from torch import optim 
 5 from torch.nn import functional as F
 6 
 7 
 8 path_lr = r'E:DataSetsSRTestTestImagesLR001.jpg'
 9 path_hr = r'E:DataSetsSRTestTestImagesHR001.jpg'
10 img_lr = plt.imread(path_lr) 
11 img_hr = plt.imread(path_hr) 
12 
13 Dl = torch.load('dictionaries/Dic_LR').reshape([16,16,2560])
14 Dh = torch.load('dictionaries/Dic_HR').reshape([16,16,2560])   
15 def img_SR(LR_img, lambda1=0.5,lambda2=1,lambda3=1,lambda4=1,epoch=100):
16   LR_img = cv2.resize(LR_img,(256,256),cv2.INTER_CUBIC)
17 
18   LR_img_YUV = cv2.cvtColor(LR_img,cv2.COLOR_RGB2YCrCb)
19 
20   LR_img_Y = torch.tensor(LR_img_YUV[:,:,0],device='cuda',requires_grad=False)*1.  
21   SR_img_Y = torch.zeros([256,256],device='cuda',requires_grad=False)*1.  
22   
23   alpha_array = [] 
24   for i in range(21):
25     al = []
26     for j in range(21): 
27       al.append(torch.normal(0,1,[2560],device='cuda',requires_grad=True)) 
28     alpha_array.append(al)
29   
30   def SRcompat_loss(sr_patch, i, j, mean):
31     loss = 0
32     if i > 0:
33       loss += torch.mean((SR_img_Y[12*i:12*i+4,12*j:12*j+16] - (sr_patch[:4]+mean))**2)
34     if j > 0:
35       loss += torch.mean((SR_img_Y[12*i:12*i+16,12*j:12*j+4] - (sr_patch[:,:4]+mean))**2)
36     return loss
37 
38   #按顺序优化SR各个区块  
39   for i in range(21):
40     for j in range(21): 
41       alpha = alpha_array[i][j]
42       opt = optim.RMSprop([alpha])
43       mea = torch.mean(LR_img_Y[i*12:i*12+16,j*12:j*12+16])
44       lr_patch = LR_img_Y[i*12:i*12+16,j*12:j*12+16] - mea
45       for k in range(1, epoch): 
46         L1 = torch.mean(torch.abs(alpha))
47         L2 = torch.mean((torch.matmul(Dl,alpha) -  lr_patch)**2)
48         L3 = SRcompat_loss(torch.matmul(Dh,alpha), i, j,mea) 
49         L4 = torch.mean((torch.matmul(Dh,alpha) - lr_patch)**2)#额外的下采样一致性
50         Loss = lambda1 * L1 + lambda2 * L2 + lambda3 * L3 + lambda4 * L4
51         if k%500 ==0:
52           print(k,Loss)
53           for l in opt.param_groups:  
54             l['lr'] *= 0.5
55         opt.zero_grad()
56         Loss.backward()
57         opt.step()
58         if Loss < 1:
59           print(k,Loss)
60           break 
61       SR_img_Y[12*i:12*i+16,12*j:12*j+16] = torch.matmul(Dh,alpha).detach()+mea
62     plt.imshow(SR_img_Y.detach().cpu()/255,cmap='gray')
63     plt.show() 
64   LR_img_YUV[:,:,0] = SR_img_Y.cpu()
65   return LR_img_YUV
66 sr = img_SR(img_lr,1,1,1,0.5,5000) 
67  
68 img_sr = cv2.cvtColor(sr,cv2.COLOR_YCrCb2RGB)
69 plt.imshow(img_sr)
70 plt.show() 
71 #%%  
72 fig = plt.figure(figsize=(15,15))
73 ax1,ax2,ax3 = fig.add_subplot(131),fig.add_subplot(132),fig.add_subplot(133)
74 ax1.imshow(cv2.resize(img_lr,(256,256))),ax1.set_title('LR',fontsize=20)
75 ax2.imshow(img_sr),ax2.set_title('SR',fontsize=20)
76 ax3.imshow(img_hr),ax3.set_title('HR',fontsize=20)
77 plt.show()

　　迭代了好久，结果如下：

　　几乎和LR图像一模一样，原因应该就是没有用论文中提到的使用反向投影法进行精炼，以及使用$F$操作。看来只用稀疏表示是不太行的。

另一种推理方式

　　由于上述优化有先后顺序，后面的区块可能会损失精读来“迁就”前面已获得的SR区块，让重叠部分一致。因此论文在完成这一步后又加了一个全局一致性的约束，来调整已获得的SR图像，就是用所谓的反向投影法，但是写得很不明白。因此，我试验了一种直接对所有区块同时进行优化的方法。

　　首先定义图像相邻区块的重叠为8像素宽，也就是16x16区块的一半。则SR相邻区块的兼容性可以通过创建两个“补丁”图来约束，如下图所示：

　　即同时优化三个SR图像，第一张是最终的SR结果$x$，第二张用于约束$x$横向相邻区块之间的匹配度，第三张用于约束$x$纵向相邻区块之间的匹配度。也就是说，第一张图像相邻区块各取一半拼接成的图块要与第二、三图像中对应的区块一致。

　　综上，对于放大后的测试LR图像$y$，分别去掉左右、上下边缘的8像素宽、高的图块，获得用于匹配约束的$y_r,y_c$。然后分别定义相应的表示$alpha,alpha_r,alpha_c$。根据以上定义，各个参数的规模如下（LR图像区块被展开为向量的形式）：

egin{equation} left{ egin{aligned} &alphain R^{2560 imes (16 imes 16)}\ &alpha_rin R^{2560 imes (16 imes 15)}\ &alpha_cin R^{2560 imes (15 imes 16)}\ &yin R^{(16 imes 16) imes(16 imes 16)}\ &y_rin R^{(16 imes 16) imes(16 imes 15)}\ &y_cin R^{(16 imes 16) imes(15 imes 16)}\ &D_lin R^{(16 imes 16) imes2560}\ &D_hin R^{(16 imes 16) imes 2560}\ end{aligned} ight. end{equation}

　　则稀疏性约束、重建一致性约束、区块匹配度约束和最终的优化式如下：

egin{equation} egin{aligned} &L_{spars} = |alpha|_1+|alpha_r|_1+|alpha_c|_1 \ &L_{recon} = |D_lalpha - y|_1+|D_lalpha_r - y_r|_1+|D_lalpha_c - y_c|_1 \ &L_{comp} = |P_1D_halpha - D_halpha_r|_1+|P_2D_halpha - D_halpha_c|_1\ &minlimits_{alpha,alpha_r,alpha_c}Loss = lambda_1 L_{spars}+lambda_2L_{recon}+lambda_3L_{comp} end{aligned} end{equation}

　　$P_1,P_2$表示将SR图像映射到相应的匹配约束图像的操作，$lambda_1,lambda_2,lambda_3$用于平衡三个约束的占比。使用$L1$范数是因为它能生成噪声更少的图像。另外，为了防止梯度过大，代码中计算的各项范数会除以元素数量，公式中没有标明。代码如下：

  1 #%%
  2 import torch,cv2
  3 import matplotlib.pyplot as plt
  4 from torch import optim
  5 from torch.utils.tensorboard import SummaryWriter 
  6 import numpy as np
  7 
  8 
  9 path_lr = r'E:DataSetsSRTestTestImagesLR003.jpg'
 10 path_hr = r'E:DataSetsSRTestTestImagesHR003.jpg'
 11 img_lr = plt.imread(path_lr) 
 12 img_hr = plt.imread(path_hr) 
 13 
 14 Dl = torch.load('dictionaries/Dic_LR')
 15 Dh = torch.load('dictionaries/Dic_HR')
 16 def img_SR(LR_img, lambda1=0.5,lambda2=1,lambda3=1,epoch=100): 
 17   LR_img = cv2.resize(LR_img,(256,256),interpolation=cv2.INTER_CUBIC)
 18   LR_img_YUV = torch.tensor(cv2.cvtColor(LR_img,cv2.COLOR_RGB2YCrCb),dtype=torch.float32) 
 19   LR_img_Y = torch.tensor(LR_img_YUV[:,:,0],device='cuda').clone()   
 20   LR_img_Y_mean = torch.zeros_like(LR_img_Y,device='cuda')    
 21   for i in range(16):
 22     for j in range(16):
 23       LR_img_Y_mean[16*i:16*i+16, 16*j:16*j+16] += torch.mean(LR_img_Y[16*i:16*i+16, 16*j:16*j+16]) 
 24 
 25   LR_Y_sub_mean = LR_img_Y - LR_img_Y_mean
 26   LR_Y_sub_mean_r = LR_Y_sub_mean[:,8:-8]
 27   LR_Y_sub_mean_c = LR_Y_sub_mean[8:-8,:]  
 28  
 29   def img2patches(img):  
 30     patch_r = int(img.shape[0]/16)
 31     patch_c = int(img.shape[1]/16)
 32     patches = torch.zeros([16*16, patch_r*patch_c],device='cuda') 
 33     for i in range(patch_r): 
 34       for j in range(patch_c): 
 35         patches[:,i * patch_c + j] = torch.flatten(img[i*16:i*16+16,j*16:j*16+16])  
 36     return patches
 37   def patches2img(patches,row,col):
 38     img = torch.zeros([row*16, col*16],device='cuda')
 39     for i in range(row):
 40       for j in range(col):
 41         img[i*16:i*16+16, j*16:j*16+16] = patches[:,i*col+j].reshape([16,16])
 42     return img
 43    
 44   alpha = torch.normal(0,1,[2560,16*16],requires_grad=True,device='cuda') 
 45   alpha_r = torch.normal(0,1,[2560,16*15],requires_grad=True,device='cuda') 
 46   alpha_c = torch.normal(0,1,[2560,15*16],requires_grad=True,device='cuda') 
 47   y = img2patches(LR_Y_sub_mean)
 48   y_r = img2patches(LR_Y_sub_mean_r)
 49   y_c = img2patches(LR_Y_sub_mean_c) 
 50  
 51   opt = optim.RMSprop([alpha,alpha_r,alpha_c])
 52 
 53   writer = SummaryWriter('InferLogs2/')
 54   for i in range(1, epoch):
 55     l_alpha = torch.mean(torch.abs(alpha))
 56     l_alpha_r = torch.mean(torch.abs(alpha_r))
 57     l_alpha_c = torch.mean(torch.abs(alpha_c))
 58     L1 = l_alpha + l_alpha_r + l_alpha_c
 59 
 60     l_rec1 = torch.mean(torch.abs(torch.matmul(Dl,alpha)-y))
 61     l_rec2 = torch.mean(torch.abs(torch.matmul(Dl,alpha_r)-y_r))
 62     l_rec3 = torch.mean(torch.abs(torch.matmul(Dl,alpha_c)-y_c))
 63     L2 = l_rec1 + l_rec2 + l_rec3
 64 
 65     l_comp1 = torch.mean(torch.abs(patches2img(torch.matmul(Dh,alpha),16,16)[:,8:-8] - patches2img(torch.matmul(Dh,alpha_r),16,15)))
 66     l_comp2 = torch.mean(torch.abs(patches2img(torch.matmul(Dh,alpha),16,16)[8:-8,:] - patches2img(torch.matmul(Dh,alpha_c),15,16)))
 67     L3 = l_comp1 + l_comp2
 68 
 69     Loss = lambda1 * L1 + lambda2 * L2 + lambda3 * L3
 70 
 71     opt.zero_grad()
 72     Loss.backward()
 73     opt.step()
 74 
 75     writer.add_scalar('L1',L1,i)
 76     writer.add_scalar('L2',L2,i)
 77     writer.add_scalar('L3',L3,i)
 78     writer.add_scalar('Loss',Loss,i)
 79     if i % 50 == 0:
 80       print(i, Loss)
 81       img = patches2img(torch.matmul(Dh,alpha),16,16).detach()+LR_img_Y_mean
 82       plt.imshow(img.cpu()/255)
 83       plt.show()
 84       img = patches2img(torch.matmul(Dl,alpha),16,16).detach()+LR_img_Y_mean
 85       plt.imshow(img.cpu()/255)
 86       plt.show()
 87     if i % 500 == 0:
 88       for i in opt.param_groups:
 89         i['lr'] *= 0.5 
 90 
 91   LR_img_YUV[:,:,0] = patches2img(torch.matmul(Dh,alpha),16,16).detach()+LR_img_Y_mean 
 92   sr = cv2.cvtColor(np.array(LR_img_YUV.numpy(),dtype=np.uint8),cv2.COLOR_YCrCb2RGB) 
 93   return np.array(sr,dtype=np.uint8)
 94 
 95 SR_img = img_SR(img_lr, lambda1=100,lambda2=1,lambda3=1,epoch=1800) 
 96 
 97 fig = plt.figure(figsize=(15,15))
 98 ax1,ax2,ax3 = fig.add_subplot(131),fig.add_subplot(132),fig.add_subplot(133)
 99 ax1.imshow(img_lr),ax1.set_title('LR',fontsize=20)
100 ax2.imshow(SR_img/255),ax2.set_title('SR',fontsize=20)
101 ax3.imshow(img_hr),ax3.set_title('HR',fontsize=20)
102 plt.show()

　　不过实际效果也并不是很好：

总结

　　综上，单纯使用稀疏表示做SR效果并不如人意。因为论文中还加了其它的方式和技巧来减少重建图像的噪声，而我这个实验没有加入，并且，论文的实验是放大3倍，区块大小为3x3，我这里与其并不相同，所以没能重现出论文的效果。另外，还可能是优化算法的缘故，论文使用的是凸优化（可能有某种方式算出解析解，但我看到L1范数就放弃了），我则是梯度下降。

　　主要还是不想再做了。。。论文作者没有给出源代码，自己敲代码加写博客用了4天时间，想看其它论文了。

论文地址

　　Image Super-Resolution via Sparse Representation