Elliptical Monogenic Wavelets for the analysis and processing of color images Raphaël Soulard and Philippe Carré

Abstract-This paper studies and gives new algorithms for image processing based on monogenic wavelets. Existing greyscale monogenic filterbanks are reviewed and we reveal a lack of discussion about the synthesis part. The monogenic synthesis is therefore defined from the idea of wavelet modulation, and an innovative filterbank is constructed by using the Radon transform. The color extension is then investigated. First, the elliptical Fourier atom model is proposed to generalize the analytic signal representation for vector-valued signals. Then a color Riesz-transform is defined so as to construct color elliptical monogenic wavelets. Our Radon-based monogenic filterbank can be easily extended to color according to this definition. The proposed wavelet representation provides efficient analysis of local features in terms of shape and color, thanks to the concepts of amplitude, phase, orientation, and ellipse parameters. The synthesis from local features is deeply studied. We conclude the article by defining the color local frequency, proposing an estimation algorithm.

Index Terms-wavelet transform, color, analytic signal, monogenic signal, monogenic wavelets.

I. INTRODUCTION

Through the last 15 years, the analytic signal and its 2D generalizations have greatly improved wavelet representations [START_REF] Chan | Coherent multiscale image processing using dual-tree quaternion wavelets[END_REF], [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF], [START_REF] Selesnick | The dual-tree complex wavelet transform -a coherent framework for multiscale signal and image processing[END_REF], [START_REF] Unser | Multiresolution monogenic signal analysis using the Riesz-Laplace wavelet transform[END_REF] by naturally embedding 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. The signal is encoded with more sparsity than with standard wavelets, which globally improves most waveletbased applications. The seemingly most appropriate 2D extension for image analysis is the monogenic representation [START_REF] Felsberg | Low-level image processing with the structure multivector[END_REF], of which a few wavelet counterparts have been proposed [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF], [START_REF] Olhede | The monogenic wavelet transform[END_REF], [START_REF] Unser | Multiresolution monogenic signal analysis using the Riesz-Laplace wavelet transform[END_REF] and applied [START_REF] Belaid | Phase-based level set segmentation of ultrasound images[END_REF], [START_REF] Joshi | The segmentation of colorectal MRI images[END_REF], [START_REF] Li | Robust stereo image matching using a two-dimensional monogenic wavelet transform[END_REF], [START_REF] Maltaverne | Motion estimation using the monogenic signal applied to ultrasound elastography[END_REF], [START_REF] Zhang | Monogenic-LBP: A new approach for rotation invariant texture classification[END_REF].

Existing monogenic wavelet transforms are reviewed in section II, and it turns out that from the numerical point of view, the topic is still at its very first steps.

Our main concern is that there is no satisfactory synthesis method from the monogenic wavelet analysis, whereas this is required for a large category of image processing applications, such as denoising, enhancement, compression and watermarking. While studying the invertibility and the synthesis monogenic wavelets, the information encoded in coefficients will be highlighted.

In this paper, we define a clear synthesis method from monogenic coefficients through the idea of wavelet modulation. An original numerical algorithm of reversible monogenic XLIM (UMR CNRS 7252) University of Poitiers, France. filterbank is derived in section III. The scheme is based on the Radon transform and 1D undecimated analytic wavelets. The reconstruction is exact and the synthesis wavelets are modulated according to monogenic coefficients.

The second part of the paper investigates the color extension in section IV. Analyzing color data is essential for a lot of applications but the Fourier analysis and wavelet representations have been almost exclusively limited to scalar signals so far. In practice, processing of color images is most often based on a suboptimal use of scalar tools, either by only considering luminance, or by working separately on each color channel ("marginal" schemes [START_REF] Bovik | Handbook of Image and Video Processing[END_REF]).

We propose to mathematically handle the dimensionality of color signals by starting back to the univariate case. The elliptical oscillation model is proposed to generalize the analytic signal to vector-valued signals. This method inspired from [START_REF] Lilly | Modulated oscillations in three dimensions[END_REF] provides a new clearly interpretable color phase concept.

The color definition is then extended to 2D. We define a directional color Riesz transform -based on the structure tensor formalism -to build the elliptical monogenic signal model. We end up with a color wavelet transform, whose implementation is easily done in parallel by directly using our Radon-based algorithm.

The proposed wavelet representation provides efficient multiscale representation of local color features in terms of amplitude, phase, orientation, and ellipse parameters. Elliptical monogenic features open up good applicative prospects by coherently embedding AM-FM and multiscale representations, as well as differential geometry, into a unified tool for color images. In particular, the color phase feature is a totally new concept.

The synthesis from coefficients is deeply studied, and two reconstruction methods are finally given so that most known wavelet-based processing application are possible.

In section V, the elliptical monogenic wavelets are successfully used to estimate local frequency from color images, by using the multiscale amplitude and phase data. This experimentally validates the color phase concept in particular and the proposed filterbank in general, opening the way towards higher-level image analysis applications for future work.

For making this research reproducible, the source code that produces the figures and the experiments presented in this paper has been made available publicly [START_REF] Soulard | MATLAB source code for the implementation and use of Elliptical Monogenic Wavelets[END_REF].

Following our previous publication [START_REF] Soulard | Vector extension of monogenic wavelets for geometric representation of color images[END_REF] on color extension of monogenic wavelets, the original contributions of the present article are the new filterbank algorithm, the synthesis method, the enriched color signal model (elliptical) and the local frequency estimation.

Notations :

2D coordinates: x = (x, y) ω = (ω 1 , ω 2 ) ∈ R 2 Euclidean norm: x = x 2 + y 2 Complex algebra: z = {z} + j {z} = |z|e j arg(z) Fourier transform: F{f }(ω) = f (t)e -jωt dt

II. THE GREYSCALE MONOGENIC REPRESENTATION

A. Generalization of the analytic signal An analytic signal s A is a multi-component signal associated to a real signal s to analyze. The definition is classical for scalar 1D signals [START_REF] Gabor | Theory of communication[END_REF], [START_REF] Ville | Théorie et applications de la notion de signal analytique[END_REF]: s A = s + jHs, where the Hilbert transform is defined by:

{Hs} (t) = p.v. 1 π t s(t -τ ) dτ F ←→ -j sgn(ω) ŝ(ω) (1)
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.

The growing interest on this tool within the image community is due to the possible interpretation of amplitude, phase and frequency in terms of local geometric shape.

Among the different existing 2D extensions [START_REF] Bülow | Hypercomplex spectral signal representation for the processing and analysis of images[END_REF], [START_REF] Felsberg | The monogenic signal[END_REF], [START_REF] Hahn | Multidimensional complex signals with single-orthant spectra[END_REF], the monogenic signal [START_REF] Felsberg | Low-level image processing with the structure multivector[END_REF] seems to be the most advanced definition. It is rotation invariant and generalizes the local phase concept intuitively through a directional model. Given a 2D real signal s, the associated monogenic signal s M is 3vector valued (instead of complex-valued in the 1D case) and must be taken in spherical coordinates:

s M =   s {Rs} {Rs}   =   A cos ϕ A sin ϕ cos θ A sin ϕ sin θ   ( 2
)
where Rs is the complex-valued Riesz transform of s:

{Rs}(x) = p.v. τ 1 +jτ 2 2π τ 3 s(x -τ ) dτ F ←→ ω 2 -jω 1 ω ŝ(ω) (3)
(see also [START_REF] Larkin | Natural demodulation of two-dimensional fringe patterns: I. general background to the spiral phase quadrature transform[END_REF], [START_REF] Nabighian | Toward a three-dimensional automatic interpretation of potential field data via generalized Hilbert transforms: Fundamental relations 780786[END_REF]). From the analysis point of view, the monogenic signal is formed by the three following features:

Amplitude: A = s 2 +|Rs| 2 Orientation: θ = arg{Rs} ∈ [-π; π[ 1D Phase: ϕ = arg{s + j|Rs|} ∈ [0; π] (4)
A direct link between the angles θ and ϕ and the geometric local structure of s has been shown [START_REF] Felsberg | Low-level image processing with the structure multivector[END_REF]. The signal is so expressed like an "A-strong" 1D structure with orientation θ. The quantity ϕ is analogous to the 1D local phase and indicates if the structure is rather a line or an edge. As a drawback, non-directional structures are not handled. Evolutions of the monogenic model [START_REF] Fleischmann | Image analysis by conformal embedding[END_REF], [START_REF] Wietzke | The signal multi-vector[END_REF] taking into account more complex shapes are out of the scope of this paper.

From a signal processing viewpoint, such an AM/FM representation is accordingly well suited for narrow-band signals.

It is natural to embed it in a multiscale transform.

B. Monogenic Wavelets: a review

Defining monogenic wavelet transforms is a recent topic. The main idea is to combine parallel filterbanks whose underlying wavelet functions form a Riesz triple. This is how, in the 1D case, the "dual-tree" algorithm [START_REF] Selesnick | The dual-tree complex wavelet transform -a coherent framework for multiscale signal and image processing[END_REF] combines good properties of Gabor/Morlet analytic wavelets [START_REF] Grossman | Decomposition of Hardy functions into square integrable wavelets of constant shape[END_REF] (invariance of amplitude) with the sparsity of bi-orthogonal filterbanks -achieving small redundancy and reversibility. However, its 2D extension is not monogenic and fails to extract the actual phase of the signal.

The monogenic curvelet frames [START_REF] Storath | Directional multiscale amplitude and phase decomposition by the monogenic curvelet transform[END_REF] -in the spirit of the 2D Gabor transform [START_REF] Daugman | Two-dimensional spectral analysis of cortical receptive field profiles[END_REF] -use a set of directional subbands. Such a combination of the intrinsic orientation analysis of the monogenic framework with a multi-oriented scheme is out of the scope of this paper.

Continuous monogenic wavelets were proposed in [START_REF] Olhede | The monogenic wavelet transform[END_REF], but does not allow any reconstruction of input signal from coefficients. Yet reversibility is necessary for use in applications such as denoising and compression, for which wavelets are known to be efficient. We will rather focus on wavelet frames and reversible filterbanks.

The Riesz-Laplace monogenic wavelets [START_REF] Unser | Multiresolution monogenic signal analysis using the Riesz-Laplace wavelet transform[END_REF] are constructed upon a real-valued primary wavelet frame relying on a dyadic pyramid, whose mother wavelet ψ is defined in the context of polyharmonic splines. The frame is the extension of a basis obtained by adding a specific subband regression algorithm at the synthesis side [START_REF] Unser | The pairing of a wavelet basis with a mildly redundant analysis via subband regression[END_REF]. The "Riesz part" is obtained from the Riesz transform of ψ, generating an other similar frame (complex-valued). 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.

The choice of the authors is to provide the "minimallyredundant wavelet counterpart of Felsberg's monogenic signal". This transform is approximately isotropic and has a parameter γ controlling the number of vanishing moments.

Exactly isotropic tight monogenic wavelet frames can be found in [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF], [START_REF] Unser | A unifying parametric framework for 2D steerable wavelet transforms[END_REF], where the wavelet is defined in the Fourier domain by a radial profile.

In these works, the algorithm uses FFT implementation of the filterbank. Note that what we call "monogenic filterbank" is in fact a pair of filterbanks -each one being reversibleforming a multiscale monogenic analysis as 3-valued coefficients. The synthesis of the input signal from these coefficients is now discussed.

C. The synthesis issue

We propose to clearly define how the input image should be coherently synthesized from the monogenic analysis. This is not trivial because both the primary and the Riesz-part filterbanks have been separately defined as redundant and perfectly reversible (see e.g. [START_REF] Chenouard | 3D steerable wavelets in practice[END_REF], [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF], [START_REF] Unser | A unifying parametric framework for 2D steerable wavelet transforms[END_REF], [START_REF] Unser | Multiresolution monogenic signal analysis using the Riesz-Laplace wavelet transform[END_REF]). The spatial redundancy of each one is sometimes exploited through the concept of wavelet frames and/or 'subband regression' [START_REF] Unser | The pairing of a wavelet basis with a mildly redundant analysis via subband regression[END_REF]. However, in the monogenic setting, an additional redundancy occurs by the use of two such filterbanks, each one being basically self-sufficient in terms of reconstruction.

For applications involving wavelet-based image synthesis, the modified monogenic features (A , ϕ , θ ) can be easily converted back to Cartesian terms with equation (2) -providing separately the wavelet coefficients for the primary and the Riesz-part filterbanks. However, one does not know if it is better to reconstruct either from the primary part only, or from the Riesz-part, or from a combination of both.

The literature on monogenic representation almost exclusively relies on analysis applications: stereo analysis [START_REF] Li | Robust stereo image matching using a two-dimensional monogenic wavelet transform[END_REF], motion estimation [START_REF] Maltaverne | Motion estimation using the monogenic signal applied to ultrasound elastography[END_REF], segmentation [START_REF] Belaid | Phase-based level set segmentation of ultrasound images[END_REF], [START_REF] Joshi | The segmentation of colorectal MRI images[END_REF], texture classification [START_REF] Zhang | Monogenic-LBP: A new approach for rotation invariant texture classification[END_REF]. A particular case is image demodulation [START_REF] Unser | Multiresolution monogenic signal analysis using the Riesz-Laplace wavelet transform[END_REF], that produces a 'reconstructed' image. The result is in fact equal to the amplitude A, which is part of the analysis, so it does not need any synthesis from the filterbanks. In [START_REF] Chenouard | 3D steerable wavelets in practice[END_REF], which studies the Riesz-wavelets with examples of inverse problems, the authors present a special monogenic setting requiring to jointly use an additional isotropic filterbank. However, this setting is used for contour detection only, the other applications being based on the Riesz-wavelets.

In all these references, no global monogenic synthesis method is given. Exceptions can be found in [START_REF] Unser | A unifying parametric framework for 2D steerable wavelet transforms[END_REF] and [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF].

In [START_REF] Unser | A unifying parametric framework for 2D steerable wavelet transforms[END_REF], the monogenic setting is a particular case of multicomponent generalized Riesz-based steerable wavelet transform, where the primary part is realized with a 0-th order Riesz transform. As a special case of a general inversion scheme, the monogenic synthesis is defined as the sum of the two reconstructions (up to some constant factors). However, this is neither discussed nor applied. A more developed presentation of monogenic synthesis can be found in [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF].

1) Reconstruction from primary part: In [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF], after processing, it is proposed to update the i-th scale primary coefficients w i according to the new amplitude and phase by w i (x) = A i (x) cos(ϕ i (x)) and then to use the synthesis side of the primary filterbank only. Now we give a practical example where this method is not satisfactory. Suppose that an image containing sharp edges is denoised by thresholding the amplitude of its monogenic representation. Whether there is an edge around x 0 for example, w i has a zero-crossing structure -since the primary wavelet acts like a Laplacian; whereas Rw i has a peak-like structure -since the Riesz-wavelet acts like a gradient. As a result, the amplitude A i shows a peak centered on x 0 . The problem is that thresholding will probably cancel out the whole neighborhood around x 0 , leaving only one non-zero coefficient in the Riesz-part Rw i (x 0 ), and nothing in the primary part (w i (x 0 ) = 0). When using only the primary part like in [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF], the information necessary to retrieve the edge is lost. Actually, the isotropic wavelet is optimal for peaks and lines, while the Riesz wavelet is optimal for edges. Contrary to the method from [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF], we believe that both filterbanks should be used so as to properly retrieve contours. More precisely, the synthesis method must be defined such that one monogenic coefficient (A, ϕ, θ) is able to modulate the synthesis wavelet, according to its local amplitude, phase and orientation.

2) Our proposition: To this end, we propose a direct analogy with the Fourier theory. Recall that the Fourier transform synthesizes any real-valued 2D signal by a sum of plane waves of the form:

a ω (x) = A cos(ω T x + ϕ) (5)
Orientation and frequency of a ω are determined by the coordinate ω, while amplitude and phase are encoded by the Fourier coefficient f (ω) = Ae jϕ . We want to point out that the modulation of the plane wave by A and ϕ is mathematically achieved by the summation of odd and even parts:

a ω (x) = { f (ω)} cos(ω T x) -{ f (ω)} sin(ω T x) = A cos(ϕ) cos(ω T x) -A sin(ϕ) sin(ω T x) (6)
In the monogenic wavelet case, the primary filterbank computes the "real part" of the monogenic coefficient, and is responsible for the "even" (isotropic) part of the monogenic wavelet. On the other hand, the Riesz filterbank computes the 'imaginary part' (as two components) and controls the orientation of the "odd" (directional) part of the wavelet. This analogy justifies the synthesis method to consist in a summation of reconstructed signals from both the primary and the Riesz filterbanks. This way, the implicit linear combination of even isotropic wavelets and odd directional wavelets will actually produce the desired shape: a wavelet whose amplitude, phase and orientation are encoded in the monogenic coefficient (A, ϕ, θ). Note that a strictly analogous summation occurs at the core of the Riesz filterbank as well, to steer the directional wavelet by a linear combination of its x and y components (although it may be "hidden" by the use of complex algebra like in eq. ( 3)).

In the 1D case, most works on analytic wavelets -including the dual-tree filterbank -intrinsically define the synthesis within the complex algebra, by the same odd/even summation we just proposed. This needed to be clarified in the monogenic case for which complex numbers cannot be used the same way. Note that our proposition turns out to be similar to the definition from [START_REF] Unser | A unifying parametric framework for 2D steerable wavelet transforms[END_REF].

The review of existing monogenic wavelet transforms reveals that the literature is still at its very first steps. The filterbank design in [START_REF] Unser | A unifying parametric framework for 2D steerable wavelet transforms[END_REF] uses all the degrees of freedom for the tuning of the analysis filter to be as isotropic as possible, ending up with highly anisotropic synthesis filters that could not be used for proper monogenic synthesis. We propose to explore an alternative strategy in the next section.

III. RADON TRANSFORM BASED MONOGENIC FILTERBANK

We propose in this section an innovative monogenic wavelet transform algorithm based on the following motivations. The monogenic framework simplifies to a 1D definition in the Radon domain [START_REF] Brackx | On generalized Hilbert transforms and their interaction with the Radon transform in clifford analysis[END_REF]; discrete Radon transform algorithms are available in the literature [START_REF] Carré | Discrete analytical ridgelet transform[END_REF]; and the design of analytic 1D wavelets is well know compared to 2D. In particular, the dual-tree algorithm [START_REF] Selesnick | The dual-tree complex wavelet transform -a coherent framework for multiscale signal and image processing[END_REF] is an interesting implementation of mildly redundant analytic wavelets. We propose to investigate a Radon based monogenic filterbank design.

In the first part we present the theoretical relation between Radon transform, analytic and monogenic wavelets. Next we present the numerical algorithm for the Radon transform, and then for the monogenic transform. Note that this section only deals with greyscale images, color will be treated right after.

A. Monogenic wavelets in Radon domain

Given any function s(x, y), its Radon transform is defined by the projection along direction ρ as :

s ρ (t) = R s( τ sin ρ + t cos ρ , -τ cos ρ + t sin ρ )dτ (7)
The Radon domain is known to connect the 1D and 2D worlds, as suggested by the Fourier slice theorem:

s ρ (ω) = s(ω cos ρ, ω sin ρ) (8)
It shows that 1D Fourier transform in the Radon domain is like 2D Fourier transform in the space domain. This property is used in image reconstruction from projection methods. The convolution property stems from it. Given two signals s(x, y), h(x, y) and their Radon transforms s θ (t), h θ (t), we have:

(s θ * h θ )(t) = (s * * h) θ (t) (9)
where * denotes 1D convolution and * * denotes 2D convolution. This is convenient to perform isotropic filtering. By applying a 1D symmetric filter h 1 on all Radon projections of a 2D signal, an isotropic 2D filtering is actually performed:

s θ (t) ← (s θ * h 1 )(t) ⇔ ŝ(ω) ← ŝ(ω) ĥ1 ( ω ) with h 1 (t) = h 1 (|t|) (10)
The resulting 2D filter has a radial frequency response identical to ĥ1 (ω) 1 . So Radon domain filtering allows a direct control of the frequency response from 1D filters.

In addition, the Riesz transform turns out to be simplified in the Radon domain, as a weighted Hilbert transform [START_REF] Brackx | On generalized Hilbert transforms and their interaction with the Radon transform in clifford analysis[END_REF]:

{Rs} ρ (t) = e jρ {Hs ρ }(t) (11)
Considering that a monogenic wavelet transform essentially consists in some isotropic band-pass filtering combined with the Riesz transform, the properties given above strongly suggest that it can reduce to using 1D band-pass filters and the Hilbert transform i.e. an analytic wavelet transform. This way we can focus on a simpler 1D filterbank design, while letting the isotropy issue be handled by the Radon transform. Given any symmetric 1D analytic wavelet:

ψ A (t) = ψ(t) + jHψ(t) with ψ(t) = ψ(|t|) (12)
let the monogenic wavelet ψ M be defined by:

ψ M (x) =   ψ iso (x) {Rψ iso (x)} {Rψ iso (x)}   Radon -→   ψ(t) cos(θ)Hψ(t) sin(θ)Hψ(t)   (13)
with the following radial frequency response:

ψ iso (ω) = ψ( ω ) (14)
From this construction, let us derive the numerical algorithm. 1 There is no such closed-form relation in the space domain.

B. Discrete Radon transform

The discretization of the Radon transform is difficult to achieve. Most methods in the literature have been devised for computerized tomography or to approximate the continuous formula. None of them, however, were specifically designed to be invertible transforms for discrete images. In previous work [START_REF] Carré | Discrete analytical ridgelet transform[END_REF] we have developped a numerical Radon transform with exact reconstruction based on discrete analytical geometry. This approach has straightforward n-D extension and computational simplicity. The present filterbank is based on it.

Thanks to the Fourier slice theorem (eq. ( 8)), we know that Radon projections can be obtained from inverse 1D FFT's applied on all slices of the 2D FFT of the input signal.

The main idea is to treat each direction θ by a discrete line, having central symmetry and forming a good approximation of the corresponding Euclidean straight line. We want the line to cross both the origin and a given "target point" p = (p 1 , p 2 ) ∈ Z 2 , tied to orientation θ = arg(p 1 + jp 2 ). The work in [START_REF] Carré | Discrete analytical ridgelet transform[END_REF] proposes to write this line L as a set of discrete 2D coordinates, according to the following definition:

L τ p = (k 1 , k 2 ) ∈ Z 2 |p 2 k 1 -p 1 k 2 | ≤ τ 2 ( 15
)
where τ is the arithmetical thickness parameter. Thickness provides control over the connectivity of the line, and will determine the redundancy factor of our Radon transform. In this paper, we will work with 8-connected "naive" lines, by setting τ = max(|p 1 |, |p 2 |). We actually tried "thicker" lines ("pythagorean", "supercover" [START_REF] Carré | Discrete analytical ridgelet transform[END_REF]) but besides increased redundancy, we experimentally noticed no significant difference. The algorithm generating the discrete path works iteratively with 1-pixel shifts from p to (0, 0). At each step, the shift is selected among the 8 possible directions from eq. ( 15). Path from (0, 0) to -p is then processed by central symmetry. This way, we ensure that coordinates are well ordered with respect to their radial frequency ω .

We can now extract all Fourier slices to perform Radon transform. Due to the desired symmetry, the dimension must be odd. We consider an input signal of size

[2K 1 +1]×[2K 2 + 1], and its 2D FFT ŝ(k) with k ∈ [-K 1 ..K 1 ] × [-K 2 ..K 2 ].
To ensure exact reconstruction, the set of lines must cover the whole lattice. Thus the set of orientations is defined so that every boundary point

p ∈ [-K 1 ..K 1 ] × [-K 2 ..K 2 ] with |p 1 | = K 1 or |p 2 | = K 2 (16
) together with its symmetric -p, is the target point for one line L τ p . We end up with 2K 1 +2K 2 +1 Fourier slices covering the whole 2D FFT and generating a discretization of θ in [0; π[.

After line extraction, a 1D inverse FFT provides the final Radon projection.

Finally, the reconstruction has to be defined. After having transformed the Radon projections back in the 1D Fourier domain, the 1D Fourier coefficients must be set back into the 2D spectrum. The same discrete 2D paths used at the extraction stage are used again to match Radon domain coordinate to Fourier domain. A choice has to be made since most 2D Fourier samples are represented through several Fourier slices (redundancy). The solution proposed in [START_REF] Carré | Discrete analytical ridgelet transform[END_REF] is a pointwise averaging of the multiple contributions from the different slices.

Although this method is satisfactory in the general case, it is not adapted for the monogenic filterbank because the Riesz transform needs to be implemented by a function depending on θ. In this case, for one Fourier sample, most of the contributions are tied to values of θ that do not fit well the angular coordinate of the Fourier sample, resulting in wrong values. Therefore the contributions must not be averaged.

We propose an alternative reconstruction method in this paper. Every Fourier sample ŝ(k) will be extracted only once, from the best possible slice s θ (k). This slice will be chosen to have θ value as close as possible to the exact orientation of the Fourier sample: θ ref = arg(k 1 + jk 2 ). Note that when the Radon transform is not modified, both reconstruction methods are identical, but the present one is required for the monogenic filterbank. The last inversion step is simply the inverse 2D FFT. The whole reconstruction algorithm is numerically exact up to machine's precision.

C. Discrete Monogenic wavelets

Now we have to choose a 1D analytic wavelet transform to be applied in the Radon domain. The first idea is to use the well established dual-tree wavelet transform [START_REF] Selesnick | The dual-tree complex wavelet transform -a coherent framework for multiscale signal and image processing[END_REF], that combines good approximation of Hilbert pairs of filters, and efficient filterbank implementation. However, this method would downsample the Radon projections at each scale, which would produce half-length projections without reducing the number of orientations. The discrete Radon algorithm is not compatible with such resizing of the data, so the inverse Radon transform could not be used. We will then rather focus on undecimated wavelet transforms.

This choice presents the advantage of releasing strong constraints on filter design, and potentially improving the signal analysis. This also implies increased redundancy, nonorthogonal decomposition, and that we found no existing work considering a Hilbert pair of undecimated filterbanks. Literature on analytic wavelets is either focused on subsampled IIR filterbanks (Phaselets [START_REF] Gopinath | The phaselet transform-an integral redundancy nearly shift-invariant wavelet transform[END_REF], B-Splines [START_REF] Chaudhury | Construction of Hilbert transform pairs of wavelet bases and Gabor-like transforms[END_REF]), or on continuous definitions, like the generalized Morse wavelets [START_REF] Lilly | Generalized morse wavelets as a superfamily of analytic wavelets[END_REF] and the well known Gabor/Morlet wavelets [START_REF] Grossman | Decomposition of Hardy functions into square integrable wavelets of constant shape[END_REF]. A known decimated method may be adapted by simply removing the sub-sampling steps but it would limit us to less optimal filters. Therefore, we propose the customized filterbank design outlined on figure 1.

s i-1 fi gi ghi s i (t) w i (t) w hi (t) f i g i g hi l + 1 2 l + s i-1
Fig. 1. Undecimated analytic 1D filterbank iterated for i ∈ {1..J}. Input is s = s 0 , wavelet coefficients are w i + jw h i and the approximation is s J .

The perfect reconstruction of input signal is ensured by:

s i-1 = s i-1 ⇔ fi fi + ( ĝi ĝi + ĝhi ĝhi )/2 = 1 (17)
To achieve the Radon domain monogenic analysis explained in section III-A, the primary low-pass (resp. high-pass) filters f i and fi (resp. g i and gi ) must be symmetric in the space domain, and the Hilbert part (anti-symmetric) must verify ghi = Hg i . For convenient interpretation of wavelet coefficients, we choose to use the same filters on the analysis and synthesis stages. With these constraints, the perfect reconstruction condition [START_REF] Fleischmann | Image analysis by conformal embedding[END_REF]

reduces to | fi | 2 + |ĝ i | 2 = 1.
We propose the following definition of the filters in the Fourier domain, based on a Gaussian low-pass filter:

fi (t) = f i (t) F ←→ fi (ω) = exp -ω 2 σ 2 i /2 gi (t) = g i (t) F ←→ ĝi (ω) = 1 -fi (ω) 2 ghi (t) = -g hi (t) F ←→ -j sgn(ω) ĝi (ω) (18)
Like in [START_REF] Chaudhury | Construction of Hilbert transform pairs of wavelet bases and Gabor-like transforms[END_REF], we opt for Fourier domain definition and FFT implementation, allowing for high order IIR filters. In this case the Hilbert transform is approximated by the truncation ω ∈] -π; π]. Note that odd-length signals will be preferred (this is always verified with our symmetric discrete Radon projections) to avoid singularity at Nyquist frequency ω = ±π.

We choose the dyadic dilation rule

σ i = 2 i-2 so that f i (t) = 1 2 f i-1 ( t 2 )
(idem for g i and g hi ). Subbands are normalized by the energy of equivalent filters at each scale. We refer the reader to the source code for further details [START_REF] Soulard | MATLAB source code for the implementation and use of Elliptical Monogenic Wavelets[END_REF]. Now we can define the final monogenic filterbank. The combination of discrete Radon transform, analytic filterbank and θ weighting according to eq. ( 11) is outlined on the figure 2. Note that FFT-based 1D filtering allows to work directly on Fourier slices and skip all the 1D FFTs.

The analysis stage basically consists in Radon domain analytic wavelet transform followed by cos(θ)-sin(θ) weighting and inverse Radon transform.

At the synthesis stage, back in the Radon domain, the two Riesz components are multiplied again by cos(θ)-sin(θ) and summed so that we get back the imaginary part of the analytic wavelet transform. With respect to discussion at section II-C, where we suggest to sum monogenic components after filtering, we must point out that the orientation information of each monogenic coefficient will here still properly modulate the orientation of the synthesis wavelet. The reason is that the steering will be actually done in the last inverse Radon transform. It is shown on figure 3 that the values of the monogenic coefficient actually modulate the wavelet shape in terms of phase and orientation. It is also visible that our Radon-based monogenic wavelet is not exactly isotropic. This is due to the Fourier slices having varying Nyquist frequencies from π at xand yaxes to π/ √ 2 at diagonals. The 1D filters are applied regardless of this bias, which intrinsically dilates

ϕ = 0 θ = 0 ϕ = π 2 θ = 0 ϕ = π 2 θ = π 3 ϕ = π 3 θ = π 3 ϕ = π 2 θ = -π 6
s i-1 FFT 2D discrete slice extract. iFFT 1D Radon transform f i g i -g h i Analyt. WT cos θ sin θ inverse Radon inverse Radon inverse Radon inverse Radon s i w i w x,i w y,i Radon Radon Radon Radon cos θ sin θ j + f i g i g h i j + 1 2 j + Inverse analytic WT FFT 1D discrete slice fitting iFFT 2D Inverse Radon transform s i-1
Fig. 2. Radon based monogenic 2D filterbank (iterated for i ∈ {1..J}). Input image is s(x) = s 0 (x). Monogenic wavelet coefficients at i-th scale are (w i (x), w x,i (x), w y,i (x)). The low-frequency approximation is s J (x).

their frequency response around diagonals, in the 2D FFT, resulting in a slightly cross-shaped wavelet. This points out the tight link between the Radon algorithm and isotropywhich has always been a crucial point in the discrete world.

Improvements may be possible but out of the scope of this paper, while the present result is perceptually coherent with respect to phase and orientation, which is satisfactory.

D. Discussion

Note that our filterbank is close to the algorithm of complex ridgelets [START_REF] Chen | Image denoising with complex ridgelets[END_REF]. In the ridgelet analysis, the goal is to construct oriented wavelets. Its complex extension is aimed at gaining translation invariance, by using analytic 1D wavelets in the Radon domain, resulting in Radon domain wavelet coefficients. In the present work, the constructed wavelets are the multiscale extension of an AM-FM model. Their definition requires the subsequent use of the inverse Radon transform, to set the coefficients back in the space domain. This particularity is necessary for the interpretation of the data as a monogenic analysis. For example, the 2D phase concept does not exist with the complex ridgelets, whereas it is natural in the monogenic case.

Compared to the existing monogenic pyramids [START_REF] Held | Steerable wavelet frames based on the Riesz transform[END_REF], [START_REF] Unser | Multiresolution monogenic signal analysis using the Riesz-Laplace wavelet transform[END_REF], our undecimated algorithm benefits from theoretical simplicity, and especially from a clear synthesis method, achieving coherent modulation of the wavelets. Let us recall that the representation is translation and rotation invariant. The n-D extension can be straightforwardly done by extending the discrete line concept like in [START_REF] Helbert | 3-D discrete analytical ridgelet transform[END_REF] and considering n -1 orientation angles. Now our greyscale monogenic filterbank is defined, let us extend the whole concept to the tri-valued bi-dimensional signals that are color images.

IV. ELLIPTICAL MODEL FOR COLOR MONOGENIC REPRESENTATION

In a previous work [START_REF] Soulard | Color monogenic wavelets for image analysis[END_REF], following the lead from [START_REF] Demarcq | The color monogenic signal: Application to color edge detection and color optical flow[END_REF], we proposed a color monogenic analysis based on Cauchy-Riemann equations in the Clifford algebra. However, from the "image processing" point of view, this construction does not provide a clear physical interpretation of the data.

We have next carried out an alternative construction [START_REF] Soulard | Vector extension of monogenic wavelets for geometric representation of color images[END_REF] based on differential geometry to handle the color orientation concept. This time, the amplitude and phase were clearly identified, and a "color axis" concept emerged. Yet, the synthesis method has to be improved, since only the primary part of the filterbank was used (see discussion in II-C).

This paper goes one step further toward color signal processing, starting from the fundamental concept of Fourier atom to construct a full color signal model. The interpretation of the features is improved thanks to a true vector approach of the basic signal processing concepts, and the synthesis method in terms of modulated color wavelets is achieved.

Note that we would think premature to consider advanced color spaces at this stage, so we restrict ourselves to Euclidean ones -like the usual RGB space. The study of perceptually meaningful color spaces and of the dependencies between colour channels is an important prospect.

A. Fourier atom

Redefining Fourier analysis for multidimensional signals s : R p → R q is an open mathematical topic, usually involving quaternions or Clifford algebra (see [START_REF] Batard | Spinor Fourier transform for image processing[END_REF], [START_REF] Hitzer | Quaternion and Clifford Fourier Transforms and Wavelets[END_REF] for example). From the physical point of view, this literature does not yet provide sufficiently intuitive interpretation of Fourier coefficients. We propose a pragmatic approach.

As a construction basis, it seems reasonable to consider separate Fourier analyzes on R, G and B channels. The color signal is then expressed like a sum of sine waves that we will call Fourier atoms. In the 1D case, we define the color Fourier atom a(t) as:

a(t) =   a R (t) a G (t) a B (t)   =   A R cos(ωt + ϕ R ) A G cos(ωt + ϕ G ) A B cos(ωt + ϕ B )   (19)
For any fixed frequency ω, the parameters of a(t) are given by three Fourier coefficients:

â = A R e jϕ R A G e jϕ G A B e jϕ B T (20)
An example of color Fourier atom is given on figure 4a.

Here is our starting point for a non-marginal analysis of color signals. As is, those three amplitudes and three phases do not intuitively convey the oscillation's behavior in the "color world". We may prefer some unified amplitude and phase concepts, as well as color features. The atom a(t) must then be studied as a point-trajectory within the 3D color space -a "color path" -as illustrated on figure 4b.

It can be shown [START_REF] Lilly | Modulated oscillations in three dimensions[END_REF] that a(t) always draws an ellipse in the color space. We propose to study the ellipse parameters to physically handle the atom's behavior in a non-marginal way.

[ A R A G A B ] = [ .7 .18 .36 ] [ ϕ R ϕ G ϕ B ] = [-.63π -.1π -.03π ] 2π -1 2 0 1 2 E t

B. The elliptical model

New "Fourier-like features" can be carried out by studying the color Fourier atom as an ellipse. First, the amplitude concept can intuitively be extended to the ellipse's size. By considering the apogees as sinusoid's maxima, the phase can then be defined as the time-relative position of the apogees, with respect to a(0), as illustrated in figure 4c.

The apogees can be quickly identified analytically from the analysis of a(t) , whose minimum and maximum values must be equal to the least radius r -and the greatest radius r + (see figure 4c). It can be verified that:

a(t) 2 = A 2 -|â T â| 2 + |â T â| cos 2 ωt + arg(â T â) 2 (21)
where the amplitude A is defined by

A = A R 2 + A G 2 + A B 2
(Amplitude) [START_REF] Helbert | 3-D discrete analytical ridgelet transform[END_REF] This writing explicitely provides the extrema of a(t) 2 :

r 2 -= (A 2 -|â T â|)/2 r 2 + = (A 2 +|â T â|)/2 (
23) By defining the phase as

ϕ = arg(â T â)/2 (Phase) (24)
the point a(-ϕ/ω) is one of the two apogees, as expected, since its norm equals r + . The classical scalar case is retrieved with â = Ae jϕ , r -= 0, r + = A.

The quantity âT â can be viewed as a combination of the three Fourier coefficients, so that ϕ results as a weighted average of the separate phases. The modulus |â T â| is an ellipse feature in itself, when expressed through the "linearity":

λ = |â T â| A 2 = r 2 + -r 2 - r 2 + + r 2 - ∈ [0; 1] (Linearity) (25)
The linearity λ is as close to 1 as the three separate phases are similar (modulo π). In this case the ellipse has a thin shape, and the Fourier atom oscillates between two main colors by crossing the origin. When λ is small, the ellipse is round and draws a richer color path. Now the ellipse's position must be identified. We propose to define the apogee by:

a + = a -ϕ ω =   A R cos(ϕ R -ϕ) A G cos(ϕ G -ϕ) A B cos(ϕ B -ϕ)   = r +   cos α 1 cos α 2 sin α 1 cos α 2 sin α 2   ( 26
)
We consider the direction of a + as the main color axis, encoded by α 1 and α 2 . The angles are computed by appropriate arc-tangent operations. Then the perigee is defined by:

a -= a -ϕ-π 2 ω =   A R sin(ϕ R -ϕ) A G sin(ϕ G -ϕ) A B sin(ϕ B -ϕ)   = R α1 R α2   0 r -cos α 3 r -sin α 3   (27)
with the rotation matrices:

R α1 =   cos α 1 -sin α 1 0 sin α 1 cos α 1 0 0 0 1   R α2 =   cos α 2 0 -sin α 2 0 1 0 sin α 2 0 cos α 2   (28)
The perigee's direction is encoded by the last parameter α 3 . Its computation requires to use the rotation by R T α2 R T α1 that sets a + on the x-axis and a -in the y-z plane. α 3 is the angle between the y and z coordinates of the rotated perigee R T α2 R T α1 a -. The back conversion from (A, λ, ϕ, α 1 , α 2 α 3 ) first consists in getting back r + and r -:

r + = A (1 + λ)/2 r -= A (1 -λ)/2 (29)
and rotating them to retrieve a + and a -from equations ( 26)- [START_REF] Joshi | The segmentation of colorectal MRI images[END_REF]. The original separate Fourier coefficients are equal to:

â = (a + + j.a -)e jϕ (30)
Let us finally give angle wrapping details. At this stage, two equivalent ellipse feature sets are possible for one given Fourier atom, according to the equivalence:

(A λ ϕ α 1 α 2 α 3 ) ≡ (A λ ϕ + π α 1 + π -α 2 -α 3 ) (31)
In order to have a unique encoding, we propose to do

if cos(ϕ R -ϕ) < 0 then ϕ ← ϕ + π (32)
before the computation of a + and a -. This corresponds to arbitrarily choosing the apogee being in the half of the color space where the red channel is positive. The advantage is that two oscillations having similar color properties will be quickly identified as having similar values of α i , even if they are out of phase each other. We finally have (α 1 , α 2 ) ∈ [-π/2; π/2] 2 and (ϕ,

α 3 ) ∈ [-π; π] 2 .
The proposed ellipse parameters are strongly inspired from [START_REF] Lilly | Modulated oscillations in three dimensions[END_REF], where the elliptical model has been studied to analyze trajectories of particles in the field of seismology, and our definitions of A, λ and ϕ are in fact identical.

However, the approach is different with respect to the encoding of the ellipse's position as well as the calculus. In [START_REF] Lilly | Modulated oscillations in three dimensions[END_REF], the method relies on the computation of a normal vector, orthogonal to the plane containing a(t). Then two specific rotation angles ("azimuth" and "zenith") allow to study the ellipse in the x-y plane, so that r + , r -, ϕ and a last "precession" angle are defined. We believe that when working in a color space, the choice of these particular angles is not justified. In contrast, our definition provides a notion of color axis which is analogous to the so-called color axis proposed in [START_REF] Soulard | Vector extension of monogenic wavelets for geometric representation of color images[END_REF], as pointing toward two opposite main colors involved in the oscillation. When λ = 1, the color axis from present paper and this from [START_REF] Soulard | Vector extension of monogenic wavelets for geometric representation of color images[END_REF] are equal. When the oscillation is more complex (λ < 1), the color axis of present paper encodes a still meaningful main color direction, while this from [START_REF] Soulard | Vector extension of monogenic wavelets for geometric representation of color images[END_REF] becomes less relevant because its is based on a simpler model. It should also be noted that the ellipse parameters computation from [START_REF] Lilly | Modulated oscillations in three dimensions[END_REF] is unstable when λ ≈ 1. The reason is that the normal vector is ill-defined when λ = 1. In contrast, our method relying on the computation of the apogee and perigee, is numerically more efficient. The only natural consequence when λ = 1 is that α 3 is irrelevant because the perigee is null (r -= 0). Finally, we believe that our approach provides a simpler computation of the conversion between the Fourier coefficients and the ellipse features.

Thanks to the elliptical model, the three Fourier analyzes can be converted into intuitive amplitude and phase completed by 4 colorimetric features λ, α 1 , α 2 , α 3 in a non-marginal way. We can now define the corresponding color analytic signal.

C. Color analytic signal

The analytic signal (see II-A) makes a pointwise estimation of amplitude and phase, under the hypothesis that the signal to analyze locally resembles a Fourier atom a(t) = A cos(ωt + ϕ). The extracted local information has the form of a Fourier coefficient Ae jϕ . In the case of an input 3-valued signal

s(t) = [ s R (t) s G (t) s B (t) ] T (33)
we define the color extension by the use of independent Hilbert transforms on R, G and B channels, followed by the conversion to elliptical parameters. The color analytic signal is then

s A (t) = [ s R (t) s G (t) s B (t) Hs R (t) Hs G (t) Hs B (t) ] T ⇒ [ A(t) ϕ(t) λ(t) α 1 (t) α 2 (t) α 2 (t) ] T (
34) Now the ellipse parameters are time varying. We plot an example of oscillating color 1D signal with varying properties, together with the estimated amplitude and phase on figure 5. The synthetic color path is shown in the RGB space as a 3d plot on the left-hand side, as three red, green and blue 1d plots on the upper graph, and as a pixel map on the middle. Colors are generated by centering the values of s(t) in the RGB space around the middle grey. The estimated amplitude A(t) is plotted as a dotted curve in the upper graph. It visually well conveys the local energy of the oscillation. The estimated instantaneous phase ϕ(t) is plotted in the lower graph. Because the frequency of input signal is supposed to evolve smoothly, the phase is expected to be locally linear, which is confirmed here. The slope of ϕ is locally according with the frequency of s. Note the singularity of ϕ where amplitude equals 0, which is usual with polar representations, and has no consequence as long as A is used as some significance information.

The linearity λ(t) (not plotted) also gives the right values:

• λ ≈ 0.5 for the first oscillation (blue-yellow), whose trajectory in the RGB space is visibly "round", as can be seen on the left-hand 3D plot; • λ ≈ 1 for the second oscillation (pink-green), whose trajectory is far much linear. The other colorimetric parameters α 1 , α 2 and α 3 are constant in the first part, and slowly varying in the second part, according to the color patterns visible on the pixel-map. The interpretation of (α 1 , α 2 ) is similar to this of the color axis proposed in [START_REF] Soulard | Vector extension of monogenic wavelets for geometric representation of color images[END_REF]. This axis points towards the two opposite main colors involved in the oscillation. More generally, the interpretation of color angles is still an open issue, and we expect to get more intuition about it by working on the 2D extension.

The reader should notice that the phase estimation is well behaved through the four quadrants, and not altered by the variation of the other parameters. The quality of the estimation of ϕ is a direct benefit from the use of elliptical modeling. With respect to the phase estimation from our previous work [START_REF] Soulard | Vector extension of monogenic wavelets for geometric representation of color images[END_REF], this phase data is numerically more stable.

Let us now combine this new formalism with the Riesz transform to carry out the elliptical monogenic representation.

D. Color Riesz transform

We have studied above how to analyse a color signal by combining it with a phase-shifted version of its color channels (its Hilbert transform) in the 1D case, in perfect analogy with the analytic signal formalism. The extension to images must consider the monogenic formalism, according to which the signal is combined with its Riesz transform ("RT"). For the extension to be well driven, we have to consider the directional RT, defined for any scalar 2D signal s(x) by: R θ s(x) = cos(θ) {Rs(x)} + sin(θ) {Rs(x)} (35)

= |Rs(x)| cos( θ -arg(Rs(x)) ) ( 36
)
Here the RT acts like a pair of steerable filters, of which a linear combination can provide any rotation of them. The RT is retrieved by combining them along xand yaxes2 Rs = R 0 s + jR π 2 s. The directional RT is equivalent to a

s s 3 {Rs 3 } {Rs 3 } -π/2 π/2 θ + R θ+ s 3 0 1 A 0 π/2 ϕ 0 1 λ -π/2 π/2 α 3 -π/2 π/2 α 1 -π/2 π/2 α 2
Fig. 6. Color monogenic wavelet transform at third scale (Subband s 3 ). Oscillating color signals are centered and normalized in the RGB cube. Amplitude A and linearity λ are shown as normalized greyscale bitmaps. Angles θ + , ϕ, α 1 , α 2 and α 3 are displayed as hue in the HSV color space, while saturation is controlled by A so as to whiten irrelevant values (The linearity is also used for α 3 ). The angles ϕ and α 3 are wrapped. The smoothing kernel h for calculus of θ + is a Gaussian filter with σ = 2.

Hilbert transform in the case of strictly directional 2D signals (a.k.a. "simple signal" or "intrinsically 1D"):

s θ1 (x) = s 1 (cos(θ 1 )x 1 + sin(θ 1 )x 2 ) (37) ⇒ R θ1 s θ1 (x) = {Hs 1 }(cos(θ 1 )x 1 + sin(θ 1 )x 2 ) (38)
This expresses well that the monogenic signal is like the 1D analytic signal, locally steered toward its main direction. We can now feed the elliptical model with the directional RT as phase-shifted versions of s. The last step is to choose a direction, toward which the RT will be steered.

In the grayscale case, the direction is naturally arg(Rs). Here, we propose to use the one that maximizes the color variation. This mathematical problem is solved in terms of Euclidean color distances, by the well known color gradient from Di Zenzo [START_REF] Zenzo | A note on the gradient of a multi-image[END_REF], [START_REF] Sapiro | Anisotropic diffusion of multivalued images with applications to color filtering[END_REF]. Thanks to a strong theoretical link between the RT and the gradient operator [START_REF] Köthe | Riesz-transforms versus derivatives: On the relationship between the boundary tensor and the energy tensor[END_REF], [START_REF] Unser | Multiresolution monogenic signal analysis using the Riesz-Laplace wavelet transform[END_REF], we have already proposed in [START_REF] Soulard | Vector extension of monogenic wavelets for geometric representation of color images[END_REF] a Riesz-based color structure tensor, analogous to Di Zenzo's gradient. It is defined by:

T R (s) = T R R + T G R + T B R = T 11 T 12 T 12 T 22 ( 39
)
with

T C R = h * {Rs C } 2 {Rs C } {Rs C } {Rs C } {Rs C } {Rs C } 2 ( 40
)
and h is a smoothing filter. The main orientation of color local structures is equal to the orientation of the eigenvector tied to the largest eigenvalue of T R :

θ + = arg(T 11 -T 22 + 2jT 12 )/2 ( 41
)
This method is used in the grayscale case in [START_REF] Unser | Multiresolution monogenic signal analysis using the Riesz-Laplace wavelet transform[END_REF] to regularize the orientation with h. Now we use θ + to steer the directional RT so as to perform amplitude/phase analysis in a unique local direction. The phase-shifted signal is defined by:

R T s(x) = R θ+(x) s R (x) R θ+(x) s G (x) R θ+(x) s B (x) (42)
The equation ( 34) can now be fed by the three channels of the signal s on the one hand, and by those of the phase-shifted version R T s on the other hand. Let us finally use it to define the elliptical monogenic representation.

E. Elliptical monogenic wavelets

The elliptical monogenic representation of any color input 2D signal s is given by: 22) to ( 27)

s M =         s R (x) s G (x) s B (x) Rs R (x) Rs G (x) Rs B (x)               T R -→ θ + (x) →   R θ+ s R R θ+ s G R θ+ s B   E eqs. (
¡ ¡ ¡ ¡         A(x) ϕ(x) λ(x) α 1 (x) α 2 (x) α 3 (x)         (43)
The conversion to ellipse parameters is possible thanks to the steering of RTs toward the local main orientation θ + based on the Riesz-part tensor constructed from s M . Note that the construction basically uses three parallel Riesz transforms, before the non-marginal conversion of all the components. The wavelet representation can therefore be based on three parallel monogenic filterbanks. The color decomposition is defined by the use of greyscale monogenic wavelets on every color channel followed by the conversion of output coefficients according to eq. [START_REF] Soulard | MATLAB source code for the implementation and use of Elliptical Monogenic Wavelets[END_REF]. In this paper we use the filterbank defined in section III.

Let us comment the illustration of the wavelet transform at the third scale given on figure 6. We first comment the upper row. The subband s 3 seems well filtered in an isotropic way, with no visibly favored direction. The Riesz part actually shows xand yphase-shifting, and the tensor orientation θ + is according to our perception, and equally efficient on lines and edges thanks to the smoothing kernel. The Riesz-part is well combined as R θ+ , where directional structures are clearly phase-shifted in the proper direction (See how the blue/yellow edge is curved to follow the input yellow line).

Let us now observe the lower row, containing the ellipse parameters. The amplitude and phase are visibly coherent. Color lines and edges are equally detected by the amplitude A, that shows a maximum at their centers. The phase is wrapped in one quadrant of the trigonometric circle to highlight the lineedge information. We actually find 0 (red) and π/2 (blue) on line-and edge-like structures respectively, which confirms that ϕ encodes contour type information. The linearity of local color ellipse λ is near 1 almost everywhere. This means that the chosen test image mainly contains "simple" color oscillations, involving two main colors at their apogees in the color space. For example, the yellow line in test image is encoded in the subband by a simple variation between blue and yellow The information conveyed by all these features is further explored below in the context of synthesis. Note that this result confirms the advantage of having defined a numerically stable computation of ellipse parameters, instead of using the method in [START_REF] Lilly | Modulated oscillations in three dimensions[END_REF], which is ill-defined when λ = 1. Another advantage is that the angles (α 1 , α 2 ) can be interpreted as the main direction in the color space, pointing towards the two opposite colors around a contour (see also [START_REF] Soulard | Vector extension of monogenic wavelets for geometric representation of color images[END_REF]). Note how the two angles show locally constant values around contours of a same color.

Note that features appear independent of each other, suggesting a complementary encoding of the visual information. Before going into the synthesis method for this wavelet transform, let us briefly discuss some related work.

Related work: The color quaternion Fourier transform

It should be noted that elliptical color paths had already been noticed in the context of color quaternion Fourier transform by Ell and Sangwine [START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF]. Yet, those ellipses are different from these of present paper, because they arise from a very different combination of real and imaginary parts of red, green and blue Fourier transforms. The "modulus" and "phase" from [START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF] are then different quantities than those presented here. We think that the quaternion framework -intrinsically 4D -is impeding to handle well the required 6 ellipse parameters. Moreover, the ellipses from [START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF] are parameterized by a globally constant color "axis", responsible for the particular combination of color channels, whereas in this paper, the analysis in nonparametric.

The prior work by McCabe et. al. [START_REF] Mccabe | Theory of spatiochromatic image encoding and feature extraction[END_REF] should also be cited, because it revealed elliptical paths from the use of a classical Fourier transform on image's chrominance.

Let us now study image synthesis from our monogenic coefficients. We expect monogenic synthesis to provide more interpretative elements about the features.

F. Synthesis

Because this research field is motivated by the global idea of improving existing wavelet transforms, we believe that the famous applications involving reconstruction -compression, contrast enhancement, denoising, watermarking... -must not be put aside, contrary to most works about analytic/monogenic wavelets. The synthesis wavelet modulation process proposed in section II-C must now be extended to the elliptical feature set, so that color and shape information can be transferred to the synthesis atoms.

Thanks to the parallel structure of the model, reconstruction from Cartesian components is straightforwardly defined by the parallel use of synthesis monogenic filterbanks. As shown on figure 7, the final synthesis method relies on the way we choose to coherently retrieve the Cartesian coordinates after having processed the elliptical monogenic features. At

s R s G s B MWT MWT MWT         s R i s G i s B i Rs R i Rs G i Rs B i         E eq. ( 43
)           A i ϕ i λ i α 1,i α 2,i α 3,i θ +,i           E ?         s R i s G i s B i Rs R i Rs G i Rs B i         iMWT iMWT iMWT s R s G s B
Fig. 7. Color monogenic analysis and synthesis filterbank. The ? symbol stands for "How to get coherent Cartesian coordinates back from the monogenic features". this stage of the work, the color monogenic signal does not perform a "split of identity" i.e. we have no reversible polarlike representation of the whole 9D feature set s M . Actually, we are close to it, since the ellipse parameters provide a 6D split of identity of primary and directional Riesz part. But the orientation information of the Riesz part -reduced to θ + -is not fully included in the model. Fortunately, coherent synthesis is possible in some way. We here propose two particular methods. The first one is direct and exact but limited to amplitude-based processing, the second one allows modification of all ellipse parameters but with approximate reconstruction.

1) Amplitude-based processing: Many wavelet-based image processing methods such as wavelet shrinkage (denoising) and contrast enhancement only modify the amplitude of coefficients. Since the amplitude A is defined as the Euclidean norm of the monogenic vector, applying any weighting k to A is equivalent to applying this weighting on every of its Cartesian components:

A ← kA ⇔ s M ← k [s R s G s B Rs R Rs G Rs B ] T (44)
This way, the effect of any processing of amplitude is properly propagated to the Cartesian monogenic representation and then to the synthesis image. For example, a wavelet-shrinkage color image denoising scheme can be carried out with:

k(x) = 1 if A(x) > T 0 else ( 45
)
where T is a given threshold. An experimental result is shown on figure 8, featuring a noisy version of the "sailing" picture, and the denoised result. This experiment is the clear illustration that the color monogenic wavelet transform proposed in this paper is able to coherently synthesize images. More precisely, at those coordinates x i where amplitude is significant, the 9 Cartesian components of s M (x i ) are preserved. The parallel monogenic filterbanks are then fed with the full information to properly combine odd and even wavelet shapes and produce the corresponding color modulation. To highlight the contribution of embedding the color dimension in our model, the corresponding marginal approach is illustrated for comparison. It consists in applying an amplitude-based greyscale monogenic thresholding on every color channel separately (The threshold has been experimentally adapted to compare as fairly as possible). It can be seen that the marginal approach is less stable with respect to false colors. See the red artifacts around the sail's number, and the green artifacts around the contour of the black-red sail. They are produced by some large noisy wavelet coefficients, which happen to be better regularized in a true color approach. With the marginal method, the unbalanced processing of the color channels generates a strong quantification-like bias in the encoding of the color information. Note how this rounding of color features in the RGB space produces "primary color" artifacts, here red and green, and visually annoying. In contrast, the proposed approach either discards or preserves all channels as a whole, which tends to preserve accurate local color information. Note that the development of a complete denoising experiment would require a probabilistic model aimed at discriminating the information from the noise, which is out of the scope of this paper.

This wavelet transform is in the spirit of sparse representations, because we expect the visual information to be succinctly encoded by few significant coefficients, and closely to the human visual system. Given the rich color and shape description contained in one only monogenic coefficient, we may even think that a good synthesis could still be obtained by keeping only those whose amplitude is a local maxima, like in [START_REF] Mallat | Characterization of signals from multiscale edges[END_REF] (This topic is a future work).

In order to use the elliptical monogenic wavelet representation, we have to further understand the link between features and visual elements. Let us now study the shape of synthesis wavelets and their modulation by the elliptical features, through an alternative synthesis method.

2) The pure directional case: The full color monogenic filterbank proposed in this paper is tied to a very large dictionary of color wavelet shapes, parameterized in geometry and color by a 9D feature set. Nevertheless, according to the directional model upon which it is based, its proper interpretation requires Original Reconst.

Original Reconst.

Original Reconstruction Fig. 9. Pure directional reconstruction of the sailing image (original on figure 8) and the face image (original on figure 6), and zoomed parts.

a unique local orientation, encoded by θ + . In some cases, s M can be such that the Riesz parts have different orientations across color channels, resulting in a biased interpretation of ellipse parameters. We propose to ignore the multi-oriented cases by assuming:

θ + (x) = arg (Rs R (x)) = arg (Rs G (x)) = arg (Rs B (x)) ( 46
)
Such assumption is usual in vector differential geometry, for example with Di Zenzo's gradient. This restriction provides a new synthesis method where the directional Riesz part is steered back to θ + with Rs C = e jθ+ R θ+ s C for C ∈ {R, G, B}. We obtain:

        A(x) ϕ(x) λ(x) α 1 (x) α 2 (x) α 3 (x)         θ + (x) T eqs. (22) to (27) c   R θ+ s R R θ+ s G R θ+ s B   E ×e jθ+ E r r j               s R (x) s G (x) s B (x) Rs R (x) Rs G (x) Rs B (x)         (47)
The dictionary is reduced to the subset of pure directional wavelets, so the reconstruction is not exact in general. It is still a good approximation, since we experimentally noticed almost no visible artifact with our test images, as shown on figure 9. The only visible differences are very mild and occur around the sharp synthetic edges from image face. It turns out that most local structures actually fit the directional model, so they can be properly retrieved.

The advantage is that the synthesis wavelets can now be Examples of modulated synthesis wavelets are given on figure 10. This experiment confirms that the modulation is properly achieved in terms of phase, orientation and color content.

The first column illustrates the special subset of isotropic wavelets, for which ϕ = 0 (the structure is symmetric) and λ = 1 i.e. the elliptical color path reduces to a line between two opposite colors. For these cases, the orientation is not used and set to an arbitrary value. We can see that the produced wavelets are approximately isotropic and actually feature two opposite colors, explicitly indicated in the color space by the color axis (α 1 , α 2 ). Since λ = 1, the angle α 3 has no influence.

The other columns imply directional wavelets, for which ϕ = π/2 indicates an edge-like shape. For the second and third columns, we can verify the correspondence between θ + and the perceptible orientation of the colored wavelets. here again, only two main colors are involved because the ellipse reduces to a line (λ = 1).

The colorimetric pattern becomes richer when the elliptical color path is rounder (λ < 1), as can be seen in the last two columns. In these cases, α 3 encode an "orthogonal color axis" containing the information of the additional colors present in the oscillation. The smaller λ, the more present these new colors. In the first row for example, some brown has been added to the lime-purple oscillation. These wavelets may not be involved for retrieving simple edges but rather for color textures. This suggests using the linearity feature λ to discriminate contours (λ ≈ 1) from textural points (λ ≈ 0).

This alternative synthesis scheme may be used for image regularization, where not only the amplitude could be thresholded, but the other elliptical features could also be regularized. For example, the phase can be simplified under some hypothesis about its local linearity, or the color angles can be made locally constant etc.

Finally, let us quickly deal with the open issue of multioriented cases. By handling three independent orientations instead of just θ + , it is possible to produce multi-oriented synthesis wavelets, that do not fit the directional model, like illustrated in the figure 11. In the greyscale computer vision (0, π/2, 0) (0, π/3, 2π/3) (0, π/5, 2π/5) literature, these shapes are referred to as "intrinsically 2D" and include corners, junctions, crossings etc. A prospect of this work is to extend the model so that it includes these shapes as well. This could involve geometric notions of isotropy and coherency -well known in the structure tensor literature [START_REF] Tschumperlé | Vector-valued image regularization with PDEs: A common framework for different applications[END_REF] and may be linked to the conformal monogenic representation [START_REF] Fleischmann | Image analysis by conformal embedding[END_REF], [START_REF] Wietzke | The signal multi-vector[END_REF]. The challenge is to preserve the phase interpretation, which has been so far strictly tied to directional models. Such a research lead could end up with a 9D split of identity associated to a unique modulating reconstruction scheme.

Now the synthesis has been defined and the elliptical monogenic features explained, let us conclude this article in a color image analysis context, through the definition and the estimation of the local color frequency.

V. LOCAL COLOR FREQUENCY ESTIMATION

It is well known in the signal processing literature that the local frequency ν can be basically obtained by differentiating the local phase ϕ. This topic is numerically challenging and several methods can be used [START_REF] Huang | On instantaneous frequency[END_REF]. In the 2D world of images, the topic has a long history in the AM-FM literature [START_REF] Bovik | Handbook of Image and Video Processing[END_REF] where local frequency is interpreted as some granularity feature.

The monogenic representation presented in this paper gives direct access to local phase in a multiscale fashion. This is the first time to our knowledge that phase data can be obtained from a color image. This motivates the proposition of a local frequency estimation method based on our filterbank. The first issue is numerical phase differentiation.

A. Differentiating the phase

After having decomposed any color image into elliptical monogenic subbands thanks to the filterbank of section IV-E, we need to compute the spatial derivative of ϕ(x).

An analytic method is given with the greyscale monogenic filterbank from [START_REF] Unser | Multiresolution monogenic signal analysis using the Riesz-Laplace wavelet transform[END_REF]. This consists in using an other version of the filterbank, where the Riesz part filters are tuned to correspond to partial derivatives of the Riesz transform. The coefficients can then be combined to analytically find the local frequency. In our color case, extending the analytical formula of the derivative of ϕ(x) is more intricate, so we adopt another approach: the finite central difference method.

The phase is locally understood toward the direction θ + , and the differentiation must of course be applied in this same direction. In the discrete world, we can only compute horizontal, vertical or diagonal differences. To steer the calculation properly, our continuous orientation data θ + must be quantized so as to select the neighbor "left" ϕ l and "right" ϕ r samples to be compared to the local value ϕ(x). We propose the following quantization formula

θ q = (π/4) round (θ + /(π/4)) (48)
giving 4 possible directions: horizontal (0), vertical (π/2 and -π/2 are merged) and the two diagonals (±π/4). For example, if θ + lies in [-π 8 ; π 8 ]; then the wavefront is approximately parallel to the y-axis, and the phase evolves mainly along the x-axis. So ϕ(x 1 , x 2 ) has to be differentiated with ϕ l = ϕ(x 1 -1, x 2 ) and ϕ r = ϕ(x 1 + 1, x 2 ).

The circularity of angle data can be handled with complex exponentials, and the local frequency defined by: ν = arg (exp(ϕ r )/ exp(ϕ l ))

1 2 / cos(θ + -θ q ) (49)
The right side term compensates the bias due to orientation quantization, according to a local linear model of ϕ. Note that for diagonal cases (θ q = ±π/4), ν must also be divided by √ 2. It shall also be noted that in practice, the phase difference is not directly estimated between ϕ r and ϕ l , but computed as the average of the phase differences on both sides, which prevents some phase wrapping. For more computation related details, see the source code [START_REF] Soulard | MATLAB source code for the implementation and use of Elliptical Monogenic Wavelets[END_REF].

B. Multiscale estimation

The advantage of working in the wavelet domain is that we can locally select the "best scale" -which is the one containing the natural frequency of the signal i.e. showing the largest energy. The pointwise best scale k(x) is defined by:

k(x) = arg max i {A i (x)} (50)
In order to improve the estimation's robustness, the final algorithm is then as follows:

• Process color monogenic wavelet transform;

• Differentiate the phase ϕ i (x) at all scales i to compute local frequency ν i (x); • Compute the local best scale k(x); • Apply 3 × 3 median filtering to k(x) so that the used scale is locally constant for numerical stability;

• Extract the local frequency ν(x) = ν k(x) (x)
The performance of this method is illustrated on figure 12. This first test is done on a synthetic image, featuring explicit local frequencies for various examples of colorimetric and geometric parameters. We expect our method to measure constant values of ν within every spatial sub-square, as shown on the second graphic. The measure (third graphic) is strikingly close to theoretical values, with some expected and not damaging instability around borders. The estimation method is satisfactory for low frequencies down to ν = 0.35 (period of about 18 pixels, lower-left wave) and up to ν = 2.5 (≈ 2.5 pixels, upper-right wave). Efficiency is regardless of orientation and color features, which experimentally confirms the soundness of the method.

The multiscale property of our method is expected to bring robustness with respect to noise. The figure 12 shows a noisy version of the same test image, and the estimated local frequency (fifth graphic). It is visible that the result is stable, which is confirmed by the good quality of estimation obtained by a simple further regularization of ν with a 5 × 5 median filter (sixth graphic). The advantage of using a multiscale local frequency estimation is therefore confirmed here.

We finally test it on a natural image as shown on figure 13. This image features a perceptually constant unique frequency, ν oscillating horizontally. The difficulty is that colorimetric properties are various -blue, green and yellow. We also expect measurements not to be as regular as in the synthetic case, in particular due to richer frequency content, and randomness of real data. The estimated ν(x) has been regularized by a 11×11 median filtering.

The figure 13 shows the estimated frequency (center graphic), together with a histogram of the measured ν(x) (right graphic). Apart from a mild random variation due to the original image itself, we can see that ν still lies within a small interval around ν = 0.51 rad./pixel. This value corresponds to a period of around 12 pixels, which is well in accord with the visible periodicity of the pattern in this image.

This experiment shows that our elliptical monogenic representation actually provides sound data, according to an oscillatory model in the spirit of the classical Fourier analysis, with a intuitive extension of the amplitude and phase concepts. This opens the way towards many possible applications.

Texture segmentation of course could benefit from the color frequency feature, as well as the local color linearity λ which may discriminate complex color patterns from simple color contours. The ellipse angles α 1 , α 2 and α 3 can also provide simple discrimination. Phase data has also been shown useful for texture classification in [START_REF] Soulard | Quaternionic wavelets for texture classification[END_REF]. Note that the development of a full segmentation method based on the monogenic features would require the design and normalization of optimal feature vectors, as well as the modelling of meaningful distances, which is out of the scope of this paper.

The elliptical monogenic features make new invariant descriptors of local color and shape for keypoint detection. The invariance properties and the signal flavor of this work could make it a competitive color alternative to the famous SIFT detector [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], probably closer to the human visual system.

Stereo vision and motion estimation are also possible applications, since the phase data can be used to detect sub-pixel accurate shifts (as already done in the greyscale case).

VI. CONCLUSION

This article investigates the numerical representation of color images by monogenic wavelet transforms.

Monogenic representations have been gaining much interest for providing rich physical meaning with respect to usual wavelet transforms. This family of sparse representations suffers from a lack of synthesis algorithms as well as color extensions. Applications have been so far limited to greyscale image analysis only, putting aside the famous wavelet-based denoising and compression as well as any color handling.

To overcome this limit, we have proposed in this paper a clear synthesis method, that we exploited through the definition of a whole new monogenic filterbank. Our numerical algorithm has an explicit synthesis part, and is then able to produce synthesis wavelets being coherently modulated by the amplitude, phase and orientation encoded in one monogenic coefficient. The construction is based on a discrete Radon transform and an original 1D analytic undecimated wavelet transform. This is the first time that the synthesis of the 2D signal from its monogenic wavelet coefficients is presented, thanks to the concept of modulated synthesis wavelets.

The second part of the paper tackles the color extension. The monogenic representation is already multi-dimensional in the scalar case, and its extension to vector-valued signals is not trivial. We start in 1D by defining the color analytic signal, revealing the fundamental elliptical oscillation and the color phase concept. The representation is then brought to 2D thanks to a directional approach, making use of the color structure tensor to build the directional Riesz transform. The elliptical monogenic representation is carried out as a non-marginal extension of the monogenic framework, tied to a directional color elliptical oscillation model. It achieves local analysis of images in terms of color and geometry features, including amplitude, phase, elliptical parameters and orientation.

The monogenic filterbank derived in the first part is used to implement the new color wavelet transform. We observe a good separation of geometrical and colorimetric information through the various local features. The synthesis is again deeply studied, and two methods are proposed. Exact reconstruction is possible in the case of amplitude-based processing, allowing for a variety of applications such as wavelet-based denoising and enhancement. For the price of an approximate reconstruction -which is experimentally not visible -a full color wavelet modulation is achieved from the elliptical feature set, opening the way to new kinds of image processing tasks such as phase based regularization or any automatic local modification of graphical elements, by tuning their elliptical properties.

Finally, the color monogenic representation is applied to estimate from a color image the local frequency information, in a multiscale fashion. The method appears to be robust against additive noise, and is successfully used on synthetic and natural images.

The elliptical monogenic wavelet transform provides a ready-to-use multiscale pointwise low-level analysis combining classical computer vision methods and new kinds of information, while being efficiently designed for color images. This suggests many prospective applications.

We think that a SIFT-like keypoint detector could be based on this filterbank, and bring better connections with the human visual system, as well as better handling of the color information. Segmentation may also benefit from the advanced encoding of color contours and texture-like information (as linearity and frequency for example). Image reconstruction from wavelet maxima seems an interesting lead, given that significant feature sets are located on amplitude's maxima, and certainly contain a sparse encoding of the neighborhoods. This could be used for denoising. Prospects of this work also include the definition of a full 9D split-of-identity for complete reversibility. Finally, perceptually meaningful color spaces and dependencies between colour channels have to be considered. VII. THANKS R. Soulard works with CRESCEN2O. We want to thank the anonymous reviewers for their valuable comments.

Fig. 3 .

 3 Fig. 3. Wavelets reconstructed from one unique monogenic coefficient with different values of ϕ and θ (A = 1, i = 6).

Fig. 4 .

 4 Fig. 4. Color Fourier atom. Colors of a(t) are generated by centering around middle grey. ã(t) is the rotated version of a(t) by (α 1 , α 2 , α 3 ).

Fig. 5 .

 5 Fig.5. Color 1D signal s(t) and its estimated amplitude A and phase ϕ.

Fig. 8 .

 8 Fig. 8. Example of color monogenic wavelet shrinkage based on amplitude. The noise is Gaussian with σ = 0.1 and hard thresholding is done through 7 scales. Marginal approach: A R , A G and A B are processed separately with T = 0.4. Proposed approach: the elliptical amplitude A is processed with T = 0.58.

α 1 = 3 π 7 α 2 = 0 α 3 = 3 π 7 α 1 = -2 π 7 α 2 Fig. 10 .

 72717210 Fig. 10. Examples of modulation of color monogenic wavelets; synthetized from the Radon based filterbank at the 6th scale. Note that the contrast of the first column has been enhanced.

Fig. 11 .

 11 Fig. 11. Examples of multi-oriented color wavelets. λ = 1, ϕ = π/2, α 1 = 0.13π, α 2 = -0.06π. Values of (θ R , θ G , θ B ) are given under each one.

Fig. 12 .

 12 Fig. 12. Estimation of the local frequency ν(x) in radians/pixel. The noisy version is obtained by adding Gaussian noise with σ = 0.1.

Fig. 13 .

 13 Fig. 13. From left to right: Natural image, estimated frequency ν with 5 × 5 median filtering, histogram of ν.

The exact same analysis is well known about the gradient operator, when studied as the solution of a maximum variation problem[START_REF] Jähne | Digital Image Processing 6th Edition[END_REF].