Wavelets and LPG-PCA for Image Denoising

Abstract

In this chapter, a new image denoising approach is proposed. It combines two image denoising techniques. The first one is based on a wavelet transform (WT), and the second one is a two-stage image denoising by PCA (principal component analysis) with LPG (local pixel grouping). In this proposed approach, we first apply the first technique to the noisy image in order to obtain the first estimation version of the clean image. Then, we estimate the noise-level from the noisy image. This estimate is obtained by applying the third technique of noise estimation from noisy images. The third step of the proposed approach consists in using the first estimation of the clean image, the noisy image, and the estimate of the noise-level as inputs of the second image denoising system (LPG-PCA). A comparative study of the proposed technique and the two others denoising technique (one is based on WT and and the second is based on LPG-PCA), is performed. This comparative study used a number of noisy images, and the obtained results from PSNR (peak signal-to-noise ratio) and SSIM (structural similarity) computations show that the proposed approach outperforms the two other denoising approaches (the first one is based on WT and the second one is based on LPG-PCA).

1. Introduction

In the image acquisition process, the noise will be inevitably introduced so denoising is a necessary step for ameliorating the image quality. As a primary low-level image processing, noise suppression has been extensively studied, and numerous denoising approaches have been proposed, from the earlier frequency domain denoising approaches and smoothing filters [1] to the lately developed wavelet [2, 3, 4, 5, 6, 7, 8, 9, 10, 11], curvelet- [12] and ridgelet- [13] based approaches, sparse representation [14] and K-SVD approaches [15], shape-adaptive transform [16], bilateral filtering [17, 18], nonlocal mean-based techniques [19, 20], and nonlocal collaborative filtering [21]. With the quick development of modern digital imaging devices and their increasingly broad applications in our daily life, there are rising necessities of new denoising techniques for higher quality of image. The WT (wavelet transform) [22] proved its effectiveness in noise cancelation [2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. This transform decomposes the input signal into multiple scales which represent different time-frequency components of the original signal. At each scale, some operations, such as statistical modeling [4, 5, 6] and thresholding [2, 3], can be applied for canceling noise. Noise reduction is performed by transforming back the processed wavelet coefficients into spatial domain. Late development of WT-based denoising techniques includes ridgelet- and curvelet-based techniques [12, 13] for line structure conservation. Despite WT proved its effectiveness in denoising, it uses a fixed wavelet basis (with translation and dilation) for representing image. However, for natural images, a rich amount of different local structural patterns exists and therefore cannot be well represented by using just one fixed wavelet basis. Consequently, WT-based techniques can generate many visual artifacts in the denoising output. To overcome the problem of WT, in [23], Muresan and Parks proposed a spatially adaptive principal component analysis (PCA)-based denoising technique, which computes the locally fitted basis for transforming the image. In [15, 16], Elad and Aharon proposed K-SVD-based denoising approach and sparse and redundant representation by training a highly over-complete dictionary. In [16], Foi et al. applied a shape-adaptive discrete cosine transform (DCT) to the neighborhood, which can attain very sparse representation of the image and consequently lead to efficient denoising. All these approaches proved better denoising performance than classical WT-based denoising techniques. The NLM (nonlocal means) schemes used a very different philosophy from the above approaches in noise cancelation. The NLM idea can be traced back to [24], where the similar image pixels are averaged according to their intensity distance. Similar ideas were used in the bilateral filtering schemes [17, 18], where both the spatial and intensity similarities are exploited for pixel averaging. The NLM denoising framework was well established in [19]. In the image, each pixel is estimated as the weighted average of all the pixels and the weights are determined by the similarity between the pixels. This scheme was improved in [20], where the pair-wise hypothesis testing was used in the NLM estimation. Inspired from the success of NLM schemes, Dabov et al. [21] proposed a collaborative image denoising scheme by sparse 3D transform and patch matching. They look for similar blocks in the image by using block matching and grouped these blocks into a 3D cube. Then, a sparse 3D transform was applied to that cube, and noise was canceled by Wiener filtering in the transformed domain. The so-called BM3D approach attains remarkable denoising results, yet its implementation is a little complex. Lei Zhang et al. [25] have presented an efficient PCA-based denoising approach with local pixel grouping (LPG). PCA is a classical de-correlation technique in statistical signal processing, and it is pervasively used in dimensionality reduction and pattern, etc. [26]. The original dataset is transformed into PCA domain, and only the different most significant principal components are conserved. Consequently, trivial information and noise can be eliminated. In [23], a PCA-based scheme was proposed for image denoising by using a moving window for computing the local statistics, from which the local PCA transformation matrix was estimated. This technique applies PCA directly to the noisy image without data selection, and much residual noise and visual artifacts appear in the denoised image. In the LPG-PCA-based technique, Lei Zhang et al. [25] modeled a pixel and its nearest neighbors as a vector variable. The training samples of this variable are chosen by grouping the pixels with similar local spatial structures to the underlying one in the local window. With such an LPG technique, the local statistics of the variables can be accurately calculated so that the image edge structures can be well conserved after shrinkage in the PCA domain for noise suppression. The LPG-PCA scheme proposed in [25] has two stages where the first stage yields an initial image estimation by eliminating the most of the noise and the second stage will further refine the first stage output. The two stages have the same procedures except for the noise-level parameter. Since the noise is significantly reduced in the first stage, the LPG precision will be much improved in the second stage so that the final denoised image is visually much better. When compared with WT which uses a fixed basis function for decomposing the noisy image, the proposed LPG-PCA approach is a spatially adaptive image representation so that it can better characterize the image local structures. When compared with BM3D and NLM approaches, the LPG-PCA-based technique proposed in [25] can use a relatively small local window to group the similar pixels for PCA training, yet it yields competitive results with state-of-the-art BM3D algorithm. In this paper we propose a new image denoising approach which combines the dual-tree discrete wavelet transform (DT-DWT)-based denoising approach [12] and the two-stage image denoising technique by PCA with local pixel grouping (LPG) [25]. To evaluate this proposed technique, we have compared it to the two techniques (the DT-DWT-based denoising technique [12] and LPG-PCA [25]). This comparison is based on PSNR and SSIM computation. In the rest of this paper, we first will deal with PCA. Then, we will be interested in DT-DWT [12]. After that we will deal with noise-level estimation from the noisy image proposed in [27, 28]. Then, we will present the two denoising techniques proposed in [12, 25]. After that we will detail the proposed image denoising technique, and finally we will give results and evaluation.

2. Principal component analysis (PCA)

Let

$$X=\left[\begin{array}{lll}{x}_1^1& x_1^2\dots & x_1^n\\ x_2^1& x_2^2\dots & x_2^n\\ \begin{array}{l}.\\ x_n^1\end{array}& \begin{array}{l}.\\ x_n^2\end{array}& \begin{array}{l}.\\ x_n^n\end{array}\end{array}\right]$$E1

$$X_i=\left[x_1^1{x}_1^2\dots x_1^n\right]$$E2

representing the sample vector of $x_i$. The mean value of $X_i$ is computed as follows:

$${\mu }_i=\frac{1}{n}\sum _{j=1}^n{X}_i\left(j\right).$$E3

And then, the sample vector is centralized as follows:

$${\overline{X}}_i=X_i-{\mu }_i=\left[{\overline{x}}_i^1{\overline{x}}_i^2\dots {\overline{x}}_i^n\right]$$E4

with ${\overline{x}}_i^j=x_i^j-{\mu }_i$. Accordingly, the centralized matrix of $X$ is expressed as follows:

$$\overline{X}=X_i-{\mu }_i={\left[X_i^T{X}_2^T\dots {\overline{x}}_m^T\right]}^T.$$E5

Finally, the covariance matrix of the centralized dataset is computed as follows:

$$\mathrm{\Omega }=\frac{1}{n}\overline{X}{\overline{X}}^T.$$E6

The PCA aim consists in finding an orthonormal transformation matrix $P$ in order to de-correlate $\overline{X}$, i.e., $\overline{Y}=P\overline{X}$, so that the matrix of covariance of $\overline{Y}$ is diagonal. Since the covariance matrix $\mathrm{\Omega }$ is symmetrical; therefore, it can be expressed as follows:

$$\mathrm{\Omega }=\varphi \mathrm{\Lambda }{\varphi }^T$$E7

where $\mathrm{\Lambda }=\mathit{diag}\left\{{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_m\right\}$ is the diagonal eigenvalue matrix with ${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_m$ and $\varphi =\left[{\varphi }_1{\varphi }_2\dots {\varphi }_m\right]$ represents the $m\times m$ orthonormal eigenvector matrix. The terms ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_m$ and ${\varphi }_1,{\varphi }_2,\dots ,{\varphi }_m$ are, respectively, the eigenvalues and the eigenvectors of $\mathrm{\Omega }$. By setting the matrix P as follow:

$$P={\varphi }^T,$$E8

$\overline{X}$ can be de-correlated, i.e., $\mathrm{\Lambda }=\frac{1}{n}{\overline{Y}\overline{Y}}^T$ and $\overline{Y}=P\overline{X}$. An interesting property of PCA is that it fully de-correlates the original dataset $\overline{X}$. In general, the signal energy will concentrate on a small subset of the PCA transformed dataset, whereas the noise energy will evenly spread over the whole dataset. Consequently, the noise and signal can be better distinguished in the domain of PCA.

3. LPG-PCA denoising algorithm

3.1. Modeling of spatially adaptive PCA denoising

In [25] and in previous literature, the noise $\upsilon$ degrading the original image $I$ is supposed to be white and additive with standard deviation $\sigma$ and zero mean, and the noisy image,$I_{\upsilon }$, is expressed as follows:

$$I_{\upsilon }=I+\upsilon$$E9

Both noise $\upsilon$ and image $I$ are supposed to be uncorrelated. The purpose of image denoising consists in estimating the clean image $I$ from $I_{\upsilon }$, and the estimate is denoted by $\widehat{I}$. The latter is expected to be as close as possible to the original image, $I$. Two quantities describe an image pixel. Those quantities are its intensity and the spatial location. However, the image local structure is represented as a set of neighboring pixels at different intensity levels. As most of the semantic information of an image is conveyed by its edge structures, edge conservation is highly required in denoising of this image. In [25], the pixel and its nearest neighbors were modeled as a vector variable and perform denoising on the vector instead of the single pixel. According to Figure 1, for an underlying pixel to be denoised, Lei Zhang et al. [25] set a $K\times K$ window centered on it, and denoted by $x={\left[x_1,\dots ,x_m\right]}^T,m=K^2$, the vector contains all components within the window.

Figure 1.

Illustration of the modeling of LPG-PCA-based denoising [25].

As the observed image is the original image degraded by the noise, they denote the noisy vector of $x$ by [25]:

$$x_{\upsilon }=x+\upsilon$$E10

$\mathrm{x}$ where $\upsilon ={\left[{\upsilon }_1\dots {\upsilon }_m\right]}^T$, $x_k^{\upsilon }=x_k+{\upsilon }_k,k=1,\dots ,m$, and $x_{\upsilon }=\left[x_1^{\upsilon }\dots x_m^{\upsilon }\right]$. For estimating $x$ from the noisy vector, $x_{\upsilon }$, they are viewed as (both noiseless and noisy) vector variables so that one can use the statistical techniques such as PCA. For canceling the noise from the noisy vector $x_{\upsilon }$ by using $\mathrm{PCA}$, a set of training samples of $x_{\upsilon }$ is needed so that the covariance matrix of $x_{\upsilon }$ and therefore the PCA transformation matrix can be computed. For this aim, Lei Zhang et al. [25] have used an $L\times L\left(L>K\right)$ training block centered on $x_{\upsilon }$ in order to find the training samples, as illustrated in Figure 1. The simplest manner consists in taking the pixels in each possible $K\times K$ block within the $L\times L$ training block as the samples of noisy variable $x_{\upsilon }$. In this way, for each component $x_k^{\upsilon }$ of $x_{\upsilon }$, there are in total ${\left(L-K+1\right)}^2$ training samples. Though, there can be very different blocks from the given central $K\times K$ block in the $L\times L$ training window, taking all the $K\times K$ blocks as the training samples of $x_{\upsilon }$ will lead to inaccurate estimation of the matrix of covariance of $x_{\upsilon }$, which subsequently leads to inaccurate estimation of the PCA transformation matrix and finally results in much residual noise. Consequently, selecting and grouping the training samples that are similar to the central $K\times K$ block are required before image denoising by applying the PCA transform.

3.2. Local pixel grouping (LPG)

Grouping the training samples similar to the central $K\times K$ block in the $L\times L$ training window is certainly a problem of classification, and therefore different grouping techniques such as correlation based matching, block matching, K-means clustering, etc. can be used based on different criteria. The block matching-based technique may be the simplest but very efficient one, and it is used in [25] for LPG. There are totally ${\left(L-K+1\right)}^2$ possible training blocks of $x_{\upsilon }$ in the $L\times L$ training window. We will denote ${\overline{x}}_0^{\upsilon }$ in the column sample vector which contains the pixels in the central $K\times K$ block, and denoted by ${\overline{x}}_i^{\upsilon },i=1,2,\dots ,{\left(L-K+1\right)}^2-1$, the sample vectors correspond to the other blocks. Let ${\overline{x}}_i$ and ${\overline{x}}_0$ be, respectively, the associated noiseless sample vectors of ${\overline{x}}_i^{\upsilon }$ and ${\overline{x}}_0^{\upsilon }$. It can be simply computed that

$$e_i=\frac{1}{m}\sum _{k=1}^m{\overline{x}}_0^{\upsilon }\left(k\right)-{\overline{x}}_i^{\upsilon }{\left(k\right)}^2\approx \frac{1}{m}\sum _{k=1}^m{\overline{x}}_0\left(k\right)-{\overline{x}}_0{\left(k\right)}^2+2{\sigma }^2$$E11

In Eq. (11), the fact that noise $\upsilon$ is white and uncorrelated with signal is used. With Eq. (11), if we have the following condition

$$e_i<T+2{\sigma }^2$$E12

where $T$ designates a preset threshold, then we select ${\overline{x}}_i^{\upsilon }$ as a sample vector of $x_{\upsilon }$. Assume that $n$ sample vectors of $x_{\upsilon }$ are selected including the central vector ${\overline{x}}_0^{\upsilon }$. For the expression convenience, these sample vectors are denoted as ${\overline{x}}_0^{\upsilon },{\overline{x}}_1^{\upsilon },\dots ,{\overline{x}}_{n-1}^{\upsilon }$. The noiseless counterparts of those vectors are denoted as ${\overline{x}}_0,{\overline{x}}_1,\dots ,{\overline{x}}_{n-1}$, accordingly. Then, the training dataset for $x_{\upsilon }$ is constituted by.

$$X_{\upsilon }=\left[{\overline{x}}_0^{\upsilon },{\overline{x}}_1^{\upsilon },\dots ,{\overline{x}}_{n-1}^{\upsilon }\right]$$E13

The noiseless counterpart of $X_{\upsilon }$ is designated as $X=\left[{\overline{x}}_0,{\overline{x}}_1,\dots ,{\overline{x}}_{n-1}\right]$. To insure the existence of enough samples in calculating the PCA transformation matrix, $n$ could not be too small. Practically speaking, it will be used in denoising at least $c\cdot m$ training samples of $x_{\upsilon }$, with $c=8˜10$. That is to say that in case of $n<c\cdot m$, we will use the best $c\cdot m$-matched samples in PCA training. Often, the best $c\cdot m$-matched samples are robust for estimating the local statistics of image, and this operation makes the algorithm more stable for computing the PCA transformation matrix. The problem now is how to estimate from the noisy data $X_{\upsilon }$, the noiseless dataset X. Once this dataset X is estimated, the central block and therefore we can extract the central underlying pixel. Such procedure is applied to each pixel, and then the entire image $I_{\upsilon }$ can be denoised. The LPG-PCA-based denoising is detailed in [25], and the denoising refinement in the second stage will be detailed in the next part of this paper.

3.3. Denoising refinement in the second stage

Most of the noise will be suppressed by employing the denoising procedures described in [25]. However, there is still much visually unpleasant residual noise in the denoised image. Figure 2 shows an example of image denoising where (a) is the original image Cameraman, (b) the noisy version of it with $\mathit{\text{PSNR}}=22.1\mathrm{dB}$$$ and $\sigma =20$, and $\left(a\right)$ is the denoised image with $\mathit{\text{PSNR}}=29.8\mathrm{dB}$ by employing the LPG-PCA technique proposed in [25]. Despite the remarkable improvement of $\mathit{\text{PSNR}}$, one can still see much residual noise in the denoising output. There are mainly two reasons for the residual noise. First, because of the strong noise in the original dataset $X_{\upsilon }$, the covariance matrix ${\mathrm{\Omega }}_{{\overline{x}}_{\upsilon }}$ is much noise degraded, which leads to estimation bias of the PCA transformation matrix and therefore deteriorates the denoising performance; second, the strong noise in the original dataset will also lead to LPG errors, which therefore results in estimation bias of the covariance matrix ${\mathrm{\Omega }}_{{\overline{x}}_{\upsilon }}$ or ${\mathrm{\Omega }}_{\overline{x}}$. Consequently, it is essential to further process the denoising output for a better image denoising. As the noise has been much canceled in the first round of LPG-PCA denoising, the LPG correctness and the estimation of ${\mathrm{\Omega }}_{{\overline{x}}_{\upsilon }}$ or ${\mathrm{\Omega }}_{\overline{x}}$ can be much ameliorated with the denoised image. Consequently, the LPG-PCA denoising procedure for the second round for enhancing the denoising results.

Figure 2.

(a) Original image Cameraman, (b) corresponding noisy image ($\mathit{\text{PSNR}}=22.1\mathrm{dB}$), (c) denoised image after the first round of the technique proposed in [25] ($\mathit{\text{PSNR}}=29.8\mathrm{dB}$), and (d) denoised image after the second round of the proposed technique ($\mathit{\text{PSNR}}=30.1\mathrm{dB}$) [25].

According to this figure, we remark that the visual quality is much ameliorated after the second round of refinement. As shown in Figure 3, in the second round of LPG-PCA denoising technique [25], the noise-level should be updated.

Figure 3.

Flowchart of the two-stage LPG-PCA denoising technique proposed in [25].

Denote by $\widehat{I}$, the denoised version of the noisy image in the first stage. The $\widehat{I}$ can be expressed as

$$\widehat{I}=I+{\upsilon }_s$$E14

where ${\upsilon }_s$ is the residual noise in the denoised image. The level estimation of ${\upsilon }_s$ is denoted by ${\sigma }_s=\sqrt{E\left[{\upsilon }_s^2\right]}$ and inputs it to the second round of LPG-PCA denoising algorithm. In [25], ${\sigma }_s$ is estimated based on the difference between $I_{\upsilon }$ and $\widehat{I}$. Let

$$\tilde{I}=I_{\upsilon }-\widehat{I}=\upsilon -{\upsilon }_s$$E15

We have $E\left[{\tilde{I}}^2\right]=E\left[{\upsilon }^2\right]+E\left[{\upsilon }_s^2\right]-2E\left[\upsilon \cdot {\upsilon }_s\right]={\sigma }^2+{\sigma }_s^2-2E\left[\upsilon \cdot {\upsilon }_s\right]$. The ${\upsilon }_s$ can be seen as the smoothed version of noise $\upsilon$, and it mainly contains the low-frequency component of $\upsilon$. Let $\tilde{\upsilon }=\upsilon -{\upsilon }_s$ be their difference, and $\tilde{\upsilon }$ mainly contains the high-frequency component of $\upsilon$. There is $E\left[\upsilon \cdot {\upsilon }_s\right]=E\left[\tilde{\upsilon }\cdot {\upsilon }_s\right]+E\left[{\upsilon }_s^2\right]$. Generally, compared to $E\left[{\upsilon }_s^2\right]$, $E\left[\tilde{\upsilon }\cdot {\upsilon }_s\right]$ is much smaller, and we can obtain the following approximation: $E\left[\upsilon \cdot {\upsilon }_s\right]\approx E\left[{\upsilon }_s^2\right]={\sigma }_s^2$. Thus, from $E\left[{\tilde{I}}^2\right]={\sigma }^2+{\sigma }_s^2-2E\left(\upsilon \cdot {\upsilon }_s\right)$, we obtain

$${\sigma }_s^2\approx {\sigma }^2-E\left[{\tilde{I}}^2\right]$$E16

In practice, ${\upsilon }_s$ will include not only the residual noise but also the estimation error of noiseless image $I$. Consequently, in the implementation [25], of Lei Zhang et al. let

$${\sigma }_s=C_s\sqrt{{\sigma }^2-E\left[{\tilde{I}}^2\right]}$$E17

where $C_s$ is a constant satisfying $C_s<1$. In [25], Lei Zhang et al. found experimentally that setting $C_s$ around $0.35$ can lead to satisfying denoising results for most of the testing images. Figure 2d shows the denoised image ($\mathit{\text{PSNR}}=30.1\mathit{dB}$) after the second round of the LPG-PCA denoising technique [25]. Although the $\mathit{\text{PSNR}}$ is not too much ameliorated on this image, we can remark clearly that the visual quality is much ameliorated by efficiently eliminating the residual noise obtained from the first round of denoising.

4. The proposed image denoising technique

As previously mentioned, in this chapter, a new image denoising technique is proposed. It combines two denoising approaches. The first one is a dual-tree discrete wavelet (DT-DWT)-based denoising method [12], and the second one is a two-stage image denoising by PCA with LPG [25]. This proposed technique consists at the first step in applying the first denoising approach [12] to the noisy image in order to obtain the first estimation of the clean image (the cleaned image). Then, we estimate the level of noise corrupting the clean image. The cleaned image, the noisy image, and the noise-level are used for applying the second approach which is two-stage image denoising by PCA with LPG [25]. Figure 4 illustrates the block diagram of the proposed technique.

Figure 4.

The block diagram of the proposed image denoising technique.

According to this figure, the first step of the proposed image denoising technique consists in applying the first denoising approach based on DT-DWT [12] to the noisy image, $I_b$, in order to obtain a first estimate of the clean image, $I_d$, and then estimates the noise-level, $\upsilon$, from $I_b$. The noisy images $I_b$, $I_d$, and $\upsilon$ constitute the inputs of the second image denoising system proposed in [25, 27]. The output of this system and the overall proposed one are the final denoised image, $I_{d1}$. In the image denoising system (LPG-PCA denoising) proposed in [25, 27], Lei Zhang et al. have used the clean image, $I$, and the noise-level, $\upsilon$, as the inputs of this system [27]. However, only the noisy image, $I_b$, is available, and for this raison, we have used in our proposed technique the denoising approach based on DT-DWT [12] in order to obtain a cleaned image, $I_d$, which is then used as a clean image, $I$. This clean image is one important input of the denoising system proposed by Lei Zhang et al. [27]. In the following two subsections, we will be interested in the first image denoising approach based on DT-DWT [12] and the technique of noise-level estimation proposed in [28, 29], from the noisy image, $I_b$.

5. The Hilbert transform

The Hilbert transform of a signal corresponds in Fourier plane to a filter with complex gain, $-i\mathit{sign}\left(\gamma \right)$ [30]. This is corresponding to an impulse response $\mathit{vp}\left(\frac{1}{\pi t}\right)$ where $\mathit{vp}$ is the principal value in Cauchy sense [30]. The analytic signal is then constructed as follows:

$$z\left(t\right)=x\left(t\right)+\mathit{iH}\left\{x\left(t\right)\right\}=x\left(t\right)+\frac{i}{\pi }\mathit{vp}\underset{-\infty }{\overset{+\infty }{\int }}\frac{x\left(s\right)}{t-s}\mathit{ds}$$E18

This analytic signal has only positive frequencies. The Hilbert transform of a real signal is also real. Instead of considering the Hilbert transform of the wavelet (which is defined through the associated filters), we can consider the Hilbert transform of the signal, and the analysis is performed with initial wavelet because we have $⟨f,H{\psi }_{a,t}⟩=⟨Hf,{\psi }_{a,t}⟩$ [30]. The latter equality is justified by the fact that the Hilbert transform is considered as a linear filter [30]. Therefore, we have the following scheme: let $X\left(n\right)$ be the signal to be analyzed with real wavelet by using the Mallat algorithm in order to obtain the wavelet coefficients, $d_1\left(j,k\right)$. Then, we analyze $\mathit{HX}\left(n\right)$ with the same wavelet, and we obtain the coefficients $d_2\left(j,k\right)$. Then, we construct the complex coefficients: $d_{\mathit{\text{complex}}}=d_1\left(j,k\right)+i{d}_2\left(j,k\right)$. As follows, the magnitude of those coefficients is named Hilbert magnitude. The drawbacks of this method are as follows: The support of the Hilbert transform of a wavelet having a compact support is infinite. There is a computing disadvantage because the cost of two wavelet transforms plus the Hilbert Transform. Theoretically speaking, it is possible to limit the drawback of the support of the Hilbert transform of the wavelet by using an approximate of the Hilbert transform. However, this approximation cannot be optimized for all scales [30]. One solution of this problem has been proposed by Kingsbury: the dual tree [30].

6. Dual-tree complex wavelet transform

The dual tree complex wavelet (DT-CWT) permits to make signal analysis by using two different trees of DWT, with filters selected in such manner to obtain approximately a signal decomposition using analytic wavelet [30]. Figure 5 shows a tree of DT-CWT, using two different filter banks: h₁ and g₁ are high-pass filters of the first and second trees, and h₀ and g₀ are low-pass filters of the same two trees [30]. The first tree gives the coefficients of the real part, d_r(j,k), and the second tree gives those of the imaginary part, d_i(j,k). After that, we construct the complex coefficients d_complex(j,k) = d_r(j,k) + i d_i(j,k). The magnitude of those coefficients is named dual-tree magnitude [30].

Figure 5.

Dual tree of real filters for the Q-shift DT-CWT, giving real and imaginary parts of complex coefficients from tree (a) and tree (b), respectively [30].

This Q-shift dual-tree complex wavelet transform (Figure 5) is in 1D. Synthesis of the filters adapted to this structure has been performed by many research works. Particularly, Kingsbury [30] proposed some filters named Q-shift. In [30], some filters are employed, and their utilization is equivalent to the signal analysis by wavelets illustrated in Figure 6.

Figure 6.

Q-shift wavelet obtained with filters Antonini [30].

We can see in this figure that the wavelet corresponding to the imaginary part tree is very near to the Hilbert transform of the wavelet corresponding to the real part tree [30]. Finally, the utilization of this structure requires an operation of prefiltering; it means that the filters used in the first step are not the same as those used in the next step. The advantages of this method compared to the simple Hilbert transform (Section 5) are [30]:

• A lower computation cost ( Two DWT),

• An approximate of the Hilbert transform, is optimized for each scale,

• The possibility of an exact reconstruction is preserved.

The principal drawback of the DT-CWT is the non-possibility of the use of the well-known wavelets of the DWT (Daubechies wavelet, Spline, etc.) and therefore the non-possibility to choose the number of vanishing moments (all the Q-shift filter gives wavelets with two vanishing moments).

6.1. 2D DT-CWT

To explain how the DT-CWT produces oriented wavelets, consider the 2D wavelet $\psi \left(x,y\right)=\psi \left(x\right)\cdot \psi \left(y\right)$ associated with the row-column implementation of the wavelet transform, where $\psi \left(x\right)$ is a complex wavelet (approximately analytic) and is expressed as follows [31]:

$$\psi \left(x\right)={\psi }_h\left(x\right)+i{\psi }_g\left(x\right).$$E19

Therefore, we obtain the following expression of $\psi \left(x,y\right)$:

$$\begin{array}{c}\psi \left(x,y\right)=\left[{\psi }_h\left(x\right)+i{\psi }_g\left(x\right)\right]\left[{\psi }_h\left(y\right)+i{\psi }_g\left(y\right)\right]={\psi }_h\left(x\right){\psi }_h\left(y\right)-{\psi }_g\left(x\right){\psi }_g\left(y\right)+\\ i\left[{\psi }_g\left(x\right){\psi }_h\left(y\right)+{\psi }_h\left(x\right){\psi }_g\left(y\right)\right]\end{array}$$E20

The following idealized diagram (Figure 7) illustrates the Fourier spectrum support of this complex wavelet [31].

Figure 7.

Idealized diagram illustrating the Fourier spectrum support of the complex wavelet, $\psi \left(x,y\right)=\psi \left(x\right)\cdot \psi \left(y\right)$ [31].

Since the (approximately) 1D wavelet spectrum is supported on just one side of the frequency axis, the complex 2D wavelet ($\psi \left(x,y\right)$) spectrum is supported in just one quadrant of the 2D frequency plane. That is why the complex 2D wavelet, $\psi \left(x,y\right)$, is oriented. If the real part of this complex wavelet is taken, then the sum of two separable wavelets is obtained:

$$\text{Re}\mathit{al}\mathit{\text{Part}}\left\{\psi \left(x,y\right)\right\}={\psi }_h\left(x\right){\psi }_h\left(y\right)-{\psi }_g\left(x\right){\psi }_g\left(y\right).$$E21

Since the real function spectrum should be symmetric with respect to the origin, then the spectrum of this real wavelet is supported in two quadrants of the 2D frequency plane (Figure 8).

Figure 8.

Idealized diagram illustrating the support of spectrum of this real wavelet, $\text{Re}\mathit{al}\mathit{\text{Part}}\left\{\psi \left(x,y\right)\right\}$ [31].

Unlike the real separable wavelet, the support of the spectrum of this real wavelet has not the checkerboard artefact and consequently this real wavelet (illustrated in the second panel of Figure 11), is oriented at $-{45}^{\circ }$. It is deserving mentioning that this construction is depending on $\psi \left(x\right)={\psi }_h\left(x\right)+i{\psi }_g\left(x\right)$ being (approximately) analytic or equivalently on ${\psi }_g\left(x\right)$ being approximately the Hilbert transform of ${\psi }_h\left(x\right)$ (${\psi }_g\left(x\right)\approx H\left\{{\psi }_h\left(x\right)\right\}$). Note that ${\psi }_h\left(x\right){\psi }_h\left(y\right)$ is the sub-band HH of a separable 2D real wavelet transform implemented employing the filters $\left\{h{}_0\left(n\right),h_1\left(n\right)\right\}$. The term ${\psi }_g\left(x\right){\psi }_g\left(y\right)$ is also the sub-band HH which is obtained from the application of a real separable wavelet transform. The latter is implemented by employing the filters $\left\{g{}_0\left(n\right),g_1\left(n\right)\right\}$. To have a real 2D wavelet oriented at $+{45}^{\circ }$, we consider now the complex 2D wavelet ${\psi }_2\left(x,y\right)=\psi \left(x\right)\overline{\psi }\left(y\right)$ where $\overline{\psi }\left(y\right)$ is the complex conjugate of $\psi \left(y\right)$ and, as previously mentioned, $\psi \left(x\right)$ is approximately the analytic wavelet, ${\psi }_h\left(x\right)+i{\psi }_g\left(x\right)$. Therefore, we have

$$\begin{array}{l}{\psi }_2\left(x,y\right)=\left[{\psi }_h\left(x\right)+i{\psi }_g\left(x\right)\right]\left[{\psi }_h\left(y\right)-i{\psi }_g\left(y\right)\right]\\ ={\psi }_h\left(x\right){\psi }_h\left(y\right)+{\psi }_g\left(x\right){\psi }_g\left(y\right)+i\left[{\psi }_g\left(x\right){\psi }_h\left(y\right)+{\psi }_h\left(x\right){\psi }_g\left(y\right)\right]\end{array}$$E22

The support in the 2D frequency plane of this complex wavelet spectrum is illustrated in Figure 9.

Figure 9.

The idealized diagram in 2D frequency plane of the spectrum of this complex wavelet [31].

As above, the spectrum of the complex wavelet, ${\psi }_2\left(x,y\right)$, is supported in just one quadrant of the 2D frequency plane. If the real part of this complex wavelet is taken, then we have

$$\text{Re}\mathit{al}\mathit{\text{Part}}\left\{{\psi }_2\left(x,y\right)\right\}={\psi }_h\left(x\right){\psi }_h\left(y\right)+{\psi }_g\left(x\right){\psi }_g\left(y\right).$$E23

The spectrum of which is supported in two quadrants of the 2D frequency plane as illustrated in Figure 10.

Figure 10.

Idealized diagram in 2D frequency plane of the spectrum $\text{Re}\mathit{al}\mathit{\text{Part}}\left\{{\psi }_2\left(x,y\right)\right\}$ [31].

Again, neither the wavelet nor the spectrum of this real wavelet has the spectrum of the checkerboard artifact. This real 2D wavelet is oriented at $+{45}^{\circ }$ as illustrated in the fifth panel of Figure 11. To have four more oriented real 2D wavelets, one can repeat this procedure on the complex wavelets expressed as follows: $\varphi \left(x\right)\psi \left(y\right)$, $\psi \left(x\right)\varphi \left(y\right)$, $\varphi \left(x\right)\stackrel{\rightharpoonup }{\psi }\left(y\right)$, and $\psi \left(x\right)\overline{\varphi }\left(y\right)$ where we have

$$\psi \left(x\right)={\psi }_h\left(x\right)+i{\psi }_g\left(x\right)$$E24

$$\varphi \left(x\right)={\varphi }_h\left(x\right)+i{\varphi }_g\left(x\right)$$E25

Figure 11.

Typical wavelets associated with the real oriented 2D dual-tree wavelet transform. Top row illustrates the wavelets in the space domain: bottom row illustrates the (idealized) support of the Fourier spectrum of each wavelet in 2D frequency plane. The absence of the checkerboard phenomenon is observed in both frequency and spatial domains.

By taking the real part of each of these wavelets, one can obtain four real oriented 2D wavelets. Moreover, the two already obtained in Eqs. (21) and (23). Precisely, we have six wavelets expressed as follows:

$${\psi }_i\left(x,y\right)=\frac{1}{\sqrt{2}}\left({\psi }_{1,i}\left(x,y\right)-{\psi }_{2,i}\left(x,y\right)\right)$$E26

$${\psi }_{i+3}\left(x,y\right)=\frac{1}{\sqrt{2}}\left({\psi }_{1,i}\left(x,y\right)+{\psi }_{2,i}\left(x,y\right)\right)$$E27

For $i=1,2,3,$ the two separable 2-D wavelet bases are expressed as follow:

$${\psi }_{1,1}\left(x,y\right)={\varphi }_h\left(x\right)\psi {}_h\left(y\right),{\psi }_{2,1}\left(x,y\right)={\varphi }_g\left(x\right)\psi {}_g\left(y\right),$$E28

$${\psi }_{1,2}\left(x,y\right)={\psi }_h\left(x\right)\varphi {}_h\left(y\right),{\psi }_{2,2}\left(x,y\right)={\psi }_g\left(x\right)\varphi {}_g\left(y\right),$$E29

$${\psi }_{1,3}\left(x,y\right)={\psi }_h\left(x\right)\psi {}_h\left(y\right),{\psi }_{2,3}\left(x,y\right)={\psi }_g\left(x\right)\psi {}_g\left(y\right),$$E30

The normalization factor $1/\sqrt{2}$ is used only so that the sum/difference operation constitutes an orthonormal operation. In Figure 11 the six real oriented wavelets derived from a pair of typical wavelets satisfying ${\psi }_g\left(x\right)\approx H\left\{{\psi }_h\left(x\right)\right\}$ are illustrated. Compared to separable wavelets, these six non-separable wavelets succeed in isolating different orientations. Each of these six wavelets are aligned with a specific direction. Moreover, no checkerboard effect appears. In addition, they cover more distinct orientations than the separable wavelets obtained from the application of DWT. Moreover, since the sum/difference operation is orthonormal, the wavelet set is obtained from integer translates and dyadic dilations from a frame [31].

7. The technique of Noise-level estimation

In many image processing applications, the noise-level is an important parameter. For example, the performance of an image denoising technique can be much degraded due to the poor noise-level estimation. The most available denoising techniques simply supposed that the noise-level is known that largely prevents them from practical employment. Furthermore, even with the given true noise-level, those denoising techniques still cannot achieve the best performance, precisely for scenes with rich texture. Xinhao Liu et al. [28, 29] have proposed a technique of patch-based noise-level estimation, and they suggested that the noise-level parameter should be tuned according to the complexity of the scene. Their approach [28, 29] includes the process of selecting low-rank patches without high-frequency components from a single noisy image. Then, the noise-level was estimated from the selected patches employing principal component analysis. Because the exact noise-level does not always provide the best performance for non-blind denoising. Experiments prove that both the stability and precision are superior to the state-of-the-art noise-level estimation technique for different noise-levels and scenes.

8. Evaluation criteria

In this section, we will evaluate the three techniques which are the proposed image denoising techniques: the first image denoising approach based on DT-CWT [12] and the second denoising approach and the two-stage image denoising by principal component analysis with local pixel grouping [25]. This evaluation is based on the computation of PSNR and SSIM which are detailed in [32].

9. Results and discussion

In this work, we have applied the proposed image denoising technique, the first image denoising technique based on DT-CWT [12] and the second denoising technique and the two-stage image denoising by principal component analysis with local pixel grouping [25], on a number of digital images such as “House,” “Lena,” and “Cameraman.” These images are degraded by additive white noise with different values of noise-level, $\sigma$. PSNR and SSIM values obtained from the application of the three mentioned techniques on the noisy images are listed in Table 1.

Technique	The first image denoising technique based on DT-DWT [12]	Two-stage image denoising by principal component analysis with local pixel grouping [25]: first stage	Two-stage image denoising by principal component analysis with local pixel grouping [25]: second stage	The proposed technique
House (σ = 10)	34.7138 (0.8778)	35.4 (0.9003)	35.6 (0.9012)	36.1223 (0.9130)
House (σ = 20)	31.6671 (0.8253)	31.8 (0.8084)	32.5 (0.8471)	33.0828 (0.8677)
House (σ = 30)	29.8494 (0.7877)	29.3 (0.7225)	30.4 (0.8185)	31.2095 (0.8393)
House (σ = 40)	28.5744 (0.8084)	27.3(0.6243)	28.9 (0.7902)	29.7344 (0.8084)
Lena (σ = 10)	33.6767(0.9170)	33.6 (0.9218)	33.7 (0.9243)	34.0765 (0.9271)
Lena (σ = 20)	30.0002 (0.8539)	29.5 (0.8346)	29.7 (0.8605)	30.5415 (0.8765)
Lena (σ = 30)	27.9859 (0.8016)	27.1 (0.7441)	27.6 (0.8066)	28.3595 (0.8292)
Lena (σ = 40)	26.6364 (0.7585)	25.4 (0.6597)	26.0 (0.7578)	26.8566 (0.7882)
Cameraman (σ = 10)	32.7481 (0.8989)	33.9 (0.9261)	34.1 (0.9356)	33.6141 (0.9241)
Cameraman (σ = 20)	28.9990 (0.8175)	29.8 (0.8320)	30.1 (0.8902)	29.7184 (0.8575)
Cameraman (σ = 30)	27.1022 (0.7641)	27.3 (0.7395)	27.8(0.8558)	27.8174 (0.8151)
Cameraman (σ = 40)	25.7866 (0.7241)	25.5 (0.6393)	26.2 (0.8211)	26.4954 (0.7826)
Monarch (σ = 10)	32.9907 (0.9369)	34.0 (0.9522)	34.2 (0.9594)	34.0698 (0.9553)
Monarch (σ = 20)	29.1114 (0.8811)	29.6 (0.8859)	30.0 (0.9202)	30.0384 (0.9145)
Monarch (σ = 30)	27.0058 (0.8346)	27.0 (0.8071)	27.4 (0.8769)	27.7209 (0.8735)
Monarch (σ = 40)	25.5973 (0.7950)	25.2 (0.7267)	25.9 (0.8378)	26.0832 (0.8293)
Peppers (σ = 10)	33.4942 (0.9056)	33.4 (0.8909)	33.3 (0.8943)	33.7904 (0.9189)
Peppers (σ = 20)	29.8124 (0.8424)	29.9 (0.8177)	30.1 (0.8413)	30.5252 (0.8743)
Peppers (σ = 30)	27.7810 (0.7924)	27.5 (0.7332)	27.9 (0.7973)	28.4765 (0.8356)
Peppers (σ = 40)	26.4045 (0.7507)	25.9 (0.6447)	26.7(0.7648)	26.9883 (0.8013)
Paint (σ = 10)	32.5488 (0.9165)	33.5 (0.9280)	33.6 (0.9311)	33.3567 (0.9276)
Paint (σ = 20)	28.5980 (0.8416)	26.8 (0.7467)	29.5 (0.8683)	29.4699 (0.8648)
Paint (σ = 30)	26.6067 (0.7817)	26.8 (0.7467)	27.2 (0.8088)	27.2540 (0.8077)
Paint (σ = 40)	25.2968 (0.7330)	25.0 (0.6590)	25.6 (0.7569)	25.6389 (0.7560)

Table 1.

PSNR (dB) and SSIM results of the denoised images for the different techniques.

These obtained results (Table 1) show clearly that the proposed technique outperforms the denoising technique based on DT-CWT proposed in [12] and the denoising approach based on LPG-PCA [25]. Figures 12–15 show four examples of image denoising using the proposed technique.

Figure 12.

(a) Clean image (Cameraman.tif), (b) Noisy image with, (c) The denoised image by the proposed technique (the first stage) and denoised image by the proposed technique (the second stage).

Figure 13.

(a) Clean image (Monarch.tif), (b) Noisy image with, (c) The denoised image by the proposed technique (the first stage) and denoised image by the proposed technique (the second stage).

Figure 14.

(a) Clean image (Lena.tif), (b) Noisy image with, (c) The denoised image by the proposed technique (the first stage) and denoised image by the proposed technique (the second stage).

Figure 15.

(a) Clean image (Peppers.tif), (b) Noisy image with, (c) The denoised image by the proposed technique (the first stage) and denoised image by the proposed technique (the second stage).

These figures show that the noise corrupting the original images is sufficiently suppressed. Moreover, the proposed technique permits to obtain denoised images with good perceptual quality. In each of these figures, the image (c) is obtained after the first denoising stage in the proposed technique. In this image (c), some noise is still existing, whereas it is considerably reduced into the image (d) obtained after the second denoising step. In the following subsection, we will give the results obtained by applying the proposed technique, the LPG-PCA-based denoising technique [25, 27] and the DT-DWT-based denoising one to a number of grayscale images. Those results are in terms of SNR and MSE and are listed in Table 2.

Technique	The first image denoising technique based on DT-DWT [12]	Two-stage image denoising by principal component analysis with local pixel grouping [25]	The proposed technique
House (σ = 10)	SNR = 78.00, MSE = 21.96	SNR = 79.41, MSE = 15.88	SNR = 79.44, MSE = 15.75
House (σ = 20)	SNR = 74.95, MSE = 44.29	SNR = 76.37, MSE = 31.97	SNR = 76.38, MSE = 31.92
House (σ = 30)	SNR = 73.14, MSE = 67.31	SNR = 74.50, MSE = 49.21	SNR = 74.50, MSE = 49.21
House (σ = 40)	SNR = 71.86, MSE = 90.28	SNR = 73.02, MSE = 69.16	SNR = 73.02, MSE = 69.14
Lena (σ = 10)	SNR = 74.67, MSE = 27.88	SNR = 75.28, MSE = 24.17	SNR = 75.28, MSE = 24.19
Lena (σ = 20)	SNR = 70.99, MSE = 65.02	SNR = 71.53, MSE = 57.40	SNR = 71.55, MSE = 57.19
Lena (σ = 30)	SNR = 68.97, MSE = 103.39	SNR = 69.35, MSE = 94.87	SNR = 69.37, MSE = 94.36
Lena (σ = 40)	SNR = 67.62, MSE = 141.07	SNR = 67.85, MSE = 134.09	SNR = 67.87, MSE = 133.35
Cameraman (σ = 10)	SNR = 75.33, MSE = 34.53	SNR = 76.19, MSE = 28.29	SNR = 76.23, MSE = 28.06
Cameraman (σ = 20)	SNR = 71.58, MSE = 81.88	SNR = 72.30, MSE = 69.38	SNR = 72.33, MSE = 68.80
Cameraman (σ = 30)	SNR = 69.68, MSE = 126.72	SNR = 70.39, MSE = 107.48	SNR = 70.45, MSE = 106.00
Cameraman (σ = 40)	SNR = 68.36, MSE = 171.56	SNR = 69.07, MSE = 145.72	SNR = 69.14, MSE = 143.51
Monarch (σ = 10)	SNR = 74.94, MSE = 32.65	SNR = 76.02, MSE = 25.47	SNR = 76.01, MSE = 25.55
Monarch (σ = 20)	SNR = 71.06, MSE = 79.78	SNR = 71.99, MSE = 64.45	SNR = 71.98, MSE = 64.53
Monarch (σ = 30)	SNR = 68.96, MSE = 129.56	SNR = 69.67, MSE = 109.89	SNR = 69.68, MSE = 109.62
Monarch (σ = 40)	SNR = 67.55, MSE = 179.20	SNR = 68.01, MSE = 161.05	SNR = 68.03, MSE = 160.25
Peppers (σ = 10)	SNR = 76.07, MSE = 29.08	SNR = 76.65, MSE = 25.43	SNR = 76.63, MSE = 25.56
Peppers (σ = 20)	SNR = 72.39, MSE = 67.89	SNR = 73.10, MSE = 57.61	SNR = 73.12, MSE = 57.43
Peppers (σ = 30)	SNR = 70.36, MSE = 108.38	SNR = 71.05, MSE = 92.34	SNR = 71.07, MSE = 92.02
Peppers (σ = 40)	SNR = 68.98, MSE = 148.80	SNR = 69.57, MSE = 130.09	SNR = 69.58, MSE = 129.58

Table 2.

SNR (dB) and MSE results of the denoised images for the different techniques.

Those results show that the proposed technique outperforms the two other techniques (the LPG-PCA-based denoising technique [25, 27] and the DT-DWT-based denoising one [12]). In fact the proposed techniques are the highest values of SNR and lowest values of MSE.

10. Conclusion

In this chapter, a new image denoising technique is proposed. It combines two denoising approaches. The first one is a dual-tree discrete wavelet transform (DT-DWT)-based denoising technique, and the second one is a two-stage image denoising by principal component analysis with local pixel grouping (LPG-PCA). The first step of this proposed technique consists in applying the first approach to the noisy image in order to obtain a first estimate of the clean image. Then, we estimate the level of noise corrupting the original image. This estimation is performed by using a method of noise estimation from noisy images. The third step of the proposed technique consists in using this first clean image estimation, the noisy image, and this noise-level estimate as inputs of the second image denoising system (LPG-PCA-based image denoising) in order to obtain the final estimation of the clean image. A comparative study is performed between the proposed image denoising technique and two others denoising approaches where the first is based on DT-DWT and the second is based on LPG-PCA. This study is based on PSNR and SSIM computations, and the obtained results show that the proposed technique outperforms the two other denoising approaches. We also computed SNR (Signal to Noise Ratio) and MSE (Mean Square Error) and the obtained results also show that the proposed technique outperforms the others techniques.

Acknowledgments

We would like to thank all the people who contributed in some way to this work which was supported by the CRTEn (Center of Research and Technology of Energy) of Borj Cedria, Tunisia, and the Ministry of Higher Education and Scientific Research.