Smoothing filters are a class of linear and nonlinear operators that gains significant attentions in image denoising recently. The formulations of these operators are simple: Consider a noisy observation y ∈ R n of a clean image z ∈ R n corrupted by additive i.i.d. Gaussian noise. A smoothing filter is a matrix W ∈ R n×n that generates a denoised estimate z as
where D def = diag {W 1} is a diagonal matrix for normalization so that each row of D -1 W sums to unity. The formulation in (1) is very general, and many denoising algorithms are smoothing filters, e.g., Gaussian filter [1], bilateral filter [2], non-local means (NLM) [3], locally adaptive regression kernel (LARK) [4], etc. Note that some of these filters are linear (e.g., Gaussian filter) whereas some are nonlinear (e.g. non-local means). There are interesting graph-theoretic interpretations of the smoothing filters [5]- [10], and there are also fast algorithms to compute the smoothing filters [11]- [16].
While smoothing filters work well for many denoising problems, it was observed in [17]- [19] that their performance can be further improved by modifying (1) as
where D c def = diag W T 1 is a diagonal matrix that normalizes the columns of W , and D r def = diag W D -1 c 1 is a diagonal matrix that normalizes the rows of W D -1 c . In other words, we modify (1) by introducing a column normalization step before applying the row normalization.
Before discussing the technical properties of (1) and ( 2), we first provide some numerical results to demonstrate an interesting phenomenon. In Figure 1, we crop the center 100 × 100 region of 10 standard clean images. We generate noisy observations by adding i.i.d. Gaussian noise of standard deviation σ = 20/255 to each clean image. These noisy images are then denoised by (1) and (2), respectively. The weight matrix W is chosen as the one defined in the non-local means (NLM) [3]. To ensure fair comparison, we choose the best parameter h r , the range parameter in NLM, for both methods. The patch size is set as 5 × 5 and the neighborhood search window is set as 21 × 21. The experiment is repeated for 20 independent Monte-Carlo trials to average out the randomness caused by different realizations of the i.i.d. Gaussian noise.
The results of this experiment are shown at the bottom of Figure 1. It is perhaps a surprise to see that (2), which is a simple modification of (1), improves the PSNR by more than 0.23 dB on average. Another puzzling observation is that if we repeatedly apply the column-row normalization, the PSNR does not always increase as more iterations are used. Figure 2 presents the result. In this experiment, we fix the NLM parameter h r to its optimal value when using the column-row normalization, i.e., (2). For 5 out of the 10 images we tested, the PSNR values actually drop after the first column-row normalization.
The above experiment piqued our curiosity and led us to a basic question: Why would the column-row normalization improve the denoising performance? Insights 0.71σ 0.77σ 0.77σ 0.77σ 0.71σ 0.71σ 0.79σ 0.77σ 0.77σ 0.79σ PSNR Improvement +0.10 +0.34 +0.23 +0.17 +0.29 +0.07 +0.33 +0.28 +0.05 +0.37 a better understanding of the underlying mechanism could potentially lead to a more systematic procedure that can generalize the operations in (2) and further improve the denoising performance. The goal of this paper is to address this issue and propose a new algorithm.
To the best of our knowledge, the above performance gain phenomenon was first discussed by Milanfar in [17], where it was shown that if we repeatedly apply the columnrow normalization we would obtain an iterative procedure called the Sinkhorn-Knopp balancing algorithm [20], [21] (or Sinkhorn-Knopp, for short.) As illustrated in Algorithm 1, Sinkhorn-Knopp is a simple algorithm that repeatedly applies the column and row normalization until the smoothing filter converges. For example, the smoothing filter defined by ( 2) is the result of applying Sinkhorn-Knopp for one iteration.
Sinkhorn-Knopp has many interesting properties. First,
when Sinkhorn-Knopp converges, the converging limit is a doubly-stochastic matrix -a symmetric non-negative matrix with unit column and row (also called a symmetric smoothing filter.) A doubly stochastic matrix has all of its eigenvalue's magnitudes bounded in [0, 1] so that repeated multiplications always attenuate the eigenvalues [22]. Moreover, the estimate z formed by a doubly stochastic matrix is admissible in the sense that no other estimates are uniformly better [23].
To explain the performance gain, Milanfar [17] considered a notion called "effective degrees of freedom", defined as
where z j is the jth pixel of the estimate and y j is the jth pixel of the input. df measures how an estimator trades bias against variance. Larger values of df imply a lower bias but higher variance. Milanfar argued that the overall mean squared error, which is the sum of the bias and the variance, is reduced because one can prove that symmetric smoothing filters have high effective degrees of freedom. However, effective degrees of freedom is not easy to interpret. It will be more useful if we can geometrically describe the actions to which the columnrow normalization are applying.
The present paper is motivated by our wish to further understand the mechanism behind the performance gain phenomenon. Our approach is to study an Expectation-Maximization (EM) algorithm for learning a Gaussian mixture model (GMM.) By analyzing the E-step and the M-step of the EM algorithm, we find that the actions of the symmetrization is a type of data clustering. This observation echoes with a number of recent work that shows ordering and grouping of non-local patches are key to high-quality image denoising [24]- [26]. There are two contributions of this paper, described as follows.
First, we generalize the symmetrization process by reformulating the denoising problem as a maximum-a-posteriori (MAP) estimation under a Gaussian mixture model. We show that the original smoothing filter in (1), the one-step Sinkhorn-Knopp in (2), and the full Sinkhorn-Knopp (i.e., iterative applications of Sinkhorn-Knopp until convergence) are all subroutines of the EM algorithm to learn the GMM. By showing that each method supersedes its preceding counterpart, we provide a possible explanation for the performance gain phenomenon.
Second, based on the analysis of the GMM, we propose a new denoising algorithm called the GMM symmetric smoothing filter (GSF). We show that GSF does not only subsume a number of smoothing filters, but also has a performance similar to some state-of-the-art denoising algorithms. We will discuss implementation and parameter selections for the GSF algorithm.
This paper is an extension of a conference article presented in [19]. In [19], the linkage between the GMM and the MAP denoising step was not thoroughly discussed. Specifically, the MAP was formulated as a weighted least squares step but the weights were indirectly obtained through a by-product of the EM algorithm. In this paper, we show that the quadratic cost function in the weighted least squares is indeed a surrogate function for solving the MAP problem. Therefore, minimizing the cost function of the weighted least squares is equivalent to minimizing an upper bound of the MAP objective function.
The rest of the paper is organized as follows. First, we provide a brief introduction to GMM and the EM algorithm in Section II. In Section III we discuss the generalization of different symmetrizations using the EM algorithm. The new GSF is discussed in Section IV and experimental results are shown in Section V. We conclude in Section VI.
Throughout this paper, we use n to denote the number of pixels in the noisy image, and k to denote the number of mixture components in the GMM model. To avoid ambiguity, we call a mixture component of a GMM as a cluster. Clusters are tracked using the running index i ∈ {1, . . . , k}, whereas patches (or pixels) are tracked using the running index j ∈ {1, . . . , n}. Without loss of generality, we assume that all pixel intensity values have been normalized to the range [0, 1].
We use bold letters to denote vectors and the symbol 1 to denote a constant vector of all ones. The vector x j ∈ R 2 represents the two-dimensional spatial coordinate of the jth pixel, and the vector y j ∈ R d represents a d-dimensional patch centered at the jth pixel of the noisy image y. For a clean image z, the jth patch is denoted by z j . A scalar y j ∈ R refers
to the intensity value of the center pixel of y j . Therefore, a d-dimensional patch y j (assume d is an odd number) is a vector y j = [y j-(d-1)/2 , . . . , y j , . . . , y j+(d-1)/2 ] T . To extract y j from the noisy image y, we use a linear operator P j ∈ R d×n such that y j = P j y. For some smoothing filters, the spatial coordinate x j is used together with a patch z j (or y j ). Therefore, for generality we define a generalized patch p j ∈ R p , which could be x j , z j (or y j ), or a concatenation of both:
The results presented in this paper are applicable to smoothing filters W taking the following form:
where κ i def = (2π) p |Σ i | is a normalization constant, and N (•) denotes a p-dimensional Gaussian:
with mean µ i ∈ R p and covariance matrix Σ i ∈ R p×p . We note that (3) is general and covers a number of widely used filters as shown in Table I.
Example 1 (Standard NLM): For the standard NLM [3], we have p j = y j . The ith mean vector is µ i = y i , i.e., the noisy patch is the mean vector. The ith covariance matrix is
where h r is the NLM parameter.
In practice, to reduce computational complexity, a search window Ω is often introduced so that neighboring patches are searched within Ω. This is equivalent to multiplying an indicator function to the weight as
where I {x ∈ Ω} = 1 if x ∈ Ω, and is zero otherwise.
Example 2 (Spatially Regulated NLM): As an alternative to the hard thresholding introduced by the indicator function in (6), one can also consider the spatially regulated NLM [17]. In this case, we can define p j = x T j , y T j T . The mean vector µ i and the covariance matrices will consist of two parts:
where h s and h r are the spatial and the range parameters, respectively. This leads to the weight
It is not difficult to see that the spatially regulated NLM coincides with the standard NLM as h s → ∞ and |Ω| → ∞.
In fact, the former uses a soft search window and the latter uses a hard search window. From our experience, we find that the spatially regulated NLM typically has better performance than the standard NLM when the best choices of h s and Ω are used in either case. Therefore, in the rest of this paper we will focus on the spatially regulated NLM.
The Gaussian mixture model (GMM) plays an important role in this paper. Consider a set of n generalized patches P def = {p 1 , . . . , p n } where p j ∈ R p . We say that P is generated from a Gaussian mixture model of k clusters if p j is sampled from the distribution
where π i ∈ R is the weight of the ith cluster, µ i ∈ R p is the mean vector, and Σ i ∈ R p×p is the covariance matrix. For the purpose of analyzing smoothing filters in this paper, we shall assume that the covariance matrices Σ i are fixed according to the underlying smoothing filter. For example, in spatially regulated NLM we fix the covariance matrices as Σ i = Σ NLM . When Σ i 's are fixed, we denote
Learning the model parameters Θ from P is typically done using the Expectation-Maximization (EM) algorithm [27]. The EM algorithm consists of two major steps: the expectation step (E-step) and the maximization step (M-step). The E-step is used to compute the conditional expected log-likelihood, often called the Q-function. The M-step is used to maximize the Q-function by seeking the optimal parameters Θ. The algorithm iterates until the log-likelihood converges. Since the EM algorithm is widely used, we skip the introduction and refer readers to [27] for a comprehensive tutorial. The EM algorithm for learning a GMM is summarized in Algorithm 2.
In this section, we discuss how various symmetric smoothing filters are generalized by the EM algorithm. We begin by discussing how the GMM can be used for denoising.
Algorithm 2 EM Algorithm for Learning a GMM with a known covariance matrix Σ [27].
, and the number of clusters k.
for i = 1, . . . , k, and set t = 0. while not converge do E-step: Compute, for i = 1, . . . , k and j = 1, . . . , n γ
Update Counter: t ← t + 1. end while
We first specify the denoising algorithm. Given the noisy image y, we formulate the denoising problem by using the maximum-a-posterior (MAP) approach:
where the first term specifies the data fidelity with λ as a parameter. The second term, which is a sum of overlapping patches (thus are dependent), is called the expected patch loglikelihood (EPLL) [25]. Note that EPLL is a general prior that uses the expectation of the log-likelihood of overlapping patches. It is not limited to a particular distribution for f (p j | Θ), although in [25] the GMM was used. Note also that the patch p j in ( 13) is extracted from the optimization variable z. Thus, by minimizing over z we also minimize over p j . Substituting ( 9) into ( 13), we obtain a GMM-based MAP denoising formulation
The computational challenge of ( 14) is the sum of exponentials inside the logarithm, which hinders closed-form solutions. In [25], Zoran and Weiss proposed to use a half quadratic splitting method to alternatingly minimize a sequence of subproblems, and select the mode of the GMM in each subproblem.
Our solution to handle the optimization problem in ( 14) is to use a surrogate function under the Majorization-Maximization framework [28]. The idea is to find a surrogate function
If we can find such function h(p j , p ′ j | Θ), then the minimization in ( 13) can be relaxed to minimizing an upper bound
For the GMM in ( 9), Zhang et al. [29] proposed one convenient surrogate function h(p j , p ′ j | Θ) whose expression is given in Lemma 1.
The function
where
is a surrogate function that satisfies ( 15) and ( 16).
Proof: See Appendix A of [29].
This variable can be chosen as p ′ j = [x T j , y T j ] T , i.e., the generalized patch using the noisy image y. Thus, γ ij is independent of p j .
Substituting the result of Lemma 1 into (17), we observe that (17) can be rewritten as
where we dropped terms involving p ′ j as they are independent to the optimization variable z. Note that ( 20) is a quadratic minimization, which is significantly easier than (14).
For the case of spatially regulated NLM, it is possible to further simplify (20) by recognizing that p j involves both spatial coordinate x j and patch z j . Therefore, by expressing p j = [x T j , z T j ] T and by defining µ i using ( 7), (20) becomes
.
In this minimization, we observe that since x j is not an optimization variable, it can be eliminated without changing the objective function. By using a patch extract operator P j , i.e., z j = P j z, we obtain the minimization
where we have absorbed the NLM parameter 2h 2 r (the diagonal term of Σ) into λ. Problem (P 1 ) is the main optimization of interest in this paper. We call (P 1 ) the GMM Symmetric Smoothing Filter (GSF). In the literature, there are other GMM based denoising algorithms. Their connections to GSF will be discussed in Section IV .
We now discuss the role of symmetrization by studying (P 1 ) and the EM algorithm for learning the GMM. To keep track of the iterations of the EM algorithm, we use the running index t denotes the iteration number. For consistency, we will focus on the spatially regulated NLM.
In spatially regulated NLM, the ith pixel of the denoised image is
with W ij defined in (8). To obtain ( 21) from (P 1 ), we consider the following choices of parameters
In this case, since the quadratic objective in (P 1 ) is separable, the ith term becomes
which coincides with (21) as
It is important to study the conditions in (22). First, the patch extractor P j extracts a pixel z j from z. This step is necessary because NLM is a weighted average of individual pixels, even though the weights are calculated using patches. Accordingly, the mean vector µ i is also a pixel y i so that it matches with the optimization variable z i . From a clustering perspective, we can interpret µ i = y i as having one cluster center for one pixel, i.e., every pixel has its own cluster. Clearly, this is a suboptimal configuration, and we will address it when we present the proposed method. We also note that in (22), the parameter λ is 0. What it means is that the data fidelity term is not used and only the EPLL prior is required to obtain the NLM result.
The most interesting observation in (22) is the choice of γ ij . In fact, the γ ij in ( 22) is only the numerator of ( 19) by assuming π i = 1/k. If we let Θ i = (π i , µ i ) be the ith model parameter of a GMM, then the physical meaning of ( 22) is a conditional probability of observing p j given that we pick the ith cluster, i.e., P(p j | Θ i ). In contrast, ( 19) is a posterior probability of picking the ith cluster given that we observe p j , i.e., P(Θ i | p j ). We will discuss this subtle difference more carefully in the next subsection when we discuss the one-step Sinkhorn-Knopp.
In one-step Sinkhorn-Knopp, we recognize that the ith pixel of the denoised image is
This result can be obtained from (P 1 ) by letting
Unlike the conditions posted by the original smoothing filter in (22), one-step Sinkhorn-Knopp defines γ ij according to (19), or with some substitutions that yields (25).
As mentioned in the previous subsection, the γ ij defined in ( 22) is a conditional probability whereas that in ( 25) is a posterior probability. The posterior probability, if we refer to the EM algorithm (see Algorithm 2), is in fact the E-step with initializations π
For the original filter, the E-step is not executed because γ ij is defined through (22). Therefore, from a clustering perspective, it is reasonable to expect a better denoising result from one-step Sinkhorn-Knopp than the original smoothing filter because the clustering is better performed.
To further investigate the difference between the conditional probability P(p j | Θ i ) and the posterior probability P(Θ i | p j ), we adopt a graph perspective by treating each patch p j as a node, and γ ij as the weight on the edge linking node i and node j [5]. In the original smoothing filter, the conditional probability P(p j | Θ i ) causes a "majority vote" effect, meaning that the parameter Θ i has a direct influence to every patch {p j } in the image. Therefore, if Θ i has many "weak friends", the sum of these "weak friends" can possibly alter the denoising result which would have been better obtained from a few "good friends".
In contrast, the one-step Sinkhorn-Knopp uses the posterior probability P(Θ i | p j ). From Bayes rule, the posterior probability is related to the conditional probability by
Since P(Θ i ) = π i and π i = 1/k, we see that P(Θ i | p j ) is the ratio of P(p j | Θ i ) and P(p j ). P(p j ) measures the popularity of p j . Thus, if p j is a popular patch (i.e., it is a "friend" of many), then the normalization P(p j | Θ i )/P(p j ) is a way to balance out the influence of p j . Interpreted in another way, it is equivalent to say that if Θ i has many "weak friends", the influence of these "weak friends" should be reduced. Such intuition is coherent to many NLM methods that attempt to limit the number of nearby patches, e.g., [30], [31].
The definite answer to whether there is performance gain due to one-step Sinkhorn-Knopp is determined by the likelihood of obtaining "weak friends". This, in turn, is determined by the image content and the NLM parameter h r which controls the easiness of claiming "friends". For example, if an image contains many textures of a variety of content and if h r is large, then it is quite likely to obtain "weak friends". In this case, one-step Sinkhorn-Knopp will improve the denoising performance. On the other hand, if an image contains only a constant foreground and a constant background, then most patches are "good friends" already. Applying the one-step Sinkhorn-Knopp could possibly hurt the denoising performance. Figure 3 shows an example. To further justify the above claim that the performance gain is caused by the likelihood of obtaining "weak friends", we consider the 10 images in Figure 1 by plotting the PSNR gain as a function of h r . Our hypothesis is that for small h r , the performance gain should be small or even negative because it is difficult for a patch to find its "friends". When h r is large, the performance gain should become significant because many "weak friends" will be balanced out by the one-step Sinkhorn-Knopp. As shown in Figure 4, this is in fact the case: For h r lies between 0 and certain threshold (around 0.65σ where σ is the noise standard deviation), PSNR gain is always zero or negative. When h r increases, PSNR gain becomes positive. The result is consistent for all 10 images we tested.
The full Sinkhorn-Knopp (Algorithm 1) is an iterative algorithm that repeatedly applies the one-step Sinkhorn-Knopp until convergence. To analyze the full Sinkhorn-Knopp algorithm, we first recognize that the algorithm can be effectively described by two steps:
where the first equation is a row normalization and the second equation is a column normalization. This pair of normalization can be linked to the EM algorithm in the following sense. First, in the EM algorithm we fix the mixture weight π (t) i as π
(t) i = 1/k for all clusters i = 1, . . . , k and all iteration numbers t = 1, . . . , t max . This step is necessary because the full Sinkhorn-Knopp does not involve mixture weights. Intuitively, setting π (t) i = 1/k ensures that all clusters have equal probability to be selected.
When
is fixed, the EM algorithm in Algorithm 2 has only two steps: Update of γ (t) ij in (10) and update of µ (t) i in (12). Inspecting this pair of equations, we observe that (10) appears a column normalization whereas (12) appears a row normalization. However, since full Sinkhorn-Knopp does not have a mean vector µ (t) i , we have to modify the EM algorithm in order to link the two. To this end, we modify ( 12) by defining a sequence β
and modify (10) by updating of γ
Under this setting, ( 27)-( 28) becomes exactly (26).
It is important to understand the difference between the original EM algorithm using (10)-( 12) and the modified EM algorithm using ( 27)- (28). In the M-step of the modified EM algorithm, the mean vector µ (t) i is absent. However, µ (t) i is still updated, though implicitly, because
where (a) follows from (27). Therefore, β
ij are the coefficients that form µ (t) i through a linear combination of p j 's.
In the E-step of the modified EM algorithm, since µ (t) i is absent, one cannot compute the conditional probability N (p j | µ (t) i , Σ) and hence the posterior probability in (10). To resolve this issue, we recognize that since
ij is the coefficient for the mean vector as shown in (29) The fact that the full Sinkhorn-Knopp does not resemble a complete EM algorithm offers some insights into the perfor-mance gain phenomenon. Recall from the results in Figure 2, we observe that the first Sinkhorn-Knopp iteration always increases the PSNR except the artificial images we discussed in Section III.C. This can be interpreted as that the clustering is properly performed by the EM algorithm. However, as more Sinkhorn-Knopp iterations are performed, some images show reduction in PSNR, e.g., images 3, 4, 8, 10. A close look at these images suggests that they contain complex texture regions that are difficult to form few but distinctive clusters. In this case, the approximation of
ij is weak and hence the denoising performance drops.
To summarize our findings, we observe that the performance of the normalization is related to how the EM algorithm is being implemented. A summary of these findings is shown in Table II. For all the three algorithms we considered: the original filter, the one-step Sinkhorn-Knopp, and the full Sinkhorn-Knopp algorithm, the EM algorithm is not completely performed or in-properly configured. For example, setting k = n causes excessive number of clusters and should be modified to k < n; the MAP parameter λ is always 0 and should be changed to a positive value to utilize the data fidelity term in (P 1 ); the E-step and the M-step are not performed as it should be. Therefore, in the following section we will propose a new algorithm that completely utilizes the EM steps for problem (P 1 ).
The proposed denoising algorithm is called the Gaussian Mixture Model Symmetric Smoothing Filter (GSF). The overall algorithm of GSF consists of two steps:
• Step 1: Estimate the GMM parameter µ (r) i and γ ij from the noisy image the by EM algorithm.
• Step 2: Solve Problem (P 1 ), which has a closed form solution. In the followings we discuss how GSF is implemented.
First of all, we recall that since (P 1 ) is a quadratic minimization, it is possible to derive a closed form solution by considering the first order optimality condition, which yields a normal equation
where the vector w j is defined as
Equations ( 30)-( 31) has a simple interpretation: The intermediate vector w j is a weighted average of the mean vectors {µ (r) i } k i=1 . These {w j } n j=1 represent a collection of (denoised) overlapping patches. The operation P T j on the right hand side of (30) aggregates these overlapping patches, similar to the TABLE II: Generalization and comparisons using EM algorithm for learning GMM with a known Σ.
One-step Full GSF Filter [3] Sinkhorn-Knopp [18] Sinkhorn-Knopp [21] (Proposed) No. Clusters
aggregation step in BM3D [32]. The addition of λy regulates the final estimate by adding a small amount of fine features, depending on the magnitude of λ.
In order to use ( 30)-( 31), we must resolve two technical issues related to the EM algorithm and Problem (P 1 ): (i) How to determine λ; (ii) How to determine k.
Ideally, λ should be chosen as the one that minimizes the mean squared error (MSE) of the denoised image. However, in the absence of the ground truth, MSE cannot be calculated directly. To alleviate this difficulty, we consider the Stein's Unbiased Risk Estimator (SURE) [33], [34]. SURE is a consistent and unbiased estimator of the MSE. That is, SURE converges to the true MSE as the number of observations grows. Therefore, when there are sufficient number of observed pixels (which is typically true for images), minimizing the SURE is equivalent to minimizing the true MSE.
In order to derive SURE for our problem, we make an assumption about the boundary effect of P j .
Assumption 1: We assume that the patch-extract operator {P j } n j=1 satisfies the following approximation:
We note Assumption 1 only affects the boundary pixels and not the interior pixels. Intuitively, what Assumption 1 does is to require that the boundary pixels of the image are periodically padded instead of zero-padded. In the image restoration literature, periodic boundary padding is common when analyzing deblurring methods, e.g., [35].
Under Assumption 1, we can substitute (32) into (30) and take the matrix inverse. This would yield
where Then, we can derive the SURE of z as follows.
Proposition 1: Under Assumption 1, the SURE of z(λ) is
where σ 2 def = 1 n uy 2 , and
where e j ∈ R d is the jth standard basis. Proof: See Appendix B. The SURE given in ( 35) is a one-dimensional function in λ. The minimizer can be determined in closed-form.
Corollary 1: The optimal λ that minimizes SURE(λ) is
Proof: (35) is a differentiable function in λ. Therefore, the minimizer can be determined by considering the first order optimality and set the derivative of SURE(λ) to zero. The projection operator max(•, 0) is placed to ensure that λ * ≥ 0.
Example 3: To demonstrate the effectiveness of SURE, we show a typical MSE and a typical SURE curve of a denoising problem. In this example, we consider a 128 × 128 image "Baboon", with noise standard deviation of σ = 30/255. The non-local means parameters are h r = σ and h s = 10. The number of clusters is k = 50, and the patch size is 5 × 5. The results are shown in Figure 5, where we observe that the SURE curve and the true MSE curve are very similar. In fact, the minimizer of the true MSE is λ = 8.0080 with a PSNR of 24.5143dB whereas the minimizer of SURE is λ = 7.9145 with a PSNR of 24.5141dB.
Remark 2: Careful readers may notice that in (36), we implicitly assume that γ ij is independent of y j . This implicit assumption is generally not valid if γ ij is learned from y. However, in practice, we find that if we feed the EM algorithm with some initial estimate (e.g., by running the algorithm with λ = 0), then the dependence of γ ij from y j becomes negligible.
The number of clusters k is another important parameter. We estimate k based on the concept of cross validation [36].
Our proposed cross-validation method is based on comparing the estimated covariance with Σ NLM . More specifically, we compute the estimated covariance
where
i ] T is the mean returned by the EM algorithm, and
is the converged weight. Then, we compute the ratio of the deviation
Ideally, if Σ i = Σ NLM , then by (39) we have δ i (k) = 1. However, if the ith estimated Gaussian component has a radius significantly larger than h r (or, h s for the spatial components), then the covariance Σ i would deviate from Σ NLM and hence δ i (k) > 1. Conversely, if the ith estimated Gaussian component has a radius significantly smaller than h r , then we will have δ i (k) < 1. Therefore, the goal of the cross validation is to find a k such that δ i (k) is close to 1.
To complete the cross-validation setup, we average δ i (k) over all k clusters to obtain an averaged ratio
The parenthesis (k) in (40) emphasizes that both δ(k) and δ i (k) are functions of k. With (40), we seek the root k of the equation δ(k) = 1.
The root finding process for δ(k) = 1 can be performed using the secant method. Secant method is an extension of
Fig. 6: Illustration of the secant method. Given k (a) and k (b) , we compute k (c) according to the slope defined by the line linking δ (a) and δ (b) .
Input: k (a) and k (b) such that δ (a) > 1 and δ (b) < 1.
. end if end while the bisection method in which the bisection step size (i.e., 1/2) is now replaced by an adaptive step size determined by the local derivative of the function. Let k (a) and k (b) be two number of clusters, and δ (a) and δ (b) be the corresponding cross-validation scores, i.e., δ (a) = δ(k (a) ). If δ (a) > 1 and δ (b) < 1, the secant method computes the new k as
If δ(k (c) ) > 1, then we replace k (a) by k (c) ; Otherwise, we replace k (b) by k (c) . The process repeats until the |k (a) -
A pictorial illustration of the secant method is shown in Figure 6. A pseudo code is given in Algorithm 3.
Example 4: To verify the effectiveness of the proposed cross validation scheme, we consider a 128 × 128 "House" image with noise σ = 60/255. The patch size is 5×5, h r = σ, and h s = 10. Figure 7 shows the PSNR value of the denoised image and the corresponding cross validation score δ(k) as a function of k. For this experiment, the maximum PSNR is achieved at k = 144, where PSNR = 26.0257dB. Using the cross-validation score δ(k), we find that δ(k) is closest to 1 when k = 130. The corresponding PSNR value is 25.9896dB, which is very similar to the true maximum PSNR.
In this section, we present additional simulation results to evaluate the proposed GSF.
We consider 10 testing images, each of which is resized to 128 × 128 (so n = 16384) for computational efficiency. The noise standard deviations are set as σ ∈ {20/255, 40/255, 60/255, 80/255, 100/255}. Several existing denoising algorithms are studied, namely the NLM [3], One-step Sinkhorn-Knopp [18], BM3D [32], EPLL [25], Global image denoising (GLIDE) [12], NL-Bayes [37], and PLE [38]. The parameters of the methods are configured as shown in Table III.
For NLM and One-step Sinkhorn-Knopp (One-step, in short), we use the spatially regulated version due to its better performance over the standard NLM. We implement the algorithms by setting the patch size as 5×5 (i.e., d = 25). The parameters are h s = 10 and h r = σ √ d. The full Sinkhorn-Knopp algorithm is implemented using GLIDE [12], where the source code is downloaded from the author's website 1 . Default settings of GLIDE are used in our experiment.
For the proposed GSF, we keep the same settings as NLM except for the intensity parameter h r where we set h r = σ. The omission of the factor √ d is due to the fact that each Gaussian component is already a d-dimensional multivariate distribution. It is therefore not necessary to normalize the distance y iy j 2 by the factor d. For BM3D, EPLL, NL-Bayes, we downloaded the original source code from the author's website 2,3,4 . For PLE, we modified an inpainting version of the source code provided by the authors. Default settings of these algorithms are used.
Among these methods, we note that EPLL is an external denoising algorithm where a Gaussian mixture is learned from a collection of 2 million clean patches. All other methods (including GSF) are single image denoising algorithms.
Configuration NLM [3] Patch size 5 × 5, h s = 10, h r = σ √ d One-step [18] Patch size 5 × 5, h s = 10, h r = σ √ d GSF (Ours) Patch size 5 × 5, h s = 10, h r = σ GLIDE [12] Default settings. Pilot estimate uses NLM. NL-Bayes [37] Base mode. Default settings. PLE [38] Default settings. Default initializations. BM3D [32] Default settings. EPLL [25] Default settings. External Database.
The overall results of the experiment are shown in Table VI. We first compare the PSNR values of GSF with NLM, Onestep and full Sinkhorn-Knopp.
In Table IV we show the average PSNR over the 10 testing images. In this table, we observe that on average One-step has a higher PSNR than NLM by 0.12dB to 1.12dB, with more significant improvements at low noise levels. This implies that the "grouping" action by the column normalization becomes less influential when noise increases. Moreover, if we compare GSF with NLM and One-step, we observe that the PSNR gain is even larger. Even at a high noise level (e.g., σ = 80/255 or σ = 100/255), the average gain from NLM is 2.5dB or more. Besides studying the trend of PSNR as a function of σ, it is also interesting to compare the PSNR when we increase the spatial parameter h s . In Table V, we show the PSNR improvement when we use different h s ∈ {5, 10, 20, 50, 100} for a 128 × 128 image. The results show that when h s increases, the PSNR improvement also increases. One reason is that in (7), the spatial parameter h s controls the diagonal bandwidth of the smoothing filter W . That is, a small h s leads to a banded diagonal W with small bandwidth. In the limit when h s → 0, W will become a diagonal matrix, and hence is immune to any column normalization. Therefore, the effectiveness of the column normalization in the One-step depends on how large h s is.
The full Sinkhorn-Knopp algorithm is implemented using GLIDE [12]. GLIDE consists of multiple steps: It first determines the weight matrix, followed by a full Sinkhorn-Knopp algorithm that symmetrizes the weight matrix. Then, it incorporates an estimator to optimally determine the number of non-zero eigenvalues and the power of eigenvalues of the smoothing filter. GLIDE can use any denoising result as its pilot estimate. For the fairness of the experiment we follow the default setting of GLIDE and use the standard NLM as the
Since the proposed GSF uses a Gaussian mixture model, we compare it with three other algorithms that also use Gaussian mixture models. These algorithms are the NL-Bayes [37], the piecewise linear estimator (PLE) [38] and the EPLL [25].
Methodologically, there are some important differences between GSF, NL-Bayes, PLE and EPLL. NL-Bayes has a grouping procedure that groups similar patches, a process similar to BM3D. The grouped patches are used to estimate the empirical mean and covariance of a single Gaussian. In GSF, there is no grouping. The other difference is that the denoising of NL-Bayes is performed by a conditional expectation over the single Gaussian. In GSF, the denoising is performed over all clusters of the GMM. Experimentally, we observe that GSF and NL-Bayes have similar performance, with NL-Bayes better in the low noise conditions and GSF better in the high noise conditions. One possible reason is that the grouping of NL-Bayes uses the standard Euclidean distance as a metric, which is not robust to noise.
PLE first estimates the model parameters using the noisy image. Then, for every patch, the algorithm selects a single Gaussian. The denoising is performed by solving a quadratic minimization and is performed for each patch individually. The algorithm iterates by improving the model parameters and the denoised estimate until convergence. GSF does not have this iterative procedure. Once the GMM is learned, the denoising is performed in closed-form. The other difference is that PLE requires a good initialization and is very sensitive to the initialization. Experimentally, we find that GSF performs better than PLE using a MATLAB code provided the authors of PLE. In this MATLAB code, the initialization was originally designed for image inpainting at a particular image resolution. Because of the sensitivity of PLE to the initializations, its performance on denoising is not particularly good. With a better initialization, we believe that PLE would improve. However, even so the gap between GSF and PLE will unlikely be significant because PLE performs worse that BM3D and EPLL.
The closest comparison to GSF is EPLL as both algorithms solve a whole-image MAP minimization with a Gaussian mixture model. To evaluate the difference between the two algorithms, we consider an experiment by feeding the noisy patches the EM algorithm to learn a GMM. The patch size is fixed at 5 × 5, and the number of clusters is fixed as k = 100. We repeat the experiment by inputting the denoised result of BM3D and the oracle clean image into the EM algorithm.
From Table VIII, we observe that EPLL with a noisy input performs poorly. The reason is that the original EPLL trains the GMM from 2 million clean patches. When feeded with noisy images, the GMM trained becomes a non-informative prior distribution. Moreover, in EPLL the GMM involves (π i , µ i , Σ i ) whereas in GSF the GMM only involves (π i , µ i ). This is a significant reduction in the number of model parameters. When feeded with only a single image, there is insufficient training sample for EPLL.
Another observation from Table VIII is that the performance of EPLL depends heavily on the quality of the GMM. For example, if we use the result of BM3D as a pilot estimate for learning the GMM, the performance of EPLL is similar to the oracle case where we use the clean image. However, using BM3D as a pilot estimate is not a plausible approach because by running BM3D alone we can get an even higher PSNR (See Table VI). This result further shows the effectiveness of the proposed GSF for single image denoising.
Finally, we discuss the complexity and limitations of the proposed GSF. GSF is a one-step denoising algorithm when γ ij and µ (r) i are known. However, learning the GMM using the EM algorithm is time-consuming, and the complexity depends on the number of clusters k. In addition, since k needs to be estimated through a cross-validation scheme, the actual complexity also depends on the number of cross-validation steps. To provide readers an idea of how k changes with other system parameters, we conduct two experiments.
In Table IX we show the number of clusters returned by the cross-validation scheme as we increase the noise level. As shown, the number of clusters increases when noise level reduces. This result is consistent with our intuition: As noise reduces, the grouping of patches becomes less important. In the limit when the image is noise-free, every patch will become its own cluster center. Therefore, one limitation of GSF is that for low-noise images the computing time could be very long. However, GSF is still a useful tool as its simple structure offers new insights to denoising. Now, if we fix the noise level but change the image size, the complexity of GSF also varies. In Table X, we show the number of clusters as a function of image size. As a reference we also show the PSNR values of GSF and that of BM3D. The result in Table X indicates that the number of clusters increases with the image size. In the table, we also observe that BM3D performs worse than GSF for small images, but becomes better as image size increases.
Remark 3 (Subtraction of mean): It is perhaps interesting to ask whether it is possible to subtract the mean before learning the GMM, as it could reduce the number of clusters. However, from our experience, we find that this actually degrades the denoising performance. If the GMM is learned from a collection of zero-mean patches, the denoising step in (P 1 ) can only be used to denoise zero-mean patches. The mean values, which are also noisy, are never denoised. This phenomenon does not appear in EPLL (in which the GMM has a zero-mean) because the means are iteratively updated. We followed the same approach to iteratively update the means. However, we find that in general the denoising performance is still worse than the original GMM with means included. Further exploration on this would likely provide more insights into the complexity reduction issue.
Motivated by the performance gain due to a column normalization step in defining the smoothing filters, we study the origin of the symmetrization process. Previous studies have shown that the symmetrization process is related to the Sinkhorn-Knopp balancing algorithm. In this paper, we further showed that the symmetrization is equivalent to an EM algorithm of learning a Gaussian mixture model (GMM). This observation allows us to generalize various symmetric smoothing filters including the Non-Local Means (NLM), the one-step Sinkhorn-Knopp and the full Sinkhorn-Knopp, and allows us to geometrically interpret the performance gain phenomenon.
Based on our findings, we proposed a new denoising algorithm called the Gaussian mixture model symmetric smoothing filters (GSF). GSF is a simple modification of the denoising framework by using the GMM prior for the maximum-aposteriori estimation. Equipped with a cross-validation scheme which can automatically determine the number of clusters, GSF shows consistently better denoising results than NLM, One-step Sinkhorn-Knopp and full Sinkhorn-Knopp. While GSF has slightly worse performance than state-of-the-art methods such as BM3D, its simple structure highlights the importance of clustering in image denoising, which seems to be a plausible direction for future research.
We thank Dr. Guosheng Yu and Prof. Guillermo Sapiro for sharing the code of PLE [38]. We also thank the anonymous reviewers for the constructive feedback that helps to significantly improve the paper.
Given an estimator z of some observation y, the SURE is defined as
Substituting ( 33) into (42) yields
where σ 2 def = 1 n uy 2 . So it remains to determine div( z).
GLIDE: https://users.soe.ucsc.edu/ ∼ htalebi/GLIDE.php