Quaternion Wavelet Analysis and Application in Image Denoising

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 L2 R2 , proves and presents quaternion wavelet’s scale basis function and wavelet basis function concepts in spatial scale space L2 R2;H , and studies quaternion wavelet transform structure. Finally, the quaternionwavelet 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.


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 twodimensional 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 dualtree 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: g x, y f x, y ε x, y , where f x, y , g x, y , and ε x, y represent the clean image, the observed noisy image, and the white Gaussian noise with variance σ 2 n .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.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 L 2 R 2 , and then in the space L 2 R 2 ; 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.

2D Hilbert Transform Concept
Definition 2.1 see 8 .Set f x, y ∈ L 2 R 2 , then f x, y along the x-axis and the y-axis and along the x, y axis of the Hilbert transform f Hx x, y , f Hy x, y , and f Hxy x, y , respectively: the corresponding the frequency domain is f Hxy u, v − sgn u sgn v f u, v .

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 {q q r iq i jq j kq k | q r , q i , q j , q k ∈ R}.Three imaginary i, j, k satisfy the following rules: And define q as the quaternion.
The conjugate of quaternion q is given by q Re q − Im q q r − iq i − jq j − kq k and the module of quaternion q is given by |q| qq q 2 r q 2 i q 2 j q 2 k .In addition, quaternion q can also express as: q |q|e iϕ e jθ e kψ , where |q| is the modulus of q, ϕ, θ, ψ are the-three phase angles which is uniquely defined within the range ϕ, θ, ψ 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 f x, y is a real two-dimensional signals, then quaternion analytic signal can be defined as f q x, y f x, y if Hx x, y jf Hy x, y kf Hxy x, y .where f Hx x, y , f Hy x, y , f Hxy x, y are the f x, y along the x-axis and the y-axis and along the x, y axis of the Hilbert transform.

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.

Quaternion Wavelet's Scale Basis and Wavelet Basis
Proof.According to the theorem of the paper in 23 , we know that standard orthogonal basis in the frequency domain is equivalent to , 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 {V j } j∈z be a 1D orthogonal multiresolution analysis, the corresponding scale function and wavelet function are ϕ h x and ψ h x , respectively; note that H x ϕ h x ϕ h y ϕ g x ϕ h y H x ϕ h x ϕ h y are said to function as ϕ h x ϕ h y along x-axis direction Hilbert transform, and the corresponding mark is similar H y ϕ h x ϕ h y ϕ h x ϕ g y ,
H xy ϕ h x ϕ h y ϕ g x ϕ g y .

3.1
So ) are the standard orthogonal base of the space V 0 V 0 and the space V 0 V 0 , where V 0 HV 0 .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 x, y and x, y axis direction's Hilbert transform, so by Theorem 3.2 one knows {ϕ In fact, by Theorem 3.2 and the above analysis, it is easy to prove the following.
Theorem 3.3.Let {V j } j∈z be a 1D L 2 R orthogonal multiresolution analysis, ϕ h x , ψ h x are scale function and wavelet function of corresponding respectively.Note Φ q x, y ϕ h x ϕ h y iϕ g x ϕ h y jϕ h x ϕ g y kϕ g x ϕ g y , Φ q j,k,m x, y ϕ h,j,k x ϕ h,j,m y iϕ g,j,k x ϕ h,j,m y jϕ h,j,k x ϕ g,j,m y kϕ g,j,k x ϕ g,j,m y ,

3.2
where is the basis function of quaternion wavelet scale space.And called Φ q x, y is quaternion wavelet scale function in L 2 R 2 ; H ; called {Φ Further one has the following.Theorem 3.4.Let {V j } j∈z be a 1D orthogonal L 2 R multiresolution analysis, and ϕ h x , ψ h x are the scale function and wavelet function of corresponding, respectively; note that Ψ q,1 x, y ϕ h x ψ h y iϕ g x ψ h y jϕ h x ψ g y kϕ g x ψ g y , Ψ q,2 x, y ψ h x ϕ h y iψ g x ϕ h y jψ h x ϕ g y kψ g x ϕ g y , Ψ q,3 x, y ψ h x ψ h y iψ g x ψ h y jψ h x ψ g y kψ g x ψ g y .

3.3
The shift and expand of Ψ q,1 x, y , Ψ q,2 x, y and Ψ q,3 x, y form of function family x, y ψ h,j,k x ψ h,j,m y iψ g,j,k x ψ h,j,m y jψ h,j,k x ψ g,j,m y kψ g,j,k x ψ g,j,m y k,m∈z ,

3.4
where ψ h,j,k x 2 −j/2 ψ h 2 −j x − k , j, k, m ∈ Z, the corresponding mark is similar.And Ψ q,1 x, y , Ψ q,2 x, y , Ψ q,3 x, y as the quaternion wavelet basis functions in L 2 R 2 ; H space; that {Ψ are discrete quaternion wavelet functions in L 2 R 2 ; H) space.

One can get concept of the quaternion wavelet transform.
Mathematical Problems in Engineering 7 Definition 3.5.For for all f x, y ∈ L 2 R 2 ; H , make 3 is called the discrete quaternion wavelet transform of f x, y .

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 Φ q x, y is the quaternion scale function, and Ψ q,1 x, y , Ψ q,2 x, y , Ψ q,3 x, y are the quaternion wavelet basis functions, if quaternion scale function and wavelet function's corresponding real component are taken out to form a matrix:

ϕ h y ϕ h x ψ h y ψ h x ϕ h y ψ h x ψ h y ϕ g x ϕ h y ϕ g x ψ h y ψ g x ϕ h y ψ g x ψ h y ϕ h x ϕ g y ϕ h x ψ g y ψ h x ϕ g y ψ h
x ψ g y ϕ g x ϕ g y ϕ g x ψ g y ψ g x ϕ g y ψ g x ψ g y Then each row of the matrix G 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 L 2 R 2 , 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 L 2 R 2 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 h 0 and h 1 are low-pass and high-pass filter of real wavelet; g 0 and g 1 are low-pass and high-pass filter, corresponding to Hilbert transform of h 0 and h 1 , respectively; h 0 and h 1 are synthesis filter.

The Magnitude and Phase Representation of Image
An image via quaternion wavelet transform, the coefficients can constitute a matrix Q:

LL ϕ h x ϕ h y LH ϕ h x ψ h y HL ψ h x ϕ h y HH ψ h x ψ h y LL ϕ g x ϕ h y LH ϕ g x ψ h y HL ψ g x ϕ h y HH ψ g x ψ h y
LL ϕ h x ϕ g y LH ϕ h x ψ g y HL ψ h x ϕ g y HH ψ h x ψ g y LL ϕ g x ϕ g y LH ϕ g x ψ g y HL ψ g x ϕ g y HH ψ g x ψ g y ⎞ ⎟ ⎟ ⎠ . 3.7 Imaginary tree (i) Imaginary tree (j)

Each column of matrix Q corresponding to a quaternion q 1 LL ϕ h x ϕ h y iLH ϕ h x ψ h y jHL ψ h x ϕ h y kHH ψ h x ψ h y ,
q 2 LL ϕ g x ϕ h y iLH ϕ g x ψ h y jHL ψ g x ϕ h y kHH ψ g x ψ h y , q 3 LL ϕ h x ϕ g y iLH ϕ h x ψ g y jHL ψ h x ϕ g y kHH ψ h x ψ g y , q 4 LL ϕ g x ϕ g y iLH ϕ g x ψ g y jHL ψ g x ϕ g y kHH ψ g x ψ g y .

3.8
The signal phase information reflects important local features; according to the concept of quaternion q |q|e iϕ e jθ e kψ , 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 Assume that the size of original image is m × n and the size of magnitude and phase matrix F is 2 2−j m × 2 2−j n in the jth scale QWT decomposition. Figure 2 shows an example of the magnitude and phase representation of QWT coefficients.

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.

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 where α σ, β σ −1 Γ 3/β / Γ 1/β 1/2 , C σ, β βα σ, β / 2Γ 1/β .And Γ t 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.

Bayesian Threshold
In common cases, the clear image is corrupted by zero mean additive white Gaussian noise, and the degradation model is g x, y f x, y ε x, y , 4.2 where w, m, and n are the observed QWT coefficients, the original QWT coefficients, and the noise QWT coefficients, respectively.The goal of denoising is estimate m from w, make m as closer as possible to the original coefficients m.According to Bayesian theory, we can get the threshold T 16 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 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.
In order to achieve better denoising effect, we amended formula 4.4 , multiplied by a parameter γ, and obtained 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 T 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.Here σ w 2 is the noisy wavelet coefficient magnitude's variance.We will estimate the variance σ w 2 by a square-shaped neighborhood window W k that is centered at |w j |:

Denoising Algorithm
where N 2 is the size of W k .Repeated experiments show that W k take 5×5 square-shaped neighborhood window will be more suitable.Thus, we can get the estimation for σ m 2 : 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 σ n 2 and σ m 2 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.

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 σ n 2 .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 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.

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 Mathematical Problems in Engineering 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.

+c
One level QWT decomposition on the Barbara image

Figure 1 :
Figure 1: The decomposition and reconstruction of quaternion wavelet.

Figure 2 :Figure 3 :
Figure 2: The quaternion wavelet transform of image Barbara.a -e Original image, magnitude, and the 3 terms of phase ϕ, θ, ψ .

Figure 4 :
Figure 4: The histogram of magnitudes subtracting its mean of boat image.

Figure 5 :
Figure 5: The average PSNR varies with γ of multiple images in different levels of noise variance.

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
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, Mihc ¸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 , and non-Gaussian model, but parameters' estimation is more complex in image denoising process.

Table 1 :
Performance comparison of different denoising methods PSNR .
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.