Image denoising based on global image similar patches searching and HOSVD to patches tensor

Nonlocal self-similarity of natural image is an essential property for image processing. But how to measure the similarity between different patches and how to better exploit the similarity of patches are two crucial problems for image denoising. In this article, we establish a novel image denoising method based on a global image similar patches searching method and the similar patches tensor high-order singular value decomposition theory. Particularly, in order to find the reasonable similar patches to a reference, a variant of Gaussian mixture model (GMM) global similar patches searching method is proposed, and the graph process unit-based GMM model training method is performed to speed up the training process. Furthermore, the k-means and the local inter-class searching are used to improve on the similarity measure between patches. To better exploit the similarity of patches to image denoising, we rearrange similar patches to a tensor for each reference, and an iterative adaptive weighted tensor low-rank approximation method is established to perform image denoising. Experimental results clearly show that the proposed method is comparable to many recent denoising algorithms from the PNSR, SSIM and NCC viewpoint. For higher noise level (e.g., sigma=50\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$sigma=50$$\end{document}), the performance of our proposed algorithm evinces a 0.06-unit, 0.05-unit and 0.12-unit improvement in the PNSR, SSIM and NCC score, respectively, when compared to state-of-the-art.

(NNS) that based on Euclidean distance (ED), that is the square error of corresponding two patches, is one of widely used patches searching method, and the so-called BM3D algorithm [13] is one of its classical application. NNS method performs similar patches searching in a local window of image, thus is a local searching method. However, patches with significant structures, for example round edges or corners in an image, do not have a repeating pattern within a local window, which indicates that only using NNS to find similar patches is not an optimal method.
In order to better capture the structural similarity patches in natural images, various similar patches searching methods have been established in the literature. In [14], a socalled Needle patch representation method is proposed to improve on the reliability of patch-matching when the image quality deteriorates, and the patch representation Needle consists of small multi-scale versions of the patch and its immediate surrounding region. With the development of recent innovations in training deep convolutional neural network, a deep learning local image patches matching methods that based on Triplet and Siamese networks using a combination of triplet loss and global loss has been proposed in [15]. In [16], a 2-channel based neural network is established to learn a general image patches similarity function directly from raw image pixels.
Patch priors that give high likelihoods for patches will yield better patches restoration performance [17]. Recently, using learned specific probabilistic distribution of patches from noise-free image dataset to globally search and recover similar patches from a degraded image has become a hot research issue in image processing. In principle, the probabilistic distribution of patches can be modeled using arbitrary distributions, but most commonly, Gaussian mixture model (GMM) is often used to model patches priors [17], even though the literature [18] argues that a generalized Gaussian mixture model is a better fit for image patches prior modeling. In fact, patches priors are learned over small image patches by GMM will make computational tasks such as learning, inference and likelihood estimation much easier than working with whole images. GMM also has been shown to be a powerful tool for patches classification and for similar patches matching. In addition, the learned means, covariance matrices and mixing weights over all patches can greatly improve the accuracy for searching similar patches.
The expected patch log likelihood (EPLL) algorithm [17] based on GMM employs a global prior to search patch such that every selected patch is more like the given local prior. In [19], by integrating the external patch prior and internal self-similarity into one framework, the learned GMM (external patch prior) is used to guide the clustering of patches from the input degraded images, and then a patches matrix low-rank approximation process (internal self-similarity) is used to estimate the latent subspace for image recovery. In [20], a patch searching method that clusters similar patch candidates into patch groups using GMM-based clustering is proposed, and the selected patch groups that contain the reference patch are used to image denoising.
However, the simple GMM clustering method generally puts some nonsimilar patches, especially geometric nonsimilar patches, into a group which indicates this kind of similarity measure between patches is not accuracy. In other words, we need a more refined clustering method to determine which patches are similar or nonsimilar, which will motivate us to refine the GMM patches clustering method. To this end, in this paper, we exploit two additional operations to refine the GMM patches cluster: firstly, we employ a simple K-means methods to refine the patches classification by using the mean intensity of each patch; secondly, for each reference patch in a certain group, we gradually expand the search radius until the number of similar patches contained in this circle satisfy our needs. The aim to do this is because the more similar images patches are grouped together (in the patches matrix sense [19] or the patches tensor used in this paper), the better the patches denoising effect will be. To better exploit the spatial geometric structural information of similar patches, we apply the 3D patches tensor X , as shown in Fig. 1, to organize the similar patches instead of vectorizing each patches to form a similar patches matrix as in [19], then low-rank tensor approximation method based on HOSVD is used to perform image denoising. An iterative adaptive weighted core tensor thresholding algorithm is proposed to achieve low-rank tensor approximation. The reason using this algorithm are twofold, firstly, the tensor nuclear norm [21] (in fact is the ℓ 1 norm penalization) will result in significantly biased estimation to the values in core tensor, thus cannot achieve a reliable image recovery. Secondly, although the nonconvex penalty, such as ℓ p with (0 ≤ p < 1) , can also provide unbiased solution, it is computationally cumbersome. As a surrogate to nonconvex penalty [22], the weighted ℓ 1 norm not only can ameliorate the bias to the larger values in core tensor, but also is low computational complexity.
The initial motivation of this work come from two observations: firstly, tensor organization of image patches can preserve more structural information of image than matrix form of vectorized image patches exploited in the literature; secondly, the more accurate the distance/similarity measure between image patches is, the better the image denoising effect achieved by low-rank image patches tensor approximation method will be. These observations inspire us to propose an image similar patches clustering method, denoted as GKGLM, i.e., "GMM" + "k-means" + "inter-class geometry location similarity patches searching. " The advantages of these variants of GMM are the hybrids globally patches searching (i.e., GMM) and the locally patches searching (i.e., the k-means and inter-class geometry location similarity patches searching). Furthermore, consider the high computational complexity of training process of GMM using expectation maximization (EM) method, we propose a novel GPU based GMM training method using the back propagation and gradient descent method. Meanwhile, an iterative weighted image patches tensor approximation is proposed to perform image denoising, which the key advantage of this iterative method is the unbiased estimation to larger entry values in core tensor that are the main feature of images. See Fig. 2 for the flowchart of this paper. j = 1, · · · , J k = 1, · · · , K i = 1, · · · , I The main contributions of this paper arise from the following analysis: 1. A global similar patches searching method is proposed based on GMM patches clustering method, the k-means method and a local searching method within patches group are also being integrated into this framework to improve the accuracy of this patches clustering method, which will preserve more geometric structural similar information of patches. 2. A GPU based parallel training method is proposed to accelerate the GMM modeling training solved by EM algorithm. 3. In order to preserve more spatial geometric structural image information and fine details information image during denoising, an iterative adaptive weighted patches tensor approximation method is proposed.
The rest of the paper is organized as follows. In Sect. 2, we provide a brief introduction of the GMM and discuss how to search similar patches of the reference patch, including the algorithm of training GMM by the expectation maximization (EM) algorithm. In Sect. 3, we introduce the definition of three-order tensor, the properties of the corresponding core tensor and HOSVD algorithm. The iterative weighted image patches tensor approximation algorithm for image denoising is discussed in Sect. 4. We present numerical experiments and results in Sect. 5.The conclusion is presented in Sect. 6. 2 Similar patches matching by deep cluster based GMM training

GMM patches prior modeling
As being beforementioned, the probabilistic distribution of patches can be modeled using arbitrary distributions; but most commonly, GMM is used to learn patches priors. We take advantage of GMM as a probabilistic model to learn patch priors from noise-free patches set; then, noisy image patches clustering guided by the learned Gaussian mixture priors. GMM makes global searching for similar patches of a reference patch in the whole noise image. The probability of a random patch from clean images sets will be defined as follows: where Θ = (w 1 , . . . , w K , θ 1 , . . . , θ K ) is a set of parameters; w k denotes the weight that probability function p(x i |θ k ) ( k = 1, . . . K ) contributes to p(x i |Θ) and K k=1 w k = 1 ; θ k = (µ k , Σ k ) are mean value and covariance matrix, respectively, in density function p(x i |θ k ) . One of the most important feature of the patches is that it is close to a lowdimensional manifold for many natural images. Therefore, to more accurately measure the distance between different patches, we replace the frequently used Euclidean distance (ED) with Mahalanobis distance (MD) to measure the patches difference used in patches cluster problem, i.e., p( , where c is a constant. As shown the patches correct matching rates results in Table 1 from [19], even though the correct matching rates of ED is better than that of MD in image homogenous region, the MD overwhelmingly surpass that of ED in the structural and textural region, which indicates that using Mahalanobis distance to measure the difference between image patches is more reasonable than using Euclidean distance.

Parameters estimation using EM algorithm based on GPU
To accelerate image denoising algorithm, we exploit a pre-trained GMM for the first step to cluster the similar patches. That is the parameters in GMMs modeling in (1) will be learned from training image patches set using the expectation maximization (EM) algorithm. However, most traditional implementation of EM algorithm is high computational complexity, thus cost a large amount of time in training processing. On the other hand, because the most time-consuming process in training is the calculation of distance between different patches, thus, it is easy to consider to apply parallel computing to accelerate GMM training with GPU. Table 1 Patches correct matching rates (ED vs. MD) [19] Bold value indicates the best results combpared between MD and ED It is well known that an EM based GMM training methods include these two steps, i.e., expectation step and maximization step, respectively. Expectation step in EM algorithm is a kind of evaluation and cluster processing, that is the model parameters learned in the maximization step will be checked by evaluating the difference between the last and the current clustering results. And the maximization step is the parameters estimation processing, which is used to fit the new cluster results by using maximum likelihood estimation (MLE) method.
In the following, we establish a deep cluster based GMM training to accelerate the EM training algorithm. Particularly, in the parameters estimation step (maximization step), we first establish a binary cross-entropy (BCE) loss function based on the difference between the last clustering results Y k and the current clustering results Ŷ in the testing step (i.e., the Expectation step), then, the back-propagation algorithm and the stochastic gradient descent (SGD) method based on GPU are used to estimate the parameters of GMM modeling. See more details in Algorithm (1). We train the GMM prior model from a set of 2 million overlapping patches that randomly sampled from [23] with their DC removed. Allowing for overlapping when sampling patches is important because otherwise there would appear blocking artifact. See Fig. 3 for the GMM model learned flowchart.

Similar patches searching by deep cluster based GMM training
Based on the previously trained GMM patches model, given a reference patch, we propose a "GMM" + "k-means" + "inter-class geometry location similarity patches searching" method to refine the patches similarity matching. Given a noisy image Y, we firstly split it into a set of overlapping patches, then we denote P i Y by the i-th patch of Y. The likelihood of the i-th patch belongs to the kth ( k = 1, . . . , K ) cluster under parameter Θ is defined: where Σ k is a covariance matrix of the learned GMMs, σ is the noise standard deviation and I is an identity matrix. Maximizing the likelihood (6) to determine Gaussian component that makes the likelihood maximum, then we denote it to be the k th Gaussian component and assign the ith noisy patch into the k cluster. This is the GMM based noisy patches globally cluster, and all patches of the noisy image Y are split into different several clusters. As we have been emphasized in Sect. 2.1, there exists a deficiency of GMM based clustering method, i.e., in image homogenous region, the Euclidean distance based clustering method is better than the GMM based patches clustering method. Motivated by this observation, to improve the accuracy of patches matching, we employ a simple K-means method to each GMM based cluster to refine patches classification. After that, if there exists a class with the number of patches in it is less than 10, we will merge the patches in this class to the other classes most likely. Next, we explain how to choose number of similar patches of a reference patch in each class to form a patches tensor. We firstly determine which class each reference patch is in, then we select similar patches for each reference from this class by calculating the location distance from reference patch to each patch in this class, i.e., the distance from the upper-left pixel of reference patch and the patches in this class. Secondly, the first T least distance corresponding to the patches is selected to form a similar patches tensor. We call this process to be the inter-class geometry location similarity patches searching method.
. Thus, our similarity patches searching method can utilize the global information and the local location information of image. For convenience, we denote this clustering method as GKGLM, i.e., "GMM" + "k-means" + "inter-class geometry location similarity patches searching. " At last, we summarize the GKGLM algorithm in Algorithm 2.

Rank of patches tensor and HOSVD
Third-order patches tensors have column, row and tube fibers; Figs. 4 and 5 show the horizontal, lateral and frontal slides of a third-order. Low-rank matrix approximation methods are to measure the low-rank structure of the observed matrix X by minimizing its rank. Similarly, the latent tensor data can be approximated from low-rank version of observed tensor measurements. But, the low-rank decomposition to a multidimensional tensor usually is ticklish, and generally, there exist at least three different rank definitions of multidimensional tensor relative to different low-rank tensor decompositions, and it is hard to find a tight convex relaxation to the nonconvex rank function of multidimensional tensor. Two often used rank definitions of tensor are the CANDECOMP/PARAFAC (CP) rank and Tucker rank.

Definition 1
The CP rank of a tensor X ∈ R I 1 ×···×I N is defined as the minimum number of rank-1 decomposition, where V i = a (i) 1 ⊗ · · · ⊗ a (i) k , symbol ⊗ denotes vector outerproduct, a (i) j ∈ R n j is a vector and c i is the decomposition coefficient [24].

Definition 2
Tucker rank of a tensor X ∈ R I 1 ×···×I N is denoted by rank T (X ) and is a N-dimensional vector expressed as follows: where rank(X (i) ) denotes the matrix rank of ith-mode matricization(see the following section) X i to tensor X for i = 1, 2, . . . , N.
Based on these two definitions of tensor rank, tensor has a similar nuclear norm like the matrix nuclear norm, i.e., tensor nuclear norm (TNN), which is one extension of matrix nuclear norm [21].
The tensor trace norm (TTN) is similar to TNN, it is developed in [25] where α i > 0 and N i=1 α i = 1. TNN and TTN, in essence, merely is a convex combination of trace norms of the nth-mode unfolding X (i) . A novel tensor nuclear norm ||X || TNN is proposed in the literature [26] and is defined as the sum of the singular values of all frontal slices. This tensor nuclear norm has also been proved to be a norm and be the tightest convex relaxation to ℓ 1 norm of tensor multilinear rank.
A HOSVD to a tensor X ∈ R I 1 ×···×I N can be expressed as a core tensor S ∈ R r 1 ×···×r N and matrices multiplication shown as follows, where U (i) ∈ R I i ×r i (1 ≤ i ≤ R) is the left singular orthogonal matrix calculated by performing SVD on matrix X i . Figure 6 shows a HOSVD to a third-order tensor. The inverse HOSVD is defined by: where (·) T denotes the matrix transpose operation.
In some tensor decompositions, we need to unfold or flatten the tensor into matrix, a.k.s. matricization. There are N mode matricization to a Nth-order tensor. The nth-mode matricization of a tensor X ∈ R I 1 ×I 2 ×···×I N is denoted by X (n) ∈ R I n ×(I 1 I 2 ...I n−1 I n+1 ...I N ) (n = 1, 2, . . . , N ) is a matrix. We also denote functions unfold n (·) and fold n (·) to be the nth-mode unfolding and folding operation, that is, X (n) = unfold n (X ) and X = fold n (X (n) ).
There are many different arrangement ways of the nth-mode matricization, but in practice, permutation of elements doesn't affect the result of calculation [27]. See Fig. 7 for a third-order tensor matricization. The nth-mode product is denoted by X × n U where X ∈ R I 1 ×I 2 ×···×I N is a tensor and U ∈ R J ×I n is a matrix. We can obtain a tensor Y = X × n U of size I 1 × I 2 × · · · × I n−1 × J × I n+1 × · · · × I N . The nth-mode product can be expressed by matrix multiplication using fold n (·) as follows At last, the Frobenius norm of tensor X is denoted by �X � F = ( i 1,...,i N |x i 1,...,i N | 2 ) 1/2 ; 2-norm of vector v is denoted by v 2 . Fig. 6 HOSVD to third-order tensor

Low-rank patches tensor approximation
In this section, we establish a low-rank patches tensor approximation method that begin by an example. Consider a tensor X ∈ R 3×2×3 and the corresponding core tensor S with three horizontal slices, shown in Fig. 8. Three associated singular values matrices Σ (i) , (i = 1, 2, 3) obtained by SVD to the corresponding nth-mode matricization of tensor X , i.e., X (i) , (i = 1, 2, 3) , are listed as follows Three nth-mode matricization of core tensor S shown as follows, From this example, we can find that the quadratic sum of the ith row of S (n) equals to the square of the ith singular values of the corresponding Σ (n) , which proves one property, that is, S 2 F = X 2 F . For instance, 7.2051 2 = 51.9134 exactly equals (−7.1846) 2 + 0.0615 2 + 0.2235 2 + (−0.4906) 2 . Not only that, but all the Σ (i) 2 F , (i = 1, 2, 3) are identical, and equal to �S� 2 F = �X � 2 F , i.e., 92. That is, we have illustrated the following properties: denotes the lth singular value of the nth-mode matricization of tensor X. Based on these observations, we propose a novel low-rank patches tensor approximation method, i.e., by directly penalizing the adaptive weighted singular values of core tensor obtained by HOSVD to patches tensor.

Image denoising based on iterative low-rank patches tensors approximation algorithm
First, we introduce the objection function for image denoising in this paper.
where Y and X denote noise image and the latent clean image, respectively, W denotes adaptive weights, and ||̥(X)|| * denotes the proposed low-rank patches tensor nuclear norm (15). ̥(X) denotes tensor formulation and HOSVD to image X. Low rank is a common assumption to natural image in image processing, so we usually use ℓ 1 or nuclear norm to approximate rank function. ℓ 1 or nuclear norm regularization can lead to sparse solution, and their convexity is another important reason why they have become so popular in signal processing field. However, both of them will result in significantly biased estimates since the convexity, thus, cannot achieve a reliable solution. In comparison, as a nonconvex penalty, such as the ℓ q ( 0 ≤ q < 1 ), smoothly clipped absolute deviation (SCAD) or minimax concave (MC) penalty is superior to ℓ 1 or nuclear norm regularization because it can ameliorate the bias problem. The advantages of nonconvex penalty regularization have been verified in many applications, and in fact, a nonconvex penalty can yield significantly better performance than a convex penalty regularization. On the other hand, nonconvex regularization based sparse and low-rank recovery are of becoming considerable interest in recent years, also partly attributes to the recent progress in nonconvex and nonsmooth optimization algorithm theory. Motivated by these observations, we exploit iterative adaptive weights to the regularization term. We claim that such an iterative adaptive weighted scheme is equivalent to a kind of nonconvex penalty to the core tensors, where the singular values are assigned different penalty weights.
Proximity operator plays a central role in tackling nonconvex and nonsmooth inverse problem, and it is a highly efficient first-order algorithms which can scale well for highdimensional problems. For a proper and lower semi-continuous penalty function P (·) , where > 0 is a threshold parameter, consider the following scalar proximal projection: As P (·) is separable, the proximity operator of a vector t = [t 1 , . . . , t n ] ∈ R n , denoted by proxo (t) can be computed in an element-wise manner as: For commonly used proximity operators, refer to Table 2. For a patch tensor M , we denote a generalized penalty for low-rank promotion to its core tensor by P (·) , which is defined as : where i is the entry values in core tensor of tensor M.
Exploiting the adaptive weights ℓ 1 -norm we can better preserve essential characteristic of image, and the nonconvex objective function can be solved easily by using the soft-thresholding operator. For each GKGLM patches tensor X k , k = 1 . . . K , K is the number of reference patches, we have the following minimization where W k j is the adaptive weight, and k j is the value in core tensor S k of tensor X k . It is worth note that the adaptive weighted ℓ 1 -norm penalized optimization problem is nonconvex due to its varying weights W k j . But, if the weight W k j is nonincreasing assigned to the increasing core tensor absolute values | k j | , the penalized optimization problem is a convex function, the conclusion has been proved in the literature. In this paper, we are not to tackle the noise image directly, but to perform denoising for each tensor of being constituted by similar patches of reference patches. So, the objective function can be written as the following subproblems.
where X s y represents the kth reference patch tensor of noisy image Y, X s x is latent clear patches tensor corresponding to X s y , s j is the value of core tensor of noisy patches tensor X s y , and W s j is the adaptive weight and assigned to | s j | . It is obvious that the optimization problem can be solved easily by using the soft-thresholding operator. After achieving (17) proxp (t) = arg min Table 2 Regularization penalties and the corresponding proximity operator ( > 0)

Penalty name Penalty formulation Proximity operator
Soft-thresholding P (x) = |x| proxp (t) = sign(t) max{|t| − , 0} Hard thresholding X s x for each given noisy tensor X s y , we will put the new patches back into the image to constitute a new one, and an aggregation procedure is need. Until now, we finish the nth image denoising and obtain the clean image X n .

Adaptive weight setting to W s j
Adding residual (Y −X n ) of nth iteration back to the (n + 1) th step denoised image Y n+1 is the essence of the iterative algorithm. For a better denoising performance, we formulate the Y n+1 as follows: where n is the outer iteration number, k is the inner iteration number, and η is the relaxation parameter. Since, we want to add the residual back to the denoised image, the variance of noise remaining can be estimated by where σ is noise variance of Y, γ is a scaling factor used to control the re-estimation of noise variance, k × l is the number of all pixels in Y.
A advantages of adaptive-thresholding is to preserve large coefficients while filters small coefficients, so the features in image can be survived through thresholding. Then, we set the weight corresponding to each j to be: where N is the number of image patches of each tensor S s y , σ n can be computed by 23, ε > 0 is a sufficient small positive parameter for avoiding dividing by zero. We apply W s j as adaptive weight to solve the optimization problem 21 by the soft-thresholding proximal operator shown in Table 2, thus solution τ s j to the jth element of S s x can be expressed by as follows

Aggregation
When we put these patches back into the image, different patches tensor will be assigned different weights. Therefore, we use the weighted averaging to give more weight to patches tensors with less noise and less weight to those with much noise. Specifically, we use the following equation to define weight. where N is the number of patches in current tensor, p is the patch size and I is the number of thresholded elements in core tensor of this patches tensor. The result image can be calculated by the following formulation where N x,y are all GKGLM patches tensor overlapping in position (x, y), J (i) x,y are all patches in ith tensor that overlap in position (x, y) and Ω i,j is a pixel in the ith tensor, the jth tensor in position (x, y). The proposed algorithm is summarized in Algorithm 3, we denote the algorithm as GKGLM-Tensor.

Experiments
In this section, we conduct image denoising experiments on several widely used natural images, for example, House, Peppers, Barbara, Monarch and C.man. And we compare the proposed GMMH algorithm with several state-of-the-art denoising methods, including NNM [28], BM3D [13], EPLL [17], LSSC [12], NCSR [29], SAIST [30], WNNM [31], BM3D bst that is the BM3D based on boosted patch searching method [20] denoising algorithms, SLR [23] and LR-GSC [32]. As a common experimental setting in the literature, we add white Gaussian noise with zero mean and standard deviation σ to the noise-free images and to test the performance of competing denoising methods. The MATLAB source code of our proposed algorithm can be downloaded at https:// github. com/ Kazuk iAmak awa/ HOSVD_ denoi se_ patch.

Implementation details
Our presented image denoising method contains two stages: the searching image similar patches stage and the denoising stage. In the denoising stage, there exist four parameters need to be set: γ , β , the patch size p, the number M of patches in a tensor X . These parameters pick different values when the standard variation σ of noise is different (shown in Table 3). Intuitively, if you want to exploit more global image information of image, number M of patches in a tensor should be set as big as possible.

Results and discussion
We evaluate the competing methods from peak signal to noise ratio (PSNR), structural similarity index measurement (SSIM), normalized cross-correlation coefficient (NCC) and visual quality viewpoints.

PSNR
In  Table 4, it can be seen that our proposed algorithm achieves much better PSNR results than NNM. Secondly, our proposed algorithm has higher PSNR values than BM3D, LSSC, EPLL and NCSR, BM3D bst , SLR and is only slightly inferior to WNNM and LR-GSC for low noise level.  With the noise level increasing, the performance of our proposed algorithm will be improved significantly, which can be verified for σ = 50 noise level.

SSIM & NCC
In Table 5, we report the average SSIM and NCC results on three noise levels σ = 20, 40, 60 for image House, Peppers, Barbara, Monarch and C.man. We only compare our methods to BM3D, BM3D bst , EPLL, WNNM, SLR and LR-GSC, which outperforms the other algorithms. It can be seen from the SSIM viewpoint that our proposed method is comparable to these methods, and particularly, our method is superior to other method for σ = 50 noise level. The NCCNCC (normalized cross-correlation) shows the similarity of the denoised image with original clear image. The range of NCC [−1, 1] , and the closer the 1 is, the more similar the two image are, the closer the −1 indicates that the more dissimilar the two images. The NCC results in Table 5 show that there exists positive correlation between SSIM and NCC, which indicates that the higher SSIM value is, the bigger NCC value is.

Visual quality
Considering that human subjects are the ultimate judge of the image quality, the visual quality of denoised image is essential to evaluate a denoising algorithm. Figures 9 and 10 show the denoised images of Pepper and Monarch by the competing method. We obtain that WNNM has strong ability to image denoising by Table 4, so we just compare the visual effect of WNNM with our proposed algorithm. It can be seen that the proposed algorithm in this paper is likely to generate artifacts than WNNM, but our algorithm preserves edge more better. In summary, comparing to state-of-the-art denoising algorithm, our proposed method has performed a comparable job in image denoising, especially in high noise level.

Conclusion
In this paper, a new denoising approach based on image internal self-similarity prior, external patch priors and the low-rank patches tensor property of natural image is established. We search similar patches to each reference patch to form a tensor X (k) , and we decompose the tensor X (k) by HOSVD, then we use soft-threshold to tackle core tensor S (k) . In order to improve the performance of denoising, we have trained GMM by EM based on deep cluster. To enhance the accuracy for searching similar patches, a globally searching similar patches methods within the whole image is established and a geometric searching method is used to further improve the accuracy for finding similar patches.