1. Introduction

Prior models play a key role in the area of wavelet-based Bayesian image denoising [1]. Over the last two decades, a number of models have been developed to characterize the statistical behaviors of natural image wavelet coefficients. The widely-used univariate statistical models include the Generalized Laplacian model [2, 3], Alpha-stable model [4, 5], exponential model [6], Bessel K Form density model [7, 8, 9] and mixture of two distributions [10, 11]. These univariate models are suitable to characterize an empirical statistic of natural images: the filter response of natural images is always sharply peaked at zero with heavy tails [12]. Recently, in order to describe the joint statistical behavior among natural image wavelet coefficients, some multivariate models also have been proposed such as the multivariate Laplacian model [13, 14, 15, 16], multivariate exponential model [17], multivariate Alpha-stable model [18], multivariate Bessel K Form density model [7], multivariate scale mixture of Gaussian (GSM) model [19], multivariate normal inverse Gaussian model [20], and multivariate Elliptically Contoured Distribution (ECD) model [21].

Based on various prior models such as those mentioned above, many Bayesian image denoising algorithms have been developed [2, 4, 5, 9, 10, 11, 13, 14, 16, 17, 18, 19, 20, 21, 22]. Due to the fact that these algorithms use different prior models, as well as different Bayesian inference technologies (such as the Minimum Mean-squared Error (MMSE) estimate or the Maximum a posterior (MAP) estimate), they always present very different forms. Interestingly, although most existing natural image wavelet coefficient models have obviously different analytical forms, they can be categorized in the same class - they essentially belong to the GSM models [23]. Thus, a natural question is: can we develop a unified denoising framework for these models? Such a unified framework, if available, will be very useful not only for the image denoising application itself but also in the research of natural image statistic. For example, using the unified framework, one can fairly compare the performance of these models for image denoising, thus understands which model can characterize the statistical behaviors of natural image wavelet coefficients the best and why. In addition, when a new GSM-type prior model is available, to use it for image denoising, one only needs to simply put it into the unified denoising framework rather than to develop a new estimator based on its specific analytical form.

In this paper, we will develop a novel image denoising algorithm by combining the EM scheme and the properties of the Scale Mixture of Gaussian (GSM) models. The developed algorithm provides a unified image denoising framework for the GSMmodels that have an analytical expression form. The algorithm is of a simple iterative form and only uses the derivative information of a prior model. When a new prior model is available, it can be easily adopted for image denoising with the developed algorithm. As examples, several classical and recently-proposed statistical models of wavelet coefficients of natural images are tested in the denoising framework and some new results are obtained.

The outline of the rest of the paper is as follows. Section 2 gives the derivation of the unified image denoising algorithm. Examples are given that test several classical and recently-proposed statistical models in Section 3. Section 4 reports experimental results and provides a sufficient discussion. Concluding remarks are summarized in Section 5.

2. Unified denoising algorithm: derivation

Assume that noisy image wavelet coefficients have been partitioned into M small neighborhoods, and each neighborhood includes d coefficients. Let Y={Y ₁,…,Y_M} denotes the M observed d-dimensional noisy data and each Y_i, 1≤i≤M, satisfies the model

where X_i, 1≤i≤M, is the corresponding noise-free wavelet coefficients and ε~N(0⃗, ρ_ε) is the additive Gaussian white noise. Suppose that the noise-free wavelet coefficient X_i, 1≤i≤M, are i.i.d random variables each having a GSM distribution $X_i\stackrel{d}{=}\sqrt{z_i}U,$ , where U~N(0⃗, ρ) and z_i is a hidden positive scalar multiplier. Let Z={z ₁, …, z_M} and z_i, 1≤i≤M, be statistically independent. Thus, we have P(Z)= $P\left(Z\right)=\underset{i=1}{\overset{M}{\Pi }}P\left(z_i\right)$ .

In [24], J. M. Bioucas-Dias proposed a highly effective image deconvolution algorithm based on the univariate GSM models. The algorithm employed a new generalized Expectation Maximization (EM) scheme, where the missing variables were the scalar multipliers of the GSM models and the maximization step of the EM scheme was replaced with a linear stationary second-order iterative method. We will use the standard EM scheme to estimate X={X ₁, …,X_M} under the observed data Y={Y ₁, …,Y_M} in this paper. Similarly to the technology adopted in [24], we will set the missing variables to be the scalar multipliers of the GSM models. As we will show below, combined with several special properties of the GSM models, the EM scheme will lead to a very simple iterative image denoising algorithm that is suitable for both univariate and multivariate GSM models.

The standard EM scheme works iteratively by alternating two steps: the E-step (expectation step) and the M-step (maximization step). In our setting, Z is the missing (or hidden) variable and (Y,Z) is the so-called complete data. The E-step finds the expected value of the loglikelihood logP(X,Y,Z) given the observed data Y and the current estimate X ^t={X^t ₁, …,X^t_M},

Using the Bayesian rule, we have P(X,Y,Z)=P(Y|Z,X)P(Z,X)=P(Y|X)P(X|Z)P(Z). As a result, Eq. (2) can be rewritten as

where C(Z,Y,X ^t) is a function of Z,Y and X ^t. In addition, since

$P\left(Z|\mathbf{Y},{\mathbf{X}}^t\right)=\frac{P(\mathbf{Y}|{\mathbf{X}}^t,Z)P\left({\mathbf{X}}^t|Z\right)P\left(Z\right)}{P\left(\mathbf{Y}|{\mathbf{X}}^t\right)P\left({\mathbf{X}}^t\right)}=\frac{P\left({\mathbf{X}}^t|Z\right)P\left(Z\right)}{P\left({\mathbf{X}}^t\right)}=P\left(Z|{\mathbf{X}}^t\right),$

we get

$Q(\mathbf{X},{\mathbf{X}}^t)=\log p\left(\mathbf{Y}|\mathbf{X}\right)+\underset{Z}{\int }\log p\left(\mathbf{X}|Z\right)P\left(Z|{\mathbf{X}}^t\right)\mathrm{dz}+C(Z,\mathbf{Y},X^t).$

An important property of a GSM model is that it is Gaussian when conditioned on the positive scalar multiplier [19]. Hence, we have X_i|z_i~N(0⃗,z_iρ), 1≤i≤M, and

$P\left(\mathbf{X}\mid Z\right)=\underset{i=1}{\overset{M}{\Pi }}P\left(X_i\mid z_i\right)={\left(2\pi \right)}^{-\frac{d}{2}\times M}\underset{i=1}{\overset{M}{\Pi }}{\mid z_i\rho \mid }^{-\frac{1}{2}}\exp [-\frac{1}{2}\sum _{i=1}^M{X}_i^T{\left(z_i\rho \right)}^{-1}{X}_i].$

It follows that

The M-step of the EM scheme is to maximize the expectation Eq. (3) we computed in the E-step above. Let $\Omega \left(X_i^t\right)=\underset{z_i}{\int }{z}_i^{-1}P\left(z_i|X_i^t\right)d{z}_i$ . We have

Hence, we get

or

Another important property of a GSM model is that its probability density function can be expressed as a function of the quadratic form of a random vector [21]. The property means P(X^t_i)=f(r) if we let r=(X^t_i)^T ρ ^-1 X^t_i. It follows that $\Omega \left(X_i^t\right)=-\frac{2}{f\left(r\right)}\frac{\mathrm{df}}{\mathrm{dr}}$ (see the Appendix). Now we can rewrite Eq. (4) as

So far, we have obtained an iterative denoising formula Eq. (5), which only uses the derivative information of a GSM model, namely, $\frac{d}{\mathrm{dr}}\log f\left(r\right)$ . During the derivation above, we utilized two important properties of the GSM models. The first property is that a GSM random variable is Gaussian when conditioned on the scalar multiplier. This property was used in the E-step. The second property is that the probability density function of a GSM random variable can be written as a function of its quadratic form. This is true, since the GSM model is in fact a special subfamily of the Elliptically Contoured Distribution Family (ECDF) [21]. This property was used in the M-step. Considering that most existing models of wavelet coefficients of natural images belong to the GSM models, Eq. (5) hence provides a unified framework for the image denoising application. To verify Eq. (5), we will use it to derive two classical denoising algorithms. One is the Wiener filtering and the other is the Bishrinkage function proposed in [13]. Let t→+∞, Eq. (5) reduces to

We first check the Wiener filtering case. Assume X_i obeys the Gaussian distribution, namely, X_i~N(0→,ρ). Since $P\left(X_i\right)=f\left(r\right)={\left(2\pi \right)}^{-\frac{d}{2}}{\mid \rho \mid }^{-\frac{1}{2}}\exp (-\frac{1}{2}r)$ , we have $\frac{d}{\mathrm{dr}}\log f\left(r\right)=-\frac{1}{2}$ . Equation (6) hence becomes

which is just the classical Wiener filtering for Gaussian data contaminated by Gaussian white noise.

In [13], Sendur and Selesnick developed a highly effective denoising algorithm referred to as the Bishrinkage. To develop the algorithm, they used a bivariate model

where X_i=(x ₁,x ₂)^T and $\rho =\left(\begin{array}{cc}{\sigma }^2& 0\\ 0& {\sigma }^2\end{array}\right).$ .

The Bishrinkage function in [13] was derived from aMAP estimator. It is easy to show that we can also obtain the Bishrinkage function using Eq. (6). Now we have $\frac{d}{\mathrm{dr}}=\log f\left(r\right)=-\frac{\sqrt{3}}{2}{r}^{-\frac{1}{2}}$ . Substituting it into Eq. (6), we can obtain the Bishrinkage function

where (∙)₊ is defined as ${\left(x\right)}_{+}=\{\begin{array}{c}2, \mathrm{if}x<0\\ x, \mathit{otherwise}\end{array}$ .

The bivariate model Eq. (8) can be extended to multivariate cases and the resulting multivariate models belong to the Elliptically Contoured Distribution Family (ECDF) [21]. In [21], we developed an effective image denoising algorithm using the Maximum a Posterior (MAP) estimator based on these extended multivariate models. Interestingly, the algorithm developed in [21] has the same form as Eq. (6), which is the limit case of Eq. (5). Of course, it is not surprising that the EM scheme and the MAP estimator can obtain the same result, if we realize that the EM scheme is essentially a method to solve MAP (or ML) problems. We would like to emphasize that the iterative formula Eq. (5) converges due to the property of EM scheme itself. Hence, the derivation above, in fact, provides an explanation of the convergence of the algorithm developed in [21].

Next, we will simply discuss the numerical implementation of the iterative formula Eq. (5). For more details please see [21]. Let S be the symmetric square root of the positive definite matrix ρ_ε (i.e., ρ_ε=SS^T), and {Q, Λ} be the eigenvector/eigenvalue expansion of the matrix S ^-1 ρS^T, namely QΛQ^T=S ^-1 ρS^T (we note here Q^T=Q ^-1). It follows that

Hence we can rewrite Eq. (5) as

where ${\left[\frac{{\lambda }_i}{{\lambda }_i-2\frac{d}{\mathrm{dr}}\log f\left(r\right)}\right]}_1^d$ denotes the diagonal matrix $\mathrm{diag}[\frac{{\lambda }_1}{{\lambda }_1-2\frac{d}{\mathrm{dr}}\log f\left(r\right)},\dots ,\frac{{\lambda }_d}{{\lambda }_d-2\frac{d}{\mathrm{dr}}\log f\left(r\right)}]$ , and λ_i, 1≤i≤d, are the diagonal elements of Λ, namely, Λ=diag(λ ₁,…,λ_d. In addition, using S, Q and Λ, we can write r as

where V_i, 1≤i≤d, denotes the i-th element of the d-dimension vector Q^T S ^-1 X^t_i. Now we can obtain the estimate of X_i,1≤i≤M, through the iterative formula Eq. (11).

3. Unified denoising algorithm: examples

As mentioned above, most existing statistical models for natural image wavelet coefficients belong to the GSM models. In this section, we will test several classical and recently-proposed GSM prior models for the image denoising application using the unified framework developed in Section 2. According to the iterative formula Eq. (11), we only need to write a GSM prior model as a function of the quadratic form and compute the derivative of its logarithm with respect to r. We summarize the unified denoising algorithm as follows.

Unified Denoising Algorithm:

3.1. Examples for univariate prior models

In this subsection, we will focus on the usage of several univariate models (namely, d=1). Note in this case, we have r=σ ^-2 x ².

1). Univariate Gaussian Model

Although the 1-D Gaussian model is not suitable to characterize the heavy-tailed marginal statistical behavior of the wavelet coefficients of natural images, we still test this model for analysis and comparison purpose. The univariate Gaussian model is of the form

where $f_G\left(t\right)=\exp \left(-\frac{r}{2}\right).$ . Obviously, we have

2). Univariate Generalized Laplacian Model

To model the heavy-tailed property of natural image wavelet coefficients, a widely-used univariate model is the two parameter Generalized Laplacian (GL) model of the form [2, 3]

where $\Psi (s,p)=\frac{2s}{p}\Gamma \left(\frac{1}{p}\right)$ is a normalization term and Γ(x)=∫^∞ ₀ t ^x-1 e ^-t dt is the Gamma function.

The GL model is zero-mean and symmetric, and can capture the large kurtosis property of the 1-D marginal statistic of natural image wavelet coefficients [2, 3]. The typical values of p for the wavelet coefficients of natural images are in the range [0.5 1.0], corresponding to kurtosis values in the range [6, 25.2] [25]. The parameters {s,p} are directly related to the second and fourth moment, σ² and m ₄, of the random variable x. Specifically, the second and fourth moments of a generalized Laplacian signal corrupted by additive Gaussian white noise are [2]

Using Eq. (16), it is convenient to compute {s, p} from the noisy wavelet coefficients.

In order to use the iterative formula Eq. (11), we rewrite Eq. (15) as

where $\alpha ={\left(\frac{\sigma }{s}\right)}^p$, $\beta =\frac{p}{2}$ and f_GL(r)=exp(-αr^β). Obviously, we have

3). Univariate Bessel K Form Model

For p>0 and s>0, the probability density function of the univariate Bessel K Form model can be expressed as [7, 8]

where K _ν is the modified Bessel function of the second kind defined as

${\mathbf{K}}_v\left(zx\right)=\frac{\Gamma \left(v+\frac{1}{2}\right){\left(2x\right)}^v}{\sqrt{\pi }\Gamma \left(p\right)}\underset{0}{\overset{\infty }{\int }}\frac{\cos \left(\mathrm{zt}\right)\mathrm{dt}}{{\left(t^2+x^2\right)}^{v+\frac{1}{2}}}\left(\mathrm{Re}\left(v\right)>-\frac{1}{2},z>0,\mid \arg x\mid <\frac{\pi }{2}\right),$

p and s are the shape and scale parameters, respectively.

The probability density function of the univariate Bessel K Form model is unimodal, symmetric around the mode. When p<1, P(x) is sharply peaked and has heavy tails. The parameters {s,p} are directly related to the second and fourth cumulant, σ² and κ, of the random variable x. Specifically, we have σ²=ps, $\kappa =\frac{3}{p}+3$ [8, 9]. Since adding a Gaussian process to a zero-mean process only modifies its variance but not the other cumulants, we can obtain the estimate of parameters {s,p} from the noisy wavelet coefficients using the relation

where $\tilde{\sigma }$ ² and $\tilde{\kappa }$ denote the sample variance and the fourth cumulant of noisy wavelet coefficients, respectively. Using the relation $s=\frac{{\sigma }^2}{p}$ , we can rewrite Eq. (19) as

where $f_{\mathrm{BK}}\left(r\right)=r^{\frac{p}{2}-\frac{1}{4}}{\mathbf{K}}_{p-\frac{1}{2}}\left(\sqrt{2\mathrm{pr}}\right)$ . Since $\frac{d}{\mathrm{dz}}\left[z^v{\mathbf{K}}_v\left(z\right)\right]=-z^v{\mathbf{K}}_{v-1}\left(z\right)$ [26], we get

Now, we have

4). Asymptotic Univariate Bessel K Form Model

Although the Bessel K Form model Eq. (19) describes the marginal statistical behavior of natural image wavelet coefficients very well [8, 9], its functional form is not easy to work with for some applications. The Bessel K Form model with the parameters p and s can be approximated well by the function (that we call asymptotic Bessel K Form model) [8]

For denoising applications, parameters of the asymptotic Bessel K Form model can be estimated using Eq. (20) directly from the noisy wavelet coefficients. Now we can rewrite Eq. (24) as

Since σ²=ps, we get

where $f_{\mathrm{AABK}}\left(r\right)=r^{\frac{p-1}{2}}\exp \left(-\sqrt{2pr}\right)$ . Obviously, we have

5). Laplacian Model (or Bi-exponential Model)

The probability density function of the Laplacian model is of the form

where p=2s. The tails of the Laplacian model are intermediate in heaviness between the Student t-distribution and the Gaussian distribution. Since for the Laplacian model, σ²=Var(x)=2s, we have

where $f_{\mathrm{BE}}\left(r\right)=\exp (-\sqrt{2r})$ . Obviously, we have

3.2. Examples for multivariate prior models

It is well-known that a residual dependency structure always remains among wavelet coefficients, although the orthonormal wavelet was found to be fairly well-decorrelated. Recently, some multivariate models have been proposed to characterize the statistical dependence among natural image wavelet coefficients [7, 13, 15, 16, 17, 18, 19, 20, 21]. Generally, image denoising algorithms based on multivariate models can produce better denoising performance than those based on univariate models. In this subsection, we will focus on the usage of the multivariate GSM models that have an analytical expression form to the image denoising application under the unified denoising framework.

1). Multivariate Gaussian Model

The multivariate Gaussian model is of the form

where $f_G\left(r\right)=\exp \left(-\frac{1}{2}r\right)$ . The same as the univariate Gaussian case, we have

2). Multivariate Laplacian Model

The multivariate Laplacian model is of the form

where q(X)=λX^T ρ ^-1 X. Themodel is useful for representingmultivariate, sparsely distributed data, with mutually dependent components [15]. Equation (33) can be rewritten as

where $f_{\mathrm{ML}}\left(r\right)=\frac{{\mathbf{K}}_{\frac{d}{2}-1}\left(\sqrt{2r}\right)}{r^{\frac{d}{4}-\frac{1}{2}}}$ . Using Eq. (22), we can obtain

3). Multivariate Elliptically Contoured Distribution Model

In [21], we extended a generalized version of the bivariate model Eq. (8) to multivariate cases and obtained several new multivariate models that belong to the Elliptically Contoured Distribution Family (ECDF). The generalized bivariate model is of the form

where X(x ₁, x ₂)^T. Obviously, we have

The extended 4-D model is of the form

where X=(x ₁, x ₂, x ₃, x ₄)^T. Similarly, we have

The extended 10-D model is of the form

where X=(x ₁, …, x ₁₀)^T. Again, we have

4). Multivariate Exponential Model

Although the Elliptically Contoured Distribution models developed in [21] are highly effective for the image denoising application, how they fit the empirical statistics of the wavelet coefficients of natural images remains unknown. In [17], we built several new multivariate models that not only match well with the observed statistics of the wavelet coefficients of natural images, but also have analytical forms simpler than those in [21]. These models are of exponential forms and their parameters were determined by directly fitting the 1-D histogram of the quadratic form of wavelet coefficients.

The exponential models can be expressed as [17]

where X=(x ₁, …, x_d)^T and f_{MEM_d}(r)=exp(-a ₂ r ^a3). For the 2-D case, we have a ₁=-2.1, a ₂=6.8 and a ₃=0.17. For the 4-D case, we have a ₁=-3.2, a ₂=6.3 and a ₃=0.22. For the 9-D case, we have a ₁=-4.0, a ₂=5.6 and a ₃=0.26. For the 10-D case, we have a ₁=-4.5, a ₂=5.5 and a ₃=0.3. From Eq. (42), we obtain

4. Experiment results and discussion

In this section, we will test and compare the performance of prior models demonstrated in Section 3 for the image denoising application under the unified denoising framework.

We use several standard test images, namely, Barbara, Lena, Boat and Pepper for our experiments. These images are available at Portilla’s homepage http://www.io.csic.es/PagsPers/JPortilla/denoise/software/index.htm. The first three images are of size 512×512, and the last one is of size 256×256. To study the dependency of our denoising algorithm on different noise levels, these test images are contaminated by additive Gaussian white noise with different standard deviation, 10, 20 and 30. The denoising performance is evaluated both in terms of PSNR and visually. A 4-scale orthonormal wavelet transform with the Symlet-8 filter is used here. The PSNR of the restoration images versus the full range of the input noise level are listed in Table 1 and Table 2 (in both tables, 1×1+1, 3×1+1, 3×3 and 3×3+1 denote the used neighborhood structures. For example, 3×+1 denotes the used neighborhood structure containing a 3×3 local window plus a parent coefficient of the central coefficient of the window). For comparison, we also list in both tables the PSNR due to several other image denoising algorithms including the hard-threshold procedure (3σ_n) (Hard Threshold), the MMSE estimator based on the GL model proposed by Simoncelli et. al (GL_MMSE) [2], the GSM-based MMSE algorithm developed by Portilla et. al (GSM_MMSE) [19] and the Wiener filtering in the wavelet domain (Wiener).

In our experiments, we use the Successive SubstitutionMethod for the iterative formula Eq. (11) and the initialization is set to be zero (for more details about the Successive Substitution Method please see the Appendix in [21]). Figure 1(a–d) display the PSNR of the solution versus the number of iterations for the Pepper image with different noise level (20 and 30) when employing the asymptotic univariate Bessel K Form model and the multivariate exponential model (using a 3×3 neighborhood structure), respectively. Although the Successive Substitution Method is only of linear convergence, these examples clearly show the unified iterative algorithm can converge quickly. From a practical viewpoint, four or five iterations are enough to get a stable solution. For the Pepper image, when using the asymptotic univariate Bessel K Form model, the algorithm only runs 0.35 seconds for 4 iterations and 0.51 seconds for 10 iterations; when using the multivariate exponential model (using a 3×3 neighborhood structure), the algorithm only runs 1.28 seconds for 4 iterations and 2.64 seconds for 10 iterations (The platform we used is as follows: Hardware Architecture: PC; Operating System: Windows XP; Processor: 3.2 GHz; Memory Size: 2G; Application software: Matlab 7.1).

From Table 1 and Table 2, we can observe some interesting phenomena:

1). The classical Wiener filtering in the wavelet domain produces almost the same performance as the unified iterative algorithm using the Gaussian model (Gaussian_UA). One can check for this the univariate case and the multivariate case in both tables. Figure 1(e,f) display two convergence plots for the multivariate Gaussian model case (using a 3×3 neighborhood structure) for the Pepper image at noise level 20 and 30, respectively. Interestingly, only after the first iteration, the algorithm has converged to the optimal solution (in Section 2, we have noted the equivalence relation between the unified algorithm and the Wiener filtering (see Eq. (7)), which is the optimal estimate under the minimum L ₂ error sense for the Gaussian data contaminated by the additive Gaussian white noise). This verifies the validity of our unified algorithm for the image denoising application.

 

Fig. 1. PSNR of the solution versus the number of iterations for the Pepper image: For the asymptotic univariate Bessel K Form model (a) at noise level 20 and (b) at noise level 30; For the multivariate exponential model (3×3 neighborhood structure) (c) at noise level 20 and (d) at noise level 30; For the multivariate Gaussian model (3×3 neighborhood structure) (e) at noise level 20 and (f) at noise level 30

Download Full Size | PDF

2). It is well known that the univariate Laplacian model is not very suitable to model the statistics of natural image wavelet coefficients since its tails are not heavy enough. However, we surprisingly find in Table 1 and Table 2, for the univariate case, all models except the Gaussian model produce a similar performance. In most cases, the univariate Laplacian model even produces results that are a little better. It has been widely-accepted in the literature that the

Table 1. Comparison of denoising performance using different prior models in terms of PSNR for Barbara and Lena image

View Table | View all tables in this article

band-adaptive algorithms commonly perform better than the non-band-adaptive ones. It is more surprising if we realize that other models are band-adaptive — they adjust their parameters for every wavelet coefficient band, whereas the univariate Laplacian model is not. This can be explained using a relation in our denoising setting. The unified denoising algorithm operates on the quadratic form r=X^T ρ ^-1 X rather than the wavelet coefficients themselves. In the univariate case, r can be written as $r=\frac{x^2}{{\sigma }^2}$ . It follows that the expectation $E\left(r\right)=\frac{E\left(x^2\right)}{{\sigma }^2}=1$ , which, we would like to emphasize, is not band-dependent. Obviously, some non-band-adaptive algorithms can obey the relation better than the band-adaptive ones. Thus, they maybe perform better than the band-adaptive ones (if the non-band-adaptive algorithms also can obey some higher order statistic relations such as $\mathrm{Var}\left(r\right)=E\left(r^2\right)-E^2\left(r\right)=\frac{E\left(x^4\right)}{{\sigma }^4}-1$ , of course). From Table 1 and Table 2, we also can find the multivariate Laplacian model leads to a very good denoising performance (note that no parameters need to be adjusted in this case).

3). In addition, we can find that for the univariate GL model, the unified algorithm and the MMSE algorithm in [2] have a very similar performance. Generally, the PSNR difference

Table 2. Comparison of denoising performance using different prior models in terms of PSNR for Boat and Pepper image

View Table | View all tables in this article

between both algorithms is only about 0.1dB. This again implies the validity of the unified algorithm.

To compare the visual effect of different denoising algorithms, some restoration images for the Pepper image with noise standard deviation 20 are shown in Fig. 2 and Fig. 3. As usual, the hard threshold procedure produces many artifacts that blemish the restoration image seriously (see Fig. 2(a)). Interestingly, in the univariate case, the visual effect due to the GL_MMSE algorithm in [2], the unified denoising algorithm using the GL model (GL_UA), the Bessel K Form model (Bessel K Form_UA), the asymptotic Bessel K Form model (Asy. Bessel K Form_UA) and the Laplacian model (Laplacian_UA) are very similar to each other. In Fig. 3, we can find that the visual effect due to the 3×3 Wiener filtering in the wavelet domain and the unified denoising algorithm using the multivariate Gaussian model (Gaussian_UA) are almost the same. Obviously, the unified denoising algorithm based on the multivariate Laplacian model (Laplacian_UA), the multivariate Elliptically Contoured Distribution model (MECDF_UA) and the multivariate exponential model (MEM_UA) all produces a satisfactory visual effect. In order to make it clear, we display crops of denoised Boat images at noise level 20 in Fig. 4. Among them, the multivariate Laplacian model produces the best visual effect.

 

Fig. 2. Restoration images for visual comparison using different univariate prior models on the Pepper image (with additive Gaussian white noise of standard deviation 20, PSNR=22.08dB) (a) Hard threshold, PSNR=26.84dB; (b) MMSE estimator based on the GL model proposed by Simoncelli et. al (GL_MMSE) in [2], PSNR=28.24dB; (c) Unified denoising algorithm based on the Generalized Laplacian model (GL_UA), PSNR=28.31dB; (d) Unified denoising algorithm based on the Bessel K Form model (Bessel K Form_UA), PSNR=28.12dB; (e) Unified denoising algorithm based on the Asymptotic Bessel K Form model (Asy. Bessel K Form_UA), PSNR=28.07dB; (f) Unified denoising algorithm based on the Laplacian model (Laplacian_UA), PSNR=28.28dB

Download Full Size | PDF

 

Fig. 3. Restoration images for visual comparison using different multivariate prior models on the Pepper image (with additive Gaussian white noise of standard deviation 20, PSNR=22.08dB) (a) 3×3 Wiener filtering in the wavelet domain, PSNR=27.84dB; (b) Unified denoising algorithm based on the Gaussian model (Gaussian_UA: 3×3), PSNR=27.85dB; (c) GSM-based MMSE algorithm developed by Portilla et. al in [19] (GSM_MMSE: 3×3+1), PSNR=29.45dB; (d) Unified denoising algorithm based on the Laplacian model (Laplacian_UA: 3×3+1), PSNR=29.43dB; (e) Unified denoising algorithm based on the Elliptically Contoured Distribution model (MECDF_UA: 3×3+1), PSNR=29.19dB; (f) Unified denoising algorithm based on the exponential model (MEM_UA: 3×3+1), PSNR=29.31dB

Download Full Size | PDF

 

Fig. 4. Crop of restoration images for visual comparison using different multivariate prior models on the Boat image (with additive Gaussian white noise of standard deviation 20, PSNR=22.10dB) (a) GSM-based MMSE algorithm developed by Portilla et. al in [19] (GSM_MMSE: 3×3+1), PSNR=29.62dB; (b) Unified denoising algorithm based on the Laplacian model (Laplacian_UA: 3×3+1), PSNR=29.76dB; (c) Unified denoising algorithm based on the Elliptically Contoured Distribution model (MECDF_UA: 3×3+ 1), PSNR=29.58dB; (d) Unified denoising algorithm based on the exponential model (MEM_UA: 3×3+1), PSNR=29.67dB

Download Full Size | PDF

5. Conclusion

An image denoising algorithm based on natural image prior models has been developed by combining the EM scheme and the properties of the GSM models. It not only has a simple iterative form but also is computationally tractable. Due to the fact that most prior models of natural image wavelet coefficients belong to the GSM models, the developed algorithm can be viewed as a unified image denoising framework. We have tested several classical and recently-proposed prior models under the unified algorithm. We found that the developed algorithm based on the multivariate Laplacian model can produce restoration results with high quality, not only in terms of the PSNR but also visually. Considering that the main purpose of the paper is to show how the algorithm is obtained and how to use it, we did not test extensively existing prior models for natural images. We will give more test results in future. It is our hope that the developed algorithm is helpful to those who will use new natural image models for the image denoising application.