Quaternion Wavelet Analysis and Application in Image Denoising

Abstract

The quaternion wavelet transform is a new multiscale analysis tool. Firstly, this paper studies the standard orthogonal basis of scale space and wavelet space of quaternion wavelet transform in spatial 𝐿²(𝑅²), proves and presents quaternion wavelet’s scale basis function and wavelet basis function concepts in spatial scale space 𝐿²(𝑅²;H), and studies quaternion wavelet transform structure. Finally, the quaternion wavelet transform is applied to image denoising, and generalized Gauss distribution is used to model QWT coefficients’ magnitude distribution, under the Bayesian theory framework, to recover the original coefficients from the noisy wavelet coefficients, and so as to achieve the aim of denoising. Experimental results show that our method is not only better than many of the current denoising methods in the peak signal to noise ratio (PSNR), but also obtained better visual effect.

1. Introduction

Wavelet analysis is a rapidly developing branch of mathematics since 1980s; its research has just been unfolding. As a mathematical tool, wavelet transform is a major breakthrough of the Fourier transform and Fourier transform window known to many people; it has good time-frequency features and multiple resolution and wavelet analysis theory has become one of the most useful tools in signal analysis, image processing, pattern recognition and other fields. In image processing, the basic idea of the wavelet transform is to decompose image multiresolution; the original image is decomposed into different space and different frequency sub-image, and then coefficients of sub-image are processed. Commonly used wavelet transforms are real discrete wavelet transform and complex wavelet transform and so on. The real discrete wavelet and complex wavelet transform have two generals shortcoming; first, the real discrete wavelet transform signal small shift will produce the energy of wavelet coefficient distribution change; second, dual-tree complex wavelet although overcame the first problem, it can generate signal phase ambiguity when represented two-dimensional image’s features. While the quaternion wavelet transform is a new multiscale analysis image processing tool, it is based on the Hilbert two-dimensional transform theory, which has approximate shift invariance and can well overcome the above drawbacks [1–3].

At present, quaternion wavelet research is divided into two branches, one is based on quaternion numerical function multiresolution analysis theory of quaternion wavelet, using a single tree structure, the earliest in 1994, Mitrea gave quaternion wavelet form concept [4]; in 2001, Traversoni used real wavelet transform and complex wavelet transform by quaternion Haar kernel and proposed discrete quaternion wavelet transform [5] and gave some applications in image processing; in 2004, He and Yu used matrix value function multiresolution analysis structure for consecutive quaternion wavelet transform [6] and gave some properties; in 2010, Bahri constructed discrete quaternion wavelet transform through complex matrix function and proved some basic properties [7]; in 2011, Bahri et al. introduced through quaternion wavelet admissibility conditions, systematically extended the consecutive wavelet transform concept to consecutive quaternion wavelet concept [8], and proved the reconstruction theorem and continuous quaternion wavelet basic properties. But these are mainly the concepts and properties of promotion, because its filters’ structure and implementation are difficulties, it has not made any progress in application at present. Another branch is based on Bulow quaternion analytic signal [9]; using real filter and dual-tree structure to achieve the quaternion wavelet transform, the filter has the advantages of simple structure, relatively easy, and there was quaternion signal application background; this forms the research focus as this paper focuses on this branch. In 2006, Corrochano through quaternion Gabor filter constructed quaternion wavelet transform (QWT) [2] and discussed the QWT properties and wavelet pyramid algorithm, pointing out DWT is without phase, CWT only has one phase, while QWT can provide three phases and puting forward the image multiresolution disparity estimation method based on the theory of QWT; in 2008, based on the dual-tree complex wavelet, Chan et al. used the concepts and properties of Bulow quaternion analytic signal and quaternion Fourier transform, constructed (dual-tree) quaternion wavelet transform, and worked out the meaning of three phases, two of which represent the image of local displacement information, another as image texture feature, which can be used to estimate the image of the local geometric features [3]; in 2010, Xu et al. used quaternion wavelet transform’s amplitude and phase method and applied it to the face recognition [10], also obtained certain result; in 2011, Soulard and Carré applied the quaternion wavelet transform to image texture analysis [11] and proved the feasibility of this method.

In most cases, the image corruption is commonly modeled by a zero-mean additive white Gaussian random noise leading to the following additive degradation model: 𝑔(𝑥,𝑦)=𝑓(𝑥,𝑦)+𝜀(𝑥,𝑦), where 𝑓(𝑥,𝑦), 𝑔(𝑥,𝑦), and 𝜀(𝑥,𝑦) represent the clean image, the observed noisy image, and the white Gaussian noise with variance 𝜎²_𝑛. The traditional image denoising methods are mainly the Wiener filter, median filter, and weighted median filter, and so forth, although these methods restrain the noise, they lose a lot of image detail information and often cause image blur and generate a ringing phenomenon. Wavelet analysis has good properties of time-frequency localization and multiresolution denoising ability, which can well reserve the detail of the image, thus becoming the main tool for image denoising. The earliest wavelet denoising method is proposed by Mallat and Hwang in 1992 based on modulus maximum denoising method [12]; its principle is observed through different scales of the wavelet coefficients’ modulus maximal of regularity, removed by the noise generated by the module maximum, with the remaining modulus maximal recovery signal. In 1994, Donoho put forward to adopt universal threshold wavelet shrinkage method [13]; the basic idea is comparison between different scale coefficients’ module a certain threshold and obtained the de-noised signal by the inverse transform. Although the implementation of thresholding method is simple, it does not take into account the correlation between wavelet coefficients. Due to that the wavelet coefficients have strong correlation, researchers put forward many prior statistical models of the wavelet coefficients, so as to get better denoising results. Then most of the research is based on a statistical model of the noise reduction method, the key technology of this kind method is to choose an accurate prior probability distribution model for the wavelet coefficients. In 1998, Crouse et al. suggested hidden Markov model (HMM) [14] in wavelet transform domain; this model can well capture interscale wavelet coefficients correlation, but it needs to get the parameters through iterative calculation which cost much time. In 1999, Mihçak et al. presented LAWML and LAWMAP methods [15], as the methods adopted local adaptive coefficients distribution model, so it obtained better denoising results. However, retained too much small wavelet coefficients, severe burr phenomenon appeared in the reconstructed image. In 2000, Chang et al. used generalized Gaussian distribution simulates intrascale wavelet coefficients’ relationship [16], the model can well reflect the wavelet coefficient correlation between scales, and obtained certain results. In 2002, Şendur and Selesnick considered the correlation of interscale wavelet coefficients and put forward BiShrink [17] and Local—BiShrink [18] denoising methods and its main utilize bivariate non-Gaussian distribution to manifest the correlation of wavelet coefficients between different scale and as a prior distribution. Experiments show that denoising results based on bivariate non-Gaussian distribution model can achieve very good results. In 2003, Portilla et al. [19] proposed the GSM model which can indicate marginal distribution of coefficients and embody the correlation of neighborhood coefficients. Later, many effective image denoising methods based on this model combining different transforms are obtained [20, 21]. In 2005, Cho and Bui proposed the multivariate generalized Gaussian distribution model [22], which adjusts different parameters and can include Gaussian, generalized Gaussian, and non-Gaussian model, but parameters’ estimation is more complex in image denoising process.

This paper gives and proves the properties of the Hilbert transform and the standard orthogonal basis of scale space and wavelet space of quaternion wavelet transform in spatial 𝐿²(𝑅²), and then in the space 𝐿²(𝑅²;H) it gives quaternion wavelet basis function and scale function concepts, further gives quaternion wavelet transform concept, and also studies quaternion wavelet transform’s structure and filter structure. Its application is applied in image denoising; generalized Gauss distribution is used to model the QWT coefficients’ magnitude, and then using the Bayesian minimum mean square error estimation method, to recovery the original coefficients’ magnitude from the noisy wavelet coefficients’ magnitude, and so as to achieve the purpose of denoising. Experimental results show that our method is better than many other current denoising methods.

2. Quaternion Analytical Signal

2.1. 2D Hilbert Transform Concept

Definition 2.1 (see [8]). Set 𝑓(𝑥,𝑦)∈𝐿²(𝑅²), then 𝑓(𝑥,𝑦) along the 𝑥-axis and the 𝑦-axis and along the 𝑥,𝑦 axis of the Hilbert transform 𝑓_H𝑥(𝑥,𝑦), 𝑓_H𝑦(𝑥,𝑦), and 𝑓_H𝑥𝑦(𝑥,𝑦), respectively: 𝑓_H𝑥=1𝑥,𝑦𝜋_𝑅𝑓𝜉₁,𝑦𝑥−𝜉₁𝑑𝜉₁, (2.1) the corresponding frequency domain is 𝑓_H𝑥(𝑢,𝑣)=−𝑗sgn(𝑢)𝑓(𝑢,𝑣); 𝑓_H𝑦=1𝑥,𝑦𝜋_𝑅𝑓𝑥,𝜉₂𝑦−𝜉₂𝑑𝜉₂, (2.2) the corresponding frequency domain is 𝑓_H𝑦(𝑢,𝑣)=−𝑗sgn(𝑣)𝑓(𝑢,𝑣); 𝑓_H𝑥𝑦=1𝑥,𝑦𝜋²_𝑅²𝑓𝜉₁,𝜉₂𝑥−𝜉₁𝑦−𝜉₂𝑑𝜉₁𝑑𝜉₂, (2.3) the corresponding the frequency domain is 𝑓_H𝑥𝑦(𝑢,𝑣)=−sgn(𝑢)sgn(𝑣)𝑓(𝑢,𝑣).

2.2. Quaternion Analytical Signal

The quaternion is proposed by W. R. Hamilton in 1843; quaternion is a complex promotion, and it can be regarded as a special Clifford algebra.

Set H={𝑞=𝑞_𝑟+𝑖𝑞_𝑖+𝑗𝑞_𝑗+𝑘𝑞_𝑘∣𝑞_𝑟,𝑞_𝑖,𝑞_𝑗,𝑞_𝑘∈𝑅}. Three imaginary 𝑖,𝑗,𝑘 satisfy the following rules: 𝑖²=𝑗²=𝑘²=−1,𝑖𝑗=−𝑗𝑖=𝑘,𝑗𝑘=−𝑘𝑗=𝑖,𝑘𝑖=−𝑖𝑘=𝑗. (2.4) And define 𝑞 as the quaternion.

The conjugate of quaternion 𝑞 is given by 𝑞=Re(𝑞)−Im(𝑞)=𝑞_𝑟−𝑖𝑞_𝑖−𝑗𝑞_𝑗−𝑘𝑞_𝑘 and the module of quaternion 𝑞 is given by |𝑞|=𝑞𝑞=𝑞²_𝑟+𝑞²_𝑖+𝑞²_𝑗+𝑞²_𝑘. In addition, quaternion 𝑞 can also express as: 𝑞=|𝑞|𝑒^𝑖𝜑𝑒^𝑗𝜃𝑒^𝑘𝜓, where |𝑞| is the modulus of 𝑞, (𝜑,𝜃,𝜓) are the-three phase angles which is uniquely defined within the range (𝜑,𝜃,𝜓)∈[−𝜋,𝜋]×[−𝜋/2,𝜋/2]×[−𝜋/4,𝜋/4][9].

So 𝐿²(𝑅²;H)={𝑓∶𝑅²→H,∫_𝑅2|𝑓(𝑥)|²𝑑²𝑥<∞}, forall𝑓,𝑔∈𝐿²(𝑅²;H), define the inner product ⟨𝑓,𝑔⟩=∫_𝑅2𝑓(𝑥)𝑔(𝑥)𝑑²𝑥, it is not difficult to prove 𝐿²(𝑅²;H) is a Hilbert space.

Quaternion wavelet transform is based on quaternion analytic signal, and we give the concept.

Definition 2.2 (see [8]). Set 𝑓(𝑥,𝑦) is a real two-dimensional signals, then quaternion analytic signal can be defined as 𝑓^𝑞(𝑥,𝑦)=𝑓(𝑥,𝑦)+𝑖𝑓_H𝑥(𝑥,𝑦)+𝑗𝑓_H𝑦(𝑥,𝑦)+𝑘𝑓_H𝑥𝑦(𝑥,𝑦).
where 𝑓_H𝑥(𝑥,𝑦),𝑓_H𝑦(𝑥,𝑦),𝑓_H𝑥𝑦(𝑥,𝑦) are the 𝑓(𝑥,𝑦) along the 𝑥-axis and the 𝑦-axis and along the 𝑥,𝑦 axis of the Hilbert transform.

3. Quaternion Wavelet Transform

This section is based on the study in [2, 3], in-depth study of quaternion wavelet theory, which proved and presented the correlative properties and concepts of scale basis and wavelet basis of quaternion wavelet, further, provided discrete quaternion wavelet transform concept, and studied quaternion wavelet transform structure and filter structure.

3.1. Quaternion Wavelet’s Scale Basis and Wavelet Basis

Lemma 3.1. Let {𝜑(𝑥−𝑘)}_𝑘∈𝑍 be the space of 𝑉⊂𝐿²(𝑅) standard orthogonal basis, set 𝜑_𝑔(𝑥)=H(𝜑(𝑥)) (i.e., that 𝜑(𝑥) of Hilbert transform), so {𝜑_𝑔(𝑥−𝑘)}_𝑘∈𝑍 is the standard orthogonal basis of space 𝑉=H𝑉⊂𝐿²(𝑅).

Proof. According to the theorem of the paper in [23], we know that standard orthogonal basis in the frequency domain is equivalent to: ∑_𝑘∈𝑍|𝜑(𝜔+2𝑘𝜋)|²=1, by the concept and property of Hilbert transform, the establishment of conclusion is right.

According to the lemma, it is not difficult to prove the following.

Theorem 3.2. Let {𝑉_𝑗}_𝑗∈𝑧 be a 1D orthogonal multiresolution analysis, the corresponding scale function and wavelet function are 𝜑_ℎ(𝑥) and 𝜓_ℎ(𝑥), respectively; note that H_𝑥(𝜑_ℎ(𝑥)𝜑_ℎ(𝑦))=𝜑_𝑔(𝑥)𝜑_ℎ(𝑦) (H_𝑥(𝜑_ℎ(𝑥)𝜑_ℎ(𝑦)) are said to function as (𝜑_ℎ(𝑥)𝜑_ℎ(𝑦)) along 𝑥-axis direction Hilbert transform, and the corresponding mark is similar H_𝑦𝜑_ℎ(𝑥)𝜑_ℎ𝑦=𝜑_ℎ(𝑥)𝜑_𝑔𝑦,H_𝑥𝑦𝜑_ℎ(𝑥)𝜑_ℎ𝑦=𝜑_𝑔(𝑥)𝜑_𝑔𝑦. (3.1) So {𝜑_𝑔(𝑥−𝑙)𝜑_ℎ(𝑦−𝑘)}_{𝑙,𝑘∈𝑍} (or {𝜑_ℎ(𝑥−𝑙)𝜑_𝑔(𝑦−𝑘)}_{𝑙,𝑘∈𝑍}) are the standard orthogonal base of the space V₀⨂𝑉₀ and the space 𝑉₀⨂𝑉₀, where 𝑉₀=H𝑉₀.
We know that, in real orthogonal wavelet form 2D tensor product of wavelet, having a scale function and form of standard orthogonal basis in scale space, there are three wavelet functions form of standard orthogonal basis in wavelet space. For a function, it can be made along the 𝑥,𝑦 and 𝑥,𝑦 axis direction’s Hilbert transform, so by Theorem 3.2 one knows {𝜑_𝑔(𝑥−𝑙)𝜑_ℎ(𝑦−𝑘)}_{𝑙,𝑘∈𝑍} (or {𝜑_ℎ(𝑥−𝑙)𝜑_𝑔(𝑦−𝑘)}_{𝑙,𝑘∈𝑍}; total of 4) is the standard orthogonal basis of scale space 𝐿²(𝑅²), {𝜑_𝑔(𝑥−𝑙)𝜓_ℎ(𝑦−𝑘)}_{𝑙,𝑘∈𝑍}({𝜑_ℎ(𝑥−𝑙)𝜓_𝑔(𝑦−𝑘)}_{𝑙,𝑘∈𝑍}; total of 12) is the standard orthogonal base of wavelet space 𝐿²(𝑅²).
In fact, by Theorem 3.2 and the above analysis, it is easy to prove the following.

Theorem 3.3. Let {𝑉_𝑗}_𝑗∈𝑧  be a 1D 𝐿²(𝑅) orthogonal multiresolution analysis, 𝜑_ℎ(𝑥),𝜓_ℎ(𝑥) are scale function and wavelet function of corresponding respectively.
Note Φ^𝑞𝑥,𝑦=𝜑_ℎ(𝑥)𝜑_ℎ𝑦+𝑖𝜑_𝑔(𝑥)𝜑_ℎ𝑦+𝑗𝜑_ℎ(𝑥)𝜑_𝑔𝑦+𝑘𝜑_𝑔(𝑥)𝜑_𝑔𝑦,Φ^𝑞_{𝑗,𝑘,𝑚}𝑥,𝑦=𝜑_{ℎ,𝑗,𝑘}(𝑥)𝜑_{ℎ,𝑗,𝑚}𝑦+𝑖𝜑_{𝑔,𝑗,𝑘}(𝑥)𝜑_{ℎ,𝑗,𝑚}𝑦+𝑗𝜑_{ℎ,𝑗,𝑘}(𝑥)𝜑_{𝑔,𝑗,𝑚}𝑦+𝑘𝜑_{𝑔,𝑗,𝑘}(𝑥)𝜑_{𝑔,𝑗,𝑚}𝑦, (3.2) where 𝜑_{ℎ,𝑗,𝑘}(𝑥)=2^−𝑗/2𝜑_ℎ(2^−𝑗𝑥−𝑘),𝑗,𝑘,𝑚∈𝑍, corresponding mark is similar.
So {Φ^𝑞_{𝑗,𝑘,𝑚}(𝑥,𝑦)}_{𝑘,𝑚∈𝑍} is the basis function of quaternion wavelet scale space. And called Φ^𝑞(𝑥,𝑦) is quaternion wavelet scale function in 𝐿²(𝑅²;H); called {Φ^𝑞_{𝑗,𝑘,𝑚}(𝑥,𝑦)}_{𝑘,𝑚∈𝑍} is discrete quaternion wavelet scale function in 𝐿²(𝑅²;H).

Further one has the following.

Theorem 3.4. Let {𝑉_𝑗}_𝑗∈𝑧 be a 1D orthogonal 𝐿²(𝑅) multiresolution analysis, and 𝜑_ℎ(𝑥),𝜓_ℎ(𝑥) are the scale function and wavelet function of corresponding, respectively; note that Ψ^𝑞,1𝑥,𝑦=𝜑_ℎ(𝑥)𝜓_ℎ𝑦+𝑖𝜑_𝑔(𝑥)𝜓_ℎ𝑦+𝑗𝜑_ℎ(𝑥)𝜓_𝑔𝑦+𝑘𝜑_𝑔(𝑥)𝜓_𝑔𝑦,Ψ^𝑞,2𝑥,𝑦=𝜓_ℎ(𝑥)𝜑_ℎ𝑦+𝑖𝜓_𝑔(𝑥)𝜑_ℎ𝑦+𝑗𝜓_ℎ(𝑥)𝜑_𝑔𝑦+𝑘𝜓_𝑔(𝑥)𝜑_𝑔𝑦,Ψ^𝑞,3𝑥,𝑦=𝜓_ℎ(𝑥)𝜓_ℎ𝑦+𝑖𝜓_𝑔(𝑥)𝜓_ℎ𝑦+𝑗𝜓_ℎ(𝑥)𝜓_𝑔𝑦+𝑘𝜓_𝑔(𝑥)𝜓_𝑔𝑦. (3.3) The shift and expand of  Ψ^𝑞,1(𝑥,𝑦),  Ψ^𝑞,2(𝑥,𝑦) and Ψ^𝑞,3(𝑥,𝑦) form of function family Ψ^𝑞,2_{𝑗,𝑘,𝑚}𝑥,𝑦=𝜓_{ℎ,𝑗,𝑘}(𝑥)𝜑_{ℎ,𝑗,𝑚}𝑦+𝑖𝜓_{𝑔,𝑗,𝑘}(𝑥)𝜑_{ℎ,𝑗,𝑚}𝑦+𝑗𝜓_{ℎ,𝑗,𝑘}(𝑥)𝜑_{𝑔,𝑗,𝑚}𝑦+𝑘𝜓_{𝑔,𝑗,𝑘}(𝑥)𝜑_{𝑔,𝑗,𝑚}𝑦_{𝑘,𝑚∈𝑧},Ψ^𝑞,1_{𝑗,𝑘,𝑚}𝑥,𝑦=𝜑_{ℎ,𝑗,𝑘}(𝑥)𝜓_{ℎ,𝑗,𝑚}𝑦+𝑖𝜑_{𝑔,𝑗,𝑘}(𝑥)𝜓_{ℎ,𝑗,𝑚}𝑦+𝑗𝜑_{ℎ,𝑗,𝑘}(𝑥)𝜓_{𝑔,𝑗,𝑚}𝑦+𝑘𝜑_{𝑔,𝑗,𝑘}(𝑥)𝜓_{𝑔,𝑗,𝑚}𝑦_{𝑘,𝑚∈𝑧},Ψ^𝑞,3_{𝑗,𝑘,𝑚}𝑥,𝑦=𝜓_{ℎ,𝑗,𝑘}(𝑥)𝜓_{ℎ,𝑗,𝑚}𝑦+𝑖𝜓_{𝑔,𝑗,𝑘}(𝑥)𝜓_{ℎ,𝑗,𝑚}𝑦+𝑗𝜓_{ℎ,𝑗,𝑘}(𝑥)𝜓_{𝑔,𝑗,𝑚}𝑦+𝑘𝜓_{𝑔,𝑗,𝑘}(𝑥)𝜓_{𝑔,𝑗,𝑚}𝑦_{𝑘,𝑚∈𝑧}, (3.4) where 𝜓_{ℎ,𝑗,𝑘}(𝑥)=2^−𝑗/2𝜓_ℎ(2^−𝑗𝑥−𝑘), 𝑗,𝑘,𝑚∈𝑍, the corresponding mark is similar.
And Ψ^𝑞,1(𝑥,𝑦),Ψ^𝑞,2(𝑥,𝑦),Ψ^𝑞,3(𝑥,𝑦) as the quaternion wavelet basis functions in 𝐿²(𝑅²;H) space; that {Ψ^𝑞,1_{𝑗,𝑘,𝑚},Ψ^𝑞,2_{𝑗,𝑘,𝑚},Ψ^𝑞,3_{𝑗,𝑘,𝑚}}_{𝑗,𝑘,𝑚∈𝑧} are discrete quaternion wavelet functions in 𝐿²(𝑅²;H) space.
One can get concept of the quaternion wavelet transform.

Definition 3.5. For for all 𝑓(𝑥,𝑦)∈𝐿²(𝑅²;H), make 𝑐_{𝑗,𝑘,𝑚}=𝑓𝑥,𝑦,Φ^𝑞_{𝑗,𝑘,𝑚},𝑑𝑥,𝑦_{𝑖𝑗,𝑘,𝑚}=𝑓𝑥,𝑦,Ψ^𝑞,𝑖_{𝑗,𝑘,𝑚}𝑥,𝑦𝑖=1,2,3,𝑗,𝑘,𝑚∈𝑍 (3.5)𝑑^𝑖_{𝑗,𝑘,𝑚}(𝑖=1,2,3) is called the discrete quaternion wavelet transform of 𝑓(𝑥,𝑦).

3.2. Quaternion Wavelet Transform’s Structure

The above discussion shows that quaternion wavelet transform by using four real discrete wavelet transforms, the first real discrete wavelet corresponding quaternion wavelet real part, the other real discrete wavelets are formed by the first real discrete wavelet transform by Hilbert transform, corresponding to the three imaginary parts of quaternion wavelet, respectively. In the last section, we know that Φ^𝑞(𝑥,𝑦) is the quaternion scale function, and Ψ^𝑞,1(𝑥,𝑦),Ψ^𝑞,2(𝑥,𝑦),Ψ^𝑞,3(𝑥,𝑦) are the quaternion wavelet basis functions, if quaternion scale function and wavelet function's corresponding real component are taken out to form a matrix: 𝜑𝐺=^ℎ(𝑥)𝜑_ℎ𝑦𝜑_ℎ(𝑥)𝜓_ℎ𝑦𝜓_ℎ(𝑥)𝜑_ℎ𝑦𝜓_ℎ(𝑥)𝜓_ℎ𝑦𝜑_𝑔(𝑥)𝜑_ℎ𝑦𝜑_𝑔(𝑥)𝜓_ℎ𝑦𝜓_𝑔(𝑥)𝜑_ℎ𝑦𝜓_𝑔(𝑥)𝜓_ℎ𝑦𝜑_ℎ(𝑥)𝜑_𝑔𝑦𝜑_ℎ(𝑥)𝜓_𝑔𝑦𝜓_ℎ(𝑥)𝜑_𝑔𝑦𝜓_ℎ(𝑥)𝜓_𝑔𝑦𝜑_𝑔(𝑥)𝜑_𝑔𝑦𝜑_𝑔(𝑥)𝜓_𝑔𝑦𝜓_𝑔(𝑥)𝜑_𝑔𝑦𝜓_𝑔(𝑥)𝜓_𝑔𝑦. (3.6) Then each row of the matrix 𝐺 corresponds to the one real wavelet of quaternion wavelet, the first column corresponding to quaternion wavelet scale function, the other columns are quaternion wavelet’s three wavelet functions corresponding to horizontal, vertical and diagonal three subbands. In the space 𝐿²(𝑅²), there are four standard orthogonal real wavelet bases, by wavelet frame and the concept of 2D real wavelet, we know that quaternion wavelet base in 𝐿²(𝑅²) form a tight frame, frame bound is 4.

Figure 1 shows the decomposition and reconstruction of quaternion wavelet transform and an example of a QWT decomposition, where ℎ₀ and ℎ₁ are low-pass and high-pass filter of real wavelet; 𝑔₀ and 𝑔₁ are low-pass and high-pass filter, corresponding to Hilbert transform of ℎ₀ and ℎ₁, respectively; ℎ₀ and ℎ₁are synthesis filter.

(a) Decompose structure

(b) Reconstruct structure

(c) One level QWT decomposition on the Barbara image

Figure 1 

The decomposition and reconstruction of quaternion wavelet.

3.3. The Magnitude and Phase Representation of Image

An image via quaternion wavelet transform, the coefficients can constitute a matrix 𝑄: 𝑄=𝐿𝐿_{𝜑ℎ(𝑥)𝜑ℎ(𝑦)}𝐿H_{𝜑ℎ(𝑥)𝜓ℎ(𝑦)}H𝐿_{𝜓ℎ(𝑥)𝜑ℎ(𝑦)}HH_{𝜓ℎ(𝑥)𝜓ℎ(𝑦)}𝐿𝐿_{𝜑𝑔(𝑥)𝜑ℎ(𝑦)}𝐿H_{𝜑𝑔(𝑥)𝜓ℎ(𝑦)}H𝐿_{𝜓𝑔(𝑥)𝜑ℎ(𝑦)}HH_{𝜓𝑔(𝑥)𝜓ℎ(𝑦)}𝐿𝐿_{𝜑ℎ(𝑥)𝜑𝑔(𝑦)}𝐿H_{𝜑ℎ(𝑥)𝜓𝑔(𝑦)}H𝐿_{𝜓ℎ(𝑥)𝜑𝑔(𝑦)}HH_{𝜓ℎ(𝑥)𝜓𝑔(𝑦)}𝐿𝐿_{𝜑𝑔(𝑥)𝜑𝑔(𝑦)}𝐿H_{𝜑𝑔(𝑥)𝜓𝑔(𝑦)}H𝐿_{𝜓𝑔(𝑥)𝜑𝑔(𝑦)}HH_{𝜓𝑔(𝑥)𝜓𝑔(𝑦)}. (3.7) Each column of matrix 𝑄 corresponding to a quaternion 𝑞₁=𝐿𝐿_{𝜑ℎ(𝑥)𝜑ℎ(𝑦)}+𝑖𝐿H_{𝜑ℎ(𝑥)𝜓ℎ(𝑦)}+𝑗H𝐿_{𝜓ℎ(𝑥)𝜑ℎ(𝑦)}+𝑘HH_{𝜓ℎ(𝑥)𝜓ℎ(𝑦)},𝑞₂=𝐿𝐿_{𝜑𝑔(𝑥)𝜑ℎ(𝑦)}+𝑖𝐿H_{𝜑𝑔(𝑥)𝜓ℎ(𝑦)}+𝑗H𝐿_{𝜓𝑔(𝑥)𝜑ℎ(𝑦)}+𝑘HH_{𝜓𝑔(𝑥)𝜓ℎ(𝑦)},𝑞₃=𝐿𝐿_{𝜑ℎ(𝑥)𝜑𝑔(𝑦)}+𝑖𝐿H_{𝜑ℎ(𝑥)𝜓𝑔(𝑦)}+𝑗H𝐿_{𝜓ℎ(𝑥)𝜑𝑔(𝑦)}+𝑘HH_{𝜓ℎ(𝑥)𝜓𝑔(𝑦)},𝑞₄=𝐿𝐿_{𝜑𝑔(𝑥)𝜑𝑔(𝑦)}+𝑖𝐿H_{𝜑𝑔(𝑥)𝜓𝑔(𝑦)}+𝑗H𝐿_{𝜓𝑔(𝑥)𝜑𝑔(𝑦)}+𝑘HH_{𝜓𝑔(𝑥)𝜓𝑔(𝑦)}. (3.8) The signal phase information reflects important local features; according to the concept of quaternion 𝑞=|𝑞|𝑒^𝑖𝜑𝑒^𝑗𝜃𝑒^𝑘𝜓, each quaternion wavelet coefficient can be represented using the magnitude and phase; thus, we can obtain the amplitude and phase matrix of quaternion coefficients (hereinafter referred to as magnitude and phase matrix) 𝐹=𝐿𝐿_|𝑞1|𝐿H_|𝑞2|H𝐿_|𝑞3|HH_|𝑞4|𝐿𝐿_𝜑1𝐿H_𝜑2H𝐿_𝜑3HH_𝜑4𝐿𝐿_𝜃1𝐿H_𝜃2H𝐿_𝜃3HH_𝜃4𝐿𝐿_𝜓1𝐿H_𝜓2H𝐿_𝜓3HH_𝜓4. (3.9) Assume that the size of original image is 𝑚×𝑛 and the size of magnitude and phase matrix 𝐹 is (2^2−𝑗𝑚)×(2^2−𝑗𝑛) in the 𝑗th scale QWT decomposition. Figure 2 shows an example of the magnitude and phase representation of QWT coefficients.

(a) Original image

(b) Magnitude

(c) Phase 𝜑

(d) Phase 𝜃

(e) Phase 𝜓

Figure 2 

The quaternion wavelet transform of image Barbara. (a)–(e) Original image, magnitude, and the 3 terms of phase (𝜑,𝜃,𝜓).

3.4. Quaternion Wavelet Transform Filters Design

In order to calculate the coefficients of QWT, looking from the structure, the quaternion wavelet filters’ system is similar to dual-tree complex wavelet, quaternion wavelet filters’ coefficients is quaternion, and it is realized by using the dual-tree algorithm, using an analytic quaternion wavelet bases in order to satisfy the Hilbert transform. Quaternion wavelet filters are dual-tree filters, each filters’ subtree part comprises 2 analysis filters and 2 synthesis filters, respectively. In this paper, we use the quaternion wavelet with near-symmetric orthogonal “Farrs” filters [24] at level 1 and Q-shift dual-tree filters [25] at higher levels.

4. Bayesian Image Denoising in Quaternion Wavelet Domain

4.1. Generalized Gaussian Distribution (GGD) Model

We find that the QWT coefficients of natural image is mainly distributed in the near to zero, and the two sides have a long tail, the traditional Gauss distribution is not accurate modeling of the QWT coefficients’ distribution. Through a large number of examples research, generalized Gaussian distribution (GGD) model can be used as the prior model of QWT coefficients in the high frequency subbands. The GGD [16] is 𝐺𝐺_𝜎,𝛽(𝑥)=𝐶𝜎,𝛽exp−𝛼𝜎,𝛽|𝑥|^𝛽,−∞<𝑥<∞,𝜎>0,𝛽>0, (4.1) where 𝛼(𝜎,𝛽)=𝜎⁻¹[(Γ(3/𝛽))/(Γ(1/𝛽))]^1/2,𝐶(𝜎,𝛽)=(𝛽𝛼(𝜎,𝛽))/(2Γ(1/𝛽)).

And Γ(𝑡) is the gamma function, the parameter 𝜎 is the wavelet coefficients’ standard variance, and 𝛽 is shape parameter. 𝛽 in 1, GGD model becomes Laplace distributions, when 𝛽 in 2 is Gauss distributions. Figure 3 shows the finest horizontal subbands coefficients’ probability histograms of Lena image (perform 5 levels decomposition) and corresponding GGD curves, where 𝜎 are 1.0918, 1.1924, 1.2936, and 1.3922 and 𝛽 are 0.7519, 0.8042, 0.8521, and 0.8011, respectively. As you can see from Figure 3, the GGD curves and the probability histograms are better fitted, which proved again that GGD distribution can accurately simulate the probability distribution of subbands coefficients.

(a) Real part

(b) Imaginary part 𝑖

(c) Imaginary part 𝑗

(d) Imaginary part 𝑘

Figure 3 

Probability histogram and corresponding GGD distribution for Lena image.

4.2. Bayesian Threshold

In common cases, the clear image is corrupted by zero mean additive white Gaussian noise, and the degradation model is 𝑔𝑥,𝑦=𝑓𝑥,𝑦+𝜀𝑥,𝑦, (4.2) where 𝑓(𝑥,𝑦), 𝑔(𝑥,𝑦), and 𝜀(𝑥,𝑦) represent the clean image, the observed noisy image, and the white Gaussian noise with variance 𝜎²_𝑛. Image via QWT, get 𝑤=𝑚+𝑛, (4.3) where 𝑤, 𝑚, and 𝑛 are the observed QWT coefficients, the original QWT coefficients, and the noise QWT coefficients, respectively. The goal of denoising is estimate 𝑚 from 𝑤, make 𝑚 as closer as possible to the original coefficients 𝑚. According to Bayesian theory, we can get the threshold 𝑇 [16] 𝜎𝑇=^2𝑛𝜎. (4.4) The formula (4.4) is under the frame of Bayesian theory, and it is suitable for the use of soft thresholding function for image denoising. Soulard and Carré [11, 26] deeply studied the amplitude and three phases of quaternion wavelet transform coefficients. In the image coding experiments, they only code the QWT magnitude that can obtain better visual effect than discrete wavelet transform, so they concluded that the QWT phase contains far more information than the DWT sign. Moreover, since the real and imaginary parts of QWT coefficients are not shift invariant individually but the magnitudes are. Based on the study of the distribution of magnitudes of a large number of image quaternion wavelet coefficients, statistics found that the histogram of magnitudes subtracting its mean can approximate for GGD model. Although this is an approximation, the performance of our algorithm suggests it is an acceptable assumption. Figure 4 shows the histogram of magnitudes subtracting its mean of the finest horizontal subbands of boat image (perform 4 levels decomposition), and the distribution of other natural images are similar. Therefore, we assume the magnitudes of the QWT coefficients are corrupted by additive Gaussian noise.

Figure 4 

The histogram of magnitudes subtracting its mean of boat image.

In order to achieve better denoising effect, we amended formula (4.4), multiplied by a parameter 𝛾, and obtained 𝜎𝑇=𝛾^2𝑛𝜎. (4.5) For each subbands of QWT domain, the threshold is too low, the denoising effect is not obvious, the threshold is high, and QWT coefficient excessive “killed” would lose more detail, from formula (4.5) known that optimal threshold 𝑇 mainly depends on the parameter 𝛾. Experiments show that although each image with different noise levels optimal parameter 𝛾 is disaffinity, but the parameter 𝛾 is concentrated in a small range. Figure 5 shows the optimal 𝛾 in the range between 0.6 and 1.1, and noise variance is bigger, the value 𝛾 is greater. In this paper, the denoising algorithm used 5 levels QWT decomposition; in order to better remove noise, we used the combined form of 𝛾 value, namely, the value of 𝛾 reduces gradually in the scales from coarser to finer.

Figure 5 

The average PSNR varies with 𝛾 of multiple images in different levels of noise variance.

4.3. Denoising Algorithm

First we estimate the variance of noise 𝜎_𝑛². It can be estimated by the finest subbands coefficients’ magnitude 𝜎_𝑛=Median𝑤_𝑖,𝑗,0.6745 (4.6) where 𝑤_𝑖,𝑗 is the finest subbands coefficient and |𝑤_𝑖,𝑗| denote its magnitude. Then calculate the marginal variance 𝜎_𝑚², because of 𝜎_𝑤²=𝜎_𝑚²+𝜎_𝑛². (4.7) Here 𝜎_𝑤² is the noisy wavelet coefficient magnitude’s variance. We will estimate the variance 𝜎_𝑤² by a square-shaped neighborhood window 𝑊(𝑘) that is centered at |𝑤_𝑗|: 𝜎_𝑤²=1𝑁²_{𝑤𝑗∈𝑊(𝑘)}𝑤_𝑗², (4.8) where 𝑁² is the size of 𝑊(𝑘). Repeated experiments show that 𝑊(𝑘) take 5×5 square-shaped neighborhood window will be more suitable. Thus, we can get the estimation for 𝜎_𝑚²: 𝜎_𝑚²=𝜎_𝑤²−𝜎_𝑛²₊. (4.9) Below, we summarize the main steps taken by the proposed image denoising method in this paper.

Step 1. Perform the 5 levels QWT decomposition on the noisy image.

Step 2. Calculate the magnitude and phase of each subbands coefficients, and using (4.6) and (4.9) calculate 𝜎_𝑛² and 𝜎_𝑚² values in each subbands.

Step 3. Compute noise-free coefficients’ magnitude using the multiplied threshold and soft thresholding function.

Step 4. Perform inverse transform on the estimated magnitude and the original phase to get QWT coefficients, and conduct the inverse QWT with estimated coefficients, and get the denoised image.

4.4. Experimental Results and Discussion

In the experiment, we have tested various denoising method for a representative set of standard 8-bit grayscale images such as Lena, Barbara, Boat, Flinstones (size 512×512), and House, Mit, and Cameraman (size 512×512) images, corrupted by simulated additive Gauss white noise with zero mean, variance 𝜎_𝑛². We compared the proposed method with several other state-of-the-art outstanding denoising methods, namely, Donoho’s hard threshold (“db8” wavelet) [13], the LAWML [15] in dual-tree complex wavelet domain, SURE-LET [27], non-Gaussian bivariate distribution model in PDTDFB [28] domain (P-Bivariate), and we evaluate denoising performance combining with PSNR and visual effects in this paper. The PSNR is defined as PSNR=20log₁₀255,MSE (4.10) where MSE=(1/𝑁²)(𝑋−𝑌)², and 𝑋, 𝑌 are the original image and the denoised image, respectively, 𝑁² is image size.

The PSNR values are listed in Table 1, and denoised images are given in Figures 6, 7, 8, and 9.

Table 1 

Performance comparison of different denoising methods (PSNR).

(a) Noisy image

(b) Hard

(c) LAWML

(d) SURE-LET

(e) P-Bivariate

(f) Our method

Figure 6 

Denoised image using different denoising methods (𝜎𝑛=20, Lena).

(a) Noisy image

(b) Hard

(c) LAWML

(d) SURE-LET

(e) P-Bivariate

(f) Our method

Figure 7 

Denoised image using different denoising methods (𝜎𝑛=20, Barbara).

(a) Noisy image

(b) Hard

(c) LAWML

(d) SURE-LET

(e) P-Bivariate

(f) Our method.

Figure 8 

Denoised image using different denoising methods (𝜎𝑛=20, Flintstones).

(a) Noisy image

(b) Hard

(c) LAWML

(d) SURE-LET

(e) P-Bivariate

(f) Our method

Figure 9 

Denoised image using different denoising methods (𝜎𝑛=20, House).

From the experimental data and the denoised images, it can be seen that the proposed method almost provides the highest objective data. The proposed method has obvious advantage over hard thresholding shrinkage, and the PSNR has been greatly improved, compared with LAWML and SURE-LET which also has a certain degree of increase; the denoised Barbara image of P-Bivariate method obtained higher PSNR values than our method; it is mainly due to that PDTDFB has multidirection that it can well represent more texture images. From that Figures 6–9, we can see that hard threshold denoised images’ fuzzy phenomenon is very serious; LAWML and SURE-LET methods cannot effectivly remove the noise, and there is still a lot of noise left in the denoised images; denoised images add a lot of false information using the P-Bivariate denoising method; due to that the QWT phase can compensate the loss of information, thus the proposed denoising method retains more detail information and obtains better visual effect.

5. Conclusions

Quaternion wavelet transform is established based on the quaternion algebra, quaternion Fourier transform, and Hilbert transform; using four real discrete wavelet transform (DWT), the first real discrete wavelet corresponding to quaternion wavelet real part, the other real discrete wavelet is obtained by the first real discrete wavelet transform’s Hilbert transform, corresponding to quaternion wavelet three imaginary part, respectively, the four real wavelet composed of quaternion analytic signal. It can be understood as the improved real wavelet and complex wavelet’s promotion, which have approximate shift invariance, abundant phase information, and limited redundancy and so forth, while still retaining the traditional wavelet time-frequency localization ability, filters design using Hilbert transform pair of dual-tree structure, and it is easy to be realized.

This paper mainly studies some of the concepts and properties of quaternion wavelet transform, gives quaternion wavelet scale and wavelet functions, and applies the quaternion wavelet in image denoising, puts forward Bayesian denoising method based on quaternion wavelet transform, considering wavelet coefficient’s correlation, and generalized Gaussian distribution is used to model the probability distribution function of wavelet coefficients’ magnitude and the best range of the Bayesian thresholding parameter is found out. The experimental results show that our method both in visual effect and PSNR are better than many current denoising methods.

Acknowledgments

The authors are grateful for the support received from the National Natural Science Foundation of China (NSFC no. 11172086), and this work is supported by the Natural Science Foundation of Anhui (no. 11040606M06).