# Quaternion Wavelet Analysis and Application in Image Denoising

Authors: Ming Yin; Wei Liu; Jun Shui; Jiangmin Wu
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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 spatialL2(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 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.


---

## Body

## 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: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 σn2. 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 spatialL2(R2), and then in the space L2(R2;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]).
Setf(x,y)∈L2(R2), then f(x,y) along the x-axis and the y-axis and along the x,y axis of the Hilbert transform fHx(x,y), fHy(x,y), and fHxy(x,y), respectively:
(2.1)fHx(x,y)=1π∫Rf(ξ1,y)x-ξ1dξ1,
the corresponding frequency domain is fHx(u,v)^=-jsgn(u)f(u,v)^;
(2.2)fHy(x,y)=1π∫Rf(x,ξ2)y-ξ2dξ2,
the corresponding frequency domain is fHy(u,v)^=-jsgn(v)f(u,v)^;
(2.3)fHxy(x,y)=1π2∬R2f(ξ1,ξ2)(x-ξ1)(y-ξ2)dξ1dξ2,
the corresponding the frequency domain is fHxy(u,v)^=-sgn(u)sgn(v)f(u,v)^.


### 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.SetH={q=qr+iqi+jqj+kqk∣qr,qi,qj,qk∈R}. Three imaginary i,j,k satisfy the following rules:
(2.4)i2=j2=k2=-1,ij=-ji=k,jk=-kj=i,ki=-ik=j.
And define q as the quaternion.The conjugate of quaternionq is given by q-=Re(q)-Im(q)=qr-iqi-jqj-kqk and the module of quaternion q is given by |q|=qq-=qr2+qi2+qj2+qk2. In addition, quaternion q can also express as: q=|q|eiφejθekψ, where |q| is the modulus of q, (φ,θ,ψ) are the-three phase angles which is uniquely defined within the range (φ,θ,ψ)∈[-π,π]×[-π/2,π/2]×[-π/4,π/4][9].SoL2(R2;H)={f:R2→H,∫R2|f(x)|2d2x<∞}, forallf,g∈L2(R2;H), define the inner product 〈f,g〉=∫R2f(x)g(x)-d2x, it is not difficult to prove L2(R2;H) is a Hilbert space.Quaternion wavelet transform is based on quaternion analytic signal, and we give the concept.Definition 2.2 (see [8]).
Setf(x,y) is a real two-dimensional signals, then quaternion analytic signal can be defined as fq(x,y)=f(x,y)+ifHx(x,y)+jfHy(x,y)+kfHxy(x,y).
wherefHx(x,y),fHy(x,y),fHxy(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.


## 2.1. 2D Hilbert Transform Concept

Definition 2.1 (see [8]).
Setf(x,y)∈L2(R2), then f(x,y) along the x-axis and the y-axis and along the x,y axis of the Hilbert transform fHx(x,y), fHy(x,y), and fHxy(x,y), respectively:
(2.1)fHx(x,y)=1π∫Rf(ξ1,y)x-ξ1dξ1,
the corresponding frequency domain is fHx(u,v)^=-jsgn(u)f(u,v)^;
(2.2)fHy(x,y)=1π∫Rf(x,ξ2)y-ξ2dξ2,
the corresponding frequency domain is fHy(u,v)^=-jsgn(v)f(u,v)^;
(2.3)fHxy(x,y)=1π2∬R2f(ξ1,ξ2)(x-ξ1)(y-ξ2)dξ1dξ2,
the corresponding the frequency domain is fHxy(u,v)^=-sgn(u)sgn(v)f(u,v)^.


## 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.SetH={q=qr+iqi+jqj+kqk∣qr,qi,qj,qk∈R}. Three imaginary i,j,k satisfy the following rules:
(2.4)i2=j2=k2=-1,ij=-ji=k,jk=-kj=i,ki=-ik=j.
And define q as the quaternion.The conjugate of quaternionq is given by q-=Re(q)-Im(q)=qr-iqi-jqj-kqk and the module of quaternion q is given by |q|=qq-=qr2+qi2+qj2+qk2. In addition, quaternion q can also express as: q=|q|eiφejθekψ, where |q| is the modulus of q, (φ,θ,ψ) are the-three phase angles which is uniquely defined within the range (φ,θ,ψ)∈[-π,π]×[-π/2,π/2]×[-π/4,π/4][9].SoL2(R2;H)={f:R2→H,∫R2|f(x)|2d2x<∞}, forallf,g∈L2(R2;H), define the inner product 〈f,g〉=∫R2f(x)g(x)-d2x, it is not difficult to prove L2(R2;H) is a Hilbert space.Quaternion wavelet transform is based on quaternion analytic signal, and we give the concept.Definition 2.2 (see [8]).
Setf(x,y) is a real two-dimensional signals, then quaternion analytic signal can be defined as fq(x,y)=f(x,y)+ifHx(x,y)+jfHy(x,y)+kfHxy(x,y).
wherefHx(x,y),fHy(x,y),fHxy(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.


## 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{φ(x-k)}k∈Z be the space of V⊂L2(R) standard orthogonal basis, set φg(x)=H(φ(x)) (i.e., that φ(x) of Hilbert transform), so {φg(x-k)}k∈Z is the standard orthogonal basis of space V~=HV⊂L2(R).Proof.
According to the theorem of the paper in [23], we know that standard orthogonal basis in the frequency domain is equivalent to: ∑k∈Z|φ(ω+2kπ)|^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{Vj}j∈z be a 1D orthogonal multiresolution analysis, the corresponding scale function and wavelet function are φh(x) and ψh(x), respectively; note that Hx(φh(x)φh(y))=φg(x)φh(y) (Hx(φ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
(3.1)Hy(φh(x)φh(y))=φh(x)φg(y),Hxy(φh(x)φh(y))=φg(x)φg(y).
So {φg(x-l)φh(y-k)}l,k∈Z (or {φh(x-l)φg(y-k)}l,k∈Z) are the standard orthogonal base of the space V0⨂V~0 and the space V~0⨂V0, where V~0=HV0.
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 thex,y and x,y axis direction’s Hilbert transform, so by Theorem 3.2 one knows {φg(x-l)φh(y-k)}l,k∈Z (or {φh(x-l)φg(y-k)}l,k∈Z; total of 4) is the standard orthogonal basis of scale space L2(R2), {φg(x-l)ψh(y-k)}l,k∈Z({φh(x-l)ψg(y-k)}l,k∈Z; total of 12) is the standard orthogonal base of wavelet space L2(R2).
In fact, by Theorem3.2 and the above analysis, it is easy to prove the following.Theorem 3.3.
Let{Vj}j∈z  be a 1D L2(R) orthogonal multiresolution analysis, φh(x),ψh(x) are scale function and wavelet function of corresponding respectively.
Note(3.2)Φq(x,y)=φh(x)φh(y)+iφg(x)φh(y)+jφh(x)φg(y)+kφg(x)φg(y),Φj,k,mq(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),
where φh,j,k(x)=2-j/2φh(2-jx-k),j,k,m∈Z, corresponding mark is similar.
So{Φj,k,mq(x,y)}k,m∈Z is the basis function of quaternion wavelet scale space. And called Φq(x,y) is quaternion wavelet scale function in L2(R2;H); called {Φj,k,mq(x,y)}k,m∈Z is discrete quaternion wavelet scale function in L2(R2;H).Further one has the following.Theorem 3.4.
Let{Vj}j∈z be a 1D orthogonal L2(R) multiresolution analysis, and φh(x),ψh(x) are the scale function and wavelet function of corresponding, respectively; note that
(3.3)Ψ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).
The shift and expand of  Ψq,1(x,y),  Ψq,2(x,y) and Ψq,3(x,y) form of function family
(3.4){Ψj,k,mq,2(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)Ψj,k,mq,2(x,y)+kψg,j,k(x)φg,j,m(y)Ψj,k,mq,2}k,m∈z,{Ψj,k,mq,1(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)Ψj,k,mq,2(x,y)+kφg,j,k(x)ψg,j,m(y)Ψj,k,mq,2}k,m∈z,{Ψj,k,mq,3(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)Ψj,k,mq,2(x,y)+kψg,j,k(x)ψg,j,m(y)Ψj,k,mq,2}k,m∈z,
where ψh,j,k(x)=2-j/2ψh(2-jx-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 L2(R2;H) space; that {Ψj,k,mq,1,Ψj,k,mq,2,Ψj,k,mq,3}j,k,m∈z are discrete quaternion wavelet functions in L2(R2;H) space.
One can get concept of the quaternion wavelet transform.Definition 3.5.
For for allf(x,y)∈L2(R2;H), make
(3.5)cj,k,m=(f(x,y),Φj,k,mq(x,y)),dj,k,mi=(f(x,y),Ψj,k,mq,i(x,y))(i=1,2,3,j,k,m∈Z)dj,k,mi(i=1,2,3) is called the discrete quaternion wavelet transform of f(x,y).


### 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Φ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:
(3.6)G=(φh(x)φ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 L2(R2), 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 L2(R2) form a tight frame, frame bound is 4.Figure1 shows the decomposition and reconstruction of quaternion wavelet transform and an example of a QWT decomposition, where h0 and h1 are low-pass and high-pass filter of real wavelet; g0 and g1 are low-pass and high-pass filter, corresponding to Hilbert transform of h0 and h1, respectively; h~0 and h~1are synthesis filter.The decomposition and reconstruction of quaternion wavelet.
(a)
Decompose structure
(b)
Reconstruct structure
(c)
One level QWT decomposition on the Barbara image


### 3.3. The Magnitude and Phase Representation of Image

An image via quaternion wavelet transform, the coefficients can constitute a matrixQ:
(3.7)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)).
Each column of matrix Q corresponding to a quaternion
(3.8)q1=LLφh(x)φh(y)+iLHφh(x)ψh(y)+jHLψh(x)φh(y)+kHHψh(x)ψh(y),q2=LLφg(x)φh(y)+iLHφg(x)ψh(y)+jHLψg(x)φh(y)+kHHψg(x)ψh(y),q3=LLφh(x)φg(y)+iLHφh(x)ψg(y)+jHLψh(x)φg(y)+kHHψh(x)ψg(y),q4=LLφg(x)φg(y)+iLHφg(x)ψg(y)+jHLψg(x)φg(y)+kHHψg(x)ψg(y).
The signal phase information reflects important local features; according to the concept of quaternion q=|q|eiφejθekψ, 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)
(3.9)F=(LL|q1|LH|q2|HL|q3|HH|q4|LLφ1LHφ2HLφ3HHφ4LLθ1LHθ2HLθ3HHθ4LLψ1LHψ2HLψ3HHψ4).
Assume that the size of original image is m×n and the size of magnitude and phase matrix F is (22-jm)×(22-jn) in the jth scale QWT decomposition. Figure 2 shows an example of the magnitude and phase representation of QWT coefficients.The quaternion wavelet transform of image Barbara. (a)–(e) Original image, magnitude, and the 3 terms of phase(φ,θ,ψ).
(a)
Original image
(b)
Magnitude
(c)
Phaseφ
(d)
Phaseθ
(e)
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.


## 3.1. Quaternion Wavelet’s Scale Basis and Wavelet Basis

Lemma 3.1.
Let{φ(x-k)}k∈Z be the space of V⊂L2(R) standard orthogonal basis, set φg(x)=H(φ(x)) (i.e., that φ(x) of Hilbert transform), so {φg(x-k)}k∈Z is the standard orthogonal basis of space V~=HV⊂L2(R).Proof.
According to the theorem of the paper in [23], we know that standard orthogonal basis in the frequency domain is equivalent to: ∑k∈Z|φ(ω+2kπ)|^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{Vj}j∈z be a 1D orthogonal multiresolution analysis, the corresponding scale function and wavelet function are φh(x) and ψh(x), respectively; note that Hx(φh(x)φh(y))=φg(x)φh(y) (Hx(φ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
(3.1)Hy(φh(x)φh(y))=φh(x)φg(y),Hxy(φh(x)φh(y))=φg(x)φg(y).
So {φg(x-l)φh(y-k)}l,k∈Z (or {φh(x-l)φg(y-k)}l,k∈Z) are the standard orthogonal base of the space V0⨂V~0 and the space V~0⨂V0, where V~0=HV0.
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 thex,y and x,y axis direction’s Hilbert transform, so by Theorem 3.2 one knows {φg(x-l)φh(y-k)}l,k∈Z (or {φh(x-l)φg(y-k)}l,k∈Z; total of 4) is the standard orthogonal basis of scale space L2(R2), {φg(x-l)ψh(y-k)}l,k∈Z({φh(x-l)ψg(y-k)}l,k∈Z; total of 12) is the standard orthogonal base of wavelet space L2(R2).
In fact, by Theorem3.2 and the above analysis, it is easy to prove the following.Theorem 3.3.
Let{Vj}j∈z  be a 1D L2(R) orthogonal multiresolution analysis, φh(x),ψh(x) are scale function and wavelet function of corresponding respectively.
Note(3.2)Φq(x,y)=φh(x)φh(y)+iφg(x)φh(y)+jφh(x)φg(y)+kφg(x)φg(y),Φj,k,mq(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),
where φh,j,k(x)=2-j/2φh(2-jx-k),j,k,m∈Z, corresponding mark is similar.
So{Φj,k,mq(x,y)}k,m∈Z is the basis function of quaternion wavelet scale space. And called Φq(x,y) is quaternion wavelet scale function in L2(R2;H); called {Φj,k,mq(x,y)}k,m∈Z is discrete quaternion wavelet scale function in L2(R2;H).Further one has the following.Theorem 3.4.
Let{Vj}j∈z be a 1D orthogonal L2(R) multiresolution analysis, and φh(x),ψh(x) are the scale function and wavelet function of corresponding, respectively; note that
(3.3)Ψ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).
The shift and expand of  Ψq,1(x,y),  Ψq,2(x,y) and Ψq,3(x,y) form of function family
(3.4){Ψj,k,mq,2(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)Ψj,k,mq,2(x,y)+kψg,j,k(x)φg,j,m(y)Ψj,k,mq,2}k,m∈z,{Ψj,k,mq,1(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)Ψj,k,mq,2(x,y)+kφg,j,k(x)ψg,j,m(y)Ψj,k,mq,2}k,m∈z,{Ψj,k,mq,3(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)Ψj,k,mq,2(x,y)+kψg,j,k(x)ψg,j,m(y)Ψj,k,mq,2}k,m∈z,
where ψh,j,k(x)=2-j/2ψh(2-jx-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 L2(R2;H) space; that {Ψj,k,mq,1,Ψj,k,mq,2,Ψj,k,mq,3}j,k,m∈z are discrete quaternion wavelet functions in L2(R2;H) space.
One can get concept of the quaternion wavelet transform.Definition 3.5.
For for allf(x,y)∈L2(R2;H), make
(3.5)cj,k,m=(f(x,y),Φj,k,mq(x,y)),dj,k,mi=(f(x,y),Ψj,k,mq,i(x,y))(i=1,2,3,j,k,m∈Z)dj,k,mi(i=1,2,3) is called the discrete quaternion wavelet transform of f(x,y).


## 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Φ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:
(3.6)G=(φh(x)φ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 L2(R2), 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 L2(R2) form a tight frame, frame bound is 4.Figure1 shows the decomposition and reconstruction of quaternion wavelet transform and an example of a QWT decomposition, where h0 and h1 are low-pass and high-pass filter of real wavelet; g0 and g1 are low-pass and high-pass filter, corresponding to Hilbert transform of h0 and h1, respectively; h~0 and h~1are synthesis filter.The decomposition and reconstruction of quaternion wavelet.
(a)
Decompose structure
(b)
Reconstruct structure
(c)
One level QWT decomposition on the Barbara image


## 3.3. The Magnitude and Phase Representation of Image

An image via quaternion wavelet transform, the coefficients can constitute a matrixQ:
(3.7)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)).
Each column of matrix Q corresponding to a quaternion
(3.8)q1=LLφh(x)φh(y)+iLHφh(x)ψh(y)+jHLψh(x)φh(y)+kHHψh(x)ψh(y),q2=LLφg(x)φh(y)+iLHφg(x)ψh(y)+jHLψg(x)φh(y)+kHHψg(x)ψh(y),q3=LLφh(x)φg(y)+iLHφh(x)ψg(y)+jHLψh(x)φg(y)+kHHψh(x)ψg(y),q4=LLφg(x)φg(y)+iLHφg(x)ψg(y)+jHLψg(x)φg(y)+kHHψg(x)ψg(y).
The signal phase information reflects important local features; according to the concept of quaternion q=|q|eiφejθekψ, 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)
(3.9)F=(LL|q1|LH|q2|HL|q3|HH|q4|LLφ1LHφ2HLφ3HHφ4LLθ1LHθ2HLθ3HHθ4LLψ1LHψ2HLψ3HHψ4).
Assume that the size of original image is m×n and the size of magnitude and phase matrix F is (22-jm)×(22-jn) in the jth scale QWT decomposition. Figure 2 shows an example of the magnitude and phase representation of QWT coefficients.The quaternion wavelet transform of image Barbara. (a)–(e) Original image, magnitude, and the 3 terms of phase(φ,θ,ψ).
(a)
Original image
(b)
Magnitude
(c)
Phaseφ
(d)
Phaseθ
(e)
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
(4.1)GGσ,β(x)=C(σ,β)exp{-[α(σ,β)|x|]β},-∞<x<∞,σ>0,β>0,
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.Probability histogram and corresponding GGD distribution for Lena image.
(a)
Real part
(b)
Imaginary parti
(c)
Imaginary partj
(d)
Imaginary partk


### 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)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 σn2. Image via QWT, get
(4.3)w=m+n,
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]
(4.4)T=σn2σ.
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
(4.5)T=γσn2σ.
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.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σn2. It can be estimated by the finest subbands coefficients’ magnitude
(4.6)σn=Median(|wi,j|)0.6745,
where wi,j is the finest subbands coefficient and |wi,j| denote its magnitude. Then calculate the marginal variance σm2, because of
(4.7)σw2=σm2+σn2.
Here σw2 is the noisy wavelet coefficient magnitude’s variance. We will estimate the variance σw2 by a square-shaped neighborhood window W(k) that is centered at |wj|:
(4.8)σw2=1N2∑|wj|∈W(k)(|wj|)2,
where N2 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 σm2:
(4.9)σm2=(σw2-σn2)+.
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 σn2 and σm2 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 (size512×512), and House, Mit, and Cameraman (size 512×512) images, corrupted by simulated additive Gauss white noise with zero mean, variance σn2. 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
(4.10)PSNR=20log10(255MSE),
where MSE=(1/N2)(X-Y)2, and X, Y are the original image and the denoised image, respectively, N2 is image size.The PSNR values are listed in Table1, and denoised images are given in Figures 6, 7, 8, and 9.Table 1
Performance comparison of different denoising methods (PSNR).
σ
Donoho’s HT
LAWML
SURE-LET
P-Bivariate
Our method
10
30.05
34.59
34.56
34.77
35.18
Lena
20
27.29
31.32
31.37
31.72
32.10
30
25.85
29.28
29.56
29.89
30.18
10
26.02
32.28
32.16
33.16
33.34
Barbara
20
23.38
28.81
27.96
29.75
29.36
30
22.38
26.83
25.82
27.77
27.07
10
27.42
32.11
32.91
32.72
33.15
Boat
20
25.08
29.04
29.47
29.66
29.88
30
23.83
27.18
27.63
27.83
27.92
10
24.57
30.54
31.17
30.51
31.45
Flinstones
20
21.72
27.17
27.40
27.28
28.12
30
19.95
25.02
25.29
25.35
26.04
10
30.07
33.18
34.27
34.28
34.75
House
20
27.09
29.98
30.90
30.89
31.76
30
25.49
27.96
28.96
28.94
29.88
10
26.09
30.92
32.07
31.63
32.87
Mit
20
22.78
27.08
27.88
27.77
28.62
30
20.95
24.91
25.60
25.60
26.22
10
26.87
31.47
32.35
31.81
33.08
Cameraman
20
24.05
27.94
28.51
28.27
29.13
30
22.54
26.00
26.48
26.38
27.05Denoised image using different denoising methods (σn=20, Lena).
(a)
Noisy image
(b)
Hard
(c)
LAWML
(d)
SURE-LET
(e)
P-Bivariate
(f)
Our methodDenoised image using different denoising methods (σn=20, Barbara).
(a)
Noisy image
(b)
Hard
(c)
LAWML
(d)
SURE-LET
(e)
P-Bivariate
(f)
Our methodDenoised image using different denoising methods (σn=20, Flintstones).
(a)
Noisy image
(b)
Hard
(c)
LAWML
(d)
SURE-LET
(e)
P-Bivariate
(f)
Our method.Denoised image using different denoising methods (σn=20, House).
(a)
Noisy image
(b)
Hard
(c)
LAWML
(d)
SURE-LET
(e)
P-Bivariate
(f)
Our methodFrom 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 Figures6–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.


## 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
(4.1)GGσ,β(x)=C(σ,β)exp{-[α(σ,β)|x|]β},-∞<x<∞,σ>0,β>0,
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.Probability histogram and corresponding GGD distribution for Lena image.
(a)
Real part
(b)
Imaginary parti
(c)
Imaginary partj
(d)
Imaginary partk


## 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)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 σn2. Image via QWT, get
(4.3)w=m+n,
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]
(4.4)T=σn2σ.
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
(4.5)T=γσn2σ.
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.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σn2. It can be estimated by the finest subbands coefficients’ magnitude
(4.6)σn=Median(|wi,j|)0.6745,
where wi,j is the finest subbands coefficient and |wi,j| denote its magnitude. Then calculate the marginal variance σm2, because of
(4.7)σw2=σm2+σn2.
Here σw2 is the noisy wavelet coefficient magnitude’s variance. We will estimate the variance σw2 by a square-shaped neighborhood window W(k) that is centered at |wj|:
(4.8)σw2=1N2∑|wj|∈W(k)(|wj|)2,
where N2 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 σm2:
(4.9)σm2=(σw2-σn2)+.
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 σn2 and σm2 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 (size512×512), and House, Mit, and Cameraman (size 512×512) images, corrupted by simulated additive Gauss white noise with zero mean, variance σn2. 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
(4.10)PSNR=20log10(255MSE),
where MSE=(1/N2)(X-Y)2, and X, Y are the original image and the denoised image, respectively, N2 is image size.The PSNR values are listed in Table1, and denoised images are given in Figures 6, 7, 8, and 9.Table 1
Performance comparison of different denoising methods (PSNR).
σ
Donoho’s HT
LAWML
SURE-LET
P-Bivariate
Our method
10
30.05
34.59
34.56
34.77
35.18
Lena
20
27.29
31.32
31.37
31.72
32.10
30
25.85
29.28
29.56
29.89
30.18
10
26.02
32.28
32.16
33.16
33.34
Barbara
20
23.38
28.81
27.96
29.75
29.36
30
22.38
26.83
25.82
27.77
27.07
10
27.42
32.11
32.91
32.72
33.15
Boat
20
25.08
29.04
29.47
29.66
29.88
30
23.83
27.18
27.63
27.83
27.92
10
24.57
30.54
31.17
30.51
31.45
Flinstones
20
21.72
27.17
27.40
27.28
28.12
30
19.95
25.02
25.29
25.35
26.04
10
30.07
33.18
34.27
34.28
34.75
House
20
27.09
29.98
30.90
30.89
31.76
30
25.49
27.96
28.96
28.94
29.88
10
26.09
30.92
32.07
31.63
32.87
Mit
20
22.78
27.08
27.88
27.77
28.62
30
20.95
24.91
25.60
25.60
26.22
10
26.87
31.47
32.35
31.81
33.08
Cameraman
20
24.05
27.94
28.51
28.27
29.13
30
22.54
26.00
26.48
26.38
27.05Denoised image using different denoising methods (σn=20, Lena).
(a)
Noisy image
(b)
Hard
(c)
LAWML
(d)
SURE-LET
(e)
P-Bivariate
(f)
Our methodDenoised image using different denoising methods (σn=20, Barbara).
(a)
Noisy image
(b)
Hard
(c)
LAWML
(d)
SURE-LET
(e)
P-Bivariate
(f)
Our methodDenoised image using different denoising methods (σn=20, Flintstones).
(a)
Noisy image
(b)
Hard
(c)
LAWML
(d)
SURE-LET
(e)
P-Bivariate
(f)
Our method.Denoised image using different denoising methods (σn=20, House).
(a)
Noisy image
(b)
Hard
(c)
LAWML
(d)
SURE-LET
(e)
P-Bivariate
(f)
Our methodFrom 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 Figures6–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.



---

*Source: 493976-2012-10-24.xml*