Adaptive Hybrid Bivariate Double Density Discrete and Complex Wavelet for Image Denoising

Gamal Fahmy^1*, Mamdouh F. Fahmy², Omar Fahmy³

¹Faculty of Computer Science and Engineering, Galala University, Al Galala, Egypt.
²Electrical Engineering Department, Assiut University, Assiut, Egypt.
³Electrical Engineering Department, Badr University, Cairo, Egypt.

DOI: 10.4236/jcc.2023.112004                

Abstract

Image denoising is an important step in eliminating any noise impact in any image transmission process. Recently we presented two approaches for Bivariate based image denoising. They were Double Density Discrete Wavelet Transform (DD DWT) and Double Density Dual Tree Complex Wavelet Transform (DD CWT). In both techniques we decomposed noisy images with either DD DWT or DD CWT decompositions and then applied the Bivariate based denoising technique for noise removal. In this paper we propose an adaptive hybrid technique for Bivariate based image denoising that is based on the synthesis of DD-DWT bands or DD-CWT bands but with different weights, to deliver enhanced image features with less denoising impact especially around image edges, which is the most effected by noisy transmission channels. This proposed technique has been also enhanced by edge sharpening and Eigen analysis, as two separate stages. Simulation result comparisons have been performed between the proposed hybrid band adaptive DD-DWT and DD-CWT technique and the two primary techniques DD-DWT, DD- CWT, as well as other superior literature techniques such the original bivariate denoising technique with both original Complex Wavelet Transform and Double Density decompositions. This work in specific compares between Double Density DWT and Double Density CWT decompositions, proposes new filter design that suits each of them and proposes a hybrid technique between as will be shown.

Keywords

Complex Wavelet, Image Deniosing, Double Density Wavelets

1. Introduction

There has always been an amount of unavoidable noise contaminated in any communication channel. Eliminating or avoiding this noise amount is inevitable. Hence, denoising became crucial to improve quality and eliminate degradations of received data in any receiver system.

Several image denoising algorithms were proposed during the last two decades. They range from the frequency-based filtering techniques [1] [2] , and the wavelet transform-based techniques [3] . This is in addition to other transform-based techniques that were not competitive enough.

Typically, in wavelet-based techniques, the signal/image is represented in frequency domain as a few numbers of large coefficients as well as large number of almost zero coefficients. The adopted denoising strategy is based on thresholding small coefficients, where the noise effect is high, while keeping or modifying large coefficients, where the noise effect is low. In [4] [5] , a complete review for image denoising was presented. It concluded that superior performance for denoising can be achieved with wavelet based denoising, rather than other transform based denoising techniques or convex optimization based denoising approaches.

A Bivariate Shrinkage technique was first adopted in wavelet image denoising in [6] [7] , where the correlation between parent and children was exploited for better noise removal performance. This was performed by utilizing the Maximum Likelihood Estimation between noisy wavelet coefficients and their adjacent at the next coarser scale. Complex Wavelet Transform achieved better denoising performance for this Bivariate technique as in [7] . In [8] , the double density discrete wavelet transform filter bank structure was first presented for denoising purposes. Performance for image denoising was outstanding due to the fact that double density wavelet bands were shift invariant and phase oriented. Denoising was also performed by simple thresholding of double density bands. The Discrete Wavelet Transform DWT has been historically is known for being an efficient time-frequency representation/analysis, however, it has the following drawbacks: 1) Shift Variance; 2) Oscillations at edges and sharp corners. 3) Aliasing; and 4) Lack of orientation identifications [9] . The Complex Wavelet Transform CWT [9] [10] [11] [12] , was presented to overcome these short comings as it provided wavelet functions that have almost single side band frequency response. The one dimensional CWT, 1-D CWT, was first presented in [10] [11] [12] [13] as two real DWT trees connected in parallel. A complete design analysis and review of the filter sets of CWT has been presented in [9] [14] [15] . The filter sets also ensured that the upper and lower trees form a Hilbert transform pair approximately in all scales. The two-dimensional CWT, 2-D CWT is implemented as four 2-D DWT connected in parallel. The Bivariate technique was later adopted in image denoising through utilizing CWT, as in [7] . One dimensional Double Density DWT, 1-D DD-DWT, was implemented using one low pass filter $H_0\left(z\right)$ and two high pass filters $H_1\left(z\right)$ & $H_2\left(z\right)$ as shown in Figure 1, [8] . Similarly, 2-D DD-DWT is represented as two 1-D DD-DWT in raw and column manners, respectively. The DD-DWT still suffered from directional ambiguity; hence the double density complex wavelet transform that is briefly presented in [9] , has been recently proposed to overcome this draw back. In [16] , a similar analysis was presented however it did not propose any filter design procedure as presented in this work. In this paper, we first propose a novel two-dimensional double density Discrete Wavelet Transform 2-D DD-DWT filter structure for denoising purposes, then we propose another two-dimensional double density Complex Wavelet Transform 2-D DD-CWT filter structure and utilize it in denoising purposes. These filters in both techniques satisfy the perfect reconstruction property as well as the alias free condition. Next, we utilize these filters in conjunction with the Bivariate Shrinkage based denoising technique. We then propose our unique hybrid technique that combines both 2-D techniques and produces an enhanced denoising technique as will be shown in our simulations. We propose two scenarios to enhance this hybrid technique first through edge sharpening and then Eigen analysis manipulation. Simulation comparisons show that these two scenarios can achieve sizeable enhancement of the final denoised image. Some early results of this work was presented in [9] , but at much less scale and scope. The proposed 2-D DD-DWT filter structure is presented in section II.A, then the enhanced 2-D DD-CWT filter structure is presented in section II.B. The bivariate shrinkage is presented in section III, and our final hybrid bivariate technique that combines DD-DWT and DD-CWT is presented in section IV, along with two scenarios that are proposed for further enhancement. Simulation comparisons for all proposed techniques are illustrated in section V.

2. The Proposed DD DWT CWT Design

2.1. DD DWT

Double Density DWT wavelet

The one-dimensional double-density 1-D DD-DWT, first presented in [8] , is implemented by recursively applying a three-channel analysis filter bank, instead of the historic two-channel analysis filter bank, Figure 1. Also, the Inverse double-density DD DWT is obtained by applying the synthesis filter bank for the initial three channel analysis set. The 2-D double-density DD-DWT is consequently implemented by applying the filters $H_0\left(z\right)$ & $H_1\left(z\right)$ & $H_2\left(z\right)$ first to the rows, then to the columns of an image. This would result in nine 2-D sub bands, where one of them is the 2-D lowpass scaling filter, and the other eight make up remaining 2-D wavelet filters, Figure 2.

The proposed design to achieve a 3-channel perfect reconstruction filter bank set is detailed as follows: Let $H_0\left(z\right)$ be an N-tap low pass filter with K zeros at $z=-1$ . These K zeros are placed to insure smoothness. We then would want to construct two N-tap high pass filters $H_1\left(z\right)$ & $H_2\left(z\right)$ that satisfy the following perfect reconstruction (PR) and Alias Free (AF) conditions:

Figure 1. Double density discrete wavelet transform.

Figure 2. Two dimensional 2-D DD-DWT.

PR (condition):

$H_0\left(z\right)H_0\left(z^{-1}\right)+H_1\left(z\right)H_1\left(z^{-1}\right)+H_2\left(z\right)H_2\left(z^{-1}\right)=C{z}^{-\left(N-1\right)}$ , (1)

where C is a constant

AC (condition):

$H_0\left(-z\right)H_0\left(z^{-1}\right)+H_1\left(-z\right)H_1\left(z^{-1}\right)+H_2\left(-z\right)H_2\left(z^{-1}\right)=0$ (2)

We propose to satisfy these two conditions in two steps:

In step 1: For a given $H_0\left(z\right)$ , assume an arbitrary $H_1\left(z\right)$ , then scale $H_1\left(z\right)$ such that its norm is less than $\sqrt[2]{2-{\left|H_0\left(z\right)\right|}^2}$ . We then construct $H_2\left(z\right)$ to meet AC condition using a root finding technique as follows:

Let

$Y\left(z\right)=H_0\left(-z\right)H_0\left(z^{-1}\right)+H_1\left(-z\right)H_1\left(z^{-1}\right)=-Y_a\left(-z\right)Y_a\left(z^{-1}\right)$ . (3)

This means that $Y\left(z\right)$ has roots at $z_k$ & ( $-\frac{1}{z_k}$ ), $k=1,2,\cdots ,N-1$ .

Let

$H_2\left(z\right)=z^{-\left(N-1\right)}{Y}_a\left(-z^{-1}\right)$ .

where $Y_a\left(z\right)$ is constructed through root grouping, either by maximum phase, minimum phase or mixed phase factorization. This would satisfy AC condition.

In step 2: we update $H_1\left(z\right)$ as follows,

$\begin{array}{l}{H}_{1u}\left(z\right)\ast H_{1u}\left(z^{-1}\right)=Y_1\left(z\right)=2\\ =C{z}^{-\left(N-1\right)}-H_0\left(z\right)H_0\left(z^{-1}\right)-H_2\left(z\right)H_2\left(z^{-1}\right)\end{array}$

This would allow us to obtain $H_{1u}\left(z\right)$ using a root finding technique of $Y_1\left(z\right)$ again either using maximum, minimum or mixed phase factorization. Then we update $H_1\left(z\right)=H_{1u}\left(z\right)$ and reiterate step 1 until conversion is obtained. It is worth mentioning that, while the filter $H_0\left(z\right)$ has K zeros at $z=-1$ ; Both $H_1\left(z\right)$ & $H_2\left(z\right)$ must has M zeros at $z=1$ where $M\le K$ . Finally, because of the factorization of $Y\left(z\right)$ , the construction of $H_2\left(z\right)$ and subsequently $H_1\left(z\right)$ is non-unique.

Example: In this example, we choose a 6-tap minimal-length Debauches low-pass filter $H_0\left(z\right)$ with $K=4$ . In our simulation, we choose $M=2$ to construct an arbitrary $H_1\left(z\right)$ , length 4. We express $H_1\left(z\right)={\left(1-z^{-1}\right)}^M{X}_1\left(z\right)$ and similarly for $H_2\left(z\right)$ . The unknown parameters $X_1\left(z\right)X_2\left(z\right)$ proceeds as described above. For this arbitrary $H_1\left(z\right)$ , $Y\left(z\right)$ is factorized as in Equation (3). Only six iteration steps are needed to obtain the exact solution. Design coefficients of the upper tree $H_0,H_1,H_2$ are listed below in Table 1.

We note here that the 2-D DD-DWT still suffers from phase ambiguities that is typical in any basic Discrete wavelet Transform DWT. The next utilized Complex Wavelet Transform CWT has been proposed to overcome this phase ambiguity issue.

2.2. DD CWT

Double Density CWT Wavelet

In this section we present our novel 2-D DD-CWT filter bank structure. We first note that the design of the regular CWT filter structure bank $H_{0u},H_{0d}$ , which is the low pass filter of the upper and lower (down) trees must be a half band filter as noted in [9] [12] [13] . This is in addition to another restriction of a group delay. The group delay of $H_{0d}\left(z\right)/H_{0u}\left(z\right)$ must approximate one sample for the first level of decomposition and half sample for the succeeding levels [9] [12] [13] . Figure 3 shows the CWT system, where filters $H_{0u},H_{0d}$ are designed as $h_0,g_0$ . Further designed examples are listed in [8] [10] .

It has been shown that to fulfill this Hilbert pair relation; the low-pass filters $\left\{h_{0,u}\left(n\right),h_{0,d}\left(n\right)\right\}$ of the upper and lower tree satisfy the following constraint

$h_{0,d}\left(n\right)=h_{0,u}\left(n-0.5\right)\Rightarrow {\psi }_d\left(t\right)=\mathcal{H}\left\{{\psi }_u\left(t\right)\right\}$

where ${\psi }_u\left(t\right)$ & ${\psi }_d\left(t\right)$ are the wavelet functions of the upper and lower trees, respectively and $\mathcal{H}$ denotes Hilbert transform. Moreover, this condition also implies ${\varphi }_d\left(n\right)={\varphi }_u\left(n-0.5\right)$ where $\varphi \left(n\right)$ represents the corresponding discrete scaling functions. This half sample shift implies that in case of multi-level decomposition, integer translates of ${\varphi }_d\left(n\right)$ should lie midway between integer translates ${\varphi }_u\left(n\right)$ . Thus in order to ensure that the upper and lower wavelet functions form a Hilbert pair at every decomposition level $1,\cdots ,r$ , the first stage and the succeeding stages filters of the upper and lower should be chosen to satisfy the condition [4]

$h_{0,d}^1\left(n\right)=h_{0,u}^1\left(n-1\right),$

$h_{0,d}^j\left(n\right)=h_{0,u}^j\left(n-0.5\right),j=2,3,\cdots ,r$ (4)

where r is no. of decomposition levels and j is the stage index, as the lattice tree structure in Figure 4.

Table 1. Design coefficients of proposed DD-DWT.

Figure 3. CWT Filter bank structure.

Figure 4. CWT Filter bank structure-one tree.

In [4] [5] [6] , two techniques were described to fulfill this half sample delay requirements, namely the maximum flat delay [6] and the quarter phase delay, [4] . However, these approaches satisfy the prescribed delay of 1 & 0.5 over the whole frequency band irrespective of the magnitude response of $H_{\left\{0,u\right\}}(\; z\; )$

For our proposed double density CWT filter structure, the proposed one-dimensional DD CWT is designed as follows:

Construct $H_{0u}\left(z\right)$ & $H_{0d}\left(z\right)$ to satisfy the half band property and group delay constraints as follows:

Assume $H_u\left(z\right)$ to be an N-tap FIR, N = even, of the form

$H_u\left(z\right)={\alpha }_0+{\alpha }_1{z}^{-1}+{\alpha }_2{z}^{-2}+\cdots +{\alpha }_{N-1}{z}^{-\left(N-1\right)}$

where the unknown α’s are determined to satisfy the following two conditions:

1) The function $R\left(z\right)=H_u\left(z\right)H_u\left(z^{-1}\right)$ must be a half band function as in PR systems.

2) The group delay of the rational function $F\left(z\right)=\frac{H_d\left(z\right)}{H_u\left(z\right)}=\frac{z^{-\left(N-1\right)}{H}_u\left(z^{-1}\right)}{H_u\left(z\right)}$ should approximate the delay $\tau =1$ for the first level of the decomposition bank and 0.5 for the succeeding stages in a least squares sense over a fraction $r_x$ of the pass-band of the low pass filter $H_u\left(z\right)$ , respectively [15] .

Then we construct the filters $\left[H_{1u},H_{2u}\right]$ of the upper tree according to the following analysis:

$\begin{array}{l}{H}_{1u}\left(z\right)H_{1u}\left(z^{-1}\right)+H_{2u}\left(z\right)H_{2u}\left(z^{-1}\right)=H_{0u}\left(-z\right)H_{0u}\left(-z^{-1}\right)\\ H_{1u}\left(-z\right)H_{1u}\left(z^{-1}\right)+H_{2u}\left(-z\right)H_{2u}\left(z^{-1}\right)=-H_{0u}\left(-z\right)H_{0u}\left(z^{-1}\right)\end{array}$ (5)

similar relations are applied for $H_{1d}\left(z\right)$ , $H_{2d}\left(z\right)$ . Hence, the design of $H_{1u}\left(z\right),H_{2u}\left(z\right)$ proceeds as follows

1) Construct $H_0\left(z\right)$ , from regular filter design [13] , keeping in mind that it is a half band filter;

2) Take an initial $H_1\left(z\right)$, that is scaled such that its frequency response is less than that of $H_0\left(-z\right)$ at all z point;

3) Construct $H_2\left(z\right)$ from the root finding and perfect reconstruction property, Equation (4);

4) Check the alias free property, and update $H_1\left(z\right)$ ;

5) Reiterate until conversion is obtained.

Design coefficients of $H_0\left(z\right)$ & $H_1\left(z\right)$ & $H_2\left(z\right)$ , $\tau =0.5$ , and 1, for the upper and lower tree, respectively, are listed.

This shows that DD-CWT is designed as a double density 3-channel decomposition for each of $H_{0u},H_{0d}$ , i.e. the low pass filter of the upper and lower (down) trees.

Figure 4 shows the 1-D DD CWT for the upper trees, while Figure 5 shows the completer set of upper and lower trees for the 1-D DD CWT structure.

This design of a 1-D DD-CWT filter along with its coefficients along with the 1-D DD-DWT design methodology presented in previous section is our primary contribution of this work (Table 2).

Figure 5. CWT Filter bank structure-two trees.

(a) (b)

Table 2. (a) Design coefficients of proposed DD-CWT τ = 1; (b) Design coefficients of proposed DD-CWT τ = 0.5.

The 2-D Double-Density Complex Wavelet Transform i.e. 2-D DD-CWT is consequently implemented by applying the filters $H_0\left(z\right)$ & $H_1\left(z\right)$ & $H_2\left(z\right)$ of $H_{0u}$ , first to the rows, then to the columns of an image, for every stage. This would result in nine 2-D sub bands, for every stage, where one of them is the 2-D lowpass scaling filter, and the other eight make up remaining 2-D wavelet filters. For a 3 stage 2-D DD-CWT decomposition we would have a resulting total of 25 2-D subbands as will be detailed in the simulation result section.

3. Double Density Bivariate Denoising Technique

Classical wavelet based denoising techniques, are based on thresholding wavelet coefficients. Different techniques are used to determine these thresholding levels. They are mainly based on thresholding the wavelet decomposition of the noisy image at every sub band by a specific thresholding parameter. The problem can be formulated as: Given the noisy wavelet coefficient $w_n$ , it is required to recover the clean wavelet coefficient w where $w_n=w+n$ , n is the associated independent noise. This is a Maximum Likelihood Estimation (MLE) problem, as solution can be found as,

$w=\max \left(p\left(w|w_n\right)\right)=\max \left(p\left(n\right)p\left(w\right)\right)$ , p is the pdf distribution. In case of zero mean Gaussian noise, $p\left(n\right)$ can be formulated in terms of its variance ${\sigma }_n^2$ . In this case, the variance can be estimated using the empirical formula ${\sigma }_n^2=\text{median}\left(w_n\right)/0.6745$ .

As far as w, it has been observed that the pdf of the wavelet coefficients of natural images approximates Laplacian distribution, [7] .

In [7] , it has been observed that there exist strong dependencies between neighbor wavelet coefficients, such as between parent coefficients at a coarser scale and its adjacent children coefficients at a finer scale. To check this dependence, we construct the CWT system of Lena image using the CWT system of section 2. The contour plots suggest the validity of the empirical circular joint pdf formula of [7] , Figure 6(b) i.e.

$p\left(w_1,w_2\right)=\frac{1}{2\pi {\sigma }^2}{\text{e}}^{-\frac{3}{\sigma }\sqrt{w_1^2+w_2^2}}$

where $w_1,w_2$ are the parent and children wavelet coefficients, respectively.

In this paper, we estimate the variance ${\sigma }^2=\min \left({\sigma }_1^2,{\sigma }_2^2\right)$ , where ${\sigma }_1^2$ and ${\sigma }_2^2$ are variances of $w_1$ and $w_2$ , respectively. Thus, in order to de-noise a noisy image, and in view of the assumption of near circular joint pdf distribution between adjacent scale wavelet coefficients, the thresholded children coefficient $w_{n_2}^{th}$ is given by, [7]

$w_{n_2}^{th}=w_{n_2}\cdot \frac{{\left(\sqrt{w_{n_1}^2+w_{n_2}^2}-\frac{\sqrt{3}{\sigma }_n^2}{\sigma }\right)}_{+}}{\sqrt{w_{n_1}^2+w_{n_2}^2}}$

where the soft thresholding function ${\left(x\right)}_{+}=x$ if $x>0$ and zero otherwise.

(a)(b)

Figure 6. (a) Histogram of Joint pdf. (b) Contour plot of Joint pdf.

Thus, the proposed de-noising scheme amounts to thresholding the wavelet coefficients of the real and imaginary wavelets of the upper and lower trees. In this scheme, the noise variance ${\sigma }_n^2$ is accurately estimated through estimating the pdf of the detail coefficients of the first level wavelet decomposition of the noisy image. The de-noising scheme is summarized as follows:

1) For a prescribed number of decomposition levels n and a prescribed number of vanishing moments K, determine the first and succeeding stages filters of the upper and lower trees of the dual tree DWT, as described in section 2.

2) Initially, for the first scale of the upper and lower trees, evaluate the real and imaginary parts of the complex wavelets $w_{n_2}$ and its adjacent parent $w_{n_1}$ at the coarser scale. Interpolate by 2 the parent coefficient $w_{n_1}$ .

3) In order to estimate $\sigma$ , ${\sigma }_1^2$ and ${\sigma }_2^2$ have to be estimated. They are estimated as the peak powers $E_1,E_2$ of the coefficients $w_{n_1},w_{n_2}$ . Then, ${\sigma }_1^2=E_1-{\sigma }_n^2$, ${\sigma }_2^2=E_2-{\sigma }_n^2$ , and ${\sigma }^2=\min \left({\sigma }_1^2,{\sigma }_2^2\right)$ . Threshold the children’s coefficients.

Repeat steps 2,3 until all children coefficients are scanned. In the last scale, threshold the parent coefficient $w_{n_1}$ as well.

In this paper, we apply this bivariate shrinkage technique to denoise 2D DD-DWT and 2D DD-CWT wavelet packets. Bivariate shrinkage technique is based on strong dependency between noisy children wavelet coefficients $x_1^n$ and their corresponding noisy parent coefficients $x_2^n$ at coarser scale. Figure 6 illustrates this joint pdf relation for an image contaminated with AWGN with noise variance ${\sigma }_n^2$ . This means that to denoise the image; we have to maximize the joint conditional probabilities

$p_x=p\left(x_1,x_2|x_1^n,x_2^n\right)$ (6)

where $x_1,x_2$ denote the clean wavelet ceofficients.

Using Bays rule; it turns out that

$p_x=p\left(n_1,n_2\right)\ast p\left(x_1,x_2\right)$ (7)

where $n_1=x_1^n-x_1,n_2=x_2^n-x_2$ are two independent Gaussian noise with variance ${\sigma }_n^2$ .

In order to obtain a closed form solution for $x_1,x_2$ we express the joint $p\left(x_1,x_2\right)$ using the empirical formulation as in [7]

$p\left(x_1,x_2\right)=\frac{3}{2\pi {\sigma }^2}\ast {\text{e}}^{-\frac{\sqrt{3}}{\sigma }\sqrt{x_1^2+x_2^2}}$ (8)

where ${\sigma }^2=\sqrt{{\sigma }_1^2\ast {\sigma }_2^2}$ ; ${\sigma }_1^2,{\sigma }_2^2$ denote the variances of $x_1,x_2$ respectively. This is verified by the near circular performance of the noiseless joint pdf distribution shown in Figure (3). Maximization of Equation (7) yields

$x_k=\frac{{\left(\sqrt{{\left(x_k^n\right)}^2+{\left(x_{k+1}^n\right)}^2}-\frac{\sqrt{3}{\sigma }_n^2}{\sigma }\right)}_{+}}{\sqrt{x_k^2+x_{k+1}^2}}\cdot x_k^n$ (9)

We note here that all DD-CWT or DD-DWT decompositions presented in the work adopts this bivariate denoising scheme between its parent and children bands.

4. Enhanced Scenarios for the Proposed Hybrid Bivariate DD-DWT and DD-CWT Denoising

In this section we propose two scenarios to enhance the performance of denoising, then we summarize how to achieve the highest performance by fusing different bands from different decompositions structures, i.e. CWT, DWT, and combine them with optimized factors to synthesize an enhanced denoising image.

4.1. Edge Sharpening

We here propose to enhance the noisy image edges by processing it through a Laplacian 2-D filter with an optimized factor and add the result to the original noisy image. This Laplacian 2-D filter would have a sharpening factor that is also optimized. The output image from this Edge sharpening stage would be as the following equation:

$w_{\text{updated}}=\alpha Lap\left(w\right)+w$

where w is original noisy CWT or DWT wavelet band, $w_{\text{updated}}$ represents the final band coefficients, and $\alpha$ is the optimization factor.

4.2. Eigen Analysis

We here propose to decompose the noisy image DD-CWT or DD-DWT with a singular value decomposition. In this decomposition we produce a diagonal matrix D, with nonnegative diagonal elements in decreasing order, and unitary matrices U and V, all of the same size of the input image band, according to following equation

$\left[U\cdot D\cdot V^t\right]=SVD(\; w\; )$

where w is original noisy CWT or DWT wavelet band, D is the diagonal matrix that will be rescaled in decreasing order, and $w_{\text{updated}}$ is reconstructed according to the following synthesis equation

$U^{t}DV=w_{\text{updated}}$

We then select the highest D matrix diagonal element and scale it to an updated denoised value and then reconstruct the wavelet band. This updated denoised value is obtained through some empirical analysis and also by denoised learning scenarios.

4.3. Hybrid DD-DWT and DD-CWT Denoising

Finally in this section we proposed to decompose the denoised enhanced images with either DD-CWT or DD-DWT by a further DD-CWT structure for each image, where each band in each denoised image is multiplied by an optimized factor and then synthesize this CWT decomposition. This optimization process that is performed for all the CWT bank factors and is aimed at reducing edge energy in the final output denoised image. CWT optimization is performed according to the following equation.

$w_{\text{hybridoptimized}}={\alpha }_1{w}_1+{\alpha }_2{w}_2+\cdots +{\alpha }_i{w}_i+\cdots +{\alpha }_{\text{last}}{w}_{\text{last}}$

where $w_{\text{hybridoptimized}}$ is final hybrid optimized wavelet band and each i original band is multiplied with an ${\alpha }_i$ multiplication factor.

We note here that this further denoising enhancement by a factorized optimized DD-CWT structure would represent a hybrid denoising technique that combines merits of both DD-DWT and DD-CWT in a fusion manner, especially if it was originally denoised with DD-DWT. The ${\alpha }_i$ factors would make this factor adaptive as we can select some high density bands to get higher weight in the final denoised image. This fusion hybrid technique would achieve significant improvement in the final PSNR or SSIM values as will be shown in next section. We justify this improvement to the optimization process that selects the highest energy bands.

5. Simulation Examples

In this section we first show the performance of 2D DD-DWT denoising after explaining its experimental procedure. We then show the performance of 2D DD-CWT that in some circumstances achieves better PSNR results. Then we show how each of our proposed two enhancement scenarios, Edge sharpening and Eigen Analysis, would enhance the final denoised image. We finally illustrate how our proposed hybrid DWT and CWT fusion denoising methodology would achieve the most enhancement performance.

We note here that the enhancement in denoising is measured in terms of reduced energy around edges. This was also compared with PSNR values of the output denoised image compared with original clean image, before any noise attack. This PSNR calculation with the clean image could be unpractical in many denoising applications where the clean image is not available for comparisons, but it is only mentioned just to verify that enhancement in PSNR values is consistent with edge energies minimizations that is performed in all proposed techniques.

In these simulations; ${\sigma }_n^2$ is estimated using the pdf technique of the first wavelet detail coefficient as described in [17] . Simulations of several noisy images verify that this choice yields the highest denoised image

The proposed filter design $H_0\left(z\right)$ has $K=4$ zeros at $z=-1$ ; $H_1\left(z\right)$ is constructed with $M=2$ zeros at $z=1$ using the maximum phase technique.

For space limitations, we have two 256 × 256 Cameraman & Lena images contaminated with zero mean Gaussian noise AWGN with different variance ${\sigma }_n^2$ . The number of decomposition levels is 3 for either DWT or CWT. Figures 7-10 shows samples of our Denoising results, with a noisy image example ${\sigma }_n^2=0.1$ , 0.05, or 0.2.

Figure 7. Example of Noisy image and before Denoising, ${\sigma }_n^2=0.05$ .

Figure 8. Denoising performance from [7] and then proposed.

Figure 9. Example of Noisy image and before Denoising, ${\sigma }_n^2=0.01$ .

Figure 10. Denoising performance from [7] and then proposed.

The proposed 2D-DD DWT Bivariate denoising technique is implemented on 2D noisy images as follows:

1) Given the 2D-DD-DWT filters $H_0\left(z\right)$ , $H_1\left(z\right)$ and $H_2\left(z\right)$ as in section 2; decompose the noisy image $I_n$ through wavelet packet structure using these filters.

2) For each of the 8 decomposed wavelet coefficients in all sub bands levels $x_k^n$ , $k=1,2,\cdots ,J$ levels; apply Bivariate shrinkage technique between $x_k^n$ and its corresponding at coarse scale $x_{k+1}^n$ . Expand $x_{k+1}^n$ by 2 to have the same size as $x_k^n$ . Estimate ${\sigma }_1^2,{\sigma }_2^2$ for both $x_k,x_{k+1}$ .

3) Reconstruct the denoised image $\hat{I}$ as in section 2.

The proposed 2D-DD CWT Bivariate denoising technique is implemented on 2D noisy images as follows:

4) Decompose the input Noisy image by its rows through the proposed 2-D DD CWT decomposition structure in Figure 5 for 3 stages, for each stage we would have 1 lowpass scaling subband and 8 bandpass higher frequency bands. The scaling subband is the only subband that gets further decomposed to the next stage in a wavelet lattice structure. This would result in total of 25 2-D subbands.

5) Repeat step one but through columns of the Noisy image.

6) The Bivariate shrinkage technique proposed in section 3 is applied between each scaling subband and its children for all stages.

7) Reconstruct the denoised image by rows and columns in a reverse manner of the decomposition.

Table 3 & Table 4 tabulate the denoising PSNR results when the noisy images are processed through our two proposed enhancing stages of Edge Sharpening and Eigen Analysis.

Results with Edge sharpening

Table 3. PSNR performance in dB for cameraman.

Results with Eigen Analysis

Table 4. PSNR performance in dB for cameraman.

Table 5 & Table 6 tabulate the denoising PSNR results when the noisy images are contaminated with different ${\sigma }_n^2$ . These results show the superiority of the proposed denoising technique DD-DWT, underlined over other DD-DWT [8] , DD-CWT proposed, CWT [7] technique.

In our final hybrid DD-DWT and DD-CWT fusion technique we decomposed the image first through either DD-DWT or DD-CWT decomposition for denoising, then we further enhance the denoising by a factorized optimized DD-CWT process. Subband multiplication factors are our main optimization variables in our simulation results as in Table 7, for Cameraman.

Table 5. PSNR performance in dB for cameraman.

Table 6. PSNR performance in dB for Lena image.

Table 7. Results from Hybrid DD-DWT/CWT denoising.

6. Discussion and Conclusion

In this paper, we proposed in more details the usage of DWT or CWT decompositions in image denoising with the adoption of Double density analysis. We proposed two methodologies to enhance these decompositions for either DWT or CWT scenarios. We also presented an adaptive hybrid technique for Bivariate based image denoising that is based on the synthesis of DD-DWT bands or DD-CWT bands but with different weights, to deliver enhanced image features with less denoising impacts. Simulation results have shown that the DD_DWT bivariate shrinkage achieves the best performance of all the denoising schemes considered. From Equations (3), (4) it is clear that there is a plenty of solutions for H₂(z) that satisfy the alias free conditions, yet simulations have shown that the maximum phase solution yields the optimum denoising performance.

Acknowledgements

This work has been funded mainly from Alexander Von Humboldt foundation, Germany, 2019. Early results of this work was presented in [10] at much less scope and scale.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Fernandes, F., van Spaendonck, R., Coates, M.J. and Burrus, S. (2000) Directional Complex Wavelet Processing. In: Aldroubi, A., Laine, A.F. and Unser, M.A., Eds., Wavelet Applications in Signal and Image Processing VIII, Vol. 4119, SPIE, Bellingham, 536-546. https://doi.org/10.1117/12.408642
[2]	Fernandes, F., van Spaendonck, R.L.C. and Burrus, C.S. (2003) A New Framework for Complex Wavelet Transforms. IEEE Transactions on Signal Processing, 51, 1825-1837. https://doi.org/10.1109/TSP.2003.812841
[3]	Donoho, D.L. (1995) De-Noising by Soft-Thresholding. IEEE Transactions on Information Theory, 41, 613-627. https://doi.org/10.1109/18.382009
[4]	Selesnick, I.W. and Bayram, I. (2014) Sparse Signal Estimation by Maximally Sparse Convex Optimization. IEEE Transactions on Signal Processing, 62, 1078-1092. https://doi.org/10.1109/TSP.2014.2298839
[5]	Ding, Y. and Selesnick, I.W. (2015) Artifact-Free Wavelet Denoising: Non-Convex Sparse Regularization, Convex Optimization. IEEE Signal Processing Letters, 22, 1364-1368. https://doi.org/10.1109/LSP.2015.2406314
[6]	Sendur, L. and Selesnick, I.W. (2002) Bivariate Shrinkage Functions for Wavelet-Based Denoising Exploiting Inter-Scale Dependency. IEEE Transactions on Signal Processing, 50, 2744-2756. https://doi.org/10.1109/TSP.2002.804091
[7]	Fahmy, M.F. and Fahmy, O.M. (2018) Efficient Bivariate Image Denoising Technique Using New Orthogonal CWT Filter Design. IET Image Processing, 12, 1354-1360. https://doi.org/10.1049/iet-ipr.2017.1117
[8]	Selesnick, I.W. (2004) The Double-Density Dual-Tree DWT. IEEE Transactions on Signal Processing, 52, 1304-1314. https://doi.org/10.1109/TSP.2004.826174
[9]	Selesnick, I.W. (2002) The Design of Approximate Hilbert Transform Pairs of Wavelet Bases. IEEE Transactions on Signal Processing, 50, 1144-1152. https://doi.org/10.1109/78.995070
[10]	Fahmy, G., Fahmy, M.F. and Fahmy, O. (2021) Bivariate Double Density Discrete Wavelet for Enhanced Image Denoising. 2021 3rd Novel Intelligent and Leading Emerging Sciences Conference (NILES), Giza, 23-25 October 2021, 308-311. https://doi.org/10.1109/NILES53778.2021.9600554
[11]	Fahmy, G., Fahmy, M.F. and Fahmy, O.M. (2017) Micro-Movement Magnification in Video Signals Using Complex Wavelet Analysis. IET Image Processing, 11, 986- 993. https://doi.org/10.1049/iet-ipr.2017.0049
[12]	Selesnick, I.W., Baraniuk, R.G. and Kingsbury, N. C. (2005) The Dual-Tree Complex Wavelet Transform. IEEE Signal Processing Magazine, 22, 123-151. https://doi.org/10.1109/MSP.2005.1550194
[13]	Selesnick, I.W. (2001) Hilbert Transform Pairs of Wavelet Bases. IEEE Signal Processing Letters, 8, 170-173. https://doi.org/10.1109/97.923042
[14]	Fahmy, M.F. and Fahmy, O.M. (2016) An Enhanced Denoising Technique Using Dual Tree Complex Wavelet Transform. 2016 33rd National Radio Science Conference (NRSC), Aswan, 22-25 February 2016, 205-211. https://doi.org/10.1109/NRSC.2016.7450829
[15]	Fahmy, O.M., Fahmy, G. and Fahmy, M.F. (2018) A New Video Magnification Technique Using Complex Wavelets with Radon Transform Application. Signal, Image and Video Processing, 12, 1505-1512. https://doi.org/10.1007/s11760-018-1306-9
[16]	Bhonsle, D. and Dewangan, S. (2012) Comparative Study of Dual-Tree Complex Wavelet Transform and Double Density Complex Wavelet Transform for Image Denoising Using Wavelet-Domain. International Journal of Scientific and Research Publications, 2.
[17]	Fahmy, M.F. and Fahmy, G. (2014) Exponential Spline Perfect Reconstruction Decomposition with Applications in Compression and Denoising. Signal, Image and Video Processing, 8, 1111-1120. https://doi.org/10.1007/s11760-014-0640-9