# Maximizing Nonlocal Self-Similarity Prior for Single Image Super-Resolution

Authors: Jianhong Li; Kanoksak Wattanachote; Yarong Wu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Prior knowledge plays an important role in the process of image super-resolution reconstruction, which can constrain the solution space efficiently. In this paper, we utilized the fact that clear image exhibits stronger self-similarity property than other degradated version to present a new prior called maximizing nonlocal self-similarity for single image super-resolution. For describing the prior with mathematical language, a joint Gaussian mixture model was trained with LR and HR patch pairs extracted from the input LR image and its lower scale, and the prior can be described as a specific Gaussian distribution by derivation. In our algorithm, a large scale of sophisticated training and time-consuming nearest neighbor searching is not necessary, and the cost function of this algorithm shows closed form solution. The experiments conducted on BSD500 and other popular images demonstrate that the proposed method outperforms traditional methods and is competitive with the current state-of-the-art algorithms in terms of both quantitative metrics and visual quality.


---

## Body

## 1. Introduction

The technology of single image super-resolution (SISR) has been widely used in many fields, such as medical imaging, remote sensing and digital home etc., which refers to estimating the corresponding high-resolution (HR) image just according to an input low resolution (LR) one through software computing way. Due to the one-to-many mapping relationship between the LR and HR images captured from the same scene, the process from LR image to HR one is obviously a typically ill-posed problem. For obtaining HR image with sharp edges and fine details, researchers generally introduce some priors explicitly or implicitly to impose on the inverse imaging process; thus the quality of the reconstructed HR image is closely related to the introduced priors. The common priors for SISR include edge prior [1], gradient profile prior [2], sparse prior [3, 4], etc.Except a handful of SISR algorithms designed for particular types of images [5–7], natural images are the main target for most SISR algorithms. Generally speaking, there are mainly three categories of SISR approaches including interpolation-based methods, reconstruction-based methods, and learning-based methods.Interpolation-based methods assumed that the structures of one natural image are piece-wise smoothed [8]. According to this assumption, these methods take a kernel interpolator or base function to approximate the pixel values in HR grid and then transform the discrete digital image into continuous surface for estimating the unknown pixel values at specific position. This category is simple, efficient, even real-time processing because they do not discriminate edge or texture from smoothed region. The drawback of this category is that they cannot restore the high frequency information which was lost during imaging procedure and will lead to visible blur and noise along the edge in the reconstructed results to degrade its perception.The reconstruction-based methods usually incorporate priors explicitly [2, 9, 10], which require not only making the reconstructed HR image to satisfy the prior but also promoting the image degenerated from the reconstructed HR image by imaging model as similar as possible to the input one. The priors used in reconstruction-based methods are obtaining from edge statistical information of the desired HR images [2, 10]; hence this category is good at reconstructing edges and can suppress noises and artifacts near edges effectively. However, when the super-resolved magnification factor is larger, the prior term provides less and less effective information, which leads to be constrained for the ability to restore high frequency information, and exhibits unpleasing artifacts in the reconstructed texture and detailed regions.Learning-based methods take advantage of an external HR image dataset to induce a mapping function from LR to HR image, which contains some implicit priors. For estimating the HR image, we need to substitute the LR image into the mapping function as a variable. Freeman et al. [11] designed a Markov Network model with HR and LR image patches, for each LR patch extracted from the input LR image, which selected the ‘best proper’ HR patch from some candidates as the reconstructed result. In essence, the whole reconstructed HR image was stitched by these ‘best proper’ HR patches from the training set, which makes noise and blur generated easily near the stitched position. Because of the nearest neighbors searching step, time-consuming is another drawback for this algorithm. Chang et al. [12] introduced an idea called ‘Local Linear Embedding’ (LLE) originated from manifold learning. For each LR image patch, itsk nearest neighbors found in training set were used for fitting it. The fitting coefficients between the LR image patch and the corresponding HR one were assumed to be identical. Thus the unknown HR patches can be estimated by weighted combination of the HR versions sourced from thek nearest neighbors. This algorithm can successfully suppress the noise appeared near the stitched position and improve the reconstructed result to some extent. However, the algorithm needs to define the number of nearest neighborskin advance, which make it cannot be adapted according to the image content. Therefore, some noise and blur can be observed along the salient edge in the reconstructed result. Yang et al. [3, 4] considered SISR as sparse representation problem and constructed two over-completed dictionaries for LR and HR image jointly. For each patch extracted from the input LR image, the algorithm selected as little as possible ‘atoms’ in the LR dictionary to sparsely represent it and transferred the sparse coefficients to the HR dictionary to estimate the corresponding HR patch. Actually, this algorithm solved the problem of ‘fixedk’ in [12] and significantly improved the quality of the reconstructed results. However, the condition for obtaining the accurate sparse coefficients is greatly rigorous, which is relaxed to choose ‘atoms’ by an approximately random style, and this will make the discriminative atoms leak into the fitting step to generate more errors unavoidably. In order to choose some specific atoms when sparse represent a patch, Dong et al. [13] proposed an adaptive sparse domain selection algorithm, which preclassify the training set and construct a subdictionary for each class. For each super-resolved patch, it would select one subdictionary according to the distant between the patch and these subdictionaries to sparsely represent it, which can alleviate the problem of existing noise along salient edges to some extent. Timofte et al. [14] combined neighbor embedding and sparse representation to find the nearest atom in the over-completed dictionary and take advantage of thek nearest neighbors of this atom in the dictionary to build a matrix operator. The matrix operator can speedily map the LR patch into HR space. Based on the above theory, Timofte et al. [15] additionally explored the original training images to further improve the algorithm’s efficiency significantly. With the development of computer technology and the increasing enhancement of its ability to deal with data, deep learning as a resurgent technology is very popular in SISR [16–22]. However, Graphic Processing Unit (GPU) is a necessary accelerating installation for deep learning algorithm to optimize the weighted coefficients among neural network nodes by training a large number of HR images. Another weakness of deep learning technique is that prior knowledge is not effortless to be incorporated during the training process. For further improvement, the quality of the reconstructed HR image is restricted. As a summary, learning-based SISR algorithm exhibits huge potential for preserving edges and restoring details. However, this category relies too heavily on external dataset, which is a huge challenge for normal computer with the increasing complexity of SISR algorithm. The necessary number of HR images also shows an exponential rising tendency. Meanwhile, little research has been done on HR training images themselves, just to stack examples to enlarge the size of training set. For a specific super-resolved image, it is hardly to explain what roles an image from training set will play up to now. The only subjective conclusion is that an HR training image will supply more high frequency information when its content is more similar to the input one.As far as image content is concerned, the input LR image itself is one of the closest ones to the corresponding HR image. Therefore, a kind of special learning-based algorithm called self-similarity based SISR has been presented [23–27]. The self-similarity of an image refers to a phenomenon when we observe an image from small patches perspective; each small patch will recur at other position or other scales with higher probability [23–29]. The self-similarity based SISR algorithms with the help of self-similarity property take the input image itself as training example to estimating the HR image. Glasner et al. [23] proposed a unified framework which considers the similar patches from the same scale as the different views from the same scenario and build the training patch pairs with similar patches from different scales and then integrated multiframe super-resolution technology with learning-based approach to magnify the LR input image. For some specific patches, there could exist obvious differences between the input patch and its nearest neighbors, which will generate blur and noise near the edges inevitably or introduce incorrect high frequency components in detailed regions. Concerning this issue, Freedman et al. [24] and Yang et al. [25] pointed out, respectively, that almost all patches will recur themselves in a small window at the same position in the input or lower scale image. This property can make the consuming time for finding nearest neighbors reduced significantly. However, due to the limited size of the training set, the edges of reconstructed HR image look too sharp to be natural. For solving the problem of small training size, He et al. [30] constructed Gaussian process regression (GPR) training set to train parameters in a small window. However, this algorithm did not explore the relationship among different windows that leads to obvious noise and artifacts existed in the result. For the same purpose, Huang et al. [31] performed specific geometric deformation on the input image patch which not only extend the research space but also make the nearest neighbors matching more accurate. Due to the geometric deformation leading to distortion of the patch, there are visible noise and blur in the reconstructed HR image especially in complex texture region.Self-similarity based SISR algorithm generally constructs the training set just with the help of self-similarity property of LR image. However, the self-similarity property of the unknown HR image has not been noticed. In this paper, we consider the self-similarity of the unknown HR image as a prior and propose a new SISR algorithm. The main contributions are as follows.First, we utilize different scales of the input LR image to construct a joint Gaussian mixture model. Then an interesting prior is discovered by deviation with the help of this model. Combined with this discovery, we deliver a novel SISR algorithm, which is competitive to the state-of-the-art algorithms in both the quality of the reconstructed image and time cost.Second, our proposed algorithm needs two steps: training and reconstructing. In training step, we train a joint Gaussian mixture model that contains a small number of Gaussian components. In reconstructing step, our algorithm can avoid nearest neighbor searching which is time-consuming and liable to generate errors.At last, the presented prior in our algorithm can be self-adapted according to the input LR image, which make our algorithm more robust and can be uncomplicated expand into other fields such as image denoising, image restoration and so on.The rest of this paper is organized as follows. Section2 explains the prior term in SISR cost function. Section 3 details our proposed algorithm. The corresponding experiments and analysis are all exhibited in Section 4. Section 5 devotes to the conclusion of our work in this paper.


## 2. Prior in Cost Function

Imaging process will be disturbed by some degradation factors such as relative motion, focus on misalignment, atmosphere turbulence, etc. which is difficult to find an ideal mathematical model to explain. Actually, researchers generally describe the imaging procedure with a linear system model for SISR problem(1)Y=AX+n,whereY is the input LR observation,X is the unknown ideal HR image,A is the degraded matrix which should include filtering and down-sampling operations, andn is additive white Gaussian noise generated in imaging process. For an observed LR image, there are infinite number of solutions that can gratify the imaging model, consequently the problem of SISR is a representative ill-posed problem. For obtaining a satisfying solution, researchers often incorporate some priors to regularize the solution space. In terms of prior modality, there are primarily two ways to introduce them, explicitly, and implicitly, as shown in Figure 1. The explicit priors generated by summarizing external HR images are intuitive and comprehensible for their concrete expression formula. The SISR cost function combined with some explicit prior can be casted as follows: (2)argminX⁡Y-AX2+λρX,where Y-AX2 is the fidelity term, ρX is the prior, and λ is a balance coefficient to tradeoff between the two terms. Natural images, as a whole, vary in content and structure, which are difficult to establish with a unified model to represent. However, when researching them from local region perspective (such as image patch with size 5×5, 7×7 …), the content of each image patch becomes primitive and shows strong regularity, which makes modeling painless. Consequently, the prior term designed from the view of image patches can be written as(3)ρX=-log⁡∏i=1NpPiX=-∑i=1Nlog⁡pPiX,whereN is the total number of those extracted patches,Pi is an extracted matrix which takes charge of extractingi-th patch from HR imageX, and p(xi) is the prior ofi-th patch, where concrete expression of it will vary with the selected prior [28, 29]. On the whole, the formulated problem in formula (2) is nonconvex and difficult to solve directly. More typically, an auxiliary variable zi is introduced for eachxi, and then formula (3) can be solved by ‘Half Quadratic Splitting’ algorithm [32–37].Figure 1
Two styles of training prior (the above is implicit style, which have no specific prior representation, and the below is explicit style, and the concrete prior expression is trained by HR images.).


## 3. Maximizing Nonlocal Self-Similarity Prior for SISR

### 3.1. Problem Formulation

Self-similarity is an inherent property of natural image itself [23–25, 28, 29, 33]. Through experiments displayed in Figure 2, we can notice that what the clearer, sharper, and more detailed image are, what the more obvious this property is. The image patch extracted arbitrarily from the clear image has an extremely similar appearance to the nearest neighbor found in its lower scale (as shown in Figure 2(a)). However, in the image with noise (as shown in Figure 2(b)) or blur (as shown in Figure 2(c)) enforced, there are big differences between the selected patch and its nearest neighbors found in the corresponding lower scales. Based on this phenomenon, we believe that the clearer the image is, the stronger the self-similarity property manifests. In other words, on the premise that the reconstructed HR result is good enough, one image patch arbitrarily extracted from it will find similar patches at its lower scale with higher probability (in this paper, the lower scale refers to the input LR image). And the following is the concrete realization of this idea.Figure 2
Illustration of cross scale self-similarity property in natural image. (a) Small patches extracted from sharp HR image are similar to their nearest neighbors in lower scale. (b) and (c) In contrast, in a noise (b) or blurry (c) image, there are big differences between small patches extracted from degraded HR image and their nearest neighbors in lower scale. The patch positions can be found well in zoomed PDF.
(a)
Clean image (b)
Noise image (c)
Blurred imageSupported by self-similarity, we use the input LR image to construct the training and testing dataset. As shown in Figure3, first, upsample both the input LR imageY and its downsampling versionZ, s times to obtain the corresponding results X′ and Y′, andX is the unknown HR image for estimating. As a matter of fact,Y and Y′ can be considered as the degraded LR image version ofX and X′ separately. It can be noticed according to the self-similarity property, for patchesxi and xi′ separately extracted fromX and X′ that there should exist recurrences in the corresponding LR imagesY and Y′ with the highest probability. Second, we can formulate this problem as follows: extracting all patches formY and Y′, then combining them to obtain the training set yjT,yj′Tj=1M, whereyj and yj′ arej-th patches extracted fromY and Y′ separately.M is the total number of patches extracted fromY or Y′, andT represents the transpose operator on a vector. Next, we combine the patchesxi and xi′ which are thei-th patches extracted formX and X′ separately. According to the self-similarity theory, when xiT,xi′TT make recurrence in training set yjT,yj′Tj=1M with the highest probability, the residual core task is how to estimate the HR patchxi.Figure 3
The framework of constructing the training set.


### 3.2. Nonlocal Self-Similarity Prior

Recently, GMM has been a popular and valid tool used in image processing [38] due to the fact that it can smoothly approximate the density distribution of arbitrary shape and is characterized as remarkable efficiency. Specifically, this property shows excellent performance in image modeling based on image patches [39, 40]. In this paper, for training set yjT,yj′Tj=1M (yj∈Ra2and yj′∈Ra2, where a is an integer that represents the size of extracted patch a×a), we introduce GMM to approximate its probability density function, as shown in formula (4),(4)pyjT,yj′TT=∑k=1KπkNyjT,yj′TT∣μk,Σk,whereK is the number of Gaussian components in GMM, πk,μk,Σk is the weighted coefficient, mean vector and covariance matrix ofk-th Gaussian component separately, according to the definition of GMM, ∑k=1Kπk=1, and Nz∣μk,Σk representsk-th Gaussian component, and its concrete expression is (5)Nz∣μ,Σ=12πn/2Σ1/2exp⁡-12z-μTΣ-1z-μ.We take the expectation maximization (EM) algorithm to training the parameters πk,μk,Σk, as illustrated in Algorithm 1. The EM algorithm can be divided into two steps including E-step and M-step, which are executed alternatively. E-step refers to estimating the probability that thej-th data pointzj in the training dataset is generated by thek-th Gaussian component if the model parameters are already known, which can be described as(6)γj,k=πkNzj∣μk,Σk∑i=1KπiNzi∣μi,Σi.Thus, the probability of the whole dataset generated byk-th Gaussian component can be casted as (7)Nk=∑j=1Nγj,k.Next, the model parameters πk,μk,Σk of GMM can be estimated by maximum likelihood in M-step.(8)πk=NkN,(9)μk=1Nk∑j=1Nγj,kzj,(10)Σk=1Nkγj,kzj-μkzj-μkT.Because each example in training set was combined by two patches extracted fromY and Y′ separately, the mean vector μk and covariance matrix Σk can be rewritten as (11)μk=μkyμky′,(12)Σk=ΣkyΣkyy′Σky′yΣky′.Thus the probability of testing example xiT,xi′T in the probability density function of the training set can be represented as(13)pxi,xi′∣πk,μk,Σkk=1K=∑k=1KπkNxi,xi′∣μk,Σk.Thus, we discovered that the probability of patch vector xi satisfied a specific Gaussian distribution(14)pxi∣xi′,πk,μk,Σkk=1K~Nxi∣μxi,Σxi,where μxi and Σxi can be denoted respectively as mean vector and covariance matrix for specifying this Gaussian distribution. (15)μxi=∑k=1Kφkxi′μky+Σkyy′Σky′-1xi′-μky′,(16)Σxi=∑k=1Kφk2xi′Σky-Σkyy′Σky′-1Σky′y,And the coefficient φkxi′ appeared in formula (15) and formula (16) is (17)φkxi′=πkNxi′∣μky′,Σky′∑i=1KπiNxi′∣μiy′,Σiy′.Then the prior in this paper can be define as(18)ρX=-∑i=1Nlog⁡pPiX=-∑i=1Nlog⁡NPiX∣μxi,Σxi,whereNis the number of patches extracted from the unknown HR imageX andPi is an extracted matrix to extracti-th patch fromX. According to the prior expression, we can conclude that the recurrence probability of each patch in unknown HR image is measured, and we call this prior as nonlocal self-similarity prior.Algorithm 1:Estimate the parameters of GMM.
Input: small image Y, low frequency version Y′ Output: the parameters πk,μk,Σk 1: Function π k , μ k , Σ k = Training(Y,Y′) 2:Partitioning Y and Y′ into overlapped patches yjj=1n and yj′j=1n; 3:Combining the corresponding patches into one vector yjT,yj′TT; 4:Training the parameters πk,μk,Σk using EM algorithm; 5:  Return π k , μ k , Σ k; 6: end Function


### 3.3. Image Super-Resolution Reconstruction

In this paper, our algorithm requires the reconstructed HR result to be as similar as possible to the input LR image after the degeneration model processed; meanwhile the result also should satisfy the proposed nonlocal self-similarity prior. Therefore, the cost function of our algorithm is designed as(19)argminX⁡12Y-AX2-λ∑i=1Nlog⁡NPiX∣μxi,Σxi,where the first is fidelity term, the second is prior term, and λ is balance coefficient between the two terms. The cost function can be further simplified as (20)argminX⁡12Y-AX2+λ∑i=1Nlog⁡Σxi+12PiX-μxiTΣxi-1PiX-μxi.The solution of (20) is relatively simple to solve; we directly take the derivative of (20) and set the derivative to zero(21)ATAX-ATY+λ∑i=1NPiTΣxi-1PiX-∑i=1NPiTΣxi-1μxi=0.The final estimation of the HR image can be obtained by further sorting the expression(22)X=ATA+λ∑i=1NPiTΣxi-1Pi-1ATY+λ∑i=1NPiTΣxi-1μxi.The above idea can be implemented by Algorithm2 to estimate the whole HR image, and the complete framework of the proposed method is illustrated in Figure 4 to exhibit how to super-resolve a color image.Algorithm 2:SISR using nonlocal self-similarity prior.
Input: small image Y, magnification factor s Output:HR image X 1: X ′ = u p s a m p l i n g Y , s; 2: Y ′ = u p s a m p l i n g d o w n s a m p l i n g Y , s , s; 3:Partitioning Y,Y′ into overlapping patches yjj=1n and yj′j=1n; 4: π k , μ k , Σ k = T r a i n i n g Y , Y ′; 5: For each patch xi′ extracted form X′; 6:  For k = 1 : K 7:Computing φkxi′ using (17); 8:  End For 9:Computing μxi using (15); 10:Computing Σxi using (16); 11: End For 12:ComputingX using (22); 13: return XFigure 4
Framework of our proposed algorithm.


## 3.1. Problem Formulation

Self-similarity is an inherent property of natural image itself [23–25, 28, 29, 33]. Through experiments displayed in Figure 2, we can notice that what the clearer, sharper, and more detailed image are, what the more obvious this property is. The image patch extracted arbitrarily from the clear image has an extremely similar appearance to the nearest neighbor found in its lower scale (as shown in Figure 2(a)). However, in the image with noise (as shown in Figure 2(b)) or blur (as shown in Figure 2(c)) enforced, there are big differences between the selected patch and its nearest neighbors found in the corresponding lower scales. Based on this phenomenon, we believe that the clearer the image is, the stronger the self-similarity property manifests. In other words, on the premise that the reconstructed HR result is good enough, one image patch arbitrarily extracted from it will find similar patches at its lower scale with higher probability (in this paper, the lower scale refers to the input LR image). And the following is the concrete realization of this idea.Figure 2
Illustration of cross scale self-similarity property in natural image. (a) Small patches extracted from sharp HR image are similar to their nearest neighbors in lower scale. (b) and (c) In contrast, in a noise (b) or blurry (c) image, there are big differences between small patches extracted from degraded HR image and their nearest neighbors in lower scale. The patch positions can be found well in zoomed PDF.
(a)
Clean image (b)
Noise image (c)
Blurred imageSupported by self-similarity, we use the input LR image to construct the training and testing dataset. As shown in Figure3, first, upsample both the input LR imageY and its downsampling versionZ, s times to obtain the corresponding results X′ and Y′, andX is the unknown HR image for estimating. As a matter of fact,Y and Y′ can be considered as the degraded LR image version ofX and X′ separately. It can be noticed according to the self-similarity property, for patchesxi and xi′ separately extracted fromX and X′ that there should exist recurrences in the corresponding LR imagesY and Y′ with the highest probability. Second, we can formulate this problem as follows: extracting all patches formY and Y′, then combining them to obtain the training set yjT,yj′Tj=1M, whereyj and yj′ arej-th patches extracted fromY and Y′ separately.M is the total number of patches extracted fromY or Y′, andT represents the transpose operator on a vector. Next, we combine the patchesxi and xi′ which are thei-th patches extracted formX and X′ separately. According to the self-similarity theory, when xiT,xi′TT make recurrence in training set yjT,yj′Tj=1M with the highest probability, the residual core task is how to estimate the HR patchxi.Figure 3
The framework of constructing the training set.


## 3.2. Nonlocal Self-Similarity Prior

Recently, GMM has been a popular and valid tool used in image processing [38] due to the fact that it can smoothly approximate the density distribution of arbitrary shape and is characterized as remarkable efficiency. Specifically, this property shows excellent performance in image modeling based on image patches [39, 40]. In this paper, for training set yjT,yj′Tj=1M (yj∈Ra2and yj′∈Ra2, where a is an integer that represents the size of extracted patch a×a), we introduce GMM to approximate its probability density function, as shown in formula (4),(4)pyjT,yj′TT=∑k=1KπkNyjT,yj′TT∣μk,Σk,whereK is the number of Gaussian components in GMM, πk,μk,Σk is the weighted coefficient, mean vector and covariance matrix ofk-th Gaussian component separately, according to the definition of GMM, ∑k=1Kπk=1, and Nz∣μk,Σk representsk-th Gaussian component, and its concrete expression is (5)Nz∣μ,Σ=12πn/2Σ1/2exp⁡-12z-μTΣ-1z-μ.We take the expectation maximization (EM) algorithm to training the parameters πk,μk,Σk, as illustrated in Algorithm 1. The EM algorithm can be divided into two steps including E-step and M-step, which are executed alternatively. E-step refers to estimating the probability that thej-th data pointzj in the training dataset is generated by thek-th Gaussian component if the model parameters are already known, which can be described as(6)γj,k=πkNzj∣μk,Σk∑i=1KπiNzi∣μi,Σi.Thus, the probability of the whole dataset generated byk-th Gaussian component can be casted as (7)Nk=∑j=1Nγj,k.Next, the model parameters πk,μk,Σk of GMM can be estimated by maximum likelihood in M-step.(8)πk=NkN,(9)μk=1Nk∑j=1Nγj,kzj,(10)Σk=1Nkγj,kzj-μkzj-μkT.Because each example in training set was combined by two patches extracted fromY and Y′ separately, the mean vector μk and covariance matrix Σk can be rewritten as (11)μk=μkyμky′,(12)Σk=ΣkyΣkyy′Σky′yΣky′.Thus the probability of testing example xiT,xi′T in the probability density function of the training set can be represented as(13)pxi,xi′∣πk,μk,Σkk=1K=∑k=1KπkNxi,xi′∣μk,Σk.Thus, we discovered that the probability of patch vector xi satisfied a specific Gaussian distribution(14)pxi∣xi′,πk,μk,Σkk=1K~Nxi∣μxi,Σxi,where μxi and Σxi can be denoted respectively as mean vector and covariance matrix for specifying this Gaussian distribution. (15)μxi=∑k=1Kφkxi′μky+Σkyy′Σky′-1xi′-μky′,(16)Σxi=∑k=1Kφk2xi′Σky-Σkyy′Σky′-1Σky′y,And the coefficient φkxi′ appeared in formula (15) and formula (16) is (17)φkxi′=πkNxi′∣μky′,Σky′∑i=1KπiNxi′∣μiy′,Σiy′.Then the prior in this paper can be define as(18)ρX=-∑i=1Nlog⁡pPiX=-∑i=1Nlog⁡NPiX∣μxi,Σxi,whereNis the number of patches extracted from the unknown HR imageX andPi is an extracted matrix to extracti-th patch fromX. According to the prior expression, we can conclude that the recurrence probability of each patch in unknown HR image is measured, and we call this prior as nonlocal self-similarity prior.Algorithm 1:Estimate the parameters of GMM.
Input: small image Y, low frequency version Y′ Output: the parameters πk,μk,Σk 1: Function π k , μ k , Σ k = Training(Y,Y′) 2:Partitioning Y and Y′ into overlapped patches yjj=1n and yj′j=1n; 3:Combining the corresponding patches into one vector yjT,yj′TT; 4:Training the parameters πk,μk,Σk using EM algorithm; 5:  Return π k , μ k , Σ k; 6: end Function


## 3.3. Image Super-Resolution Reconstruction

In this paper, our algorithm requires the reconstructed HR result to be as similar as possible to the input LR image after the degeneration model processed; meanwhile the result also should satisfy the proposed nonlocal self-similarity prior. Therefore, the cost function of our algorithm is designed as(19)argminX⁡12Y-AX2-λ∑i=1Nlog⁡NPiX∣μxi,Σxi,where the first is fidelity term, the second is prior term, and λ is balance coefficient between the two terms. The cost function can be further simplified as (20)argminX⁡12Y-AX2+λ∑i=1Nlog⁡Σxi+12PiX-μxiTΣxi-1PiX-μxi.The solution of (20) is relatively simple to solve; we directly take the derivative of (20) and set the derivative to zero(21)ATAX-ATY+λ∑i=1NPiTΣxi-1PiX-∑i=1NPiTΣxi-1μxi=0.The final estimation of the HR image can be obtained by further sorting the expression(22)X=ATA+λ∑i=1NPiTΣxi-1Pi-1ATY+λ∑i=1NPiTΣxi-1μxi.The above idea can be implemented by Algorithm2 to estimate the whole HR image, and the complete framework of the proposed method is illustrated in Figure 4 to exhibit how to super-resolve a color image.Algorithm 2:SISR using nonlocal self-similarity prior.
Input: small image Y, magnification factor s Output:HR image X 1: X ′ = u p s a m p l i n g Y , s; 2: Y ′ = u p s a m p l i n g d o w n s a m p l i n g Y , s , s; 3:Partitioning Y,Y′ into overlapping patches yjj=1n and yj′j=1n; 4: π k , μ k , Σ k = T r a i n i n g Y , Y ′; 5: For each patch xi′ extracted form X′; 6:  For k = 1 : K 7:Computing φkxi′ using (17); 8:  End For 9:Computing μxi using (15); 10:Computing Σxi using (16); 11: End For 12:ComputingX using (22); 13: return XFigure 4
Framework of our proposed algorithm.


## 4. Experiments and Analyses

In this section, we verify the efficiency of our nonlocal self-similarity prior super-resolution algorithm by subjective visual observation and objective parameters comparisons. We firstly elaborate the parameters setting which is necessary in implementing our algorithm and then give the comparison results including the quality of the reconstructed HR images and their consuming time. At last, the performance of our algorithm is further discussed.


### 4.1. Experiment Setting

We regard BSD500 dataset and some commonly used test images in image processing field as the ideal HR images and all these HR images were downloaded from the websitehttps://www2.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/resources.html and http://people.csail.mit.edu/mrub/retargetme/download.html. In order to simulate the image degradation model to generate LR images for testing input, these HR images were downsampleds times by bicubic method and then filtered by Gaussian filtering with mean value of 0 and variance of 0.5. Finally, we input these LR images into our proposed algorithm to upsample them s times. As the self-similarity between two different scales of one image will be weaken with the decreasing scale [20, 25, 26, 30, 35], when the magnification factor is larger, zooming the input LR image to the target size directly will not achieve the best performance for the reconstructed HR image. Aiming at this problem, we adopted a stepwise super-resolving method, and the magnification factor is1.25 for each time. For the input color image, we firstly transformed this image fromRGB space toYUV space, because channelY representing brightness is sensitive to human eyes. The channelY is processed by our algorithm. TheU andV channels representing color information are relatively sluggish to human eyes. And theU andV channels can be magnified to the target size directly by bicubic method. Finally, the three reconstructed channel images are merged and transformed fromYUV space toRGB space for displaying. In addition, the parameters involved in our algorithm are set as λ=0.5 (as in formula (22)), a=7 (the patch size extracted form image is 7×7), and K=8. These two parameters λ andK are determined empirically and will be explained latter. In the experiment, we noticed that, for almost all the Patches xi′i=1N, their corresponding weights φkxi′k=1K were almost sparse, such as the example shown in Figure 5. As a consequence, for reducing the amount of calculation, formula (23) was used and φkxi′ was set to(23)φkxi′=φkxi′,φkxi′≥K-10φkxi′≺K-1.Figure 5
Lenna example to illustrate the sparsity ofφkxi′. Left is image lenna and its three randomly extracted patches. The right is the value of φkxi′ whereK is set to 8.


### 4.2. Experiment Results

In order to verify the effectiveness of our proposed algorithm, we chose some state-of-the-art algorithms and conducted a series of comparison experiments. The selected algorithms for comparison include ScSR [4], Glasner [23], GPR [30], SPM [41], SelfExSR [31], A+ [15], SRCNN [17], and FSRCNN [18], where Glasner, GPR, SelfExSR, and our proposed are self-similarity based super-resolution algorithms. The codes of all these algorithms except Glasner were downloaded from the authors own websites separately. The code of Glasner was achieved by ourselves, and the results generated by this code are very similar to those image displayed in the original paper [23]. All these codes were designed by Matlab software. All these algorithms were executed in inter (R) core (TM) i7-5600 CPU @ 2.60 GHz, 8.0 GM cache and windows 7 professional, 64 bit operating system with Matlab 2015b. As shown in Figures 6, 7, and 8, these three groups of images are the comparison results of these algorithms that are magnified by 2 times, 3 times, and 4 times separately. Observing these three groups of images carefully, we can notice that the algorithms of ScSR, Glasner, GPR, and SPM did not preserve the edges well and made noise, artifacts, and distortion along the salient edges in the reconstructed HR images. Meanwhile, the generated noise and artifacts become more obvious with the increasing of the image magnification factor in these methods. On the contrary, the algorithms of selfExSR, A+, SRCNN, FSRCNN, and our proposed have a good effect on the edge and texture region in the images, and with the increasing of magnification factor, there is no significant noise and artifacts generated near the edges or in the texture region. Furthermore, by observing the smallest windows marked in the reconstructed results, it can be found that our method has the smallest deformation near the edge, which is closest to the original image and shows excellent texture processing capability. In order to further illustrate the performance of this algorithm, we randomly selected 10 images from ideal HR image dataset, as shown in Figure 9, which were super-resolved by these selected state-of-the-art algorithms and our proposed with 3 and 4 times separately. Then, the parameters of peak single to noise ratio (PSNR) and structure similarity index (SSIM) of each reconstructed image were calculated. PSNR and SSIM are two parameters to measure the image quality under laboratory conditions. Generally speaking, the higher PSNR and SSIM values represent the better quality of the reconstructed results. These PSNR and SSIM values of the reconstructed images are shown in Tables 1 and 2, where for each reconstructed image, the first line is PSNR and the second line is SSIM. By observing the two tables, it can be concluded that the quality of our proposed algorithm is superior to others in a majority of results.Table 1
Quality evaluation results of 10 test images with super resolution 3×.
Image ScSR Glasner GPR SPM SelfExSR A+ SRCNN FSRCNN Our Fence 28.51 28.12 27.72 27.77 29.06 28.69 28.99 29.17 29.41 0.87 0.80 0.79 0.85 0.88 0.88 0.89 0.88 0.91 Building 28.79 27.72 27.57 28.55 30.83 29.49 29.39 30.36 30.47 0.93 0.85 0.84 0.90 0.95 0.94 0.93 0.94 0.95 Baby 36.99 32.23 32.24 36.55 37.20 37.14 37.14 37.20 37.55 0.96 0.92 0.92 0.96 0.96 0.96 0.96 0.96 0.97 Barbara 26.59 25.80 25.59 25.95 26.33 26.72 26.70 26.29 26.69 0.85 0.77 0.73 0.83 0.85 0.86 0.86 0.85 0.88 Lena 34.88 34.37 33.19 33.69 35.16 35.05 35.14 35.48 35.21 0.93 0.86 0.86 0.92 0.93 0.93 0.93 0.94 0.93 Tiger 28.94 28.22 27.43 28.39 29.00 29.06 29.12 29.28 29.48 0.80 0.70 0.67 0.76 0.80 0.80 0.80 0.81 0.83 Kola 34.01 33.93 32.11 33.25 33.98 34.07 34.17 34.27 34.43 0.93 0.83 0.83 0.91 0.93 0.93 0.93 0.94 0.95 Butterfly 31.47 31.92 30.70 30.87 31.83 32.05 32.28 32.63 32.58 0.96 0.89 0.91 0.93 0.96 0.97 0.96 0.96 0.96 Parrot 32.69 31.36 30.28 32.20 33.02 33.08 33.39 34.02 34.24 0.95 0.90 0.90 0.94 0.95 0.95 0.95 0.95 0.96 Zebra 30.76 29.77 29.25 30.71 31.67 31.01 31.74 31.17 31.33 0.95 0.83 0.82 0.93 0.95 0.95 0.96 0.95 0.95 Average 31.36 30.34 29.61 30.79 31.51 31.64 31.71 31.99 32.14 0.91 0.84 0.82 0.89 0.92 0.92 0.92 0.93 0.93Table 2
Quality evaluation results of 10 test images with super resolution 4.
Image ScSR Glasner GPR SPM SelfExSR A+ SRCNN FSRCNN Our Fence 26.83 26.00 26.84 26.98 27.65 26.94 27.59 27.78 27.86 0.74 0.66 0.71 0.75 0.76 0.74 0.76 0.78 0.79 Building 26.38 26.53 26.06 26.13 27.03 27.27 27.24 27.84 27.88 0.74 0.64 0.73 0.75 0.81 0.78 0.78 0.80 0.81 Baby 33.20 33.86 33.04 33.62 34.81 34.73 34.82 34.89 34.97 0.91 0.86 0.92 0.93 0.93 0.93 0.93 0.87 0.94 Barbara 23.82 23.59 23.78 23.83 24.04 24.03 24.05 24.51 25.84 0.70 0.65 0.71 0.72 0.73 0.72 0.73 0.75 0.77 Lena 30.06 29.94 29.00 29.85 31.90 31.76 31.87 31.60 31.79 0.85 0.76 0.85 0.87 0.91 0.89 0.89 0.86 0.89 Tiger 26.43 26.42 26.58 26.35 26.79 26.81 26.89 26.28 27.11 0.63 0.62 0.62 0.65 0.65 0.65 0.65 0.63 0.68 Kola 28.04 25.31 28.04 28.06 28.63 28.68 28.74 28.10 28.92 0.75 0.63 0.74 0.78 0.77 0.77 0.76 0.73 0.78 Butterfly 25.42 25.77 25.50 24.66 28.77 29.40 29.53 29.66 29.55 0.80 0.81 0.89 0.75 0.93 0.97 0.92 0.96 0.95 Parrot 27.85 27.64 27.37 27.84 29.18 28.93 29.28 29.52 29.59 0.88 0.88 0.87 0.89 0.91 0.90 0.90 0.93 0.94 Zebra 26.51 25.95 25.12 27.04 28.87 28.46 28.70 28.55 28.94 0.81 0.81 0.80 0.84 0.87 0.86 0.86 0.84 0.88 Average 27.45 27.10 27.13 27.44 28.77 28.70 28.87 28.87 29.25 0.78 0.73 0.78 0.79 0.83 0.82 0.82 0.82 0.84Figure 6
Super-resolved image (2×) of ‘temple’.Figure 7
Super-resolved image (3×) of ‘penguin’.Figure 8
Super-resolved image (4×) of ‘Buddha’.Figure 9
Ten test images for effectiveness validation.In practical, the number of Gaussian componentsKin GMM has to be set manually in our algorithm. In order to reveal the relationship between Gaussian component and reconstructed results, we magnified the 10 testing images with 3 times when the number of Gaussian components is set to be 3,4…,11,12, respectively, and then computed the average PSNR and SSIM values, and the results are shown in Figure 10. As can be seen from Figure 10, when K=8, both the PSNR and SSIM values achieved their maximum separately, which mean the image quality is best at this point. How to set the parameter λ is a key factor affecting the quality of the reconstructed high-resolution image which balances between the fidelity term and the prior term. We randomly selected 10 high-resolution images from the training set to be downsampled 3 times as the inputs, then set the parameter λ as 0.1,0.2…,1.0, and keep other parameters to be fixed to super-resolve them separately and finally compute the average PSNR and SSIM for each λ about the 10 images. The result was illustrated below in Figure 11. We can discover from Figure 11 that the parameters both PSNR and SSIM were maximized at λ=0.5. For verifying the effectiveness of our algorithm from statistical perspective, we randomly selected 200 images from BSD500 as testing images firstly which were super-resolved by these state-of-the-art algorithms and ours with 3 times. Then, calculate and restore these PSNR values for each reconstructed HR image and then count the HR results which entered into the top three maximums for each algorithm. The process of SSIM parameter is the same as PSNR. The statistical results were shown in Figure 12 where these algorithms were listed horizontally, and the number of best quality images for each algorithm was displayed vertically. It is effortless to find that the best image numbers of SRCNN, FSRCNN, and ours are more than half of the testing dataset, and our algorithm is the most one.Figure 10
The relationship between the number of Gaussian componentK and the quality of reconstructed images (when K=8, the quality of the reconstructed image is best).Figure 11
The relationship between the balance weightλ and the quality of reconstructed images (when λ=0.5, the quality of the reconstructed image is best).Figure 12
Statistical experiment with BSD500.At last, we compared the consuming time of our algorithm with others. We also selected 200 images from BSD500 randomly, super-resolved them with those algorithms, and then recorded the average PSNR and the consuming time for each algorithm. The average PSNR and consuming time for each algorithm were shown in Figure13. The X-coordinate represents the PSNR and the Y-coordinate represents the consuming time in Figure 13. The conclusion is that the average PSNR of our algorithm is the largest one, and its consuming time lies at midlevel which is less than GPR and self-EXSR.Figure 13
The comparison of consuming time and PSNR among state-of-the-art algorithms and ours.


### 4.3. Experimental Analysis

From the above experimental results, we can find that our algorithm has certain advantages in the quality of reconstructed image compared with other state-of-the-art algorithms, and the consuming time is at the intermediate level. By literature retrieval, we noticed that there are some formal resemblances between the algorithms [32, 33, 38, 40, 42] and our proposed, but they are different in nature. The algorithms [32, 38, 40, 42] achieved the task of image restoration with the help of GMM. For an image patch to be restored, the assumption that this image patch is generated by just some Gaussian component with highest probability is necessary, and then to estimate the clear version of this patch by wiener filtering. However, our algorithm did not need such assumption, and the cost function of our algorithm has a closed form solution. In paper [33], a prior that ‘clearer image has stronger self-similarity’ is used to deblur the natural image. This prior proposed in [33] is similar to the nonlocal self-similarity algorithm proposed in this paper, but the expressions of the two priors are completely different. The prior proposed in [33] directly calculated the distance between the blurred image patch and its nearest neighbor found in the lower scale of the input image. The result of this estimation is essentially the weighted combination of these nearest neighbors of the blurred patch. However, our proposed algorithm is to estimate the HR image patch with maximum value of the probability density function. Meanwhile, the algorithm in [33] still needs the step of searching nearest neighbors which is time-consuming and inaccurate; however, the algorithm in this paper can avoid this step. The algorithm presented in [40] is formally closest to our proposed algorithm; both paper [40] and our proposed algorithm resorted to probability density function trained with LR/HR image patch pairs to estimate HR image. However, the algorithm presented in [40] still needs the restricted assumption and the wiener filtering solution used in papers [32, 38, 42] to estimate. More importantly, the wiener filtering solution is equivalent to the mean vector μxiof the specific Gaussian distribution prior of the algorithm in this paper, and the information of its covariance matrix Σxi is not used in paper [40], which reduces the accuracy of the results. Based on the experimental results presented in paper [40], our algorithm is far more efficient because it exploits the information of covariance matrix Σxi.The time complexity of one algorithm is an important index to evaluate its quality. The algorithm in this paper is mainly composed of two parts: model training and super-resolution. EM algorithm was used to train the Gaussian mixture model, and the time complexity isO(n2), wheren is the number of low resolution and high-resolution patch pairs. Since the algorithm in this paper only uses the input image itself for training, the number of image pairs will not be too large to consume unacceptable time. The super-resolution part would take up most of the consuming time of the algorithm, since a large number of matrix inversion and matrix multiplication operations were necessary in this step, and the time complexity of this step is O(n3). However, the weight φkxi′ is sparse, which make a large number of μxi and Σxi set to be 0 directly when φkxi′ is very small or even 0. PiTΣxi-1 was used in both of the two terms of (22) which can be computed one time and then restored the result for the second term. As a whole, the total consuming time by the proposed algorithm is acceptable.


## 4.1. Experiment Setting

We regard BSD500 dataset and some commonly used test images in image processing field as the ideal HR images and all these HR images were downloaded from the websitehttps://www2.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/resources.html and http://people.csail.mit.edu/mrub/retargetme/download.html. In order to simulate the image degradation model to generate LR images for testing input, these HR images were downsampleds times by bicubic method and then filtered by Gaussian filtering with mean value of 0 and variance of 0.5. Finally, we input these LR images into our proposed algorithm to upsample them s times. As the self-similarity between two different scales of one image will be weaken with the decreasing scale [20, 25, 26, 30, 35], when the magnification factor is larger, zooming the input LR image to the target size directly will not achieve the best performance for the reconstructed HR image. Aiming at this problem, we adopted a stepwise super-resolving method, and the magnification factor is1.25 for each time. For the input color image, we firstly transformed this image fromRGB space toYUV space, because channelY representing brightness is sensitive to human eyes. The channelY is processed by our algorithm. TheU andV channels representing color information are relatively sluggish to human eyes. And theU andV channels can be magnified to the target size directly by bicubic method. Finally, the three reconstructed channel images are merged and transformed fromYUV space toRGB space for displaying. In addition, the parameters involved in our algorithm are set as λ=0.5 (as in formula (22)), a=7 (the patch size extracted form image is 7×7), and K=8. These two parameters λ andK are determined empirically and will be explained latter. In the experiment, we noticed that, for almost all the Patches xi′i=1N, their corresponding weights φkxi′k=1K were almost sparse, such as the example shown in Figure 5. As a consequence, for reducing the amount of calculation, formula (23) was used and φkxi′ was set to(23)φkxi′=φkxi′,φkxi′≥K-10φkxi′≺K-1.Figure 5
Lenna example to illustrate the sparsity ofφkxi′. Left is image lenna and its three randomly extracted patches. The right is the value of φkxi′ whereK is set to 8.


## 4.2. Experiment Results

In order to verify the effectiveness of our proposed algorithm, we chose some state-of-the-art algorithms and conducted a series of comparison experiments. The selected algorithms for comparison include ScSR [4], Glasner [23], GPR [30], SPM [41], SelfExSR [31], A+ [15], SRCNN [17], and FSRCNN [18], where Glasner, GPR, SelfExSR, and our proposed are self-similarity based super-resolution algorithms. The codes of all these algorithms except Glasner were downloaded from the authors own websites separately. The code of Glasner was achieved by ourselves, and the results generated by this code are very similar to those image displayed in the original paper [23]. All these codes were designed by Matlab software. All these algorithms were executed in inter (R) core (TM) i7-5600 CPU @ 2.60 GHz, 8.0 GM cache and windows 7 professional, 64 bit operating system with Matlab 2015b. As shown in Figures 6, 7, and 8, these three groups of images are the comparison results of these algorithms that are magnified by 2 times, 3 times, and 4 times separately. Observing these three groups of images carefully, we can notice that the algorithms of ScSR, Glasner, GPR, and SPM did not preserve the edges well and made noise, artifacts, and distortion along the salient edges in the reconstructed HR images. Meanwhile, the generated noise and artifacts become more obvious with the increasing of the image magnification factor in these methods. On the contrary, the algorithms of selfExSR, A+, SRCNN, FSRCNN, and our proposed have a good effect on the edge and texture region in the images, and with the increasing of magnification factor, there is no significant noise and artifacts generated near the edges or in the texture region. Furthermore, by observing the smallest windows marked in the reconstructed results, it can be found that our method has the smallest deformation near the edge, which is closest to the original image and shows excellent texture processing capability. In order to further illustrate the performance of this algorithm, we randomly selected 10 images from ideal HR image dataset, as shown in Figure 9, which were super-resolved by these selected state-of-the-art algorithms and our proposed with 3 and 4 times separately. Then, the parameters of peak single to noise ratio (PSNR) and structure similarity index (SSIM) of each reconstructed image were calculated. PSNR and SSIM are two parameters to measure the image quality under laboratory conditions. Generally speaking, the higher PSNR and SSIM values represent the better quality of the reconstructed results. These PSNR and SSIM values of the reconstructed images are shown in Tables 1 and 2, where for each reconstructed image, the first line is PSNR and the second line is SSIM. By observing the two tables, it can be concluded that the quality of our proposed algorithm is superior to others in a majority of results.Table 1
Quality evaluation results of 10 test images with super resolution 3×.
Image ScSR Glasner GPR SPM SelfExSR A+ SRCNN FSRCNN Our Fence 28.51 28.12 27.72 27.77 29.06 28.69 28.99 29.17 29.41 0.87 0.80 0.79 0.85 0.88 0.88 0.89 0.88 0.91 Building 28.79 27.72 27.57 28.55 30.83 29.49 29.39 30.36 30.47 0.93 0.85 0.84 0.90 0.95 0.94 0.93 0.94 0.95 Baby 36.99 32.23 32.24 36.55 37.20 37.14 37.14 37.20 37.55 0.96 0.92 0.92 0.96 0.96 0.96 0.96 0.96 0.97 Barbara 26.59 25.80 25.59 25.95 26.33 26.72 26.70 26.29 26.69 0.85 0.77 0.73 0.83 0.85 0.86 0.86 0.85 0.88 Lena 34.88 34.37 33.19 33.69 35.16 35.05 35.14 35.48 35.21 0.93 0.86 0.86 0.92 0.93 0.93 0.93 0.94 0.93 Tiger 28.94 28.22 27.43 28.39 29.00 29.06 29.12 29.28 29.48 0.80 0.70 0.67 0.76 0.80 0.80 0.80 0.81 0.83 Kola 34.01 33.93 32.11 33.25 33.98 34.07 34.17 34.27 34.43 0.93 0.83 0.83 0.91 0.93 0.93 0.93 0.94 0.95 Butterfly 31.47 31.92 30.70 30.87 31.83 32.05 32.28 32.63 32.58 0.96 0.89 0.91 0.93 0.96 0.97 0.96 0.96 0.96 Parrot 32.69 31.36 30.28 32.20 33.02 33.08 33.39 34.02 34.24 0.95 0.90 0.90 0.94 0.95 0.95 0.95 0.95 0.96 Zebra 30.76 29.77 29.25 30.71 31.67 31.01 31.74 31.17 31.33 0.95 0.83 0.82 0.93 0.95 0.95 0.96 0.95 0.95 Average 31.36 30.34 29.61 30.79 31.51 31.64 31.71 31.99 32.14 0.91 0.84 0.82 0.89 0.92 0.92 0.92 0.93 0.93Table 2
Quality evaluation results of 10 test images with super resolution 4.
Image ScSR Glasner GPR SPM SelfExSR A+ SRCNN FSRCNN Our Fence 26.83 26.00 26.84 26.98 27.65 26.94 27.59 27.78 27.86 0.74 0.66 0.71 0.75 0.76 0.74 0.76 0.78 0.79 Building 26.38 26.53 26.06 26.13 27.03 27.27 27.24 27.84 27.88 0.74 0.64 0.73 0.75 0.81 0.78 0.78 0.80 0.81 Baby 33.20 33.86 33.04 33.62 34.81 34.73 34.82 34.89 34.97 0.91 0.86 0.92 0.93 0.93 0.93 0.93 0.87 0.94 Barbara 23.82 23.59 23.78 23.83 24.04 24.03 24.05 24.51 25.84 0.70 0.65 0.71 0.72 0.73 0.72 0.73 0.75 0.77 Lena 30.06 29.94 29.00 29.85 31.90 31.76 31.87 31.60 31.79 0.85 0.76 0.85 0.87 0.91 0.89 0.89 0.86 0.89 Tiger 26.43 26.42 26.58 26.35 26.79 26.81 26.89 26.28 27.11 0.63 0.62 0.62 0.65 0.65 0.65 0.65 0.63 0.68 Kola 28.04 25.31 28.04 28.06 28.63 28.68 28.74 28.10 28.92 0.75 0.63 0.74 0.78 0.77 0.77 0.76 0.73 0.78 Butterfly 25.42 25.77 25.50 24.66 28.77 29.40 29.53 29.66 29.55 0.80 0.81 0.89 0.75 0.93 0.97 0.92 0.96 0.95 Parrot 27.85 27.64 27.37 27.84 29.18 28.93 29.28 29.52 29.59 0.88 0.88 0.87 0.89 0.91 0.90 0.90 0.93 0.94 Zebra 26.51 25.95 25.12 27.04 28.87 28.46 28.70 28.55 28.94 0.81 0.81 0.80 0.84 0.87 0.86 0.86 0.84 0.88 Average 27.45 27.10 27.13 27.44 28.77 28.70 28.87 28.87 29.25 0.78 0.73 0.78 0.79 0.83 0.82 0.82 0.82 0.84Figure 6
Super-resolved image (2×) of ‘temple’.Figure 7
Super-resolved image (3×) of ‘penguin’.Figure 8
Super-resolved image (4×) of ‘Buddha’.Figure 9
Ten test images for effectiveness validation.In practical, the number of Gaussian componentsKin GMM has to be set manually in our algorithm. In order to reveal the relationship between Gaussian component and reconstructed results, we magnified the 10 testing images with 3 times when the number of Gaussian components is set to be 3,4…,11,12, respectively, and then computed the average PSNR and SSIM values, and the results are shown in Figure 10. As can be seen from Figure 10, when K=8, both the PSNR and SSIM values achieved their maximum separately, which mean the image quality is best at this point. How to set the parameter λ is a key factor affecting the quality of the reconstructed high-resolution image which balances between the fidelity term and the prior term. We randomly selected 10 high-resolution images from the training set to be downsampled 3 times as the inputs, then set the parameter λ as 0.1,0.2…,1.0, and keep other parameters to be fixed to super-resolve them separately and finally compute the average PSNR and SSIM for each λ about the 10 images. The result was illustrated below in Figure 11. We can discover from Figure 11 that the parameters both PSNR and SSIM were maximized at λ=0.5. For verifying the effectiveness of our algorithm from statistical perspective, we randomly selected 200 images from BSD500 as testing images firstly which were super-resolved by these state-of-the-art algorithms and ours with 3 times. Then, calculate and restore these PSNR values for each reconstructed HR image and then count the HR results which entered into the top three maximums for each algorithm. The process of SSIM parameter is the same as PSNR. The statistical results were shown in Figure 12 where these algorithms were listed horizontally, and the number of best quality images for each algorithm was displayed vertically. It is effortless to find that the best image numbers of SRCNN, FSRCNN, and ours are more than half of the testing dataset, and our algorithm is the most one.Figure 10
The relationship between the number of Gaussian componentK and the quality of reconstructed images (when K=8, the quality of the reconstructed image is best).Figure 11
The relationship between the balance weightλ and the quality of reconstructed images (when λ=0.5, the quality of the reconstructed image is best).Figure 12
Statistical experiment with BSD500.At last, we compared the consuming time of our algorithm with others. We also selected 200 images from BSD500 randomly, super-resolved them with those algorithms, and then recorded the average PSNR and the consuming time for each algorithm. The average PSNR and consuming time for each algorithm were shown in Figure13. The X-coordinate represents the PSNR and the Y-coordinate represents the consuming time in Figure 13. The conclusion is that the average PSNR of our algorithm is the largest one, and its consuming time lies at midlevel which is less than GPR and self-EXSR.Figure 13
The comparison of consuming time and PSNR among state-of-the-art algorithms and ours.


## 4.3. Experimental Analysis

From the above experimental results, we can find that our algorithm has certain advantages in the quality of reconstructed image compared with other state-of-the-art algorithms, and the consuming time is at the intermediate level. By literature retrieval, we noticed that there are some formal resemblances between the algorithms [32, 33, 38, 40, 42] and our proposed, but they are different in nature. The algorithms [32, 38, 40, 42] achieved the task of image restoration with the help of GMM. For an image patch to be restored, the assumption that this image patch is generated by just some Gaussian component with highest probability is necessary, and then to estimate the clear version of this patch by wiener filtering. However, our algorithm did not need such assumption, and the cost function of our algorithm has a closed form solution. In paper [33], a prior that ‘clearer image has stronger self-similarity’ is used to deblur the natural image. This prior proposed in [33] is similar to the nonlocal self-similarity algorithm proposed in this paper, but the expressions of the two priors are completely different. The prior proposed in [33] directly calculated the distance between the blurred image patch and its nearest neighbor found in the lower scale of the input image. The result of this estimation is essentially the weighted combination of these nearest neighbors of the blurred patch. However, our proposed algorithm is to estimate the HR image patch with maximum value of the probability density function. Meanwhile, the algorithm in [33] still needs the step of searching nearest neighbors which is time-consuming and inaccurate; however, the algorithm in this paper can avoid this step. The algorithm presented in [40] is formally closest to our proposed algorithm; both paper [40] and our proposed algorithm resorted to probability density function trained with LR/HR image patch pairs to estimate HR image. However, the algorithm presented in [40] still needs the restricted assumption and the wiener filtering solution used in papers [32, 38, 42] to estimate. More importantly, the wiener filtering solution is equivalent to the mean vector μxiof the specific Gaussian distribution prior of the algorithm in this paper, and the information of its covariance matrix Σxi is not used in paper [40], which reduces the accuracy of the results. Based on the experimental results presented in paper [40], our algorithm is far more efficient because it exploits the information of covariance matrix Σxi.The time complexity of one algorithm is an important index to evaluate its quality. The algorithm in this paper is mainly composed of two parts: model training and super-resolution. EM algorithm was used to train the Gaussian mixture model, and the time complexity isO(n2), wheren is the number of low resolution and high-resolution patch pairs. Since the algorithm in this paper only uses the input image itself for training, the number of image pairs will not be too large to consume unacceptable time. The super-resolution part would take up most of the consuming time of the algorithm, since a large number of matrix inversion and matrix multiplication operations were necessary in this step, and the time complexity of this step is O(n3). However, the weight φkxi′ is sparse, which make a large number of μxi and Σxi set to be 0 directly when φkxi′ is very small or even 0. PiTΣxi-1 was used in both of the two terms of (22) which can be computed one time and then restored the result for the second term. As a whole, the total consuming time by the proposed algorithm is acceptable.


## 5. Conclusion

In this paper, the property that the self-similarity of clear natural images is stronger than other versions is considered as a prior and integrated into the image super-resolution reconstruction framework. With the help of training joint GMM using LR and HR patch pairs, the prior is represented as a special Gaussian distribution. The algorithm in this paper has a closed solution, which does not require sophisticated training and time-consuming nearest neighbor searching. Therefore, the quality of reconstructed HR image with our proposed is better than other state-of-the-art algorithms. In the following future work, we hope to combine the image sparsity with the prior presented in this paper to further improve the quality of the reconstructed images.



---

*Source: 3840285-2019-02-28.xml*