Block thresholding image denoising with dual-tree complex wavelet transform

Abstract

As a new tool of signal and image analysis, the dual-tree complex wavelet transform (DTCWT) has been found very useful for a lot of applications in machine learning and pattern recognition. In this paper, we present a method for denoising of images corrupted with additive white Gaussian noise, which employs DTCWT on noisy image to obtain complex coefficients with properties of approximate shift invariance and directional selection. Furthermore, our method adopts block thresholding scheme in denoising procedure, which can empirically choose the optimal block size and threshold at each resolution level by minimizing Stein’s unbiased risk estimate. Experimental results on several benchmark images show that our method outperforms most conventional term by term or block thresholding based methods for eliminating different levels of noise. Moreover, its denoising result is compared to state-of-the-art denoising methods with finer structures preservation.

Keywords

Image denoising undecimated dual-tree complex wavelet transform block thresholding

1 Introduction

Denoising is often used as a crucial pre-processing step in image processing, which aims to preserve important structures of the image while removing noise. So far, wavelet denoising methods with thresholding or shrinkage scheme have received a lot of attention from many researchers. Such methods perform well under a number of applications since the interest information about a signal is essentially contained in relatively small number of large wavelet coefficients and the rest of small wavelet coefficients can be set to zero. The classical thresholds include Universal threshold [1], SureShrink [2] and BayesShrink [3]. However, these methods sometimes exhibit visual artifacts since they assume that the wavelet coefficients are independent, which is not always realistic [4, 5]. Therefore, wavelet thresholding rules considering statistical dependencies among wavelet coefficients were present to solve the problem. Wang et al. [6] proposed an adaptive shrinkage denoising method by using neighbourhood characteristics, which can obtain better results than Universal threshold and SureShrink based methods. Sendur et al. [7] proposed non-Gaussian bivariate distributions by considering the dependencies between the coefficients and their parents, and the corresponding nonlinear threshold functions were derived from a Bayesian estimation based model. Experimental results validated the effectiveness of the bivariate shrinkage rules. Chen and Han [8] proposed a new threshold by capturing the dependency of inter-scale wavelet coefficients, whose values are adaptive to the changing statistics of images while maintaining the denoising efficiency. Chen et al. [9] proposed a thresholding image denoising method by considering a small neighbourhood around the customized wavelet coefficient.

Block thresholding proposed by Cai et al. [10, 11] is another scheme considering the influence of neighbourhood wavelet coefficients for denoising. In such method, the wavelet coefficients are considered within overlapping blocks. If the sum of squared wavelet coefficients in a block is larger than a given threshold, the coefficients should be retained; otherwise all of them are set to zero. Block thresholding can achieve the exact optimal rates of convergence over a range of Besov bodies. The superiority of block thresholding method has been proved via the minimax approach under the L^p risk [12, 13]. The adaptivity of block thresholding, however, depends on the choice of block size and threshold level [14]. Cai and Zhou [15] proposed a data-driven block thresholding method, which empirically chooses the block size and threshold at each resolution level by minimizing Stein’s unbiased risk estimate. The procedure is easy to implement and sharply adaptive over a class of Besov bodies. Experimental results show that it has significant advantages over the most conventional wavelet thresholding methods with fixed block sizes. Chesneau et al. [16] also proposed a similar Stein based block denoising method, and compared it with state-of-the-art denoising methods on a large set of test images. Moreover, the combination block thresholding rules were proposed by Autin et al. [17, 18], which can generate larger maxisets and achieve better denoising performance.

Note that the discrete wavelet transform (DWT) has been most commonly used in thresholding denoising methods mentioned above. However, DWT is not shift invariant and lacks phase information. To overcome the problem, Kingsbury proposed a new form of DWT, called dual-tree complex wavelet transform (DTCWT), which can provide properties of approximate shift invariance and directional selection while preserving computational efficiency [19 –21]. Since the magnitudes of complex coefficients are more strongly dependent among interscale and intrascale neighborhoods, substituting DTCWT for DWT can improve the performance of thresholding denoising method [22]. Several successful DTCWT based signal and image denoising methods can be referred to [23 –25].

Main contributions of the paper include: (a) A Stein’s unbiased risk estimate (SURE) based block thresholding image denoising method has been proposed. Shift invariance wavelet coefficients on six directional subbands have been obtained by using DTCWT. After that, the block size and threshold at each subband can be determined by minimizing SURE. (b) The proposed method has been tested on several benchmark images, and compared with some state-of-the-art methods for effectiveness. The rest of this paper is organized as follows. Section 2 reviews the theory of 2-D DTCWT. Section 3 describes SURE based block thresholding in DTCWT domain, and presents the corresponding image denoising method. Section 4 analyzes the experimental results and section 5 addresses the conclusions.

2 2-D dual-tree complex wavelet transform

The DTCWT calculates the complex transform of a signal using two separate DWT decompositions (tree a and tree b). If the filters used in one tree are specifically designed different from those in the other, it is possible for the first DWT to produce the real coefficients and the second DWT to produce the imaginary coefficients. Figure 1 shows the block diagram for a 3-level DTCWT decomposition with a different set of filters at each stage.

Note that good shift invariance of DTCWT requires proper filter design. In this paper, (13, 19)-tap biorthogonal filter set is adopted at the first decomposition level, and 14-tap linear-phase Q-Shift filter set [21] is adopted at the other decomposition levels. Given two signals with 5 samples shift (see Fig. 2a), the real coefficients computed by using DWT with Daubechies14 filters and the magnitude of the complex coefficients computed by using DTCWT filters mentioned above are shown in Fig. 2b and c respectively. The total energies at scale 5 of the two signals are 0.1979 and 0.2292 for DWT, whereas the total energies at scale 5 of the two signals are 0.9435 and 0.9446 for DTCWT, which is nearly constant.

In 2-D DTCWT, an image signal f (x, y) is decomposed by using a complex scaling function and six complex wavelet functions as follow

$\begin{matrix} f (x, y) & = & \sum_{l \in Z^{2}} A_{j_{0}, l} φ_{j_{0}, l} (x, y) + \\ \sum_{k \in α} \sum_{j = 1}^{j_{0}} \sum_{l \in Z^{2}} D_{j, l}^{k} φ_{j, l}^{k} (x, y) \end{matrix}$ (1) where A_j₀,l and $D_{j, l}^{k}$ are scaling coefficients and wavelet coefficients respectively. φ_j₀,l (x, y) denotes the scaling function and $φ_{j, l}^{k} (x, y)$ denotes six wavelet functions which are oriented at k ∈ α = {±15°, ± 45°, ± 75°}. The impulse response of the six wavelets is shown in Fig. 3.

Considering the 2-D wavelet φ (x, y) = φ (x) φ (y) associated with the row-column implementation of the wavelet transform, if φ_g (t) is approximately the Hilbert transform of φ_h (t), i.e., φ_g (t) ≈ H {φ_h (t)}, one can obtain the following six wavelets [22]:

$φ_{i} (x, y) = \frac{1}{\sqrt{2}} (φ_{1, i} (x, y) - φ_{2, i} (x, y))$ (2) $φ_{i + 3} (x, y) = \frac{1}{\sqrt{2}} (φ_{1, i} (x, y) + φ_{2, i} (x, y))$ (3)

for i = 1, 2, 3, where the two separable 2-D wavelet bases are defined as follows: $φ_{1, 1} (x, y) = φ_{h} (x) φ_{h} (y), φ_{2, 1} (x, y) = φ_{g} (x) φ_{g} (y)$ (4) $φ_{1, 2} (x, y) = φ_{h} (x) φ_{h} (y), φ_{2, 2} (x, y) = φ_{g} (x) φ_{g} (y)$ (5) $φ_{1, 3} (x, y) = φ_{h} (x) φ_{h} (y), φ_{2, 3} (x, y) = φ_{g} (x) φ_{g} (y)$ (6)

In fact, the 2-D DTCWT is implemented by taking sum/difference of two separable wavelet filter banks. Since each of the above six wavelets are aligned along a specific direction, the 2-D DTCWT can capture more image information than the conventional wavelet transform.

3 The proposed block thresholding based image denoising

We first decompose a noisy image by using the complex scaling function and six complex wavelet functions defined in Equations 2–6. Suppose we obtain the complex wavelet coefficients $D_{j, l}^{k}$ , where j denotes the resolution level (1 ≤ j ≤ j₀, j₀ is the number of decomposition level), and k denotes the orientation (k ∈ α = {±15°, ± 45°, ± 75°}). For the sake of simplicity, we index the wavelet coefficients using one number, writing d_l in place of $D_{j, l}^{k}$ , where l = 1, 2, …, N.

Considering the problem of estimating an image corrupted by additive white Gaussian noise of variance $σ_{n}^{2}$ , it can be turned into the problem into the wavelet domain: estimate y_l given d_l = y_l + η_l, $η_{l} \sim N (0, σ_{n}^{2})$ . Now we estimate y = (y₁, y₂, … y_N) based on the coefficient observations d = (d₁, d₂, … d_N) under the average mean squared error as follow $R (\hat{y}, y) = \frac{1}{N} \sum_{l = 1}^{N} E ({\hat{y}}_{l} - y_{l})^{2}$ (7)

Assume that n (n ≤ N) is the number of 2-D DTCWT coefficients at given resolution level and oriented direction, and it divisible by L for simplicity. Let m = n/L be the number of blocks, where L is the block size. Then the sum of squared empirical wavelet coefficients in the i-th block, denoted as B_i, can be calculated as follow $S_{i}^{2} = \sum_{{\underline{d}}_{i} \in B_{i}} {\underline{d}}_{i}^{2}$ (8)

where ${\underline{d}}_{i} = (d_{(i - 1) L + 1}, \dots, d_{iL}) \subseteq d$ .

In our block thresholding procedure, if $S_{i}^{2}$ is larger than given threshold λ, all the coefficients in B_i are retained; otherwise all the coefficients in B_i are set to zero. For a given block B_i the threshold estimator can be defined as follow ${\underline{y}}_{i} i (λ, L) = (1 - λ / S_{i}^{2})_{+} \cdot {\underline{d}}_{i} i, i = 1, 2, \dots, m$ (9)

Furthermore, a weakly differentiable function, denoted as $g ({\underline{d}}_{i})$ , can be defined as follow $g ({\underline{d}}_{i}) = (1 - λ / S_{i}^{2})_{+} {\underline{d}}_{i} - {\underline{d}}_{i} = {\begin{matrix} - \frac{λ}{S_{i}^{2}} {\underline{d}}_{i} λ < S_{i}^{2} \\ - {\underline{d}}_{i} λ \geq S_{i}^{2} \end{matrix}$ (10)

According to Stein’s theory, one can obtain $\begin{matrix} E_{{\underline{y}}_{i}_{i}} | | {\underline{y}}_{i} i (λ, L) - {\underline{y}}_{i} i i | |_{2}^{2} = E_{{\underline{y}}_{i} i} (SURE ({\underline{d}}_{i} i, λ, L)) \\ = L + | | g ({\underline{d}}_{i} i) | |_{2}^{2} + 2 \frac{\partial g ({\underline{d}}_{i})}{\partial {\underline{d}}_{i}} \end{matrix}$ (11)

Therefore, the total risk can be calculated as follow $E_{y} | | \hat{y} (λ, L) - y | |_{2}^{2} = E_{y} (SURE (d, λ, L))$ (12)

where $SURE (d, λ, L) = \sum_{i = 1}^{m} SURE ({\underline{d}}_{i}, λ, L)$ is an unbiased estimate of risk.

Thus, the optimal block size and threshold, denoted as L^S and λ^S, can be derived by minimizing SURE. $(λ^{S}, L^{S}) = \underset{λ, L}{arg} min (SURE (d, λ, L))$ (13)

Now, our proposed block thresholding image denoising method with dual-tree complex wavelet transform can be described as follows.

Perform 2-D DTCWT on the noisy image according to Equations (2–6), and obtain the complex coefficients at all directional subbands and resolution level,

At each resolution level and oriented direction, select the block size and threshold according to Equation (13),

Estimate all complex coefficients obtained in step (a) by using the hybrid Block James-Stein rule [15],

Perform inverse 2-D DTCWT to obtain the denoised image.

4 Experimental results

In this section, we perform extensive experiments on several benchmark images (see Fig. 4) to validate the effectiveness of our proposed method. All experiments are performed on the platform of Lenovo M7100 with Dual core 2.83 GHz CPU and 4G RAM.

4.1 Decomposition level of DTCWT

We first determine the optimal number of decomposition level for DTCWT through experiments. All six original images are added with white Gaussian noise of mean zero and variance 900 (σ_n = 30). The PSNR results of the proposed denoising method with different DTCWT decomposition level are show in Fig. 5. Taking Baboon as an example, the PSNR is 30.39 at decomposition level 1, and the PSNR improves to 30.77 at decomposition level 2. The PSNR is 32.05 at decomposition level 3, which obtains a significant improvement. The PSNR decreases to 31.60 at decomposition level 4, and further decreases to 31.37 and 31.12 respectively at decomposition levels 5 and 6. The similar results can be obtained on Barbara, Lena, Peppers, Boat and Hill. It is illustrated that the proposed method can obtain the best denoising performance when the number of DTCWT decomposition level is three.

4.2 Comparison of thresholding scheme

In this sub-section, the performance of four thresholding procedures, i.e., SureShrink [2], Bayesshrink [3], DenBlock [11] and our proposed block thresholding procedure are compared for DTCWT based image denoising. The former two procedures compute thresholds through term by term scheme, whereas the latter two procedures compute thresholds through block scheme. In the experiment, the parameters of SureShrink, Bayesshrink, and DenBlock are set as configured by their authors. While implementing our proposed method, the threshold is calculated at each resolution level, and used to process the high frequency DTCWT coefficients of images. The image denoising results of DTCWT+SureShrink, DTCWT+Bayesshrink, DTCWT+DenBlock, and our proposed method are recorded in Table 1. Note that all six original images are added with white Gaussian noise of mean zero and variance 900 (σ_n = 30), and the DTCWT decomposition level of all methods are set to 3. As is shown in the table, our proposed method indeed achieves the best PSNR results on all of the test images. This is because our proposed method determines the block size and threshold at each DTCWT resolution level by minimizing SURE. DTCWT+DenBlock obtains slightly lower PSNR since it chooses the fixed block size for thresholding, which lack of local adaptive. As for the two term by term thresholding based methods, DTCWT+ Bayesshrink outperforms DTCWT+ SureShrink on most cases, since it aims to achieve the global minimax rate of convergence, while considering the local density in a Besov class.

4.3 Comparison of conventional methods with different noisy level

In this sub-section, we implement our proposed method on the six test images varying the level of Gaussian white noise (σ_n = 10, 20, 30, 40, 50). Four other conventional denoising methods are performed for comparison: (a) the thresholding method using contourlet transform [26], (b) the multiple local thresholding method using curvelet transform [27], (c) the interscale dependency based thresholding method using DTCWT [28], (d) the thresholding method using ridgelet transform and DTCWT [29]. All methods estimate the noise variance by median estimation and implement three-level decomposition on input image to give a fair comparison.

Table 2 lists the corresponding results of the five denoising methods on six test images. It can be seen from the table that methods of [26] and [27] yield relatively lower PSNR values. This is because these two methods assume that the coefficients in transform domains are independent, which is not suitable for describing the features of most natural images. Since the method of [26] uses a single value as the global threshold, its denoising performance is worse than that of method of [27]. Compared to the former two methods, method of [28] can obtain better PSNR values except on very low noise level (σ_n = 10). This may be due to the fact that it considers the dependencies between wavelet coefficients and their parents in detail for threshold selection, and implements DTCWT to keep the shift invariant property of input images. Our proposed method is nearly 3.99 dB better than method of [26], nearly 1.67 dB better than method of [27], and nearly 0.84 dB better than method of [28] for Baboon at the average of all noise levels. The similar results can be obtained on the other five test images. It is validated that the effectiveness by combining block thresholding and dual-tree complex wavelet for image denoising. Especially when the noise level is high, the advantage of our proposed method is more obviously (see Figs. 6 and 7 for example). Note that although method of [29] yields the best PSNR values, the whole process is the most complicated since it need to partition the image into a number of overlapping squares and compute the 2-D FFT of the image. Our proposed method, however, can approximate the performance of method of [29] with less execution time.

Visual results for the six test images are shown in Figs. 8–13. Although contourlet and curvelet have certain advantages over the classical wavelet, methods of [26] and [27] introduce lower denoising performance due to the diverse and irregular patterns contained in the natural images. Obviously, methods of [28] and [29], and our method can better recover the periodic textures at different orientations in Baboon and Barbara, and provide quite a good visual quality due to the direction selection provided by DTCWT. Our method is better than method of [28] in smoothing the Lena and Pepper, which shows that the block thresholding is a good choice for preserving the smooth regions in denoising procedure. While comparing method of [29] with our method from the results on Baboon, Barbara, Lena and Peppers, we cannot easily say which of them is the winner.

5 Conclusions

This paper presents an image denoising method based on DTCWT and block thresholding. Since DTCWT can use six wavelets aligned along different directions to capture more image information than the classical wavelet transform, we first employ DTCWT on the noisy image to obtain complex coefficients with properties of approximate shift invariance and directional selection. After that, block thresholding is applied in shrinkage phase so as to fully consider the influence of neighbourhood wavelet coefficients at each subband, and the optimal block size and threshold can be empirically choosed by minimizing Stein’s unbiased risk estimate. Compared to conventional term by term or block thresholding methods, such as SureShrink, Bayesshrink and DenBlock, experimental results indicate that our method can achieve comparable better denoising effectiveness. Furthermore, our method outperforms most state-of-the-art denoising methods both in terms of signal-to-noise ratio and visual quality. Future works include combining statistical model into our proposed denoising method.

Footnotes

Acknowledgments

We are grateful to the referees for their valuable comments. This work is financially supported by the National Natural Science Foundation of China (61363050, 61272077, 61563037), the Natural Science Foundation of JiangXi Province (20142BDH80026), and the Ph.D. Programs Foundation of Nanchang Hangkong University (No. EA201504141).

References

Donoho

D.L.

and Johnstone

I.M.

, Ideal spatial adaptation by wavelet shrinkage, Biometrika81 (1994), 425–455.

Donoho

D.L.

and Johnstone

I.M.

, Adapting to Unknown Smoothness via Wavelet Shrinkage, Journal of the American Statistical Association90 (1999), 1200–1224.

Chang

S.G.

and Yu

, Adaptive wavelet thresholding for image denoising and compression, IEEE Transactions on Image Processing9 (2000), 1532–1546.

Liu

and Moulin

, Information-theoretic analysis of interscale and intrascale dependencies between image wavelet coefficients, IEEE Transactions on Image Processing10 (2001), 1647–1658.

Zhang

, Feng

, Wang

and Liu

, Image denoising via 2D dictionary learning and adaptive hard thresholding, Pattern Recognition Letters34 (2013), 2110–2117.

Shengqian

, Yuanhua

and Daowen

, Adaptive shrinkage denoising using neighbourhood characteristic, Electronics Letters38 (2002), 502–503.

Sendur

and Selesnick

I.W.

, Bivariate shrinkage with local variance estimation, IEEE Signal Processing Letter9(12) (2002), 438–441.

Chen

and Han

, Adaptive wavelet threshold for image denoising, Electronics Letters41(10) (2005), 586–587.

Chen

G.Y.

, Bui

T.D.

and Krzyżak

, Image denoising with neighbour dependency and customized wavelet and threshold, Pattern Recognition38(1) (2005), 115–124.

10.

Cai

T.T.

and Silverman

B.W.

, Incorporating Information on Neighboring Coefficients into Wavelet Estimation, Sankhya Series B63(2) (2001), 127–148.

11.

Cai

and Chicken

, Block thresholding for density estimation: Local and global adaptivity, Journal of Multivariate Analysis95(1) (2005), 76–106.

12.

Chicken

, Nonparametric regression with sample design following a random process, Communications in Statistics: Theory and Methods36(10) (2007), 1915–1936.

13.

Chesneau

, Wavelet block thresholding for samples with random design: Minimax approach under the Lp risk, Electronic Journal of Statistics1 (2007), 331–346.

14.

Chicken

and Cai

, Block thresholding for density estimation: Local and global adaptivity, Journal of Multivariate Analysis95(1) (2004), 76–106.

15.

Cai

T.T.

and Zhou

H.H.

, A data driven block thresholding approach to wavelet estimation, The Annals of Statistics37(2) (2009), 569–595.

16.

Chesneau

, Fadili

J.M.

and Starck

J.-L.

, Applied and Computational Harmonic Analysis28(1) (2010), 67–88.

17.

Autin

, Freyermuth

J.M.

and Sachs

R. von

, Ideal, denoising, within a family of tree-structured wavelet estimators, Electronic Journal of Statistics5 (2011), 829–855.

18.

Autin

, Freyermuth

J.M.

and von

, Sachs, Combining thresholding rules: A new way to improve the performance of wavelet estimators, Journal of Nonparametric Statistics24(4) (2012), 905–922.

19.

Kingsbury

N.G.

, A dual-tree complex wavelet transform with improved orthogonality and symmetry properties, International Conference on Image Processing (2000), 375–378.

20.

Kingsbury

N.G.

, Complex wavelets for shift invariant analysis and filtering of signals, Applied and Computational Harmonic Analysis10(3) (2001), 234–253.

21.

Kingsbury

N.G.

, Design of q-shift complex wavelets for image processing using frequency domain energy minimization, International Conference on Image Processing2003.

22.

Selesnick

I.W.

, Baraniuk

R.G.

and Kingsbury

N.C.

, The dual tree complex wavelet transform, IEEE Signal Processing Magazine22(6) (2005), 123–151.

23.

and Lu

C.C.

, A complex wavelet domain Markov model for image denoising, IEEE International Conference on Image Processing (2003), 365–368.

24.

Wang

, He

and Zi

, Enhancement of signal denoising and multiple fault signatures detecting in rotating machinery using dual-tree complex wavelet transform, Mechanical Systems and Signal Processing24(1) (2009), 119–137.

25.

Abdullah

H.N.

, SAR image denoising based on dual-tree complex wavelet transform, Journal of Engineering and Applied Sciences3(7) (2008), 587–590.

26.

Satheesh

and Prasad

KVSVR

, Medical image denoising using adaptive threshold based on contourlet transform, Advanced Computing2(2) (2011), 52–58.

27.

Ali

A.-D.

, Swami

P.D.

and Singhai

, Modified curvelet thresholding algorithm for image denoising, Journal of Computer Science6(1) (2010).

28.

Sendur

and Selesnick

I.W.

, Bivariate shrinkage functions for wavelet-based denoising exploiting interscale dependency, IEEE Transactions on Signal Processing50(11) (2002), 2744–2756.

29.

Chen

G.Y.

and Kégla

, Image denoising with complex ridgelets, Pattern Recognition40(2) (2007), 578–585.