Wavelet-Based Adaptive Image Denoising with Edge Preservation

Report 2 Downloads 81 Views
WAVELET-BASED ADAPTIVE IMAGE DENOISING WITH EDGE PRESERVATION Charles Q. Zhan and Lina J. Karam Telecommunications Research Center Department of Electrical Engineering Arizona State University Tempe, AZ 85287-7206 {charles.zhan,karam}@asu.edu ABSTRACT This paper presents a state-of-the-art adaptive wavelet-based denoising method with edge preservation. More specifically, a redundant discrete dyadic wavelet transform (DDWT) is performed on the noisy image to get the wavelet frame decomposition at different scales. Based on the Lipschitz regularity theory, correlation analysis across scales is performed to detect the significant coefficients from the signal and the insignificant coefficients from the noise for each subband. Different denoising techniques are applied to the significant coefficients and insignificant coefficients separately, based on different statistical models. Unlike most of the existing image denoising methods, the proposed method is able to not only shrink but also increase the magnitude of the noisy wavelet coefficients. Simulation results show that the proposed method has a remarkably superior ability to preserve the edge information and to achieve better visual quality. 1. INTRODUCTION Many types of noise can be introduced when processing or transmitting images. The problem of image denoising is to recover an image that is ”cleaner” than its noisy observation. Wavelet denoising techniques received wide attention mainly due to the breakthrough work by Donoho and Johnstone [1]. It has been shown that wavelet denoising with VisuShrink thresholding possesses nearly optimal rates of convergence over a broad family of Besov spaces [1]. Wavelet denoising with VisuShrink thresholding has been shown to be very successful for 1-D signals. However, its performance is not satisfactory for 2-D images [2]. To improve the performance of wavelet-based image denoising, several methods have been proposed in recent years. It is found that denoising based on the undecimated discrete wavelet transform (UDWT) can help in significantly reducing artifacts in the areas with discontinuities [3]. While the objective of denoising is to remove the noise from a

0-7803-7750-8/03/$17.00 ©2003 IEEE.

signal, it is also very important that the edges in the image are not blurred by the denoising operation. The Lipschitz regularity theory is widely used to detect the edge and non-edge wavelet coefficients, based on the dyadic discrete wavelet transform (DDWT) [4, 5, 6]. Several waveletbased methods, sometimes categorized as ”denoising from singularity detection”, have been reported in the literature [5, 7, 8]. Another area that attracted a lot of attention is adaptive wavelet-based denoising [9, 10, 11, 12]. This kind of methods often assume a general type of statistical model for the image wavelet coefficients. For each wavelet coefficient, the parameters of the statistical model is calculated and used to estimate the clean wavelet coefficient value. A new spatially adaptive image denoising method is presented in this paper. Correlation analysis is applied to classify the DDWT coefficients as significant and insignificant coefficients. Different statistical models and denoising techniques are applied to the significant and insignificant coefficients. The paper is organized as follows. The problems of the existing denoising methods are pointed out in Section 2. Section 3 presents the new proposed denoising method. Simulation results are given in Section 4. 2. PROBLEMS WITH EXISTING APPROACHES Let G(k) and Fˆ (k) denote the noisy and denoised wavelet coefficients value, respectively. The wavelet shrinkage type of methods [1, 11, 10] estimate Fˆ (k) from G(k) as follows: Fˆ (k) = sign(G(k))(|G(k)| − λ(k)),

(1)

where λ(k) is the thresholding value for each wavelet coefficient. The MMSE estimation type of methods [9] estimate Fˆ (k) from G(k) as follows: Fˆ (k) =

σk2

σk2 G(k), + σn2

(2)

where σk2 and σn2 are the variances of the estimated wavelet coefficient and noise, respectively.

ICIP 2003

From (1) and (2), it is quite obvious that the estimation can only shrink the magnitude of the noisy wavelet coefficients. If the noisy wavelet coefficient value has a smaller magnitude than the clean wavelet coefficient, the estimation can only make it worse. The estimators resulting from (1) and (2) are guaranteed to be biased. For the singularity detection type of methods [5, 7, 8], we have: ½ G(k), if G(k) is significant coefficient ˆ F (k) = (3) 0, otherwise While this type of methods can give unbiased estimations for the detected significant coefficients, the noise variances are not reduced for those coefficients.

Fig. 1. 2-D DDWT decomposition structure. 3.2. Robust significant coefficients detection

3. PROPOSED DENOISING APPROACH We propose a new adaptive wavelet-based with edge preservation (AWEP) denoising method that can give an unbiased estimation for the significant coefficients, while at the same time reducing the noise variances. This method can be summarized as follows: The DDWT is performed on the noisy image across different scales. Correlation analysis is applied to the detailed subbands to separate the significant coefficients from the insignificant coefficients. The significant coefficients are used to denote those wavelet coefficients that are important for the edge reconstruction; the insignificant coefficients are used to denote the rest of the wavelet coefficients, which are mainly contributed by noise. The edge direction of each significant coefficient is calculated by using the neighboring wavelet coefficients. We model each significant wavelet coefficient Gsig (k) with N (mk , σk2 ). The parameters mk , σk are calculated using the neighboring wavelet coefficients along the edge direction. We model each insignificant wavelet coefficient Ginsig (k) with N (0, σk2 ). MMSE estimation is applied to get a denoised wavelet coefficient value. Inverse DDWT is then performed to get the denoised image. 3.1. Discrete dyadic wavelet transform By applying the DDWT to the noisy image g(u, v), two detailed-subband images W 1 g(u, v, s) and W 2 g(u, v, s) and one low-frequency subband image W 3 g(u, v, s) are obtained at each scale s = 2j . The DDWT can be implemented using a filter bank structure. Fig. 1 shows the one-level DDWT of a 2-D image signal g(u, v). Only W 3 g(u, v, s) will be used for decomposition at the next scale. A 4-level DDWT decomposition is performed in the proposed method. As an example, Fig. 3(b) shows the noisy Lena image g(u, v). The three figures in the left column in Fig. 2 show the DDWT decomposition results W 1 g(u, v, s) with s = 1, 2, 4.

The significant coefficients detection method used as part of the proposed AWEP denoising technique is a modified version of the correlation method introduced in [7]. It can be summarized in the following five steps: 1) For each of the detailed subbands W l g(u, v, 2j ), l = 1, 2 and j = 0, 1, 2, a correlation measure is computed as follows: corrl g(u, v, 2j ) = W l g(u, v, 2j ) · W l g(u, v, 2j+1 ) (4) 2) Scaling is then performed to corrl g(u, v, 2j ) as: kW l g(u, v, 2j )k kcorrl g(u, v, 2j )k (5) 3) For each point (u, v), label the coefficient as a significant coefficient if corrl g(u, v, 2j ) > |W l g(u, v, 2j )|, also set corrl g(u, v, 2j ) = 0 and W l g(u, v, 2j ) = 0 in further calculations. 4) Calculate the power in the remaining 2 wavelet coefficients W l g(u, v, 2j ). Use σn,j to denote the variance of the noise wavelet coefficients at scale 2j . If P l j 2 2 u,v W g(u, v, 2 ) /M > σn,j and the iteration is still less than a pre-defined number (10 used in this paper), go back to Step 2; otherwise, proceed to Step 5. 5) Calculate the standard deviation σc of the remaining nonzero values of corrl g(u, v, 2j ); label the coefficients as significant if corrl g(u, v, 2j ) > 3σc . All the wavelet coefficients go through additional postprocessing steps, taking advantage of the geometrical properties of the image. For example, if an insignificant coefficient is surrounded by significant coefficients, then it will be reclassified as a significant coefficient and vice versa. In the three middle column figures in Fig. 2, the detected significant coefficients of W 1 g(u, v, s) are shown as white pixels. corrl g(u, v, 2j ) = corrl g(u, v, 2j ) ·

3.3. Robust edge direction calculation Based on the two detailed subband components W 1 g(u, v, s) and W 2 g(u, v, s), we use the following method to robustly

nonzero-mean Gaussian random variables for the significant coefficients. In this paper, that set of wavelet coefficients is defined to be the set of coefficients along the edge direction of the current considered significant coefficient inside a local window (a 9 by 9 local window is used in our implementation). Since the coefficient position is discrete, the edge direction might fall in between samples, in which case both the two neighboring coefficients are selected. Let the set Ns represent the selected coefficients along the edge direction, and let Ms be the number of the selected coefficients in the set Ns . For the current significant coefficient W l (u0 , v0 , 2j ), mk and σk are calculated as follows: mk =

1 Ms

X

W l g(u, v, 2j )

(u,v)∈Ns





1 σk2 = max 0, Ms Fig. 2. Intermediate results of AWEP denoising method for the noisy 512 by 512 Lena image shown in Fig. 3 (b). Left column: W 1 g(u, v, s); middle column: detected significant coefficients of W 1 g(u, v, s); right column: edge direction of significant coefficients. From top to bottom: s = 1, 2, 4. calculate the edge direction for each significant coefficient. First, a local window is defined around the current coefficient’s position, a window size of 5 by 5 is used in our implementation. Then, the absolute values of the sum of W 1 g(u, v, s) and W 2 g(u, v, s) at the current scale inside the local window are calculated, and are denoted as sum1 and sum2 , respectively. The angle (direction) is calculated as: µ ¶ sum1 Angle(u, v, s) = sign · arctan , (6) sum2 where the sign of the angle is set to 1 if W 1 g(u, v, s) and W 2 g(u, v, s) have a majority coefficients inside the local window with matching signs, i.e. sign(W 1 g(u, v, s)) = sign(W 2 g(u, v, s)); otherwise the sign is set as -1. In the three right column figures in Fig. 2, the edge directions of the detected significant coefficients of W 1 g(u, v, s) are plotted as image intensity values. In these figures, the pixels with similar intensities are the pixels with similar directions. 3.4. Gaussian model parameters calculation and MMSE estimation The first step in calculating the Gaussian model N (mk , σk2 ) parameters mk and σk is to select a set of neighboring wavelet coefficients for each significant coefficient. Those coefficients should exhibit the behavior as a realization of a

(7)

X

2  (W l g(u, v, 2j ) − mk )2 − σn,j

(u,v)∈Ns

(8) 2 where σn,j is the variance of the noise wavelet coefficients at scale s = 2j . For the insignificant coefficients, heavier smoothing can be applied. Define Nn to be the set containing all the insignificant coefficients in a bigger local window (17 by 17 used in our implementation). Let Mn be the number of the coefficients in Nn , σk is calculated as follows:   X 1 2  σk2 = max 0, , (W l g(u, v, 2j ) − mn )2 − σn,j Mn (u,v)∈Nn

(9) P where mn = M1n (u,v)∈Nn W l g(u, v, 2j ), mn is often a small number. The MMSE estimation of each significant and insignificant wavelet coefficient can be derived based on the models N (mk , σk2 ) and N (0, σk2 ), respectively, as follows: Fˆsig (k) = mk +

σk2 (Gsig (k) − mk ) σk2 + σn2

(10)

σk2 (Ginsig (k)) σk2 + σn2

(11)

Fˆinsig (k) =

4. SIMULATION RESULTS By implementing the proposed AWEP denoising method, good performance can be achieved in both the MSE sense and in terms of visual quality. Best performance can be achieved for those images containing a lot of strong edges. Table 1 compares the MSE performance of the proposed AWEP denoising method with existing denoising methods for the 512 by 512 Lena image. We generate noisy images by adding zero-mean white noise sequences with different

(b) Noisy image, σn = 20, SNR=16.44

(a) Original image

(c) Denoised image using AWEP, SNR=26.27

Fig. 3. Denoising example of the 512 by 512 Lena image.

Table 1. MSE performance comparison for the 512 by 512 Lena image. σn 10 12.5 15 20 25 AWEQ 23.7 37.6 52.1 67.5 SSMM 21.5 36.0 SANS 27.5 32.7 44.1 56.5 SAOE 20.9 26.7 32.5 44.1 56.6 AWEP (ours) 20.8 26.5 31.2 41.7 52.7 variances σn2 . The MSE results of the existing methods have been taken from [12], in which Li and Orchard summarized some of the best and newest adaptive denoising methods in the literature, which are: 1) The method proposed by Mihcak et al. [13], denoted as AWEQ(Adaptive Wiener EQ). 2) The method proposed by Liu and Moulin [10], denoted as SSMM(Scale-Space Mixture Modeling). 3) The method proposed by Chang et al. [11], denoted as SANS(Spatially Adaptive). 4) The method proposed by Li and Orchard [12], denoted as SAOE (spatially adaptive denoising under overcomplete expansion). Fig. 3(a) & 3(b) show the 512 by 512 original and noisy Lena images with σn = 20, respectively. Fig. 3(c) shows the denoised Lena image using our method. Compared to the original clean image in Fig. 3(a), it can be seen that the visual quality of the denoised image is excellent. 5. REFERENCES [1] D. L. Donoho and I. M. Johnstone, “Ideal spatial adaptation via wavelet shrinkage,” Biometrika, vol. 81, pp. 425–455, 1994. [2] P. Moulin and J. Liu, “Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors,” IEEE Tran. Info. Theory, vol. 45, no. 3, pp. 909–919, 1999.

[3] R. R. Coifman and D. L. Donoho, “Translation-invariant denoising,” in Wavelet and Statistics, A. Antoniadis and G. Oppenheim, Eds. 1995, pp. 125–150, Springer. [4] S. Mallat and S. Zhong, “Characterization of signals from multiscale edges,” IEEE Trans. Patt. Anal. Machine Intell., vol. 14, no. 7, pp. 710–732, 1992. [5] S. Mallat and W. L. Hwang, “Singularity detection and processing with wavelets,” IEEE. Tran. Inform. Theory, vol. 38, no. 2, pp. 617–643, 1992. [6] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, 2nd edition, 1999. [7] Y. Xu, J. B. Weaver, D. M. Healy, and J. Lu, “Wavelet transform domain filters: a spatially selective noise filtration techniques,” IEEE Trans. Image Proc., vol. 3, pp. 747–758, 1994. [8] T. Hsung, D. P. Lun, and W. Siu, “Singularity detection and processing with wavelets,” IEEE Trans. Info. Theory, vol. 47, no. 11, pp. 3139–3142, 1999. [9] M. K. Mihcak, I. Kozintsev, and K. Ramchandran, “Lowcomplexity image denoising based on statistical modeling of wavelet coefficients,” IEEE Signal. Proc. Lett., vol. 8, no. 12, pp. 300–303, 1999. [10] J. Liu and P. Moulin, “Image denoising based on scale-space mixture modeling of wavelet coefficients,” in ICIP, 1999, pp. 386–390. [11] S. G. Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” in ICIP, 1998. [12] X. Li and M. T. Orchard, “Spatially adaptive image denoising under overcomplete expansion,” in ICIP, 2000, pp. 300– 303. [13] M. K. Mihcak, I. Kozintsev, and K. Ramchandran, “Local statistical modeling of wavelet image coefficients and its application,” in ICASSP, 1999, pp. 3253–3256.