Image de-noising using wavelets

Image de-noising using wavelets J. N. Ellinas, T. Mandadelis, A. Tzortzis, L. Aslanoglou T.E.I. of Piraeus, Department of Electronic Computer Systems ...
7 downloads 0 Views 842KB Size
Image de-noising using wavelets J. N. Ellinas, T. Mandadelis, A. Tzortzis, L. Aslanoglou T.E.I. of Piraeus, Department of Electronic Computer Systems

Abstract In this paper, we propose an adaptive method of image de-noising in the wavelet subband domain. This approach is based on threshold estimation for each subband of the wavelet decomposition of a noise-contaminated image, by considering that the subband coefficients have a Generalized Gaussian Distribution (GGD). Under this framework, the proposed technique estimates the threshold level by applying a robust median estimator on either all the detail coefficients or every detail subband of each decomposition level. Hence, it improves the performance of the Matlab 2-D denoising function, which defines a global threshold for all the subbands of a Discrete Wavelet Transform (DWT). Also, the effectiveness of the criteria, which are the same with those of Matlab, that define the threshold levels is evaluated. The experimental evaluation of our proposition shows that it removes noise significantly and more effectively than the existed Matlab function. Keywords: Discrete Wavelet Transform; Image denoising. Περίληψη Σε αυτή την εργασία, προτείνεται μία μέθοδος για την απομάκρυνση του θορύβου από μία εικόνα, η οποία είναι προσαρμοζόμενη και εφαρμόζεται στις υποζώνες ενός μετασχηματισμού κυματιδίων. Η προτεινόμενη μέθοδος βασίζεται στην ανάλυση μίας εικόνας με τον διακριτό μετασχηματισμό κυματιδίων (DWT) και την εφαρμογή ενός κατωφλίου για την απόρριψη μέρους των συντελεστών του μετασχηματισμού. Η εκτίμηση του κατωφλίου πραγματοποιείται θεωρώντας ότι οι συντελεστές έχουν Gaussian κατανομή και με τη βοήθεια ενός εκτιμητή, ο οποίος εφαρμόζεται είτε σε όλους τους συντελεστές ή στους συντελεστές κάθε υποζώνης για κάθε επίπεδο ανάλυσης του μετασχηματισμού. Ακόμα στη παρούσα εργασία γίνεται εκτίμηση της απόδοσης των κριτηρίων, τα οποία χρησιμοποιούνται και στο Matlab, για τον προσδιορισμού του επιπέδου κατωφλίου. Τα πειραματικά αποτελέσματα δείχνουν ότι ο προτεινόμενος τρόπος απομάκρυνσης του θορύβου από μία εικόνα είναι πολύ πιο αποδοτικός από την υπάρχουσα συνάρτηση του Matlab, η οποία χρησιμοποιεί ένα ενιαίο κατώφλι για όλους τους συντελεστές του μετασχηματισμού. 1. Introduction An image is often corrupted by noise during its acquisition or transmission. The de-noising process is to remove the noise while retaining and not distorting the quality of the processed image. The traditional way of image de-noising is filtering. Recently, a lot of research about non-linear methods of signal de-noising has been developed.

1

These methods are mainly based on thresholding the Discrete Wavelet Transform (DWT) coefficients, which have been affected by additive white Gaussian noise. Since the work of Donoho and Johnstone [1]-[4], there has been a lot of research on the way of defining the threshold levels and their type (i.e. hard or soft threshold). Matlab wavelet toolbox includes functions for 1-D or 2-D de-noising [5], which are based on Donoho’s proposed thresholds. Nevertheless, in the 2-D case, there is no option for the selection of a threshold criterion and the threshold is not level dependent. In this work, a new 2-D Matlab function is developed, employing the same thresholding criteria as in Matlab, which has the option to define three different methods of threshold level estimation. According to the proposed methods, the noisecontaminated image is decomposed by a DWT. The transform coefficients within the subbands are modeled as independent identically distributed (i.i.d) random variables with Generalized Gaussian Distribution (GGD). The threshold is estimated under three different approaches and the coefficients are killed or remain unchanged or shrinked, depending on the type of thresholding (i.e. hard or soft). The first method estimates the threshold level by a median estimator, which implements the noise standard deviation from the coefficients of the diagonal subband of the first level (i.e. HH1) and is called global. The second method refers to a median estimator, which is applied on all the detail coefficients of each level, so is level dependent. Eventually, the third approach employs a median estimator which is applied on the horizontalvertical-diagonal detail coefficients of each subband, so is detail dependent. The experimental evaluations of the proposed adaptive methods of image de-noising shows that both the subjective and objective performance is superior to Matlab’s 2-D de-noising function. This paper is organized as follows. Section 2 is a brief review of the discrete wavelet transform. In section 3, the concept of wavelet thresholding is developed. Section 4 explains the proposed method of de-noising based on wavelet decomposition. Experimental evaluation is performed in section 5 and finally conclusions are given in section 6. 2. The discrete wavelet transform (DWT) The mathematical approach to the discrete wavelet transform (DWT) is based on the fact that a function f(t) can be linearly represented as: f (t ) = ∑ a kψ k (t )

(1)

k

where ak are the analysis coefficients and ψk the analyzing functions, which are called basis functions, if the above analysis is unique. If the basis functions are orthogonal, that is,

ψ k (t ),ψ l (t ) = ∫ψ k (t )ψ l (t )dt = 0 for k ≠ l

(2)

the coefficients can be estimated from the following equation: a k = f (t ),ψ k (t ) = ∫ f (t )ψ k (t ) dt

(3)

2

where f(t) is given from (1). For example, the orthogonal analysis functions for the Fourier transform are the sin(kω0t) and cos(kω0t). In general, a 2-D signal may be transformed by DWT as: f (t ) = ∑∑ α j ,kψ j ,k (t ) k

(4)

j

where αj,k and ψj,k are the transform coefficients and basis functions respectively. Equation (4) is the inverse transform, given αj,k and ψj,k. Therefore, a function f(t) may be represented by transform coefficients, which are estimated from the internal product of that function with an orthogonal basis function. Inversely, the desired function may be reconstructed from these coefficients and the basis function. These basis functions are called wavelets [6], [9], [10]. Another consideration of the wavelets is the subbnad coding theory or multiresolution analysis [11]. The signal passes successively through pairs of low pass and high pass filters, the analysis filters, which produce the transform coefficients. These coefficients, if pass successively through the synthesis filters, may reproduce the initial signal at the decoder’s side. The synthesis filters have a specific relationship with the analysis filters in order to provide perfect reconstruction and constitute the Perfect Reconstruction Bank (PRB). The orthogonal relation among the filters is their most basic property, as the produced coefficients must preserve the energy of the initial signal (Parseval’s theorem).

(Downsampling) (a)

(Upsampling) (Decomposition)

(b)

(Reconstruction)

(c)

Fig. 1 Analysis and synthesis of a 1-D signal with subband coding. (a) DWT and Inverse DWT of 1-D signal; (b) DWT tree; (c) Wavelet packet tree.

3

Fig. 1 shows the analysis of a 1-D signal with the subband coding theory, which is equivalent to wavelet decomposition. The low pass filters represent the “approximation” of the signal or its dc component and the high pass filters represent the “details” or its high frequency components. The successive analysis of the low pass component only, is called wavelet decomposition, Fig. 1(b), whereas the analysis of both the low and high pass components is called wavelet packet decomposition, Fig.1 (c). The samplers are used so that the total number of the produced coefficients to be the same as the number of signal samples. Hence, an input signal S may be equivalently analysed as: S=A3+D3+D2+D1 or S=A2+D2+D1 or S=A1+D1 Similarly, by using wavelet packet decomposition, the signal may be analysed as: S=A1+AAD3+DAD3+ADD3+DDD3 The decomposition tree is called “wavelet tree” and is transferred with its coefficients to the decoder for the reconstruction of the initial signal. The wavelet transform provides an appropriate basis for image handling because of its beneficial features. The main assets of the wavelet transform are: • The ability to compact most of the signal’s energy into a few transformation coefficients, which is called energy compaction. • The ability to capture and represent effectively low frequency components (such as image backgrounds) as well as high frequency transients (such as image edges). • The variable resolution decomposition with almost uncorrelated coefficients. • The ability of a progressive transmission, which facilitates the reception of an image at different qualities. In that sense, the existence of small coefficients is more likely to be due to the noise contamination, whereas the large coefficients contain significant image details. Hence, the small magnitude coefficients may be thresholded without affecting the large ones and therefore the quality of the image. The thresholding techniques are simple non-linear techniques that eliminate all the subband coefficients that their magnitude is under a certain threshold. The remaining coefficients remain either unaffected, which is called hard-thresholding or modified, which is called soft thresholding. The soft thresholding techniques eliminate the coefficients with magnitude less than the pre-specified threshold and shrink the rest of them. The reconstruction of the “clean” image, after the thresholding process, is performed with the inverse wavelet transform. The quality of the reconstructed image, which will contain some noise and may be distorted, is measured either subjectively by an optical evaluation or objectively by the Peak Signal to Noise Ratio (PSNR). 3. Wavelet thresholding Suppose x={xij , i=1,2,…,M and j=1,2,…,N}is an image of M x N pixels, which is corrupted by independent and identically distributed (i.i.d) zero mean, white Gaussian noise nij with standard deviation σn. The noise signal can be denoted as nij ~ N (0,σn2). This noise may corrupt the signal in a transmission channel. The observed, noise contaminated, image is y={yij , i=1,2,…,M and j=1,2,…,N}. Therefore, the noised image can be expressed as: 4

yij = xij + nij

(5)

The object of a de-noising process is to estimate image x from the noised image y, so that the Mean Square Error (MSE) to be minimum. Let W and W -1 denote the two dimensional DWT and its inverse respectively. Then, the original signal, its noised version and the noise have a matrix form in the transform domain that includes the subband coefficients. X=W x , Y=W y , V=W n

(6)

Fig. 2 shows the two level DWT of a 2-D signal, which consists of the subbands LL2 (low frequency or approximation coefficients), HL2 (horizontal details), LH2 (vertical details), HH2 (diagonal details) and the first level details HL1, LH1, HH1 [6], [8].

LL2

HL2 HL1

LH2

HH2

LH1

HH1

Fig. 2 Two level Discrete Wavelet Transform Therefore equation (5), in the spatial domain, becomes in the transform domain as follows: Y=X+V

(7)

where X, Y and V are the transform domains of the original image, its noised version and the noise respectively. The orthogonal property of the transform insures that the

5

noise in the transform domain is also of Gaussian nature. The de-noising algorithms, which are based on thresholding, suggest that each coefficient of every detail subband is compared to a threshold level and is either retained or killed if its magnitude is greater or less respectively. The approximation coefficients are not submitted in this process, since on one hand they carry the most important information about the image and on the other hand the noise mostly affects the high frequency subbands. The type of the threshold is either hard or soft. Fig. 3 indicates the two types of thresholding, which can be expressed analytically as follows. y = x if | x | > T Hard threshold:

(8) y = 0 if | x | < T

Soft threshold:

y = sign (x) ( | x | - T)

(9)

where x is the input signal, y is the signal after threshold and T is the threshold level.

y

0 (a)

x

0

0

(b)

(c)

Fig. 3 Threshold types: (a) Original signal; (b) Hard; (c) Soft The hard type does not affect the coefficients that are greater than the threshold level, whereas the soft type causes shrinkage to these coefficients. In the present work, both types of threshold are evaluated but hard thresholding may create abrupt artifacts because of its discontinuous nature. The reconstructed image is a de-noised estimate of x, which is produced by the inverse DWT.

xˆ = W −1Yˆ

(10)

where Yˆ consists of the thresholded subbands of the noised image. The threshold level is estimated by various methods called thresholding criteria, which are based on the minimization of the averaged squared error.

6

arg min[

1 N

∑i (Yˆi − X i )2 ]

(11)

where Xi and Yˆi are all the detail subbands coefficients of the original image and the noised image after thresholding respectively. 4. The proposed method of de-noising

The noise-corrupted image, which is expressed as y = x + σnn, where σn is the standard deviation of noise, according to equation (1), is subjected to a DWT. The threshold level is usually estimated for σn = 1. Therefore, when σn is different than unity there should be a threshold rescaling. Estimating the noise standard deviation in one of the following ways may perform this rescaling of the threshold level. a) A robust median estimator of the highest subband diagonal coefficients (i.e. HH1) estimates the noise standard deviation [2]. ^

σn =

median

(

Yij

0.6745

)

(12)

where Yij represents the coefficients of HH1 subband. The threshold level is rescaled according to the thresholding criterion, which is selected and is applied in every coefficient of all the subbands. This case is mentioned as global threshold. b) The same median estimator is employed in order to determine the standard deviation of noise for each level separately. In this case Yij represents all the detail coefficients of the corresponding level (i.e. coefficients of the horizontal, vertical and diagonal subbands). The threshold level is rescaled according to the selected thresholding criterion and is applied on every coefficient of the corresponding level. Hence, the number of the threshold levels is equal to the number of decomposition levels. This case is mentioned as level dependent threshold. c) The same median estimator is applied on the horizontal, vertical and diagonal detail coefficients of each level. The matrix Yij represents the corresponding detail coefficients of each level. Therefore, the number of the threshold levels is equal to the number of decomposition levels times the decomposition subbands of every level, which are three in the 2-D case. This case is mentioned as detail dependent threshold. The proposed de-noising algorithm is summed up to the following steps: • A four level DWT transforms the noise-corrupted image. • Estimate the noise standard deviation with one of the above three proposed methods by employing a thresholding criterion. • For each subband (except the low pass or approximation subband) apply hard or soft threshold to the subband coefficients. • Reconstruct the image by employing the inverse DWT.

7

5. Experimental results

The experimental evaluation is performed on two gray scale images like “Lena” and “Barbara” of size 512 X 512 pixels at different noise levels. The wavelet transform employs Daubechies’s least asymmetric compactly supported wavelet with eight vanishing moments [7] at four levels of decomposition. The objective quality of the reconstructed image is measured by: PSNR = 10 log10

2552 mse

(13)

where mse is the mean square error between the original (i.e. x) and the de-noised image (i.e. xˆ ) with size M x N: mse =

M N

1 [ x(i, j ) − xˆ (i, j )]2 MxN i =1 j =1

∑∑

(14)

The thresholding criteria, which are used for the evaluation of the proposed de-noising methods are the “rigrsure”, “heursure”, “sqtwolog”, “minimaxi”, as they are mentioned in bibliography. Firstly, the three proposed methods are evaluated. Tables 1, 2 and 3 show the performance of each method for soft thresholding and for each of the above mentioned thresholding criteria. It is apparent that the “rigrsure” thresholding criterion for all the proposed methods provides better results and for this criterion the global median estimator outperforms the level median estimator which in turn provides better results than the detail median estimator. Table 4 illustrates the comparative performance of the best global estimator, which employs the best “rigrsure” thresholding criterion, for soft and hard thresholding. These results prove that soft thresholding outperforms hard thresholding from 2.6 dB to 4 dB and for that reason it is almost always used in image de-noising applications. Table 1. Performance of the global median estimator

σ rigrsure heursure sqtwolog minimaxi rigrsure heursure sqtwolog minimaxi

Noise Standard Deviation 10 15 20 25 PSNR (dB) - LENA 33.06 31.02 29.67 28.58 33.05 30.77 29.21 28.41 27.36 25.88 24.84 24.07 28.64 27.08 25.97 25.15 PSNR (dB) - BARBARA 30.54 28.18 26.68 25.59 30.54 27.69 26.45 25.50 23.68 22.53 21.86 21.40 24.90 23.55 22.72 22.16

30

27.82 27.75 23.44 24.46 24.71 24.67 21.06 21.74

8

Table 2. Performance of the level dependent median estimator Noise Standard Deviation 10 15 20 25 PSNR (dB) - LENA 30.06 28.78 27.87 rigrsure 31.46 29.67 27.37 26.88 heursure 31.15 22.86 22.33 21.98 sqtwolog 23.34 24.09 23.46 23.01 minimaxi 24.67 PSNR (dB) - BARBARA 26.77 25.58 24.79 rigrsure 28.58 25.99 24.50 23.78 heursure 27.13 20.36 20.08 19.92 sqtwolog 20.80 21.31 20.88 20.63 minimaxi 21.90 σ

30

27.19 26.43 21.70 22.66 24.06 23.31 19.80 20.44

Table 3. Performance of the detail dependent median estimator

σ rigrsure heursure sqtwolog minimaxi rigrsure heursure sqtwolog minimaxi

Noise Standard Deviation 10 15 20 25 PSNR (dB) - LENA 30.33 29.37 27.97 27.30 30.32 28.72 27.66 27.13 22.32 22.00 21.68 21.48 23.20 22.85 22.45 22.22 PSNR (dB) - BARBARA 27.78 26.48 25.29 24.32 27.78 26.14 25.14 24 20.25 19.98 19.79 19.63 21.15 20.74 20.43 20.19

30

26.58 26.53 21.26 21.94 23.59 22.46 19.52 19.99

Table 4. Comparative performance of the best de-noising method (global median estimator and “rigrsure” thresholding criterion) for hard or soft thresholding

σ hard soft hard soft

Noise Standard Deviation 10 15 20 25 PSNR (dB) - LENA 30.39 27.46 25.87 24.47 33.06 31.02 29.67 28.58 PSNR (dB) - BARBARA 29.09 26.16 24.07 22.71 30.54 28.18 26.68 25.59

30

23.96 27.82 21.68 24.71

9

34 Matlab 2-D denoising function unique median estimator level median estimator detail median estimator

32

PSNR (dB)

30

28

26

24

22 10

12

14

16

18 20 22 Standard Deviation

24

26

28

30

Fig. 4 Comparative performance of the proposed de-noising methods

Fig. 4 demonstrates the performance of the three proposed methods, which employ soft thresholding and the “rigrsure” thresholding criterion, versus the Matlab 2-D denoising function for the test image “Lena”. All the proposed methods present significant better performance, which ranges from 3 dB to 6 dB for the whole examined range, than the Matlab 2-D de-noising function. Finally, Fig. 5 shows the effect of the compared methods to a noise contaminated image with noise of standard deviation 30 (strong interference). It is obvious that the subjective optical quality of the de-noised image is the best for the global median estimator, becomes slightly worse for the level estimator, the detail estimator and much worse for the Matlab de-noising function. It is quite apparent that the Matlab’s de-noising function creates image blurring and some annoying artifacts, Fig. 5(c), whereas the proposed methods reduce significantly these effects and give to the denoised image smoother feeling, Fig. 5(d, e, f). The best optical quality is presented by the global median estimator and the reason is that it produces more accurate estimate of the noise standard deviation as the diagonal coefficients of the first decomposition subband seem to have statistical properties closer to noise. The rest two of the proposed methods present slightly worse characteristics, but the subjective optical quality of the de-noised image is quite good.

10

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 5 Comparative performance of the proposed methods. (a) Original image; (b) noised image; (c) Matlab de-noising; (d) Global estimator; (e) Level estimator; (f) Detail estimator

11

6. Conclusions

In this paper, a simple global and two subband adaptive median estimators are proposed and evaluated with respect to the Matlab’s 2-D de-noising function in order to recover an image from noise contamination effectively. They are all based on the wavelet decomposition of the image and the Generalized Gaussian Distribution modeling of the subband coefficients. The threshold estimation is either detail or subband level dependent. The proposed methods are tested for various thresholding criteria and use soft thresholding to provide smoothness and better edge preservation, avoiding the discontinuity character of the hard thresholding methods. The experimental evaluation showed that the proposed methods have far better performance than the Matlab 2-D de-noising function that uses global thresholding.

12

References

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]

D.L. Donoho, De-noising by soft thresholding, IEEE Trans. on Info. Theory, pp. 933-936, 1993. D.L. Donoho, I.M. Johnstone, Ideal spatial adaptation via wavelet shrinkage, Biometrica, vol. 81, pp. 425-455, 1994. D.L. Donoho, I.M. Johnstone, Adapting to unknown smoothness via wavelet shrinkage, Journal of American statistical assoc., vol. 90, no. 432, pp. 12001224, 1995. D.L. Donoho, I.M. Johnstone, Wavelet shrinkage: Asymptopia, J.R. Stat. Soc., series B, vol. 57, no. 2, pp. 301-369, 1995. Matlab Wavelet Toolbox User’s Guide, MathWorks, 1996. M. Vetterli, J. Kovacevic, Wavelets and subband coding, Englewood Cliffs, NJ, Prentice Hall, 1995. I. Daubechies, Ten lectures on wavelets, Proceedings CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM, vol. 61, 1992. J.N.Ellinas, M.S.Sangriotis, Modern Techniques of Image Compression, 2nd Conference of « Archipelagos Technologies », Piraeus, April 2002. C.S.Burrus, R.A.Gopinath, H.Guo, “Introduction to Wavelets and Wavelet Transforms”, Prentice Hall, 1998, pp. 2-18. G. Strang, T. Nguyen, “Wavelets and Filter Banks”, Wellesley,1997. S. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation”, IEEE Trans. PAMI, vol. 11, no. 7, pp. 674-693.

13

Suggest Documents