Kévin Polisano
email: kevin.polisano@univ-grenoble-alpes.fr

Marianne Clausel
email: marianne.clausel@univ-lorraine.fr

Valérie Perrier
email: valerie.perrier@univ-grenoble-alpes.fr

Laurent Condat
email: laurent.condat@gipsa-lab.grenoble-inp.fr

Riesz-based orientation of localizable Gaussian fields

Keywords: Fractional fields, H-sssi fields, anisotropy function, Riesz analysis, structure tensor, orientation vector, localizable fields, tangent fields

come

Introduction

Anisotropic images, admitting different characteristics along a considered direction, are ubiquitous in many areas as computer vision [START_REF] Peyré | Oriented patterns synthesis[END_REF], image processing [START_REF] Jennane | Analyse de champs browniens fractionnaires anisotropes, XXVIIIème colloque GRETSI[END_REF], and hydrology [START_REF] Benson | Aquifer operator scaling and the effect on solute mixing and dispersion[END_REF]. A major issue is then the definition of a suitable concept of local anisotropy.

A widely used approach, in the image processing community, consists in defining directionality properties of an image by means of its Riesz transform [START_REF] Felsberg | The monogenic signal[END_REF]. Several characteristics can then be derived from the knowledge of the Riesz transform of an image: its local orientation, which is roughly speaking the dominant direction at a given point and the structure tensor [START_REF] Bigun | Optimal orientation detection of linear symmetry[END_REF] whose rank is related to the local dimensionality of the image. This approach has proved to be successful for many applications such as classification or texture retrieval [START_REF] Depeursinge | Lung texture classification using locally-oriented Riesz components[END_REF]. Recently, this framework has been extended to the case of superimposed patterns. An extension of the synchrosqueezing method to the bidimensionnal setting, based on wavelet analysis, has been proposed in [START_REF] Clausel | The Monogenic Synchrosqueezed Wavelet Transform: A tool for the Decomposition/Demodulation of AM-FM images[END_REF].

In many cases, the analyzed anisotropic image is related to some physical phenomena, that can be well modeled using a stochastic approach. Anisotropic random fields are then naturally involved in the modeling of medical images [START_REF] Jennane | Analyse de champs browniens fractionnaires anisotropes, XXVIIIème colloque GRETSI[END_REF] or in spatial statistics [START_REF] Davies | Fractal analysis of surface roughness by using spatial data (with discussion)[END_REF]. In such situations, the Riesz framework is not so easy to apply. The main difficulty lies in giving a rigourous definition of the Riesz transform of a random field. Indeed [START_REF] Stein | Singular Integrals Differentiability Properties of Functions[END_REF], the Riesz transform of a function is well defined if it belongs to L p for some p > 1, which is not the case for the sample paths of many classical random fields, like Fractional Fields widely used to model random textures. The nonlocal character of the Riesz transform then prevents any definition based on a restriction of the considered random field to a compact set.

To overcome all these difficulties, we choose to use a Riesz-based approach, adapting the work of [START_REF] Bigun | Optimal orientation detection of linear symmetry[END_REF][START_REF] Felsberg | Low-Level Image Processing with the Structure Multivector[END_REF] about Riesz-analysis of anisotropic images characteristics. In [START_REF] Felsberg | Low-Level Image Processing with the Structure Multivector[END_REF], the authors defined a structure tensor and an orientation of the neighborhood of an image, derived from the Riesz framework. Numerical experiments have put in evidence the effectiveness of the approach and especially the fact that one clearly recovers the anisotropic features of the image. From the theoretical point of view, the orientation and the structure tensor depend on the chosen analysis function. We show that, surprisingly, considering the very general case of localizable Gaussian fields, the anisotropic characteristics of a random field become intrinsic: neither the structure tensor nor the orientation vector depend on the analysis function, but only on the anisotropy encoded in the density function of its tangent fields.

Our paper is organized as follows. In Section 2, we first recall some basic facts about the Riesz transform and its use for defining anisotropic features of an image. Then in Section 3, we define the Riesz-based notion of the orientation and the structure tensor in the Gaussian selfsimilar with stationary increments case. We prove in Theorem 1 that these two characteristics are intrinsic in the sense that they depend only on the anisotropic properties of the analyzed random field. Section 4 is then devoted to the extension of all these notions to the localizable case. We then provide two classes of Gaussian models with prescribed orientation. For sake of clarity, we postponed all proofs in Section 5.

Classical tools in directionality analysis of images

In this section, we give some background about two classical tools for analyzing the anisotropy properties of an oriented texture: the local orientation and the structure tensor. We first recall in Section 2.1 the usual definitions based on the Riesz transform introduced in [START_REF] Felsberg | The monogenic signal[END_REF].

Local orientation of an image and structure tensor

The classical notion of local orientation of a texture is based on the Riesz transform. The Riesz transform Rf of any f ∈ L 2 (R 2 ) is defined in the Fourier domain 1 as

Rf = R 1 f R 2 f with R 1 f (ξ) = -j ξ 1 ξ f (ξ), R 2 f (ξ) = -j ξ 2 ξ f (ξ) , ∀ξ = (ξ 1 , ξ 2 ) ∈ R 2 .
The main properties of R [START_REF] Stein | Singular Integrals Differentiability Properties of Functions[END_REF][START_REF] Unser | Multiresolution Monogenic Signal Analysis Using the Riesz-Laplace Wavelet Transform[END_REF] are summarized in the two following propositions. The first ones concern the invariance with respect to dilations, translations, and the steerability property (relation with the rotations).

Proposition 1. The Riesz transform commutes both with the translation, and the dilation operator, that is for any f ∈ L 2 (R 2 ), a > 0 and b ∈ R 2 , one has

RD a f = D a Rf with D a f = f (a -1 •) ,
and

RT b f = T b Rf with T b f = f (• -b) .
Proposition 2. The Riesz transform is steerable, that is, for any f ∈ L 2 (R 2 ) one has

R θ (Rf ) = R -1 θ R(R θ f ) = cos θ R 1 (R θ f ) + sin θ R 2 (R θ f ) -sin θ R 1 (R θ f ) + cos θ R 2 (R θ f ) ,
where R θ f = f (R -θ •) is the rotation operator by the angle θ, and R -θ = cos θ sin θ -sin θ cos θ .

is the matrix of the spatial rotation of angle -θ.

The Riesz transform is also a unitary and componentwise antisymmetric operator on L 2 (R 2 ). Proposition 3. For any i ∈ {1, 2}, the i-th component of the Riesz transform R i is an antisymmetric operator, namely for all f, g ∈ L 2 (R 2 ) we have

R i f, g L 2 (R 2 ) = -f, R i g L 2 (R 2 ) . and R 1 2 + R 2 2 = -I, which imply in particular that Rf, Rg L 2 (R 2 ) = R 1 f, R 1 g L 2 (R 2 ) + R 2 f, R 2 g L 2 (R 2 ) = f, g L 2 (R 2 ) .
Using the Riesz transform one can also define the so-called orientation of an image [START_REF] Felsberg | The monogenic signal[END_REF][START_REF] Yang | Phase Derivative of Monogenic Signals in Higher Dimensional Spaces[END_REF].

Definition 1 (Local orientation of a deterministic function). Let f ∈ L 2 (R 2 , R * + ). Assume that a.e. Rf = 0. One then defines ∀a.e. x ∈ R 2 , n(x) = Rf (x) Rf (x) ,
the local orientation of f at point x.

The Riesz transform can be viewed as a smooth version of the gradient operator. To characterize the degree of directionality of f at some point x, a classical tool, widely used in the image processing community, is the so-called structure tensor involving the gradient [START_REF] Bigun | Optimal orientation detection of linear symmetry[END_REF][START_REF] Harris | A combined corner and edge detector[END_REF][START_REF] Jahne | Practical handbook on image processing for scientific and technical applications[END_REF], which has been revisited through the Riesz transform by [START_REF] Felsberg | Low-Level Image Processing with the Structure Multivector[END_REF]. Let us begin with the 2 × 2 matrix:

J f (x) = Rf (x)Rf (x) T = R 1 f (x) 2 R 1 f (x)R 2 f (x) R 2 f (x)R 1 f (x) R 2 f (x) 2 . ( 1
)
This matrix is symmetric, positive definite, of rank one and then admits 0 as eigenvalue and also 1 with the associated eigenvector Rf (x), that is the local orientation, since

J f (x)Rf (x) = (Rf (x)Rf (x) T )Rf (x) = Rf (x)(Rf (x) T Rf (x)) = Rf (x) .
In practice, we define the local orientation vector on a point x ∈ R 2 as the vector minimizing the distance to the set of the Riesz vectors Rf (x ′ ) on a neighborhood of x ∈ R 2 , defined by a positive window W (xx ′ ), that is we are looking for the vector n maximizing the quantity

W (x -x ′ )(Rf (x ′ ) T n) 2 dx ′ or equivalently n T J W f (x) n, with J W f (x) = (W * J f )(x) , (2)
which is the matrix J f (x) filtered by the positive windows function W to form the structure tensor [START_REF] Bigun | Optimal orientation detection of linear symmetry[END_REF][START_REF] Harris | A combined corner and edge detector[END_REF][START_REF] Jahne | Practical handbook on image processing for scientific and technical applications[END_REF]. This matrix is symmetric and positive definite, then admits two nonnegative eigenvalues λ 1 (x) and λ 2 (x). It is easy to see that the local orientation n(x) maximizing n T J W f (x) n is always an eigenvector of the matrix J W f (x) associated to its largest eigenvalue λ 1 (x). The following coherency index provides a degree of directionality at any point [START_REF] Bigun | Optimal orientation detection of linear symmetry[END_REF]:

χ f (x) = λ 1 (x) -λ 2 (x) λ 1 (x) + λ 2 (x) ∈ [0, 1] .
The case χ f (x) ≈ 1 corresponds to an almost one dimensional image at x, whereas the case χ f (x) = 0 may correspond to different situations, such as isotropy or existence of a corner. Note that, all these quantities depend on the chosen windows W . The aim of next sections is to adapt this framework to the case of random Gaussian fields. We first shall consider the case of self-similar Gaussian fields admitting stationary increments in Section 3, and in Section 4 to the more general classe of localizable Gaussian fields.

Riesz-based orientation of self-similar Gaussian fields with stationary increments

We observe first that all the definitions of the previous section cannot be extended directly to the case of random fields. Indeed, there is a difficulty to circumvent, which is the definition of the Riesz transform of a Gaussian random field, like a Brownian motion for example. The usual definition of the Riesz transform holds for L p functions. Unfortunately, the sample paths of many classical random Gaussian fields do not belong to these spaces. Defining the Riesz transform of a Gaussian field by duality (as it is often done for the transforms of distributions) is not straightforward, because the Schwartz class is not stable by the Riesz transform. Although some work on the Hilbert transform of temperate distributions [START_REF] Ishikawa | Generalized Hilbert transforms in tempered distributions[END_REF][START_REF] Koizumi | On the Hilbert transform I*[END_REF] can potentially be extended to a generalized Riesz transform adapted to random fields, we favor here to define a notion of orientation for Gaussian fields, which is easily interpretable (intuitively related to the spectral density of the field) and easy to manipulate from a computational point of view (which can be explicitly computed for usual random fields and estimated by the use of wavelets or filter bank).

In this section, we first begin with the simple case of self-similar Gaussian fields. We give in Section 3.1 some background on such fields. Then in Section 3.2, we define our notion of orientation, based on a Riesz analysis of these fields with an appropriated analysis function. Finally in Section 3.3, we give several examples of oriented self-similar Gaussian fields with stationary increments. From now on, we restrict ourselves to bidimensionnal centered real valued Gaussian fields, since our goal is to analyse anisotropic images. We shall also assume that the Gaussian field X under consideration is stochastically continuous; that is, the covariance

(x, y) → E[X(x)X(y)] ,
is a continuous function on (R 2 ) 2 .

Self-similar Gaussian fields with stationary increments

In what follows we shall focus on the special case of H-self-similar Gaussian fields admitting stationary increments (H-sssi), studied in [START_REF] Dobrushin | Gaussian and their subordinated self-similar random generalized fields[END_REF][START_REF] Benassi | Elliptic Gaussian random processes[END_REF][START_REF] Bonami | Anisotropic analysis of some Gaussian models[END_REF]. Remember that the bidimensionnal Gaussian field X is said to admit stationary increments if for any x ∈ R 2 ,

{X(x + h) -X(x)} h∈R 2 (f dd) = {X(h)} h∈R 2 ,
whereas X is said to be H-self-similar (see [START_REF] Dobrushin | Gaussian and their subordinated self-similar random generalized fields[END_REF]), for some H ∈ (0, 1) if

∀c > 0, {X(cx)} x∈R 2 (f dd) = {c H X(x)} x∈R 2 ,
where as usual

(f dd)
= means equality of finite dimensional distributions. Since X is assumed to be stochastically continuous, the self-similarity implies in particular that X(0) = 0 a.s.

We now recall, following [START_REF] Yaglom | Some classes of random fields in n-dimensional space, related to stationary random processes[END_REF], the notion of spectral measure of a Gaussian field admitting stationary increments, based on the following classical result. Proposition 4. Let X = {X(x)} x∈R 2 be a centered real-valued Gaussian field with stationary increments. Then, there exists a unique Borel measure σ X satisfying

R 2 min(1, ξ 2 ) dσ X (ξ) < ∞ ,
such that for any x, y ∈ R 2 , the covariance reads:

E(X(x)X(y)) = R 2
(e j x, ξ -1)(e -j y, ξ -1) dσ X (ξ) .

The measure σ X is called the spectral measure of the Gaussian field with stationary increments X.

In what follows we shall consider only Gaussian fields whose spectral measure σ X admits a density f X , called the spectral density of X, with respect to the Lebesgue measure: dσ X (ξ) = f X (ξ) dξ. Since X is real-valued this function is necessarily even. Such field admits an harmonizable representation:

X(x) = R 2 e j x, ξ -1 f 1/2 X (ξ) W(dξ) , (3)
where W is a complex-valued white noise. By uniqueness of the spectral density, the representation of H-self-similar Gaussian fields follows (see also [START_REF] Dobrushin | Gaussian and their subordinated self-similar random generalized fields[END_REF][START_REF] Benassi | Elliptic Gaussian random processes[END_REF][START_REF] Bonami | Anisotropic analysis of some Gaussian models[END_REF][START_REF] Cohen | Fractional fields and applications[END_REF]):

Proposition 5 (Dobrushin [START_REF] Dobrushin | Gaussian and their subordinated self-similar random generalized fields[END_REF]). Let H ∈ (0, 1) and X be a H-self-similar Gaussian field with stationary increments admitting a spectral density f X . Then f X is of the form

f X (ξ) = ξ -2H-2 C X ξ ξ ,
where C X is a positive homogeneous function defined on the sphere S 1 = {ξ ∈ R 2 , ξ = 1}. The function C X is called the anisotropy function of X.

Remark 1. The estimation problem of the anisotropy function has been addressed by Istas in [START_REF] Istas | Identifying the anisotropical function of a d-dimensional Gaussian self-similar process with stationary increments[END_REF].

We now investigate the orientation properties of a self-similar Gaussian field deformed by a linear transform.

Proposition 6. Let X be a H-self-similar Gaussian field with stationary increments admitting as spectral density f X and as anisotropy function C X . Let L be an invertible 2 × 2 real valued matrix. Define X L by X L (x) = X(L -1 x). Then X L is a H-self-similar Gaussian field admitting as

• spectral density f XL (ξ) = det(L) f X (L T ξ) , ξ ∈ R 2 ,
• anisotropy function

C XL (Θ) = det(L) L T Θ 2H+2 C X L T Θ L T Θ , Θ ∈ S 1 .
Proof. The self-similarity and stationarity properties of X L directly come from that of X and of the linearity of L. To compute the spectral density of X L , observe that:

E(X L (x), X L (y)) = Cov(X(L -1 x), X(L -1 y)) , = R 2 e j L -1 x, ξ -1 e -j L -1 y, ξ -1 f X (ξ) dξ , = R 2 e j x, (L -1 ) T ξ -1 e -j y, (L -1 ) T ξ -1 f X (ξ) dξ , = R 2 e j x, ζ -1 e -j y, ζ -1 f X (L T ζ) |det L| dζ ,
using the change of variable ζ = (L -1 ) T ξ, which directly leads by identification to the explicit expression of f XL , as well as that of C XL given in Proposition 6.

We now explain how to define in a proper way the directional characteristics of a Gaussian field X admitting stationary increments in the self-similar case, and the relation to its anisotropy function C X .

Definition of a Riesz-based orientation in the self-similar case

Our notion of Riesz-based orientation of a self-similar Gaussian field will be based on the following preliminary result, leading to a new formulation for the structure tensor.

Theorem 1. Let X be a H-self-similar Gaussian field admitting a spectral density f X whose anisotropy function C X is bounded and ψ a real isotropic wavelet with 2 zero moments and fast decay, defined by its 2-D Fourier transform ψ(ξ) = ϕ( ξ ). By analogy with the structure tensor (1)-( 2) evaluated at the origin, let us define the matrix:

J ψ X = | X, R 1 ψ | 2 X, R 1 ψ X, R 2 ψ X, R 1 ψ X, R 2 ψ | X, R 2 ψ | 2 , with the Gaussian variables X, R ℓ ψ = R 2 X(x)R ℓ ψ(x) dx.
Then, the covariance matrix of the Gaussian vector

( X, R 1 ψ , X, R 2 ψ ) T is E J ψ X = +∞ 0 |ϕ(r)| 2 r 2H+1 dr J X ,
where for any

ℓ 1 , ℓ 2 ∈ ({1, 2}) 2 , [J X ] ℓ1,ℓ2 = Θ∈S 1 Θ ℓ1 Θ ℓ2 C X (Θ) dΘ , (4)
with the notation Θ = (Θ 1 , Θ 2 ). J X is a nonnegative definite 2 × 2 matrix depending only on the anisotropy function C X , and will be called the structure tensor of X.

Proof. The assumptions on ψ insure that the function g ℓ defined by g ℓ (ξ) = R ℓ ψ(ξ)f

1/2
X (ξ) is square-integrable; that is, the stochastic integral R 2 g ℓ (ξ) W(dξ) is well defined and the following operations lead to the a.s. existence of the Gaussian variables X, R ℓ ψ , for ℓ = 1, 2:

R 2 R ℓ ψ(ξ)f 1/2 X (ξ) W(dξ) = R 2 R ℓ ψ(ξ) -R ℓ ψ(0) f 1/2 X (ξ) W(dξ) , = R 2 R 2 e j x, ξ R ℓ ψ(x) dx - R 2 R ℓ ψ(x) dx f 1/2 X (ξ) W(dξ) , = R 2 R 2 e j x, ξ -1 f 1/2 X (ξ) W(dξ) R ℓ ψ(x) dx , = X, R ℓ ψ .
The first equality follows since R ℓ ψ(0) = 0 by the moment assumption on ψ and the third equality comes from the classical stochastic Fubini Theorem (see [START_REF] Pipiras | Regularization and integral representations of Hermite processes[END_REF]), which holds since

R 2 R 2 |e j x, ξ -1| 2 f X (ξ)|ψ(x)| 2 dξ 1/2 dx = R 2 R 2 |e j x, ξ -1| 2 f X (ξ) dξ 1/2 |ψ(x)| dx < ∞ ,
by the integrability properties of ψ and the existence of the stochastic integral defining X. The covariance matrix is then easily computed: by definition of the Riesz transform, one has for any

ℓ 1 , ℓ 2 ∈ ({1, 2}) 2 , E[ X, R ℓ1 ψ X, R ℓ2 ψ ] = R 2 ξ ℓ1 ξ ℓ2 ξ 2 ψ(ξ) 2 f X (ξ) dξ .
We set ξ = rΘ with (r, Θ) ∈ R + × S 1 . Hence, remembering that ψ is isotropic, that is ψ(ξ) = ϕ(r) and using the expression of f X given in Proposition 5, we obtain

E[ X, R ℓ1 ψ X, R ℓ2 ψ ] = Θ∈S 1 +∞ r=0 1 r 2H+2 Θ ℓ1 Θ ℓ2 ϕ(r) 2 C X (Θ) r dr dΘ , = +∞ r=0 |ϕ(r)| 2 r 2H+1 dr Θ∈S 1 Θ ℓ1 Θ ℓ2 C X (Θ) dΘ .
Theorem 1 then follows.

We now define the structure tensor of X, its orientation and its coherency index.

Definition 2 (Orientation and coherency). The matrix J X defined in Theorem 1 ( 4) is called the structure tensor of X. Let λ 1 , λ 2 be its two eigenvalues. The coherency index of X is defined as

χ(X) = |λ 2 -λ 1 | λ 1 + λ 2 .
An orientation n is any unit eigenvector associated to the largest eigenvalue of J X .

Examples

We present below several examples of self-similar Gaussian fields, and every time we make explicit their structure tensor and their intrinsic orientation.

Example 1: Fractional Brownian Field (FBF)

We begin with the Fractional Brownian Field (FBF) defined in [START_REF] Lindstrøm | Fractional Brownian fields as integrals of white noise[END_REF]. This random field is the isotropic multidimensional extension of the famous Fractional Brownian Motion defined by [START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF]. Its harmonizable representation is:

X(x) = R 2 1 ξ
H+1 (e j x, ξ -1)

W(dξ) 2π . We can easily check (C X ≡ 1 2π ), setting Θ = (cos θ, sin θ) ∈ S 1 , that [J X ] 1,1 = 1 2π 2π 0 cos 2 θ dθ = 1 2 , [J X ] 1,2 = [J X ] 2,1 = 1 2π 2π 0 cos θ sin θ dθ = 0 , and
[J X ] 2,2 = 1 2π 2π 0 sin 2 θ dθ = 1 2 ,
which directly implies

J X = 1 2 0 0 1 2 , χ(X) = 0 .
Any unit vector is thus an orientation of the FBF, which is clearly consistent with its isotropic nature.

Example 2: Elementary fields (EF)

We focus here to a special case of H-sssi fields called elementary fields (EF), whose anisotropy function is a cone with an orientation α 0 and a half-width δ > 0. It admits the following harmonizable representation (α 0 ∈ (-π/2, π/2), δ > 0):

X α0,δ (x) = R 2 (e j x, ξ -1)f 1/2 X (ξ) W(dξ) , (5)
where

f X (ξ) = ξ -2H-2 C X ξ ξ , (6)
C X (Θ) = 1 4δ 1 [α0-δ,α0+δ] (arg Θ) + 1 [α0+π-δ,α0+π+δ] (arg Θ) . (7)
Let us compute its structure tensor J X , using the definition given in Theorem 1.

Remark 2. The contribution of each of the portions of the cone is identical in the integrals defining the coefficients of J X . Thus, in the following computations we only consider one portion whose contribution is doubled.

We start with the diagonal terms:

[J X ] 1,1 = Θ∈S 1 Θ 2 1 C X (Θ) dΘ = 1 2δ α0+δ α0-δ cos 2 θ dθ = 1 2 + 1 2 cos(2α 0 ) sin(2δ) 2δ .
By the relation cos 2 θ + sin 2 θ = 1, we get as well

[J X ] 2,2 = 1 2 - 1 2 cos(2α 0 ) sin(2δ) 2δ .
The last terms are computed as follows:

[J X ] 1,2 = [J X ] 2,1 = 1 2δ α0+δ α0-δ cos θ sin θ dθ = 1 2 sin(2α 0 ) sin(2δ) 2δ .
Hence the structure tensor of the elementary field is

J X = 1 2 + 1 2 cos(2α 0 ) sin(2δ) 2δ 1 2 sin(2α 0 ) sin(2δ) 2δ 1 2 sin(2α 0 ) sin(2δ) 2δ 1 2 -1 2 cos(2α 0 ) sin(2δ) 2δ .
Remark that J X diagonalizes as

J X = cos α 0 -sin α 0 sin α 0 cos α 0 1 2 + 1 2 sin(2δ) 2δ 0 0 1 2 -1 2 sin(2δ) 2δ cos α 0 -sin α 0 sin α 0 cos α 0 T .
Denoting λ 1 λ 2 the two eigenvalues of J X ,

λ 1 = 1 2 + 1 2 sin(2δ) 2δ and λ 2 = 1 2 - 1 2 sin(2δ) 2δ ,
the coherency index of X is given by

χ(X) = λ 1 -λ 2 λ 1 + λ 2 = sin(2δ) 2δ .
An orientation of the elementary field X α0,δ being a unit eigenvector associated with λ 1 , we obtain

n = cos α 0 sin α 0 .
Remark 3. This is in accordance with what we observe when performing simulations. For δ = π/2 we recover the classical FBF which is isotropic and we get χ(X) = sin(π)/π = 0. Conversely, notice that χ(X) tends to 1 when δ → 0, meaning that the coherency is strong when the cone of admissible directions is tight around the angle α 0 . Note that in the limite case δ → 0, the density function C X tends to the Dirac measure along the line arg Θ = α 0 , and the structure tensor degenerates to

J X = cos 2 α 0 cos α 0 sin α 0 cos α 0 sin α 0 sin 2 α 0 ,
which diagonalizes as

J X = R α0 1 0 0 0 R T α0 .
leading to the same orientation vector (cos α 0 , sin α 0 ) T . Notice also that in the limit case, the structure tensor is not invertible anymore.

Example 3: sum of two elementary fields

To understand how our notion of Riesz-based orientation can be adapted to the setting of multiple oriented random fields, we consider the following toy model (α 0 , α 1 ∈ (-π/2, π/2), δ > 0):

X = X α0,δ + X α1,δ .
The Gaussian field X is then the sum of two elementary fields of same regularity H and of respective directions α 0 = α 1 (as defined in Example 2 above). We assume that δ < |α 1 -α 0 |/2. This last condition implies in particular that [α 0 -δ, α 0 + δ] ∩ [α 1 -δ, α 1 + δ] = ∅ and the spectral densities have disjoint supports, that is X α0,δ and X α1,δ are supposed to be independent.

Then we have

J X = 1 + 1 2 sin(2δ)
2δ (cos(2α 0 ) + cos(2α 1 )) .

As previously the matrix J X diagonalizes as

J X = R (α0+α1)/2 1 + sin(2δ) 2δ cos(α 0 -α 1 ) 0 0 1 -sin(2δ) 2δ cos(α 0 -α 1 ) R T (α0+α1)/2 ,
where we denoted

R (α0+α1)/2 = cos α0+α1 2 -sin α0+α1 2 sin α0+α1 2 cos α0+α1 2 .
Thus, the coherency index is

χ(X) = sin(2δ) 2δ cos(α 0 -α 1 ) ,
which tends to cos(α 0 -α 1 ) when δ → 0. One can also observe that the closer α 0 and α 1 are, the more coherent the random field is and then admits a dominant orientation. An orientation of X is given by

n = cos α0+α1 2 sin α0+α1 2 .
We then recover a dominant orientation, related to half the angle of the two orientations.

Example 4: deformation of an elementary field by a linear transform

Let L be an invertible 2 × 2 matrix and α 0 ∈ (-π/2, π/2), δ > 0. Set

X L (x) = X α0,δ (L -1 x) . (8)
Thanks to Proposition 6, we have an explicit expression for the spectral density of X L :

f XL (ξ) = |det(L)| 2δ L T ξ -2H-2 1 [α0-δ,α0+δ] (arg(L T ξ)) , = |det(L)| 2δ L T ξ -2H-2 1 [tan(α0-δ),tan(α0+δ)] ((L T ξ) 2 /(L T ξ) 1 ) .
Since the matrix L is invertible, it admits a Singular Value Decomposition L = U∆V T where U, V are two orthogonal matrices and ∆ a diagonal matrix with nonnegative eigenvalues. One can then deduce the general case of an invertible matrix L from three specific ones:

L ∈ O + 2 (R), L ∈ O - 2 (
R) and L diagonal with nonnegative eigenvalues. Before deriving the general form of an orientation vector, we will consider each term of the SVD.

(i) We first consider the case where L is an orthogonal matrix of the form

L = R θ0 = cos θ 0 -sin θ 0 sin θ 0 cos θ 0 , one has f XL (ξ) = 1 2δ ξ -2H-2 1 [α0+θ0-δ,α0+θ0+δ] (arg ξ) ,
which implies that one can choose as orientation for X L the unit vector

n L = cos (α 0 + θ 0 ) sin (α 0 + θ 0 ) = L cos α 0 sin α 0 = (L -1 ) T cos α 0 sin α 0 , (9)
since any orthogonal matrix equals the transpose of its inverse.

(ii) We now deal with the case where L is an orthogonal matrix of the form

L = cos θ 0 sin θ 0 sin θ 0 -cos θ 0 = cos θ 0 -sin θ 0 sin θ 0 cos θ 0 × 1 0 0 -1 . One has f XL (ξ) = 1 2δ ξ -2H-2 1 [θ0-α0-δ,θ0-α0+δ] (arg ξ) ,
which implies that one can choose as orientation for X L the unit vector

n L = cos (θ 0 -α 0 ) sin (θ 0 -α 0 ) = L cos α 0 sin α 0 = (L -1 ) T cos α 0 sin α 0 , (10)
since as above any orthogonal matrix equals the transpose of its inverse.

(iii) We finally deal with the case where L is a diagonal matrix L = λ 1 0 0 λ 2 , with λ 1 , λ 2 > 0.

In this case observe that the condition

tan(α 0 -δ) < (L T ξ) 2 (L T ξ) 1 < tan(α 0 + δ) , is equivalent to λ1 λ2 tan(α 0 -δ) < ξ2 ξ1 < λ1 λ2 tan(α 0 + δ), that is to δ ∆ < arg ξ < δ ∆ ,
with δ ∆ = arctan( λ1 λ2 tan(α 0 -δ)) and δ ∆ = arctan( λ1 λ2 tan(α 0 + δ)). Hence,

f XL (ξ) = |det(L)| 2δ L T ξ -2H-2 1 [δ ∆ ,δ∆] (arg ξ) , C XL (Θ) = |det(L)| 2δ L T Θ -2H-2 1 [δ ∆ ,δ∆] (arg Θ) .
Now, recalling that:

[J XL ] ℓ1,ℓ2 = Θ∈S 1 Θ ℓ1 Θ ℓ2 C XL (Θ) dΘ ,
we obtain

[J XL ] 1,1 |det(L)| = 1 2δ δ∆ δ ∆ cos 2 θ C XL (cos θ, sin θ) dθ = 1 2δ δ∆ δ ∆ cos 2 θ (λ 2 1 cos 2 θ + λ 2 2 sin 2 θ) H+1 dθ , [J XL ] 2,2 |det(L)| = 1 2δ δ∆ δ ∆ sin 2 θ C XL (cos θ, sin θ) dθ = 1 2δ δ∆ δ ∆ sin 2 θ (λ 2 1 cos 2 θ + λ 2 2 sin 2 θ) H+1 dθ , [J XL ] 1,2 |det(L)| = 1 2δ δ∆ δ ∆ cos θ sin θ C XL (cos θ, sin θ) dθ = 1 2δ δ∆ δ ∆ cos θ sin θ (λ 2 1 cos 2 θ + λ 2 2 sin 2 θ) H+1 dθ .
Now, let us define (u 1 (θ), u 2 (θ)) = (cos θ, sin θ) and let us introduce the functions

f ℓ1,ℓ2 : θ → u ℓ1 (θ)u ℓ2 (θ) (λ 2 1 u 1 (θ) 2 + λ 2 2 u 2 (θ) 2 ) -H-1 , F ℓ1,ℓ2 : x → x 0 f ℓ1,ℓ2 (θ) dθ, ν : α → arctan λ 1 λ 2 tan α , G ℓ1,ℓ2 : α → F ℓ1,ℓ2 (ν(α)) .
Each term of the structure tensor writes

[J XL ] ℓ1,ℓ2 |det(L)| = G ℓ1,ℓ2 (α 0 + δ) -G ℓ1,ℓ2 (α 0 -δ) 2δ .
When the parameter δ is small, we have

G ℓ1,ℓ2 (α 0 + δ) -G ℓ1,ℓ2 (α 0 -δ) 2δ = G ′ ℓ1,ℓ2 (α 0 ) + δ 2 12 G ′′′ ℓ1,ℓ2 (α 0 ) + O(δ 4 ) , = ν ′ (α 0 ) F ′ ℓ1,ℓ2 (ν(α 0 )) + O(δ 2 ) , = ν ′ (α 0 ) f ℓ1,ℓ2 (ν(α 0 )) + O(δ 2 ) . Hence, G ℓ1,ℓ2 (α 0 + δ) -G ℓ1,ℓ2 (α 0 -δ) 2δ = λ 1 λ 2 λ 2 2 cos 2 α 0 + λ 2 1 sin 2 α 0 × u ℓ1 (ν(α 0 )) u ℓ2 (ν(α 0 )) (λ 2 1 cos 2 (ν(α 0 )) + λ 2 2 sin 2 (ν(α 0 ))) H+1 +O(δ 2 ) . Let us define C H,λ1,λ2,α0 = λ 2 1 λ 2 2 (λ 2 2 cos 2 α 0 + λ 2 1 sin 2 α 0 ) -1 (λ 2 1 cos 2 (ν(α 0 )) + λ 2 2 sin 2 (ν(α 0 ))) -H-1 .
Then, one has for small δ,

J XL = C H,λ1,λ2,α0 cos 2 (ν(α 0 )) cos(ν(α 0 )) sin(ν(α 0 )) cos(ν(α 0 )) sin(ν(α 0 )) sin 2 (ν(α 0 )) + O(δ 2 ) ,
which can be written as

J XL = C H,λ1,λ2,α0 R ν(α0) 1 0 0 0 R T ν(α0) + O(δ 2 ) .
Therefore n L = (cos ν(α 0 ), sin ν(α 0 )) T can be viewed as an approximate eigenvector of J XL associated to its largest eigenvalue, and then an orientation of X L . Finally, we remark that

ν(α 0 ) = arctan λ 1 sin α 0 λ 2 cos α 0 = arg λ 2 cos α 0 λ 1 sin α 0 = arg λ 2 0 0 λ 1 × cos α 0 sin α 0 .
Consequently, an approximate (up to δ 2 ) orientation of X L is in this case

n L = λ 2 0 0 λ 1 n λ 2 0 0 λ 1 n with n = cos α 0 sin α 0 .
Observe that λ 2 0 0 λ 1 is the adjugate matrix of L. Then, dividing the numerator and denominator of the last equation by det(L) = λ 1 λ 2 , we get

n L = L -1 n L -1 n = (L -1 ) T n (L -1 ) T n , (11)
since the diagonal matrix L -1 equals its transpose. We now gather ( 9), ( 10) and ( 11): using the existence of the SVD for every matrix, we deduce the following proposition.

Proposition 7. Let L be an invertible 2 × 2 matrix and X L the Gaussian field defined by [START_REF] Davies | Fractal analysis of surface roughness by using spatial data (with discussion)[END_REF]. Set n = cos α0 sin α0 the orientation vector of X. Then the unit vector

n L = (L -1 ) T n (L -1 ) T n ,
is an approximate (up to δ 2 ) orientation vector of X L .

Proof. Let L = U∆V T be the SVD singular decomposition of L, with U, V ∈ O 2 (R) and ∆ diagonal with nonnegative eigenvalues.

X L (x) = X(V -T ∆ -1 U -1
x) be the Gaussian field defined by [START_REF] Davies | Fractal analysis of surface roughness by using spatial data (with discussion)[END_REF]. Let decompose the three operations as follows :

X L = X • (V T ) -1 X1 •∆ -1 X2 •U -1 .
Then, since V T is an orthogonal matrix, we have from (i) and (ii) that the unit orientation vector of X 1 is

n 1 = ((V T ) -1 ) T n . Now from (iii), the unit orientation vector (up to δ 2 ) of X 2 = X 1 • ∆ -1 is n 2 = (∆ -1 ) T n 1 (∆ -1 ) T n 1 = (∆ -1 ) T ((V T ) -1 ) T n (∆ -1 ) T ((V T ) -1 ) T n .
Finally, from (i) and (ii) again, the unit orientation vector (up to δ 2 ) of

X L = X 2 • U -1 is n L = (U -1 ) T n 2 (U -1 ) T n 2 , = (U -1 ) T (∆ -1 ) T ((V T ) -1 ) T n (U -1 ) T (∆ -1 ) T ((V T ) -1 ) T n , = ((V T ) -1 ∆ -1 U -1 ) T n ((V T ) -1 ∆ -1 U -1 ) T n = (U∆V T ) T n (U∆V T ) T n = (L -1 ) T n (L -1 ) T n .

Riesz-based orientation of localizable Gaussian fields

We now extend the notion of intrinsic orientation, defined for self-similar random fields, to a much more general setting, that of localizable Gaussian fields. This will be the purpose of Section 4.2. In Section 4.3, we will apply it to two classes of model with prescribed orientation. First of all, let us recall in Section 4.1 the definition of localizable Gaussian fields.

Localizable Gaussian fields

We first recall, following [START_REF] Benassi | Elliptic Gaussian random processes[END_REF][START_REF] Falconer | Tangent Fields and the local structure of random fields[END_REF][START_REF] Falconer | The local structure of random processes[END_REF], the definition of H(x 0 )-localizable Gaussian fields.

Definition 3 (Localizable Gaussian field). We say that the random field

Y = {Y (x), x ∈ R 2 } is H(x 0 )-localizable, at x 0 ∈ R 2 with H(x 0 ) ∈ (0, 1)
and with tangent field (or local form) the non-trivial random field

Y x0 = {Y x0 (h), h ∈ R 2 } if Y (x 0 + ρh) -Y (x 0 ) ρ H h∈R 2 d -→ Y x0 (h) h∈R 2 , (12)
as ρ → 0, where d → means convergence in distribution, that is the weak convergence for stochastic processes (see [START_REF] Billingsley | Convergence of Probability Measures[END_REF]). A random field Y = {Y (x), x ∈ R 2 } is said to be localizable if for all x ∈ R 2 it is H(x)-localizable for some H(x) ∈ (0, 1).

In Theorem 3.9 and Corollary 3.10 of [START_REF] Falconer | Tangent Fields and the local structure of random fields[END_REF], Falconer proved the following result that we state in the Gaussian case. It enables to describe the whole class of possible tangent fields of a Gaussian field with continuous sample paths.

Theorem 2. Let X be a localizable Gaussian field with continuous sample paths. For almost all x 0 in R 2 the tangent field Y x0 of X at x 0 has stationary increments and is self-similar, that is for some H(x 0 ) ∈ (0, 1) and for all ρ 0,

{Y x0 (ρx), x ∈ R 2 } (f dd) = {ρ H(x0) Y x0 (x), x ∈ R 2 } . (13)
In short, a Gaussian field with continuous sample paths will have at a.e. point, a "fractal" tangent field behaving like a FBF. We now illustrate this notion considering a classical example of Gaussian field with prescribed tangent field: the Multifractional Brownian Field defined in the unidimensional setting in [START_REF] Lévy-Vehel | Multifractional Brownian Motion : definition and preliminary results[END_REF], and in the multivariate case in [START_REF] Benassi | Elliptic Gaussian random processes[END_REF][START_REF] Herbin | From N-parameter fractional Brownian motions to N-parameter multifractional Brownian motions[END_REF]. Such field is localizable at each point, with a fractional Brownian Field for tangent field.

Example 1 (Multifractional Brownian Field). Let h : R 2 → (0, 1) be a continuously differentiable function whose range is supposed to be a compact interval [α, β] ⊂ (0, 1). The Multifractional Brownian Field (MBF) with multifractional function h, is the Gaussian field defined by its harmonisable representation as follows

X h (x) = R 2 e j x, ξ -1 ξ h(x)+1 W(dξ) . (14)

Tensor structure and orientation of localizable Gaussian fields

The results of Section 3 together with Theorem 2 of section 4.1, will allow us to define the Riesz-based orientation of any localizable Gaussian field X almost everywhere. Definition 4 (Localizable field orientations). Let X be a Gaussian field with continuous sample paths. Assume that X is localizable at the point x 0 , with tangent field Y x0 , and that Y x0 is a self-similar Gaussian field with stationary increments. One then defines:

• The local anisotropy function C x0 at x 0 of the localizable Gaussian field X is the anisotropy function of its tangent field Y x0 .

• The local structure tensor J x0 at x 0 of the localizable Gaussian field X is the structure tensor of its tangent field Y x0 .

• A local orientation at x 0 of the localizable Gaussian field X is any orientation of its tangent field Y x0 .

In view of these definitions and of Theorem 2, we deduce that any localizable Gaussian field X admits a local orientation at almost every point

x 0 ∈ R 2 .
Example 2 (Local structure tensor and orientation of a MBF). The Multifractional Brownian Field X h (14) admits at each point a structure tensor proportional to the identity matrix. In particular, any unit vector is an orientation of X h . Indeed, the tangent field Y x0 of the MBF at point x 0 is a FBF of Hurst index h(x 0 ), whose structure tensor has been determined in example 3.3.1.

Two new models of localizable Gaussian fields with prescribed orientation

In this section, we will extend our previous works [START_REF] Polisano | Texture modeling by Gaussian fields with prescribed local orientation[END_REF][START_REF] Polisano | Modélisation de texture par champ gaussien à orientation locale prescrite[END_REF] and define two classes of Gaussian fields with prescribed orientation. The details about numerical aspects and synthesis of the model, as well as comparison between them, will be detailed in the companion paper [START_REF] Polisano | Simulation of oriented pattern with prescribed local orientation[END_REF]. These two models will be derived from two general classes: Generalized Anisotropic Fractional Brownian Fields (GAFBF) and Warped Anisotropic Fractional Brownian Fields (WAFBF) that we describe in Sections 4.3.1 and 4.3.2 respectively.

First model: Generalized Anisotropic Fractional Brownian Fields (GAFBF)

We introduce below the definition of Generalized Anisotropic Fractional Brownian Fields (GAFBF) which generalizes the notion of Locally Anisotropic Brownian Fields (LAFBF) introduced in [START_REF] Polisano | Texture modeling by Gaussian fields with prescribed local orientation[END_REF], and whose simulation will be studied in [START_REF] Polisano | Simulation of oriented pattern with prescribed local orientation[END_REF].

Our Gaussian field will be defined from two functions h from R 2 to [0, 1] and C from R 2 × R 2 to R + satisfying the following set of assumptions:

Assumptions (H) • h is a β-Hölder function, such that a = inf x∈R 2 h(x) > 0, b = sup x∈R 2 h(x) < 1 and b < β 1. • (x, ξ) → C(x, ξ) is bounded, that is ∀(x, ξ) ∈ R 2 × R 2 , C(x, ξ) M .
• ξ → C(x, ξ) is even and homogeneous of degree 0: ∀ρ > 0, C(x, ρξ) = C(x, ξ).

• x → C(x, ξ) is continuous and satisfies: there exists some η, with b < η 1 such that ∀x ∈ R 2 , sup z∈B(0,1)

z -2η S 1 C(x + z, Θ) -C(x, Θ) 2 dΘ A x < ∞ . (15)
Morever x → A x is bounded on any compact set of R 2 .

We now define our model, the Generalized Anisotropic Fractional Brownian Field.

Definition 5. Generalized Anisotropic Fractional Brownian Fields (GAFBF)

Let us consider h : R 2 → [0, 1] and C satisfying Assumptions (H). We then define the GAFBF as the following Gaussian field generalizing [START_REF] Polisano | Texture modeling by Gaussian fields with prescribed local orientation[END_REF][START_REF] Polisano | Modélisation de texture par champ gaussien à orientation locale prescrite[END_REF] by

X(x) def = R 2 (e j x, ξ -1) C(x, ξ) ξ h(x)+1 W(dξ) . ( 16
)
The main properties of the GAFBF X are summarized in the following propositions.

Theorem 3. The GAFBF X (16) admits at any point x 0 ∈ R 2 , a tangent field Y x0 given by

Y x0 (x) = R 2 (e j x, ξ -1)f 1/2 (x 0 , ξ) W(dξ) = R 2 (e j x, ξ -1) C(x 0 , ξ) ξ h(x0)+1 W(dξ) . (17)
In particular, for each point x 0 , the local anisotropy function of the Gaussian field X at x 0 is

C x0 : Θ → C(x 0 , Θ) 2 .
Proof. Theorem 3 is proven in Section 5.2.

Example 3. We now derive our first example of Gaussian field with prescribed orientation. The LAFBF introduced in [START_REF] Polisano | Texture modeling by Gaussian fields with prescribed local orientation[END_REF][START_REF] Polisano | Modélisation de texture par champ gaussien à orientation locale prescrite[END_REF] is a particular case of GAFBF where the function C is a localized version of the cone (7) with constant half-width δ > 0 and whose orientation may vary spatially, that is where α : R 2 → (-π/2, π/2) is now a continuously differentiable function, which is 2η-Holderian with b < η 1/2. One can verify that the so-called function C satisfies assumptions (H). Applying Theorem 3, the corresponding localizable Gaussian field defined by formula ( 16) admits a tangent field at any points x 0 ∈ R 2 and the following cone as local anisotropy:

C x0 (Θ) = 1 4δ 1 [α(x0)-δ,α(x0)+δ] (arg Θ) + 1 [α(x0)+π-δ,α(x0)+π+δ] (arg Θ) . (18)
The above example shows that the tangent field of a LAFBF is actually an elementary field. Then, using the results of Example 3.3.2 and Definition 4, we immediately deduce the following proposition:

Proposition 8. The LAFBF defined in the Example 3 admits at each point x 0 an approximate (up to δ 2 ) local orientation vector given by n = cos α(x 0 ) sin α(x 0 ) .

Warped Anisotropic Fractional Brownian Fields (WAFBF)

We now consider a second model, satisfying similar properties, in the same spirit as the approach developed in [START_REF] Perrin | Reducing non-stationary stochastic processes to stationarity by a time deformation[END_REF][START_REF] Perrin | Reducing non-stationary random fields to stationarity and isotropy using a space deformation[END_REF][START_REF] Guyon | Identification of space deformation using linear and superficial quadratic variations[END_REF] but in the case where the warped Gaussian field is a H-self-similar Gaussian field with stationary increments. Definition 6 (Warped Anisotropic Fractional Brownian Fields). Let X be a H-sssi field, with anisotropy function C X , as explicit in Proposition 5. Let Φ : R 2 → R 2 be a continuously differentiable function. The Warped Anisotropic Fractional Brownian Field (WAFBF) Z Φ,X is defined as the deformation of the elementary field X:

Z Φ,X (x) = X(Φ(x)) . ( 19
)
The aim of this section is to study the local properties of such Gaussian fields.

Proposition 9. The Gaussian field Z Φ,X defined by ( 19) is localizable at any point x 0 ∈ R 2 , with tangent field Y x0 defined as

Y x0 (x) = X(DΦ(x 0 ) x) , ∀x ∈ R 2 , ( 20
)
where DΦ(x 0 ) is the Jacobian matrix of Φ at point x 0 .

Proof. Proposition 9 is proven in Section 5.3.

Proposition 10. Let Z Φ,X be the WAFBF defined in Definition 6, from an elementary field X = X α0,δ defined by ( 5)-( 6)-( 7) and whose orientation, computed in section 3.3.2, is given by the unit vector n = (cos α 0 , sin α 0 ). In addition, we assume that the

C 1 -differentiable deformation Φ is a diffeomorphism on an open set U ⊂ R 2 .
Then, at each point x 0 ∈ U , the WAFBF Z Φ,X admits an approximate (up to δ 2 ) local orientation given by

n Z (x 0 ) = DΦ(x 0 ) T n DΦ(x 0 ) T n .
Proof. According to Definition 4, the local orientation of Z Φ,X at x 0 ∈ U is given by the one of its tangent field Y x0 . From Proposition 9, Y x0 (x) = X(DΦ(x 0 ) x), and since Φ a diffeomorphism in a neighborhood of x 0 , DΦ(x 0 ) is invertible. Proposition 7 applied to L -1 = DΦ(x 0 ) directly leads to the result.

Example 4 (Local rotation). We now illustrate this result considering the case α 0 = 0, then the orientation n = (cos α 0 , sin α 0 ) of the elementary field X is now the unit vector n = e 1 = (1, 0) T . The deformation we consider is a local rotation governed by a continuously differentiable function x → α(x). We have the following proposition:

Proposition 11. We consider the following warped field from the standard elementary field:

Z Φ,X (x) = X 0,δ (Φ(x)) , with Φ(x) = R -α(x) x = cos α(x) x 1 + sin α(x) x 2 -sin α(x) x 1 + cos α(x) x 2 ≡ Φ 1 (x) Φ 2 (x) , ( 21
)
where α : R 2 → R is a C 1 function on R 2 such that, on an open set U ⊂ R 2 , one has:

∀x 0 ∈ U, ∇α(x 0 ) ∧ x 0 = ∂α ∂x 1 (x 0 )x 0,2 - ∂α ∂x 2 (x 0 )x 0,1 = -1 . ( 22
)
Then, for each point x 0 ∈ U satisfying [START_REF] Istas | Identifying the anisotropical function of a d-dimensional Gaussian self-similar process with stationary increments[END_REF], Z Φ,X admits as local orientation vector:

n(x 0 ) = u(α(x 0 )) + u(α(x 0 )) ⊥ , x 0 ∇α(x 0 ) u(α(x 0 )) + u(α(x 0 )) ⊥ , x 0 ∇α(x 0 ) . ( 23
)
with u(α(x 0 )) = (cos(α 0 (x 0 )), sin(α 0 (x 0 )).

Proof. Since the function α is assumed to be C 1 , the deformation Φ (21) is also C 1 . Its Jacobian matrix is given by

DΦ(x) = cos α(x) + ∂ α ∂x 1 (x)Φ 2 (x) sin α(x) + ∂α ∂x2 (x)Φ 2 (x) -sin α(x) -∂α ∂x1 (x)Φ 1 (x) cos α(x) -∂α ∂x2 (x)Φ 1 (x) , whose determinant is det DΦ(x) = 1 + ∂α ∂x 1 (x)x 2 - ∂α ∂x 2 (x)x 1 .
Under the assumption [START_REF] Istas | Identifying the anisotropical function of a d-dimensional Gaussian self-similar process with stationary increments[END_REF] followed by α, the determinant on the open set U is non-zero, so Φ is a C 1 -diffeomorphism on U . Then, Proposition 9 and 10 hold, and at each point x 0 ∈ U , Z Φ,X admits as local orientation vector n(x 0 ) = DΦ(x 0 ) T e 1 / DΦ(x 0 ) T e 1 , which writes

n(x 0 ) = u(α(x 0 )) + u(α(x 0 )) ⊥ , x 0 ∇α(x 0 ) u(α(x 0 )) + u(α(x 0 )) ⊥ , x 0 ∇α(x 0 ) .
(x 1 , x 2 ) = -π 2 + x 1 (top right), (c) α(x 1 , x 2 ) = -π 2 + x 2 (bottom left), (d) α(x 1 , x 2 ) = -π 2 + x 2 1 -x 2 (bottom right).
Some examples of realizations of Z Φ,X on the domain [0, 1] 2 are displayed on Figure 1 for different values of α, fulfilling the condition [START_REF] Istas | Identifying the anisotropical function of a d-dimensional Gaussian self-similar process with stationary increments[END_REF]. Remark that the orientation vector given by ( 23) is equal to u(α(x 0 )) = (cos α(x 0 ), sin α(x 0 )) plus a term depending on the gradient of α. Consequently, we do not exactly have a prescribed orientation governed by α.

We now inverse the problem, and investigate the construction of a deformation Φ, to obtain a prescribed orientation α. To this aim, we will employ a conformal deformation, which has the particularity to preserve the angles. An important result, stated in the following proposition, is that we can prescribe the orientation α of a Gaussian field, if this orientation is supposed to be harmonic.

Proposition 12. Let Z Φ,X (x) be the Gaussian field [START_REF] Bonami | Anisotropic analysis of some Gaussian models[END_REF], warped by a conformal deformation Φ defined as follows: let α : R 2 → R be an harmonic function, and λ its harmonic conjugate function such that Ψ = λ -α is holomorphic (as a complex function, identifying R 2 with C). Define now Φ as any complex primitive of exp(Ψ), as an holomorphic function on C. Then at any point x 0 , an approximate (up to δ 2 ) local orientation of Z Φ,X is

n Z (x 0 ) = cos(α(x 0 )) sin(α(x 0 )) ,
which is exactly the orientation vector defined by the angle function α.

Proof. Firstly, the existence of λ is the classical result of the existence of an harmonic conjugate of α (see [START_REF] Stein | Singular Integrals Differentiability Properties of Functions[END_REF]). Then Ψ is holomorphic, and exp(Ψ) is holomorphic too (as the composition of holomorphic functions). In addition, since Φ is a complex primitive of exp(Ψ) as an holomorphic function on C, we have at any point:

Φ ′ (x 0 ) = exp(Ψ(x 0 )) = e λ(x0) e -iα(x0)
(as a complex function in C). Moreover, since Φ is holomorphic,

Φ ′ (x 0 ) = ∂Φ ∂x 1 (x 0 ) = -i ∂Φ ∂x 2 (x 0 ) ,
which leads to the Jacobian matrix:

DΦ(x 0 ) = exp(λ(x 0 )) cos(α(x 0 )) sin(α(x 0 )) -sin(α(x 0 )) cos(α(x 0 ))
and concludes the proof.

Example 5 (Affine orientation functions). We consider the family of harmonic functions

α(x 1 , x 2 ) = ax 1 + bx 2 + c ,
with a, b, c real constants. By the procedure of Proposition 12, we are able to construct the deformation function Φ, whose explicit formula is

Φ(x 1 , x 2 ) = exp(ax 2 -bx 1 ) a 2 + b 2 a sin(ax 1 + bx 2 + c) -b cos(ax 1 + bx 2 + c) a cos(ax 1 + bx 2 + c) + b sin(ax 1 + bx 2 + c) . (24)
Then we can verify that

DΦ(x) T e 1 = exp(ax 2 -bx 1 ) cos(ax 1 + bx 2 + c) sin(ax 1 + bx 2 + c) , n(x) = cos α(x) sin α(x) .
An example of simulation of such a prescribed local orientation field is provided in Figure 2, where the angle variations are governed by the function α(

x 1 , x 2 ) = 2x 1 -x 2 .

Proofs

This last section is devoted to the proofs of Theorem 3 and Proposition 9.

Technical lemmas

We first state and prove some lemmas that are used in the proof of Theorem 3.

Lemma 1.

Assume that h is a β-Hölder function with Lipschitz constant Λ h defined on R 2 satisfying

0 < a = inf x∈R 2 h(x) sup x∈R 2 h(x) = b < β 1 .
Then, for all x, y ∈ R 2 and for all ξ ∈ R 2 ,

ξ -h(y) -ξ -h(x) Λ h y -x β log ξ ξ -a-1 1 ξ >1 + ξ -b-1 1 ξ 1 .
Proof. Let us fix ξ ∈ R 2 and apply the mean value inequality to the function

h → ξ -h-1 = exp(-(h + 1) log ξ ) .
We obtain that

ξ -h2-1 -ξ -h1-1 |h 1 -h 2 | log ξ ξ -α-1 , with α = min(h 1 , h 2 ) if ξ > 1 and α = max(h 1 , h 2 ) if ξ < 1.
This leads to the inequality:

∀(x, y) ∈ (R 2 ) 2 , ∀ξ ∈ R 2 , ξ -h(y) -ξ -h(x) |h(y) -h(x)| log ξ × ξ -infz h(z)-1 1 ξ >1 + ξ -sup z h(z)-1 1 ξ 1 .
The holderianity of h allows to conclude.

Lemma 2. Assume that h : R 2 → [0, 1] and C : R 2 × R 2 → R + are two functions satisfying assumptions (H). Let x ∈ R 2 . Then, there exists some constant K x > 0 depending only on x such that for any (ρ, w) ∈ R + × R 2 with |ρ| 1 and w 1 we have

R 2 e j x, ξ -1 2 f 1/2 (x + ρw, ξ) -f 1/2 (x, ξ) 2 dξ K x |ρ| 2β max( w 2β , w 2η ) 1 + x 2 .
Moreover the function x → K x is bounded on any compact set.

Proof. Observe that

C(x + ρw, ξ) ξ h(x+ρw)+1 - C(x, ξ) ξ h(x)+1 = C(x + ρw, ξ) ξ h(x+ρw)+1 - C(x + ρw, ξ) ξ h(x)+1 + C(x + ρw, ξ) ξ h(x)+1 - C(x, ξ) ξ h(x)+1 .
Using the classical inequality |a -b| 2 2(|a| 2 + |b| 2 ), we deduce that

R 2 e j x, ξ -1 2 f 1/2 (x + ρw, ξ) -f 1/2 (x, ξ) 2 dξ , 2 R 2 e j x, ξ -1 2 C(x + ρw, ξ) ξ h(x+ρw)+1 - C(x + ρw, ξ) ξ h(x)+1 2 dξ +2 R 2 e j x, ξ -1 2 ξ 2h(x)+2 C(x + ρw, ξ) -C(x, ξ) 2 dξ .
To bound the latter integral R 2 e j x, ξ -1

2 ξ 2h(x)+2 C(x + ρw, ξ) -C(x, ξ) 2 dξ, we set ξ = rΘ with (r, Θ) ∈ R * + × S 1
and use the homogeneity of C. It yieds:

R 2 e j x, ξ -1 2 ξ 2h(x)+2 C(x + ρw, ξ) -C(x, ξ) 2 dξ , = R * + |e js -1| 2 s 2h(x)+1 S 1 | x, Θ | 2h(x) C(x + ρw, Θ) -C(x, Θ) 2 dΘ ds , x 2h(x) R * + |e js -1| 2 s 2h(x)+1 ds S 1 C(x + ρw, Θ) -C(x, Θ) 2 dΘ ,
where we set s = r x, Θ in the second equality. We now use condition [START_REF] Ishikawa | Generalized Hilbert transforms in tempered distributions[END_REF] of assumptions (H) with z = ρw ∈ B(0, 1). Hence,

R 2 e j x, ξ -1 2 ξ 2h(x)+2 C(x + ρw, ξ) -C(x, ξ) 2 dξ A x |ρ| 2η w 2η x 2h(x) R * + |e js -1| 2 s 2h(x)+1 ds . Then, since x 2h(x) x 2 + 1 is always valid (h(R 2 ) ⊂ [0, 1]), one has R 2 e j x, ξ -1 2 ξ 2h(x)+2 C(x + ρw, ξ) -C(x, ξ) 2 dξ B x w 2η |ρ| 2η x 2 + 1 , (25)
with

B x = A x R * + |e js -1| 2 s 2h(x)+1 ds < ∞ .
We now bound R 2 e j x, ξ -1

2 C(x+ρw,ξ) ξ h(x+ρw)+1 -C(x+ρw,ξ) ξ h(x)+1 2
dξ. Since C is bounded and by Lemma 1 we have for some A > 0 depending only on h and C:

R 2 e j x, ξ -1 2 C(x + ρw, ξ) ξ h(x+ρw)+1 - C(x + ρw, ξ) ξ h(x)+1 2 dξ Λ h,C |ρ| 2β w 2β R 2 e j x, ξ -1 2 log ξ 2 ξ -2a-2 1 ξ >1 + ξ -2b-2 1 ξ 1 dξ . Since e j x, ξ -1 min( x • ξ , 2), we directly get that R 2 e j x, ξ -1 2 C(x + ρw, ξ) ξ h(x+ρw)+1 - C(x + ρw, ξ) ξ h(x)+1 2 dξ Λ |ρ| 2β w 2β x 2 + 1 . (26)
with

Λ = Λ h,C R 2 log ξ 2 ξ -2a-2 1 ξ >1 + ξ -2b 1 ξ 1 dξ .
The conclusion then follows from ( 25) and ( 26) with K x = 2B x + 2 A. The fact that x → K x is bounded on any compact set comes from the fact that x → A x is bounded on any compact set.

Proof of Theorem 3

Let X be the Gaussian field defined by formula [START_REF] Koizumi | On the Hilbert transform I*[END_REF], and x 0 ∈ R 2 . Let Z x0 be the Gaussian field

Z x0,ρ (u) = X(x 0 + ρu) -X(x 0 ) ρ h(x0) ,
and Y x0 the H-sssi field defined by formula [START_REF] Dobrushin | Gaussian and their subordinated self-similar random generalized fields[END_REF]. We are going to prove that Y x0 is the tangent field of X at x 0 ∈ R 2 , that is

Z x0,ρ (h) h∈R 2 d -→ Y x0 (h) h∈R 2 .
as ρ → 0, in the sense of weak convergence of stochastic processes. The proof is divided in two steps :

(i) We first prove that the finite dimensional distribution of Z x0,ρ converge to those of Y x0 as ρ → 0:

(Z x0,ρ (h 1 ), . . . , Z x0,ρ (h N )) -→ (Y x0 (h 1 ), . . . , Y x0 (h N )) ,
which means the convergence of the measures of these finite dimensional random vectors on R N . The Lévy theorem insures that it is equivalent to prove the converge in term of the characteristic functions of these random vectors, which is, in the Gaussian case, equivalent to show that we have convergence with respect to the covariance:

∀(u, v) ∈ (R 2 ) 2 , lim ρ→0 E[Z x0,ρ (u)Z x0,ρ (v)] = E(Y x0 (u)Y x0 (v)) . (27)
(ii) Thereafter, we set ρ n = 1/n ∈ [0, 1] and prove that the sequence of random fields (Z n ) n∈N * def = (Z x0,ρn ) n∈N * , satisfies a tightness property, which is fulfilled if (Z n ) n∈N * satisfies the following Kolmogorov criteria (see for example [START_REF] Karatzas | Brownian Motion and stochastic calculus[END_REF] p.64):

∀T > 0, ∀u, v ∈ [-T, T ] 2 , sup n>1 E(|Z n (u) -Z n (v)| γ1 ) C 0 (T ) u -v 2+γ2 , (28)
for some positive constant C 0 (T ) which may depend on T and γ 1 , γ 2 which are universal positive constants.

Remark 4. Since Z x0,ρ (u) -Z x0,ρ (v) is a Gaussian vector, then for all γ 1 > 0 sup ρ∈(0,1) are equal up to a multiplicative constant depending only on γ 1 . The inequality (28) is satisfied by considering γ 1 > 2/γ 0 with γ 2 = γ 0 γ 1 -2. Therefore, it will be sufficient to verify that 0 γ 0 1 such that

u -v -γ0γ1 E Z x0,ρ (u) -Z x0,ρ (v)
u -v -2γ0 E Z x0,ρ (u) -Z x0,ρ (v) 2 γ1/2
∀T > 0, ∀u, v ∈ [-T, T ] 2 , sup n>1 E(|Z n (u) -Z n (v)| 2 ) C 0 (T ) u -v 2γ0 . ( 29
)
Remark 5. Since the notion of tangent field is a local notion at point x 0 , it is equivalent to determine the tangent field of X at x 0 or that of X(x) = ϕ x0 (x)X(x) at point x 0 , where ϕ is a C ∞ function which is equal to 1 on a neighborhood of x 0 (e.g. the ball B(x 0 , a) of radius a) then decreases and vanished outside a compact set (e.g. the ball B(x 0 , b)). Such bump functions can easily be constructed [START_REF] Tu | Bump functions and partitions of unity[END_REF], from a 1-D function ϕ 0 as illustrated in Figure 3 by taking ϕ x0 (x) = ϕ 0 ( xx 0 ). Then, we define

Z x0,ρ (u) = X(x 0 + ρu) -X(x 0 ) ρ h(x0) .
In step (i), computing the limit for given u and v, nothing changes since from a certain rank, ρ n is such that x 0 + ρ n u and x 0 + ρ n v fall in the ball B(x 0 , a) in which X = X.

In step (ii), with

Z n (u) -Z n (v) = ρ -h(x0) n
( X(x 0 + ρ n u) -X(x 0 + ρ n v)), the following inequality

sup n>1 E(| Z n (u) -Z n (v)| 2 ) C 0 (T ) u -v 2γ0 , (30)
requires to distinguish three cases:

• If u and v are in B(x 0 , b) c , then X(x 0 + ρ n u) = X(x 0 + ρ n v) = 0 and so (30) is satisfied.

• If u ∈ B(x 0 , b) and v ∈ B(x 0 , b) c , then Z n (u) -Z n (v) = ρ -h(x0) n X(x 0 + ρ n u) , = ρ -h(x0) n ϕ x0 (x 0 + ρ n u) -ϕ x0 (x 0 + ρ n v) X(x 0 + ρ n u) , with ϕ x0 (x 0 + ρ n u) -ϕ x0 (x 0 + ρ n v) 2 A(T )ρ 2 n u -v 2 since ϕ is C ∞ that is a Lipschitz
function and E X(x 0 + ρ n u) 2 B(T ) by continuity of the covariance function, on the compact set [-T, T ] 2 . Again, (30) is satisfied.

• It remains to deal with the case where u and v are in B(x 0 , b), that is to say we can restrict ourselves for the proof to the case where u and v are in the neighborhood of x 0 as small as we want. We will take for the purposes of the demonstration b = 1/2, in other words we will be able to restrict ourselves to the compact set [-T, T ] = [-1/2, 1/2]. Now we prove these two conditions ( 27) and ( 28).

(i) First step:

We aim at proving [START_REF] Falconer | The local structure of random processes[END_REF] that is for all (u, v) ∈ (R 2 ) 2 :

lim ρ→0 E[(X(x 0 + ρu) -X(x 0 ))(X(x 0 + ρv) -X(x 0 ))] ρ 2h(x0) = E(Y x0 (u)Y x0 (v)) . (31)
Let us take (u, v) ∈ (R 2 ) 2 and consider ρ ρ 0 such that ρ 0 u, ρ 0 v ∈ B(0, 1/2).

First observe that, by definition of X, one has

E[(X(x 0 + ρu) -X(x 0 ))(X(x 0 + ρv) -X(x 0 ))] = R 2 g u,v (x 0 , ξ) dξ , (32)
where we set

g u,v (x 0 , ξ) = e j x0+ρu, ξ -1 f 1/2 (x 0 + ρu, ξ) -e j x0, ξ -1 f 1/2 (x 0 , ξ) × e -j x0+ρv, ξ -1 f 1/2 (x 0 + ρv, ξ) -e -j x0, ξ -1 f 1/2 (x 0 , ξ) ,
and f 1/2 (x, ξ) = C(x, ξ)/ ξ h(x)+1 . We now split the integral into four terms:

E[(X(x 0 + ρu) -X(x 0 ))(X(x 0 + ρv) -X(x 0 ))] = R 2 e jρ u, ξ -1 e -jρ v, ξ -1 f 1/2 (x 0 + ρu, ξ)f 1/2 (x 0 + ρv, ξ) dξ (I 1 ) + R 2 e jρ u, ξ -1 1 -e j x0, ξ f 1/2 (x 0 + ρu, ξ) f 1/2 (x 0 + ρv, ξ) -f 1/2 (x 0 , ξ) dξ (I 2 ) + R 2 1 -e -j x0, ξ e -jρ v, ξ -1 f 1/2 (x 0 + ρu, ξ) -f 1/2 (x 0 , ξ) f 1/2 (x 0 + ρv, ξ) dξ (I 3 ) + R 2 e j x0, ξ -1 2 f 1/2 (x 0 + ρu, ξ) -f 1/2 (x 0 , ξ) f 1/2 (x 0 + ρv, ξ) -f 1/2 (x 0 , ξ) dξ (I 4 ) = I 1 + I 2 + I 3 + I 4 .
In order to prove [START_REF] Polisano | Texture modeling by Gaussian fields with prescribed local orientation[END_REF], we now investigate the behavior of each integral I 1 , I 2 , I 3 , I 4 when ρ → 0.

• Study of the first term I 1

We suppose below that ρ > 0: indeed, since ξ → C(x, ξ) is even, the case ρ < 0 derives in the same way. In the integral I 1 , we set ζ = ρξ (ζ = -ρξ if ρ < 0), dζ = ρ 2 dξ and use the explicit expression of f , then:

I 1 = ρ h(x0+ρu)+h(x0+ρv) R 2 e j u, ζ -1 e -j v, ζ -1 C(x 0 + ρu, ζ/ρ)C(x 0 + ρv, ζ/ρ) ζ h(x0+ρu)+h(x0+ρv)+2 dζ .
By homogeneity of ξ → C(x, ξ), we deduce:

ρ -2h(x0) I 1 = ρ h(x0+ρu)+h(x0+ρv)-2h(x0) × R 2 e j u, ζ -1 e -j v, ζ -1 C(x 0 + ρu, ζ)C(x 0 + ρv, ζ) ζ h(x0+ρu)+h(x0+ρv)+2 dζ .
Observe now that

ρ h(x0+ρu)+h(x0+ρv)-2h(x0) = exp log ρ [h(x 0 + ρu) + h(x 0 + ρv) -2h(x 0 )] .
Using that h is β-Hölder, we obtain

|h(x 0 +ρu)+h(x 0 +ρv)-2h(x 0 )| |h(x 0 +ρu)-h(x 0 )|+|h(x 0 +ρv)-h(x 0 )| ( u β + v β )|ρ| β .
Since β > 0 by assumption, and lim ρ→0 + ρ β log ρ = 0, we then deduce the limit:

lim ρ→0 + log ρ [h(x 0 + ρu) + h(x 0 + ρv) -2h(x 0 )] = 0 , (33)
and hence lim

ρ→0 + ρ h(x0+ρu)+h(x0+ρv)-2h(x0) = 1. It implies that lim ρ→0 + ρ -2h(x0) I 1 = lim ρ→0 + ζ 1 e j u, ζ -1 e -j v, ζ -1 C(x 0 + ρu, ζ)C(x 0 + ρv, ζ) ζ h(x0+ρu)+h(x0+ρv)+2 dζ + lim ρ→0 + ζ 1 e j u, ζ -1 e -j v, ζ -1 C(x 0 + ρu, ζ)C(x 0 + ρv, ζ) ζ h(x0+ρu)+h(x0+ρv)+2 dζ .
We now apply the Lebesgue's Dominated Convergence Theorem to each integral separately. We first bound the two integrands as follows:

∀ ζ 1, (e j u, ζ -1)(e -j v, ζ -1) C(x 0 + ρu, ζ)C(x 0 + ρv, ζ) ζ h(x0+ρu)+h(x0+ρv)+2 M 2 u v ζ 2 ζ 2(b+1) , ∀ ζ 1, (e j u, ζ -1)(e -j v, ζ -1) C(x 0 + ρu, ζ)C(x 0 + ρv, ζ) ζ h(x0+ρu)+h(x0+ρv)+2 4M 2 1 ζ 2(
ρ→0 + ρ -2h(x0) I 1 = R 2 e j u, ζ -1 e -j v, ζ -1 C(x 0 , ζ) 2 ζ 2h(x0)+2 dζ = E(Y x0 (u)Y x0 (v)) .
since the functions h and x → C(x, ζ) are continuous.

• Study of the other terms I 2 , I 3 , I 4

We now prove that the three other integrals I 2 , I 3 , I 4 are negligible with respect to the first one when ρ is small. We only detail the negligibility of I 2 , the other cases I 3 and I 4 being similar. Using the Cauchy-Schwarz inequality, we get that

I 2 R 2 e jρ u, ξ -1 2 f (x 0 + ρu, ξ) dξ 1 /2 × R 2 e j x0, ξ -1 2 f 1/2 (x 0 + ρv, ξ) -f 1/2 (x 0 , ξ) 2 dξ 1 /2 .
The analysis of the first integral has already been done in the study of I 1 taking u = v. We then obtain lim

ρ→0 ρ -2h(x0) R 2 e jρ u, ξ -1 2 f (x 0 + ρu, ξ) dξ = E[Y x0 (u) 2 ] . (34)
The bound of the second integral directly comes from Lemma 2. Since β sup x h(x), we get that lim ρ→0 ρ -2h(x0) I 2 = 0 .

The same approach also yields for I 3 and I 4 , leading to

lim ρ→0 ρ -2h(x0) I 3 = lim ρ→0 ρ -2h(x0) I 4 = 0 .
which concludes the proof of ( 31).

(ii) Second step:

We now have to prove that the sequence (Z n ) n∈N * def = (Z x0,ρn ) n∈N * satisfies (28), with ρ ≡ ρ n = 1 n ∈ [0, 1] and with u and v restricted to B(0, 1/2) from Remark 5, that in the following T = 1/2. Then, we have:

Z x0,ρ (u) -Z x0,ρ (v) = 1 ρ h(x0) R 2 e j x0+ρu, ξ -1 f 1/2 (x 0 + ρu) -e j x0+ρv, ξ -1 f 1/2 (x 0 + ρv) W(dξ) .
Hence,

E Z x0,ρ (u) -Z x0,ρ (v) 2 = 1 ρ 2h(x0) R 2 e j x0+ρu, ξ -1 f 1/2 (x 0 + ρu, ξ) -e j x0+ρv, ξ -1 f 1/2 (x 0 + ρv, ξ) 2 dξ , = 1 ρ 2h(x0) R 2 e j x0+ρu, ξ -1 f 1/2 (x 0 + ρu, ξ) -f 1/2 (x 0 + ρv, ξ) -e j x0+ρv, ξ -1 -(e j x0+ρu, ξ -1) f 1/2 (x 0 + ρv, ξ) 2 dξ , 2 ρ 2h(x0) R 2 e j x0+ρu, ξ -1 2 f 1/2 (x 0 + ρu, ξ) -f 1/2 (x 0 + ρv, ξ) 2 dξ + 2 ρ 2h(x0) R 2 e j ρ(v-u), ξ -1 2 f (x 0 + ρv) dξ .
the last inequality coming from |a -b| 2 2|a|

2 + 2|b| 2 . We now apply Lemma 2 with x = x 0 + ρu and w = vu ∈ B(0, 1). It implies the following bound for the first integral

2 ρ 2h(x0) R 2 e j x0+ρu, ξ -1 2 f 1/2 (x 0 + ρu, ξ) -f 1/2 (x 0 + ρv, ξ) 2 dξ , 2 sup y∈B(x0,2T ) K y 1 + 2 x 0 2 + 2 u 2 max v -u 2β , v -u 2η |ρ| 2β-2h(x0) , K 1 |ρ| 2β-2h(x0) max v -u 2β , v -u 2η , (35)
with

K 1 = 2 sup y∈B(x0,2T ) K y 1 + 2 x 0 2 + 4T 2 and β -h(x 0 ) positive.
To bound the second one observe that C is homogeneous w.r.t. the second variable and bounded. Set ζ = ρξ uv and deduce that

2 ρ 2h(x0) R 2 e j ρ(v-u), ξ -1 2 f (x 0 + ρv) dξ 2 u -v 2h(x0+ρu) C ∞ ρ 2h(x0+ρu)-2h(x0) × R 2 e j Θ, ζ -1 2 ζ -2a-2 1 ζ >1 + ζ -2b-2 1 ζ 1 dζ ,
with Θ = (uv)/ uv . By the same arguments than in [START_REF] Polisano | Simulation of oriented pattern with prescribed local orientation[END_REF] we have that lim ρ→0 ρ 2(h(x0+ρv)-h(x0)) = 1 , so it has a finite upper bound A 1 > 0, which is achieved on a compact by continuity

A 1 = max ρ,v (ρ, v) ∈ [0, 1] × [-T, T ] 2 , ρ 2(h(x0+ρv)-h(x0)) .
Identically we have uv 2h(x0+ρv) = uv 2(h(x0+ρv)-h(x0)) uv 2h(x0) and the first term tends to 1, then the function (ρ, u, v) → uv 2(h(x0+ρv)-h(x0)) also achieves its upper bound A 2 > 0. Thus,

2 u -v 2h(x0+ρu) C ∞ ρ 2h(x0+ρu)-2h(x0) 2A 1 A 2 C ∞ u -v 2h(x0) .
Hence, using that e j Θ, ζ -1 R 2 e j ρ(v-u), ξ -1

2 f (x 0 + ρv, ξ) dξ K 2 u -v 2h(x0) , (36)
with

K 2 = 2A 1 A 2 C ∞ R 2 min(2, ζ 2 )( ζ -2a-2 1 ζ >1 + ζ -2b 1 ζ 1 ) dζ.
Since K 1 , K 2 are two positive constants depending only on T , inequalities [START_REF] Perrin | Reducing non-stationary random fields to stationarity and isotropy using a space deformation[END_REF] and [START_REF] Guyon | Identification of space deformation using linear and superficial quadratic variations[END_REF] imply that

E Z x0,ρ (u) -Z x0,ρ (v) 2 K 2 u-v 2h(x0) 1 + K 1 K 2 max v -u 2(β-h(x0)) , v -u 2(η-h(x0)) ,
and since β -h(x 0 ) > and η -h(x 0 ) > 0, the second factor achieved its bounds on the compact set [-T, T ] 2 × [-T, T ] 2 , hence sup ρ∈(0,1)

u -v -2h(x0) E Z x0,ρ (u) -Z x0,ρ (v) 2 < ∞ .
Which proves the inequality (28) from Remark 4 with γ 0 = h(x 0 ).

The proof of the points (i) and (ii) completes the proof of the Proposition 3.

Proof of Proposition 9

Let x 0 ∈ R 2 . Since X is H-self-similar, one has

X(x) = R 2
(e j x,ξ -1)f 1/2 (ξ) W(dξ) , with f (ξ) = C X (ξ) ξ -2H-2 . Then,

Z Φ,X (x) = X(Φ(x)) = R 2
(e j Φ(x),ξ -1)f 1/2 (ξ) W(dξ) .

As for the proof of Proposition 3 in Section 5.2, we divide the following proof into two steps.

(i) First step:

Let u, v ∈ R 2 , and consider:

1 ρ 2H E Z Φ,X (x 0 + ρu) -Z Φ,X (x 0 ) Z Φ,X (x 0 + ρv) -Z Φ,X (x 0 ) = 1 ρ 2H
R 2 e j Φ(x0), ξ e j Φ(x0+ρu)-Φ(x0), ξ -1 e -j Φ(x0), ξ e -j Φ(x0+ρv)-Φ(x0), ξ -1 f (ξ) dξ , The second inequality is obtained by Cauchy-Schwarz inequality, the third by mean value inequality, the forth under the assumption that Φ which is continuously differentiable so Φ ′ K, the fifth with C = K max( u , v ). Finally, we show that G is integrable since: (ii) Second step:

R 2 G(ζ) dζ = R 2 min(2, C ζ ) 2 f (ζ) dζ , = 1 C 2 R 2 min(2, ξ ) 2 f ξ C dξ , = C 2H+2 C 2
We then prove that the convergence holds in the sense of finite dimensional distributions. To deduce Proposition 9, we follow the same way as in Step 2 of Theorem 3 a Kolmogorov criteria.

E

Z Φ,X (x 0 + ρu) -Z Φ,X (x 0 ) ρ H -Z Φ,X (x 0 + ρu) -Z Φ,X (x 0 ) ρ H Φ ′ (x) ρ(uv)

2H

,

I H I C Φ ′ 2H u -v 2H , C 0 u -v 2H .
with C 0 = I H I C Φ ′ 2H . We conclude like at the end of Step 2 of Theorem 3.

Figure 1 :

 1 Figure 1: Image texture of size 512 × 512 resulting from the simulation of the field Z Φ,X (x) = X(R -α(x) x) on [0, 1] 2 , where X is the standard elementary field with parameters H = 0.5, α 0 = 0 and δ = 0.3, for the following functions α: (a) α(x 1 , x 2 ) = -π 3 (top left), (b) α(x 1 , x 2 ) = -π 2 + x 1 (top right), (c) α(x 1 , x 2 ) = -π 2 + x 2 (bottom left), (d) α(x 1 , x 2 ) = -π 2 + x 2 1 -x 2 (bottom right).

Figure 2 :

 2 Figure 2: Image texture of size 512 × 512 resulting from the simulation of the field Z Φ,X (x) = X(Φ(x)) on [0, 1] 2 , where X is the standard elementary field with parameters H = 0.5, α 0 = 0 and δ = 0.3, and Φ is the deformation function defined by (24) corresponding to the harmonic function α(x 1 , x 2 ) = ax 1 + bx 2 + c, with (a, b) = (2, -1).

Figure 3 :

 3 Figure 3: Ilustration of a bump function ϕ 0 .

  a+1) , In the first line we used |e jt -1| 2 | sin(t/2)| |t|, and the Cauchy-Schwarz inequality applied to the R 2 -scalar product. Secondly, since a > 0 and b < 1 by assumption, we easily check that the functions ζ → ζ -2b and ζ → ζ -2(a+1) are respectively integrable on ζ 1 and ζ 1. The Lebesgue's Dominated Convergence Theorem then implies that lim

2 min( 2 ,

 22 Θ ζ ) = min(2, ζ ), 2 ρ 2h(x0)

Φ

  +ρv)-Φ(x 0 ) ρ , ζ -1 f (ζ/ρ) dζ/ρ 2 , +ρv)-Φ(x 0 ) ρ , ζ -1 ρ 2H+2 f (ζ) dζ/ρ 2 , 0 +ρu)-Φ(x 0 ) ρ , ζ -1 e -j Φ(x 0 +ρv)-Φ(x 0 ) ρ , ζ -1 f (ζ) dζ .To compute the limit of this quantity when ρ → 0, let us denote by g(ρ, ζ) the integrand of the last integral. We havelim ρ→0 g(ρ, ζ) = e j DΦ(x0)u, ζ -1 e -j DΦ(x0)v, ζ -1 f (ζ) .Now we have to bound the integrand |g(ρ, ζ)|: using the inequality |e jx -1| min(2, |x|), one has|g(ρ, ζ)| min   2, Φ(x 0 + ρu) -Φ(x 0 ) 0 + ρu) -Φ(x 0 ) ζ min 2, 1 ρ Φ(x 0 + ρv) -Φ(x 0 ) ζ f (ζ) , ′ (x) ρv ζ f (ζ) , min 2, K u ζ min 2, K v ζ f (ζ) , min(2, C ζ ) 2 f (ζ) ≡ G(ζ) .

R 2 min( 2 ,2H R 2 min( 4 , ξ 2 2 e

 222422 ξ ) 2 f (ξ) dξ , C )f (ξ) dξ < +∞ .where we have used the homogeneity of f , and Proposition 4. Hence, using the Lebesgue's Dominated Convergence Theorem, we obtainlim ρ→0 1 ρ 2H E[(Z Φ,X (x 0 + ρu) -Z Φ,X (x 0 ))(Z Φ,X (x 0 + ρv) -Z Φ,X (x 0 ))] = E[Y x0 (u)Y x0 (v)]where we denotedY x0 (u) = R j DΦ(x0)u, ξ -1 f 1/2 (ξ) W(dξ) ,which is by definition the tangent field.

2 e 2 e 2 f 2 f 1 Φ 1 C 2 I

 2222112 j Φ(x0+ρu), ξ -e j Φ(x0+ρv), ξ 2 f (ξ) dξ , j Φ(x0+ρu)-Φ(x0+ρv), ξ -1 (ξ) dξ , = R 2 e j Φ(x 0 +ρu)-Φ(x 0 +ρv) ρ , ζ -1 (ζ) dζ , (x 0 + ρu) -Φ(x 0 + ρv) ρ , Θ 2H C X (Θ) dΘ ds , Φ(x 0 + ρu) -Φ(x 0 + ρv) ρ X (Θ) dΘ . with s = r Φ(x0+ρu)-Φ(x0+ρv) ρ, Θ . Let denote by I H the first integral above and I C the second one. Then, since Φ is C 1 , one haveE Z Φ,X (x 0 + ρu) -Z Φ,X (x 0 ) ρ H -Z Φ,X (x 0 + ρu) -Z Φ,X (x 0 ) ρ H H I C ρ 2H sup [x0+ρu,x0+ρv]

where the

2D Fourier transform is defined for f ∈ L 1 (R 2 ) by f (ξ) = R 2 f (x)e -j x, ξ dx and then extended by the Plancherel theorem and by continuity arguments for f ∈ L 2 (R 2 ).

Acknowledgements

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-13-BS03-0002-01 (ASTRES).


B Jodar

G Seisson

D Hébert

I Bertron

M Boustie

L Berthe

Enhancement of a dynamic porous model considering compression-release hysteresis behavior: application to graphite

Keywords: POREQST, Pα, carbon, foam, constitutive relation

Because of their shock wave attenuation properties, porous materials and foams are increasingly used for various applications such as graphite in the aerospace industry and polyurethane (PU) foams in biomedical engineering. For these two materials, the absence of residual compaction after compression and release cycles limits the efficiency of the usual numerical dynamic porous models such as Pα and POREQST. In this paper, we suggest a simple enhancement of the latter in order to take into account the compressionrelease hysteresis behavior experimentally observed for the considered materials. The new model, named HPOREQST, was implemented into a Lagrangian hydrocode and tested for simulating plate impact experiments at moderate pressure onto a commercial grade of porous graphite (EDM3). It proved to be in far better agreement with experimental data than the original model which encourages us to pursue numerical tests and developments.

Introduction

The shock wave attenuation capability of foams and porous materials is one of the main research topics in preventing damage caused by the propagation of shock waves. The mech anisms responsible for the shock wave attenuation in such materials are not yet completely understood nor described. For instance, a porous graphite similar to those used in the aerospace industry as a constitutive element of C/C composite materials, revealed no residual compaction after static confined compression (see figure 1 and [START_REF] Seisson | Experimental and theoretical study of the damaging of graphite under dynamic loading-application to hypervelocity impacts[END_REF]). In biomedical engineering, polyurethane (PU) foams, used as a specimen to understand shock interaction with complex media, showed a similar behavior during uniaxialstrain experiments (see figure 1(b) and [START_REF] Kitagawa | Attenuation of shock waves propagating in polyurethane foams[END_REF]). This particular behavior is also expected during shock wave propagation. The concern of this paper is to pave the way for developing a numerical model able to correctly reproduce such complex mechanisms under shockandrelease loading.

To describe the dynamic behavior of porous materials, sev eral numerical models have been developed during the last decades. Among them, the simplest one is the 1D model of Thouvenin [START_REF] Thouvenin | Action d'une onde de choc sur un solide poreux[END_REF]. It can be called a geometric model since it consists in describing a porous material as an assembly of parallel slabs of dense matter. The ratio between the slab thickness and the empty spaces between them is related to the macroscopic density of the material. It originally gave good results for sustained shocks onto porous metals but was also tested for laserinduced shocks onto aluminum [START_REF] Zagouri | Laser shock compaction of porous materials[END_REF] and porous graphite [START_REF] Hébert | Plate impact experiments and simulation on porous graphite[END_REF]. Nevertheless, its validity is limited to a range of pressures where compaction is complete and where materials do not undergo phase transition since it assumes the symmetry of the Hugoniot and the release isentrope.

Conversely, nongeometric models such as those expressed by Zel'dovich and Raizer [START_REF] Zel | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF], McQueen et al [START_REF] Mcqueen | The equation of state of solids from shock wave studies High Velocity Impact[END_REF] and Charakhch'yan [START_REF] Khishchenko | Hydrodynamic simulation of converging shock waves in porous conical samples enclosed within solid targets[END_REF][START_REF] Charakhch'yan A A | Shock compression of some porous media in conical targets: numerical study[END_REF][START_REF] Charakhch'yan A A | Converging shock waves in porous media Tech[END_REF] allow 2D and 3D computations. They suppose an almost instantaneous density transition between the initial porous state and the solid state. Hence, they have the property to correctly deal with highpressure phase trans itions but cannot properly manage lowpressure ranges where compaction is uncomplete.

The Pα model, firstly introduced by Herrmann [START_REF] Herrmann | Constitutive equation for the dynamic compaction of ductile porous materials[END_REF], is certainly the most known and used porous model in the litera ture. It is also able to deal with phase transition and palliates the lacuna of the models of Thouvenin, Zel'dovich, McQueen and Charakhch'yan when compaction is uncomplete. It uses a plastic compaction curve that manages the gap between the initial porous state to the dense one. The totally compacted material is modeled by the original dense equation of state (EOS). After a partial compaction, the return to a zero pres sure state is performed following elastic surfaces. Several developments of this model were then proposed, e.g. [START_REF] Carroll | Static and dynamic pore collapse relations for ductile porous materials[END_REF][START_REF] Wünnemann | A strain based porosity model for use in hydrocode simulations of impacts and implications for transient crater growth in porous targets[END_REF], but none of them accounted for a potential hysteresis behavior.

The POREQST model expressed by Seaman et al [START_REF] Seaman | Computational representation of constitutive relations for porous materials[END_REF] is based on fundamentals similar to Pα but the compac tion curve is easier to fit on experimental data. Moreover, the elastic surfaces are defined by a Mie-Grüneisen EOS that allows computing time to be saved since the porosity is ana lytically calculated instead of with an algorithm like for Pα. It is also able to deal with negative pressures introducing a pore reopening curve and all of its parameters have a physical meaning, most of them being obtainable at the first order by quasistatic tests. Extensive investigations have already been conducted by the present authors [START_REF] Seisson | Experimental and theoretical study of the damaging of graphite under dynamic loading-application to hypervelocity impacts[END_REF][START_REF] Hébert | Plate impact experiments and simulation on porous graphite[END_REF][START_REF] Seisson | Dynamic behavior of a porous brittle material: experiments and modeling[END_REF][START_REF] Seisson | Modeling of laserdriven waterconfined shocks into porous graphite[END_REF][START_REF] Seisson | Dynamic cratering of graphite: experimental results and simulations[END_REF] on this model in order to shape a set of parameters adapted to the commer cial grade of porous graphite evoked in figure 1. Therefore POREQST has been chosen for being enhanced in order to reproduce a compressionrelease hysteresis behavior.

The next section recalls the main features of the original POREQST model before proposing a substantial modification for the representation of the compressionrelease hysteresis behavior. The third section will adjust the new model, called HPOREQST, for a commercial grade of graphite (EDM3) and discuss its limitations and benefits.

Enhancement of POREQST

Original model

The original POREQST model [START_REF] Seaman | Computational representation of constitutive relations for porous materials[END_REF] is implemented into Hésione, a Lagrangian-Eulerian hydrocode developed at CEA. It supplies constitutive relations for porous materials submitted to dynamic loading. Figure 2 presents an overview of its features in the zero internal energy plane.

2.1.1. Dense material. Let ρ s be the initial density of the dense (or solid) material at zero pressure P and internal energy e. The EOS of the solid can be any tabulated EOS such as a SESAME or a Mie-Grüneisen expressed as follows: with ρ the current density, K s the dense compressibility modulus and Γ the Grüneisen coefficient of the solid which is assumed to be constant. The two parameters of the elasticplastic behavior law are tabulated and depend on P and e:

( ) ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ⎛ ⎝ ⎜ ⎜ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ⎞ ⎠ ⎟ ⎟ ρ ρ ρ ρ ρ ρ = - - Γ - +Γ P e K e , 1 1 2 1 s s s s (1)
• = ( ) ( ) ( ) Y P e Y P f e , s s 1
: the elastic limit of the solid;

• G P e G P f e , s s 2 = ( ) ( ) ( ): its shear modulus.
f 1 and f 2 are tabulated softening functions varying between 0 for the melting internal energy and 1 for a zero energy.

( ) ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ρ ρ ρ ρ = - P e K T , , 1
i i i i f ( 2
)
with

T e K K K f e L L K K 1 , 1 , 1 f s s i s i i s 3 1 1 0 0 0 ρ α α α α = + Γ = + - = - - ( ) ( ) /
(3) where T f is a term of thermal dilatation which is constant what ever the intermediate surface. f 3 (as the following functions f i of the model) is a tabulated softening function similar to f 1 and f 2 . Moreover, we define / α ρ ρ = i s i as the material disten sion of the porous material at the intermediate density ρ i . In the specific case of the first intermediate surface, i.e. the initial elastic surface of the porous material, we have

= K K i 0 , ρ ρ = i 0
and then / α α ρ ρ = = i s 0 0 . Note that the hidden hypothesis is that the dilatation coefficient at zero pressure is the same whatever the intermediate surface. It means that the volume V is written as:

( ) α α ρ = + ∆ = Γ V V T C K 1 w ith . i t t s p s ( 4
)
As regards the deviatoric stress, the elastic limit and the shear modulus progress according to the intermediate density and the internal energy in the same way as the compressibility modulus:

Y Y f e L L Y Y 1 with 1 i s i i s 4 2 2 0 0 0 α α α α = + - = - - ( ) ( ) / ( 5
)
α α α α = + - = - - ( ) ( ) / G G f e L L G G 1 with 1 .
In the principal system, let the extrem deviatoric stresses be:

( ) ( ) = = = = S S S S max a nd min . i i i i max 1 ,3 min 1 ,3 (8)
In compression, the limitation of the intermediate surfaces is fulfilled by a compaction surface described by:

( ) ( ) ( ) σ ρ ρ = Σ e f e , c c 6 (9)
where ( ) ρ Σ c results from experimental data of a static uniaxial strain test in compression of the healthy porous material such as displayed in figure 2 (dashed line).

In traction, the limitation is undertaken by a pore re opening curve whose intersect at zero energy with the dense material surface occurs at (ρ l ; σ l ):

( ) ( ) ( ) ( ) σ ρ ρ ρ ρ σ ρ = Σ Σ = e f e T , w ith o o f o l l 7 ( 10
)
σ l being the main parameter of this surface. Note that in prac tice, ( )

( ) ( ) = = f e f e f e 3 6 7
.

2.1.3. General algorithm. Then, it is necessary to perform a series of tests (see algorithm 1) in order to establish what value of the pressure the code must use. For that purpose, we note ID the index caracterizing the behavior pattern:

• ID = 0: no change of intermediate surface;

• ID = 1: pore reopening;

• ID = 2: compaction;

• ID = 3: complete densification. [START_REF] Herrmann | Constitutive equation for the dynamic compaction of ductile porous materials[END_REF] It means that damage is activated when both of the following conditions are met: (1) the matter in a Lagrangian cell has a negative pressure located on the pore reopening surface of the POREQST model (i.e. ID = 1); (2) the intermediary distension in this cell is higher than the critical value α lim . In that case, pressure and deviatoric stress are relaxed to zero, and a damage variable D is set to 1. Then in each cell, σ l is maintained to its original value if D = 0 or set to 0 if D = 1, which is equivalent to withdrawing the local tension strength of the damaged material.

α ρ ρ α = = ⩽ ID If 1, . s i i lim

New model: H-POREQST

The original POREQST (just as the Pα model and its deriva tives) obviously predicts a residual compaction after compres sion. As written in the introduction, this is not realistic regarding the available data for such materials as graphite and polyurethane foams which present compressionrelease hysteresis cycles. This work suggests an enhancement of the original POREQST model adapted to materials that experience com pressionrelease hysteresis cycles and calls it HPOREQST (i.e. hysteresisPOREQST). For consistency with the phe nomenological framework of the original model, we intro duce a simple plane surface denominated hysteresis surface as shown in figure 3. It is only defined by the parameter σ h which refers to its intersection with the dense EOS in the zero internal energy plane and its equation is given by:

σ ρ ρ ρ σ ρ ρ ρ ρ = Σ Σ = - - ( ) ( ) ( ) ( ) ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ e f e T , w ith h h h f h hys 8 0 0 ( 12
)
where ρ h is the density of the intersection point and T f the same thermal dilatation coefficient defined by equation [START_REF] Thouvenin | Action d'une onde de choc sur un solide poreux[END_REF]. A numerical adjustment of the parameter σ h will be required to reproduce the experiments. Henceforth, pressure relaxation follows successively the intermediary and the hysteresis sur faces that involves two distinct ranges of release wave velocity. Tension is only permitted when the current density is smaller than the initial one, i.e. as long as the distension is greater than / α ρ ρ = s 0 0 . In that case, the pressure evolution will follow the negative part of the hysteresis surface. This behavior is not wellsuited for a tensile initial loading but might be appropriate for tension after a compressionrelease cycle during which the microstructure of the material may have changed. Hence, HPOREQST must be considered rele vant only for computations without gives the series of tests necessary to implement the new model (the dissimilarities with the original POREQST are in blue) thanks to a new index of behavior pattern defined as:

• ID = 4: hysteresis.
Finally, just as POREQST, our new model will encounter density anomalies in tension. To avoid this phenomenon, the damage criterion and treatment enunciated in section 2.1.4 are also used with HPOREQST. Hence, α lim is set to 1, meaning that pore reopening is never authorized.

Application to graphite and discussion

The EDM3, a polycrystalline graphite produced by sintering, is manufactured by the compagny POCO [18,[START_REF] Sheppard R G, Morgan | Properties and Characteristics of Graphite for the EDM Industry 5th edn[END_REF] [START_REF] Seisson | Experimental and theoretical study of the damaging of graphite under dynamic loading-application to hypervelocity impacts[END_REF][START_REF] Hébert | Plate impact experiments and simulation on porous graphite[END_REF][START_REF] Seisson | Modeling of laserdriven waterconfined shocks into porous graphite[END_REF] in order to investigate its behavior under quasistatic and dynamic loading. Quasistatic confined compressionrelease cycles revealed a hysteresis behavior, as displayed in figure 1(a). In the case of dynamic loadings, the release waves are expected to be governed by this phenomenon inducing two distinct ranges of velocities. Furthermore, the experiment tends to show that the density of EDM3 returns to its initial value after compression, i.e. EDM3 does not experience residual compaction.

EDM3 modeling

The main parameters used for modeling EDM3 with POREQST and HPOREQST are given in table 1. The dense EOS is the tabulated SESAME 7832 of carbon, the compaction curve is adapted from the dashed line of figure 1(a), and the intermediate surfaces are defined through the elastic moduli K and G. The deviatoric stress is constrained by the VonMises criterion and the elastic limit Y. Note that G s and Y s are tabu lated piecewiselinear parameter σ l is set to -95 MPa in order to make the intersect between the pore reopening surface and the initial intermediate surface at zero energy to be equal to -70 MPa, which is the tension strength of healthy EDM3 [START_REF] Seisson | Experimental and theoretical study of the damaging of graphite under dynamic loading-application to hypervelocity impacts[END_REF]. The limit distension α lim is set to 1, i.e. the condition of equation ( 11) is always fulfilled. Thus, the use of the pore reopening surface is totally prohibited, thereby transcribing the brittleness of the material (fracture toughness of about 1 MPa • m 1/2 [START_REF] Takahashi | Impact fracture toughness of a nuclear graphite measured by the onepoint bending method[END_REF][START_REF] Ayatollahi | Mixed mode brittle fracture of sharp and blunt vnotches in polycrystalline graphite[END_REF][START_REF] Burchell | A microstructurally based fracture model for polygranular graphites[END_REF][START_REF] Salazar | In situ observation of damage nucleation in graphite and carbon/carbon composites[END_REF][START_REF] Latella | The initiation and propagation of thermal shock cracks in graphite[END_REF][START_REF] Mirhabibi | Graphite flakecarbon composites. II: fracture behaviour, toughness, notch insensitivity and Weibull modulus[END_REF]). Last, the unique parameter of the hysteresis surface σ h is equal to 200 MPa (see section 3.2).

As shown in figure 4, this set of parameters allows a satis fying match between the test from figure 1

/ K f e T P L K f e PL 1 i s f s s 1 1 1 1 ρ ρ ρ - + - ← ( ( ) (
))/( ( ) ) ; end end up to 1 GPa with HPOREQST and the hydrocode Hésione. It was a prerequisite prior to simulating dynamic experiments.

Plate impact test

In order to appraise HPOREQST, two experiments with EDM3 samples have been performed. Both consisted in a 2 mmthickness copper plate launched at 330 m s -1 by a mono stage lightgas gun onto a 3 mmthickness plate of graphite. A PDV probe [START_REF] Strand | Using the heterodyne method to measure velocities on shock physics experiments[END_REF][START_REF] Mercier | PDV measurements of ns and fs laser driven shock experiments on solid targets[END_REF][START_REF] Hébert | A KDP equation of state for laserinduced damage applications[END_REF] recorded the rear free surface velocity of the latter. Since the experiments were perfectly repeatable only one velocity is displayed in figure 5(a) (circles).

Numerical simulations were performed with Hésione in a planar 1D Lagrangian configuration. The maximum generated stress onto the graphite target rear face was around 400 MPa. It allowed the model to be tested in the middle range of the static uniaxial compression data of figure 1.

Results are given in figure 5(a). The continuous line gives the rear surface velocity of EDM3 obtained with the original model of POREQST associated to the damage model pre sented in section 2.1.4, whilst the dotted lines display sig nals computed with the model HPOREQST with the same damage model. In both case, the curve foot is squared due to the successive rebounds of the shock elastic precursor on the rear surface of the target corresponding to the elastic initial surface of the models. Throughout an adjustment of the value of σ h was set to 200 MPa, the new model clearly improves the first peak of the freesurface velocity. Furthermore, pullback the spall ejec tion velocity satisfactorily oscillates around the experimental signal.

Discussion

The increase of the maximum free surface velocity up to the experimental level induced by HPOREQST confirms that hysteresis behavior also occurs in the dynamic regime. It is well understood thanks to figure 5(b) which gives the numer ical σu diagram of the target recorded close to the back face, u being the particle velocity and σ the longitudinal stress. The target undergoes a nonmonotonic compression (because of the rebounds of the elastic precursor on the back face) (1) up to a maximum stress of about 400 MPa [START_REF] Kitagawa | Attenuation of shock waves propagating in polyurethane foams[END_REF]. The shock break out on the rear surface involves a stress decrease (3) whose shape depends on the model. Original POREQST implies a straightlike pressure decay (plain) due to the intermediate sur faces generating a low free surface velocity around 360 m s -1 . Whereas HPOREQST implies an upper free surface velocity around 450 m s -1 coming from the twostraightpart pressure However, despite the simplicity of the modeling, the overall agreement between the experiment and the computa tion with HPOREQST is very satisfactory which paves the way for potential further improvements, when more various experimental results will be available.

Conclusion

The purpose of this paper was to enhance the POREQST model in order to be able to numerically reproduce the assumed density recovery after shock compression of certain porous materials and foams.

Hence, we introduced a new constitutive surface beside the original ones and called it the hysteresis surface that is defined through a single parameter. This simple modeling was tested for plate impact experiments on EDM3, a commercial porous graphite. It clearly improved the target free surface velocity first peak as well as the spall phenomenon and its flight average velocity, thereby confirming the existence of the hysteresis behavior under dynamic shockandrelease loading.

Finally, testing this model on more different loading cases (unsustained shock, compression up to higher pressures and release, etc) might be relevant for its improvement and vali dation. Then, the new HPOREQST model should be tested for computing hypervelocity impacts of steel spheres into graphite [START_REF] Hébert | Simulations of hypervelocity impacts into graphite[END_REF][START_REF] Seisson | Penetration and cratering experiments of graphite by 0.5 mm diameter steel spheres at various impact velocities[END_REF] where it may help to simulate the crater refill phenomenon observed experimentally. σu diagram of the experiment extracted from numerical simulations (with artifical viscosity [START_REF] Wilkins | Use of artificial viscosity in multidimensional fluid dynamic calculations[END_REF]). Comparison between POREQST and HPOREQST.

Figure 1 .

 1 Figure 1. Highlighting of the compressionrelease hysteresis phenomenon for polycrystalline porous graphite and polyurethane foam. For both materials, no residual compaction exists after compression (at least up to 1.1 GPa for EDM3 and 0.26 MPa for PU foam). (a) Uniaxial strain tests in compression (dashed) and cycled compression and release (solid) onto an isotropic porous grade of graphite. Adapted from [1]. (b) Uniaxial strain test in compression and release onto a PU foam. Extracted with permission from [2]. Copyright 2006 SpringerVerlag.

Figure 2 .

 2 Figure 2. Schematic of the original POREQST model adapted from [14] in the zero internal energy plane. The dashed arrow depicts a typical compressiontension path (without damage model).

  α α = i 0 we obtain the initial values Y 0 and G 0 corresponding to the initial porous material; when α = 1 i , we have the values of the solid material Y s and G s . Delimitation of the intermediate surfaces. Assuming a posi tive stress in compression, it follows that:

  = 1 or ID = 2 or ID = 4 then /* Deduce the new intermediate surface *

Figure 4 .

 4 Figure 4. Comparison between the quasistatic cycled test from figure 1 and its simulation with HPOREQST and Hésione using the parameters of EDM3 from table 1. As an illustration of the discussion of section 3.3, both simulations with σ = 100 h MPa (black rounds) and with σ = 200 h MPa (red squares) are given. Discontinuities in the slope of the intermediate surfaces are logically due to the elastoplastic modeling.Table 1. Main parameters used for computations with POREQST and HPOREQST performed in section 3.2.

Figure 3 .

 3 Figure 3. Schematic of the new HPOREQST model in the zero internal energy plane. The dashed arrow depicts a typical compressiontension path (without damage model). Dotted lines are unusable surface parts.

Figure 5 .

 5 Figure 5. Comparison between the original model and the enhanced one in the case of a plate impact experiment on EDM3 with spallation. (a) Experimental free surface velocity of the EDM3 target compared to the simulations results with POREQST and HPOREQST. (b)σu diagram of the experiment extracted from numerical simulations (with artifical viscosity[START_REF] Wilkins | Use of artificial viscosity in multidimensional fluid dynamic calculations[END_REF]). Comparison between POREQST and HPOREQST.

  The porous material is represented as a succession of intermediate elastic surfaces that are delimited by two other surfaces, one accounting for its compaction and the other modeling the pore reopening.

2.1.2. Porous material.

Intermediate surfaces. The first of the intermediate sur faces i is the elastic surface of the porous material. Their density at zero pressure and energy, called intermedi ate density, is denominated ρ i . Along these intermediate surfaces, the pres sure is also expressed thanks to a Mie-Grüneisen EOS:

  . Its porosity Algorithm 1. General algorithm of POREQST, where e melt is the assumedconstant melting energy by volume unit.

	begin
			← ID 0;
			If ⩾ e e melt then
				/* Use the dense EOS ← ( ) ρ P P e , s	*/
				← ID 3
			else
				/* Research of the correct pressure expression If ( ) ( ) ρ σρ -< P e S e , , i o max and ( ) ( ) ρ σρ -< P e S e , , s o max then	*/
				P	← ( ) σ ρ + e , o	S max ;
				← ID 1;
				end
				If ( ) ρ P e , s	>	( ) ρ P e , i	and ID = 0 then
				← ( ) ρ P P e , s ;
				← ID 3;
				end
				If ( ) ρ P e , s	-	S min	>	( ) σρ e , c	and ID = 0 then
				← ( ) ρ P P e , s ;
				← ID 3;
				end
				If ( ) ρ P e , i	-	S min	>	( ) σρ e , c	and ID = 0 then
				P	← ( ) σ ρ + e , c	S min ;
				← ID 2;
				end
				If ( ) ρ P e , i	-	S min	<	( ) σρ e , c	and ID = 0 then
				← ( ) ρ P P e , i ;
				← ID 0;
				end
			end
			If ID = 1 or ID = 2 then
				/* Deduce the new intermediate surface K f e T P L K f e PL 1 i s f s s 1 1 1 1 ρ ρ ρ -+ -← ( ( ) ( ))/( ( ) ) ;	*/
			end
	end
	(i.e.	1	-	/ ρ ρ s 0 ) is about 20% and it is macroscopically homo
	geneous and isotropic. Characterization tests were conducted
	by the present authors

  (a) and its simulation Algorithm 2. General algorithm of HPOREQST, where e melt is the assumedconstant melting energy by volume unit. In blue are the dissimilarities with the original POREQST. (ρ, e) -S max < σ hys (ρ, e) and σ o (ρ, e) < σ hys (ρ, e) then P ← σ hys (ρ, e) + S max ;

	begin
	← ID 0;
	If ⩾ e e melt then
	/* Use the dense EOS ← ( ) ρ P P e , s	*/
	← ID 3
	else
	/* Research of the correct pressure expression if ( ) ( ) ρ σρ -< P e S e , , i o max and ( ) ( ) ρ σρ -< P e S e , , s o max then	*/
	P	← ( ) σ ρ + e , o	S max ;
	← ID 1;
	end
	if P i ID ← 4;
	end
	if ( ) ρ P e , s	>	( ) ρ P e , i	and ID = 0 then
	← ( ) ρ P P e , s ;
	← ID 3;
	end
	if ( ) ρ P e , s	-	S min	>		( ) σρ e , c	and ID = 0 then
	← ( ) ρ P P e , s ;
	← ID 3;
	end
	if ( ) ρ P e , i	-	S min	>	( ) σρ e , c	and ID = 0 then
	P	← ( ) σ ρ + e , c	S min ;
	← ID 2;
	end
	if ( ) ρ P e ,	-	S	<	( ) σρ e ,

Table 1 .

 1 Main parameters used for computations with POREQST and HPOREQST performed in section 3.2.

	Dense graphite EOS		SESAME 7832
	Compaction curve		figure 1 a
	Pore reopening		-95 a
	parameter σ l (MPa)
	Hysteresis parameter		200
	σ h (MPa)
	Softening function f i (e)		Linear
	Initial densities (kg m -3 )	Porous ρ 0	1754 a
	Compressibility	Dense ρ s K 0	2250 c 8.58 a
	modulus (GPa)	K s	34.05 c
	Shear modulus (GPa)	G 0	5.01 a
		G s (P = 0)	20 a
		G s (P = 18 GPa)	50 a
		G s (P = 27.5 GPa) 280 a
	Elastic limit (GPa)	Y 0	0.1 b
		Y s (P = 0)	0.1 b
		Y s (P = 25 GPa)	5 b
	Limit distension α lim		1

a Data extracted from

[START_REF] Seisson | Experimental and theoretical study of the damaging of graphite under dynamic loading-application to hypervelocity impacts[END_REF]

. b Data from

[START_REF] Hébert | Simulations of hypervelocity impacts into graphite[END_REF]

. c Values deduced from the SESAME 7832 table.

Acknowledgments

The authors express their gratitude to the team of the Complexe 500 at CEA CESTA for having performed the plate impact experiments.