demodulation of AM–FM images

Appl. Comput. Harmon. Anal. 39 (2015) 450–486

Contents lists available at ScienceDirect

Applied and Computational Harmonic Analysis www.elsevier.com/locate/acha

The monogenic synchrosqueezed wavelet transform: a tool for the decomposition/demodulation of AM–FM images Marianne Clausel ∗ , Thomas Oberlin, Valérie Perrier University of Grenoble–Alpes, Laboratoire Jean Kuntzmann, UMR CNRS 5224, BP 53, 38041 Grenoble cedex 9, France

a r t i c l e

i n f o

Article history: Received 20 November 2012 Received in revised form 26 September 2014 Accepted 16 October 2014 Available online 22 October 2014 Communicated by the Editors MSC: 65T60 92C55 94A08 Keywords: Monogenic signal Wavelet transform Directional time–frequency image analysis Synchrosqueezing

a b s t r a c t The synchrosqueezing method aims at decomposing 1D functions into superpositions of a small number of “Intrinsic Modes”, supposed to be well separated both in time and frequency. Based on the unidimensional wavelet transform and its reconstruction properties, the synchrosqueezing transform provides a powerful representation of multicomponent signals in the time–frequency plane, together with a reconstruction of each mode. In this paper, a bidimensional version of the synchrosqueezing transform is defined, by considering a well-adapted extension of the concept of analytic signal to images: the monogenic signal. We introduce the concept of “Intrinsic Monogenic Mode”, that is the bidimensional counterpart of the notion of Intrinsic Mode. We also investigate the properties of its associated Monogenic Wavelet Decomposition. This leads to a natural bivariate extension of the Synchrosqueezed Wavelet Transform, for decomposing and processing multicomponent images. Numerical tests validate the effectiveness of the method on synthetic and real images. © 2014 Elsevier Inc. All rights reserved.

1. Introduction In the past few years, there has been an increasing interest in representing images using spatially varying sinusoidal waves. As the understanding of the theory advanced, amplitude- and frequency-modulation (AM–FM) decompositions have been applied in a large range of problems: for example motion estimation using a optical-flow method based on an assumption of phase-invariance [1–3], reconstruction of breast cancer images [4,5], texture analysis [6,7], or ultrasound image segmentation methods [8] based on the concept of monogenic signal and quadrature filters [9]. A survey of applications of AM–FM decompositions in medical imaging can also be found in [10].

* Corresponding author. E-mail addresses: [email protected] (M. Clausel), [email protected] (T. Oberlin), [email protected] (V. Perrier). http://dx.doi.org/10.1016/j.acha.2014.10.003 1063-5203/© 2014 Elsevier Inc. All rights reserved.

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

451

In each case, the main challenge is to decompose any input image s(x) into a sum of bidimensional AM–FM harmonics of the form s(x1 , x2 ) =

L 

s (x1 , x2 ) =

=1

L 

  A (x1 , x2 ) cos ϕ (x1 , x2 ) ,

(1)

=1

where A > 0 denotes a slowly-varying amplitude function, ϕ denotes the phase, and  = 1, · · · , L indexes the different AM–FM harmonics. With each phase function, one can associate a local frequency vector field defined as ω = ∇ϕ . Finding the components s from the bidimensional signal s is called the decomposition problem. When the decomposition is given, another subsequent problem is to determine for each component s its amplitude, phase and frequency functions involved in (1). The amplitudes can be related to the energy contained at each point of the image, whereas significant texture variations are captured in the frequency content. For example, in the case of a single component, the local frequency vectors are orthogonal to the isovalue intensity lines of an image, while their magnitudes provide a measure of local frequency content. The problem of estimating the local amplitudes, phases and frequencies of a signal of the form (1) (bidimensional or not) is called the demodulation problem. In the one dimensional case, these two problems have been widely investigated both from a theoretical and practical point of view. To address the demodulation problem in the one dimensional setting, a first and necessary step is to define, in a proper way, the concept of instantaneous amplitude, phase and frequency of a given signal. To this end, one can use the well-known Hilbert transform, defined for any f ∈ L2 (R) as  For a.e. t ∈ R,

Hf (t) = lim

ε→0

1 π

 |t−s|>ε

 f (s) ds . t−s

The analytic signal associated with f is then the complex-valued function F = (1 + iH)f . Thereafter, one defines amplitude A and phase ϕ of f as, respectively, the modulus and argument of the analytic signal F associated with f (the uniqueness of phase being ensured under smoothness assumptions on f and under some initial conditions of the form ϕ(0) = ϕ0 ). One then obtains For a.e. t ∈ R,

F (t) = A(t) eiϕ(t) .

The instantaneous frequency of f is thus the derivative of the phase ϕ [11]. The practical estimation of amplitude, phase and instantaneous frequency was the subject of numerous studies (see [12] for a survey). The “naive” method, which consists in using the Hilbert transform to estimate A and ϕ, is known to be numerically unstable [13]. Alternative methods use the wavelet transform associated with reallocation techniques: in the one dimensional case, the time–frequency localization of wavelets allows to compute robust estimations of instantaneous frequency lines in the time–frequency representation: several approaches have been developed in the 90’s, one should cite the method of “wavelet ridges” (see the two pioneer works [14] and [15]), the squeezing method introduced in [16] or the reassignment method [17]. On the other side, the decomposition problem consists in developing decompositions of any signal into a sum of “well-behaved” AM–FM components, and was widely studied. A recent attempt in this context is Empirical Mode Decomposition (EMD), initially introduced by Huang et al. in [18] and later popularized by Flandrin and his co-authors (see e.g. [19] or [20]). Basically, EMD is the output of an iterative algorithm which provides an adaptive decomposition of any signal into several AM–FM components (see [18,21] for more details). It is now used in a wide range of applications including meteorology, structural stability analysis, and medical studies [12]. In spite of its simplicity and its efficiency, this method is hard to analyze mathematically, since it is defined in an empirical way. To tackle this problem, in [22], the authors propose

452

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

and analyze an alternative method called “synchrosqueezing” derived from reassignment methods: they introduce a class of functions which can be viewed as the superposition of a reasonable number of modes,  i.e. which read  A (t) eiϕ (t) , where the amplitude A (t) of each mode is slowly varying with respect to its instantaneous frequency ω (t) = ϕ (t). For each function belonging to this class, the synchrosqueezing provides an estimation of this (unknown) decomposition. Synchrosqueezing presents many similarities with Empirical Mode Decomposition as shown in [22] (see also [23] for a comparison of these two methods or [24] for more about potential applications). In addition, since synchrosqueezing is based on the wavelet transform, it can be used to address simultaneously the decomposition and demodulation problem as done in [22]. This paper aims to generalize the synchrosqueezing to dimension 2, in order to analyze multicomponent images of the form (1). Up to now, many works have been achieved that analyze such multicomponent images. Notably, several bidimensional extensions of Empirical Mode Decomposition have been proposed in the pioneering works [25–27] and [28]. All these different versions of 2D EMD rely on the concepts of local extrema, envelopes and local mean: but in practice, such concepts are formulated differently, leading to various shapes of the 2D IMF, and finally various forms of the 2D EMD. Another drawback of the EMD is that it suffers from a lack of theoretical background. Recently, alternative methods to bidimensional EMD based on a rigorous mathematical framework have been proposed. Let us for instance mention [29] where the 2D sifting process is replaced by the resolution of a nonlinear PDE. The next step after the image decomposition into IMF, relies to compute the bidimensional characteristics of each component, namely the local amplitudes and frequencies. To define these characteristic features, several extensions of the Hilbert transform may be used, such as the directional Hilbert transform (see [30–32] and references therein), the total Hilbert transform [33], or the Riesz transform [27,34]. Despite its apparent simplicity, the use of directional versions of Hilbert transform leads to quite high computational cost for truly 2D signals, i.e. 2D signals which do not vary with respect to one orientation only. In addition, in the case of intrinsically bidimensional images, this method yields an ambiguous definition of the amplitude and phase, as stated by [30,32]. It is the reason why other multidimensional extensions of the Hilbert transform have been proposed, such as the total Hilbert transform defined in [35] or the Riesz transform used in [30,31,9]. In what follows, our bidimensional extension of the analytic signal will be based on the Riesz transform rather than on the total Hilbert transform, since it leads to easier interpretation. The approach we proposed here is in the same spirit as in [36], in the sense that we also use the synchrosqueezing reassignment step. But, as in [22], we aim to establish theoretical results about the decomposition and the reconstruction of an image into a sum of AM–FM modes from the synchrosqueezed wavelet coefficients. Note that similar theoretical reconstruction results are not achievable in the wave-packet framework of [36], as it does not allow direct reconstruction from synchrosqueezed coefficients. In addition, our concept of mode is quite different of that proposed in [36]. Indeed, in [36], the components are complex modes of the type A(x) eiϕ(x) : to our knowledge, there is no unique way of writing a bidimensional real function f in the form f (x) = Re(A(x) eiϕ(x) ). On the opposite, our monogenic synchrosqueezing equally deals with real AM–FM images or with the so-called associated Intrinsic Monogenic Mode Functions, defined in a unique way using the Riesz transform. This is similar to the one-dimension case, where the analytic wavelet indifferently processes complex modes A(t) eiϕ(t) or their associated real part. Recall now that if we are given a real valued function f ∈ L2 (R2 ), one defines its Riesz transform Rf as the following vector-valued function  Rf = where for any i = 1, 2 and a.e. x ∈ R2 one has

R1 f R2 f

 ,

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486



1 Ri f (x) = lim ε→0 π

 |x−y|>ε

453

 (xi − yi ) f (y) dy . |x − y|3

Then, the monogenic signal associated with f is  Mf =

f Rf

 .

It can also be read Mf = f + R1 f i + R2 f j, where (1, i, j, k) denotes the canonical basis of the algebra H of the quaternions (see Appendix A for some background about quaternionic calculus and Appendix B for the main properties of the Riesz transform). The monogenic signal is the bidimensional counterpart of the analytic signal, when replacing the Hilbert transform with the Riesz transform. It can be described using three main characteristics: local amplitude, local phase, and local orientation which can be defined by means of the polar form of the monogenic signal (see [9] and [37]) Mf = A eϕnθ = A(cos ϕ + nθ sin ϕ)

with nθ = cos θ i + sin θ j,

where A is a positive real valued function and ϕ and θ are two real valued functions defined in a unique way (under suitable assumptions). The function A is called the (local) amplitude of f , ϕ its (local) phase and nθ its (local) orientation. The function ω = ∇ϕ will be called the local frequency of f . To illustrate these definitions, we consider the case where f (x) = A0 cos(k · x) with k = (k1 , k2 ) and k1 > 0. Set θ0 = Arctan(k2 /k1 ). One then has  Rf (x1 , x2 ) = A0

sin(k1 x1 + k2 x2 ) cos θ0 sin(k1 x1 + k2 x2 ) sin θ0

 = A0

k sin(k · x), |k|

and  Mf (x) =

f (x) Rf (x)

 = A0 e(k·x)(cos θ0 i+sin θ0 j) .

Then for all x ∈ R2 , one can set A(x) = A0 ,

ϕ(x) = k · x,

θ(x) = θ0 = Arctan(k2 /k1 ).

The aim of the present paper is to define the synchrosqueezing of bidimensional images, in the monogenic setting, in order to address both the decomposition and the demodulation problems. To this end, we will consider monogenic signals which are superposition of monogenic waves of the form A(x) eϕ(x)nθ(x) ,

(2)

where A and θ are slowly varying functions with respect to ϕ. We will consider successively the two different problems of deconvolution and decomposition. As in [22], our method is based on the continuous wavelet transform. Therefore, in Section 2, we recall a recent extension of the continuous bidimensional wavelet transform to the monogenic setting. We then investigate in Section 3, the problem of estimating the local amplitude and frequency of a monogenic spectral line of the

454

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

form (2). Subsequently, we define in Section 4 the bidimensional synchrosqueezing and its application to the decomposition of multicomponent images. Finally, Section 5 gives some insights on the practical implementation of bidimensional synchrosqueezing, and provides numerical experiments. Proofs are postponed to Section 6 to make reading easier. 2. 2D-Wavelet transform and monogenic wavelet transform In the one dimensional setting, the synchrosqueezing method is based on the one-dimensional wavelet analysis: indeed, the wavelet transform is used both to decompose an analytic signal into several components, and also to estimate the local amplitude and frequency of each component. The main advantage of using the wavelet analysis is to avoid the use of the discrete Hilbert transform, since there is a link between the wavelet coefficients of the analytic signal associated with a real function f , and its analytic wavelet coefficients (i.e. the wavelet coefficients of f among an analytic wavelet). In the bidimensional case, it is then quite natural to wonder if this approach can be extended using bidimensional wavelet analysis, replacing the analytic signal by the monogenic signal associated with our data. Since the monogenic signal can be viewed as a 3D vector field, it can be componentwise analyzed by the usual wavelet transform: Section 2.1 recalls first some basics about the bidimensional continuous wavelet transform. Thereafter in Section 2.2, we recall the main characteristics of the monogenic wavelet analysis, which is an extension of the usual wavelet analysis introduced in [38] and [39]. We would like to point out that, for any bidimensional real image f , the wavelet coefficients of the monogenic signal Mf can be related to the monogenic wavelet coefficients of f (see Proposition 2.1 for a precise statement). Then, using monogenic wavelet analysis, one can recover the characteristics of the monogenic signal associated with our data, without estimating it. In what follows, for any (a, b, α) ∈ R∗+ × R × (0, 2π), we denote by Da the dilation operator, Tb the translation operator and Rα the rotation operator defined on L2 (R2 ) by: Da f (x) = a−1 f (x/a),

  Rα f (x) = f rα−1 x ,

Tb f (x) = f (x − b),

where rα is the usual 2 × 2 rotation matrix of angle α:  rα =

cos α sin α

− sin α cos α

 .

(3)

2.1. The usual bidimensional wavelet transform We briefly recall classical definitions about the bidimensional continuous wavelet transform. A (real or complex) function ψ ∈ L2 (R2 ) is called an admissible wavelet if it satisfies1 : 

2 |ψ(ξ)| dξ < +∞ |ξ|2

2

Cψ = (2π)

R2

(4)

As usual (see [40]), the wavelet family {ψa,α,b }(a,α,b)∈R∗+ ×(0,2π)×R2 is defined by ψa,α,b = Tb Rα Da ψ. One then defines the wavelet coefficients of any f ∈ L2 (R2 ) as follows:  cf (a, α, b) =

f (x)ψa,α,b (x) dx. R2

1

The Fourier Transform of ψ being defined by: ψ(ξ) =

1 2π

R2

ψ(x) e−iξ·x dx

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

455

If the wavelet is assumed to be isotropic, the wavelet coefficients do not depend on α. In this case, we will denote cf (a, b) = cf (a, 0, b) and ψa,b = ψa,0,b . Any square integrable function f can be recovered from its wavelet coefficients using the so-called reconstruction formula:    1 da f (x) = cf (a, α, b)ψa,α,b (x) 3 dα db, (5) Cψ a b∈R2 α∈(0,2π) a∈(0,+∞)

where this equality stands in L2 (R2 ). When the wavelet is assumed to be isotropic, we get a simpler expression: f (x) =

2π Cψ



 cf (a, b)ψa,b (x)

da db. a3

b∈R2 a∈(0,+∞)

Other reconstruction formulas are available (see e.g. [40]), among which the pointwise reconstruction formula, obtained by summing over scales only. In the isotropic case, it reads: 2π f (x) = C˜ψ

+∞  da cf (a, x) 2 a

with C˜ψ =

 ψ(ξ) dξ. |ξ|2

(6)

R2

0

As in the one dimensional case [22], this formula will play a special role in the bidimensional synchrosqueezing. For any real function f ∈ L2 (R2 ), one can also define the wavelet transform of its monogenic signal F = Mf = f + R1 f i + R2 f j, as: ⎛

cf

⎞

cF = cf + cR1 f i + cR2 f j = ⎝ cR1 f ⎠ . cR2 f 2.2. The monogenic wavelet transform We now present the continuous monogenic wavelet transform as defined in [38]. We consider a real admissible wavelet ψ and define ⎛

ψ (M )

⎞ ψ = Mψ = ⎝ R1 ψ ⎠ R2 ψ

the associated monogenic wavelet. Observe that ψ (M ) is also an admissible wavelet (taking values in R3 ), since each of its components satisfies relation (4). We define below the monogenic wavelet coefficients of a real function. (M )

Definition 2.1. Let f ∈ L2 (R2 ). The monogenic wavelet coefficients cf (M ) cf (a, α, b)

 =

(a, α, b) of f are defined by:

(M )

f (x)ψa,α,b (x) dx, R2 (M )

where for any (a, α, b) ∈ (0, +∞) × (0, 2π) × R2 , ψa,α,b = Tb Rα Da (Mψ).

(7)

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

456

(M )

Remark 2.1. For any (a, α, b), cf cf (a, α, b) +

(1) cf (a, α, b) i

+

(M )

(a, α, b) is a Clifford vector, and thus can be written as cf

(2) cf (a, α, b) j

(a, α, b) =

where we denote for i = 1, 2

(i) cf (a, α, b)

 f (x) (Tb Rα Da Ri ψ)(x) dx

= R2

(see Appendix A for more details). The next proposition states that there exists an explicit relationship between the monogenic coefficients of a real image f and the usual wavelet coefficients of its monogenic signal Mf . Proposition 2.1. Let f ∈ L2 (R2 ) and denote by F = Mf the monogenic signal associated with f . Then for any (a, α, b) ∈ (0, +∞) × (0, 2π) × R2  cF (a, α, b) =

1 0

0 −rα

 (M )

cf

(a, α, b).

(8)

Proof. We will prove that  (M )

cf

(a, α, b) =

1 0

0 −rα−1

 cF (a, α, b),

which is equivalent to (8). By definition of the monogenic wavelet transform one has: (M )

cf

 (a, α, b) =

f (x) Tb Rα Da (Mψ)(x) dx. R2

The translation-, scale-invariance and steerability properties of the Riesz transform recalled in Propositions B.3 and B.4 of Appendix B, imply that: Tb Rα Da (Rψ) = rα−1 R(ψa,α,b ). Hence 

f (x) Tb Rα Da (Rψ)(x) dx = rα−1

R2

 f (x) R(ψa,α,b )(x) dx.

(9)

R2

Since by Proposition B.5, the Riesz transform is a componentwise antisymmetric operator on L2 (R2 ), one deduces that:     f, R1 (ψa,α,b )L2 (R2 ) = − (Rf )(x)ψa,α,b (x) dx. f (x)R(ψa,α,b )(x) dx = (10) f, R2 (ψa,α,b )L2 (R2 ) R2

Eqs. (9) and (10) then directly imply relation (8).

R2

2

We now consider the case an admissible isotropic real-valued wavelet ψ. Examples include the Mexican hat (Laplacian of Gaussian) or the Morse wavelets (see [41] and Section 5). In such a case, the wavelets ψa,α,b as well as the wavelet coefficients cF (a, α, b) of the monogenic signal do not depend on the orientation α. From now, we then denote ψa,α,b and cF (a, α, b) respectively by ψa,b and cF (a, b). Thus Eq. (8) can be simplified in the following way:

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

 cF (a, b) =

1 0 0 −rα



 (M )

cf

(a, α, b) =

⎛

1 0 0 −rα



⎛

⎞

⎜ c(1) (a, α, b) ⎟ ⎝ f ⎠ (2)

⎞

cf (a, b)

cf (a, b)

457

cf (a, α, b)

(1) (2) ⎟ ⎜ = ⎝ − cos αcf (a, α, b) + sin αcf (a, α, b) ⎠ . (1) − sin αcf (a, α, b)

−

(11)

(2) cos αcf (a, α, b)

In practical situations, it will be easier to consider the case in which α = 0. We then get that  cF (a, b) =

1 0 0 −Id

 (M )

cf

(a, 0, b).

(12)

To sum up, in the isotropic case, the monogenic wavelet coefficients of the 2D real image f can be related in a very simple way to the wavelet coefficients of the monogenic signal using (12). It will be useful in the sequel since when analyzing images, one deals with a real-valued signal f and aims at recovering the wavelet coefficients cF (a, b) of its monogenic signal F without computing F . 3. Wavelet analysis of bidimensional spectral lines In this section, we tackle the problem of estimating the local amplitude and frequency of a bidimensional mode of the form A(x) eϕ(x)n(x) under slow-variation assumptions on functions A, ϕ and nθ . Note that the following reasoning can be adapted to the particular case of real modes of the form A(x) cos(ϕ(x)). This is similar to the dimension one, where the analytic wavelet transform considers indifferently complex modes eiϕ(x) or their real parts. In the one dimensional setting, the “wavelet-ridge” method has proved to be efficient to address this problem [14,15]. In [38,41] it has been extended to the bidimensional setting using monogenic wavelet analysis. Here we follow an alternative approach consisting in extending synchrosqueezing to the bidimensional context. We first define precisely the bidimensional modes that we are considering, and that we will call “Intrinsic Monogenic Mode Function” (IMMF). Then we deduce the notion of bidimensional frequency vector and state an approximation result using the monogenic wavelet transform. −σ Before any statement we introduce the functional space B˙ ∞,1 (R2 , H): −σ Definition 3.1. Let σ > 0. The homogeneous Besov space B˙ ∞,1 (R2 , H) is the collection of tempered distributions F such that:

∞ −σ = F B˙ ∞,1

0

 da  sup cF (a, b) 2−σ < +∞. a b∈R2

(13)

Remark 3.1. Alternative definitions based on finite differences can be found in [42]. We now define our concept of bidimensional modes: Definition 3.2. Let ε > 0. An Intrinsic Monogenic Mode Function (IMMF) with accuracy ε and smoothness −σ σ > 0 is a function F ∈ B˙ ∞,1 (R2 , H) of the form F (x) = A(x) eϕ(x)nθ(x)

    with nθ(x) = cos θ(x) i + sin θ(x) j,

(14)

458

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

where A > 0 and A, θ ∈ C 1 (R2 ) ∩ L∞ (R2 ), ϕ ∈ C 2 (R2 ). The function A is called the amplitude of F , whereas ϕ and nθ are called respectively the scalar phase and the local orientation of F . The function ϕ is assumed to have bounded derivatives of order 1 and 2: ∀i ∈ {1, 2},

    0 < inf2 ∂xi ϕ(x) ≤ sup ∂xi ϕ(x) < ∞.

(15)

M = max

  sup ∂x2i1 xi2 ϕ(x) < ∞.

(16)

x∈R

x∈R2

i1 ,i2 =1,2 x∈R2

Further we assume that A, θ, ∇2 ϕ are slowly varying functions with respect to ∇ϕ, namely that for i = 1, 2 ∀x ∈ R2 ,

      max ∂xi A(x), ∂xi θ(x) ≤ ε∂xi ϕ(x)

(17)

and ∀x ∈ R2 ,

    max ∂x2i1 xi2 ϕ(x) ≤ ε∇ϕ(x).

i1 ,i2 =1,2

(18)

Remark 3.2. The approximation result that we will establish in the following, will also hold when considering spectral lines of the form f (x) = A(x) cos(ϕ(x)) as done in [41], provided that functions A, ϕ and n = ∇ϕ/|∇ϕ| satisfy assumptions of Definition 3.2. Thus, the results stated hereafter for IMMFs can be transposed to the case of real modes. Remark 3.3. Observe that, in the sequel, we shall consider non-square integrable functions. In this case, the two reconstruction formulas (5) and (6) may not hold (see Appendix B of [43]). We then need an additive

+∞ assumption to ensure that the integral 0 cF (a, b) da a2 exists to give some sense to the reconstruction formula (6). Moreover, we will estimate the behavior of +∞    cF (a, b) da a2 a0 −σ for any a0 > 0. The assumption F ∈ B˙ ∞,1 (R2 , H) implies that for any a0 > 0:

∞ ∞  da  1 da     −σ . cF (a, b) 2 = sup cF (a, b) σ 2−σ ≤ a−σ sup 0 F B˙ ∞,1 a a a b∈R2 b∈R2 a0

(19)

a0

Moreover, by hypothesis, F ∈ C 1 (R2 , H) and ∇F is bounded. Using the expression of the wavelet coefficients then implies that:   cF (a, b) ≤ a2 ∇F L∞ xψL1 ,

∀b ∈ R2

which ensures that for any a0 > 0 a0  da  sup cF (a, b) 2 < ∞ 2 a b∈R 0

and finally leads to the existence of the integral

+∞ 0

cF (a, b) da a2 .

(20)

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

459

We now focus on the estimation of the local amplitude and frequency of an IMMF. To this end, we need some assumptions on the wavelet that we list below: Assumptions (W) (i) the wavelet ψ is isotropic and real-valued. (ii) The moments of order 0, 1, 2, 3 of ψ and ∇ψ are finite, namely Iα < ∞,

sup α∈{0,1,2,3}

Iα < ∞,

sup α∈{0,1,2,3}

where for any α > 0  Iα =



 |x| ψ(x) dx, α

R2

Iα

 =

  |x|α ∇ψ(x) dx.

(21)

R2

We first deal with the simple case of an IMMF of constant amplitude and phase: Lemma 3.1. Let F be an IMMF of the form F (x) = A0 e(k·x)nθ with k = (k1 , k2 ) ∈ R2 , A0 ∈ R∗+ , θ ∈ R, and let ψ be a wavelet satisfying Assumptions (W). The wavelet transform of F is given by:     cF (a, b) = A0 aψ(ak) cos(k · b) + sin(k · b)(cos θ i + sin θ j) = aψ(ak) A0 e(k·b)nθ

(22)

and one has for i = 1, 2:    ∂bi cF (a, b) = ki nθ aψ(ak) A0 e(k·b)nθ . The vectors k and nθ can then be recovered using the two following relations:  −1 k1 nθ = ∂b1 cF (a, b) × cF (a, b) ,  −1 . k2 nθ = ∂b2 cF (a, b) × cF (a, b) Proof. By definition, since the wavelet is isotropic real, one has: cF (a, b) = a−1



 F (x)ψ R2

  x−b dx = a F (au + b)ψ(u)du a R2

A simple calculation leads to:

⎛ ⎞ ⎞ cos(k.b) cos(k · (au + b))ψ(u)du ⎝ sin(k.b) cos(θ) ⎠ = aψ(ak)F cF (a, b) = aA0 ⎝ cos θ R2 sin(k · (au + b))ψ(u)du ⎠ = A0 aψ(ak) (b),

sin θ R2 sin(k · (au + b))ψ(u)du sin(k.b) sin(θ) ⎛

R2

where we have used ψ(−ak) = ψ(ak) since the wavelet is isotropic. Now, since ψ satisfies Assumptions (W), the moment I0 is finite and F ∈ L∞ (R2 ), we get: ∂bi cF (a, b) = a−1

 R2

      x−b x−b F (x)∂bi ψ dx = −a−2 F (x)(∂xi ψ) dx. a a R2

(23)

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

460

A similar calculation, using ∂ xi ψ(ξ) = iξi ψ(ξ), ∀ξ = (ξ1 , ξ2 ), and ψ isotropic, yields: ⎛

⎞ − sin(k.b) ⎝ cos(k.b) cos(θ) ⎠ = aki ψ(ak)n ∂bi cF (a, b) = A0 aki ψ(ak) θ F (b), cos(k.b) sin(θ) where we have used:   nθ F (b) = nθ cos(k · b) + nθ sin(k · b) = − sin(k · b) + nθ cos(k · b). This leads to (23).

2

Remark that, more generally, if F (x) = A0 e(k·x+α)nθ (α ∈ R): cF (a, b) = A0 aψ(ak) e(k·b+α)nθ = aψ(ak)F (b)

(24)

∂bi cF (a, b) = aki ψ(ak) nθ F (b)

(25)

The formulas (22) and (23) can be extended in the general case of an IMMF F of the form (14). Proposition 3.1. Let ψ a wavelet satisfying Assumptions (W) and F an IMMF with accuracy ε and smoothness σ > 0 of the form (14). For all (a, b) ∈ R∗+ × R2 , one has (i)    cF (a, b) = aψ a∇ϕ(b) A(b) eϕ(b)n(b) + εa2 R1 (a, b), with   √   √    R1 (a, b) ≤ I1 2A(b) + 1 ∇ϕ(b) + aI2 A(b) ∇ϕ(b) + 2M /2 + aI2 M/2 + a2 I3 A(b)M/6. (26) (ii)     ∂bi cF (a, b) = ∂bi ϕ(b)nθ(b) aψ a∇ϕ(b) A(b) eϕ(b)nθ(b) + εaR2 (a, b), with        √   R2 (a, b) ≤ A(b) a∇ϕ(b)I2 /2 + a2 M I3 /6 + 2A(b) + 1 ∇ϕ(b)I1 + aM I2 /2 .

(27)

(M )

Remark 3.4. As a consequence of this proposition, since in the isotropic case cF (a, b) and cf (a, α, b) are related by (11), the above proposition specifies the results stated in [38] and [41]. Notably, it gives a bound of the approximation error. Proof. See Section 6.1.

2

We now define the two following Clifford vectors:  −1 Λ1 (a, b) = ∂b1 cF (a, b) × cF (a, b) ,  −1 . Λ2 (a, b) = ∂b2 cF (a, b) × cF (a, b)

(28)

These two vectors will provide an approximation of the vectors ∂bi ϕ(b)nθ(b) for an IMMF F of the form (14). More precisely, we can state the following result:

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

461

Theorem 3.2. Let F be an IMMF with accuracy ε > 0 and smoothness σ > 0 of the form (14). Assume that we are given a wavelet ψ satisfying Assumptions (W) and consider (a, b) ∈ R∗+ × R2 such that   cF (a, b) ≥ εν

for some ν ∈ (0, 1/2).

Then for i = 1, 2:        Λi (a, b) − ∂b ϕ(b)nθ(b)  ≤ ε1−ν a R2 (a, b) + aR1 (a, b)∇ϕ(b) , i

(29)

where R1 , R2 have been defined respectively in (26) and (27). In addition for any a0 > 0, there exists some ε0 > 0 such that, for all (a, b) ∈ (0, a0 ) × R2 and any 0 < ε ≤ ε0 :   Λi (a, b) − ∂b ϕ(b)nθ(b)  ≤ εν . i Proof. See Section 6.2.

(30)

2

Since cos(ϕ(x)) = cos(−ϕ(x)), there is an ambiguity in the definition of the local frequency ∇ϕ since ϕ can be replaced by −ϕ. To solve this problem, practitioners usually assume that the first component of the local frequency is positive (see [10] and references therein for more details). Under this additional assumption, our method yields an estimate of the local amplitude and frequency of an IMMF: Proposition 3.3. The notations and assumptions are those of Theorem 3.2. In addition, we assume that: ∀x ∈ R2 ,

∂x1 ϕ(x) > 0.

Then for i = 1, 2, and (a, b) ∈ (0, a0 ) × R:    ∂b ϕ(b) − Λ1 (a, b) ≤ εν 1

(31)

      ∂b ϕ(b) − Λ2 (a, b) sgn Re ∂b cF (a, b)∂b cF (a, b)  ≤ εν . 2 1 2

(32)

and

Proof. See Section 6.3.

2

4. Decomposition of a multicomponent image into spectral lines In this section we will consider a more general context, where our model is now a superposition of several IMMFs (introduced in Section 3), assumed to be well separated, as defined in the following Definition 4.1. We will show in Theorem 4.1 how such functions can be decomposed into monogenic modes, using synchrosqueezing. Thereafter, for each component, the associated local amplitude and frequency will be estimated using the results of Section 3. Definition 4.1. A function F defined from R2 to H is said to be a superposition of well separated Intrinsic Monogenic Mode Components with smoothness σ > 0 up to accuracy ε > 0 and with separation d > 0, if there exists a finite integer L such that

F (x) =

L  =1

F (x),

(33)

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

462

where all the F are IMMFs with accuracy ε and smoothness σ ≥ σ > 0 of the form (14): F (x) = A (x) eϕ (x)nθ (x) , and moreover satisfy for any x,  >  and i ∈ {1, 2}:     ∂x ϕ (x) > ∂x ϕ (x), i i

(34)

      ∂x ϕ (x)nθ (x) − ∂x ϕ (x)nθ  (x)  ≥ d ∂x ϕ (x) + ∂x ϕ (x) . i i i i  

(35)

and

We denote by Aε,σ,d the class of all functions F satisfying these conditions. We now define the monogenic synchrosqueezing transform (MSST) of a function F belonging to Aε,σ,d . Definition 4.2. Assume that we are given ψ a wavelet satisfying Assumptions (W) such that   ⊂ ξ ∈ R2 ; 1 − Δ ≤ |ξ| ≤ 1 + Δ , supp(ψ) with Δ<

d . 2(1 + d)

(36)

Consider h ∈ Cc∞ (R, R) such that R h(x)dx = 1. Let δ > 0, ν ∈ (0, 1), and ε > 0. For F ∈ Aε,σ,d , the monogenic synchrosqueezing transform of F is defined for any (b, k, n) ∈ R2 × R2 × S1 (S1 being the unit sphere of R2 ) by  δ,ν SF,ε (b, k, n)

=

    1 k2 − Re(nΛ2 (a, b)) da k1 − Re(n Λ1 (a, b)) h cF (a, b) 2 h , δ δ δ a2

(37)

Aε,F (b)

with     Aε,F (b) = a ∈ R+ ; cF (a, b) > εν , and where Λ1 (a, b), Λ2 (a, b) are defined by Eq. (28). We can now state our main result: Theorem 4.1. The notations are those of Definition 4.2. Let F ∈ Aε,σ,d and ν ∈ (0, 1/(2 + 4/σ)). Provided that ε > 0 is sufficiently small, the following assertions hold: • |cF (a, b)| > εν only if there exists some  such that       (a, b) ∈ Z = (a, b), a∇ϕ (b) − 1 < Δ .

(38)

• For each  = {1, · · · , L} and each pair (a, b) ∈ Z for which |cF (a, b)| > εν , one has for i = 1, 2 and for some C > 0:  Λi (a, b) − ∂i ϕ (b)nθ

 (b)

  ≤ Cεν .

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

463

• Moreover for any  ∈ {1, · · · , L}, there exists C > 0, such that for any b ∈ R2   lim 

δ→0

1 2π C˜ψ



  δ Sf,ε (b, k, n)dkdn − A (b) eϕ (b)nθ (b)  ≤ Cεν ,



S1 B (εν ,n,b)

where we denote B (εν , n, b) = {k ∈ R2 ; maxi |ki n − ∂bi ϕ (b)nθ (b) | ≤ εν }, and where C˜ψ was introduced in (6). Remark 4.1. Note that the wave-packet synchrosqueezing of [36] does not allow direct reconstruction from the synchrosqueezed coefficients, that prevents from stating an approximation result similar to Theorem 4.1. Nevertheless, in our approach an additional assumption on the separation of modes in frequency is required, to ensure the uniqueness of the decomposition. Contrary to [36], we thus do not consider two waves with same instantaneous frequency magnitude and different orientations. Proof. See Section 6.4.

2

5. Implementation and numerical experiments This section aims to illustrate the efficiency of the monogenic synchrosqueezing transform, and its potential applications. We will first describe the computation of the discrete MSST from a N × N discrete image, then we will show several examples on artificial and real images. The Matlab code used for all the experiments and figures of this paper can be downloaded from http://www-ljk.imag.fr/membres/Thomas. Oberlin/MSST.tar.gz. 5.1. The discrete monogenic synchrosqueezing transform The first step consists in the discretization of the variables, therefore one builds discrete sets A, K, O for the scales a, the normalized frequencies (k1 , k2 ), and the orientations θ respectively, the frequency resolution in k being denoted by Δk . Note that one usually chooses a logarithmic scale for a and k, e.g. A = {2j/nv }j=0···na −1 , so that we have nv coefficients per octave. Then, we compute the discrete monogenic wavelet transform cF (a, b), and the estimations of the local frequencies Λ1 (a, b) and Λ2 (a, b) defined in Eq. (28) by writing for i = 1, 2:  ∂bi cf (a, b) =

f (x)∂bi ψa,b (x) dx, R2

the other components ∂bi cR1 f (a, b) and ∂bi cR2 f (a, b) being computed in the same manner. The monogenic synchrosqueezing transform consists in a partial reallocation of the monogenic wavelet transform according to space, frequency and orientation parameters. Given a threshold γ, and for all a ∈ A, (k1 , k2 ) ∈ K2 , θ ∈ O, one simply writes SF,γ (b, k1 , k2 , nθ ) =

log(2) nv

 ⎧ |cF (a,b)|>γ ⎪ ⎪ ⎨ Δk a∈A s.t. |k1 −Re(Λ1 (a,b)nθ )|≤ 2 1 ⎪ ⎪ Δk ⎩ |k2 −Re(Λ2 (a,b)nθ )|≤ 2 2

cF (a, b) . a

(39)

da The term log(2) anv comes from the measure a2 in Eq. (37), as we deal with logarithmic scales. Note that if we use a logarithmic scale for the frequencies like in [22], the resolution Δk depends on the value k.

464

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

This discrete MSST is relatively easy and cheap to compute, but it is not of great interest in practice because of the large number of variables which prevents both easy representation and fast processing (i.e. decomposition, demodulation). To decrease this number of variables, one can remark that if a superposition of IMMFs satisfies the separation condition of Eq. (34), then for any  >  it also satisfies |∇ϕ (x)| > |∇ϕ (x)|, so that the IMMFs are also separated in terms of the norm of the frequency vector. This leads us to define the following discrete isotropic MSST for all a ∈ A and k ∈ K: SF,γ (b, k) =

log(2) nv

  |cF (a,b)|>γ  a∈A s.t. Δ |k− |Λ1 (a,b)|2 +|Λ2 (a,b)|2 |≤ 2k

cF (a, b) . a

(40)

Intuitively, SF,γ (b, kp ) contains all the coefficients cF (a, b) whose “local isotropic frequency”  |Λ1 (a, b)|2 + |Λ2 (a, b)|2 is around the value k. The latter isotropic MSST can be particularly useful when visualizing the time–frequency distribution, as it only depends on three scalar parameters (b1 , b2 , k), that is by far easier to represent. Hence, this version will be used in the experiments hereafter. Additionally we will use the isotropic Morlet wavelet with parameters μ > 0 and σ > 0, defined in the Fourier domain as follows   2  ψˆμ,σ (ξ) = exp −π 2 σ |ξ| − μ .

(41)

5.2. Representing a synthetic 3-components signal This section illustrates the interest in using the MSST to represent multicomponent signals in 2 dimensions. Let us first define the synthetic 3-components test-signal f = f1 + f2 + f3 with ⎧ 2 2 ⎨ f1 (x, y) = e−10((x−0.5) +(y−0.5) ) sin(10π(x2 + y 2 + 2(x + 0.2y))), f (x, y) = 1.2 sin(40π(x + y)), ⎩ 2 f3 (x, y) = cos(2π(70x + 20x2 + 50y − 20y 2 − 41xy)).

(42)

We compute both its monogenic isotropic wavelet transform cF and its isotropic MSST SF . We then aim at representing |cF | (resp. |SF |), which depends on the three scalar parameters x, y and a (resp. k). We will here propose two distinct visualizations: either a 3-dimensional scatter “density” plot, or a 2-dimensional representation for x or y fixed, which looks similar to a 1D continuous wavelet representation. The following Fig. 1 shows these different representations for the monogenic wavelet transform and the MSST of signal (42), illustrating how the MSST sharpens the space-scale representation. Indeed, the volumetric visualization (Fig. 1(g)) is almost perfectly concentrated on three surfaces, corresponding to the three modes of the signal. The 2D representations (Fig. 1(f) and (h)) also shows that the MSST sharpens the representation, in the same manner as in 1D the reassignment or synchrosqueezing transforms. 5.3. Decomposition and demodulation Once the MSST has been computed, the key point to separate and reconstruct the modes is to identify the local (anisotropic) frequency of each mode |∇ϕ (b)| at each point b. The th IMMF F is then approximated by F  (b) as follows: F  (b) =

 ˆ (b)−κ≤k≤k ˆ (b)+κ k

SF,γ (b, k),

(43)

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

465

Fig. 1. (a): The 2D test-signal represented as an image. (b, c, d): Its components. (e, f): The monogenic wavelet transform |cF | represented as a density in 3D, and in 2D for y = 0.5 fixed. (g, h): The same visualizations of the MSST |SF |.

where kˆ (b) is an appropriate estimate of |∇ϕ (b)|. Obtaining this estimate is an involved problem which will not be discussed here, so we simply use the approach used in [22] based on a greedy algorithm. Our approximation F  (b) of the th IMMF is computed by summing the coefficients of the synchrosqueezed transform in the vicinity of this local frequency, to regularize the solution. In practice the width κ is set to 7, and parameter nv equals 16. To evaluate the method, we compare the mode F of the test-signal defined in Eq. (42) to the corresponding mode Fˆ obtained from the discrete MSST after extraction, by computing the normalized Mean-Squared Error (MSE):

MSE(Fˆ ) =

Fˆ − F  . F 

(44)

466

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

Fig. 2. The extracted and reconstructed modes 1, 2 and 3 for the test signal of Eq. (42). The normalized MSE are respectively 0.03, 0.06 and 0.05. We used the isotropic Morlet Wavelet ψμ,σ defined in (41) with σ = 2 and μ = 1, and removed the border ( 81 at each side for each dimension) to avoid border effects.

Fig. 3. Top: the first three modes given by EMD [28] for the test signal of Eq. (42). The corresponding MSE are 4.1, 0.99 and 1.0. Bottom: the decomposition given by EEMD with MSE 0.59, 0.74 and 0.86.

Each reconstructed IMMF is displayed on Fig. 2, where one also provides the corresponding MSE. It is clear that the modes are reconstructed with very high accuracy, although our method for ridge extraction is quite simple. To compare with, we also display on Fig. 3 the first three components computed by standard EMD, using the implementation of A. Linderhed available at http://www.aquador.vovve.net/IEMD/. We also show the decomposition provided by 2D Ensemble EMD (EEMD) [44]. Fig. 3 shows that both EMD and EEMD fail to extract the modes with high accuracy. In the EMD decomposition, the first component seems to contain f1 + f2 , whereas the second mode contains alternatively f3 and a piece of f2 : this is the well-known mode-mixing phenomenon. The decomposition provided by EEMD seems better, but the modes 2 and 3 are not separated properly. Indeed, EEMD is well-known to decrease the mode-mixing effect, but it does not improve the frequency resolution, i.e. the ability to separate two modes with close frequency content. Then, we aim at showing that the MSST remains efficient on real images or textures, even though they are not truly multicomponent. The following Fig. 4 shows the image Lena where we artificially added a fingerprint image. After a 2D synchrosqueezing transform one is able to extract the oscillating component to roughly reconstruct the fingerprint, in a similar way than EMD or EEMD. This suggests that the MSST is a compact and efficient representation for oriented textures, that should lead in the future to interesting applications in texture synthesis or segmentation. Finally, let us briefly discuss the issue of noise, which is of particular interest in practical situations. As the MSST is based on the wavelet transform, it is quite robust to noise compared to direct tools like EMD.

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

467

Fig. 4. Illustration of an application on real images. (a): Input image, superposition of Lena and a fingerprint. (b): The extracted fingerprint by the 2D MSST. The parameters are the same as for Fig. 2. (c): The first mode given by EMD. (d): The first mode computed by EEMD.

Fig. 5. Illustration of the robustness to noise. (a): Second component of Fig. 1 with white Gaussian noise (MSE = 0.14). (b): The amplitude of its MSST. (c): The reconstructed component (MSE = 0.6).

The situation is similar to the one dimensional case, where the SST remains stable for reasonable noise levels [45,46]. The following Fig. 5 shows the behavior of the MSST on a noisy version of the mode from Fig. 1. We can see the noise in the volumetric representation at high frequencies, but the ridge (the surface) associated with the signal remains clear, allowing us to reconstruct the denoised mode in Fig. 5(c) with a good accuracy. 6. Proofs 6.1. Proof of Proposition 3.1 We prove Proposition 3.1 into several steps. 6.1.1. A preliminary result We first need a preliminary result: Lemma 6.1. Assume that F is an IMMF with accuracy ε of the form (14). The following estimates hold for any y, h ∈ R2 :       A(y + h) − A(y) ≤ ε |h|∇ϕ(y) + |h|2 M , 2       cos θ(y + h) − cos θ(y) ≤ ε |h|∇ϕ(y) + |h|2 M , 2       sin θ(y + h) − sin θ(y) ≤ ε |h|∇ϕ(y) + |h|2 M , 2

(45) (46) (47)

468

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

      ∇ϕ(y + h) − ∇ϕ(y) ≤ ε |h|∇ϕ(y) + |h|2 M , 2

(48)

where M is the constant defined in (16). Proof. We only prove inequality (45), the others follows in the same way. Observe that: 1 A(y + h) − A(y) =

∇A(y + th) · h dt 0

By assumption (17), ∇A is slowly varying with respect to ∇ϕ, then:   A(y + h) − A(y) ≤ ε|h|

1

  ∇ϕ(y + th)dt

0

Applying the usual mean value theorem, we deduce that for any t:       ∇ϕ(y + th) ≤ ∇ϕ(y) + ∇ϕ(y + th) − ∇ϕ(y)      ≤ ∇ϕ(y) + t|h| max sup ∂x2i1 ,xi2 ϕ(y) i1 ,i2

y∈R2

  ≤ ∇ϕ(y) + M t|h|,

applying assumption (16) on second order partial derivatives of ϕ. To sum up, one has   A(y + h) − A(y) ≤ ε|h|

1

  ∇ϕ(y) dt + εM |h|2

0

1

    |h|2   , tdt ≤ ε |h| ∇ϕ(y) + M 2

0

which is the required result. 2 6.1.2. Proof of Proposition 3.1 We now prove Proposition 3.1. Proof of point (i) of Proposition 3.1. By definition of the wavelet coefficients of F , one has  cF (a, b) =

A(x) e

ϕ(x)nθ(x) −1

a

R2



x−b ψ a

 dx.

We now split this sum into three terms −1

cF (a, b) = a



 A(b)

e

R2 −1

+a

+ a−1

ϕ(x)nθ(b)



A(b)  R2

R2

x−b ψ a

 dx

   ϕ(x)n  x−b ϕ(x)nθ(b) θ(x) dx e −e ψ a

    x−b dx. A(x) − A(b) eϕ(x)nθ(x) ψ a

(49)

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

469

Observe now that 1 ϕ(x) = ϕ(b) + ∇ϕ(b) · (x − b) +

    ∇ϕ b + t(x − b) − ∇ϕ(b) · (x − b) dt.

(50)

0

We set u =

and use Eq. (24) with k = ∇ϕ(b). We then deduce that:

x−b a

a





−1

e

(ϕ(b)+∇ϕ(b)·(x−b))nθ(b)

R2

x−b ψ a



  dx = eϕ(b)nθ(b) aψ a∇ϕ(b) .

(51)

Combining Eqs. (75), (50) and (51) implies:    cF (a, b) − aψ a∇ϕ(b) A(b) eϕ(b)nθ(b)   

1   x−b dx = A(b) e(ϕ(b)+∇ϕ(b)·(x−b))nθ(b) enθ(b) 0 [∇ϕ(b+t(x−b))−∇ϕ(b)]·(x−b) dt − 1 a−1 ψ a R2



+ A(b)  + R2



e

ϕ(x)nθ(x)

−e

R2

ϕ(x)nθ(b)

   −1 x−b dx a ψ a

    ϕ(x)n x−b θ(x) −1 dx. A(x) − A(b) e a ψ a

Hence     cF (a, b) − aψ a∇ϕ(b) A(b) eϕ(b)nθ(b)       1  x − b  ≤ a−1 A(b)  e[ 0 [∇ϕ(b+t(x−b))−∇ϕ(b)]·(x−b)dt]nθ(b) − 1ψ dx a R2

−1

+a



A(b) 

−1

+a

R2

R2

    ϕ(x)n  x − b  ϕ(x)nθ(b)  θ(x) e −e ψ dx a

      A(x) − A(b)ψ x − b dx.   a

We now give an upper bound of each of these three terms. We bound the first one using the mean value theorem and the triangular inequality: a

−1

 A(b) R2

     [ 1 ([∇ϕ(b+t(x−b))−∇ϕ(b)]·(x−b)) dt]n x − b  θ(b)   e 0 − 1 ψ  dx a 

1

R2

0

≤ a−1 A(b)

!           ∇ϕ b + t(x − b) − ∇ϕ(b) · (x − b)dt ψ x − b  dx.   a

We now use inequality (48) with y = b and h = t(x − b). We then obtain:

a

−1



1

R2

0

A(b)

!           ∇ϕ b + t(x − b) − ∇ϕ(b) · (x − b)dt ψ x − b dx   a

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

470



1 

R2

0

≤ εa−1 A(b) ≤ εa−1 A(b)

 

R2



 !   3    x − b  |x − b| 2 2   |t||x − b| ∇ϕ(b) + M |t| dt ψ dx 2 a

     |x − b|3  x − b  |x − b|2  ∇ϕ(b) + M · ψ  dx 2 6 a

3

 a  a4 ∇ϕ(b)I2 + M I3 ≤ εA(b) 2 6

 (52)

where we set u = x−b a in the last integral. Let us now give an upper bound of the second term. Since   eϕ(x)nθ(x) − eϕ(x)nθ(b) = sin ϕ(x) (nθ(x) − nθ(b) ), one has: a−1 A(b)



R2

≤a

−1

    ϕ(x)n  x − b  θ(x) e − eϕ(x)nθ(b) ψ  dx a 

A(b) R2

≤ εa

−1

    x − b   |nθ(x) − nθ(b) |ψ  dx a

    √     x − b  2M    ψ |x − b| ∇ϕ(b) + |x − b| A(b) 2  dx 2  a R2

  √     M  ≤ εa2 A(b) 2 |u|∇ϕ(b) + a|u|2 ψ(u)du 2 R2

  √   2 3M  ≤ εA(b) 2 a ∇ϕ(b) I1 + a I2 2

(53)

where we have used inequalities (46) and (47) with y = b and h = x − b for the third inequality and set u = (x − b)/a for the fourth one. Finally, using the same approach we prove that the last term is bounded by, using (45): a

−1

 R2

          A(x) − A(b)ψ x − b  dx ≤ ε a2 ∇ϕ(b)I1 + a3 M I2 .   a 2

Combining inequalities (52, 53, 54) leads to:   2       cF (a, b) − aψ a∇ϕ(b) A(b) eϕ(b)nθ(b)  ≤ εa2 A(b) a ∇ϕ(b)I2 + a M I3 2 6   √   M 2   + εa A(b) 2 ∇ϕ(b) I1 + a I2 2     M + εa2 ∇ϕ(b)I1 + a I2 2 which gives (26). 2

(54)

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

471

Proof of point (ii) of Proposition 3.1. Since A ∈ L∞ (R2 ) and I0 , I0 are both finite, one has:  ∂bi cF (a, b) = −

A(x) e

ϕ(x)nθ (x) −2

a

R2



x−b ∂xi ψ a

 dx.

As below we split the sum into three terms −2

∂bi cF (a, b) = −a



 A(b)

e

R2

− a−2 A(b) −2



−a



R2



ϕ(x)nθ (b)

x−b ∂xi ψ a

 dx

   ϕ(x)n (x)  x−b θ dx e − eϕ(x)nθ (b) ∂xi ψ a 



A(x) − A(b) e

ϕ(x)nθ (x)

R2

x−b ∂xi ψ a

 dx.

(55)

We again use (50) and Eq. (24) with k = ∇ϕ(b) to write: −a−2



 e(ϕ(b)+∇ϕ(b)·(x−b))nθ(b) ∂xi ψ R2

x−b a



  dx = eϕ(b)nθ(b) ∂bi ϕ(b)nθ(b) aψ a∇ϕ(b) .

(56)

Hence we deduce that      ∂b cF (a, b) − ∂b ϕ(b)nθ(b) aψ a∇ϕ(b) A(b) eϕ(b)nθ(b)  i i       1 x − b  ≤ a−2 A(b) e( 0 [∇ϕ(b+t(x−b))−∇ϕ(b)]·(x−b)dt)nθ(b) − 1∂xi ψ  dx a R2

+a

+a

−2

−2

   √   x − b   A(b) 2 |nθ(x) − nθ(b) |∂xi ψ dx a  R2

R2

      A(x) − A(b)∂x ψ x − b dx.  i  a

A similar approach than in the proof of point (i) leads to the following upper bound:      ∂bi cF (a, b) − ∂bi ϕ(b)nθ(b) aψ a∇ϕ(b) A(b) eϕ(b)nθ(b)    2   a3 a    ∇ϕ(b) I2 + M I3 ≤ εA(b) 2 6   √   M + εA(b) 2 a∇ϕ(b)I1 + a2 I2 2     M + ε a∇ϕ(b)I1 + a2 I2 2 which leads to the estimate (ii) of Proposition 3.1, with:     2      √   M R2 (a, b) ≤ A(b) a ∇ϕ(b)I2 + a M I3 + 2A(b) + 1 ∇ϕ(b)I1 + a I2 . 2 6 2

2

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

472

6.2. Proof of Theorem 3.2 We can now prove Theorem 3.2. By definition for i = 1, 2:   −1 Λi (a, b) = ∂bi cF (a, b) cF (a, b) . Set B = aψ(a∇ϕ(b))(A(b) eϕ(b)nθ(b) ). Proposition 3.1 implies that:     cF (a, b) − B  ≤ εa2 R1 (a, b),

(57)

    ∂b cF (a, b) − ∂b ϕ(b)nθ (b)B  ≤ εaR2 (a, b), i i

(58)

and

where R1 (a, b) and R2 (a, b) satisfy inequalities (26) and (27) respectively. Then, for i = 1, 2:    −1 Λi (a, b) − ∂bi ϕ(b)nθ(b) = ∂bi cF (a, b) − ∂bi ϕ(b)nθ(b) B + ∂bi ϕ(b)nθ(b) B − cF (a, b) cF (a, b) . Using (57) and (58) implies that:          Λi (a, b) − ∂b ϕ(b)nθ(b)  ≤ ε aR2 (a, b) + a2 R1 (a, b)∂b ϕ(b) cF (a, b)−1 . i i By assumption |cF (a, b)| ≥ εν for all (a, b), which leads to (29). The last estimate (30) follows from the fact that for any given a0 > 0       aR2 (a, b) + a2 R1 (a, b)∇ϕ(b) < ∞.

sup (a,b)∈(0,a0 )×R

Then, since ν ∈ (0, 1/2), there exists some ε0 > 0 depending on a0 such that for any 0 < ε ≤ ε0      aR2 (a, b) + a2 R1 (a, b)∇ϕ(b) ≤ ε2ν−1 , which is equivalent to       ε1−ν aR2 (a, b) + a2 R1 (a, b)∇ϕ(b) ≤ εν . This last inequality and inequality (29) clearly implies inequality (30). 6.3. Proof of Proposition 3.3 The proof of Proposition 3.3 relies on Theorem 3.2 and on the following lemma: Lemma 6.2. Let F an IMMF with accuracy ε > 0. Provided that ε > 0 is sufficiently small, the sign of Re(∂b1 cF (a, b)∂b2 cF (a, b)) is this of ∂b1 ϕ(b)∂b2 ϕ(b). Proof. Let us first observe that Theorem 3.2 implies that for all (a, b) under consideration:   Λ1 (a, b) = ∂b1 ϕ(b)nθ(b) + O εν ,   Λ2 (a, b) = ∂b2 ϕ(b)nθ(b) + O εν .

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

473

Then   Λ1 (a, b)Λ2 (a, b) = ∂b1 ϕ(b)∂b2 ϕ(b) + O εν .

(59)

Hence, for ε > 0 sufficiently small, the sign of Λ1 (a, b)Λ2 (a, b) is this of ∂b1 ϕ(b)∂b2 ϕ(b). To get the required conclusion, we now relate Λ1 (a, b)Λ2 (a, b) and ∂b1 cF (a, b)∂b2 cF (a, b). Let us first remark that the definition of the two Clifford vectors Λ1 (a, b), Λ2 (a, b) from Eq. (28) implies that:     −1  −1  Re Λ1 (a, b)Λ2 (a, b) = Re ∂b1 cF (a, b) cF (a, b) ∂b2 cF (a, b) cF (a, b) .

(60)

We use now the fact that ∂bi cF (a, b) and cF (a, b) are both Clifford vectors. We then apply equality (A.4) with q = ∂b2 cF (a, b) and q  = (cF (a, b))−1 . Hence we get that  −1  −1 ∂b2 cF (a, b) cF (a, b) = cF (a, b) × ∂b2 cF (a, b). Combining this last equality with (60) leads to −2     Re Λ1 (a, b)Λ2 (a, b) = Re ∂b1 cF (a, b) ∂b2 cF (a, b) cF (a, b) .

(61)

Taking into account (61), Eq. (59) reads   Re(∂b1 cF (a, b)∂b2 cF (a, b)) = ∂b1 ϕ(b)∂b2 ϕ(b) + O εν , 2 |cF (a, b)| which implies that Re(∂b1 cF (a, b)∂b2 cF (a, b)) and ∂b1 ϕ(b)∂b2 ϕ(b) have the same sign for ε > 0 sufficiently small. 2 We can now prove Proposition 3.3. Proof of Proposition 3.3. Theorem 3.2 directly implies that for any i = 1, 2 and (a, b) under consideration:     Λi (a, b) − ∂b ϕ(b) ≤ εν . i Since ∂b1 ϕ(b) is assumed to be always positive, we get (31). Since ∂b2 ϕ(b) and ∂b1 ϕ(b)∂b2 ϕ(b) have the same sign, Lemma 6.2 implies (32). 2 6.4. Proof of Theorem 4.1 The proof of Theorem 4.1 consists in several steps. Lemma 6.3. Let F be a function belonging to Aε,d . For any (a, b), there can be at most one  ∈ {1, · · · , L} such that     a∇ϕ (b) − 1 < Δ.

(62)

Remark 6.1. Condition (62) is a necessary condition for having ψ(a∇ϕ  (b)) = 0. Lemma 6.3 then states there is at most one  such that ψ(a∇ϕ (b)) =  0. Then the two following sums, involved in the estimation  of cF (a, b) and ∂bi cF (a, b):

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

474

L 

  A (b) eϕ (b)nθ (b) aψ a∇ϕ (b) ,

=1

and L 

  A (b) eϕ (b)nθ (b) aψ a∇ϕ (b) ∂i ϕ (b)nθ (b) ,

=1

contain at most one term. Proof. We follow the same line as in [22]. Assume that there exists some (1 , 2 ) such that (62) is satisfied with 1 = 2 . One may suppose that 1 < 2 . One then has for j = 1, 2, a

−2

     ∂ϕj (b) 2  ∂ϕj (b) 2     < a−2 (1 + Δ)2 , (1 − Δ) <  + ∂x1  ∂x2  2

which directly implies          ∂ϕ1 (b) 2  ∂ϕ1 (b) 2  ∂ϕ2 (b) 2  ∂ϕ2 (b) 2 −2 2   +  +  +   ∂x1   ∂x2   ∂x1   ∂x2  > 2a (1 − Δ) , and          ∂ϕ2 (b) 2  ∂ϕ2 (b) 2  ∂ϕ1 (b) 2  ∂ϕ1 (b) 2 −2   +  −  −   ∂x1   ∂x2   ∂x1   ∂x2  < 4a Δ. Further, by assumption (34), for any i = 1, 2 one has           2   ∂ϕ2 (b) 2  ∂ϕ1 (b) 2  ∂ϕ2 (b) 2  ∂ϕ2 −1 (b) 2       −  ≥  −  ≥ d  ∂ϕ2 (b)  +  ∂ϕ2 −1 (b)  ,  ∂xi   ∂xi   ∂xi   ∂xi   ∂xi   ∂xi  where we applied assumption (35):      ∂xi ϕ2 (x) − ∂xi ϕ2 −1 (x) ≥ ∂xi ϕ2 (x)nθ

 (x)

     − ∂xi ϕ2 −1 (x)nθ (x)  ≥ d ∂xi ϕ2 (x) + ∂xi ϕ2 −1 (x) .

The usual inequality (u + v)2 ≥ u2 + v 2 , valid for any u, v ≥ 0 implies that     2  2      ∂ϕ2 (b) 2  ∂ϕ1 (b) 2     −  ≥ d  ∂ϕ2 (b)  +  ∂ϕ1 (b)  .  ∂xi   ∂xi   ∂xi   ∂xi  Summing over i = 1, 2 gets: 4a−2 Δ > 2da−2 (1 − Δ)2 ≥ 2da−2 (1 − 2Δ), that is   Δ > d/ 2(1 + d) , which leads to a contradiction. 2

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

475

Lemma 6.4. Let F be a superposition of IMMF with accuracy ε and separation d of the form (33). For all (a, b) ∈ R∗+ × R2 one has   L         ψ a∇ϕ (b) A (b) eϕ (b)nθ (b)  ≤ εa2 Γ1 (a, b), cF (a, b) − a     =1

with Γ1 (a, b) = I1

L  √ =1

+a

I2 2

! L      I2      2A (b) + 1 ∇ϕ (b) + a A (b) ∇ϕ (b) + M 2 =1

L 

M  + a2

=1

I3 6

L 

A (b)M

(63)

=1

(Iα being introduced in (21)). In particular • If for some 0 ∈ {1, · · · , L}     a∇ϕ (b) − 1 < Δ 0

(64)

    cF (a, b) − aψ a∇ϕ (b) A (b) eϕ0 (b)nθ0 (b)  ≤ εa2 Γ1 (a, b). 0 0

(65)

then

• If for all  ∈ {1, · · · , L}, Condition (64) is not satisfied then   cF (a, b) ≤ εa2 Γ1 (a, b).

(66)

Proof. By assumption F is of the form F =

L 

A eϕ n ,

=1

where for each , F = A eϕ n is an IMMF of accuracy ε. Proposition 3.1 applied successively to F1 , · · · , FL and the linearity of the continuous wavelet transform imply:   L       ϕ (b)nθ (b)  ψ a∇ϕ (b) A (b) e cF (a, b) − a  ≤ εa2 Γ1 (a, b).   =1

The rest of the proof follows from Remark 6.1, which states that under condition (64), the following sum L 

  A (b) eϕ (b)nθ (b) aψ a∇ϕ (b) ,

=1

reduces to A0 (b) e

ϕ0 (b)nθ

0

(b)

aψ(a∇ϕ 0 (b)), and 0 if condition (64) is never satisfied. 2

Lemma 6.5. Let F be a superposition of IMMF with accuracy ε and separation d. For all a, b such that     a∇ϕ (b) − 1 < Δ,

(67)

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

476

one has, for any i = 1, 2:     ∂b cF (a, b) − a∂i ϕ (b)nθ (b) ψ a∇ϕ (b) A (b) eϕ (b)nθ (b)  ≤ εaΓ2 (a, b), i 

(68)

with Γ2 (a, b) =

I1

" L # L $     √  ∇ϕ (b) 2A (b) + 1 + aI2 M 2 =1

+

aI2

=1

L 

A (b)

√

" L # $    2  6 2M + ∇ϕ (b) /2 + a I3 A (b) M

=1

(69)

=1

(Iα being introduced in (21)). L Proof. As in the proof of Lemma 6.4, F is of the form F = =1 A eϕ n , where each F = A eϕ n is an IMMF of accuracy ε. Proposition 3.1 applied successively to F1 , · · · , FL and the linearity of the continuous wavelet transform imply:   L        ϕ (b)n  θ  (b)   ∂i ϕ (b)nθ (b) ψ a∇ϕ (b) A (b) e   ≤ εaΓ2 (a, b). ∂bi cF (a, b) − a     =1

Remark 6.1 states that under condition (64), the following sum L 

   ∂i ϕ (b)nθ (b) ψ a∇ϕ (b) A (b) eϕ (b)nθ (b) ,

 =1 ϕ (b)nθ (b) reduces to ∂i ϕ (b)nθ (b) ψ(a∇ϕ ).  (b))(A (b) e

2

Lemma 6.6. Let F be a superposition of IMMF with accuracy ε > 0 and separation d. Let  ∈ {1, · · · , L}. For all a, b such that     a∇ϕ (b) − 1 < Δ,

(70)

      Λi (a, b) − ∂b ϕ (b)nθ (b)  ≤ aε1−ν Γ2 (a, b) + a∇ϕ (b)Γ1 (a, b) . i 0 

(71)

one has, for any i = 1, 2

Proof. The proof follows the same lines than that of Theorem 3.2, replacing the two estimates (57) and (58) on the wavelet coefficients and their partial derivatives, by (65) and (68) respectively. 2 We can now prove Theorem 4.1. Proof of Theorem 4.1. In the proof of Theorem 4.1, the following notations will be needed: amin =

1−Δ , max∈{1,···,L} supb∈R2 (|∇ϕ (b)|)

amax =

1+Δ . min∈{1,···,L} inf b∈R2 (|∇ϕ (b)|)

By assumptions on the phases ϕ of the modes of F , amin and amax are both finite and positive. For any  ∈ {1, · · · , L}, ε% > 0, b ∈ R2 and n ∈ S1 , one also defines:

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

477

' &   B (% ε, n, b) = k ∈ R2 ; maxki n − ∂bi ϕ (b)nθ (b)  ≤ ε% . i

We now show Theorem 4.1. We set ε% = εν and fix  ∈ {1, · · · , L}. Let us first prove that 

 δ,ν SF,ε (b, k, n)dkdn

lim

δ→0 S1 B (% ε,n,b)



a−2 cF (a, b)da.

= 2π

(72)

Aε,F (b)∩{a;maxi |Λi (a,b)−∂bi ϕ (b)nθ (b) |<% ε} δ,ν By definition of SF,ε , we have



 δ,ν SF,ε (b, k, n)dkdn

S1 B (% ε,n,b)



" 



= S1 B (% ε,n,b)

Aε,F (b)

 #  2 ( da k − Re(Λ (a, b) n) i i cF (a, b)δ −2 h dkdn. δ a2 i=1

We first use Fubini Theorem to interchange the order of summations in this integral. Remark that: 





S1 B (% ε,n,b) Aε,F (b)

 ≤

a

   2  ( − Re(Λ (a, b) n) k   da i i h cF (a, b)δ −2  2 dkdn  a δ i=1  



 cF (a, b)

−2 

Aε,F (b)

δ S1



−2

 2   (  ki − Re(Λi (a, b)n)  dk dn da h   δ

i=1

R2

  a−2 cF (a, b)da < +∞,

≤ 2πh2L1 (R) Aε,F (b)

where the last inequality is coming from the fact that F is a finite superposition of IMMFs F , whose wavelet coefficients satisfy (19) and (20) (see Remark 3.3 for more details). Then by Fubini Theorem: 

 δ,ν SF,ε (b, k, n)dkdn

S1 B (% ε,n,b)

 =

a

−2

"



cF (a, b)

δ S1 B (% ε,n,b)

Aε,F (b)

−2

#   2 ( ki − Re(Λi (a, b)n) dkdn da. h δ i=1

Since the integrand is bounded by:  "   −2 a cF (a, b) 



S1 B (% ε,n,b)

δ

−2

#   2  (   ki − Re(Λi (a, b) n)  dkdn  ≤ 2πh2L1 (R) a−2 cF (a, b), h  δ i=1

which belongs to L1 (Aε,F (b)), the Dominant Convergence Theorem applies. Hence,

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

478



 δ,ν SF,ε (b, k, n)dk dn

lim

δ→0 S1 B (% ε,n,b)

 =

"

a−2 cF (a, b) lim





δ→0

Aε,F (b)

 =

S1 B (% ε,n,b)

" a−2 cF (a, b) lim



#   2 ( − Re(Λ (a, b) n) k i i dk dn da δ −2 h δ i=1 

δ→0 S1 B (% ε,n,b)∩B (δR,n,b)

Aε,F (b)

#   2 ( − Re(Λ (a, b)n) k i i dk dn da δ −2 h δ i=1

where we set B  (δR, n, b) = {k ∈ R2 ; maxi |ki n − Λi (a, b)| ≤ δR}, with R the bound of the support of h. Remark now that, if maxi |Λi (a, b) − ∂bi ϕ (b)nθ (b) | > ε%, then for δ sufficiently small B (% ε, n, b) ∩ B  (δR, n, b) = ∅, and the limit above is equal to 0. On the other side, if maxi |Λi (a, b) − ∂bi ϕ (b)nθ (b) | < ε%, for δ sufficiently small: B (% ε, n, b) ∩ B  (δR, n, b) = B  (δR, n, b), and 



S1

B (δR,n,b)

lim

δ

δ→0

  = lim

δ

δ→0 S1

R2

−2

−2

  2 ( ki − Re(Λi (a, b) n) dk dn h δ i=1

  2 ( ki − Re(Λi (a, b) n) dk dn h δ i=1

= 2π,

where the last equality comes from the assumption R h = 1. One then deduces equality (72). Observe now that for any ε > 0 sufficiently small there exists some a0 ≥ amax such that −σ ν −σ a 2πF B˙ ∞,1 0 ≤ε

and

sup (a,b)∈(0,a0

)×R2

    ε1−ν a Γ2 (a, b) + a∇ϕ (b)Γ1 (a, b) ≤ εν .

(73)

Indeed, by definition of Γ1 (a, b), Γ2 (a, b) given respectively in Eqs. (63) and (69), there exists CΓ > 0 not depending on a0 such that sup (a,b)∈(0,a0

)×R2

    a Γ2 (a, b) + a∇ϕ (b)Γ1 (a, b) ≤ CΓ a40 .

This inequality implies that to prove (73), it is sufficient to find a0 such that  1/σ  1/4 amax ≤ 2πF B˙ −σ ε−ν ≤ a0 ≤ CΓ−1 ε2ν−1 . ∞,1

(74)

Since ν < 1/(2 + 4/σ), one deduces that ε(2ν−1)/4 /ε−ν/σ tends to ∞ when ε → 0. Hence, for ε sufficiently small, there exists some a0 ≥ amax satisfying (74). From now to the end of the proof, we fix such a0 .

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

479

We now split the integral in the right-hand side of (72) into two terms 

 δ,ν SF,ε (b, k, n)dk dn

lim

δ→0 S1 B (% ε,n,b)



= 2π

cF (a, b)

da a2

cF (a, b)

da . a2

Aε,F (b)∩{a∈(0,a0 ); maxi |Λi (a,b)−∂bi ϕ (b)nθ (b) |<% ε}



+ 2π Aε,F (b)∩{a≥a0 ; maxi |Λi (a,b)−∂bi ϕ (b)nθ (b) |<% ε}

(75)

We now deal with each of these two terms successively. Recall that since by assumption F is a superposition of IMMFs, then inequality (19) holds for some σ > 0 (see Remark 3.3). We then get that for any a0 satisfying (73): 

I% = 2π

cF (a, b)

Aε,F (b)∩{a≥a0 ;maxi |Λi (a,b)−∂bi ϕ (b)nθ (b) |<% ε}

da a2

+∞    cF (a, b) da ≤ 2πF  ˙ −σ a−σ ≤ ε%. ≤ 2π B∞,1 0 a2

(76)

a0

Let us now consider the integral  2π

cF (a, b)

Aε,F (b)∩{a∈(0,a0 );maxi |Λi (a,b)−∂bi ϕ (b)nθ (b) |<% ε}

da , a2

and prove that under assumption (73) ' &   Aε,F (b) ∩ a ∈ (0, a0 ), maxΛi (a, b) − ∂bi ϕ (b)nθ (b)  < ε% i       = Aε,F (b) ∩ a ∈ (0, a0 ); a∇ϕ (b) − 1 < Δ .

(77)

Let us first assume that a ∈ Aε,F (b) ∩ {a ∈ (0, a0 ); |a|∇ϕ (b)| − 1| < Δ}. Since |a|∇ϕ (b)| − 1| < Δ, Lemma 6.6 and assumption (73) imply that       maxΛi (a, b) − ∂bi ϕ (b)nθ (b)  ≤ aε1−ν Γ2 (a, b) + a∇ϕ (b)Γ1 (a, b) ≤ ε%, i

that is a ∈ {a ∈ (0, a0 ), maxi |Λi (a, b) − ∂bi ϕ (b)nθ (b) | < ε%}. Conversely, let us assume that a ∈ Aε,F (b) ∩ {a ∈ (0, a0 ); maxi |Λi (a, b) − ∂bi ϕ (b)nθ (b) | < ε%}. Since a ∈ Aε,F (b), |cF (a, b)| > ε%. Then by Lemma 6.4 there exists some  ∈ {1, · · · , L} such that |a|∇ϕ (b)| − 1| < Δ, otherwise assumption (73) and Eq. (66) would lead to a contradiction. In addition, remark that if  = , one has       maxΛi (a, b) − ∂bi ϕ (b)nθ (b)  ≥ ∂bi ϕ (b)nθ (b) − ∂bi ϕ (b)nθ (b)  − maxΛi (a, b) − ∂bi ϕ (b)nθ (b) . i

i

As above, Lemma 6.6 and assumption (73) imply that   maxΛi (a, b) − ∂bi ϕ (b)nθ (b)  ≤ ε%. i

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

480

Hence, we get that       maxΛi (a, b) − ∂bi ϕ (b)nθ (b)  ≥ d ∂bi ϕ (b) + ∂bi ϕ (b) − ε% i

≥ ε%, provided that ε is sufficiently small, namely that:  ε≤

   d  ∂bi ϕ (b) + ∂bi ϕ (b) 2

 ν1 i.e.

    d ∂bi ϕ (b) + ∂bi ϕ (b) ≥ 2˜ ε,

∀,  , ∀i.

(78)

Hence, necessarily,  =  and then a ∈ Aε,F (b) ∩ {a ∈ (0, a0 ), |a|∇ϕ (b)| − 1| < Δ}. In addition, one has 

cF (a, b)a−2 da

2π Aε,F (b)∩{a∈(0,a0 ),|a|∇ϕ (b)|−1|<Δ}





cF (a, b)a−2 da − 2π

= 2π {a∈(0,a0 ),|a|∇ϕ (b)|−1|<Δ}

cF (a, b)a−2 da,

{a∈(0,a0 ),|a|∇ϕ (b)|−1|<Δ}\Aε,F (b)

% ≤ ε%. Using this last equation and Eqs. (75) and (76), we deduce that with |I|     lim 1 δ→0 C˜ ψ

S1

  2π ≤  C˜ψ

 δ,ν SF,ε (b, k, n)dkdn B (% ε,n,b)

− A (b) e

  cF (a, b)a−2 da − A (b) eϕ (b)nθ (b) 



{a∈(0,a0 ),|a|∇ϕ (b)|−1|<Δ}

  2π +  C˜ψ

  

ϕ (b)nθ (b) 

  cF (a, b)a−2 da + ε%.



{a∈(0,a0 ),|a|∇ϕ (b)|−1|<Δ}\Aε,F (b)

Observe that the domain {a ∈ (0, a0 ), |a|∇ϕ (b)| − 1| < Δ} is bounded from below by amin which only depends on F . In addition, if a ∈ / Aε,F (b), inequality (66) is satisfied. Hence    



{a∈(0,a0 ), |a|∇ϕ (b)|−1|<Δ}\Aε,F (b) %

−2

cF (a, b)a

 a0     * )    2 −2 da ≤ ε sup a Γ1 (a, b) a da   a∈(0,a0 ) amin

) * ≤ ε sup a2 Γ1 (a, b) a−1 min a∈(0,a0 )

≤ ε1−ν

)

* sup a2 Γ1 (a, b) a−1 min

a∈(0,a0 )

≤ C ε%,

(79)

where C = (amin inf b∈R2 |∇ϕ (b)|)−1 . The last inequality comes from assumption (73) valid for ε > 0 sufficiently small. This last inequality and inequality (65) proved in Lemma 6.4 then imply that

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

    lim 2π δ→0 C˜ ψ

S1

  2π ≤  C˜ψ

 δ,ν SF,ε (b, k, n)dkdn B (% ε,n,b)

− A (b) e

  

ϕ (b)nθ (b) 

      −2 ϕ (b)nθ (b) ϕ (x)nθ (x)  aψ a∇ϕ (b) A (b) e a da − A (x) e 



{a∈(0,a0 ), |a|∇ϕ (b)|−1|<Δ}

  2π +  C˜ψ

481





   ε. εa2 Γ1 (a, b) a−2 da + (C + 1)%

{a∈(0,a0 ), |a|∇ϕ (b)|−1|<Δ}

To conclude, observe now that the first term of the right-hand side vanishes. Indeed, the wavelet is real isotropic and for any b ∈ R2 , a → ψ(a∇ϕ  (b)) is supported in (0, amax ). Since a0 ≥ amax , we then have:  {a∈(0,a0 ),|a|∇ϕ (b)|−1|<Δ}

 da  = ψ a∇ϕ (b) a

 R

 da  1 = ψ a∇ϕ (b) a 2π

 C˜ψ ψ(ξ) , dξ = 2 |ξ| 2π

R2

where C˜ψ has been defined in (6). The second term of the last inequality can be bounded using the same approach than in inequality (79). We then get Theorem 4.1. 2 Acknowledgments The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-13-BS03-0002-01 (ASTRES). Appendix A. Quaternionic calculus In this appendix, we give a short introduction to the algebra of quaternions H. More details can be found for instance in [47]. H is the real 4-D vector space spanned by {1, i, j, k}, i.e. a quaternion is of the form q = q0 +q1 i +q2 j +q3 k, where the algebra product is defined by i2 = j2 = k2 = −1 and ij = −j i = k, j k = −kj = i, k i = −i k = j. For a given quaternion q = q0 + q1 i + q2 j + q3 k, one defines: • its real part: Re(q) = q0 , • its vectorial part: Vect(q) = (q1 , q2 , q3 ) ∈ R3 . If Re(q) = 0, q = q1 i + q2 j + q3 k is called a pure quaternion. The conjugate of the quaternion q is q = q0 − q1 i − q2 j − q3 k, and its norm (or modulus) is defined by: |q| =



+ qq =

q02 + q12 + q22 + q32

For any (q, q  ) ∈ H2 , one has: qq  = q  q

    and qq   = |q|q  .

Additionally, the quaternions are a division algebra, which means that any (non-zero) quaternion admits a multiplicative inverse given by: q −1 =

q . |q|2

(A.1)

482

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

Observe that one has: q −1 = |q|q 2 = (q)−1 . The set of unit quaternion {q ∈ H; |q| = 1} will be denoted by S3 . The exponential function is defined on the quaternion algebra as exp : q → eq =

+∞   q =0

!

,

which converges since exp(|q|) converges. Using the exponential map, each unit quaternion of S3 (|q| = 1) such that Vect(q) = 0 can be written as: q = (cos ϕ + n sin ϕ) = eϕn ,

(A.2)

with n=

Vect(q) , |Vect(q)|

  cos ϕ = Re(q) and sin ϕ = Vect(q),

which is an extension of the complex exponential. Notice that the usual property eμ eν = eμ+ν is no more satisfied in general. The polar form of any quaternion q is then given by: q = |q|(cos ϕ + n sin ϕ) = |q| eϕn

(A.3)

where n is a pure unit quaternion: n = a i + b j + c k with a2 + b2 + c2 = 1, and ϕ ∈ R. If (ϕ, n) ∈ R × S3 satisfies (A.3), one says that ϕ is a scalar argument of q and n a vectorial orientation of q. A quaternion of the form q = q0 + q1 i + q2 j with (q0 , q1 , q2 ) ∈ R3 (q4 = 0) is called a Clifford vector. Observe that if q = a0 + a1 i + a2 j and q  = a0 + a1 i + a2 j with (a0 , a1 , a2 ), (a0 , a1 , a2 ) ∈ R3 then qq  = q  × q.

(A.4)

In this case any vectorial orientation of q is of the form: n = ai + bj with a2 + b2 = 1. Then there exists θ ∈ R such that n = cos θ i + sin θ j. This implies that q admits the following polar form: q = |q|(cos ϕ + sin ϕ cos θ i + sin ϕ sin θ j) = |q| eϕ(cos θ i+sin θj) .

(A.5)

If (ϕ, θ) ∈ R2 satisfies (A.5), one says that ϕ is a scalar argument of q and θ a scalar orientation of q. Observe that ϕ and θ are also Euler angles of the vector ⎛

⎞ q0 ⎝ q1 ⎠ ∈ R3 . q2 Appendix B. Hilbert and Riesz transform In this appendix we briefly recall the Hilbert transform and the analytic signal analysis, and then proceed with a natural extension in two dimension, the Riesz transform, used to define the 2-D monogenic signal in Section 1.

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

483

B.1. The 1-D Hilbert transform We first define the Hilbert transform of a 1-D real valued function. Definition B.1. Let f ∈ L2 (R, R). The Hilbert transform of f , denoted as Hf is defined in the Fourier domain as follows: for a.e. ξ ∈ R, ,(ξ) = −i sgn(ξ)f (ξ). Hf Definition B.2. The complex analytic signal associated with f ∈ L2 (R, R) is then F = f + i Hf, which reads in the Fourier domain as   F = 1 + sgn(ξ) f . We recall below some properties of the Hilbert transform and the analytic signal. Proposition B.1. The Hilbert transform is an antisymmetric operator on L2 (R, R): Hf1 , f2  = −f1 , Hf2 . Proposition B.2. Let f ∈ L2 (R, R) 1. f and Hf are orthogonal 2. F 2 = 2f 2 3. The spectrum of the analytic signal is one sided. B.2. 2-D Riesz transform We begin with the definition of the 2-D Riesz transform of a real valued image. Definition B.3. Let f ∈ L2 (R2 ). The Riesz transform of f , denoted as Rf is the vector valued function:  Rf =

R1 f R2 f

 ,

where for any i = 1, 2, Ri f is defined in Fourier domain as follows: for a.e. ξ ∈ R2 , ξi , R f (ξ). i f (ξ) = −i |ξ| We now present the key properties of R [48,39]. The first ones concern the invariance with respect to dilations, translations, and the steerability property (relation with the rotations). Proposition B.3. The Riesz transform commutes both with the translation, and the dilation operator: for any f ∈ L2 (R2 ), a > 0 and b ∈ R2 , one has: RDa f = Da Rf

and

RTb f = Tb Rf.

484

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

Proposition B.4. The Riesz transform is steerable, that is, for any f ∈ L2 (R2 ) one has Rθ (Rf ) = rθ−1 R(Rθ f ) =



cos θR1 (Rθ f ) + sin θR2 (Rθ f ) − sin θ R1 (Rθ f ) + cos θ R2 (Rθ f )

 ,

(B.1)

where rθ is the rotation matrix defined by (3). Proof. We will prove in the Fourier domain that R(Rθ f ) = rθ Rθ (Rf ), which is equivalent to Eq. (B.1). Let us first remark that for a.e. ξ ∈ R2 ,   , −1 (R θ f )(ξ) = f rθ ξ .

(B.2)

Using the definition of the Riesz transform in Fourier domain, one has, for a.e. ξ ∈ R2 :   ξ  −1  −1  R(R f )(ξ) = r r × −i ξ f r θ θ θ θ |ξ|  −1  rθ ξ  −1  ξ = −irθ f r θ |rθ−1 ξ|   = rθ (, Rf ) rθ−1 ξ  = rθ (R θ Rf )(ξ), the last relation coming from (B.2). One then deduces the required result. 2 The Riesz transform is also a unitary and componentwise antisymmetric operator on L2 (R2 ): Proposition B.5. For any i ∈ {1, 2}, the i-th component of the Riesz transform Ri is an antisymmetric operator, namely for all f, g ∈ L2 (R2 ): Ri f, gL2 (R2 ) = −f, Ri gL2 (R2 ) .

(B.3)

Since R21 + R22 = −Id, it implies in particular that: Rf, RgL2 (R2 ,R2 ) = R1 f, R1 gL2 (R2 ) + R2 f, R2 gL2 (R2 ) = f, gL2 (R2 )

(B.4)

References [1] D.J. Fleet, A.D. Jepson, Computation of component image velocity from local phase information, Int. J. Comput. Vis. 5 (1) (1990) 77–104. [2] A. Basarab, H. Liebgott, P. Delachartre, Analytic estimation of subsample spatial shift using the phases of multidimensional analytic signals, IEEE Trans. Image Process. 18 (2) (2009). [3] A. Basarab, P. Gueth, H. Liebgott, P. Delachartre, Phase-based block matching applied to motion estimation with unconventional beamforming strategies, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 56 (5) (2009). [4] M. Elshinawy, M. Chouikha, Using AM–FM image modeling technique in mammograms, in: 2003 IEEE 46th Midwest Symposium on Circuits and Systems, vol. 2, 2003, pp. 660–663. [5] J. Elshinawy, M. Zeng, S.-C.B. Lo, M. Chouikha, Breast cancer detection in mammogram with AM–FM modeling and Gabor filtering, in: 2004 7th International Conference on Signal Processing, ICSP ’04, vol. 3, 2004, pp. 2564–2567. [6] I. Kokkinos, G. Evangelopoulos, P. Maragos, Texture analysis and segmentation using modulation features, generative models, and weighted curve evolution, IEEE Trans. Pattern Anal. Mach. Intell. 31 (1) (2009) 142–157. [7] P. Tay, AM–FM image analysis using the Hilbert–Huang transform, in: Southwest Symposium on Image Analysis and Interpretation, SSIAI, 2008, pp. 13–16. [8] A. Belaid, D. Boukerroui, Y. Maingourd, J.F. Lerallut, Phase-based level set segmentation of ultrasound images, IEEE Trans. Inf. Technol. Biomed. 15 (1) (2011) 138–147. [9] M. Felsberg, G. Sommer, The monogenic signal, IEEE Trans. Signal Process. 49 (12) (2001) 3136–3144.

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

485

[10] V. Murray, M.S. Pattischis, E.S. Barriga, P. Soliz, Recent multiscale AM–FM methods in emerging applications in medical imaging, EURASIP J. Adv. Signal Process. 23 (2012) 1–14. [11] T. Qian, Analytic signals and harmonic measures, J. Math. Anal. Appl. 314 (2006) 526–536. [12] N.E. Huang, Z. Wu, S.R. Long, K.C. Arnold, K. Blank, T.W. Liu, On instantaneous frequency, Adv. Adapt. Data Anal. 1 (2009) 177–229. [13] P.V. Rodriguez, M.S. Pattichis, New algorithms for fast and accurate AM–FM demodulation of digital images, in: IEEE International Conference on Image Processing, 2005. [14] R.A. Carmona, W.L. Hwang, B. Torrésani, Characterization of signals by the ridges of their wavelet transforms, IEEE Trans. Signal Process. 45 (10) (1997) 2586–2590. [15] R.A. Carmona, W.L. Hwang, B. Torrésani, Multiridge detection and time–frequency reconstruction, IEEE Trans. Signal Process. 47 (2) (1999) 480–492. [16] I. Daubechies, S. Maes, A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models, in: A. Aldroubi, M. Unser (Eds.), Wavelets in Medicine and Biology, CRC Press, 1996. [17] F. Auger, P. Flandrin, Improving the readability of time–frequency and time-scale representations by the reassignment method, IEEE Trans. Signal Process. 43 (5) (1995) 1068–1089. [18] N.E. Huang, Z. Shen, S.R. Long, S.C. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. A 454 (1998) 903–995. [19] P. Flandrin, G. Rilling, P. Gonçalves, Empirical mode decomposition as a filter bank, IEEE Signal Process. Lett. 11 (2) (2004) 112–114. [20] G. Rilling, P. Flandrin, One or two frequencies? The empirical mode decomposition answers, IEEE Trans. Signal Process. 56 (1) (2008) 85–95. [21] G. Rilling, P. Flandrin, P. Gonçalves, On empirical mode decomposition and its algorithms, in: IEEE EURASIP Workshop on NSIP, Grado, Italy, 2003. [22] I. Daubechies, J. Lu, H.-T. Wu, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl. Comput. Harmon. Anal. 30 (2) (2011) 243–261. [23] H.-T. Wu, P. Flandrin, I. Daubechies, One or two frequencies? The synchrosqueezing answer, Adv. Adapt. Data Anal. 3 (1–2) (2011) 29–39. [24] S. Meignen, T. Oberlin, S. Mclaughlin, A new algorithm for multicomponent signal analysis based on synchrosqueezing: with an application to signal sampling and denoising, IEEE Trans. Signal Process. 60 (11) (2012) 5787–5798. [25] J.C. Nunes, Y. Bouaoune, E. Delechelle, O. Niang, P. Bunel, Image analysis by bidimensional empirical mode decomposition, Image Vis. Comput. 21 (2003) 1019–1026. [26] C. Damerval, S. Meignen, V. Perrier, A fast algorithm for bidimensional EMD, IEEE Signal Process. Lett. 12 (10) (2005) 701–704. [27] J.C. Nunes, E. Delechelle, Empirical mode decomposition: applications on signal and image processing, Adv. Adapt. Data Anal. 1 (1) (2009) 125–175. [28] A. Linderhed, Image empirical mode decomposition: a new tool for image processing, Adv. Adapt. Data Anal. 1 (2) (2009) 265–294. [29] E.H.S. Diop, R. Alexandre, L. Moisan, Intrinsic nonlinear multiscale image decomposition: a 2D empirical mode decomposition-like tool, Comput. Vis. Image Underst. 116 (2012) 102–119. [30] K.G. Larkin, D.J. Bone, M.A. Oldfield, Natural demodulation of two-dimensional fringe patterns. I. General background of the spiral phase quadrature transform, J. Opt. Soc. Amer. A 18 (8) (2001) 1862–1870. [31] K.G. Larkin, D.J. Bone, M.A. Oldfield, Natural demodulation of two-dimensional fringe patterns. II. Stationary phase analysis of the spiral phase quadrature transform, J. Opt. Soc. Amer. A 18 (8) (2001) 1871–1881. [32] M. Felsberg, Low-level image processing with the structure multi-vector, PhD thesis. Christian-Albrechts-Universitat zu Kiel, Germany, 2002. [33] X. Guanlei, W. Xiaotong, X. Xiaogang, Improved bi-dimensional EMD and Hilbert spectrum for the analysis of textures, Pattern Recognit. 42 (2009) 718–734. [34] D. Zhang, J. Pan, Y.Y. Tang, Image analysis based on an Improved Bidimensional Empirical Mode Decomposition Method, in: International Conference on Wavelet Analysis and Pattern Recognition, ICWAPR, 2010, pp. 144–149. [35] M. Bülow, G. Sommer, A novel approach to the 2D analytic signal, in: F. Solina, A. Leonardis (Eds.), Proc. CAIP’99, Ljubljana, Slovenia, 1999, pp. 25–32. [36] H. Yang, Y. Ying, Synchrosqueezed wave packet transform for 2d mode decomposition, SIAM J. Imaging Sci. 6 (4) (2013) 1979–2009. [37] Y. Yang, T. Qian, F. Sommen, Phase derivative of monogenic signals in higher dimensional spaces, Complex Anal. Oper. Theory 6 (2012) 987–1010. [38] S.C. Olhede, G. Metikas, The monogenic wavelet transform, IEEE Trans. Signal Process. 57 (9) (2009) 3426–3441. [39] M. Unser, D. Van De Ville, Multiresolution monogenic signal analysis using the Riesz–Laplace wavelet transform, IEEE Trans. Image Process. 18 (11) (2009) 2402–2418. [40] J.P. Antoine, R. Murenzi, P. Vandergheynst, S.T. Ali, Two-dimensional Wavelets and Their Relatives, Cambridge University Press, 2004. [41] S.C. Olhede, G. Metikas, Multiple multidimensional Morse wavelets, IEEE Trans. Signal Process. 55 (3) (2007) 921–936. [42] H. Triebel, Interpolation Theory, Function Spaces, Differential Operators, North-Holland, Amsterdam, 1978. [43] S. Jaffard, Y. Meyer, R.D. Ryan, Wavelets: Tools for Science and Technology, SIAM, 2001. [44] Z. Wu, N.E. Huang, X. Chen, The multi-dimensional ensemble empirical mode decomposition method, Adv. Adapt. Data Anal. 1 (03) (2009) 339–372.

486

M. Clausel et al. / Appl. Comput. Harmon. Anal. 39 (2015) 450–486

[45] G. Thakur, E. Brevdo, N.S. Fučkar, H.-T. Wu, The synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications, Signal Process. 93 (5) (2013) 1079–1094. [46] Y.-C. Chen, M.-Y. Cheng, H.-T. Wu, Non-parametric and adaptive modelling of dynamic periodicity and trend with heteroscedastic and dependent errors, J. R. Stat. Soc. Ser. B Stat. Methodol. 76 (3) (2014) 651–682. [47] A. Sudberry, Quaternionic analysis, Math. Proc. Cambridge Philos. Soc. 85 (1979) 199–225. [48] E.M. Stein, Singular Integrals Differentiability Properties of Functions, second edition, Princeton University Press, New York, 1970.