Image Denoising Based on Bivariate Distribution

The literature has shown that the performance of the de-noising algorithm was greatly influenced by the dependencies between wavelet coefficients. In this paper, the bivariate probability density function (PDF) was proposed which was symmetric, and the dependencies between the coefficients were considered. The bivariate Cauchy distribution and the bivariate Student’s distribution are special cases of the proposed bivariate PDF. One of the parameters in the probability density function gave the estimation method, and the other parameter can take any real number greater than 2. The algorithm adopted a maximum a posteriori estimator employing the dual-tree complex wavelet transform (DTCWT). Compared with the existing best results, the method is faster and more efficient than the previous numerical integration techniques. The bivariate shrinkage function of the proposed algorithm can be expressed explicitly. The proposed method is simple to implement.


Introduction
Image denoising is very important in many image processing applications. The methods of de-noising using wavelet transform have been widely proposed. Chang et al. proposed a subband adaptive method based on the generalized Gaussian distribution (GGD) model in [1]. Mihçak et al. in [2] gave the Bayesian linear least mean square error estimation (LMMSE) method by using the central square window. Portiilla et al. in [3] proposed a stochastic model by a hidden random field with an overcomplete wavelet pyramid. In [4], the method of combining dual tree complex wavelet with total variation is studied to eliminate speckle noise and Gaussian noise in ultrasonic image. In [5], a threshold estimation denoising method based on double density wavelet transform (DDWT) was proposed. Shen et al. in [6] proposed a denoising algorithm by constructing nonseparable Parseval frames. The paper in [7] used the designed tight frame wavelet filters for image de-noising.
The de-noising method can effectively improve the processing performance by considering the dependence between wavelet subbands. Senduer and Selesnick in [8,9] have developed a new non-Gaussian bivariate distribution, and a scheme was developed which considers the statistical dependency between the dual-tree complex wavelet coefficients.Achim et al. in [10] designed a bivariate Bayesian estimator using bivariate α-stable distributions. A spatially adaptive method had been introduced using the Cauchy Prior for the speckle reduction in [11,12], and developed method using Gaussian probability density function (PDF) in [13] by Bhuiyan et al., respectively. Ranjani and Thiruvengadam in [14,15] proposed a despeckling algorithm using the bivariate Cauchy probability density function (PDF) by the DTCWT subbands, respectively. In [16], a denoising speckle algorithm was proposed using a heavy tailed Levy distribution based on DTCWT for the ultrasound images. A denoising algorithm was developed using scale mixtures of Rayleigh distribution in [17]. The dual-tree complex wavelet (DTCWT) is tight frame wavelet in [18,19]. The DTCWT have advantages that orthogonal bases can not have for image de-noising [20][21][22].
In this paper, a bivariate symmetric probability density function (PDF) was proposed, which utilizes the dependence of DTCWT wavelet coefficients of each subband. A novel noise reduction algorithm has been introduced for noisy images.
The organization of this paper is as follows. Section 2 describes a noise reduction algorithm. The expression of the parameter for the distribution in each wavelet subband is also derived. The standard maximum posterior (map) estimation of undamaged wavelet coefficients is given in Section 2. The experimental results of the proposed method are discussed in Section 3. Finally, conclusions are given in Section 4.

Proposed Algorithm
The paper discusses the denoising of the image with additive white Gaussian noise. Let ω 1k represent the kth wavelet coefficient,and ω 2k represent the wavelet coefficient at the next scale. y 1k and y 2k are noisy observations of ω 1k and ω 2k , and n 1k and n 2k are noise samples.
where w k = (w 1,k , w 2,k ), y k = (y 1,k , y 2,k ), and n k = (n 1,k , n 2,k ) in the wavelet domain. To improve readability, the coefficient index k is omitted. In this paper, we proposed a joint bivariate PDF as where γ is the dispersion parameter. The joint bivariate symmetric probability density function (PDF) for the coefficients was illustrated in Figure 1. In the previous literature, most of the exponential functions were used as densities for noise removal [8,9,13,16], and the exponential functions attenuated rapidly. The density of this paper is composed of power functions when the parameter θ is given.
The density of this paper was composed of power functions when the parameter θ was given. The attenuation speed of power functions is relatively slow and the density has a thick tail. We know that most of the practical application problems have thick tails. The density function of this paper is more consistent with the actual problems and is more conducive to dealing with practical problems. The bivariate Cauchy distribution and the bivariate student's distribution are special cases of the bivariate PDF in (2).
The value of parameter θ in (2) has a relatively large influence on the experimental results. Different kinds of image corresponding θ values have different noising effects. For example, this paper used θ = 4 to reduce boat image noise. When θ = 4, the marginal probability density function of bivariate PDF f w (w 1 , w 2 ) is calculated as follows Using the following identity in [23] +∞ 0 y µ−1 The solution of integral in (3) is as where the gamma function Similarly, the marginal density of w 2 can be obtained. Because the product of two marginal densities is not equal to the joint density, this indicates that w 1 and w 2 are dependent. To motivate this model, see Figure 2. The solid line showed the fitted curve of the original high-band coefficients of boat image from the first scale, obtained by employing and the marginal PDF in (5) and DTCWT in [8], and the dispersion parameter γ values for θ = 4 were estimated by local marginal standard deviations σ rearranging (7) which are described later. The dashed line shows the fitted curve by Gaussian probability density function (PDF). Observe that the normalized histogram is well approximated by the marginal PDF of the bivariate PDF in (5) .
The independent bivariate Gaussian PDF of noise n k is.
where σ 2 n is the variance of Gaussian noise. The parameters γ for the bivariate PDF in (2) are calculated using the marginal variance σ 2 of the kth wavelet coefficient as where θ > 2. When θ ≤ 2, the marginal variance does not exist. The local marginal variance σ 2 of the kth coefficient is calculated using the central square window Λ as the following [8] where N Λ is number of points in window Λ, σ 2 n be the noise marginal variance. The estimator in [24] is used to estimate the noise variance σ 2 n σ n = median y j 0.6745 , y j ∈ subband HH.
Rearranging (7), the estimator of the parameter γ for the bivariate PDF is derived aŝ The standard maximum a posteriori (MAP) estimator of w for the observation y iŝ It is just sufficient to maximize p n (y|w)p w (w) The MAP estimator in (2) can also be calculated aŝ The probability density functions of signal and noise defined in (2) and (6) are substituted into (13), respectively, we getŵ = arg max where The solution to (14) can be obtained by solving the following equations: Rearranging (15) and (16), we get w 2 = w 1 y 2 y 1 , and by substituting w 2 into (15), it is rewritten as y 2 1 + y 2 2 w 3 1 − y 1 y 2 1 + y 2 2 w 2 1 + γ 2 y 2 1 + (θ + 2)σ 2 n y 2 1 w 1 − γ 2 y 3 1 = 0.
According to the method of [25], the real solution of the cubic equation in Equation (17) is as follows where p = −18γ 2 + 9(θ + 2)σ 2 n 54 y 2 1 + y 2 2 − 1 27 The shrinkage function is obtained as followinĝ where (g) + is defined as  (19), as shown in the figure, there was a dead zone (the dead zone is the area where the estimate is zero), that is, The denoising image was obtained by inverse DTCWT of the estimated wavelet coefficients by (19). The block diagram showing the proposed de-noising process steps was given in Figure 4.
The denoising algorithm is described as follows: (1) First implement the DTCWT to get y k .

Experimental Results
In the experiments, the proposed density is used to give the algorithm, and it is compared with other famous image restoration models such as using Gaussian probability density function (PDF) [2], using the GSM (Gaussian scale mixture) algorithm [3], using non Gaussian probability density function (PDF) [8], using the tight frame in [6]. The results show that the denoising effect of the new model is effective compared with the above methods. Table 1 show the performance of the proposed algorithm and other literature methods. The peak signal-to-noise ratio(PSNR) is used to test the performance of these algorithms. To do this, we use different noise levels σ 2 n for the Barbara, boat, Lena images. Let s, d represent the original image and the denoised image. The root mean-square (rms) error is determined by (20) where N is the number of pixels. The PSNR is determined by PSNR = 20 log 10 255 rms .
The proposed algorithm in this paper uses the filter in [8] and the DTCWT in [18,19], which are tabulated in " Proposed-Method" column of Table 1. The number of decomposition levels is 5, using 7 × 7 window size N Λ . The parameter θ of the shrinkage function in (19) affects the noising effect. As the parameter θ values are manually selected prior to running a simulation, the user may have to experiment with different values to achieve the highest increase in noise removal effect. In our experiments, the parameter θ takes 11 when the standard deviation σ n = 10, 15, and θ takes 3 when the standard deviation σ n = 20, 25, 30 for lena image. The parameter θ takes 11 for all noise levels of the Barbara image, and the parameter θ takes 4 for all noise levels of the Boat image. The PSNR is the average value of the data obtained by running 100 times the proposed algorithm . We also compare the algorithm with other methods reported in the literature. Here, we give the results compared with the other four methods. In Table 1, the PSNR values using the denoising algorithm in [2] with the exponential prior are tabulated in "Method in [2]" column, the PSNR values using the GSM (Gaussian scale mixture) algorithm described in [3] are tabulated in "Method in [3]" column, the PSNR values with the DTCWT in [8] are tabulated in "Method in [8]" column, and the PSNR values with the tight frame in [6] are tabulated in "Method in [6]" column, separately. We indeed found out that a new model improves the denoising performance . In [8], the Matlab implementation for that algorithm takes 25 s for a 512 × 512 image on 450-MHz Pentium II, the Matlab program for the proposed algorithm takes 2.51 s for a 512 × 512 image on core(TM) i3 M 350.
The denoised image using a 512 × 512 Barbara image is shown in Figure 5. The original and the noisy images are shown in Figure 5a,b. The denoised image using the method in [8] was illustrated in Figure 5c and had PSNR value of 28.52 dB, and the denoised image using the proposed method was illustrated in Figure 5d and had PSNR value of 28.7222 dB, with the Gaussian noise standard deviation 25, respectively. Other denoised image using a 512 × 512 boat image is shown in Figure 6. The original and the noisy images are shown in Figure 6a,b. The denoised image using the method in [8] was illustrated in Figure 6c and had PSNR value of 28.93 dB, and the denoised image using the proposed method was illustrated in Figure 6d and had PSNR value of 29.0215 dB, with the Gaussian noise standard deviation 25, respectively. Table 1. Peak signal-to-noise ratio values of denoised images different test images and noise levels (σ n ) of noisy. [2] Method in [6] Method in [3] Method in [8]

Conclusions
In this paper, a novel bivariate symmetric probability density function (PDF) was proposed that was more suitable for actual images than Gaussian probability density function. The closed expression for the parameter γ of the bivariate PDF is derived. Based on the dependence of dual-tree complex wavelet subbands, a new noise reduction algorithm is proposed. The parameter θ values in this algorithm can not only take integers, but also take any real number greater than 2. Three kinds of images with different texture features were used to test the effectiveness of the proposed algorithm. The experiments show that the parameter θ values may exist corresponding optimal values for noise removal with different noise levels for images, it is natural to ask which values are best suited for the task. The calculation method of the optimal parameter θ values for each noise level image are explored.
Our results indicate that the algorithm produces denoised images with less ringing artifacts at less computational complexity. The experimental results show that considering the statistical dependencies between wavelet coefficients and its other neighbors can significantly improve the noise removal effect. By improving the statistical dependencies method and considering the inter scale and intra scale dependencies of wavelet coefficients [26][27][28][29][30], the work of the paper will be further developed.
Author Contributions: Conceptualization, methodology, P.Z.; software, validation, formal analysis, investigation, resources, data curation, X.Z. and C.Z.; writing-original draft preparation, P.Z.; writing-review and editing, X.Z. and C.Z.; supervision, project administration, funding acquisition, P.Z. All authors have read and agreed to the published version of the manuscript.