Julien Toulouse

Kamal Sharkas

Eric Brémond

Carlo Adamo

Rationale for a new class of double-hybrid approximations in density-functional theory

We provide a rationale for a new class of double-hybrid approximations introduced by Brémond and Adamo [J. Chem. Phys. 135, 024106 (2011)] which combine an exchange-correlation density functional with Hartree-Fock exchange weighted by λ and second-order Møller-Plesset (MP2) correlation weighted by λ 3 . We show that this double-hybrid model can be understood in the context of the density-scaled double-hybrid model proposed by Sharkas et al. [J. Chem. Phys. 134, 064113 (2011)], as approximating the density-scaled correlation functional Ec[n 1/λ ] by a linear function of λ, interpolating between MP2 at λ = 0 and a density-functional approximation at λ = 1. Numerical results obtained with the Perdew-Burke-Ernzerhof density functional confirms the relevance of this double-hybrid model.

The double-hybrid (DH) approximations introduced by Grimme [1], after some related earlier work [2,3], are increasingly popular for electronic-structure calculations within density-functional theory. They consist in mixing Hartree-Fock (HF) exchange with a semilocal exchange density functional and second-order Møller-Plesset (MP2) correlation with a semilocal correlation density functional:

E DH xc = a x E HF x + (1 -a x )E x [n] +(1 -a c )E c [n] + a c E MP2 c , (1)
where the first three terms are calculated in a usual selfconsistent hybrid Kohn-Sham (KS) calculation, and the last perturbative term is evaluated with the previously obtained orbitals and added a posteriori. The two empirical parameters a x and a c can be determined by fitting to a thermochemistry database. For example, the B2-PLYP double-hybrid approximation [1] is obtained by choosing the Becke 88 (B) exchange functional [4] for E x [n] and the Lee-Yang-Parr (LYP) correlation functional [5] for

E c [n],
and the empirical parameters a x = 0.53 and a c = 0.27 are optimized for the G2/97 subset of heats of formation.

Another approach has been proposed in which the perturbative contribution is evaluated with normal B3LYP orbitals rather than orbitals obtained with the weighted correlation density functional (1

-a c )E c [n] [6, 7].
Recently, Sharkas, Toulouse, and Savin [8] have provided a rigorous reformulation of the double-hybrid approximations based on the adiabatic-connection formalism, leading to the density-scaled one-parameter double- * Electronic address: julien.toulouse@upmc.fr † Electronic address: carlo-adamo@chimie-paristech.fr hybrid (DS1DH) approximation

E DS1DH,λ xc = λE HF x + (1 -λ)E x [n] +E c [n] -λ 2 E c [n 1/λ ] + λ 2 E MP2 c , (2)
where E c [n 1/λ ] is the usual correlation energy functional evaluated at the scaled density n 1/λ (r) = (1/λ) 3 n(r/λ). This reformulation shows that only one independent empirical parameter λ is needed instead of the two parameters a x and a c . The connection with the original doublehybrid approximations can be made by neglecting the density scaling

E c [n 1/λ ] ≈ E c [n], (3)
which leads to the one-parameter double-hybrid (1DH) approximation [8]

E 1DH,λ xc = λE HF x + (1 -λ)E x [n] +(1 -λ 2 )E c [n] + λ 2 E MP2 c . (4)
Equation ( 4) exactly corresponds to the double-hybrid approximation of Eq. ( 1) with parameters a x = λ and a c = λ 2 . Very recently, Brémond and Adamo [12] have proposed a new class of double-hybrid approximations where the correlation functional is weighted by (1 -λ 3 ) and the MP2 correlation energy is weighted by λ 3 , instead of (1 -λ 2 ) and λ 2 , respectively. Applying this formula with the Perdew-Burke-Ernzerhof (PBE) [13] exchangecorrelation density functional, they have constructed the PBE0-DH double-hybrid approximation which performs reasonably well. In this work, we give a rationale for this class of double-hybrid approximations. For this, we start by recalling that the density-scaled correlation functional E c [n 1/λ ] tends to the second-order Görling-Levy (GL2) [14] correlation energy when the density is squeezed up to the high-density limit (or weakinteraction limit) 3)], PBE with density scaling (from the parametrizations of Ref. 9), and linear interpolation between MP2 (with PBE orbitals) and PBE [Eq. ( 9)]. The MP2 calculations for He and Be (including core excitations) have been performed with the cc-pV5Z and cc-pCV5Z basis sets [10,11], respectively. which is finite for nondegenerate KS systems. The GL2 correlation energy can be decomposed as (see, e.g., Ref. [START_REF] Engel | A Primer in Density Functional Theory[END_REF])

lim λ→0 E c [n 1/λ ] = E GL2 c , (5)
E GL2 c = E MP2 c + E ∆HF c , (6)
where

E MP2 c
is the usual MP2 correlation energy expression

E MP2 c = - 1 4 ij ab | φ i φ j ||φ a φ b | 2 ε a + ε b -ε i -ε j , (7)
with the antisymmetrized two-electron integrals φ i φ j ||φ a φ b , and E ∆HF c is an additional contribution involving the difference between the local multiplicative KS exchange potential vKS x and the nonlocal nonmultiplicative HF exchange potential vHF

x E ∆HF c = - i a φ i |v KS x -vHF x |φ a 2 ε a -ε i . (8)
In both Eqs. ( 7) and ( 8), φ k are the KS orbitals and ε k are their associated energies, and the indices i,j and a,b stand for occupied and virtual orbitals, respectively. The single-excitation contribution E ∆HF c vanishes for two-electron systems, and in most other cases is negligible [START_REF] Engel | A Primer in Density Functional Theory[END_REF], so that the GL2 correlation energy is well approximated by just the MP2 contribution (evaluated with KS orbitals),

E GL2 c ≈ E MP2 c
. This leads us to propose an approximation for E c [n 1/λ ] based on a linear interpolation formula

E c [n 1/λ ] ≈ (1 -λ)E MP2 c + λE c [n]. (9)
Plugging Eq. ( 9) into Eq. ( 2), we directly arrive at what we call the linearly scaled one-parameter double-hybrid (LS1DH) approximation

E LS1DH,λ xc = λE HF x + (1 -λ)E x [n] +(1 -λ 3 )E c [n] + λ 3 E MP2 c , (10)
with the weights (1 -λ 3 ) and λ 3 , thus giving a stronger rationale to the expression that Brémond and Adamo have proposed on the basis of different considerations.

Further insight into this approximation can be gained by rewriting Eq. ( 9) in the alternative form

E c [n 1/λ ] ≈ E MP2 c + λ E c [n] -E MP2 c , (11)
which can then be interpreted as a first-order expansion in λ around λ = 0 with E c [n] -E MP2 c approximating the third-order correlation energy correction E 2) just gives the usual one-parameter hybrid (1H) approximation with the full correlation density functional [16,[START_REF] Ernzerhof | Density Functional Theory[END_REF]

E 1H,λ xc = λE HF x + (1 -λ)E x [n] + E c [n]. ( 12
)
In this sense, the LS1DH approximation of Eq. ( 10) can be considered as a next-order approximation in λ beyond the usual hybrid approximation. Figure 1 illustrates the different approximations to the density-scaled correlation energy E c [n 1/λ ] as a function of λ for the He and Be atoms. The accurate reference curve is from the parametrization of Ref. 9. For λ = 1, it reduces to the exact correlation energy, while for λ = 0 it is the GL2 correlation energy which tends to overestimate the correlation energy. In between these two limits, it is nearly linear with λ. The PBE correlation energy without density scaling [Eq. ( 3)] is the crudest approximation to E c [n 1/λ ]. The PBE correlation energy with density scaling (taken from the parametrization of Ref. 9) gives a nearly linear curve. It is a fairly good approximation for He, but a less good approximation for Be where it underestimates the correlation energy for all λ, especially at λ = 0. This is due to the presence of static correlation in this system that is not described by the PBE functional. Finally, the linear interpolation of Eq. ( 9) between MP2 (evaluated with PBE orbitals) and PBE appears as a good approximation for both He and Be. In fact, the linear interpolation is clearly the best approximation for Be, at least with the PBE functional.

For a more comprehensive comparison of the different approximations, we have performed calculations on the AE6 and test sets [START_REF] Lynch | [END_REF] with the 1DH, DS1DH, and LS1DH double hybrids for the PBE functional, using a development version of the MOLPRO 2010 program [START_REF] Werner | Molpro, version 2010.1, a package of ab initio programs[END_REF]. The AE6 set is a small representative benchmark set of six atomization energies consisting of SiH → HS + H 2 . All the calculations are performed at the optimized QCISD/MG3 geometries [20]. We use the Dunning cc-pVQZ basis set [START_REF] Dunning | [END_REF]22]. Core electrons are kept frozen in all our MP2 calculations. Spin-restricted calculations are performed for all the closed-shell systems, and spin-unrestricted calculations for all the open-shell systems. In Fig. 2, we plot the mean absolute errors (MAE) for the two sets as a function of λ. For the AE6 set, the MAEs of the DS1DH and LS1DH approximations are quite similar, and are both much smaller that the MAE of the 1DH approximation for a wide range of λ. For the BH6 set, the LS1DH approximation gives a MAE which is significantly smaller than those of both DS1DH and 1DH for intermediate values of λ.

As regards the choice of the parameter λ, the present data on the AH6 set gives an optimal value of λ = 0.75 for the LS1DH double hybrid, with a minimal MAE of 3.59 kcal/mol, and for the BH6 set the optimal value is λ = 0.70 giving a minimal MAE of 0.73 kcal/mol. However, the MAE is not very sensitive to the value of λ around the optimal value, and Brémond and Adamo [12] have argued for using λ = 0.5 in defining the PBE0-DH approximation, using a similar argument as the one used by Becke for his "half-and-half" hybrid [23].

In summary, we have shown that the new class of double-hybrid approximations named here LS1DH [Eq. (10)] can be understood as approximating the density-scaled correlation functional E c [n 1/λ ] by a linear function of λ, interpolating between MP2 at λ = 0 and a density-functional approximation at λ = 1. Numerical results obtained with the PBE density functional confirms that the LS1DH approximation is a relevant double-hybrid model, and in fact tends to be more accurate than the DS1DH double-hybrid model [Eq. ( 2)] in which density scaling is applied to the PBE correlation functional. More generally, it can be expected that the LS1DH double-hybrid model will be more accurate than the DS1DH double-hybrid model when applied with a density-functional approximation that is inaccurate in the high-density limit. We hope that this work will help constructing other theoretically justified double-hybrid approximations.

FIG. 1 :

 1 FIG. 1: (Color online) Density-scaled correlation energy Ec[n 1/λ ] for the He (left) and Be (right) atoms as a function of λ. Accurate calculations (from the parametrizations of Ref. 9) are compared with different approximations: PBE without density scaling [Eq. (3)], PBE with density scaling (from the parametrizations of Ref. 9), and linear interpolation between MP2 (with PBE orbitals) and PBE [Eq. (9)]. The MP2 calculations for He and Be (including core excitations) have been performed with the cc-pV5Z and cc-pCV5Z basis sets[10,11], respectively.

FIG. 2 :

 2 FIG.2:(Color online) Mean absolute errors for the AE6 (left) and BH6 (right) test sets as functions of the parameter λ for the 1DH, and DS1DH, and LS1DH approximations with the PBE exchange-correlation density functional. All calculations were carried out with the cc-pVQZ basis set.

  perturbation theory. In comparison, the zeroth-order approximation E c [n 1/λ ] ≈ E MP2 c plugged in Eq. (

  4 , S 2 , SiO, C 3 H 4 (propyne), C 2 H 2 O 2 (glyoxal), and C 4 H 8 (cyclobutane). The BH6 set is a small representative benchmark set of forward and reverse barrier heights of three reactions, OH + CH 4 → CH 3 + H 2 O, H + OH → O + H 2 , and H + H 2 S

Acknowledgments

We thank A. Savin for discussions. This work was partly supported by Agence Nationale de la Recherche (ANR) via contract numbers 07-BLAN-0272 (Wademecom) and 10-BLANC-0425 (DinfDFT). Sanofi-Aventis is also acknowledged for financial support.


R Soulard

P Carré

COLOR EXTENSION OF MONOGENIC WAVELETS WITH GEOMETRIC ALGEBRA : APPLICATION TO COLOR IMAGE DENOISING

Keywords: Color Wavelets, Analytic, Monogenic, Wavelet transforms, Image analysis, Denoising

We define a color monogenic wavelet transform. This is based on the recent grayscale monogenic wavelet transform and an extension to color signals aimed at defining non-marginal tools. Wavelet based color image processing schemes have mostly been made by using a grayscale tool separately on color channels. This may have some unexpected effect on colors because those marginal schemes are not necessarily justified. Here we propose a definition that considers a color (vector) image right at the beginning of the mathematical definition so we can expect to bring an actual color wavelet transform -which has not been done so far to our knowledge. This so provides a promising multiresolution color geometric analysis of images. We show an application of this transform with a statistical modeling of coefficients for color denoising issue.

INTRODUCTION

Wavelets have been widely used for handling images for more than 20 years. It seems that the human visual system sees images through different channels related to particular frequency bands and directions; and wavelets provide such decompositions. Since 2001, the analytic signal and its 2D generalizations have brought a great improvement to wavelets [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF][START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF][START_REF] Unser | Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform[END_REF]] by a natural embedding of an AM/FM analysis in the subband coding framework. This yields an efficient representation of geometric structures in grayscale images thanks to a local phase carrying geometric information complementary to an amplitude envelope having good invariance properties. So it codes the signal in a more coherent way than standard wavelets. The last and seemingly most appropriate proposition [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform[END_REF] of analytic wavelets for image analysis is based on the monogenic signal [START_REF] Felsberg | Low-level image processing with the structure multivector[END_REF] defined with geometric algebra.

In parallel a color monogenic signal was proposed [START_REF] Demarcq | The color monogenic signal: a new framework for color image processing[END_REF] as a mathematical extension of the monogenic signal; paving the way to non-marginal color tools especially by using geometric algebra and above all by considering a color signal right at the foundation of the mathematical construction.

We define here a color monogenic wavelet transform that extends the monogenic wavelets of [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform[END_REF] to color. These analytic wavelets are defined for color 2D signals (images) and avoid the classical pitfall of marginal processing (grayscale tool used separately on color channels) by relying on a sound mathematical definition. We may so expect to handle coherent information of multiresolution color geometric structure; which would make easier any wavelet based color image processing. To our knowledge color wavelets have not been proposed so far.

We first give a technical study of analytic signal/wavelets with the intent to popularize them since they rely on non-trivial concepts of geometric algebra, complex/harmonic analysis, as well as non-separable wavelet frames. Then we describe our color monogenic wavelet transform; and finally an application to color denoising will be presented.

Notations : 2-vector coordinates : x = (x, y) , ω = (ω 1 , ω 2 ) ∈ R 2 ; k ∈ Z 2 Euclidean norm : x = x 2 + y 2
Complex imaginary number : j ∈ C Argument of a complex number : arg Convolution symbol : * Fourier transform : F

ANALYTIC SIGNAL AND 2D GENERALIZATION

An analytic signal s A is a multi-component signal associated to a real signal s to analyze. The definition is well known in the 1D case where s A (t) = s(t) + j (h * s)(t) is the complex signal made of s and its Hilbert transform (with h(t) = 1 πt ). The polar form of the 1D analytic signal provides an AM/FM representation of s with |s A | being the amplitude envelope and ϕ = arg (s A ) the instantaneous phase. This classical tool can be found in many signal processing books and is used in communications for example.

Interestingly we can also interpret the phase in terms of signal shape i.e. there is a direct link between the angle ϕ and the local structure of s. Such a link between a 2D phase and local geometric structures of images would be very attractive in image processing. That is why there were several attempts to generalize it for 2D signals; and among them the monogenic signal [START_REF] Felsberg | Low-level image processing with the structure multivector[END_REF] seems the most advanced since it is rotation invariant.

The Monogenic Signal

Without going beyond strictly needed details we here review the key points of the fundamental construction of the monogenic signal; which will be necessary to understand the color extension.

The definition of the 1D case given above can be interpreted in terms of signal processing : the Hilbert transform makes a "pure π 2 -dephasing". But such a dephasing is not straightforward to define in 2D (same issue with many 1D signal tools) so let us look at the equivalent complex analysis definition of the 1D analytic signal. It says that s A is the holomorphic extension of s restricted to the real line. But complex algebra is impeding for generalizations to higher dimensions. To bypass this limitation we can see a holomorphic function like a 2D harmonic field that is an equivalent harmonic analysis concept involving the 2D Laplace equation ∆f = 0. It so can be generalized within the framework of 3D harmonic fields by using the 3D Laplace operator ∆ 3 = δ δx + δ δy + δ δz . The whole generalization relies on this natural choice and remaining points are analogous to the 1D case (see [START_REF] Felsberg | Low-level image processing with the structure multivector[END_REF] for more details). Note that in Felsberg's thesis this construction is expressed in terms of geometric algebra but here we avoided it for simplicity's sake. Finally the 2D monogenic signal s A associated to s is the 3-vector valued signal :

s A (x) =   s(x) s r1 (x) = x 2π x 3 * s(x) s r2 (x) = y 2π x 3 * s(x)   (1)
Where s r1 and s r2 are analogous to the imaginary part of the complex 1D analytic signal. Interestingly, this construction reveals the two components of a Riesz transform :

R{s} = ( s r1 (x) , s r2 (x) ) = x 2π x 3 * s(x) , y 2π x 3 * s(x) (2)
in the same way that the 1D case exhibits a Hilbert transform. Note that we get back to a signal processing interpretation since the Riesz transform can also be viewed like a pure 2D dephasing.

In the end, by focusing on the complex analysis definition of the analytic signal we end up with a convincing generalization of the Hilbert transform. Now recall that the motivation to build 2D analytic signals arises from the strong link existing between the phase and the geometric structure. To define the 2D phase related to the Riesz transform the actual monogenic signal must be expressed in spherical coordinates that yield the following amplitude envelope and 2-angle phase :

Amplitude : A = s 2 +s 2 r1 +s 2 r2 s = A cos ϕ Orientation : θ = arg (s r1 + j s r2 ) s r1 = A sin ϕ cos θ 1D Phase : ϕ = arccos ( s A ) s r2 = A sin ϕ sin θ (3)
Felsberg shows a direct link between the angles θ and ϕ and the geometric local structure of s.

The signal is so expressed like an "A-strong" 1D structure with orientation θ. ϕ is analogous to the 1D local phase and indicates if the structure is rather a line or an edge. A direct drawback is that intrinsically 2D structures are not handled. Yet this tool found many applications in image analysis from contour detection to motion estimation (see [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform[END_REF] and references therein p. 1). From a signal processing viewpoint the AM/FM representation provided by an analytic signal is accordingly well suited for narrowband signals. That is why it seems natural to embed it in a wavelet transform that performs subband decomposition. We now present the monogenic wavelet analysis proposed in [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform[END_REF].

MONOGENIC WAVELETS

So far there is one proposition of computable monogenic wavelets in the literature [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform[END_REF]. It provides 3-vector valued monogenic subbands consisting of a rotation-covariant magnitude and this new 2D phase. This representation -specially defined for 2D signals -is a great theoretic improvement of the complex and quaternion wavelets [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF][START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF]; as well as the monogenic signal itself is an improvement of its complex and quaternion counterparts.

The proposition of [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform[END_REF] consists of one real-valued "primary" wavelet transform in parallel with an associated complex-valued wavelet transform. Both transforms are linked each other by the Riesz transform so they carry out a multiresolution monogenic analysis. We end up with 3-vector coefficients forming subbands that are monogenic.

Primary transform

The primary transform is real-valued and relies on a dyadic pyramid decomposition tied to a wavelet frame. Only one 2D wavelet is needed and the dyadic downsampling is done only at the low frequency branch; leading to a redundancy of 4:3. The scaling function ϕ γ and mother wavelet ψ are defined in the Fourier domain :

ϕ γ F ←→ 4(sin 2 ω 1 2 + sin 2 ω 2 2 ) -8 3 sin 2 ω 1 2 sin 2 ω 2 2 γ 2 ω γ (4)
ψ(x) = (-∆) γ 2 ϕ 2γ (2x) (5)
Note that ϕ γ is a cardinal polyharmonic spline of order γ and spans the space of those splines with its integer shifts. It also generates -as a scaling function -a valid multiresolution analysis. This particular construction is made by an extension of a wavelet basis (non-redundant) related to a critically-sampled filterbank. This extension to a wavelet frame (redundant) adds some degrees of freedom used by the authors to tune the involved functions. In addition a specific subband regression algorithm is used at the synthesis side. The construction is fully described in [START_REF] Unser | The pairing of a wavelet basis with a mildly redundant analysis via subband regression[END_REF].

The monogenic transform

The second "Riesz part" transform is a complex-valued extension of the primary one. We define the associated complex-valued wavelet by including the Riesz components :

ψ = - x 2π x 3 * ψ(x) + j y 2π x 3 * ψ(x) (6)
It can be shown that it generates a valid wavelet basis and that it can be extended to the pyramid described above. The joint consideration of both transforms form monogenic subbands from which can be extracted the amplitude and phase for an overall redundancy of 4:1. So far no applications of the monogenic wavelets have been proposed. In [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform[END_REF] a demonstration of AM/FM analysis is done with fine orientation estimation and gives very good results in terms of coherency and accuracy. Accordingly this tool may be rather used for analysis tasks than processing.

Motivated by the powerful analysis provided by the monogenic wavelet transform we propose now to extend it for color images. Figure 2: Color MWT scheme. Each color channel is analyzed with the primary wavelet transform symbolized by a ψ bloc and the sum "R + G + B" is analyzed with the "Riesz part" wavelet transform (ψ r1 and ψ r2 blocs).

ψ R ψ G ψ B ψ r1 ψ r2

COLOR MONOGENIC WAVELETS

We define here our proposition that combines a fundamental generalization of the monogenic signal to color with the monogenic wavelets described above. The challenge is to avoid the classical marginal definition that would be applying a grayscale monogenic transform on each of the three color channels of a color image. We believe that the monogenic signal has a favorable theoretical framework for a color extension and this is why we propose to start from this particular wavelet transform rather than from a more classical one.

The color generalization of the monogenic signal is expressed within the geometric algebra framework. This algebra is very general and embeds the complex and quaternion as subalgebras. Its elements are "multivectors" naturally linked with various geometric entities. The use of this fundamental tool is gaining popularity in the literature because it allows rewriting sophisticated concepts with simpler algebraic expressions and so paves the way to innovative ideas and generalizations in many fields.

For simplicity's sake and since anyway we would not have enough space to present the fundamentals of geometric algebra we here express the construction in classical terms; as we already did section 2. Yet we may sometimes point out some necessary specific mechanisms but we refer the reader to [START_REF] Felsberg | Low-level image processing with the structure multivector[END_REF][START_REF] Demarcq | The color monogenic signal: a new framework for color image processing[END_REF] for further details.

The Color Monogenic Signal

Starting from Felsberg's approach that is originally expressed in the geometric algebra of R 3 ; the extension proposed in [START_REF] Demarcq | The color monogenic signal: a new framework for color image processing[END_REF] is written in the geometric algebra of R 5 for 3-vector valued 2D signals of the form (s R , s G , s B ). By simply increasing the dimensions we can embed each color channel along a different axis and the original equation from Felsberg involving a 3D Laplace operator can be generalized in 5D with ∆ 5 = δ δx 1 + δ δx 2 + δ δx 3 + δ δx 4 + δ δx 5 . Then the system can be simplified by splitting it into three systems with a 3D Laplace equa-tion, reducing to applying Felsberg's condition to each color channel. At this stage appears the importance of geometric algebra since an algebraic simplification between vectors leads to a 5-vector color monogenic signal that is non-marginal. Instead of naively applying the Riesz transform to each color channel, this fundamental generalization carries out the following color monogenic signal : s A = (s R , s G , s B , s r1 , s r2 ) where s r1 and s r2 are the Riesz transform applied to s R +s G +s B . Now the color extension of Felsberg's monogenic signal is defined let us construct the color extension of the monogenic wavelets.

The Color Monogenic Wavelet Transform

We can now define a wavelet transform whose subbands are color monogenic signals. The goal is to obtain vector coefficients of the form

(c R , c G , c B , c r1 , c r2 ) such that c r1 = x 2π x 3 * (c R + c G + c B ) and c r2 = y 2π x 3 * (c R + c G + c B )
. It turns out that we can very simply use the transforms presented above by applying the primary one on each color channel and the Riesz part on the sum of the three. The five related color wavelets illustrated Fig. 1 and forming one color monogenic wavelet ψ A are :

ψ R =   ψ 0 0   ψ G =   0 ψ 0   ψ B =   0 0 ψ   (7)
ψ r1 =    x 2π x 3 * ψ x 2π x 3 * ψ x 2π x 3 * ψ    ψ r2 =    y 2π x 3 * ψ y 2π x 3 * ψ y 2π x 3 * ψ    (8)
ψ A = (ψ R , ψ G , ψ B , ψ r1 , ψ r2 ) (9)
We then get 5-vector coefficients verifying our conditions and so forming a color monogenic wavelet transform. The associated decomposition is described by the diagram Fig. 2. This provides a multiresolution color monogenic analysis made of a 5-vector valued pyramid transform. The 5 decompositions of two images are shown Fig. 3 from left to right. Each one consists of 4 juxtaposed image-like subbands resulting from 3-level decomposition. We fixed γ = 3 that gave good experimental results.

Let us look at the first 3 graymaps. These are the 3 primary transforms c R , c G and c B where white (resp. black) pixels are high positive (resp. negative) values. Note that our transform is non-separable and so provides at each scale only one subband related to all orientations. We are not subjected to the arbitrarily separated horizontal, vertical and diagonal analyses of usual wavelets. This advantage is even greater in color. Whereas marginal separable transforms show 3 arbitrary orientations within each color channel -which is not easily interpretable -the color monogenic wavelet transform provides a more compact energy representation of the color image content regardless of the local orientation. The color information is well separated through c R , c G and c B : see that blue contours of first image are present only in c B . And in each of the 3 decompositions it is clear that every orientation is equally represented all along the round contours. That is different from separable transform that privileges particular directions. The multiresolution framework makes the horizontal blue low frequency structure of second image be coded mainly in the third scale of c B .

But the directional analysis is not lost thanks to the Riesz part that completes this representation. Now look at the "2-in-1" last decomposition forming the Riesz part. It is displayed in one color map where the geometric energy c 2 r1 + c 2 r2 is encoded into the intensity (with respect to the well known HSV color space) and the orientation arg (c r1 +jc r2 ) (π) is encoded in the hue (e.g. red is for {0, π} and cyan is for ± π

2 ). This way of displaying the Riesz part well reveals the provided geometric analysis of the image.

The Riesz part makes a precise analysis that is local both in space and scale. If there is a local color geometric structure in the image at a certain scale the Riesz part exhibits a high intensity in the corresponding position and subband. This is completed with an orientation analysis (hue) of the underlying structure. For instance a horizontal (resp. vertical) structure in the image will be coded by a cyan (resp. red) intense point in the corresponding subband. The orientation analysis is strikingly coherent and accurate. See for example that color structures with constant orientation (second image) exhibit a constant hue in the Riesz part over the whole structure.

Note that low intensity corresponds to "no structure" i.e. where the image has no geometric information. It is coherent not to display the orientation (low intensity makes the hue invisible) for these coefficients since this data has no sense in those cases.

In short the color and geometric information of the image are well separated from each other and the orientation analysis is very accurate. In addition the invariance properties of the primary and Riesz wavelet transforms are kept in the color extension for a slight overall redundancy of 20 : 9 ≈ 2.2. This transform is non-marginal because RGB components is considered as well as intensity (R+ G+ B) -which involves two different color spaces.

WAVELET COEFFICIENTS STUDY FOR DENOISING ISSUE

Let us now study reduction of additive gaussian noise by wavelet domain thresholding. First we will experimentally characterize decomposition of a sole noise to identify effect of scale on distribution and subband correlation. Then we will apply color MWT thresholding to achieve color denoising and show experimental results.

Characterization of noise

With classical biorthogonal wavelets optimal thresholding involves a constant threshold over the whole transform (see citeDonoho1995). Here we have a non-orthogonal transform so we need to study how the noise term of a noisy image is transformed through color monogenic wavelet analysis.

As shown Fig. 4 we observe that decomposition of a centered gaussian noise with variance σ 2 remains centered and gaussian after decomposition -with different variances though. For a given scale, distribution of Primary part as well as real and imaginary components of the Riesz part is clearly gaussian. Experimental values for standard deviations are given in table 1. These values can also be derived analytically from definition of filters. Recall that linear filtering of a stationary random signal x by filter H with output y can be studied with power spectral densities

Ψ x = |F[x]| 2 (PSD) and autocorrelations R x (τ ) = F -1 Ψ x . In par- ticular we have Ψ y = Ψ x |H| 2 . The output variance σ 2 y is equal to R y (0) which reduces to σ 2 y = σ 2 x 4π 2 (2π,2π) (0,0)
|H(ω)| 2 dω. For example the first scale output of our filterbank is directly linked to the first stage highpass filter :

σ 2 1 = σ 2 4π 2 (2π,2π) (0,0) 4(sin 2 ω 1 2 + sin 2 ω 2 2 ) -8 3 sin 2 ω 1 2 sin 2 ω 2 2 γ 2 ω γ 2 dω (10)
Remaining coefficients are tied to equivalent filters of each filterbank output.

A slight spatial correlation is introduced by the filterbank due to its redundancy. Whereas input noise autocorrelation R = σ 2 δ is null for τ > 0; R s has some small coefficients around Noisy Color MWT thresholding τ = 0 -though with very fast decay -in a neighborhood of about 2 pixels (R s is of the same form as wavelet functions showed Fig. 1). Some correlation is also present between primary and Riesz parts of subbands. So this transform has to be used carefully if some data compression is wanted.

Let us consider modulus of "Riesz part" outputs that will be of interest in the sequel. As real and imaginay part of subbands follow centered gaussian distributions with same variance σ 2 s , the complex subband follows a 2D isotropic gaussian distribution with variance 2σ 2 s . As a result the modulus of "Riesz part" outputs follow a Rayleigh distribution f s (m) = m σ 2 s e -m 2 /2σ 2 s 1 m>0 . We can see that modeled curves Fig. 4 well correspond to measured histograms. Now that subbands are statistically modeled we can apply automatic thresholding to perform color denoising.

Thresholding

We apply a soft thresholding of color wavelet coefficients with a different threshold for each scale s. For the primary part threshold is fixed to T prim s = k * σ s , with σ s as reported in table 1 and k = 3. This choice of k is classical with usual wavelet denoising.

For the Riesz part that is complex valued we propose to threshold modulus because this data is tied to some amplitude information related to geometrical structures. But modulus follows a Rayleigh distribution so classical thresholding does not hold. A similar issue is dealt in [START_REF] Kovesi | Phase preserving denoising of images[END_REF] with modulus of gabor wavelets. Knowing that the mean (resp. variance) of a Rayleigh distribution is µ r = σ π/2 (resp. σ 2 r = σ 2 (4 -π)/2), the noise shrinkage threshold at each scale is set to be some number of standard deviations beyond the mean of the distribution -that is µ r +kσ r . Since in our case a Rayleigh distribution is also involved this definition is pertinent for our Riesz part thresholding. So we have T Riesz s = σ s ( π 2 + k 4-π 2 ) (as in [START_REF] Kovesi | Phase preserving denoising of images[END_REF] we will take k = 2). Threshold is classically estimated by a simple analysis of first scale (High frequencies) (see [START_REF] Mallat | A wavelet tour of signal processing[END_REF]). By assuming that first scale contains mainly noise plus a few significant coefficients, one can use the experimental median of coefficients to well estimate standard deviation of the gaussian noise distribution. So we have the estimate σ = median(|W 1 |)/0.6745 with W 1 being the coefficients of first scale. This estimation is done in the primary part but can be used for the Riesz part as both T prim s and T Riesz s only depends on estimate of σ s . By using reference coefficients of Table 1 we use σs = σs σ 1 σ. A result of color denoising is showed Fig. 5. This is a 4 scale thresholding of color MWT coefficients where T P rim s is used in Primary subbands and T Riesz s in the Riesz part. Thresholds are processed from estimation of noise variance described above; and low frequency subband is not modified. This first result is comparable to a classical wavelet denoising and so confirms the potential of this transform. In practice the Riesz part remains difficult to deal with and many improvements may be added.

CONCLUSION

We define a color extension of the recent monogenic wavelet transform proposed in [START_REF] Unser | Multiresolution monogenic signal analysis using the riesz-laplace wavelet transform[END_REF]. This extension is non-marginal since it takes care of considering a vector signal at the very beginning of the fundamental construction and leads to a definition basically different from the marginal approach. The use of non-separable wavelets joint with the monogenic framework allows for a good orientation analysis well separated from the color information. This color transform can be a great color image analysis tool thanks to this good separation of information through various data. A statistical modeling of coefficient for thresholding/denoising issue is given.

Although it is not marginal the color generalization has a marginal style since it reduces to apply the Riesz transform on the intensity of the image. So the geometric analysis is done without considering the color information and it would be much more attractive to have a complete representation of the color monogenic signal into magnitude and phase(s) with color/geometric interpretation.

Figure 1 :

 1 Figure 1: Space representation of the 5 color wavelets.

Figure 3 :

 3 Figure 3: Color MWT of images. The two components of the Riesz part are displayed in the same graphic with the magnitude of c r1 +jc r2 encoded in the intensity and argument (local orientation) encoded in the hue.

Figure 4 :

 4 Figure 4: Histograms (bars) of primary subbands and modulus of Riesz subbands; and PDF models (line).

Figure 5 :

 5 Figure 5: Denoising.

Table 1 :

 1 Standard deviations of subbands after decomposition of gaussian noise with variance 1 (with γ = 3).

	Scale s	1(High freq.)	2	3	4
	St. dev. σ s	1.339	1.502 1.512 1.560


Raphaël Soulard
email: raphael.soulard@univ-poitiers.fr

Philippe Carré
email: carre@sic.univ-poitiers.fr

Quaternionic Wavelets for Texture Classification

Keywords: Wavelet transform, 2D Phase, Quaternionic Wavelet Transform, Image texture analysis, Image classification

This article proposes a study of the recent quaternionic wavelet transform (QWT) from a practical point of view through a digital image analysis application. Based on a theoretic 2D generalization of the analytic signal leading to a strong 2D signal modeling, this representation uses actual 2D analytic wavelets and yields subbands having a shift-invariant magnitude and a 3angle phase, using the quaternion algebra.

Our experiment furthers the understanding of this quite sophisticated tool, and shows its practical interest by a clear improvement of a famous wavelet application : texture classification. Thanks to coherent multiscale analysis brought by the QWT we obtain better classification results than with standard wavelets in a similar process.

Introduction

Since the 90's wavelets have been widely used for the analysis of textural images [START_REF] Tuceyran | Handbook of pattern recognition and computer vision[END_REF][START_REF] Lasmar | Multivariate statistical modeling for texture analysis using wavelet transforms[END_REF]. Accordingly the human visual system sees textures through different channels related to particular frequency bands and directions; and wavelet representations offer this kind of decomposition. So by extracting simple features from each subband, the image content can be characterized separately at different scales and orientations, making an efficient texture analysis.

In 2001, the well known power of the analytic signal for signal modeling led to a new definition of wavelets : the Complex Wavelet Transform (CWT) [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF], whose subbands are analytic signals giving access to a shift-invariant amplitude envelope (magnitude). The CWT codes signals in a more coherent way than standard wavelets (DWT); which overcomes their famous shift-variance problem. Its superiority over DWT for texture analysis has been shown in [START_REF] De Rivaz | Complex wavelet features for fast texture image retrieval[END_REF][START_REF] Celik | Multiscale texture classification using dual-tree complex wavelet transform[END_REF]. But CWT phase is ambiguous for 2D signals and is almost not used in applications excepted in a very recent work [START_REF] Celik | Bayesian texture classification and retrieval based on multiscale feature vector[END_REF].

In parallel, the work in [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF] provides a strong theoretic 2D generalization of the analytic signal defined in quaternion algebra. The 3-angle local 2D phase related to it contains 2D geometric information. The Quaternionic Wavelet Transform (QWT) -first proposed in 2004 [START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF][START_REF] Bayro-Corrochano | The theory and use of the quaternion wavelet transform[END_REF]] -uses wavelets that are 2D-analytic in the sense of Bülow, providing quaternion valued subbands with shift-invariant magnitude and this new 2D phase. This representation -specially defined for 2D signals -is a great theoretic improvement of CWT yielding a coherent description of local 2D structures.

Currently, many evolutions of wavelets focus on directionality (Contourlets, Complex Wavelets, Directionlets, Adaptive lifting schemes); bringing some redundancy to fit numerous directions from a quite algorithmic point of view. QWT is a different approach to a coherent representation of local structures in images based on a strong 2D signal processing theory; and sets some redundancy in a local phase rather than in directionality.

Given the promising theoretical properties of this new transform, we propose here a global study of this sophisticated tool to further its understanding through a famous application of wavelets. We aim to show its potential and verify its practical superiority over standard DWT in a comparative texture analysis application.

The main point is the understanding of the information carried by QWT coefficients. QWT magnitude can intuitively be used like DWT but interpretation of phase is far from straightforward. With a first QWT based feature extraction algorithm this work gives an application not did yet to our knowledge and furthers practical use of QWT coefficients. We expect QWT to provide better texture analysis than DWT -with simple feature extraction methods -thanks to better separation of information and shift-invariance.

An additional experimental comparison is done with CWT because the topic is near.

After a presentation of the transform, we propose several feature extraction methods based on both DWT and QWT. To compare with literature, we also describe some CWT based feature descriptors recently published in [Celik andTjahjadi, 2009, 2011]. The use of QWT phase is then discussed as well as the combination of QWT magnitude and phase based features. We finally give experimental classification results to compare DWT and QWT in terms of recognition rates. Note that the point is to observe changes in recognition and not to get a competitive texture recognition performance. The aim is to carry out a simple process to help us studying wavelet transforms, rather than achieving excellent recognition. Types of texture better described by one or another transform will be analyzed in terms of visual/geometric content.

Quaternionic Wavelet Transform

The Quaternionic Wavelet Transform (QWT) is a perfect reconstruction 2D filterbank decomposition for grayscale images; addressing three common drawbacks of standard DWT :

• Oscillations -Using oscillating wavelets complicates the representation of a simple structure, involving several coefficients in a neighborhood to compensate each other their own oscillations, and so making local analysis difficult.

• Shift-variance -The look of a wavelet transform may substantially change by a small shift of the image; causing incoherence while most recognition processes need invariant features.

• No phase notion -Fourier phase and local phase from the analytic signal are powerful analysis tools so it is interesting to fit DWT up with a phase notion in a pattern recognition context.

The key is to embed the analytic signal modeling into the wavelet framework. The Dual-Tree CWT [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF] achieves 1D analytic wavelet analysis with a 2 times redundant perfect reconstruction filterbank overcoming the 3 drawbacks listed above. But its 2D version yields an ambiguous phase failing to efficiently describe local structures. More precisely, the phase angle encodes a local shift around pixel position -in the exact direction that is orthogonal to the intrinsic orientation of considered subband.

Since the structure may be not well aligned with wavelet orientation, this shift information is not sufficient to accurately localize the structure.

The point then is the 2D generalization of analytic signal and the geometric meaning of the 2D phase extracted from it. This is why the QWT uses the 2D analytic signal defined in [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF] within the quaternion algebra. According to Bülow, this generalization is superior to previous attempts in C -including the one implicitly related to 2D CWT -and is well adapted for 2D signals whereas standard complex analytic signal is only adapted in 1D. The extracted quaternionic local phase (3 angles) contains rich geometric information. Note that Bülow proposed quaternionic Gabor filters in 1999 [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF] and used the phase for image analysis applications.

The QWT provides a quaternionic multiresolution analysis with 2D analytic wavelets overcoming the above-cited drawbacks of DWT. Whereas DWT coefficients are real, QWT is quaternion valued i.e. 4-vectors made of one magnitude and a 3-angle phase. Information is better separated to describe more explicitly the image content. In 2004, the Rice University from Houston proposes to use their Dual-Tree algorithm to carry out a QWT with perfect reconstruction filterbanks [START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF] (that we use in this work). Note that at the same time, Bayro proposes a quaternionic Gabor pyramid [START_REF] Bayro-Corrochano | The theory and use of the quaternion wavelet transform[END_REF], windowing the quaternionic Fourier transform (QFT), but this approach has no perfect reconstruction so the use can be limited for example in image coding.

In this section we present the quaternionic 2D analytic signal upon which the QWT is based, the definition of the transform, the algorithm we used and the use of this tool.

The Quaternionic Analytic Signal

In his thesis [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF], T. Bülow showed that complex algebra C is only adapted for handling 1D signals and that 2D signals are best described by embedding signal processing tools in the more general quaternion algebra H. He defined a 'Quaternionic 2D Fourier Transform' (QFT) and a 'Quaternionic 2D Analytic Signal'.

A quaternion is a generalization of a complex number, related to 3 imaginary units i, j, k following the rules i 2 = j 2 = k 2 = -1 and ij = -ji = k; and written q = a + bi + cj + dk. A polar form q = |q|e iϕ e jθ e kψ analogous to the complex one is defined by one modulus |q| and three angles (ϕ, θ, ψ) called phase. The phase corresponds here to Euler angles in a 4D Euclidean space :

ϕ = atan( 2(cd+ab) , a 2 -b 2 +c 2 -d 2 )/2 + kπ ∈ [-π; π] θ = atan( 2(bd+ac) , a 2 +b 2 -c 2 -d 2 )/2 ∈ [-π 2 ; π 2 ] ψ = arcsin( 2(ad-bc) )/2 ∈ [-π 4 ; π 4 ]
with |q| = 1. See [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF] pp 16-21 for a detailed definition.

As a fundamental signal processing tool the analytic signal associated to a real 1D signal f (t) is well known to be constructed by setting its Hilbert transform Hf (t) in the imaginary part; being equivalent to a simple spectral operation :

f A (t) = f (t) + iHf (t) ⇔ fA (ω) =          0 if ω < 0 f (ω) if ω = 0 2 f (ω) if ω > 0
Note that the spectrum of f A is null for negative frequencies.

Modulus and argument of f A can be interpreted as instantaneous magnitude and phase. When the original signal is considered oscillating or narrowband these two pieces of data become meaningful -making this tool useful for instance in amplitude and frequency modulation.

Looking at one point of interest, we can interpret a high magnitude as a strong 'presence of some oscillation' around this point. Phase indicates the relative location of this point within the oscillation.

Bülow pointed out that instantaneous phase equivalently describes the kind of structure at a point -among four fundamental structures : 0-'positive impulse', π/2-'falling step', ±π-'negative impulse', -π/2-'rising step'. This is an important point for pattern recognition and that is why the author wanted to generalize it for 2D signals.

Several generalizations to 2D are possible and depend on properties we choose to keep. We consider 2D versions of the Hilbert transform (HT). The partial HT (H θ ) is equivalent to a 1D HT along direction θ. The total HT

(H T = H 0 H π/2
) is the combination of two partial HT's along x and y axes.

A 2D analytic signal of f can be constructed by adding H θ f as the imaginary part; what would cancel out a θ-oriented half-plane in the frequency domain.

But this method depends on an arbitrary choice of the direction -not really adapted for our issue. We can use as well H T ; what cancels out 3 quadrants of the Fourier domain, but the original signal cannot be retrieved from the analytic signal (This is the classical approach in [START_REF] Hahn | Hilbert transforms in signal processing[END_REF]).

So Bülow proposed to embed both partial and total HT's into a quaternionic analytic signal :

f A (x, y) = f (x, y) + iH 0 f (x, y) + jH π/2 f (x, y) + kH T f (x, y)
This method cancels out 3 quadrants of the (quaternionic) Fourier domain and allows retrieval of the original signal. Its polar representation provides 2D local magnitude and phase that can be used to analyze 2D signals -in a similar way as with 1D analytic signal.

The QWT uses this generalization to carry out an 'actual' 2D analytic wavelet transform.

Quaternionic Wavelets : Definition

A standard wavelet transform (DWT) provides a scale-space analysis of an image; yielding a matrix in which each coefficient is related to a 'subband'

(localization in the 2D Fourier domain) and to a position in the image. A 'subband' means both an oscillation scale (i.e. a 1D frequency band) and a spatial orientation (i.e. rather vertical, horizontal or diagonal).

For 2D signals, the QWT provides a richer scale-space analysis than DWT. Its coefficients are analytic according to Bülow's theory and contrary to DWT the magnitude is near shift-invariant. Thus it inherits the magnitude-phase local analysis from the very useful 1D analytic signal. In the one hand the QWT can be viewed like a local '2D Quaternionic Fourier

Transform' (QFT); in the other hand its subbands are '2D Quaternionic Analytic Signals' associated with bandpass filtered versions of the original signal.

Note that usual interpretation of magnitude remains analogous to 1D as it indicates the relative 'presence' of a local geometric element; whereas the local phase is now represented by 3 angles carrying a complete description of this 2D structure.

We start with real separable scaling function φ and mother wavelets ψ D , ψ V , ψ H , and we construct their analytic extensions in H :

Real function : Analytic extension :

ψ D = ψ h (x)ψ h (y) → ψ D + iH 0 ψ D + jH π/2 ψ D + kH T ψ D ψ V = φ h (x)ψ h (y) → ψ V + iH 0 ψ V + jH π/2 ψ V + kH T ψ V ψ H = ψ h (x)φ h (y) → ψ H + iH 0 ψ H + jH π/2 ψ H + kH T ψ H φ = φ h (x)φ h (y) → φ + iH 0 φ + jH π/2 φ + kH T φ
Mathematically, 2D HT's of separable functions are equivalent to 1D HT's along rows and/or columns. So we consider the 1D Hilbert pair of wavelets and scaling functions :

( ψ h , ψ g = Hψ h ) ( φ h , φ g = Hφ h )
and we have for instance :

H 0 ψ V = Hφ h (x)ψ h (y) = φ g (x)ψ h (y)
. The analytic 2D wavelets can then be written in terms of separable products :

ψ D = ψ h (x)ψ h (y)+iψ g (x)ψ h (y)+jψ h (x)ψ g (y)+kψ g (x)ψ g (y) ψ V = φ h (x)ψ h (y)+iφ g (x)ψ h (y)+jφ h (x)ψ g (y)+kφ g (x)ψ g (y) ψ H = ψ h (x)φ h (y)+iψ g (x)φ h (y)+jψ h (x)φ g (y)+kψ g (x)φ g (y) φ = φ h (x)φ h (y)+iφ g (x)φ h (y)+jφ h (x)φ g (y)+kφ g (x)φ g (y)
From a practical point of view, the mother wavelet is related to a quaternionic 2D analytic filter -computable with separable 2D filterbanks. This means the decomposition is heavily dependent on the position of the image with respect to x and y axes (rotation-variance) and the wavelet is not isotropic but the advantage is the easy computation.

Each subband from QWT can be seen as the analytic signal associated with a narrowband part of the image. Shift-invariant magnitude |q| represents elements at any space position in each frequency subband; and phase (ϕ, θ, ψ) describes the 'structure' of these elements. We discuss below the interpretation of the phase.

Dual-Tree Implementation

The QWT is implemented by the Dual-Tree algorithm [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF]] -a filterbank using a Hilbert pair as a complex 1D wavelet -originally performing a CWT as shown Fig. 2 for the 1D case. It is approximately analytic and achieves near shift-invariance with a redundancy of 2:1. These two filterbanks lead to four separable 2D filterbanks that can be combined to fit the QWT definition. Fig. 1 shows an example of such decomposition.

Analysis g 2 h 2 g 1 h 1 ↓ 2 ↓ 2 ↓ 2 ↓ 2 ¤ ¤ ¤ ¤ h h h h C C E E h 2 = H{h 1 }, g 2 = H{g 1 }. Synthesis ↑ 2 ↑ 2 ↑ 2 ↑ 2 g2 g1 h2 h1 + + + 1 2 E : Real part, : Imaginary part

Interpretation and Use

For now, literature is quite poor about QWT and the major difficulty with the use of this transform is the interpretation of the phase.

In [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF] 2D quaternionic analytic Gabor filters are defined and used in a Gabor based texture segmentation. Filtered images are 2D analytic and form a scale-space analysis of the image; from which are extracted magnitudes and local phases at each point to characterize the texture.

A shift theorem similar to this of the well known Fourier transform also stands for the quaternionic Fourier transform defined in [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF] : a shift of the image is equivalent to an offset of the two first terms of phase ϕ and θ. As a result, in a Gabor decomposition ϕ and θ describe small 2D spatial shifts of the coded structure -around the quaternionic coefficient position. Note that in 1D that shift is equivalent to structural information; but 2D structures may be more complex than lines or edges (e.g. corners, T-junctions) and cannot be described by this kind of shift coding. The third term ψ completes the structure analysis and is seen as a texture feature.

Bülow observed that where ψ is around ±π/4 the coded structure is a line or an edge oriented along a diagonal.

We found three recent references [START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF][START_REF] Bayro-Corrochano | The theory and use of the quaternion wavelet transform[END_REF][START_REF] Zhou | Quaternion wavelet phase based stereo matching for uncalibrated images[END_REF] where ϕ and θ are used in disparity estimation. As the QWT performs local QFT's; the shift theorem approximately stands for QWT so ϕ and θ code quite simply spatial shifts of structures.

In another application of [START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF] ('wedgelet' representation)

ϕ and θ are used for wedges position and ψ is used for their orientation.

We In this paper, we only discuss wavelet based methods for this application in order to compare QWT with DWT. Feature extraction via standard wavelet representation of images (DWT) has been widely used in texture analysis [START_REF] Tuceyran | Handbook of pattern recognition and computer vision[END_REF][START_REF] Arivazhagan | Texture classification using wavelet transform[END_REF][START_REF] Lasmar | Multivariate statistical modeling for texture analysis using wavelet transforms[END_REF] for the two last decades. Accordingly, the multiscale analysis provided by the DWT is well adapted to textural images. From each subband one may calculate a mean, standard deviation, energy, mean power, or more sophisticated statistical modeling. Those features -well combinedcan form a powerful texture descriptor.

The Dual-Tree Complex Wavelet Transform (CWT) -a complex extension of the DWT -is motivating by the near shift-invariance of its magnitude and the oriented aspect of its wavelets. It was used in [de Rivaz and [START_REF] De Rivaz | Complex wavelet features for fast texture image retrieval[END_REF]Celik andTjahjadi, 2009, 2011] revealing its efficiency with respect to DWT in texture analysis. The invariance of magnitude to shifts makes the extracted feature independent of the precise location of textural patterns; and so allows a better characterization. We will include CWT based feature vectors of [Celik andTjahjadi, 2009, 2011] in our experiments to compare with literature.

We here propose to further texture analysis using QWT because we keep the advantages due to shift-invariance while adding new information with the ψ phase. The QWT can be regarded as the proper 2D generalization of analytic wavelets since the quaternionic analytic signal is the proper generalization according to Bülow. We use the Dual-Tree based QWT algorithm [START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF] and for comparison a DWT with the well known CDF 9/7 wavelets. The latter are widely used and offer a good separation of frequency components; what is well adapted to texture analysis.

QWT needs special "Q-Shift" filters designed in [START_REF] Kingsbury | Complex wavelets for shift invariant analysis and filtering of signals[END_REF] to form Hilbert pairs. Three different filter sets are possible : the 9-tap, 14-tap and 18-tap Q-Shift filters. In wavelet based texture analysis long filters usually achieve better performance. In this paper we mainly use the 9-tap filters to ensure fairer comparison with 9/7 classical filters; nevertheless the other Q-Shift filters will be treated in the experimental part of the paper.

A L-level decomposition provides 3L subbands for analysis and a lowfrequency subband. We do not use the latter; as low frequency structuresalmost constant -are usually not considered to be a feature of texture.

Texture Databases

In order to provide relevant performance with respect to the huge variety of textures we can find, we use several known texture databases :

• Brodatz database [START_REF] Brodatz | Textures : A Photographic Album for Artists and Designers[END_REF] • We also use Outex TC14 database, which has 68 classes of 60 samples with different lightings but only one orientation.

Experimental Protocol

Recall that our point is the understanding of the QWT. The recognition scheme must be elementary so we can identify the influence of the choice of transform/features on the recognition rate. This is why we propose to In order to evaluate performance we split our available classes into two groups :

• Then recognition rates obtained by using a certain feature extraction method will inform about the quality of this method i.e. the power of description brought by this feature vector for texture classes. Some classes are better recognized than others so we can study what kinds of texture make DWT or QWT more suitable.

Cross validation

Since experimental rates will depend on the chosen training base, we propose to iterate the process 100 times with different selections of the training base (but fixed N A ). We obtain 100 recognition rates for each feature extraction method and each texture class. So our recognition rates will be displayed as average rates with a certain deviation, following the so-called cross validation protocol (see Fig. 3 for example). Thanks to cross validation recognition rates are very stables from an experiment to another (rates vary less than 1% experimentally).

Sample selection

Let us precise how the sample selection is done. We need to separate samples of each class 100 times differently. If we naïvely pick up N A consecutive samples as they are stored in the database, we introduce a bias due to the existing correlation between consecutive samples (e.g. in TC12 basis, consecutive samples have same orientation). So N A samples are randomly selected from each class while the test base is formed by all remaining samples. This process is repeated at each of the 100 iterations.

In the sequel, we first present different feature extraction methods from DWT/QWT magnitude and QWT phase, as well as complex wavelet features from [Celik andTjahjadi, 2009, 2011] for comparison. Then the experimental classification results related to these descriptors are given and analyzed.

Feature extraction

Feature extraction is the automatic processing of several numbers describing texture characteristics. The output is a so-called feature vector that will be considered in a feature space to compare with other texture samples.

Many wavelet based texture classifications use a global measure of energy or statistics from the magnitude; which corresponds to characterize particular bands of the Fourier domain. This type of feature extraction is well adapted to textures that can be viewed as superimposed signals located in the Fourier domain. We will study several methods. First we use magnitude based features that are analogous for both DWT and QWT. Then we propose QWT phase based features that will improve texture description;

and a combination of QWT magnitude and phase. Finally the CWT-based features of [Celik andTjahjadi, 2009, 2011] -also made of global subband measures -are described.

Magnitude based features

After having processed the norm

M ij = |q ij | of wavelet coefficients of a
given image (where i and j are discrete coordinate of a pixel); two measures in each subband are considered.

Energy :

m = 1 E i,j M 2 ij
where E is the total energy of the image minus the low-frequency energy.

This normalization makes m the relative amount of energy in a subband.

Standard deviation (St. dev.) :

m = 1 N i,j (M ij -µ) 2
where N is the number of pixels in the subband and µ is the mean.

QWT phase based features

We dealt with QWT phase in section 2.4. We will first use some phase features only and then we combine magnitude and phase information. Now, how can we use QWT phase to describe textures?

First, ϕ and θ are irrelevant because they inform about position of elements whereas we are interested in their structure. We may think of a matching between the 2D shift coding and some 2D kinds of structure -analogously to the 1D phase -but this interpretation is not clear. In addition, our experiments on ϕ and θ gave very bad results (less than 3% recognition).

So we focus on the third term ψ used by Bülow in a Gabor based texture segmentation process [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF].

In this application we assume that textures are stationary so we do not want to describe spatially precise texture patterns but rather extract global measures. However, we want to describe a part of the behavior of the phase within each subband so a measure such as a mean would not be suitable because it is too global. Experimentally, phase means1 by subband give very bad results (less than 1% recognition).

As we aim to do a fair comparison we won't carry out a high level process such as spatial measures (extrema search, connexity . . . ) that would need to be differently adapted to either DWT or QWT. To make the results as comparable as possible our application must keep simple and we would rather extract a single value from a whole subband as a global feature of the phase spatial behavior. To our knowledge, literature is quite poor about the measure of phase. The Global Phase Coherence [START_REF] Blanchet | Measuring the global phase coherence of an image[END_REF] gives a measure of image sharpness but relies on a too much restricted mathematical model with respect to the topic of the paper that is a first approach toward QWT. A relative phase distribution measure was used with some complex wavelets in [START_REF] Vo | Using phase and magnitude information of the complex directional filter bank for texture image retrieval[END_REF] but is specially defined for this type of transform. The quaternionic phase is a totally different concept than the complex one; so it would be irrelevant to use such a sophisticated measure.

It seems that the simple calculus of the standard deviation (st. dev.)

within a subband would be adapted because it describes a part of the behavior of ψ. Since ψ ∈ [-π 4 , π 4 ] (see Euler angles) we avoid the usual problems about circular data (±π discontinuity, phase wrapping) so there is no ambiguity to calculate angle differences and means.

To improve robustness the ψ-deviation can be weighted by the QWT magnitude. A high magnitude means an important presence of an element while a low value means 'no element' and also provides a numerically unstable phase. So the measure should be more representative by not considering the structure of low magnitude features. The weight function W is the magnitude of the QWT coefficients normalized so the sum within the subband is 1; and is integrated in the standard deviation formula as defined below. Here are the two phase measures we use for feature extraction :

St. dev. : m = 1 N i,j (ψ ij -µ) 2
where µ = 1 N i,j ψ ij , N is the number of pixels in the subband, and (i, j) spans a subband.

Weighted st. dev. :

m = i,j W ij (ψ ij -µ W ) 2
where µ W = i,j W ij ψ ij is the weighted mean of the subband.

Combination of QWT magnitude and phase

Let us now combine the st. dev. QWT magnitude features with the weighted st. dev. QWT phase features (Given experimental results that will be presented in the next section; st. dev. measures are retained for both magnitude and phase).

First we experimentally observed that a simple concatenation of both vectors gives very good performance -better than DWT. Here is a classical point in classification : our two features are not of the same kind. One is an amplitude (∈ R + ) and the other is an angle (∈ [0, π 2 ]). It produces a lack of coherence in the final feature vector as all that terms are seen the same way by the Euclidean distance of the k-NN algorithm. A metric often solves this heterogeneity problem but the different metrics we tried did not improve the process. Even considering the optimal metric with respect to our database; performance is not significantly higher since in practice magnitude and phase features are of the same order of magnitude. So for this paper we finally consider simple concatenation of QWT magnitude and phase features.

To compare with literature we propose to also use CWT. Even though the concept is totally different from QWT the underlying algorithm is quite similar and was used in texture analysis for example in [de Rivaz [START_REF] De Rivaz | Complex wavelet features for fast texture image retrieval[END_REF]Celik andTjahjadi, 2009, 2011].

Complex wavelets based features

These feature vectors are processed from the 6 subbands of the Dual-Tree Complex Wavelet Transform (CWT). In most cases, classical features such as entropy or energy measures as well as statistics are extracted from the magnitude of each subband. In this experiment we choose to use CWT magnitude based feature vector from [START_REF] Celik | Multiscale texture classification using dual-tree complex wavelet transform[END_REF]] that uses variance and entropy. The very recent work in [START_REF] Celik | Bayesian texture classification and retrieval based on multiscale feature vector[END_REF] involves standard deviations and energy measures from both magnitude and phase. Note that in the original paper the vector is further reduced and weighted based on statistics (PCA) -which we did not reproduce because our scheme needs to be elementary so we can identify the contribution of DWT/CWT/QWT (see discussion section 3.2).

It is important to stress that the aim of cited papers is to carry out a good recognition whereas the point of the present paper is the analysis of the QWT. This is why we do not use the whole recognition process described in [Celik andTjahjadi, 2009, 2011] but just the feature extraction bloc. We refer the reader to original publications for detailed expressions.

Now the feature extraction part is done, let us give our experimental results of texture recognition with every feature vectors described above. The next section also brings analysis of interesting types of textures.

Experimental results

We obtained many results, depending on L (depth of wavelet decomposition), chosen feature extraction, value of k, size of training base N A and used database.

Experimentally for most feature vectors best results are obtained with L = 3 and k = 3. The larger the training base the better the recognition, so in this comparative study we consider rather small training bases to put the simulation into a rather hard context. The st. dev. measure globally gives slightly better performance than energy, so we will only consider the st. dev.

feature in the sequel.

The choice of Q-Shift filters in Dual-Tree algorithm has almost no influence on recognition rates (gain is less than 1% recognition). A reason is that this modification does not affect the first scale (see [START_REF] Kingsbury | Complex wavelets for shift invariant analysis and filtering of signals[END_REF]) so 1 3 of features remain the same. So we keep 9-tap filters.

For every result presented in the paper, we keep :

• Depth is L = 3 (9 subbands) so feature vectors are of size 9 for DWT/QWT magnitude and phase features, size 18 for QWT combined magnitude and phase, size 36 for CWT magnitude from [START_REF] Celik | Multiscale texture classification using dual-tree complex wavelet transform[END_REF] and size 72 for CWT magnitude and phase from [START_REF] Celik | Bayesian texture classification and retrieval based on multiscale feature vector[END_REF].

• Number of considered neighbors is k = 3;

• Size of the training base per class is N A = 10 (resp. N A = 30 and N A = 10) for Brodatz (resp. TC12 and TC14) database;

• Magnitude measure is a standard deviation;

• 9-tap Q-Shift filters are used in the Dual-Tree algorithm.

We now present experimental results in terms of recognition rates for DWT/QWT magnitude based features, QWT phase features, combination of QWT magnitude and phase and CWT based features. And works in [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF] and [START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF] show that QWT phase can provide a powerful image analysis; making clear that the QWT is not fully exploited here. Now let us use QWT phase to improve recognition and outperform DWT. provides additional geometric information.

497

We can hardly identify some kind of textures better recognized by either

Combined QWT magnitude and phase

The final feature vector based on QWT magnitude and phase provides very good classification performance with an improvement of about 10% recognition with respect to DWT. We present global recognition rates (mean rate and standard deviation of obtained rates according to cross validation) in table 1. As we cannot present every texture class we selected interesting cases -illustrated Fig. 4 -with recognition rates for respectively DWT magnitude based feature extraction, QWT magnitude, QWT phase, and final QWT feature vector. A sample is displayed under each class of the x-axis. We hope such an illustration will give some insight about visual texture characteristics that we could link with QWT properties. Fig. 5 shows more global results between DWT and QWT to conclude about the general improvement of classification thanks to the QWT. We have about a dozen 'bad' cases but the overwhelming majority of classes are better recognized by QWT and figure 5 is very convincing about its superiority. Now we discuss some particular cases to analyze why the QWT is better or worse than DWT.

A first remark is about texture Brod. 52 that is strikingly better analyzed by the QWT phase feature vector with 100% recognition. The reason is this texture contains only horizontal, vertical and -mostly -diagonal structures that makes ψ very significant.

Most of TC14 classes for which DWT is better are visually similar (see The same remark as in section 5.1 about shift-invariance still holds for textures Brod. 52, TC14 14 and TC14 5.

Note that most textures for which the single QWT magnitude is very worse than DWT (see Fig. 3) are better represented with the final feature vector which outperforms the DWT. This result shows that QWT magnitude and phase are strikingly complementary because both give separately bad results whereas very efficient when combined. Class TC12 20 (Fig. 4) makes a good example.

Finally, this use of quaternionic representation of textures greatly improves an analogous standard wavelet based texture analysis. We proposed a wavelet based feature vector that yields an average recognition rate about 10% better than standard wavelets; by using both QWT magnitude -containing similar information than DWT but with shift-invariance -and QWT phase that provides structural information. Standard DWT does not allow using such a phase -making it a less coherent representation. With this of CWT together with shift-invariance of its magnitude can improve a DWT based analysis but the proposed QWT-phase analysis performs better.

An example of the advantage of CWT directionality is texture TC14 25 (Fig. 6). This texture class exhibits straight lines with an orientation around 75 • . This is precisely the orientation on which one of the 6 CWT subbands is aligned. The CWT is well suited in this very particular case so the recognition performance is much better than DWT/QWT.

However directional analysis brought by CWT is globally not as powerful as the 2D quaternionic phase concept of QWT. The 6 subbands of CWT offer better frequency selectivity than DWT but it still carry some similar "energy" information. On the other hand the fundamental relation between QWT ψ-phase and geometry provides a new kind of additional information that complements "energy" data. Finally our results completely confirm that quaternionic wavelets provide better geometric information than DWT and CWT for a texture analysis purpose.

Conclusion

The Quaternionic Wavelet Transform is a recent improvement of standard wavelets that aims to yield coefficients with a shift-invariant magnitude and a phase containing geometric information. This transform makes use of a A comparison is also made with complex wavelets that have been used in literature in similar works. Advantage of complex wavelets is good directionality; which is confirmed in practice by better recognition for some highly oriented textures. But its complex phase fails to provide such geometric information as this of QWT. As a result QWT outperforms complex wavelets.

Our classification algorithm is a first approach as we focus on studying QWT so there are many ways of improvements in terms of pattern recognition. In particular, spatial analysis of subbands is not used and the classification proposed in [START_REF] Arivazhagan | Texture classification using wavelet transform[END_REF] for example could be improved by using a comprehensive use of QWT phase to provide an excellent classification that could compete with state of the art algorithms.

Our future work implies the monogenic wavelet transform -an improvement of the QWT based on another generalization of the analytic signal : the monogenic signal. This transform is harder to implement but can lead to a more coherent signal analysis thanks to better properties. Monogenic magnitude is rotation invariant and the phase is easier to interpret than QWT

Figure 1 :

 1 Figure 1: Quaternionic wavelet transform of a textural image. From left to right : Original image, magnitude |q| (intensity inverted for visual convenience), ϕ ∈ [-π; π], θ ∈ [-π 2 ; π 2 ], ψ ∈ [-π 4 ; π 4 ]. The 3 terms of phase are represented in color -the hue corresponding to the angle (cyan for 0, red for ±π). Darker zones in phase correspond to negligible magnitude (making phase absurd).

Figure 2 :

 2 Figure 2: 1D analytic wavelet transform with the Dual-Tree algorithm.

  contains 111 photos called textures by the author. We use those [640×640] images to create texture classes, by cutting each one into 25 [128 × 128] sub-images considered to belong to the same class. • Outex databases (http://www.outex.oulu.fi/) were specially made for texture classifiers benchmarking. We use the TC12 database that contains 24 classes of 380 [128 × 128] samples. In each TC12 class we have different images of a same kind of texture at different orientations and under different lightings.

  use a simple k Nearest Neighbors classifier (k-NN) without any automatic feature selection or normalization that would drown the contribution of the transform and confuse the interpretation. Note that we use DWT-based and QWT-based feature vectors that will be presented in the next section. Here is the k-NN algorithm. From each image of a training base a feature vector describing the texture is extracted and labelled by its actual class number. When a test image is given, its feature vector is processed and compared in the feature space to the feature vectors of the training base in terms of Euclidean distance. According to a parameter k, the k nearest vectors are kept to find the most representative class; that is the most represented class within the k neighbors.

  A test base : N T examples the program does not know; • A training base : N A labelled examples the program knows.

Figure 3 :

 3 Figure 3: Magnitude based classification performance : cross validation results for some classes of the 3 databases.

5. 1 .

 1 Magnitude based descriptor : comparing DWT and QWT Some cross validation results for magnitude based feature extraction are illustrated Fig. 3. Recognition rates are expressed in terms of mean, standard deviation and extrema of the rates obtained for each class with the various training bases. Class labels are on the x-axis and recognition rate on the y-axis. For each class the vertical line spans all 100 obtained rates between the minimum and the maximum, the bar is centered on the average rate and is high two times the standard deviation.

5. 2 .Figure 4 :

 24 Figure 4: Cross validation results for some classes of the 3 databases. See also Fig. 3. Performance of each class is shown for the measures : DWT magnitude st. dev. (white), QWT magnitude st. dev. (light-gray), QWT phase weighted st. dev. (dark-gray) and QWT combination of magnitude and phase (black).

Figure 5 :

 5 Figure 5: Final classifier performance for the DWT feature vector and the final QWT feature vector using both magnitude and phase.

Fig. 5

 5 Fig. 5 TC14 55 and 56), almost uniform and hard to distinguish. These images have a pure stochastic behavior so the orthogonality of the DWT is better adapted than the structural analysis brought by QWT -yet this issue is open.

  2D generalization of the analytic signal : a classical powerful tool for signal analysis. This generalization due to Bülow gives a coherent local analysis of 2D signals and was naturally embedded into a wavelet framework to overcome common drawbacks of standard wavelets. Redundancy brought by QWT phase adds complete structural information about local features of images.The QWT is not straightforward to interpret and requires a quite sophisticated theoretical framework. But here we gave an application study crossing the gap between that framework and the way to use this tool; by showing its superiority over standard wavelets in a texture analysis context. With a large texture database our QWT based feature extraction gives globally 10% better recognition than DWT. This improvement is due to shift-invariance of QWT magnitude together with the use of QWT phase that contains useful structural information for texture analysis.

  Note that globally, both transforms have equivalent performance. We do not present results for all classes because they are too numerous but we selected the cases where difference of performance is maximum between DWT and QWT. These extreme cases could give clues about what type of texture is better recognized by either the DWT or the QWT.

	are present. But for other cases, this superiority is experimentally confirmed
	as the structures are better encoded by QWT magnitude than by DWT -
	giving higher recognition rates.
	To conclude about magnitude based feature extraction, global recognition
	performance is quite similar between DWT and QWT. But if we look closer at
	particular textures we see that shift-invariance of QWT improves recognition
	of textures made of 'randomly located similar structures'; and that structures
	aligned with a multiple of 45 • are better represented in QWT magnitude
	while other orientations must have their related information in the phase.
	We propose to analyze the results through particular cases to highlight
	QWT characteristics. Observing Fig. 3, we see that texture Brod. 67 con-
	tains numerous very similar and randomly placed small black disks. Shift-
	invariance of QWT magnitude is essential to code those structures with vari-
	ous positions; whereas each disk is surely quite differently coded in the DWT
	because of small space shifts. In this case QWT gives a more robust descrip-
	tion of texture -keeping everywhere a coherent representation of the disks -
	what explains the superiority of QWT for Brod. 67.
	Let us call 'common' orientations multiples of 45 • i.e. horizontal, vertical
	and diagonal. We can globally see that textures better recognized by QWT
	have many structures right aligned with 'common' orientations. Those for
	which QWT is less efficient rather contain structures oriented between them.
	According to us, this is because QWT phase contains some information
	about these intermediate orientations. So by using only magnitude, superior-
	ity of QWT is not visible at those cases where many 'intermediate' structures

Table 1 :

 1 Final recognition rates.

	498
	499	magnitude or phase but the recognition rates really show that magnitude
	500	and phase are completing one another; so it seems natural to combine both
	501	feature vectors.

The angle is 'unwrapped' in the calculus of a mean. But note that ψ does not show phase wrapping as it lies in [-π/4; π/4].

Acknowledgements

This work is supported by the ANR project VERSO -CAIMAN. We also would like to thank the anonymous reviewer for his/her valuable comments.

texture classification we validate the recent QWT from an application point of view as superior to standard DWT.

We will now present experimental results with CWT. Recall that CWT coefficients are fundamentally different from QWT coefficients in terms of analysis and interpretation of their geometric meaning. Even if both are computed by a separable Dual-Tree algorithm the combination of the 4 related filterbank outputs is different. In the QWT case the 4 outputs are seen like Cartesian components of 3 quaternion valued subbands while for CWT some sums and differences of outputs are considered to be real and imaginary parts of 6 complex subbands -see [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF]. This is why analysis provided by CWT has not much to do with QWT; and we will see in the next section that kinds of texture better recognized by QWT or CWT -as well as global performance -are not the same.

Results with complex wavelets

Results obtained with features from [Celik andTjahjadi, 2009, 2011] have to be studied carefully. Recall that we present it to compare with literature and that we considered only the feature extraction part from cited references.

So result we get are not comparable with the very high performance of original papers partly due to a clever choice of classification algorithm together with proper normalization and feature selection/normalization. In addition size of training bases may differ.

Globally, the first feature descriptor based on CWT magnitude provides better recognition than DWT but worse than final QWT based features; which is an expectable result. Second feature vector based on both CWT magnitude and phase gave very bad recognition (about 5%). We observe that in [START_REF] Celik | Multiscale texture classification using dual-tree complex wavelet transform[END_REF]. We selected cases of biggest difference of performance between DWT/QWT and CWT.

CWT phase measures are numerically unstable because no special weighting is used. In practice many absurd angles -related to the numerous zero coefficients -discredit variance end energy features processed on phase data.

In [START_REF] Celik | Bayesian texture classification and retrieval based on multiscale feature vector[END_REF] a PCA based feature selection is done before classification -what we did not do -so phase based features may be discarded by this stage. Thus we will only study magnitude based analysis from [START_REF] Celik | Multiscale texture classification using dual-tree complex wavelet transform[END_REF].

We show Fig. 6 classification results with the first feature vector [START_REF] Celik | Multiscale texture classification using dual-tree complex wavelet transform[END_REF]