經典超分辨率重建論文，基於稀疏表示。下面首先介紹稀疏表示，然后介紹論文的基本思想和算法優化過程，最后使用python進行實驗。
稀疏表示
稀疏表示是指，使用過完備字典中少量向量的線性組合來表示某個元素。過完備字典是一個列數大於行數的行滿秩矩陣，也就是說，它的列向量有無數種線性組合來表達列向量空間中的任意點。由於它的列數通常遠大於行數，可以使用占比很小的列向量來表示特定的向量，我們稱這種表示為稀疏表示。
那么如何獲得這個字典呢？它在特定的任務下有特定的取值。和煉丹類似，我們先要用大量數據來訓練這個矩陣，讓它提取出能稀疏表示這些數據的特征，進而擁有稀疏表示其它相似數據的能力。
訓練過完備字典的過程稱為稀疏編碼。設訓練數據集為矩陣$X=(x_1,x_2,...,x_n)\in R^{m\times n}$，待訓練矩陣為$A\in R^{m\times K}$，矩陣對每一數據的表示權重為$\alpha = (\alpha_1,\alpha_2,...,\alpha_n)\in R^{K\times n}$。進行如下優化：
\begin{align} \min\limits_{A,\alpha}\|\alpha\|_0\;\;\;s.t. \; \|A\alpha - X\|_2^2\le \epsilon \end{align}
容易理解，約束重建向量與原始向量差異的同時，優化表示權重$\alpha$的稀疏性。通過優化獲得所需的過完備字典后，我們就可以用它來稀疏表示新的數據了。對於某個向量$y$，我們可以進行類似的優化來獲得它在字典中的稀疏表示：
\begin{align} \min\limits_{\alpha}\|\alpha\|_0\;\;\;s.t. \; \|A\alpha - y\|_2^2\le \epsilon \end{align}
因為零范數的離散性，以上優化是NP-難問題，無法用正常的優化算法解決。所以通常會對上式進行調整，以方便優化。而在線性約束下，對一范數的優化有稀疏性（點擊鏈接理解），因此可以轉換為對一范數的優化。然后根據拉格朗日乘子法（據說如此，但我覺得這個轉換並不等價），不等式約束可以移入優化中：
\begin{align} \min\limits_{\alpha} \lambda \|\alpha\|_1 + \|A\alpha - 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\{ \begin{aligned} &D_l\alpha_i \approx y_i \\ &D_h\alpha_i \approx x_i \end{aligned} \right.$
在測試階段，我們已經有了訓練好的$D_l$和相對應的$D_h$。對於測試圖像$y_t$，首先通過優化獲得$y_t$在$D_l$中的稀疏表示$\alpha_t$，此時有$D_l\alpha_t \approx y_t$。然后用這個表示通過$D_h$映射出對應的SR圖像，即$\hat{x}_t=D_h\alpha_t$。
訓練過程
訓練過程就是訓練上述的過完備字典對。因為性能的因素，我們不可能直接對整張圖進行稀疏編碼，論文是將圖像分為方形的區塊（patch）進行編碼的。因此，用於訓練的成對數據不是整張的LR-HR圖像對，而是所有圖像分割開來的區塊對。現在把LR訓練集的所有區塊表示為$Y=(y_1,y_2,...,y_n)\in R^{M\times n}$，相應的HR訓練集區塊表示為$X = (x_1,x_2,...,x_n)\in R^{N\times n}$。如果放大倍數為$t$倍，則有$N=t^2M$。由於論文中是先將LR圖像通過Bicubic放大到與HR圖像一致，則$N=M$。
優化式很直觀：
\begin{align} \min\limits_{D_l,D_h,\alpha}\frac{1}{N}\|X - D_h\alpha\|_2^2+\frac{1}{M}\|Y - D_l\alpha\|_2^2 + \lambda \|\alpha\|_1 \end{align}
其中$D_l\in R^{M\times K},D_h\in R^{N\times K}$，分別表示待訓練的LR和HR字典，$K$表示字典的原子數；$\alpha\in R^{K\times 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$個區塊，優化式如下：
\begin{align} \min\limits_{\alpha_i} \lambda\|\alpha_i\|_1 + \|FD_l\alpha_i - Fp_i\|_2^2+ \|PD_h\alpha_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\times 16\times 1=256$。另外，對於每個圖像區塊中像素的亮度值，都會減去區塊的亮度平均值，使得模型能專注於紋理的訓練與生成，這一點也是與論文一致的。
訓練
綜上，對於$(4)$式中的各個參數，有：
\begin{equation} \left\{ \begin{aligned} &X \in R^{(16\times 16)\times (1024\times 16^2)}\\ &Y \in R^{(16\times 16)\times (1024\times 16^2)}\\ &D_h \in R^{(16\times 16)\times 2560}\\ &D_l \in R^{(16\times 16)\times 2560}\\ &\alpha \in R^{2560\times (1024\times 16^2)}\\ \end{aligned} \right. \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:\DataSets\SRTest\LR' 9 HR_path = r'E:\DataSets\SRTest\HR' 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)$式各個變量的規模如下
\begin{equation} \left\{ \begin{aligned} &\alpha_{ij}\in R^{2560\times 1},\;\;i,j=1,2,...,21\\ &p_{ij}\in R^{(16\times 16) \times 1},\;\;i,j=1,2,...,21\\ &D_l\in R^{(16\times 16) \times2560}\\ &D_h\in R^{(16\times 16) \times 2560}\\ \end{aligned} \right. \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:\DataSets\SRTest\TestImages\LR\0001.jpg' 9 path_hr = r'E:\DataSets\SRTest\TestImages\HR\0001.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圖像區塊被展開為向量的形式）：
\begin{equation} \left\{ \begin{aligned} &\alpha\in R^{2560\times (16\times 16)}\\ &\alpha_r\in R^{2560\times (16\times 15)}\\ &\alpha_c\in R^{2560\times (15\times 16)}\\ &y\in R^{(16\times 16) \times(16\times 16)}\\ &y_r\in R^{(16\times 16) \times(16\times 15)}\\ &y_c\in R^{(16\times 16) \times(15\times 16)}\\ &D_l\in R^{(16\times 16) \times2560}\\ &D_h\in R^{(16\times 16) \times 2560}\\ \end{aligned} \right. \end{equation}
則稀疏性約束、重建一致性約束、區塊匹配度約束和最終的優化式如下：
\begin{equation} \begin{aligned} &L_{spars} = \|\alpha\|_1+\|\alpha_r\|_1+\|\alpha_c\|_1 \\ &L_{recon} = \|D_l\alpha - y\|_1+\|D_l\alpha_r - y_r\|_1+\|D_l\alpha_c - y_c\|_1 \\ &L_{comp} = \|P_1D_h\alpha - D_h\alpha_r\|_1+\|P_2D_h\alpha - D_h\alpha_c\|_1\\ &\min\limits_{\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:\DataSets\SRTest\TestImages\LR\0003.jpg' 10 path_hr = r'E:\DataSets\SRTest\TestImages\HR\0003.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天時間，想看其它論文了。