Dictionary learning based noisy image super-resolution via distance penalty weight model

In this study, we address the problem of noisy image super-resolution. Noisy low resolution (LR) image is always obtained in applications, while most of the existing algorithms assume that the LR image is noise-free. As to this situation, we present an algorithm for noisy image super-resolution which can achieve simultaneously image super-resolution and denoising. And in the training stage of our method, LR example images are noise-free. For different input LR images, even if the noise variance varies, the dictionary pair does not need to be retrained. For the input LR image patch, the corresponding high resolution (HR) image patch is reconstructed through weighted average of similar HR example patches. To reduce computational cost, we use the atoms of learned sparse dictionary as the examples instead of original example patches. We proposed a distance penalty model for calculating the weight, which can complete a second selection on similar atoms at the same time. Moreover, LR example patches removed mean pixel value are also used to learn dictionary rather than just their gradient features. Based on this, we can reconstruct initial estimated HR image and denoised LR image. Combined with iterative back projection, the two reconstructed images are applied to obtain final estimated HR image. We validate our algorithm on natural images and compared with the previously reported algorithms. Experimental results show that our proposed method performs better noise robustness.


Introduction
Single image super-resolution (SR) is a classical problem in computer vision. In general, it uses signal processing techniques to recover a high resolution (HR) image from only one low resolution (LR) image. SR methods can be broadly classified into three categories: interpolationbased methods, reconstruction-based methods, and example-based methods.
Interpolation-based SR such as [1,2] has been proposed for in various applications and it demonstrates the advantage of fast computational simplicity. But they usually fail to generate fine details in discontinuous regions and often result in introducing blurring of edges and other high-frequency features in practice [3].
Reconstruction-based methods usually integrate one or more sophisticated priors such as gradient profile prior [4], edge prior [5], and total variation [6] into SR literature to estimate the missed details. Recently, sparse-based regularization [7][8][9][10] has also been shown to be a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 particularly effective for the ill-posed problems of SR. Usually, these methods achieved impressive results in preserving sharper edges and suppressing aliasing artifacts. However, the performance depends heavily upon a rational prior imposed on the up-sampled image [11].
Over the years, many example-based SR methods [12][13][14] have been proposed with demonstrated promising results and become the mainstream approaches of SR domain. The methods assume that the missing high frequency details can be estimated based on learning the mapping relationship from LR-HR patch pairs of external database and input LR patches. Two kinds of relationship models exist for these methods. One is that between LR patches and the corresponding HR patches in the database. After Freeman et al. [15] used Markov network to model the relationship, regression functions [16] are employed to exploit the relationship between HR and LR patch pairs. In addition, supervised or semi-supervised learning models are introduced into some of the algorithms [17][18][19]. Recently, a mapping of LR-HR image pairs was learned using a deep convolutional neural network [20], and has shown favorable results. D. Dai et al. [21] jointly learned a collection of regressors from LR to HR patches, which collectively yielded the smallest error for all training data. The other is that between LR example patches and input LR patches. Most of the methods [22,23] is based on Nearest Neighbor Embedding (NNE). In these methods, a fixed number of nearest neighbors are extracted from database for each input LR patch, and then the corresponding HR patches are used to estimate the output HR patch by a linear combination determined by LR patch and its neighbors. Despite the algorithms are demonstrated by successful results, they highly depend on the number of neighbors which is difficult to determine. For this problem, [24] operates on a dynamic k-nearest neighbor algorithm, where k is small for test point with highly relevant neighbors and large others. Some researchers calculate the distance between input patch and its neighbors respectively. The neighbors will be abandoned when the distance is smaller than mean value. Yang [25] exploited sparse coding to perform image SR. The algorithm assumes that LR-HR patch pairs share the same sparse coefficients with respect to their respective dictionaries which are jointly learned from a set of external training images. It can be considered as neighbor embedding in sparse domain without choosing the number of neighbors. Since then, sparse coding is applied to SR problem [21][22][23], and achieves impressive results. Zeyde [26] used dimensionality reduction and orthogonal matching pursuit for sparse representation to improve efficiency. S. Wang [27], proposed a semi-coupled dictionary learning model, under which a pair of dictionaries and a mapping function describing the relationship between sparse coefficients of LR-HR patch pairs will be simultaneously learned. In [28], kernel ridge regression is employed to connect sparse coefficients of LR-HR patch pairs. Kaibing Zhang [29] determine the relationship between LR image patches and HR image patches by assuming that LR image patches and HR image patches are share the same sparse coefficients. R. Timofte et al. [30] proposed a fast image SR method called anchored neighbourhood regression (ANR) which learns sparse dictionaries and regressors anchored to dictionary atoms. This algorithm is faster, while making no compromise on quality. R. Timofte et al. [31] then produced an improved variant of ANR. The study in [31] enhanced these features and anchored regressors for ANR. Instead of learning the regressors on the dictionary, their method uses the full training material. It obtained improved quality, and became the fastest method indisputably. S. Gu [32] proposed a convolutional sparse coding based SR method to address consistency issue. In addition, researches show that image structures tend to repeat themselves within and across scales. [33][34][35] exploits the self-similarity of structures in nature image and extracts the database directly from the LR input image instead of the external database. Good reconstruction quality relies on much additional memory and running time to build counterparts across different scales in a recursive scheme. Therefore, its application is limited.
Although the algorithms can results in better performance, most of the SR algorithms including other learning-based methods assume that the input LR image is noise-free. Such assumption is not in accord with real applications. The algorithms are less robustness to noisy image SR. So another challenge is the super-resolution for noisy images. While compared with SR on clear LR input images, less attention has been paid to develop effective SR algorithms for noisy ones. J. Xie [36] first employs an adaptively regularized Shock filter to tackle the jagged noise, and then perform SR for depth image. The disadvantage of such scheme is that the artifacts can be created in denoising process and magnified in super-resolution process. Therefore, researchers started on simultaneously denoising and super-resolution. In [37], LR training images are magnified by a TV regularization model with a constraint before dictionaries training stage. However, the level of noise dealt with the method is small. Furthermore, it focuses on magnification only. Based on the current research status, we devote to design an algorithm to complete SR and denoising in the same framework to deal with noisy image patches.
Sparse representation makes the signal energy only concentrated in a few atoms. Because of the special nature, some sparse coding based SR algorithms such as [25] show certain robustness to noisy image. In addition, sparse representation has been successfully employed in image denoising [38,39], image restoration [40,41] and other processing [42,43]. The dictionary plays an important role in the sparse representation process. A predefined analytical dictionary (e.g., wavelet dictionary, Gabor dictionary) make the coding fast and explicit, but it is less effective to model the complex local structures of natural images. A synthesis dictionary (e.g., K-SVD dictionary) can be learned from example natural images and has more expensive computation but can better model complex image local structures [44]. In recent years, lots of dictionary learning methods have been proposed and achieved obvious performance. Feng et al. [45] propose to learn jointly the projection matrix for dimensionality reduction and the discriminative dictionary for face representation. Zhang et al. [46] propose a semisupervised label consistent dictionary learning framework for machine fault classification. Inspired by these, we introduce sparse theory to our research. The synthesis procedure is illustrated in Fig 1. The input LR image and example images are firstly cropped into patches. The example images are noise-free. Then the features of example patch pairs are extracted, which will be learned for dictionary pair. For each input LR patch, according to its features, it is easy to achieve simultaneously similar dictionary atom pairs ðu h i ; u l i Þ finding and calculating distance b i between input LR patch and its similar atoms. Next, combined with the input LR image patch feature, LR dictionary atom u l i and distance b i are used to compute weight ω i . After the weight is computed, we can obtain estimated HR image patch and denoised LR image patch from u h i o i . Put all the estimated HR patches into an estimated HR image, which is computed by averaging in overlapping regions. In the same way, we obtain the denoised LR image from all the denoised LR patches. At last, combined with the iterative back projection (IBP), the estimated HR image and the denoised LR image are applied to obtain the final output HR image.
The contributions can be summarized as follows.
(1) Different from the conventional methods, the proposed algorithm can process noisy image, and present for simultaneously image superresolution and denoising. Furthermore, in the training stage of our method, LR example images are noise-free. For different input LR images, even if the noise variance varies, the dictionary pair does not need to be retrained.
(2) The core idea of our proposed method is that the estimated HR patch is weighted average of similar HR example patches. To reduce computational cost for finding similar patches from millions of examples, example patches are replaced by the learned sparse dictionary which makes the signal energy only concentrate in few atoms.
(3) Penalty function is applied to least squares regression regularized by l 2 -norm for modeling weight. It makes the objective function treat each similar atom unequally. The function is determined by the similarity between input LR patch and its similar atom of LR dictionary. When the similarity is strong, we make the penalty small, which forces large weight at the same time. Conversely, when the similarity is weak, we make the penalty large, which forces small or zero weight at the same time.
(4) LR example patches subtracted mean pixel value are used for training dictionary rather than just their gradient features like other literatures such as [25]. In the training stage, for each LR example patch, we first subtract its mean pixel value, then connect it to its corresponding HR example patch into a single vector. All the new vectors are used as new HR examples to learn HR dictionary. Thus, the HR dictionary represents textures of HR example patches, but also that of LR example patches which are noise-free. Therefore, in the reconstruction stage, the HR dictionary can also be used to recover denoised input LR patches. This is different from conventional learning methods. Combined with iterative back projection (IBP), the denoised LR patches are applied to enhance robustness to noise.
The remainder of this paper is organized as follows. The proposed algorithm is presented in detail in Section 2. Experimental results and comparisons are demonstrated in section 3. Section 4 concludes this paper.

The proposed method
Firstly, let us recall the image degradation model which is shown in Eq (1). Given an observed LR image Y 2 R M that is a degraded version of a HR image X 2 R N of the same scene Where, G s is the down-sampling operator with scaling factor s; H is the blurring operator; v is the noise. It is the task of SR reconstruction to recover X from Y as accurate as possible. It is considered that the image is noise-free by conventional SR methods.

Example database
From the example images fI h 1 ; I h 2 ; :::; I h N g, LR images fI l 1 ; I l 2 ; :::; I l N g are first obtained, which are considered as noise-free ones. For each image I h j , its corresponding LR image I l j is determined

Distance penalty weight model
For the super-resolution, given a LR image Y L , which is generated form HR image X H by Eq (1), the task is to recover the unknown X H from Y L with the help of example patch pairs. The algorithm is performed with patch for the unit. Similar to [25], Y L is firstly divided into overlapping patches Where, y l i is the vectorized LR image patch of size The estimated vectorized HR image X H can be represented as Where,x h i is the estimated HR image patch of size ffiffiffi ffi w p Â ffiffiffi ffi w p . According to Eq (1), the relationship can be described by Where,v i is the noise. We assume that it is Gaussian noise with zero-mean and variance σ 2 . Thus, it become the purpose of super-resolution to estimate HR image patch x h i from input LR image patch y l i . As we known, for each x h i , it can be approximated by HR example patches through weighted average, which have similar structures. Therefore, based on this core idea, the problem in this method is to find the similar patches of x h i in database and to calculate the weight. Due to the repetition of local structures of images, a subset of patches Where, weight vector is ω i = [ω i1 , ω i2 , . . ., ω ij , . . ., ω ik ] T , k is the number of the patch pairs in this subset ðu h i ; u l i Þ. There are many methods to determine the weight, such as set the weights to be inversely proportional to the distance between patches. These methods relying on number of similar patches heavily, and cannot suppress noise. Now, we discuss a new weight model in details.
According to the degradation model Eqs (1) and (7), we have From Eq (8), we can obtain Where, v i is assumed as Gaussian noise with zero-mean and variance σ 2 . Thus, Where, ε i is related to σ 2 . We can see that the LR patch y l i can be represented by the same weight vector ω i over u l i , with an error ε i . That is to say, we can get the weight from input LR image patch and similar LR example patches with a controlled error.
Based on the above discussions, We formulate the weight solution as a least squares regression regularized by l2-norm: From Eq (12), the objective function treats the patches u l i equally. It is not flexible to obtain accurate weights for the input patch. Motivated by this, we introduce distance penalty to the least square problem Where, Á denotes a point wise vector product, between y l i and each similar example patch in u l i . When the similarity between u l ij and y l i is strong, we make the b ij small, which forces large ω ij at the same time. Conversely, when the similarity is weak, we make b ij large, which forces small or zero ω ij at the same time. It is simply determined by the squared Euclidean distance.
Eq (13) can be written as Where, λ is a regularization parameter. According to Eq (10), we have Where, γ is a positive constant. So we set λ = γσ 2 , when σ 6 ¼ 0. Thus, the main task in reconstruction stage is to find the patches u l i from p l , which is similar to y l i and compute the weight. Squared Euclidean distance can be adopted in to quantify the similarity. The corresponding u h i is assumed to have similar structures with x h i . But it is uneasy to find similar patches for each input patch from millions of example patch pairs. It will take lots of time for the repetitive computation. Sparse dictionary make the signal energy only concentrate in few atoms, and some sparse coding based SR algorithm [25] show certain robustness to noisy image, so that we use a learned sparse dictionary instead of examples. We find similar patch pairs ðu h i ; u l i Þ from dictionary atom pairs, meaning ðu h i ; u l i Þ 2 ðD h ; D l Þ. Two dictionaries D h and D l are trained to have the same sparse coding for each HR and LR patch pair. Similar to Yang [25] and Chang [22], we subtract the mean pixel value for each HR example patch, so that the dictionary D h represents image textures rather than absolute intensities. In the reconstruction stage, the mean value for each estimated patch is then predicted by its LR version. Also we employ first-and second-order derivatives as the feature extraction for LR example patches to train. Thus, D l represents the gradient feature of images rather than absolute intensities. The four filters used here are: In addition, to enhance robustness to noise, we also subtract mean pixel value for each LR example patch, and connect the LR example patch to its corresponding HR example patch into a single vector, which is also used to learn D h . Thus, dictionary D h represents textures of HR example patches, but also that of LR example patches which are noise-free. In the reconstruction stage, the D h can also be used to recover denoised input LR patches. This is different from conventional learning methods.
From above, the training set is obtained by Where, (p h ,p l ) is original HR-LR patch pairs in Eq (3), " p h i is the mean value of p h i , " p l i is the mean value of p l i , F(Á) is the operator to get four gradient vectors by Eq (16) and connect the four vectors into a single vector.
The set (P H ,P L ) is used to jointly train the dictionaries as Where, N and M are the vector dimensions of P H and P L , respectively. To solve the problem easily, Eq (18) can be rewritten as The minimization of Eq (19) is a typical patch-based sparse problem. Many methods can be used to solve it. Yang [25] proposed the framework and acquired good results. However, it takes a large amount of time to solve this sparse model. Zeyde [26] improve the execution speed by dimensionality reduction on the patches through PCA and Orthogonal Matching Pursuit for the Sparse coding. For sparse dictionaries learning, we use the approach of Zeyde [26].
Gradient features(see Eq (16)) of LR example patches are used to learn LR dictionary. D l represents the image gradient feature and u l i 2 D l . Therefore, the weight model is rewritten bŷ s:t: Where,ω i is the weight. This problem Eq (20) is l 2 -norm constraint. We solve it forω i by taking @L Where, The final optimal weight is obtained by rescaling it so that P k j¼1ôij ¼ 1.

Reconstruction
Based on the above discussions, for each input y l i , we start by extracting its gradient features and finding k similar atom pairs ðu h i ; u l i Þ. Because the dictionary atoms are learned basis vectors, we find the similar atoms based on the correlation between the LR dictionary atoms and input LR patch rather than the Euclidean distance. Now, we describe how to compute the correlation.
Fðy l i Þ can be represented by dictionary D l ¼ ½d l 1 ; d l 2 ; :::; d l j ; :::; d l nd (d l j is the LR dictionary atom, nd is dictionary size) Where, β = [β 1 , β 2 , . . ., β j , . . ., β nd ], β j is the correlation between d l j and Fðy l i Þ. Eq (23) shows that every dictionary atom makes its own contribution to representing the input patch. The contribution of the j th atom d l j can be evaluated by β j . In other words, β j is a measurement of the similarity between the input patch and the j th dictionary atom. We consider that the larger the β j , the larger scale of similarity between input patch Fðy l i Þ and dictionary atom d l j ; and a small β j means that there is little similarity. We can solve β by Thus, ðD l Þ T Fðy l i Þ could return the correlation. In Eq (20), we use distance b i as the penalty. When the similarity between Fðy l i Þ and d l j is strong, we make the b ij small, which forces largê o ij at the same time. Conversely, when the similarity is weak, we make b ij large, which forces small or zeroô ij at the same time. Therefore, we use the reciprocal of β j to compute the penalty. The atom pairs corresponding to the maximal k correlation coefficients constitute Where, Sort(a, num) is a function returning num top biggest values of vector a, abs(.) is absolute value operation. The scheme can achieve simultaneously similar atoms finding and distance computing. If σ = 0, after finding similar atoms, we set b i = 1.
After this, we can easily obtain the weightω i by Eq (20) and u h iωi . According to section 2.2, the reconstructed vector u h iωi represents the estimated HR patch and the denoised LR patch correspondent to y l i . And the estimated patch and the denoised patch are subtracted mean pixel value. Based on this, we havê Where,x h i is the estimation of x h i ,ŷ l i is the denoised patch of y l i , C 1 2 R w 1 is an all-one column vector, C 2 2 R w 2 is an all-one column vector, w 1 is the size ofx h i , w 2 is the size ofŷ l i , E(Á) is the mean evaluation operator.
Noise here is assumed as zero-means, so We can see that the noise has little effect on image mean. The mean ofŷ l i andx h i could be estimated by the mean of y l i . Eq (26) can be written bŷ Put all estimated patchesx h i into a HR imageX H , which is computed by averaging in overlapping regions. In the same way, we obtain a denoised imageŶ L fromŷ l i . In order to strengthen the reconstruction constraint Eq (1), we compute the final estimated HR image X Ã by The iterative back-projection (IBP) method [32] is used to solve this optimization problem Where, X Ã t is the estimate of the HR image at the t th iteration, " s denote up-scaling by factor s, p is a symmetric Gaussian filter.
The entire SR process is summarized as Algorithm 1.

Experiments
In this section, we will show the robustness of the proposed algorithm to noise and compare the state-of-the-art methods [20,22,25,26,31,32]. In the training stage, we used 77 standard natural images as training set. For testing, we used Set5 [20,31], Set14 [20,31] and B100 [20,31] to evaluate the performance of upscaling factors ×2, ×3 and ×4, respectively. Set5 and Set14 contain 5 and respectively 14 images for super-resolution evaluation. B100 contains 100 testing images of Berkeley Segmentation Dataset called BSDS300.
All LR images (training or test images) are generated from the original HR images. Firstly, the original HR images are directly blurred and down-sampled. The MATLAB function "imresize" is used here to complete the process. The function "imresize" involved a smooth filtering before down-sampling. Similar to [7], the noise is generated by MATLAB function "randn", and σ times noise is added to the blurred and down-sampled test images. It should be noted that LR example images for training dictionary are noise-free. For color images used in experiments, SR algorithms are performed only on luminance channel, because humans are more sensitive to illuminant changes. Therefore, we first changes channels into YCbCr ones and then apply our method to the Y channel. We interpolate the color layers (Cb, Cr) using bicubic interpolation.

Parameters
In this section, we analyze the main parameters of our algorithm. The standard settings we use are Set5 [20,31] database, dictionary size 1024, γ = 0.08 and k = 24 for upscaling factor ×2, k = 8 for upscaling factor ×3, ×4. Peak signal-to-noise ratio (PSNR) and reconstruction time were used as the objective criteria.
3.1.1 Regularization parameter. γ is a key regularization parameter of our method. Here, we validate the effectiveness of using different γ, and choose an appropriate one. The results of Set5 are shown in Fig 2. Experimental setting is dictionary size 1024 and k = 24 for upscaling factor ×2, k = 8 for upscaling factor ×3, ×4. We can see that the curves are not monotonic, and PSNR peaks at γ = 0.08. For different datasets, the optimal γ is slightly different (0.06 of Set14 and B100 compared to 0.08 of Set5) for reconstruction quality. The results of Set14 and B100 are shown in S1-S6 Figs. Therefore, we suggest determining γ to be around 0.08 in practice.
Here, in all of our following experiments, we set γ as 0.08 for convenience.

Dictionary size.
In this experiments, dictionary size is varied from 32 up to 2048, while the training samples are extracted from the same training images previously mentioned. In Fig 3, we present the results that show the relation between our method's performance and the dictionary size when γ = 0.08 and k = 24 for upscaling factor ×2, k = 8 for upscaling factor ×3, ×4. Actually, noise has little effect on reconstruction time. So we only show the reconstruction time when σ = 10. We can see that the larger we learn the dictionary, the better reconstruction quality becomes. However, this comes with a higher computational cost. The result is the same as that of [25,47]. Other datasets Set14 and B100 can also achieve similar results. The results of Set14 and B100 are shown in S7-S12 Figs. In practice, we suggest choosing the appropriate dictionary size as a tradeoff between reconstruction quality and computation. Dictionary size here is 1024 in our following experiments.    show the reconstruction time when σ = 10. We can see that k = 24 is best for reconstruction quality when upscaling factor is ×2. The PSNR peaks at k = 8 when upscaling factor is ×3 or ×4. Moreover, average reconstruction time increases distinctly as k increases. It is due to the fact that by having a larger k, the computation of matrix inversion in Eq (21) increases. Other datasets Set14 and B100 can also achieve similar results. The results of Set14 and B100 are shown in S13-S18 Figs. Therefore, in resource-limited systems, a reasonable selection of k depends on the tradeoff between reconstruction quality and computational time. We will use k = 24 when upscaling factor is ×2, k = 8 when upscaling factor is ×3 or ×4 in our further experiments.
3.1.4 Patch size and overlap. Intuitively, using a too large or too small patch size tends to produce a smooth or unwanted artifact as noticed also in [25,29] and a larger overlapping leads to a better SR results [25]. Therefore, patch size is set as 6×6, 6×6 and 8×8 for upscaling factor ×2, ×3 and ×4, respectively, and overlap is set as 4, 3 and 4 for upscaling factor ×2, ×3 and ×4, respectively.

Performance evaluation
In this section we analyze the performance of our algorithm in quantitative and qualitative comparison with the state-of-the-art methods including NE [22], SCSR [25], Zeyde [26], A+ [31], SRCNN [20], and CSC [32]. We also show the reconstruction times of the algorithms. The code of the compared method was downloaded from the authors' homepage. Peak signalto-noise ratio (PSNR) and structural similarity (SSIM) were used as the objective criteria. The parameters are analyzed in the previous section. Besides the patch size and overlap(see section 3.1.4), the other parameter are unified (γ = 0.08, dictionary size = 1024, k = 24 for upscaling factor ×2, k = 8 for upscaling factor ×3 and ×4). Tables 1-3 list the PSNR and SSIM comparisons. When σ = 0, the approach CSC [32] achieves the best performance. But it is not in accord with real application. When σ 6 ¼ 0, as repeatedly shown, the results demonstrate the superiority of our proposed algorithm over other approaches on Set5, Set14 and B100. The average PSNR of the recent method CSC [32] is 0.24 dB (Set14, upscaling factor ×4, σ = 5) and 7.4 dB (Set5, upscaling factor ×2, σ = 20) behind our method. Compared with CSC, for dataset B100, the average PSNR improvement is

Quality.
Noisy image super-resolution from the minimum 0.52 dB (upscaling factor ×4, σ = 5) to the maximum 6.18 dB (upscaling factor ×2, σ = 20). In addition, our method improves on average 3.62 dB (Set5, upscaling factor ×2, σ = 20) over the next top robustness method SCSR [25]. Figs 5-8 provide a visual assessment. We can see that our method gets similar quality performance as the top methods it was compared to when σ = 0, and it has the strongest robustness.

Reconstruction time.
Average reconstruction time of test images in Set5 was compared when σ = 10. Actually, noise has little effect on test results. The experiments were conducted on the same computer. The results are summarized in Table 4. The reconstruction

Effect of IBP
Combined with iterative back projection (IBP), the denoised LR patches are applied to improve SR performance in our algorithm. According to [47], IBP has an important role to improve SR performance. But if the input is a noisy image, the model of IBP will propagate the noise to the HR image. Experimental results show that if we use IBP algorithm directly on the Noisy image super-resolution input LR image, the performance will become worse. The results are listed in Table 5. The iteration number of IBP here is 20. From this comparison, we can see that the superiority of our method is obvious. Other datasets Set14 and B100 can also achieve similar results. The results of Set14 and B100 are shown in S1 Table.

Effect of distance penalty
Distance penalty is applied to model the weight. To check the effect of the penalty for improving SR performance, we perform our method with and without the penalty respectively on We can see that our method with distance penalty obtains better performance and the superiority of our method with distance penalty is obvious. Other datasets Set14 and B100 can also achieve similar results. The results of Set14 and B100 are shown in S19-S24 Figs.

Conclusion
In this research, we proposed an algorithm of noisy image super-resolution based on sparse representation. For the problem of noisy image super-resolution, most of the existing methods will become less effective because they assume that the input LR image is noise-free. The Noisy image super-resolution proposed algorithm can achieve simultaneously image super-resolution and denoising. For different input LR images, even if the noise variance varies, the dictionary pair does not need to be retained. The core idea of the proposed algorithm is that HR image patch is reconstructed through weighted average of similar HR example patches. In particular, atoms of learned sparse dictionary are used to compute the weight and reconstruct HR patch instead of example patches. This strategy can reduce time computation and suppress noise. In addition, LR example patches subtracted mean pixel value are also used to learn dictionary rather than just their gradient features, which will help IBP to further improve the SR performance. The experimental results show that our method performs better noise robustness.  Noisy image super-resolution