P.-O Mattei

SPECTRUM OF A FLUID-LOADED VIBRATING PLATE: THE MULTIPLE RESONANCE PHENOMENON

It was recently observed in a numerical study on a high order perturbation method under heavy fluid loading that a loaded vibrating plate results, not only in the classical frequency shift of the in vacuo single resonance (in both the real part because of the fluid added mass and the imaginary part because of energy lost by radiation), but also in an increase in the number of the resonance. As a result of the loading, a single in vacuo resonance of the structure is transformed into a multiple resonance. Here we show that this phenomenon is a refinement of the Sanchez's classical result where it was established, using asymptotic analysis, that in the case of a light loading conditions the scattering frequencies of a fluid loaded elastic structure (ie the resonance frequencies) are nearly the real eigenfrequencies of the elastic body alone and the complex scattering frequencies of the fluid with a rigid solid. A theoretical explanation of the multiple resonances is given using classical results on theory of entire functions. It is established that every single in vacuo resonance of a simply supported rectangular plate is transformed into an infinite number of resonances under fluid-loading condition.

Introduction

Accurately describe fluid structure interactions increases the computational cost of the various numerical methods used. Most numerical methods are leading to the resolution of linear systems of simultaneous equations. When the loading conditions are taken into account, the size of the matrices involved increases considerably, and they become full and frequency dependent. Any method of reducing these drawbacks is therefore most welcome; and one of the best methods available is that based on the perturbation approach. Asymptotic analysis can be conduced when the loading is light, as in the case of a metal structure surrounded by air. This approach consists in introducing a small parameter ǫ, which is the ratio between surface mass density and fluid volume density. Using perturbation expansions, one can then construct an approximate solution, based on the in vacuo eigenmodes, which is not only easier to calculate than the exact solution but also leads to a better understanding of the phenomena involved [START_REF] Filippi | The rôle of the resonance modes in the response of a fluid-loaded structure[END_REF].

The light fluid loading perturbation method (involving high density structures in contact with a small density fluid) was recently extended [START_REF] Mattei | High-Order Perturbation Expansion for the Spectral Analysis of Fluid-Loaded Vibrating Structure[END_REF] to cases where the perturbation parameter becomes large (which corresponds to a light and/or thin structure in contact with a high density fluid). During this study multiple resonance processes were observed where a single resonance mode can have several resonance frequencies. The spectrum of the operator can no longer be described in terms of the discrete spectrum based on a countable series of resonance mode/resonance frequency pairs. This typically nonlinear phenomenon is closely linked to the concept of non-linear frequency modes ( [START_REF] Dazel | Frequency non-linear modes[END_REF]), where the non linearity depends on the time parameters (such as those involved in porous materials) rather than on the more usual geometrical parameters. This frequency non-linearity is related to problems which are written in the form O(ω)U = S, where u is the unknown factor and S is the source term. O(ω) is an operator which depends non linearly on the angular frequency ω, as in the case of vibroacoustics, where the coupling depends non-linearly on the frequency via the Helmholtz equation Green's kernel.

One shows that this behavior is a general property of at least fluid-loaded simply supported rectangular plates. The main results obtained in our previous study [START_REF] Mattei | High-Order Perturbation Expansion for the Spectral Analysis of Fluid-Loaded Vibrating Structure[END_REF] are recalled in paragraph 2. A theoretical description of multiple resonance in the case of a simply supported plate loaded with a non vanishing density fluid is given in paragraph 3. It is then established, using classical results on the connection between the order of an entire function and the distribution of its zeros, that any resonance mode in this structure will have an infinite number of resonance frequencies. Numerical examples of a steel plate in contact with air or water are given paragraph 4. The conclusions and possible extensions of this work are presented in paragraph 5.

Statement of the problem

Consider a finite elastic structure occupying a domain Σ, the behavior of which is described by a differential operator A (in the case treated in this communication it is a simple bi-Laplacian operator). This structure is loaded with a perfect fluid at rest extending to infinity. The equation of the fluid-loaded structure is given by the classical integrodifferential equation

AU(M, ω) -ρ s hω 2 U(M, ω) -ǫ Σ U(M ′ , ω)G(M; M ′ , ω)dM ′ = F (M) (1)
where ǫ = ρ f /ρ s h is a small parameter in the case of a metal structure (volume mass density ρ s and thickness h) in contact with air (volume mass density ρ f ). G(M; M ′ , ω) is the Green's function of the Neumann problem in the Helmholtz equation outside the surface Σ. F (M) is the source. To facilitate reading, the weak formulation of the problem is introduced. One defines the radiation impedance

β(U, V, ω) = Σ Σ U(M ′ )G(M; M ′ , ω)V * (M)dMdM ′ ,
where U(M) and V (M) are functions defined on the surface Σ, a(U, V * ) stands for the potential elastic energy of the structure, U, V * is the usual inner product and ρ s hω 2 U, V * is the kinetic energy. Then the weak formulation reads: find U(M) such that for every

V (M) one has a(U, V * ) -Λm (ω) ( U, V * -ǫβ(U, V * , ω)) = F, V * . ( 2
)
2.1 Resonance modes and eigenmodes for a fluid-loaded structure

The eigenmodes Ũm (M, ω) and eigenvalues Λm (ω) (or equivalently the eigenpulsations ωm (ω) given by Λm (ω) = ρ s hω 2 m (ω)) are the non zero solutions of the homogeneous weak formulation [START_REF] Dazel | Frequency non-linear modes[END_REF], that is find Ũm (M, ω) and ωm (ω) such that for every V (M)

a( Ũm , V * ) -ρ s hω 2 m (ω) Ũm , V * -ǫβ( Ũm , V * , ω) = 0, (3)
β(U, V, ω), that characterizes the energy loss by radiation, depends on the frequency. Then the eigenmodes, well suited to deal with harmonic regime, are complex and depend on the frequency. One introduces the resonance modes, that allow to compute the transient regime in a natural way; they are complex and do not depend on the frequency [START_REF] Filippi | The rôle of the resonance modes in the response of a fluid-loaded structure[END_REF]. The resonance modes Ûm (M) and resonances pulsations ωm , given by Λm = ρ s hω 2 m , are the non zero solutions of the following problem: find Ûm (M) and ωm such that for every V (M)

a( Ûm , V * ) -ρ s hω 2 m Ûm , V * -ǫβ( Ûm , V * , ωm ) = 0. ( 4
)
These two kind of modes, that are identical when the loading is not took into account, have a strong link as each resonance pulsation is the solution of the resonance equation

ω2 m (ω m ) = ω2 m , (5)
and each resonance mode is an eigenmode calculated at the corresponding resonance pulsation.

Ûm (M) = Ũm (M, ωm ). ( 6
)
It is worth noting that in the eigenmodes are easier to compute than the resonance modes.

High order perturbation method

When the loading parameter ǫ is small, the eigenmodes and eigenvalues can be calculated using a perturbation method. The following perturbation expansion is introduced into the weak formulation (3)

Ũm (M) = Ũ(0) m (M) + ǫ Ũ(1) m (M) + • • • + ǫ s Ũ(s) m (M) + • • • Λm = Λ(0) m + ǫ Λ(1) m + • • • + ǫ s Λ(s) m + • • •
where Ũ(0) m (M) and Λ(0) m are the eigenmodes and eigenvalues of the elastic structure in vacuo. To facilitate the reading, one poses:

β mn (ω) = β Ũ(0) m , Ũ (0)
* n , ω . It can then be easily established that one has to the order 1 in ǫ:

Λ(1) m (ω) = Λ(0) m β mm (ω),
to the order 2 in ǫ:

Λ(2) m (ω) = Λ(0) m (β mm (ω)) 2 + p =m Λ(0) m Λ(0) m - Λ(0) p β mp (ω) ,
and so on to the higher orders. Next, as previously shown [START_REF] Mattei | High-Order Perturbation Expansion for the Spectral Analysis of Fluid-Loaded Vibrating Structure[END_REF], by grouping the terms and by a re-summation technique, one shows that the high order perturbation expansion is given by

Λm (ω) ≈ Λ(0) m 1 -ǫβ mm (ω) , (7)
which remains valid even for higher values of ǫ. It is worth noting that, compared to the classical perturbation expansion of the eigenvalues given by Λm (ω) ≈ Λ(0) m (1 + ǫβ mm (ω)), the high-order approximation [START_REF] Mangiarotty | Acoustic Radiation Damping of Vibrating Structure[END_REF] is particularly interesting as it gives very precise results for the same numerical cost.

The calculation of the eigenmodes and eigenvalues does not requires any particular effort and, for a given ω, each eigenpulsation is the solution of the following eigenpulsation equation

ω2 m (ω) ≈ ω(0)2 m (1 -ǫβ mm (ω)) . (8)
Equation [START_REF] Mattei | A two-dimensional Tchebycheff collocation method for the study of the vibration of a baffled fluid-loaded rectangular plate[END_REF] shows that the calculation of the various eigenpulsations do not requires any computational cost but the evaluation of multiple integrals β mm (ω).

The resonance equation

By introducing in equation ( 5) the high order expansion of the eigenpulsation given by equation ( 8), the resonance equation ( 5) is approximated as

ω2 m (1 -ǫβ mm (ω m )) ≈ ω(0)2 m . (9)
This shows that the calculation of the resonance pulsations is completely different than that of the eigenpulsations since β mm (ω) has a rather complicated dependence on the frequency and solving equation ( 9) requires the search for the zeros of non-convex complex function which is a challenging task. A numerical example given in the paper [START_REF] Mattei | High-Order Perturbation Expansion for the Spectral Analysis of Fluid-Loaded Vibrating Structure[END_REF] has shown a curious phenomenon of multiple resonance. That particular example has showed that the resonance frequency of a particluar mode is not unique at least two resonance frequencies could exist simultaneously for a given geometrical configuration.

To have a better understanding of this non-usual phenomenon the roots of the resonance equation ( 9) are now studied in a more general manner using classical results in the framework of the distribution of zero of entire functions [START_REF] Ya | Distribution of Zeros of Entire Functions[END_REF][START_REF] Ya | Lectures on Entire Functions[END_REF]. In the following the resonance equation ( 9) is written as f mm (z):

f mm (z) = z 2 (1 -ǫβ mm (z)) -z 2 0m , (10)
where z = ωm and z 0m = ω(0) m . The zeros of f mm (z) are the resonance pulsations.

3 Multiple resonances

Methodology

Now all the problem focuses on the existence and the counting of the roots of f mm (z). More or less, the theory of the zeros of entire functions (analytic functions in the whole complex plane) is based on the Great Picard's Theorem [START_REF] Saks | Analytic Functions[END_REF] which states that if f (z) has an essential singularity at a point z w then on any open set containing z w , f (z) takes on all possible values, with at most one exception, infinitely often. The exception is needed as showed by the function exp(1/z) that has an essential singularity at z = 0 but never attains 0 as a value. The second point is that every entire function is either a polynomial or has an essential singularity at infinity. Now if one can prove that the resonance function f mm (z) is an entire function, that is not a polynomial and that the value 0 is not the exception needed by the Picard's theorem, then the resonance function f mm (z) has an infinite number of zeros. Then for every mode, an infinite number of resonance exists. Now to characterize complex function f (z) as an entire function, the simplest way is to study its development as a power series of the form f

(z) = n=∞ n=0 c n z n . When lim n→∞ n √ c n = 0, f (z) is an entire function.
To show that it is not a polynomial, one can study its order. Let us recall that if f (z) is an entire function, the function M(r) = max |z|=r |f (z)| will increase indefinitely together with r. The rate of growth of M(r) can be characterized by comparing it with function exp(r). By definition, the order (of growth) ρ is the smallest number such that M(r) ≤ exp r ρ+η for every η > 0 and for r > r 0 (η). From this, it can be easily seen that the order is given by the formula ρ = lim sup r→∞ ln ln M(r)/ ln r ; as examples, a polynomial in z is of order 0, and the exponential function is of order 1. The order ρ of an entire function depends on the asymptotic behavior of the coefficients of its series expansion. One has 1/ρ = lim n→∞ inf(ln

(1/|c n |)/(n ln n)) ([11]); or in other words, if c n ∝ 1/(n!) α then ρ = 1/α.
One can simplify again the problem under consideration. The resonance function can be written as

f mm (z) = z 2 -z 2 0m -ǫz 2 β mm (z), as a sum of a polynomial (z 2 -z 2 0m
) of order 0 and a function which is the product of a polynomial (ǫz 2 ) also of order 0 and the radiation impedance β mm (z) which order is to be determined. The distribution of zeros of f mm (z) can be deduced from that of β mm (z) as long as ǫ has a non zero value. If β mm (z) is a polynomial, the resonance function is also a polynomial with a finite number of roots. If β mm (z) is an entire function (in the case considered below it is of order 1) with essential singularity at infinity, then f mm (z) is also an entire function with essential singularity at infinity and by the Picard's theorem f mm (z) can takes on all possible values at infinity, with at most one exception, infinitely often.

The main results of this paper is, as shown below on the particular example of a simply supported plate, that the radiation impedance and the resonance function are entire functions of exponential type (that is of order one and normal type [START_REF] Ya | Distribution of Zeros of Entire Functions[END_REF]) that both have an infinite number of roots.

The resonance function of a simply supported baffled plate as an entire function of the frequency

Now let us consider a thin simply supported rectangular plate (dimensions a × b with thickness h) made of isotropic material with Young's modulus E, Poisson's coefficient ν and density ρ s . In a vacuum, the eigenmodes are given by Ũ(0) mn (x, y) = 2/ √ ab sin mπx/a sin nπy/b and the eigenpulsations by ω 2 and c is the speed of waves in the fluid. In the following, the pulsation dependence in ω is changed into the equivalent wavenumber dependence in k = ω/c. Since the in vacuo modes involve a double index, the notation for the radiation impedance is simplified as β(

(0) mn = π 2 D/(ρ s h) [(m/a) 2 + (n/b) 2 ] with D = Eh 3 /12(1-ν 2 ). When the plate is baffled, the radiation impedance β( Ũ(0) mn , Ũ(0) * mn , ω) is given as β( Ũ(0) mn , Ũ(0) * mn , ω) = - ab π 1 0 1 0 1 0 1 0 sin mπx sin nπy sin mπx ′ sin nπy ′ e ı ω c d d dxdydx ′ dy ′ , ( 11
) where d = a 2 (x -x ′ ) 2 + b 2 (y -y ′ )
Ũ(0) mn , Ũ(0) * mn , ω) = β mn (k). The resonance function is denoted by f mn (k) = k 2 (1 -ǫβ mn (k)) -k 2 0mn , with k 0mn = ω (0) mn /c.
It is worth noting that β mn (k) and f mn (k) decrease in the upper complex plane (with ℑ(k) > 0, the term exp(ıkd) has an exponential decrease when |k| → ∞) and conversely grow in the lower complex plane, then all the roots of these two functions are located in the lower complex half plane. As the time dependency had been chosen as exp(-ıωt), this ensures that the movement of the plate remains finite over time [START_REF] Filippi | The rôle of the resonance modes in the response of a fluid-loaded structure[END_REF].

By making a change of variable [START_REF] Mangiarotty | Acoustic Radiation Damping of Vibrating Structure[END_REF], it can be easily shown that β mn (k) is given by

β mn (k) = - 4ab π 1 0 1 0 G m (X)G n (Y ) e ıkD D dXdY, (12)
where

D = √ a 2 X 2 + b 2 Y 2 .
The functions G m (x) and G n (y) are given by single analytic integrals

G m (x) = 1 0 sin mπ(x+x ′ ) sin mπx ′ dx ′ with G m (x) = (sin mπx + mπ(1 -x) cos mπx) /(2mπ).
By making a change of variables from rectangular to polar [START_REF] Habault | A numerical method for the computation of the resonance frequencies and modes of a fluid-loaded plate: application to the transient response of the system[END_REF], with 0 < θ α = arctan(b/a) < π/2, the radiation impedance ( 12) is written as

β mn (k) = - 4 π θα 0 J 1mn (k, θ)dθ - 4 π π/2 θα J 2mn (k, θ)dθ, (13)
where J 1mn (k, θ) and J 2mn (k, θ) are defined, with

H mn (R, θ) = G m R cos θ a G n R sin θ b
, as

J 1mn (k, θ) = a cos θ 0 H mn (R, θ)e ıkR dR, J 2mn (k, θ) = b sin θ 0 H mn (R, θ)e ıkR dR. (14)
By expanding the exponential as a power series, the two integrals ( 14) are given as

J 1mn (k, θ) = ℓ=∞ ℓ=0 (ık) ℓ ℓ! a cos θ 0 H mn (R, θ)R ℓ dR, J 2mn (k, θ) = ℓ=∞ ℓ=0 (ık) ℓ ℓ! b sin θ 0 H mn (R, θ)R ℓ dR.
As H mn (R, θ)R ℓ is a continuous integrable function of R in ]0, a/ cos θ[ or in ]0, b/ sin θ[ and of θ in ]0, π/2[, the mean value theorem ensures that there exist

R 1 ∈]0, a/ cos θ[, R 2 ∈]0, b/ sin θ[, θ 1 ∈]0, θ α [ and θ 2 ∈]θ α , π/2[ such that θα 0 J 1mn (k, θ)dθ = ℓ=∞ ℓ=0 (ık) ℓ ℓ! H mn (R 1 , θ 1 )R ℓ 1 , π/2 θα J 2mn (k, θ)dθ = ℓ=∞ ℓ=0 (ık) ℓ ℓ! H mn (R 2 , θ 2 )R ℓ 2 .
One then obtains β mn (k) as an entire function of the variable k:

β mn (k) = - 4 π H mn (R 1 , θ 1 ) ℓ=∞ ℓ=0 (ıkR 1 ) ℓ ℓ! + H mn (R 2 , θ 2 ) ℓ=∞ ℓ=0 (ıkR 2 ) ℓ ℓ! . ( 15
)

Order of the radiation impedance

To calculate the order of the impedance radiation, one studies the comportment of the function M β mn (r) = max |k|=r |β mn (k)| for r → ∞. Let us begins with successive integrations by parts of the integrals given by equations (14). Introducing the results given in the appendix leads to

J 1mn (k, θ) = -1 4ık + h 2mn (θ) (ık) 3 - h 1 3mn (θ)e ık a cos θ -h 0 3mn (θ) (ık) 4 + a cos θ 0 ∂ 4 H mn (R, θ)
∂R 4 e ıkR (ık) 5 dR, (16)

J 2mn (k, θ) = -1 4ık + h 2mn (θ) (ık) 3 - h 2 3mn (θ)e ık a cos θ -h 0 3mn (θ) (ık) 4 + b sin θ 0 ∂ 4 H mn (R, θ)
∂R 4 e ıkR (ık) 5 dR, (17

)
As one studies the comportment of max |k|→∞ |β mn (k)|, in the equations ( 16) and ( 17), the first two power terms in k are negligible compared to the exponential ones when |k| → ∞; similarly, only the first exponential term has a contribution when |k| → ∞. Then

M β mn (r) = max |k|=r 1 |k| 4 θα 0 h 1 3mn (θ)e ık a cos θ dθ + π 2 θα h 2 3mn (θ)e ık b sin θ dθ . ( 18
)
Now by expanding the exponentials in equation (18) as power series, one obtains

M β mn (r) = max |k|=r 1 |k| 4 ℓ=∞ ℓ=0 (ık) ℓ ℓ! θα 0 h 1 3mn (θ) a cos θ ℓ dθ + π 2 θα h 2 3mn (θ) b sin θ ℓ dθ . (
a cos θ ℓ dθ = h 1 3mn (θ 1 ) a cos θ 1 ℓ , 0 < θ 1 < θ α , (20)
π 2 θα h 2 3mn (θ) b sin θ ℓ dθ = h 2 3mn (θ 2 ) b sin θ 2 ℓ , θ α < θ 2 < π/2. ( 21
) Let us denote c 1 = a/ cos θ 1 and c 2 = b/ sin θ 2 . It is easy to see a ≤ c 1 ≤ √ a 2 + b 2 and b ≤ c 2 ≤ √ a 2 + b 2 . M β mn (r) is then written as M β mn (r) = max |k|=r 1 |k| 4 h 1 3mn (θ 1 ) ℓ=∞ ℓ=0 (ıkc 1 ) ℓ ℓ! + h 2 3mn (θ 2 ) ℓ=∞ ℓ=0 (ıkc 2 ) ℓ ℓ! = max |k|=r 1 |k| 4 h 1 3mn (θ 1 )e ıkc 1 + h 2 3mn (θ 2 )e ıkc 2 . (22)
The previous result shows that the function M β mn (r), characterizing the growth of the radiation impedance, is the sum of two exponential functions divided by a polynomial of the fourth degree. Its order is easy to obtain from classical results on the sum and product of entire functions (see for example [START_REF] Saks | Analytic Functions[END_REF], chapter VII paragraph 6). There is two possibilities for characterizing the growth of M β mn (r) for r sufficiently large: the first is that if σ denotes a finite positive number and C a positive constant one has M β mn (r) ≈ C exp(σr) which implies that the radiation impedance is of order 1, the second is that M β mn (r) = 0 which implies that the radiation impedance is of order 0. From the expressions of the functions h 1 3mn (θ) and h 2 3mn (θ) given in the appendix, it is obvious that for a = b and for m = n one cannot have h 1 3mn (θ 1 )e ıkc 1 = -h 2 3mn (θ 2 )e ıkc 2 . If, for example c 1 > c 2 , then

M β mn (r) = max |k|=r 1 |k| 4 h 1 3mn (θ 1 )e ıkc 1 . (23)
If a = b and m = n care must be taken; in this case θ α = π/4 and a little algebra shows that

π 4 0 h 1 3mm (θ) a cos θ ℓ dθ = π 2 π 4 h 2 3mm (θ) a sin θ ℓ dθ. ( 24
)
Then the maximum of β mm (k) is

M β mm (r) = max |k|=r 2 |k| 4 h 1 3mm (θ 1 )e ıkc 1 . (25)
Then β mn (k) is an entire function of order 1 and the resonance function

f mn (k) = k 2 (1 -ǫβ mn (k)) -k 2
0mn is an also an entire function of order 1 for all ǫ = 0.

Distribution of zeros of the radiation impedance and of the resonance function

In addition to the order, the relations ( 22), ( 23) and ( 25) give a more refined information on the distribution of zeros of the radiation impedance. An entire function that grows as exp(σr), with 0 < σ < ∞, is an entire function of order one and normal type σ. Such functions are known as entire function of exponential type (EFET) [START_REF] Ya | Lectures on Entire Functions[END_REF]. The resonance function and the radiation impedance are EFET. They both possess an additional property on the distribution of their zeros. They are functions that belong to the Cartwright class C

[6]. A EFET f (z) belongs to the class C if +∞ -∞ log + |f (t)| 1 + t 2 dt < ∞, (26)
where the function log + |w| is defined as log + |w| = max(0, log |w|) or equivalently as log + |w| = 1/(2π)

As H mn (r, θ) is continuous integrable function of r and θ the two integrals that appear in the right member of equation ( 27) are finite. β mn (k) is bounded for every k real and belongs to class C. Showing that the resonance function

f mn (k) = k 2 (1 -ǫβ mn (k)) + k 2 0mn
belongs to class C requires a more refined estimation as, because of the k 2 term, it is unbounded at infinity. Nevertheless, as β mn (k) is bounded for every k, the resonance function is also bounded for every finite real value k, then for every finite value A,

+A -A log + |f mn (t)| 1 + t 2 dt < ∞. ( 28
)
It remains to estimate

±∞ ±A log + |f mn (t)| 1+t 2
dt. As seen in the previous paragraph, the radiation impedance for real values of k → ∞ has the following asymptotic approximation

β mn (k) ≈ 1 (ık) 4 θα 0 h 1 3mn (θ)e ık a cos θ dθ + π 2 θα h 2 3mn (θ)e ık b sin θ dθ . ( 29
)
The two integrals I 1 (k) = θα 0 h 1 3mn (θ)e ık a cos θ dθ and I 2 (k) = π 2 θα h 2 3mn (θ)e ık b sin θ dθ that appear in this equation can both be approximated for k → ±∞ by the stationary phase method [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers I[END_REF]:

I 1 (k) ≈ k→±∞ h 1 3mn (0) π 2ka e ıka+ı π 4 = (-1) m m 2 π 2 2a 3 π 2ka e ıka+ı π 4 , I 2 (k) ≈ k→±∞ h 2 3mn π 2 π 2kb e ıkb+ı π 4 = (-1) n n 2 π 2 2b 3 π 2kb e ıkb+ı π 4 .
Now the radiation impedance equation has the following approximation for k → ±∞

β mn (k) ≈ 1 (ık) 4 (-1) m m 2 π 2 2a 3 π 2ka e ıka+ı π 4 + (-1) n n 2 π 2 2b 3 π 2kb e ıkb+ı π 4 , (30)
when k → ±∞, β mn (k) tends to zero and the resonance function can be approximated by

f mn (k) ≈ k 2 . For k → ±∞ one has log + |f mn (k)| = log |k 2 |.
Then for every value of

A +∞ A log + |f mn (t)| 1 + t 2 dt ≈ 2 +∞ A log t 1 + t 2 dt = 2 -∞ -A log |t| 1 + t 2 dt = 2 1 + log A A + O(A -2 ) < ∞.
(31) Combining results of equations ( 28) and (31) shows that [START_REF] Ya | Lectures on Entire Functions[END_REF]) the density of the set of zeros n ± (r, α)/r has the limit lim r→∞ n ± (r, α)/r = σ/π. This result is a very strong one, since not only it ensures that the resonance function possesses an infinite number of zeros but also that "almost all" zeros lie in the neighborhood of the real axis.

+∞ -∞ log + |f mn (t)| 1+t 2 dt < ∞.
It is worth noting that the zeros are counted with their multiplicity order. Then the zeros of the resonance function appears as conjugate pairs and then must be counted twice. Actually, the resonance function

f mn (k) = k 2 (1 -ǫβ mn (k)) + k 2 0mn can be ap- proached for |k| → ∞ by f mn (k) = k 2 (1 -ǫ2h 1 3mm (θ 1 )/k 4 e ıkc 1 + k 2 0mn
, then if k is a solution of f mn ( k) = 0, it is easy to see that -k⋆ , the opposite of its conjugate satisfies f mn (-k⋆ ) = 0.

Numerical examples

In this section, some numerical examples of the multiple resonance phenomenon are presented. As it is impossible to obtain a closed formula for J 1mn (k, θ) and J 2mn (k, θ) given by equations ( 14), a numerically efficient approximation of these functions has been constructed. As shown below, this approximation preserves the order properties of the involved functions. Since one always has 0 < θ α < π/2, the functions to be integrated are expanded with respect to the angular variable θ to the order 1. This gives

J 1mn (k, θ) ≈ j 1m (k) = a ıa 3 k 3 + 2m 2 π 2 -ıakm 2 π 2 -(-1) m 2e ıak m 2 π 2 4(ak -mπ) 2 (ak + mπ) 2 , (32)
J 2mn (k, θ) ≈ j 2n (k) = b ıb 3 k 3 + 2n 2 π 2 -ıbkn 2 π 2 -(-1) n 2e ıbk n 2 π 2 4(bk -nπ) 2 (nk + nπ) 2 . ( 33
)
After reporting the two expansions (32) and (32) in the radiation impedance given in equation ( 13), the following approximation of β mn (k), denoted by βmn (k), is obtained

βmn (k) = - 4 π (j 1m (k)θ α + j 2n (k)(π/2 -θ α )). (34)
It is worth noting that the accuracy of the approximation of β mn (k) decreases as the mode order increases. However, at low orders (say m ≤ 3, n ≤ 3) it remains sufficient. It is easy to obtain

j 1m (k = ±mπ/a) = a ±3ı + mπ 16mπ , and
j 1m (k = 0) = a 1 -(-1) m 2m 2 π 2 , (35)
j 2n (k = ±nπ/b) = b ±3ı + nπ 16nπ , and
j 2n (k = 0) = b 1 -(-1) n 2n 2 π 2 . ( 36
)
Then as j 1m (k) and j 2n (k) are entire functions of k, βmn (k) is also an entire function of k. It still remains to estimate the order of βmn (k). For given m and n and for |k| → ∞, one has

lim |k|→∞ j 1m (k) = -(-1) m m 2 π 2 e ıak 2a 3 k 4 and lim |k|→∞ j 2n (k) = -(-1) n n 2 π 2 e ıbk 2b 3 k 4 , then lim |k|→∞ βmn (k) = 2π (-1) m m 2 a 3 θ α e ıak + (-1) n n 2 b 3 π 2 -θ α e ıbk 1 k 4 . (37)
It is worth noting that the approximation of the radiation impedance given in equation (34) or in equation (37) presents a similar comportment than the exact one: a sum of exponential functions divided by a polynomial of the fourth degree; also the lower the mode order, the better the approximation βmn (k) becomes. To confirm that the approximation keeps the order of βmn (k) to one, classical results on the sum and product of entire functions are used again. First it is worth noting that the entire function (-1) m m 2 /a 3 θ α e ıak + (-1) n n 2 /b 3 (π/2 -θ α ) e ıbk cannot be identically zero for any value of a, b, m or n when |k| → ∞. To see this more clearly, one can deal separately with the cases a = b and a = b.

The case a = b. One has θ α = π/4 and βmn (k

) ≈ π 2 2a 3 ((-1) m m 2 + (-1) n n 2 ) e ıak 1 k 4
becomes the product of the exponential function e ıak of order 1 and the inverse of a degree 4 polynomial in k, which is of order 0. βmn (k) is therefore an entire function of order 1 and of type a.

The case a = b. Let us take for example a > b. Then with |k| → ∞ one has max e ıak ≫ max e ıbk ; it is worth noting that this is true in the half plane ℑ(k) < 0, whereas in the half plane ℑ(k) > 0, βmn (k) tend rapidly to zero. In the half plane ℑ(k) < 0, one therefore has βmn (k) ≈ 2π a 3 (-1) m m 2 θ α e ıak 1 k 4 which is also the product of the exponential function e ıak of order 1 and the inverse of a degree 4 polynomial in k, which is of order 0. βmn (k) is therefore an entire function of order 1 and of type a. It is worth noting that if b > a its type is b. In other words, one can assume that σ = max(a, b).

n f 11 (K) K = 1 3.1407 , n f 12 (K) K = 1 3.1407 , n f 23 (K) K = 1 3.14268 , n f 22 (K) K = 1 3.14268
This confirms that lim K→∞ n(K, α)/K = 2n f mn (K)/K = 2 max(a, b)/π. Similar results are obtained if for the fluid water is chosen instead of air or for various plate dimensions. All these numeric values are a validation of the theoretical results obtained in the paragraph 3. The resonance relation of the fluid loaded plate has an infinite number of zeros. Therefore, while in a vacuum, each mode of a simply supported plate will have only one resonance frequency (and its opposite), under loading condition each mode will have an infinite number of resonance frequencies. This result is confirmed in figures 1 and 2. In these various contour plots, the dark spots that appear correspond to the resonances of the simply supported plate, that is the roots of fmn (k) = 0. In figures 1, the material of the plate is steel with E = 200 Gpa, ν = 0.3 and ρ s = 7800 Kg/m 3 . In figure 2 the plate is made of aluminium E = 70 Gpa, ν = 0.3 and ρ s = 2700 Kg/m 3 . The dimensions of the two plates under consideration are the same and are given by 1 m × 0.7 m × 8.2 cm. While the spectra of both plates are similar in air (excepted in a higher imaginary part of all resonances of the aluminium plate) since they both have more or less same ratio E/ρ s and more or less the same spectrum in vacuo, when loaded by water the difference between the two spectra is obvious. It is worth noting that the multiple resonance phenomenon makes classical modal analysis rather difficult. This could be seen in figure 2 where it is easy to se that there is two resonances very close together (around ℜ(k) ≈ 18 and ℑ(k) ≈ -3). While in air (or in water for the steel plate) choosing the resonance close to the real axis and neglect the resonances with high imaginary part (corresponding to rapidly vanishing waves) is natural, in water there is, in the present case, two resonances (and then two resonance modes) that have equal importance. Now let us return to the initial question of the existence or non-existence of multiple resonances in fluid-loaded plates. In every boundary condition, one has to look for the resonance pulsations ωmn such that ω2 mn (1 -ǫβ mn ωmn (ω))ω(0)2 mn = 0. Although it has not been proved that such a comportment is true under various boundary conditions, it seems reasonable to assume that under boundary conditions of all kinds every resonance mode in a plate has an infinite number of resonance frequencies as ǫ = 0 (the radiation impedance given by relation [START_REF] Sanchez-Hubert | Vibration and Coupling of Continuous Systems[END_REF] seems to be an EFET for a large variety of functions G m and G n ).

CONCLUSIONS

The loading of a rectangular plate with a high density fluid transforms each resonance into an infinite number of resonances. These results open up new prospects. In particular, it would be interesting to obtain similar estimates on the vibratory modes occurring in clamped plates or mechanical structures of other kinds, although these may not be easy to obtain. Another question which arises is how to determine under what conditions moderately coupled systems (such as a steel plate in contact with water) will have these properties. Another important aspect worth investigating is the validity of modal expansions in the case of resonance modes having two or more resonances frequencies. It would also be interesting to develop experimental device which would make it possible to test the occurrence of the behavior described here.

A Appendix

The function H mn (r, θ) is given by Only a few derivatives with respect to r of H mn (r, θ) up to order 3 taken at the points r = 0, r = a/ cos θ and r = b/ sin θ are not zero. They are given by :

2π0π θα 0 a cos θ 0 HH 4 π θα 0 a cos θ 0 |H0

 040 log |wexp(ıθ)|dθ. A function f (t) bounded on the real axis satisfies log + |f (t)| ≤ |f (t)|. To prove that f (z) belongs to the class C it is sufficient to prove that |f (t)| is bounded for every t ∈ IR and belongs to class C. The radiation impedance given by equation (13) can be majored, for k ∈ IR, by: |β mn (k)| = 4 mn (r, θ)e ıkr dr dθ + mn (r, θ)e ıkr dr dθ ≤ mn (r, θ)| dr dθ + |H mn (r, θ)| dr dθ .

  Then both impedance radiation and resonance relation belong to the Cartwright class C. For 0 < α < π let us denote by n + (r, α) the number of zeros of the function f mn (z) in the sector {z : |z| ≤ r, | arg(z)| ≤ α} and n -(r, α) the number of zeros in {z : |z| ≤ r, |πarg(z)| ≤ α}. By the theorem of Cartwright and Levinson (theorem 1, page 127 of

Figure 1 :

 1 Figure 1: Contour plots of |k 2 (1 -ǫ βmn (k)) -k 2 0mn | in the mode m = 2, n = 3. Steel plate in contact with air ǫ ≈ 0.002 (left) and water ǫ ≈ 1.5 (right).Both radiation impedance βmn (k) and resonance relation f mn (k) belong to the Cartwright class C and then all, but a finite number of, their zeros must lie in the neighborhood of the real axis.As shown above the type σ of βmn (k) is given by σ = max(a, b). Then the type of f mn (k) is given the type of k 2 (1 -ǫβ mn (k)) + k 2 0mn ; and for ǫ = 0 it is given by σ = max(a, b). By the Cartwright and Levinson theorem, the density of resonances lim r→∞ (n + (r, α)/r + n -(r, α)/r) is then given by 2 max(a, b)/π. To verify this, the simplest way is to count the zeros of f mn (z) using the argument theorem on the disk C with radius K centered on the origin. It is clear from equations (35) and (36) that βmn (k) and f mn (k) have no pole inside C. The number of its zeros inside this disk, in the case of the mode (m -n), is denoted by n f mn (K). It is given by the integral n f mn (K) = 1/(2ıπ) C f ′mn (k)/ f mn (k)dk with k = Ke ıθ . This integral is computed numerically without difficulty. As shown before, the zeros occur in combined pairs of negative imaginary part, then each zero must be counted twice; then the density of the zeros given by (n(K, α)/K = n + (K, α)/K + n -(K, α)/K has the asymptotic value lim K→∞ n(K, α)/K = 2n f mn (K)/K.For example, let us take a steel (E = 200 Gpa, ν = 0.3 and ρ s = 7800 Kg/m 3 ) simply supported rectangular plate in contact with air (with density ρ f = 1.2 Kg/m 3 and sound wave celerity c = 340 m/s) and let the plate dimensions be a = 1 m, b = 0.7 m and h = 8.2 cm. It can easily be established that in the mode m = 1, n = 1, one has n f11 (12) = 2, n f11 (18) = 4 and n f11 (30) = 8. In the mode m = 2, n = 3, one has n f23 (19) = 1, n f23 (20) = 3, n f23 (30) = 7. As the radius K of the disk C becomes large one obtains for the density n f mn (K)/K, with e.g. K = 5000

Figure 2 :

 2 Figure 2: Contour plots of |k 2 (1-ǫ βmn (k))-k 2 0mn | in the mode m = 2, n = 3. Aluminium plate in contact with air ǫ ≈ 0.005 (left) and water ǫ ≈ 4.5 (right).

H

  mn (r, θ)

∂ 2 H 2 r=0 = - π 2 4 m 2 cos 2 θ a 2 + n 2 sin 2 θ b 2 =m m 2 π cos 2 θ 2a 3 bn b cos θ nπ cos anπ tan θ b + sin anπ tan θ b -anπ cos anπ tan θ b sin θ = h 1 3mn (θ) ∂ 3 Hn n 2 πsin 2 θ 2ab 3 m a sin θ mπ cos bmπ cot θ a + sin bmπ cot θ a -bmπ cos θ cos bmπ cot θ a = h 2 3mn

 22422132 mn (r, θ) ∂r h 2mn (θ) ∂ 3 H mn (r, mn (r, (θ)

  19) As h 1 3mn (θ) and a/ cos θ are continuous integrable functions for θ ∈ [0, θ α ] and h 2 3mn (θ) and b/ sin θ are continuous integrable functions for θ ∈ [θ α , π/2]. One applies the mean value theorem for the two integrals that appear in equation (19): there exist two values θ 1 and θ 2 such that

	θα
	0	h 1 3mn (θ)

Email address: mattei@lma.cnrs-mrs.fr

Acknowledgments: this work was partially supported by the Agence Nationale de la Recherche (ANR) thru grant ANR-06-BLAN-0081-01.


Sébastien Bougleux
email: sebastien.bougleux@unicaen.fr

Gabriel Peyré
email: peyre@ceremade.dauphine.fr

Laurent Cohen
email: cohen@ceremade.dauphine.fr

Image Compression with Anisotropic Geodesic Triangulations

We propose a new image compression method based on geodesic Delaunay triangulations. Triangulations are generated by a progressive geodesic meshing algorithm which exploits the anisotropy of images through a farthest point sampling strategy. This seeding is performed according to anisotropic geodesic distances which force the anisotropic Delaunay triangles to follow the geometry of the image. Geodesic computations are performed using a Riemannian Fast Marching, which recursively updates the geodesic distance to the seed points. A linear spline approximation on this triangulation allows to approximate faithfully sharp edges and directional features in images. The compression is achieved by coding both the coefficients of the spline approximation and the deviation of the geodesic triangulation from an Euclidean Delaunay triangulation. Numerical results show that taking into account the anisotropy improves the approximation by isotropic triangulations of complex images. The resulting geodesic encoder competes well with wavelet-based encoder such as JPEG-2000 on geometric images.

Introduction

Image compression requires the extraction of meaningful features of an image to minimize the rate distortion curve. A popular tool is the wavelet transform, that offers a sparse representation of natural images. Wavelet coefficients are mapped to a binary code using an adaptive arithmetic encoder such as in the JPEG-2000 standard, see for instance [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Because of their square support, wavelets are however sub-optimal to extract geometric features such as edges or directional textures. Geometric compression. Anisotropic representations such as curvelets [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF] or bandlets [START_REF] Pennec | Sparse geometric image representations with bandelets[END_REF] allows one to improve over the wavelet approximation, but their extension to practical coders remains difficult. Another attractive representation is the design of an adaptive triangulation to approximate the image with linear splines whose supports follow the geo-metric features of the image. Optimal adaptation of triangulations is a NP-hard problem [START_REF] Agarwal | Surface Approximation and Geometric Partitions[END_REF], and one has to use approximate greedy schemes to design the sampling and the triangulation layout. Triangulation simplification. Methods from computer graphics based on mesh simplification are efficient for surface approximation [START_REF] Garland | Surface simplification using quadric error metrics[END_REF]. Similar methods based on thinning together with Euclidean Delaunay triangulations have been applied to image compression [5]. Using an Euclidean metric might be problematic in highly anisotropic regions where an optimal triangulation deviates from an Euclidean Delaunay triangulation. Triangulation refinement. Mesh adaptation using vertex insertion is popular for the resolution of PDEs whose solutions might exhibit singular features [8]. The first attempt to use anisotropic metric for Delaunay refinement [START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF] makes use of an approximate Riemannian geodesic metric. Contributions. We propose to define a truly Riemannian geodesic sampling, and use the anisotropic Delaunay triangulation for image approximation. This combines in a common framework several important features from previous works, including fast Delaunay insertions, anisotropic sampling [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF][START_REF] Moenning | Fast marching farthest point sampling[END_REF], anisotropic triangulations [12,[START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF][START_REF] Peyré | Geodesic remeshing using front propagation[END_REF]2,[START_REF] Bougleux | Anisotropic geodesics for perceptual grouping and domain meshing[END_REF]. We show that these properties are indeed crucial to improve over the sate of the art for the compression of geometric images.

Approximation by Triangulation and Anisotropy

Image approximation by triangulations

It is possible to approximate an image, expressed as a function f ∈ L 2 (Ω), using m linear spline functions defined on a triangulation (V, F), where V = {v i } i=n i=1 ⊂ Ω denotes the vertex set and F ⊂ {1, . . . , m} 3 the triangle set. Each triangle of F is represented as a triplet of vertex indicies. The resulting approximated image is given by

f m = i a i ϕ i , (1)
where ϕ i is linear on each triangle of F, and ϕ i (v j ) = 1 if i = j and 0 otherwise. The coefficients a i are computed such that f m is the othogonal projection of f on the space V generated by the ϕ i , e.g.

f m = argmin g∈V f -g . (2)
This leads to solve the linear system Φa = F , where Φ i,j = ϕ i , ϕ j and F i = f, ϕ i .

To minimize the error ff m relative to the triangulation (V, F), it is necessary to optimize the position of the vertices and the connections between these vertices. We must put more points near areas of strong image gradient, and use strongly anisotropic triangles oriented in the direction of contours. For a cartoon image, which is a C 2 function outside C 2 contours, such an optimization of the triangulation leads to an error in ff m = O(m -1 ), see for instance [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise C 2 singularities[END_REF]. In this case, the side length ∆ 1 and ∆ 2 of the triangles must verify ∆ 2 ≈ ∆ 2 1 . We propose an efficient triangulation algorithm adapted to this problem. This latter result is far better than the wavelet m-terms approximation, based on isotropic triangulations, which generates an error in O(m -1/2 ).

Anisotropic metric

The desired anisotropy of a triangulation (V, F) can be represented by a tensor field H : Ω → R 2×2 . At each point x ∈ Ω, the local metric H(x) characterizes the shape (its size and its orientation) of the mesh elements surrounding x. It can be defined as a symmetric positive definite matrix, and diagonalized as H(x) = λ 1 (x)e 1 (x)e 1 (x) T + λ 2 (x)e 2 (x)e 2 (x) T , [START_REF] Bougleux | Anisotropic geodesics for perceptual grouping and domain meshing[END_REF] where the eigenvalue fields λ 1 and λ 2 satisfying 0 < λ 1 (x) ≤ λ 2 (x), and e 1 , e 2 are the orthogonal unoriented eigenvector fields.

Each triangle of (V, F) conforming to H must be inscribed in an ellipsoid centered at x 0

(x 0 -x) T H(x 0 )(x 0 -x) ≤ t m ,
where t m depends on the number m of triangles in F. Triangle sides are thus aligned with e 1 (x 0 ), and they satisfy the ratio

∆ 1 /∆ 2 ≈ λ 2 (x 0 )/λ 1 (x 0 ).
For a C 2 function f , the ratio λ 2 (x 0 )/λ 1 (x 0 ) can be taken as the ratio |µ 2 (x 0 )/µ 1 (x 0 )| of the eigenvalues of the Hessian of f at x 0 . This improves the constant in the approximation error ff m = O(m -1 ) [START_REF] Mirebeau | Greedy bisection generates optimally adapted triangulations[END_REF]. However, this estimate is quite unstable and cannot be used to treat noisy images, oriented textures, or cartoon images with non-smoothed contours.

So, we propose to use the structure tensor [START_REF] Kass | Analyzing oriented patterns[END_REF] which is computed from the image gradient by

T σ (x) = G σ * T 0 (x), (4)
where G σ is a Gaussian kernel applied to each component of the 1-rank tensor T 0 (x) = (∇ x f ) T ∇ x f . The parameter σ allows a robust estimation in the presence of noise. We have chosen σ = 5/n for an image of size n × n.

From the diagonalization T (x) = µ 1 (x)e 1 (x)e 1 (x) T + µ 2 (x)e 2 (x)e 2 (x) T , we define the triangulation tensor as

H(x) = (ǫ + µ 1 (x)) α e 1 (x) + (ǫ + µ 2 (x)) α e 2 (x), (5)
where ǫ > 0 controls the isotropic adaptivity of the triangulation (the variation of density, see [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF]), and α > 0 controls the anisotropic adaptivity. A metric with a high value of ǫ is almost constant, and a metric with a low value of α is almost isotropic. We have fixed ǫ to 10 -5 and the influence of the parameter α is adapted to the type of images to process.

Geodesic distance

As explained in Section 2.2, at each point x 0 ∈ Ω , the local metric H(x 0 ) defines the size of triangles around x 0 , which must be of constant size t m for the modified metric

x 0 -x H(x0) = (x 0 -x) T H(x 0 )(x 0 -x).
We compute a point set that achieves this condition by sampling them uniformly according to a global metric d H (x, y) between any pair of points (x, y) ∈ Ω × Ω.

The length of a piecewise smooth curve γ : [0, 1] → Ω, according to the metric H, is measured by

L H (γ) = 1 0 γ′(t) H(γ(t)) dt, = 1 0 γ′(t) T H(γ(t))γ′(t)dt.
A curve γ reaching a point γ(t) = x of Ω, with speed γ′(t), has a shorter local length if γ′(t) is approximately collinear to e 1 (x) rather than to another direction. The geodesic distance between two points x and y of Ω is the length of the shortest curve joining x and y

d H (x, y) = min γ∈P (x,y) L(γ), (6)
where P (x, y) = {γ : γ(0) = x ∧ γ(1) = y} denotes the set of piecewise smooth curves joining the points x and y. A curve γ, which satisfies L H (γ) = d H (x, y), is a geodesic. Given a point set S = {x k } k ⊂Ω, its minimal action map is defined as

U S (x) = min x k ∈S d H (x, x k ). (7)
An important theoretical result is if H is a smooth field, U S is the unique viscosity solution of the following Hamilton-Jacobi equation

∇ x U S H(x) -1 = 1, with U S (x k ) = 0, ∀k. (8)
Figure 1. Examples of geodesic distance map US from a central point S = {x0}, together with a few geodesic curves. The metric H is computed from the image f on the left using the structure tensor ( 5) with an increasing value of ǫ and α = .75, so that the metric becomes more isotropic from left to right.

The numerical solution of this equation can be computed on a regular grid of N = n × n pixels. In the case of an isotropic metric H(x) = λ(x)Id 2 , the Fast Marching (FM) algorithm, introduced by Sethian [START_REF] Hjelle | Triangulations and Applications[END_REF], allows to efficiently solve a discretization of (8) in O(N log N ) operations. Several extensions of the FM algorithm have been proposed to solve (8) for a generic metric, see for instance [18]. Figure 1 shows examples of geodesics distance and minimal paths.

Anisotropic Geodesic Meshes

Computing a triangulation (V, F) conformed to a tensor field H is a difficult optimization problem. The conformity constraint on V imposes that the vertices must be uniformely distributed according to the geodesic distance

d H . This implies that d H (v i , v j ) is approximately constant if (v i , v j ) is an edge of a face of F. A uniform distribution of the vertices satisfies ∀x ∈ Ω, U V (x) ≤ δ, (9)
where δ > 0 controls the density of the vertex set. Such a sample set of Ω can be computed using a farthest point sampling strategy.

Anisotropic farthest point sampling

The farthest point strategy (FPS) has been introduced by Eldar et al. [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF] for image sampling. It is an iterative method which starts with an initial sample set V ⊂ Ω. Each iteration places a new sample at the farthest point p ∈ Ω from the previously computed samples

p = argmax q∈Ω U V (q), (10)
until condition (9) is reached. This sampling strategy is summed up by the following algorithm:

Algorithm: FPS(V, ǫ) 1. Find the point p ∈ Ω satisfying (10). 2. If U V (p) < ǫ, then exit, otherwise set V ← V ∪ {p}.
3. Update the distance map U V , and goto 1.

Step 3 only requires a local update of the distance map in the region {x\U {q} (x) ≤ U V-{q} (x)}. This makes the whole seeding process becoming faster as new points are inserted.

In the case of the Euclidean distance, the FPS is linked to the incremental construction of Voronoi diagrams. As pointed out by Eldar et al. [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF], the farthest point from the set V corresponds to the farthest vertex of the Voronoi diagram of V, which is to the center of the largest ball not containing any point of V. This nice property allows to compute efficiently the farthest point from the Voronoi diagram. However, it does not extend to the case of geodesic distances, for which the farthest point does not necessarily correspond to a vertex of the geodesic Voronoi diagram. Even if this drawback is not taken into account, computing geodesic Voronoi diagrams is not possible with classical computational geometry techniques (see for instance [12,[START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF]).

To sample images with the FPS, Eldar et al. [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF] proposed to assign a weight to the vertices of the (Euclidean) Voronoi diagram with an application-dependent function, and to select the next point to include in V as the Voronoi vertex having the maximum weight. Later, Peyré and Cohen [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF], and Moenning and Dodgson [START_REF] Moenning | Fast marching farthest point sampling[END_REF], proposed independently to approximate the geodesic distance with the fast marching method, leading to a fast marching FPS (FM-FPS).

With an anisotropic fast marching algorithm, the proposed extension of the FM-FPS to anisotropic metrics is obvious. Nevertheless, it should be useful for the resolution of problems which need a step of point selection.

Using the structure tensor (4) as a metric, the anisotropic FM-FPS concentrates the samples near the most important image structures, see Figure 2 for a comparison with the isotropic FM-FPS (computed with H(x) = µ 1 (x)Id 2 ).

In the sequel, we present a modified version of the FPS algorithm which generates triangulations of images conformed to a metric H.

Anisotropic farthest point meshing

We have modified the farthest point strategy so that the faces of the triangulation (V, F) are aligned with the direction of anisotropy e 1 . This is imposed by defining F from the dual complex of the anisotropic Voronoi diagram of the sample set V.

Anisotropic Voronoi diagram. From the anisotropic geodesic distance d H , the anisotropic Voronoi cell V i (V) generated by a point x i ∈ V is defined as the set of points p ∈ Ω for which the geodesic between p and x i is shorter than the geodesics between p and the other points of V

V i (V) = {x ∈ Ω : d H (x, x i ) ≤ d H (x, x j ), ∀s j ∈ V}.
The union of all Voronoi cells, for i = 1, . . . , n, partition Ω to form the bounded anisotropic Voronoi diagram of V, noted V (V). This diagram is also known as the Riemannian Voronoi diagram [12].

Anisotropic Delaunay complex.

The dual complex of V (V) is the anisotropic Delaunay triangulation of (V, noted D(V). Contrary to the (Euclidean) Delaunay triangulation, its anisotropic version D(V) does not necessarily define a valid planar triangulation:

-Some vertices can be isolted (connected to only one other vertex of D(V)).

-Some edges can overlap.

Also, when V is the result of the FPS, D(V) is not necessarily a triangulation of Ω. Indeed, the boundary ∂Ω should always be represented in D(V).

Farthest point meshing. To construct a valid planar triangulation of Ω according to the density constraint (9), we initialize the set V with at least the four corner points of ∂Ω, such that the boundary of D(V) corresponds to ∂Ω.

Then the FPS is modified as follows.

The farthest point x k (see [START_REF] Kass | Analyzing oriented patterns[END_REF]) is inserted in

V if V k (V) ∩ ∂Ω = ∅.
Otherwise, there exists a boundary edge

(x i , x j ) ∈ D(V) which satisfies V k (V) ∩ (x i , x j ) = ∅.
In this case, the edge is said to be encroached by x k . To ensure that ∂Ω is always represented in D(V), the encroached edge is split into two sub-edges (x i , x k ) and ( x k , x j ), where

x k = argmin x∈(xi,xj )∩V k (V) d H (x, x k ).
Once the FPS with boundary constraints has been computed, the proposed algorithm iteratively inserts points to V while the triangulation is not planar:

-If a vertex x i ∈ D(V) is isolated, it is connected to one vertex x j ∈ D(V). In this case, the point xi = argmin x∈Vi(V)∩Vj (V) d H (x, x i ) is inserted to V. -If a triangle (x i , x j , x l ) is inverted, its dual Voronoi vertex V i,j,l = V i (V) ∩ V j (V) ∩ V l (V) is inserted to V.
When there is no inverted triangle, then no edge can overlap.

The farthest point meshing method is summed up by the following algorithm:

Algorithm: FPM(V, ǫ) 1. Find the point x k ∈ Ω satisfying (10). 2. If U V (x k ) < ǫ goto 5. 3. If x k encroaches a boundary edge (x i , x j ) ∈ D(V) Split (x i , x j ) into (x i , x k ) and ( x k , x j ) V ← V ∪ { x k }, update U V and goto 1. 4. V ← V ∪ {x k }, update U V and goto 1. 5. If ∃x i ∈ D(V) isolated, then x k ← xi and goto 3. 6. If ∃ inverted triangle (x i , x j , x l ) ∈ D(V)
x k ← V i,j,l and goto 3.

This algorithm is similar to the anisotropic Delaunay reffinement meshing algorithm of Labelle and Shewchuk [START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF]2]. The main difference is we use a truly anisotropic geodesic distance, meanwhile they use the modified metric (3). This is also the strategy adopted by Bougleux et al. [START_REF] Bougleux | Anisotropic geodesics for perceptual grouping and domain meshing[END_REF], but in a different context. This allows to obtain better results, in particular for triangulations with a few number of vertices, where the approximation with (3) can be important.

Image Approximation and Compression

Spline approximation

Once the image f has been triangulated with the anisotropic farthest point meshing method, its spline approximation f m , given by (1), is progressively refined as the value of m increases. Figure 5 (left column) shows examples of triangulations for two values of m and α = 0.75.

As shown in Figure 4, we found the error ff m minimal for α ≈ 0.75, with f m the image of Figure 5 (bottom left), and this for a wide range of m. Following an experimental study, α ≈ 0.75 seems to be an acceptable choice for a wide range of natural images. From now on, that is the value we are going to use in the paper.

Figure 3 shows the evolution of SNR, according to m, for an anisotropic Delaunay triangulation. It is compared to the SNR of the Euclidean Delaunay triangulation of the sample set generated with α = 0.75. This shows the importance to compute both an appropriate sampling and a suitable connectivity to obtain a better approximation.

Quantification

To encode the approximation f m , we quantify the position of the vertices of V with a setp T v , as well as the coefficients a i with a step T a . Thus, we obtain rounded integer values vi = [a i / T v ] and āi = [a i / T a ]. The coefficients ā and v are encoded with an arithmetic encoder on R v + R a bits. The decoder can compute the decompressed positions V = { v i } i and coefficients by v i = ( vi + sign( vi ) / 2)T v and a i = ( āi + sign( āi ) / 2)T a . The decompressed approximation is thus given by

f m = i a i ϕ i ,
where ϕ i is the spline basis assocated to the vertex v i .

The triangulation connectivity must also be encoded, on R c bits. In order to reduce this number R c , we transform the anisotropic Delaunay triangulation ( V, F) to the Euclidean Delaunay triangulation of V, which can be constructed by the encoder and the decoder. We use the local optimiza-tion procedure of Lawson (see for instance [9]) to compute this transformation. Given an interior edge (v i , v j ) ⊂ F incident to the triangles (v i , v j , v k ) and (v i , v j , v l ), this edge is replaced by the edge (v k , v l ) if the triangles are not in the Euclidean Delaunay triangulation. The number of bits needed to encode F is reduced to R c = K log 2 (A), where K is the number of substitutions and A the number of edges in F.

Other image compression methods, based on triangulations, have already been proposed in litterature, see for instance [5]. These methods begin with a dense triangulation which is decimated to reduce the approximation error. Our approach is different, it progressively refines the triangulation, which is very fast for low values of m. Moreover, it uses a Delaunay triangulation based on an anisotropic geodesic distance, where the elongation of triangles easily adapts to the image structures, contrary to the Euclidean Delaunay triangulation used for instance in [5].

The choice of T v and T a can be optimized to minimize the curve of coding rate ff m according to the number of bits R = R v + R a + R c . To simplify matters, we have used R v = 10 -2 and R a = f ∞ /100. Figure 5 (left) compares compressions given by our method and JPEG-2000 (see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]). This shows that for low coding rate, and for images with sufficient geometric structures, our method improves the current standard.

Conclusion and Perspectives

We have presented a new image compression method based on anisotropic geodesic Delaunay triangulations. The anisotropy of images is exploited to generate triangulations through a farthest point sampling approach. Images are then approximated by linear splines on the triangulations. Finally, the compression is achieved by coding both the coefficients of the spline approximation and the deviation between the geodesic Delaunay triangulation and the Euclidean one. By experimentation, we have shown that the resulting geodesic encoder competes well with waveletbased encoder such as JPEG-2000.

Since anisotropic geodesic distances can also be computed on meshes with fast marching methods, the proposed approach can be extended, with some modifications, to compress meshes.

Figure 2 .

 2 Figure 2. Results of the FM-FPS algorithm for 2000 samples.

Figure 3 .

 3 Figure 3. For α = 0.75, SNR according to m for the anisotropic Delaunay triangulation (solid) and for the Euclidean Delaunay triangulation (dotted).

Figure 4 .

 4 Figure 4. For m = 230, 260, 300, SNR according to α (bottom).


Javier Ripollés-Melchor

Ane Abad-Motos

Maurizio Cecconi

Rupert Pearse

Samir Jaber

Karem Slim

Nader Francis

Antonino Spinelli

Jean Joris

Orestis Ioannidis

Association between use of enhanced recovery after surgery protocols and postoperative complications in colorectal surgery in Europe: The EuroPOWER international observational study

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.