Hermine Biermé
email: hermine.bierme@univ-tours.fr

Agnès Desolneux
email: agnes.desolneux@ens-paris-saclay.fr

The anisotropy of 2D or 3D Gaussian random fields through their Lipschitz-Killing curvature densities

Keywords: 62M40, 62H11; 60G15, 60D05 Lipschitz-Killing curvatures, Gaussian random fields, Image analysis, Spatial statistics, Central limit theorems, Coarea formula

We are interested here in modeling and estimating the anisotropy of 2D and 3D Gaussian random fields through the geometry of their excursion sets. In order to do this, we use Lipschitz-Killing curvatures of the level sets as functions of the levels and see them as generalized processes for which we are able to obtain a joint functional Central Limit Theorem. For 2D and 3D stationary Gaussian fields we provide explicit formulas for the Lipschitz-Killing curvature densities from which we can deduce geometrical equivalent of second spectral moments and anisotropy ratios that allow the estimation of the anisotropy of the underlying Gaussian field.

Introduction

Since the seminal works of Robert Adler on the geometry of random fields [START_REF] Adler | The Geometry of Random Field[END_REF] and of Keith Worsley on topological inference in neuroimaging [START_REF] Worsley | The geometry of random images[END_REF], the study of the geometry of excursion sets of random fields has known a growing interest, with important developments especially for the Euler characteristic, used as a good approximation of the tail distribution of the supremum of stationary smooth Gaussian fields [START_REF] Worsley | A three-dimensional statistical analysis for CBF activation studies in human brain[END_REF]. In particular the so-called Gaussian kinematic formula [START_REF] Adler | Random fields and geometry[END_REF] nicely links mean geometries of an excursion set with observation window's ones with respect to the underlying statistical properties of the Gaussian field. However, this formula strongly depends on the metric induced by the stationary Gaussian field, and only mean Euler characteristic can be effectively computed in practice without any prior knowledge on the field. In contrast, in this paper we consider Lipschitz-Killing curvatures of excursion sets computed with the usual Euclidean metric in dimension d for both d = 2 and d = 3, corresponding to the usual dimensions of medical images. This allows us to rely on numerous results and algorithms developed in stochastic geometry for intrinsic volumes or Minkowski measures closely related to Lipschitz-Killing curvatures [START_REF] Rataj | Curvature measures of singular sets[END_REF].

Moreover, as in our previous paper [START_REF] Biermé | Mean Geometry for 2D random fields: level perimeter and level total curvature integrals[END_REF], we will have here a "weak" point of view. Instead of fixing a threshold level and considering the Lipschitz-Killing curvatures of the excursion set above this level, we will consider simultaneously all the levels and thanks to a change of variable formula (the coarea formula), we will be able to have a representation of these Lipschitz-Killing curvatures as integrals over the function domain. This point of view also allows us to work on the fine functional framework of generalized random processes [START_REF] Gel'fand | Generalized functions[END_REF]. In this setting, we are able to propose consistent and asymptotic (as the size of the observation window grows) Gaussian estimators of the Lipschitz-Killing curvature densities. Note that in sharp contrast with previous results, we do not assume isotropy of the Gaussian field.

Actually, a main point of interest is: how to "read" and "estimate" the anisotropy of a random field from the geometry of (some of) its excursion sets? We provide explicit parametric expressions for Lipschitz-Killing densities in both dimension d = 2 and d = 3, as well as numerical evaluations. It allows us to define new geometrical equivalent of spectral moments, related through elegant isoperimetric inequalities, as well as robust anisotropy ratios in the sense that they do not depend on the mean, nor on the standard deviation of the field that is crucial for image comparisons. As illustrated on Figure 1 we will summarize the geometry of an excursion set as a (2D here) point and this will allow us to visualize and estimate the anisotropy of the underlying Gaussian field. Two excursion sets of Gaussian random fields. The left one is isotropic and the right one is not. They have the same area, the same perimeter but different total curvature. Their geometry can be summarized as a point that is on the almond curve in the isotropic case and that is inside the almond domain in the anisotropic case.

The paper is organized as follows. In Section 2 we recall the main definitions of Lipschitz-Killing curvatures following [START_REF] Thäle | 50 years sets with positive reach-a survey[END_REF], and introduce the Integral Lipschitz-Killing curvatures of a smooth function that allow us to consider Lipschitz-Killing curvatures of its level sets as a tempered distribution (with respect to the levels). Considering a smooth stationary random field we are therefore able to define Lipschitz-Killing curvatures for which a functional joint Central Limit Theorem is established in dimension d = 2 or d = 3 under the additional assumption that the underlying field is Gaussian. In contrast with previous Gaussian results [START_REF] Estrade | A central limit theorem for the Euler characteristic of a Gaussian excursion set[END_REF][START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF] based on Itô-Wiener chaos expansion, our proof relies on preliminary general results stated under a quasi-association assumption closer to [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF], and which are postponed to the Appendix A. In Section 3, we focus on dimension d = 2 and extend the Gaussian isotropic results of [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF] to the anisotropic framework. We give explicit formulas for the Lipschitz-Killing curvature densities with respect to two important ordered geometrical equivalent of second spectral moments. This allows us to define an almond curve of anisotropy and a first ratio of anisotropy R with values in (0, 1]. We proceed in a simililar way in Section 4 for the dimension d = 3. Lipschitz-Killing curvature densities are also expressed with respect to now three important geometrical equivalent of second spectral moments. The ordering of these quantities is much more difficult to obtain than in the 2D case and strongly relies on some isoperimetric inequalities. In addition to the almond curve of anisotropy we define the heart curve of anisotropy involving the third index and get two respective ratios of anisotropy R and R G ≤ R, both with values in (0, 1]. We illustrate all this in both dimensions d = 2 and d = 3 through numerical experiments, checking the formulas for the Lipschitz-Killing curvatures, the Central Limit behavior, and the statistical estimation of the anisotropy from the geometry of an excursion set.

Lipschitz-Killing curvatures

Lipschitz-Killing curvatures : the general framework

Let us introduce some notations. In the following we will denote by L d or simply by L (when there is no ambiguity) the Lebesgue measure on R d and by H s the s-dimensional Hausdorff measure on R d . We have in particular

H d = L d on R d .
We now recall the definition of Lipschitz-Killing curvatures in the framework of smooth manifolds of R d , as they are presented in the paper of Thäle [START_REF] Thäle | 50 years sets with positive reach-a survey[END_REF]. Let M be a compact d-dimensional submanifold in R d with a C 2 smooth boundary ∂M (that is then (d-1)-dimensional). For x ∈ ∂M , let κ i (x), i = 1, . . . , d -1, denote the d -1 principal curvatures of ∂M at x. Then the Lipschitz-Killing curvatures of M are defined by ∀j = 0, . . . , d -1, C j (M ) = 1 α d-1-j ∂M σ d-1-j (κ 1 (x), . . . , κ d-1 (x)) H d-1 (dx), [START_REF] Adler | The Geometry of Random Field[END_REF] and

C d (M ) = L d (M ),
where σ k is the symmetric elementary function of order k, that is

σ k (κ 1 (x), . . . , κ d-1 (x)) = 1≤i1≤...≤i k ≤d-1 κ i1 (x) . . . κ i k (x),
with the convention that σ 0 = 1, and where α k , k ≥ 0 is the k-dimensional surface area of the unit ball of R k+1 , that is :

α 0 = 2, α 1 = 2π, α 2 = 4π, etc.
The Lipschitz-Killing curvatures have therefore the following interpretations:

• C 0 (M ) = 1 α d-1 ∂M κ 1 (x) . . . κ d-1 (x) H d-1 (dx) = χ(M )
is the Euler Characteristic of M , that is related to the integral of the Gaussian curvature on ∂M by the Gauss-Bonnet Theorem.

• C d-2 (M ) = 1 α1 ∂M (κ 1 (x) + . . . + κ d-1 (x)) H d-1 (dx). When d ≥ 3 it is equal to d-1 2π b(M )
where b(M ) is the mean breadth of M defined as the integral of the mean curvature on ∂M .

• C d-1 (M ) = 1 α0 ∂M H d-1 (dx) = 1 2 H d-1 (∂M ) is half the surface area (perimeter when d = 2) of ∂M in R d . • C d (M ) = L d (M ) is the volume of M in R d .
Let U be an open bounded set in R d . Then, the Lipschitz-Killing curvatures can be localized in U : they are then denoted by C j (M, U ) and their formula involving the principal curvatures are analogous to the ones of the C j (M ) given by [START_REF] Adler | The Geometry of Random Field[END_REF], except that the integrals on ∂M are replaced by integrals on ∂M ∩ U .

The Lipschitz-Killing curvatures are related to the volume of the r-parallel set to M through the Weyl's tube formula. This point will be discussed in the further remarks of Section 5.2.

The framework of smooth functions

Let f : R d → R be a C 2 function defined on R d , with here d ≥ 2. For t ∈ R, the excursion set (or level set) of f above the level t is given by

E f (t) := {x ∈ R d ; f (x) ≥ t}. Let U be an open bounded subset of R d , for instance U = (0, T ) d with T > 0. Let t ∈ R, and x ∈ U such that f (x) = t and ∇f (x) = 0. Then ν f (x) = ∇f (x) ∇f (x) ∈ S d-1 is the unit inner normal vector of E f (t) at x ∈ ∂E f (t).
Now in an open neighborhood of such x ∈ U we can compute Dν f using D 2 f , the Hessian matrix of f , and have

Dν f (x) = 1 ∇f (x) D 2 f (x) - 1 ∇f (x) 3 ∇f (x)(D 2 f (x)∇f (x)) t .
As the principal curvatures are the eigenvalues of the second fundamental form, it follows that the d symmetric functions of the curvatures σ d-1-j (κ 1 (x), . . . , κ d-1 (x)) may be explicitly computed with formulas involving f (x

) := f (x), ∇f (x), D 2 f (x) through σ d-1-j (κ 1 (x), . . . , κ d-1 (x)) = 1 ∇f (x) F d-1-j (f (x)), (2)
with in particular (see [START_REF] Spivak | A comprehensive introduction to differential geometry[END_REF] or [START_REF] Goldman | Curvature formulas for implicit curves and surfaces[END_REF] for instance)

F 0 (f (x)) = ∇f (x) , (3)
F 1 (f (x)) = -Tr(D 2 f (x)) + ∇f (x) t D 2 f (x)∇f (x) ∇f (x) 2 , (4)
and for d ≥ 3

F d-1 (f (x)) = (-1) d-1 ∇f (x) t adj(D 2 f (x))∇f (x) ∇f (x) d , (5)
where adj(D 2 f (x)) is the adjugate matrix of D 2 f (x) (that is the transpose of its cofactor matrix). Note that for d = 2, one can also express F 1 as

F 1 (f (x)) = - (∇f (x) ⊥ ) t D 2 f (x)∇f (x) ⊥ ∇f (x) 2 , ( 6
) with ∇f (x) ⊥ = (-∂ 2 f (x), ∂ 1 f (x)). Now, considering A = {x ∈ U ; ∇f (x) = 0}
, the coarea formula (see [START_REF] Evans | Measure theory and fine properties of functions[END_REF] p.112 for instance) states that

R H d-1 A ∩ f -1 (t) dt = A ∇f (x) dx = 0.
Hence, for a.e. t ∈ R and H d-1 almost every x ∈ ∂E f (t) ∩ U we have f (x) = t and ∇f (x) = 0. Under an additional assumption this will allow us to define for a.e. t ∈ R

C j (E f (t), U ) = 1 α d-1-j ∂E f (t)∩U 1 ∇f (x) F d-1-j (f (x))H d-1 (dx), such that t → C j (E f (t), U ) is L 1 -integrable (in
the sense that it coincides a.e. with a measurable integrable function on R) and therefore to consider it as a tempered distribution. More precisely, let us recall that the Schwartz space S = S(R) consists of infinitely differentiable functions h : R → R that are rapidly decreasing, that is, for all q ∈ N and k ∈ N,

h q,k = sup t∈R (1 + |t|) q |h (k) (t)| < ∞,
where h (k) denotes the derivative of order k. It is a real vector space and it is equipped with the topology given by the family of semi-norms • q,k , q ∈ N and k ∈ N. It is well known that S is a Fréchet space (i.e. a locally convex metrizable complete space). The space of tempered distributions S = S (R) is the topological dual of S, that is the space of continuous linear functionals on S. Now for h ∈ S we are interested in the Integral Lipschitz-Killing Curvatures of f defined by

LC f j (h, U ) := R h(t)C j (E f (t), U ) dt. ( 7
)
Proposition 1. Let f : R d → R be a C 2 function and note f (x) := f (x), ∇f (x), D 2 f (x) for x ∈ R d . Let U be a bounded open set of R d and assume that, for k ∈ {0, 1, d -1}, U |F k (f (x))|1 I ∇f (x) >0 dx < +∞, (8)
with F k given by (3), ( 4) and [START_REF] Berzin | Estimation of local anisotropy based on level sets[END_REF].

Then for j ∈ {0, d -2, d -1}, one has C j (E f (•), U ) ∈ L 1 (R). Moreover, for j ∈ {0, d -2, d -1, d}, one has C j (E f (•), U ) ∈ S with for all h ∈ S LC f j (h, U ) = C j (E f (•), U ), h = U L j (h, f (x))dx, (9)
where h ∈ S → L j (h, f (x)) ∈ R is linear continuous for all x ∈ U , and given by

L j (h, f (x)) = 1 α d-1-j h(f (x))F d-1-j (f (x))1 I ∇f (x) >0 forj = d, and
L d (h, f (x)) = R h(t)1 I f (x)≥t dt.
Proof. Let us recall the coarea formula for Lipschitz mappings (see [START_REF] Evans | Measure theory and fine properties of functions[END_REF] p.117 for instance): for any non-negative L d -measurable or any L d -integrable function g : R d → R, the function t

→ ∂E f (t)∩U g(x)H d-1 (dx) is L 1 -measurable and U g(x) ∇f (x) dx = R ∂E f (t)∩U g(x)H d-1 (dx) dt.
Hence for k ∈ {0, 1, d-1}, and h a continuous bounded function we can define measurable functions on U by

g k (x) = h(f (x)) F k (f (x)) ∇f (x) 1 I ∇f (x) >0 for x ∈ U.
Under the assumption [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF] we obtain that

R ∂E f (t)∩U |g k (x)|H d-1 (dx) dt ≤ h 0,0 U |F k (f (x))|1 I ∇f (x) >0 dx < +∞.
Moreover, using again the coarea formula with max(g k , 0) and -min(g k , 0) and subtracting we get

R ∂E f (t)∩U g k (x)H d-1 (dx) dt = U h(f (x))F k (f (x))1 I ∇f (x) >0 dx. It follows that t → ∂E f (t)∩U g k (x)H d-1 (dx) is L 1 -integrable.
g k (x)H d-1 (dx) = h(t) ∂E f (t)∩U σ k (κ 1 (x), . . . , κ d-1 (x))1 I ∇f (x) >0 H d-1 (dx) dt = h(t) ∂E f (t)∩U σ k (κ 1 (x), . . . , κ d-1 (x))H d-1 (dx) dt = α d-1-k C d-1-k (E f (t), U ) h(t), since for a.e. t ∈ R and H d-1 a.e. x ∈ ∂E f (t) ∩ U we have f (x) = t and ∇f (x) = 0. This implies (taking h = 1) that t → C j (E f (t), U ) is L 1 -integrable for j ∈ {0, d -2, d -1}
and in view of [START_REF] Biermé | Mean Geometry for 2D random fields: level perimeter and level total curvature integrals[END_REF],

LC f j (h, U ) = U L j (h, f (x))dx, with L j (h, f (x)) = 1 α d-1-j h(f (x))F d-1-j (f (x))1 I ∇f (x) >0 . Hence we have also h → L j (h, f (x)) ∈ S with |L j (h, f (x))| ≤ 1 α d-1-j h 0,0 |F d-1-j (f (x))|1 I ∇f (x) >0 .
Finally, for j = d we can also define

L d (h, f (x)) = R h(t)1 I f (x)≥t dt, satisfying |L d (h, f (x))| ≤ 2 h 2,0 , using the fact that R (1 + |t|) -2 dt = 2. Since U is bounded we also have h → LC f d (h, U ) = U L d (h, f (x)) dx ∈ S , with, by Fubini's Theorem, LC f d (h, U ) = R h(t)C d (E f (t), U ) dt.

Smooth stationary random fields and Lipschitz-Killing curvature densities

We extend now the Lipschitz-Killing curvatures integrals to the case of smooth stationary random fields.

Proposition 2. Let (Ω, A, P) be a complete probability space and (X(x)) x∈R d be a stationary real-valued random field that is almost surely (a.s.) a C 2 function on R d . Assume moreover that

∀k ∈ {0, 1, d -1}, E |F k (X(x))|1 I ∇X(x) >0 < +∞, ( 10
)
for F k given by (3), ( 4) and (5), and

X(x) = (X(x), ∇X(x), D 2 X(x)), x ∈ R d . Let U be a bounded open set of R d . Then we can define C j (E X (•), U ) as a generalized integrable process for j ∈ {0, d -2, d -1, d}, in the sense that C j (E X (•), U ) : (Ω, A) → (S , B(S ))
is measurable, where B(S ) is the Borel σ-field of the strong topology on S and, for all h ∈ S,

E |LC X j (h, U )| ≤ R |h(t)|E (|C j (E X (t), U )|) dt < +∞.
Proof. First note that since (Ω, A, P) is a complete probability space, in view of the a.s. continuity, up to setting the values of X at 0 on a negligible set, the field (h, U ) = 0 for ω ∈ N . Using the fact that (ω, x) → L j (h, X(ω, x)) is P ⊗ L d -measurable, it follows by Fubini's theorem that LC X j (h, U ) : ω ∈ Ω → LC X(ω) j

X : Ω × R d → R is a P ⊗ L d -
(h, U ) ∈ R is a well-defined real valued random variable. Therefore, (LC X j (h 1 , U ), . . . , LC X j (h m , U )) is a random vector for all m ≥ 1 and h 1 , . . . , h m ∈ S. Since the cylinder σ-field coincides with the Borel σ-field of the weak topology but also of the strong topology on S denoted by B(S ) (see Corollary 3.9 of [START_REF] Biermé | Generalized random fields and Lévy's continuity theorem on the space of tempered distributions[END_REF]) we can deduce that C j (E X (•), U ) : (Ω, A) → (S , B(S )) is measurable and therefore defines a generalized random process. To conclude for integrability, it is enough to remark that by Fubini's theorem and the coarea formula

E |LC X j (h, U )| ≤ R |h(t)|E (|C j (E X (t), U )|) dt ≤ h 0,0 α d-1-j U E |F d-1-j (X(x))|1 I ∇X(x) >0 dx < +∞,
by [START_REF] Bulinski | Limit theorems for associated random fields and related systems[END_REF], using the stationarity of X and the boundedness of U . The result for j = d follows the same lines as (ω, x) → L d (h, X(ω, x)) is P ⊗ L d -measurable and bounded by 2 h 2,0 .

Thanks to the stationarity we will focus on Lipschitz-Killing curvature densities defined for

j ∈ {0, d -2, d -1} by C X j , h = R h(t)C X j (t)dt = E LC X j (h, U ) L d (U ) = 1 α d-1-j E h(X(x))F d-1-j (X(x))1 I ∇X(x) >0 . (11)
Note that C X j ∈ L 1 (R) ⊂ S and that for a.e. t ∈ R we have

C X j (t) = E (C j (E X (t), U )) L d (U ) .
We can define in a simpler way

C X d (t) = E(C d (E X (t),U )) L d (U ) = P(X(x) ≥ t) for all t ∈ R, such that C X d ∈ S (but not in L 1 (R)) with C X d , h = R h(t)C X d (t)dt = E LC X d (h, U ) L d (U )
.

To ensure sample paths smoothness, a simple assumption for stationary Gaussian random fields may be given in term of covariance functions. More precisely, when X is a stationary Gaussian field with mean m and variance σ 2 , we denote by ρ its covariance function, given by ρ(x) = Cov(X(x), X(0)) = Cov(X(x + y), X(y)), for all x, y ∈ R d by stationarity. Then it is sufficient to make the following hypothesis, denoted by H : ρ : R d → R is a C 4 positive definite function, with det(D 2 ρ(0)) = 0 and such that there exist C > 0 and δ > 0 with

∂ k ρ(x) -∂ k ρ(0) ≤ C x δ , for all k = (k 1 , . . . , k d ) ∈ N d with |k| = d l=1 k l = 4. (H)
Corollary 1. Let (Ω, A, P) be a complete probability space and (X(x)) x∈R d be a stationary Gaussian random field whose covariance function satisfies H. Then, there exists a modification of X such that the assumptions of Proposition 2 are satisfied. Moreover for d ∈ {2, 3}, 0 ≤ j ≤ d and h ∈ S, the random variable LC X j (h, U ) is square integrable. Proof. By Proposition 2.1 of [START_REF] Cheng | The mean Euler characteristic and excursion probability of Gaussian random fields with stationary increments[END_REF], under H, there exists a modification of X such that X is C 2 a.s. Note that in this case ∇X(x) is an R d -valued Gaussian vector of covariance given by Γ = -D 2 ρ(0) and therefore ∇X(x) > 0 a.s. with E( ∇X(x) p(2-d) ) < +∞ for all p ≥ 1 such that p(d-2) < d. In view of (3), (4), or [START_REF] Berzin | Estimation of local anisotropy based on level sets[END_REF] we have therefore for all k ∈ {0, 1, d -1},

E (|F k (X(x))|) < +∞,
and (10) is satisfied. Moreover, for d ∈ {2, 3}, we also have E F k (X(x)) 2 < +∞, and LC X d-1-k (h, U ) is a square integrable random variable. The result for LC X d (h, U ) simply comes from the fact that C d (E X (t), U ) is a bounded random variable.

In this Gaussian setting, we obtain the following result, whose technical proof is postponed to Appendix A.

Theorem 2.1. Let X be a stationary Gaussian field defined on R d for d ∈ {2, 3}, whose covariance function satisfies H and assume moreover that there exists β > 12 + d such that

max |k|≤4 ∂ k ρ(x) ≤ C(1 + |x|) -β . ( 12
)
Then for U n = (0, n) d one has

L d (U n ) C j (E X (•), U n ) L d (U n ) -C X j 0≤j≤d distrib. -→ n→+∞ (B j ) 0≤j≤d ,
where (B j ) 0≤j≤d are centered Gaussian generalized random processes whose covariance functions are given for all h, h ∈ S and 0 ≤ j, k ≤ d, by

Cov B j , h , B k , h = R d Cov L j (h, X(z)), L k ( h, X(0)) dz.
Here the convergence holds in distribution with respect to the strong topology of S (R).

Let us mention that the decay assumption (12) could be relaxed into max |k|≤4 ∂ k ρ ∈ L 1 (R d ) for Gaussian fields, as assumed in [START_REF] Estrade | A central limit theorem for the Euler characteristic of a Gaussian excursion set[END_REF], but here the proof of our theorem relies on an intermediate result given under a more general assumption of quasi-associated field (not necessarily Gaussian). Now, let us be more explicit about the covariances involved in Theorem 2.1. Note first that for all h, h ∈ S,

Cov L d (h, X(z)), L d ( h, X(0)) = R R h(s) h(t)Cov(1 I X(z)≥s , 1 I X(0)≥t ) ds dt. Then, introducing for z ∈ R d , Σ st dd (z) := Cov(1 I X(z)≥s , 1 I X(0)≥t ), ( 13
) we obtain Cov L d (h, X(z)), L d ( h, X(0)) = R R h(s) h(t)Σ st dd (z)ds dt.
Furthermore, for the other covariances, we will also be able to write them as

Cov L j (h, X(z)), L k ( h, X(0)) = R R h(s) h(t)Σ st jk (z)ds dt.
Actually, by stationarity, for j < d,

E (L j (h, X(z))) = E (L j (h, X(0)) = 1 α d-1-j R h(t)E (F j (X(0))|X(0) = t) p X(0) (t)dt,
where p X(0) denotes the density of the Gaussian variable X(0), and

E L j (h, X(z))L d ( h, X(0))) = R h(t)E h(X(z))E F d-1-j (X(z)) α d-1-j 1 I X(0)≥t |X(z) dt = R 2 h(t)h(s) E F d-1-j (X(z)) α d-1-j 1 I X(0)≥t |X(z) = s p X(z) (s) ds dt.
Then we can set

Σ st jd (z) := E F d-1-j (X(z))1 I X(0)≥t |X(z) = s -E (F d-1-j (X(0))|X(0) = s) P(X(0) ≥ t) p X(0) (s) α d-1-j .
(14) And similarly, for k < d and z = 0, we set

Σ st jk (z) := E F d-1-j (X(z)) α d-1-j F d-1-k (X(0)) α d-1-k |X(0) = t, X(z) = s p (X(0),X(z)) (t, s) -E F d-1-j (X(0)) α d-1-j |X(0) = s E F d-1-k (X(0)) α d-1-k |X(0) = t p X(0) (t)p X(0) (s), (15)
denoting by p (X(0),X(z)) the joint density of (X(0), X(z)). Note that [START_REF] Cheng | The mean Euler characteristic and excursion probability of Gaussian random fields with stationary increments[END_REF] implies that ρ is an integrable function and therefore X admits a spectral density, thus implying |ρ(z)| < ρ(0) = σ 2 for all z = 0. This ensures the existence of p (X(0),X(z)) as soon as z = 0. Now, if we could justify the interchange of integrals we should also obtain

Cov B j , h , B k , h = R R h(t)h(s)Σ jk (s, t) ds dt, with Σ jk (s, t) = R d Σ st jk (z)dz,
allowing us to give a pointwise representation of the Gaussian generalized processes obtained at the limit. This is precisely the purpose of the following corollary.

Corollary 2. Under the assumptions of Theorem 2.1, assuming moreover that for J ⊂ {0, . . . , d} and V an open interval of R, we have:

(A1) The map t ∈ V → C X j (t, U n ) ∈ L 2 (Ω, A, P) is continuous for all n ≥ 1 and j ∈ J; (A2) For all ε > 0, there exists δ > 0 such that for all t, s ∈ V with |t -s| ≤ δ, for all n ≥ 1 and j ∈ J, Var C X j (t, U n ) -C X j (s, U n ) ≤ εL d (U n ); (A3) For all j, k ∈ J, for all (s, t) ∈ V × V we have R d |Σ st jk (z)|dz < +∞, and the function (s, t) ∈ V × V → Σ jk (s, t) := R d Σ st jk (z)dz is continuous. Then L d (U n ) C j (E X (t, U n )) L d (U n ) -C X j (t) j∈J,t∈V f dd -→ n→+∞ Bj (t) j∈J,t∈V
, where Bj (t)

j∈J,t∈V is a centered Gaussian process with covariance given by Cov( Bj (s), Bk (t)) = Σ jk (s, t).

Proof. For sake of simplicity, we only sketch the proof for k = j ∈ J and s = t ∈ V . Let us first choose h ∈ S non-negative with compact support such that R h(u)du = 1, as well as (h m ) m≥1 its associated approximation of the identity, given by h m (u) = mh(mu). We set τ t ȟ(u) = h(t -u) and may assume that τ t ȟm has support in V for all m ≥ 1. Then for all m, p ≥ 1,

Cov B j , τ t ȟm , B j , τ t ȟp = R R τ t ȟm (u)τ t ȟp (v)Σ jj (u, v) du dv -→ m,p→+∞
Σ jj (t, t), using Fubini's theorem and continuity, thanks to (A3). Then B j , τ t ȟm m is a Cauchy sequence of Gaussian variables in L 2 and we can set Bj (t) its Gaussian limit. Note that therefore we also

have Var Bj (t)) = Σ jj (t, t). Moreover, since C j (E X (t, U n )) ∈ L 2 (Ω, A, P) by (A1), we have E   L d (U n ) C j (E X (t, U n )) L d (U n ) -C X j (t) -L d (U n ) LC X j (τ t ȟm , U n ) L d (U n ) -C X j , τ t ȟm 2   = 1 L d (U n ) Var C j (E X (t), U n )) -LC X j (τ t ȟm , U n ) . But C j (E X (t), U n ) -LC X j (τ t ȟm , U n ) = R h(u) C j (E X (t), U n ) -C j E X t - u m , U n du. Therefore, Var C j (E X (t), U n )) -LC X j (τ t ȟm , U n ) = R 2 h(u)h(v)... Cov C j (E X (t), U n ) -C j E X t - u m , U n , C j (E X (t), U n ) -C j E X t - v m , U n dudv.
By (A2), for ε > 0 and u, v with |u| ≤ δm and |v| ≤ δm, and by Cauchy-Schwarz inequality, one has

Cov C j (E X (t), U n ) -C j E X t - u m , U n , C j (E X (t), U n ) -C j E X t - v m , U n ≤ εL d (U n ).
Hence, for m large enough such that h m has compact support in [-δ, δ] one has,

1 L d (U n ) Var C j (E X (t), U n )) -LC X j (τ t ȟm , U n ) ≤ ε.
Hence, choosing m such that we also have

Var B j , τ t ȟm -Bj (u) ≤ ε,
the results follows from Theorem 2.1, using the fact that we have

L(U n ) LC X j (τ t ȟm , U n ) L d (U n ) -C X j , τ t ȟm d -→ n→+∞ B j , τ t ȟm .
Note that for J = {d}, under the assumptions of Theorem 2.1, the random field X satisfies also the assumptions of Corollary 2 on V = R and we recover results of Theorem 2 of [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF]. Actually, z → Σ st dd (z) is continuous on R d and for all z ∈ R d , by Lemma 2 of [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF], we also have

Σ st dd (z) ≤ 1 4σ 2 |ρ(z)|, such that Σ dd (s, t) = R d Σ st dd (z)
dz is well-defined and Lebesgue's theorem allows to check (A3) on R.

For J = {d -1} and d = 2, our result corresponds to the one obtained in Theorem 1 of [START_REF] Iribarren | Asymptotic behaviour of the integral of a function on the level set of a mixing random field[END_REF], under a general mixing condition, in Theorem 3 of [START_REF] Kratz | Central limit theorems for level functionals of stationary Gaussian processes and fields[END_REF] or in Theorem 4.7 of [START_REF] Berzin | Estimation of local anisotropy based on level sets[END_REF]. For a general j ∈ {0, . . . , d}, we also refer to [START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF] who proved marginal Central limit theorems under an assumption of isotropy that we have relaxed here. We also obtain a stronger result with a joint Central limit theorem.

In the sequel we will compute explicitly Lipschitz-Killing densities C X j (t) at some level t ∈ R and will estimate them using

Cj (E X (t,U )) L d (U )
for an observation window U large enough.

Smooth stationary 2D Gaussian fields

Lipschitz-Killing curvature densities

When the dimension is d = 2, the quantities of interest are the Lipschitz-Killing curvatures C 0 , C 1 and C 2 that are here, up to multiplicative constants, the Total Curvature (related to the Euler Characteristic), the Perimeter, and the Area of the excursion sets of X in an open bounded domain U . More precisely we define for a smooth stationary field X defined on R 2 and h ∈ S

LA X (h, U ) := LC X 2 (h, U ), LP X (h, U ) := 2 LC X 1 (h, U )
and LTC X (h, U ) := 2π LC X 0 (h, U ). Now we give explicit formulas for the Lipschitz-Killing curvature densities of a stationary C 2 Gaussian field. Theorem 3.1. Let X be a stationary Gaussian field of mean m ∈ R, variance σ 2 > 0 and whose covariance function ρ satisfies Hypothesis H. We write γ 1 , γ 2 > 0 the eigenvalues of Γ = -D 2 ρ(0) and denote by Φ the standard Gaussian cumulative distribution. Then the Lipschitz-Killing curvature densities of X are given for a.e. t ∈ R by

A X (t) := C X 2 (t) = 1 -Φ t -m σ , Per X (t) := 2C X 1 (t) = γ Per σ 2 1 2 e -(t-m) 2 /2σ 2 TC X (t) := 2πC X 0 (t) = γ TC σ 2 1 √ 2π t -m σ e -(t-m) 2 /2σ 2 where γ Per = 1 2π 2π 0 γ 1 cos 2 θ + γ 2 sin 2 θ dθ 2 , (16)
and

γ TC = √ γ 1 γ 2 . ( 17
)
Therefore, the Gaussian random field X has same Perimeter density as a (weakly) isotropic Gaussian field with second spectral moment γ Per and same Total Curvature density as a (weakly) isotropic Gaussian field with second spectral moment γ TC . Moreover, we have

min(γ 1 , γ 2 ) ≤ γ TC ≤ γ Per ≤ max(γ 1 , γ 2 ),
where inequalities are strict if and only if

γ 1 = γ 2 .
Proof. Let h ∈ S. For j = 2, the result simply comes from the fact that C X 2 (t) = P(X(0) ≥ t) with X(0) ∼ N (m, σ 2 ). For j ∈ {0, 1}, by [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF], since ∇X(0) > 0 a.s. by H and since the random variable X(0) admits a probability density given by the function p X(0

) (t) = 1 σ √ 2π e -(t-m) 2 /2σ 2 , we have α 1-j C X j , h = E (h(X(0))F 1-j (X(0))) = R h(t)E (F 1-j (X(0))|X(0) = t) p X(0) (t)dt.
Hence, for a.e. t ∈ R, using (3) and ( 6),

C X 1 (t) = 1 2 E( ∇X(0) X(0) = t)p X(0) (t) C X 0 (t) = - 1 2π E (∇X(0) ⊥ ) t D 2 X(0)∇X(0) ⊥ ∇X(0) 2 |X(0) = t p X(0) (t).
In order to go ahead in the computations, we first remark that writing X = m + σY • P , for P orthogonal such that Γ = P t ∆P with ∆ = diag(γ 1 , γ 2 ), one has therefore

C X j (t) = C Y j ( t-m
σ ) and Y is centered with unit variance and admits γ1 σ 2 , γ2 σ 2 for spectral moments with (denoting by Y i and Y ij the first and second order partial derivatives of Y ):

γ 1 σ 2 = E(Y 1 (0) 2 ), γ 2 σ 2 = E(Y 2 (0) 2 ), and E(Y 1 (0)Y 2 (0)) = 0.
Since Y is stationary Gaussian, we also have that ∇Y (0) is independent from Y (0) and D 2 Y (0), and that Y (0) and D 2 Y (0) are correlated with covariance

E(Y (0)Y 11 (0)) = - γ 1 σ 2 , E(Y (0)Y 22 (0)) = - γ 2 σ 2 , and E(Y (0)Y 12 (0)) = 0.
Therefore the conditional expectations are given by

E(Y 11 (0)|Y (0) = t) = - γ 1 σ 2 t, E(Y 22 (0)|Y (0) = t) = - γ 2 σ 2 t, and E(Y 12 (0)|Y (0) = t) = 0.
Now, using that for a vector v ∈ R 2 , we have, denoting e θ = (cos θ, sin θ),

v = 1 4 2π 0 | v, e θ | dθ,
and using that for a 1D Gaussian random variable Z of mean 0 and variance σ 2 Z , we have

E(|Z|) = 2 π σ Z , we get C Y 1 (t) = 1 2 E( ∇Y (0) )p Y (0) (t) = 1 8πσ e -t 2 /2 2π 0 γ 1 cos 2 θ + γ 2 sin 2 θ dθ.
To compute C Y 0 (t), we first need to compute E(cos 2 Θ) where Θ is the angle of ∇Y (0), i.e. ∇Y (0) = ( ∇Y (0) cos Θ, ∇Y (0) sin Θ). Since Y is not assumed isotropic, the law of Θ is not necessarily uniform on [0, 2π). Indeed we have

E(cos 2 Θ) = σ 2 2π √ γ 1 γ 2 y 2 1 y 2 1 + y 2 2 e -(σy1) 2 /2γ1 e -(y2σ) 2 /2γ2 dy 1 dy 2 = 1 2π 2π 0 γ 1 cos 2 θ γ 1 cos 2 θ + γ 2 sin 2 θ dθ. ( 18
)
It follows that

E(sin 2 Θ) = 1 2π 2π 0 γ 2 sin 2 θ γ 1 cos 2 θ + γ 2 sin 2 θ dθ.
Finally,

C Y 0 (t) = - 1 2π E(sin 2 Θ)E(Y 11 (0)|Y (0) = t) + E(cos 2 Θ)E(Y 22 (0)|Y (0) = t) p Y (0) (t),
and combining everything, we get

C Y 0 (t) = t σ 2 ( √ 2π) 5 e -t 2 /2 2π 0 γ 1 γ 2 γ 1 cos 2 θ + γ 2 sin 2 θ dθ = t σ 2 ( √ 2π) 3 e -t 2 /2 √ γ 1 γ 2 .
The last equality comes from the following computation:

1 2π 2π 0 γ 1 γ 2 γ 1 cos 2 θ + γ 2 sin 2 θ dθ = 2 π π/2 0 γ 1 γ 2 γ 1 cos 2 θ + γ 2 sin 2 θ dθ = 2 π +∞ 0 γ 1 γ 2 γ 1 + γ 2 s 2 ds = √ γ 1 γ 2 , ( 19
)
where we have used a change of variable s = tan θ, and then recognized the derivative of the function s → arctan( γ2 γ1 s). The remark on the "equivalent" isotropic fields that have same Perimeter density or same Total curvature density as X comes from the fact that when γ 1 = γ 2 = γ, then the above computations yield γ Per = γ TC = γ.

To prove the inequality between γ TC and γ Per , let us first notice that

γ Per = 1 2π 2π 0 1 2 γ 1 cos 2 θ + γ 2 sin 2 θ + γ 1 sin 2 θ + γ 2 cos 2 θ dθ.
Then, for any s ∈ [0, 1], we have

1 2 γ 1 s + γ 2 (1 -s) + γ 1 (1 -s) + γ 2 s ≥ ((γ 1 s + γ 2 (1 -s))(γ 1 (1 -s) + γ 2 s)) 1/4 = γ 2 1 s(1 -s) + γ 2 2 s(1 -s) + γ 1 γ 2 (s 2 + (1 -s) 2 ) 1/4 = γ 1 γ 2 + s(1 -s)(γ 1 -γ 2 ) 2 1/4 ≥ (γ 1 γ 2 ) 1/4 ,
and this proves the inequality, by simply setting s = cos 2 θ and then integrating over θ. We also remark that this inequality is strict if and only if

γ 1 = γ 2 .
Note that we recover our previous results of [START_REF] Biermé | Mean Geometry for 2D random fields: level perimeter and level total curvature integrals[END_REF] in the (weakly) isotropic case. Moreover, in the anisotropic case, assuming that γ 2 > γ 1 and considering k = 1 -γ 1 /γ 2 , one can write

γ Per = √ γ 2 2 π π/2 0 (1 -k 2 ) cos 2 θ + sin 2 θ dθ = √ γ 2 2 π E(k),
where E(k) is the elliptic integral of the first kind, such that Per X (t) corresponds to the formula of Theorem 11.3. of [START_REF] Azais | Level sets and extrema of random processes and fields[END_REF]. We also refer to Proposition 3.4 of [START_REF] Estrade | Anisotropic Gaussian wave models[END_REF] for planar Gaussian waves.

Visualization and estimation of the anisotropy

In this section, we are interested in visualizing and estimating the anisotropy of a Gaussian field from the Lipschitz-Killing curvatures of its excursion sets. Here the anisotropy is understood in a weak way, and it is summarized through the anisotropy ratio R defined by

R := γ TC γ Per .
Thanks to Theorem 3.1, we have R ∈ [0, 1], and R = 1 if and only if the spectral moments γ 1 and γ 2 are equal, meaning that the field X is (weakly) isotropic. More precisely, we can see R as a function of α := min(γ1,γ2) max(γ1,γ2) ∈ [0, 1] since

R = γ TC γ Per = √ γ 1 γ 2 1 2π 2π 0 γ 1 cos 2 θ + γ 2 sin 2 θ dθ 2 = √ α 1 2π 2π 0 α cos 2 θ + sin 2 θ dθ 2 .
This is illustrated on the left of Figure 2 where we show the graph of the function α → R. It is an increasing function that maps [0, 1] to [0, 1].

The Almond curve of anisotropy

Inspired by the paper of Klatt et al. [START_REF] Klatt | Characterization of anisotropic Gaussian random fields by Minkowski tensors[END_REF] for the visualization of anisotropy, from the formulas for the Perimeter density and the Total Curvature density (Theorem 3.1), we see that if we set

x(t) = Per X (t) Per X (m) = e -(t-m) 2 /2σ 2 and ỹ(t) = TC X (t) Per X (m) 2 = γ TC γ Per 4 √ 2π t -m σ e -(t-m) 2 /2σ 2 ,
then the point (x(t), ỹ(t)) is on the curve

C R = {(x, y) ∈ R 2 ; y 2 + 16 π R 2 x 2 log x = 0}.
Notice that this curve is in the sub-domain x ∈ (0, 1] and that, thanks to the fact that R ≤ 1, it is inside the domain defined by the closed curve

C 1 = {(x, y) ∈ R 2 ; y 2 + 16 π x 2 log x = 0}.
This curve is "canonical" in the sense that it is independent of the mean and of the variance of the field. We call this curve the almond curve of anisotropy because of its shape, see Figure 2 right. Points on this curve correspond to R = 1, that is γ TC = γ Per and thus γ 1 = γ 2 , meaning the field X is (weakly) isotropic.

Estimating the anisotropy from one excursion set

Assume we observe an excursion set Ê of a sample of a Gaussian stationary random field X, in a window U . The question is then: how to estimate the anisotropy ratio R from this set Ê only, that is without having access to the whole field (nor its increments for instance). We ignore also at which level t the field X was thresholded, and also what was its mean m or its variance σ 2 . Now from the three observed Lipschitz-Killing curvatures densities (that are Â(E) = A( Ê)/L 2 (U ) the normalized Area of Ê, Per(E) = Per( Ê)/L 2 (U ) the normalized Perimeter of Ê, and TC(E) = TC( Ê)/L 2 (U ) the normalized Total Curvature of E), we can proceed the following way:

-Compute an "effective" threshold t, that is an empirical estimate of t-m σ , using the density area of Ê:

t = Φ -1 (1 -Â(E)).
-Then estimate Per X (m) by p m := Per(E)e t 2 /2 . -And the "effective" spectral moments via

γ Per σ 2 = (2 p m ) 2 , γ TC σ 2 = TC(E) √ 2π t e t 2 /2 .
Finally estimate the anisotropy ratio with

R = γ TC γ Per = TC(E) p m 2 √ 2π 4 t e t 2 /2 .
This estimation procedure is inspired from the one developed by Elena Di Bernardino and Céline Duval in [START_REF] Di Bernardino | Statistics for Gaussian random fields with unknown location and scale using Lipschitz-Killing curvatures[END_REF], where they introduce the notion of "effective" threshold, but where the framework is quite different since they assume the Gaussian field is isotropic.

Note that only "effective" spectral moments γ Per σ 2 and γ TC σ 2 are empirically accessible in practice when X is not assumed standard (i.e. so that σ = 1) but our anisotropic ratio R itself is empirically accessible. Moreover we can also compute a point on the almond curve (see Figure 3) by setting x =

Per(E) p m = e -t 2 /2 and ŷ = TC(E)

p m 2 = R 4 √ 2π te -t 2 /2 .
Let us remark that, assuming Ê is the realization of E X (t)∩U n for some t ∈ R and n large, under the assumptions of Corollary 2 on an open neigborhood V of t and J = {0, 1, 2}, Â(E), Per(E) and TC(E) are consistent and asymptotically Gaussian. When moreover t = m, by the delta method, it should imply the same properties for the effective spectral moments and anisotropy ratio (see Figure 8 in the additional experiments section B in the Appendix, where we also numerically check the central limit behavior of the anisotropy ratio). This restriction on t means that we shall not have Â(E) close to 1/2 = Φ(0).

On Figures 3 and4 (and also on Figure 9 in the additional experiments section B in the Appendix), we illustrate all the theoretical results of this section. The numerical computations of the Lipschitz-Killing curvatures of the excursion sets are done thanks to the Matlab "imMinkowski" library1 developed by David Legland and described in [START_REF] Legland | Computation of Minkowski measures on 2D and 3D binary images[END_REF]. As far as binary images are concerned, it will be more convenient to work with "quantiles" q ∈ [0, 1] so that the corresponding threshold levels are given by t = m + σΦ -1 (q).

Smooth stationary 3D Gaussian fields

Lipschitz-Killing curvature densities

We are here interested in the 3D framework, and we consider a real-valued random field X defined on R 3 . The quantities of interest are the Lipschitz-Killing curvatures C 0 , C 1 , C 2 and C 3 that are now related to the Total Gaussian Curvature (also related to the Euler Characteristic by the Gauss-Bonnet Theorem), the Total Mean Curvature (related to the mean breadth), the Surface Area and the Volume of an excursion set. More precisely, for a smooth stationary field X defined on R 3 and h ∈ S we denote

LV X (h, U ) := LC X 3 (h, U ), LSA X (h, U ) := 2LC X 2 (h, U ), LTMC X (h, U ) := πLC X 1 (h, U ) and LTGC X (h, U ) := 4πLC X 0 (h, U ).
Theorem 4.1. Let X be a stationary Gaussian field of mean m ∈ R, variance σ 2 > 0 and whose covariance function ρ satisfies Hypothesis H. We write γ 1 , γ 2 , γ 3 > 0 the eigenvalues of Γ = -D 2 ρ(0) and denote by Φ the standard Gaussian cumulative distribution. Then the Lipschitz-Killing curvature densities of X are given for a.e. t ∈ R by

V X (t) := C X 3 (t) = 1 -Φ t -m σ , SA X (t) := 2C X 2 (t) = γ SA σ 2 2 π e -(t-m) 2 /2σ 2 , TMC X (t) := πC X 1 (t) = γ TMC σ 2 1 √ 2π t -m σ e -(t-m) 2 /2σ 2 , TGC X (t) := 4πC X 0 (t) = γ TGC σ 2 3/2 1 π (t -m) 2 σ 2 -1 e -(t-m) 2 /2σ 2 ,
where

γ SA = 1 4π S 2 γ 1 u 2 1 + γ 2 u 2 2 + γ 3 u 2 3 H 2 (du) 2 , (20)
γ TMC = 1 2 (γ 1 + γ 2 + γ 3 ) - 1 8π S 2 γ 2 1 u 2 1 + γ 2 2 u 2 2 + γ 2 3 u 2 3 γ 1 u 2 1 + γ 2 u 2 2 + γ 3 u 2 3 H 2 (du), (21)
and

γ TGC = (γ 1 γ 2 γ 3 ) 1/3 . (22)
Therefore, the Gaussian random field X has same Surface Area density as a (weakly) isotropic Gaussian field with second spectral moment γ SA , same Total Mean Curvature density as a (weakly) isotropic Gaussian field with second spectral moment γ TMC and same Total Gaussian Curvature density as a (weakly) isotropic Gaussian field with second spectral moment γ TGC . Moreover, we have

min(γ 1 , γ 2 , γ 3 ) ≤ (γ 1 γ 2 γ 3 ) 1/3 = γ TGC ≤ γ TMC ≤ γ SA ≤ γ 1 + γ 2 + γ 3 3 ≤ max(γ 1 , γ 2 , γ 3 ).
where the inequalities are strict if and only if min(γ 1 , γ 2 , γ 3 ) < max(γ 1 , γ 2 , γ 3 ).

Proof. As in the 2D case, we may assume that m = 0, σ 2 = 1 and Γ = Diag(γ 1 , γ 2 , γ 3 ) meaning that ∀j = 1, 2, 3, γ j := Γ jj = E(X j (0) 2 ), and E(X i (0)X j (0)) = 0 when i = j,

denoting X i := ∂X ∂xi , X ij := ∂ 2 X ∂xi∂xj for 1 ≤ i, j ≤ d.
The formula for V X is straightforward. Now for j < 3, by Hypothesis H, we have ∇X(0) > 0 a.s. and E( ∇X(0) -1 ) < +∞. Since D 2 X(0) is independent from ∇X(0), recalling (3), (4) (5), it follows that the assumption (10) is satisfied. Then, for h ∈ S, by [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF], since the standard Gaussian random variable X(0) admits a probability density given by a function p X(0) , we have

α 2-j C X j , h = E (h(X(0))F 2-j (X(0))) = R h(t)E (F 2-j (X(0))|X(0) = t) p X(0) (t)dt.
Hence, for a.e. t ∈ R,

C X 2 (t) = 1 2 E( ∇X(0) X(0) = t)p X(0) (t) C X 1 (t) = 1 2π E ∇X(0) t D 2 X(0)∇X(0) ∇X(0) 2 -Tr D 2 X(0) X(0) = t p X(0) (t). C X 0 (t) = 1 4π E ∇X(0) t adj(D 2 X(0))∇X(0) ∇X(0) 3 X(0) = t p X(0) (t).
Note that since X is Gaussian, we have that ∇X(0) is independent from X(0) and D 2 X(0), and that X(0) and D 2 X(0) are correlated with covariance ∀j = 1, 2, 3, E(X(0)X jj (0)) = -γ j , and E(X(0)X ij (0)) = 0 when i = j. ∀j = 1, 2, 3, E(X(0)X jj (0)) = -γ j , and E(X(0)X ij (0)) = 0 when i = j.

Therefore the conditional expectations are given by E(X jj (0)|X(0) = t) = -γ j t, and E(X ij (0)|X(0) = t) = 0.

Moreover, we have

∂ 4 iijj ρ(0) = E(X ii (0)X jj (0)) = E(X ij (0) 2 ). Since by assumption D 2 ρ(0) = -Γ is diagonal, we have E(X ij (0)X(0)) = 0 and therefore E(X ij (0) 2 |X(0) = t) = ∂ 4 iijj ρ(0).
Using the formulas for the conditional distribution of Gaussian vectors, the covariance of (X ii (0), X jj (0)) knowing X(0) = t is given by

E(X ii (0)X jj (0)|X(0) = t) = γ i γ j t 2 + ∂ 4 iijj ρ(0) -γ i γ j .
Therefore, putting all together, we get

E(X ii (0)X jj (0) -X ij (0) 2 |X(0) = t) = γ i γ j (t 2 -1),
which doesn't require the fourth-order spectral moment anymore. We are now in position to compute the different densities. We start with C X 2 (t), which is equal to

C X 2 (t) = 1 2 E( ∇X(0) |X(0) = t)p X (t) = 1 2 E( ∇X(0) )p X(0) (t).
We use here the property that if

v ∈ R 3 , then v = 1 2π S 2 | v, u | H 2 (du).
Now, as in the 2D case, using the fact that the expectation of the absolute value of a centered Gaussian real random variable is 2 π times its standard deviation, we get, denoting u

= (u 1 , u 2 , u 3 ) ∈ S 2 , C X 2 (t) = 1 4π 2 e -t 2 /2 S 2 γ 1 u 2 1 + γ 2 u 2 2 + γ 3 u 2 3 H 2 (du) = 1 π γ SA e -t 2 /2 .
Then, for C X 1 , using the conditional expectations and the fact that ∇X(0) is independant from (X(0), D 2 X(0)), we have 2πC

X 1 (t) = tp X(0) (t) γ 1 + γ 2 + γ 3 -E γ 1 X 1 (0) 2 + γ 2 X 2 (0) 2 + γ 3 X 3 (0) 2 X 1 (0) 2 + X 2 (0) 2 + X 3 (0) 2 .
Using the law of the X i (0) and a change of variable with spherical coordinates, we get that

E γ 1 X 1 (0) 2 + γ 2 X 2 (0) 2 + γ 3 X 3 (0) 2 X 1 (0) 2 + X 2 (0) 2 + X 3 (0) 2 = 1 4π S 2 γ 2 1 u 2 1 + γ 2 2 u 2 2 + γ 2 3 u 2 3 γ 1 u 2 1 + γ 2 u 2 2 + γ 3 u 2 3 H 2 (du).

Hence, 2πC

X

1 (t) = 2γ TMC tp X(0) (t).
Finally, by independence of ∇X(0) from (X(0), D 2 X(0)), and the conditional expectations computed above, we get (where the indices i, j, k below are all distincts):

4πC X 0 (t) = p X(0) (t) 3 i=1 E X i (0) 2 (X 1 (0) 2 + X 2 (0) 2 + X 3 (0)) 3/2 E(X jj (0)X kk (0) -X jk (0) 2 |X(0) = t) = (t 2 -1)p X(0) (t) 3 i=1 E X i (0) 2 (X 1 (0) 2 + X 2 (0) 2 + X 3 (0)) 3/2 γ j γ k = (t 2 -1)p X(0) (t) × 1 (2π) 3/2 S 2 γ 1 γ 2 γ 3 (γ 1 u 2 1 + γ 2 u 2 2 + γ 3 u 2 3 ) 3/2 H 2 (du) = γ 3/2 TGC 1 π (t 2 -1)e -t 2 /2 .
The last line above comes from the following computation. Let J denote the above integral on S 2 , that can be also written as

J := 2 π/2 0 2π 0 γ 1 γ 2 γ 3 (γ 1 sin 2 ϕ cos 2 θ + γ 2 sin 2 ϕ sin 2 θ + γ 3 cos 2 ϕ) 3/2 sin ϕ dθdϕ.
Then by a change of variable s = tan ϕ, we have ds = (1+s 2 )dϕ, cos 2 ϕ = 1 1+s 2 and sin ϕ = s √ 1+s 2 . Therefore

J = 2 +∞ 0 2π 0 γ 1 γ 2 γ 3 s ((γ 1 cos 2 θ + γ 2 sin 2 θ)s 2 + γ 3 ) 3/2 dθ ds = 2 √ γ 3 2π 0 γ 1 γ 2 γ 1 cos 2 θ + γ 2 sin 2 θ dθ = 4π √ γ 1 γ 2 γ 3 ,
where we have used the computation made in the 2D case.

To prove the inequalities, let us first rewrite things using the random variable U = (U 1 , U 2 , U 3 ) ∈ R 3 that is assumed to be uniformly distributed on the unit sphere S 2 . Then

γ SA = E γ 1 U 2 1 + γ 2 U 2 2 + γ 3 U 2 3 2 , γ TMC = 1 2 E γ 1 (γ 2 + γ 3 )U 2 1 + γ 2 (γ 1 + γ 3 )U 2 2 + γ 3 (γ 1 + γ 2 )U 2 3 γ 1 U 2 1 + γ 2 U 2 2 + γ 3 U 2 3
.

By concavity of the square root function, we first have that

γ SA ≤ E(γ 1 U 2 1 + γ 2 U 2 2 + γ 3 U 2 3 ) = 1 3 (γ 1 + γ 2 + γ 3 ), because E(U 2 i ) = 1/3 (
this comes from the symmetry of the U i 's and the fact that U 2 1 +U 2 2 +U 2 3 = 1 a.s.). Note also that the inequality is strict for min(γ 1 , γ 2 , γ 3 ) < max(γ 1 , γ 2 , γ 3 ) by strict concavity of the square root function.

The inequality between γ SA and γ TGC = (γ 1 γ 2 γ 3 ) 1/3 can also be obtained in a straightforward way. Indeed, using again the concavity of the square root function, and the fact that i U 2 i = 1, we have

γ 1 U 2 1 + γ 2 U 2 2 + γ 3 U 2 3 ≥ √ γ 1 U 2 1 + √ γ 2 U 2 2 + √ γ 3 U 2 3 .
And then taking the expectation on both sides we have

E γ 1 U 2 1 + γ 2 U 2 2 + γ 3 U 2 3 ≥ 1 3 ( √ γ 1 + √ γ 2 + √ γ 3 ) ≥ (γ 1 γ 2 γ 3 ) 1/6 ,
thanks to the arithmetic-geometric inequality. Therefore

γ SA ≥ γ TGC .
The inequalities for γ TMC are more difficult to obtain. We first introduce the function

g γ : x ∈ R 3 → γ 1 x 2 1 + γ 2 x 2 2 + γ 3 x 2 3
and note that it is a twice continuously differentiable positive function on R 3 {0} with gradient given by ∇g γ (x) = 1 gγ (x) (γ 1 x 1 , γ 2 x 2 , γ 3 x 3 ). It follows that its directional gradient on the sphere at a point x ∈ S 2 is given by the orthonormal projection of ∇g γ (x) onto the tangent space at x that is precisely x ⊥ , namely

∇ S 2 g γ (x) = ∇g γ (x) -∇g γ (x), x x = ∇g γ (x) -g γ (x)x.
According to Poincaré -Wirtinger inequality (Theorem 5.4.1 of [START_REF] Groemer | Geometric applications of Fourier series and spherical harmonics[END_REF]) we have

Var (g γ (U )) ≤ 1 2 E ∇ S 2 g γ (U ) 2 .
But we remark that γ SA = E (g γ (U )) 2 and that E g γ (U )

2 = 1 3 (γ 1 + γ 2 + γ 3 ), while E ∇ S 2 g γ (U ) 2 = E ∇g γ (U ) 2 -E g γ (U ) 2 .
Then Poincaré -Wirtinger inequality rewrites as

3 2 E g γ (U ) 2 -γ SA ≤ 1 2 E ∇g γ (U ) 2 . But 3 2 E g γ (U ) 2 - 1 2 E ∇g γ (U ) 2 = 1 2 E γ 1 + γ 2 + γ 3 - γ 2 1 U 2 1 + γ 2 2 U 2 2 + γ 2 3 U 2 3 γ 1 U 2 1 + γ 2 U 2 2 + γ 3 U 2 3 = γ TMC .
Hence we have obtained γ TMC ≤ γ SA . Note also that the case of equality in Poincaré -Wirtinger inequality may only be achieved when g γ is constant on S 2 . Actually by Theorem 5.4.1 of [START_REF] Groemer | Geometric applications of Fourier series and spherical harmonics[END_REF], the equality holds if and only if g γ -E(g γ (U )) is a spherical harmonic of degree 1. This implies that gγ : R 3 → g γ (x) -E(g γ (U )) x is an homogeneous polynomial of degree 1 with constant gradient. But since, for (e i ) 1≤i≤3 the canonical basis, we have ∇g γ (e i ) = ( √ γ i -E(g γ (U )))e i , this implies that √ γ i = E(g γ (U )) and therefore min(γ 1 , γ 2 , γ 3 ) = max(γ 1 , γ 2 , γ 3 ) with g γ (x) = E(g γ (U )) x .

For the other inequality involving γ TMC , we will use the isoperimetric inequality for 3D convex bodies (see [START_REF] Osserman | The isoperimetric inequality[END_REF] for instance), that states that

SA 3 ≥ 36πV 2 ,
where SA is the surface area and V is the volume. Now, taking a 3D ellipsoid of semi-axes

√ γ 1 , √ γ 2 ,
√ γ 3 , we have that its volume is

V = 4 3 π √ γ 1 √ γ 2 √ γ 3 = 4 3 πγ TGC 3/2 ,
and that its surface area is [START_REF] Gusakova | Intrinsic volumes of ellipsoids[END_REF] SA = 2π(γ

1 + γ 2 + γ 3 ) -2πE γ 2 1 U 2 1 + γ 2 2 U 2 2 + γ 2 3 U 2 3 γ 1 U 2 1 + γ 2 U 2 2 + γ 3 U 2 3 = 4πγ TMC .
Therefore the isoperimetric inequality implies that

4 3 π 3 γ 3 TMC ≥ 36π 4 2 3 2 π 2 γ 3 TGC ,
which exactly means γ TMC ≥ γ TGC . Note also that equality in the isoperimetric inequality may only be achieved when the body is a sphere (see [START_REF] Osserman | The isoperimetric inequality[END_REF] p.1190, result (2.7) for instance), that is when min(γ 1 , γ 2 , γ 3 ) = max(γ 1 , γ 2 , γ 3 ).

Visualization and estimation of the anisotropy

We are here interested in the visualization and in the estimation of the anisotropy from an excursion set of a sample of a 3D Gaussian field, in a way similar to what was done in the 2D framework.

We define now two anisotropy ratios:

R := γ TMC γ SA and R G := γ TGC γ SA .
We have 0 ≤ R G ≤ R ≤ 1 because of Theorem 4.1, and since they are symmetric functions of the γ i 's, we can assume γ 3 = max({γ i }), and see R and R G as functions of α 1 := γ1 γ3 and α 2 := γ2 γ3 . This is illustrated on Figure 5.

The Almond and Heart curves of anisotropy

Again, as in the 2D case, we can see the Lipschitz-Killing curvature densities as points on a curve, that is now in R 3 . More precisely, considering

x(t) = SA X (t) SA X (m) = e -(t-m) 2 /2σ 2 , ỹ(t) = TMC X (t) SA X (m) 2 = γ TMC γ SA π √ π 4 √ 2 t -m σ e -(t-m) 2 /2σ 2 = R π √ π 4 √ 2 t -m σ e -(t-m) 2 /2σ 2 , and z(t) = TGC X (t) SA X (m) 3 = γ TGC γ SA 3/2 π 2 8 (t -m) 2 σ 2 -1 e -(t-m) 2 /2σ 2 = R 3/2 G π 2 8 (t -m) 2 σ 2 -1 e -(t-m) 2 /2σ 2 ,
then this defines a parametric curve in R 3 , that has as a projection on the first two coordinates an almond curve of cartesian equation

C (3d) R = {(x, y) ∈ R 2 ; y 2 + π 3 16 R 2 x 2 log x = 0},
and on the last two coordinates, the projection is a parametric curve of the form

C(3d) R,R G = {(y(s), z(s)) ∈ R 2 ; y(s) = R π √ π 4 √ 2 se -s 2 /2 , z(s) = R 3/2 G π 2 8 (s 2 -1)e -s 2 /2 , s ∈ R}.
These two curves are shown on the first line of Figure 7, where the blue curves correspond to the isotropic case (R = R G = 1) and the black ones correspond on the figures to R 0.91 and R G 0.84. These black curves are inside the domains defined by the isotropic (blue) curves.

Estimating the anisotropy from one excursion set

Again, as in the 2D case, assume we observe an excursion set Ê of a sample of a Gaussian stationary random field X, in a window U . The question is then: how to estimate the anisotropy of X from this set Ê only, that is without having access to the whole field (nor its increments for instance). We ignore also at which level t the field X was thresholded, and also what was its mean m or its variance σ 2 . Now, from the four observed Lipschitz-Killing curvatures densities (that are V(E) = V( Ê)/L 3 (U ) the normalized Volume of Ê, SA(E) = SA( Ê)/L 3 (U ) the normalized Surface Area of Ê, TMC(E) = TMC( Ê)/L 3 (U ) the normalized Total Mean Curvature of Ê and TGC(E) = TGC( Ê)/L 3 (U ) the normalized Total Gaussian Curvature of Ê), we can proceed the following way: -Compute an "effective" threshold t, that is an empirical estimate of t-m σ , using the volume of Ê:

t = Φ -1 (1 -V(E)).
-Then estimate SA X (m) by SA m := SA(E)e t 2 /2 .

-And finally estimate the anisotropy ratios via

R = TMC(E) SA(E) 2 4 √ 2 π √ π 1 t e -t 2 /2 and R G 3/2 = TGC(E) SA(E) 3 8 π 2 1 ( t 2 -1) e -t 2 .
As in the 2D case, we have performed some numerical experiments to illustrate all this. Let us mention that such numerical experiments are much more difficult to conduct than in the 2D case. The memory size and the computation time needed to proceed with a 3D volume are very high. Only "small" volumes were considered, with size 200 3 voxels. Therefore the variances in the estimation of the Lipschitz-Killing curvatures is quite high, see Figure 6. The estimation of the anisotropy ratios from one excursion set is illustrated on Figure 7 (anisotropic case). See the caption of the figures for more details. Additional experiments are also provided in the experimental section B of the Appendix.

Further remarks

Link with the kinematic formula

For sake of simplicity, let us assume that X = (X(x)) x∈R d is a stationary smooth centered Gaussian random field with unit variance, whose excursion set is given by E X (t) for some t ∈ R. By the kinematic formula of Theorem 9.4.1 of [START_REF] Schneider | Stochastic and integral geometry[END_REF], under the assumptions that Z := E X (t) is a standard set (see Definition 9.2.1), it follows (see end of p.416) that for r > 0,

E V j Z ∩ B(0, r) = d k=j c k,d-k+j j,d V k (Z)V d-k+j (B(0, r)),
where, following the notations of [START_REF] Schneider | Stochastic and integral geometry[END_REF]

, c k,d-k+j j,d = γ( k+1 2 )γ( d-k+j+1 2 ) γ( j+1 2 )γ( d+1 2 )
, with γ(a) = +∞ 0 t a-1 e -t dt for a > 0, and the intrinsic volumes V l correspond to Lipschitz-Killing curvatures C l . Moreover,

(-γ 1 x 2 1 -γ 2 x 2 2 -γ 3 x 2
3 ) with γ 1 = 0.01, γ 2 = 0.02 and γ 3 = 0.05. An example of such a Gaussian volume is shown on the left of the first line, with one excursion set on the right. Using the Matlab toolbox ImMinkowski [START_REF] Legland | Computation of Minkowski measures on 2D and 3D binary images[END_REF], we have estimated the volume densities (second line, left), the surface area densities (second line, right), the total mean curvature densities (third line, left) and the total Gaussian curvature densities (third line, right). The densities are plotted else as functions of the level t or as functions of the quantile q. The stars are the values for the 20 different samples.

following Corollary 9.4.1, we can identify V k (Z) with C X k (t) and rewrite the previous formula as

E C j Z ∩ B(0, r) = d k=j c k,d-k+j j,d V d-k+j (B(0, 1))C X k (t)r d-k+j .
Meanwhile, under the additional assumption that the joint distribution of the Gaussian vector X(0) is non degenerate, the stationary random field X fulfills the assumptions of Theorem 4.2.1 of [START_REF] Adler | Topological Complexity of Smooth Random Functions: École D' Été[END_REF] and the Gaussian kinematic formula states that (see (4.0.1) therein) with the change of variable k = d -i, where L ∇ k denotes the Lipschitz-Killing curvatures under the Riemannian metric given by Γ, the covariance matrix of ∇X(0),

E L ∇ j Z ∩ B(0, r) = d-j i=0 i + j i L ∇ i+j (B(0, r))(2π) -i/2 h i-1 (t) = d k=j d -k + j d -k L ∇ d-k+j (B(0, 1))(2π) -(d-k)/2 h d-k-1 (t)r d-k+j ,
k k -j = k k -j ω k ω k-j ω j ,
with ω k the volume of the k-dimensional unit ball of R k , that is

ω 0 = 1, ω 1 = 2, ω 2 = π, ω 3 = 4 3 π, etc.. and h -1 (t) = Φ(t), whereas h k (t) = 1 √ 2π H k (t)e -t 2 /2
are the Hermite functions of order k ∈ N. For j = 0 we obtain the mean Euler characteristic of E X (t) ∩ B(0, r) as

E χ(E X (t) ∩ B(0, r) = d k=0 L ∇ d-k (B(0, 1))(2π) -(d-k)/2 h d-k-1 (t)r d-k = d k=0 c k,d-k 0,d V d-k (B(0, 1))C X k (t)r d-k .
Hence, identifying, we get for all 0

≤ k ≤ d, C X k (t) = L ∇ d-k (B(0, 1)) V d-k (B(0, 1)) c d,k h d-k-1 (t),
where we set

c d,k = (2π) -(d-k)/2 c k,d-k 0,d . Note that c d,d = 1, c d,0 = (2π) -d/2 , and c 2,1 = (2π) -1/2 γ(1/2) 2 2 = 1 2 π 2 , while c 3,1 = (2π) -1 1 2 -1 = 1 π and c 3,2 = (2π) -1/2 1 2 -1 = 2 π .
Hence we can define geometrical second spectral moments for 0 ≤ k ≤ d -1 by

λ k = L ∇ d-k (B(0, 1)) V d-k (B(0, 1)) 2/(d-k) .
Note that assuming that Γ = diag(γ 1 , . . . , γ d ) we obtain that

λ k = V d-k (E (Γ)) V d-k (B(0,1)) 2/(d-k)
, where E (Γ) is the Ellipsoid of semi-axes √ γ 1 , . . . , √ γ d (see [START_REF] Gusakova | Intrinsic volumes of ellipsoids[END_REF] for explicit computations). Then, extended isoperimetric inequalities (see (1.1) of [START_REF] Paouris | Random ball-polyhedra and inequalities for intrinsic volumes[END_REF]) rewrite as

∀0 ≤ k ≤ d -1, λ 0 ≤ λ k ,
while the generalized Urysohn inequality (see (1.2) of [START_REF] Paouris | Random ball-polyhedra and inequalities for intrinsic volumes[END_REF]) states that

∀0 ≤ k ≤ d -1, λ k ≤ λ d-1 ,
in accordance with our results for dimension d ∈ {2, 3} with

λ 0 = γ TC , λ 1 = γ Per in dimension d = 2 and λ 0 = γ TGC , λ 1 = γ TMC , λ 2 = γ SA in dimension d = 3 .
Finally, let us end by remarking that the mean Euler characteristic of E X (t) ∩ [0, T ] d for some T > 0 is given according to Theorem 4.4.1 of [START_REF] Adler | Topological Complexity of Smooth Random Functions: École D' Été[END_REF] by

E χ(E X (t) ∩ [0, T ] d ) = d k=0 L ∇ d-k+j ([0, T ] d )(2π) -(d-k)/2 h d-k-1 (t).
and therefore, the Euler characteristic density of

E X (t) ∩ [0, T ] d satisfies lim T →+∞ 1 T d E(χ(E X (t) ∩ [0, T ] d )) = (γ 1 . . . γ d ) 1/2 (2π) -d/2 h d-1 (t) = C X 0 (t).

A functional Weyl's tube formula

Considering a compact d-dimensional submanifold in R d with a C 2 smooth boundary, according to Weyl's tube formula [START_REF] Weyl | On the volume of tubes[END_REF] also called Steiner's formula [START_REF] Federer | Curvature measures[END_REF], the Lipschitz-Killing curvatures are related to the volume of the r-parallel set (also called r-extension or r-dilation) of M by

L d (M r ) = d k=0 ω k C d-k (M )r k ,
where M r = M ⊕ B(0, r) = {x; dist(x, M ) ≤ r} is the r-parallel set to M , and where ω k is still the volume of the k-dimensional unit ball of R k . The above tube formula is valid for r small enough (more precisely smaller than the so-called reach of M , denoted by reach(M ), and that is strictly positive when M is a smooth compact submanifold). Now, this formula can be localized in the sense that if we consider an open bounded subset U ⊂ R d , the above tube formula becomes (see [START_REF] Federer | Curvature measures[END_REF] or [START_REF] Thäle | 50 years sets with positive reach-a survey[END_REF] Theorem 26), for 0 < r < reach(M ),

L d (M r,U ) = d k=0 ω k C d-k (M, U )r k , (23)
where here M r,U = {x; dist(x, M ) ≤ r and Π M (x) ∈ U }, with Π M (x) being the projection of x on M , that is the closest point to x in M which is unique when dist(x, M ) ≤ r < reach(M ). Notice that generally M r,U is not equal to M r ∩ U . Here the Lipschitz-Killing curvatures C d-k (M, U ) are localized in U and their formula involving the principal curvatures are analogous to the ones of the C d-k (M ) given in Section 2 except that the integrals on ∂M are replaced by integrals on ∂M ∩ U . These are the quantities that we have considered in the previous sections. Let f : R d → R be a smooth function, and let U be a bounded open subset of R d . For almost every t ∈ R, E f (t) = {x ∈ R d ; f (x) ≥ t} is a smooth manifold and we can consider its Lipschitz-Killing curvatures in U . We denote them by C k (E f (t), U ) for k = 0, . . . , d. As recalled by Weyl's tube formula [START_REF] Iribarren | Asymptotic behaviour of the integral of a function on the level set of a mixing random field[END_REF], we can obtain these Lipschitz-Killing curvatures as the coefficients of the powers of r in the volume of the r-parallel set to E f (t). More precisely we have, for r small enough,

L d ((E f (t)) r,U ) = d k=0 ω k C d-k (E f (t), U )r k .
In Section 2.2, we have 1) used the formulas for the C k involving the symmetric functions of the curvatures, and the explicit formulas for these curvatures in 2D or 3D; 2) integrate these with a test function h to have formulas for t h(t)C k (E f (t), U ) dt; 3) and finally these integrals were written as an integral on x ∈ U thanks to the coarea formula [START_REF] Evans | Measure theory and fine properties of functions[END_REF].

But we could have done it the other way round. Indeed, we can use directly the coarea formula on t h(t)L d ((E f (t)) r,U ) dt, without knowing the formulas for the curvatures, and in fact recovering them. More precisely, let h be a test function, and for z = (x, s)

∈ U × R + ⊂ R d+1 , let us define the functions h : U × [0, r] → R and F : U × [0, r] → R d+1 by h(z) = h(x, s) = h(f (x)) and F (x, s) = (f (x), x -s ν f (x)) ,
where we denote

ν f (x) = ∇f (x) ∇f (x) ∈ S d-1 if ∇f (x) = 0, and ν f (x) = 0 if ∇f (x) = 0.
A simple computation shows that ν f is differentiable at all x such that ∇f (x) = 0 and that

Dν f (x) = 1 ∇f (x) D 2 f (x) - 1 ∇f (x) 3 ∇f (x)(D 2 f (x)∇f (x)) t .
Now, let K be a compact ball containing U 1 (1-dilation of U ). Then for ε > 0, let us consider the following set of values

V ε = {t ∈ R ; reach(E f (t) ∩ K) ≥ ε}.
Then, since by Morse-Sard's theorem the set of critical values has Lebesgue measure 0 in R, we have that

L 1 (R \ ∪ ε>0 V ε ) = 0. Then, let us denote G ε f = f -1 (V ε ). For r < ε and t ∈ V ε , we have that E f (t) r,U \ (E f (t) ∩ U ) = {x -s ν f (x) ; x ∈ U and f (x) = t and s ∈ (0, r)}.
Therefore, by the coarea formula:

(U ∩G ε f )×[0,r] h(f (x)) det(DF (x, s)) dx ds = Vε h(t)(L d ((E f (t)) r,U ) -L d (E f (t) ∩ U )) dt.
The right-hand term is equal to

d k=1 ω k r k Vε h(t)C d-k (E f (t), U ) dt, whereas in the left-hand term we have det(DF (x, s)) = det ∇f (x) t 0 I d -sDν f (x) -ν f (x) = ∇f (x) det ν f (x) t 0 I d -s Dν f (x) -ν f (x)
.

(24) This determinant is a polynomial in s, with degree d -1 and with constant term (corresponding to s = 0) equal to ±1.

To identify the values of R h(t)C d-k (E f (t), U ) dt for each k we just need to develop the above determinant, and then let ε goes to 0.

When the dimension is d = 2, we have that

det(DF (x, s)) = ∇f (x) 1 -s (Tr(Dν f (x)) -ν f (x) t Dν f (x)ν f (x)) ,
which is positive for s small enough. Therefore, computing the integral on s between 0 and r, we get a polynomial in r, that allows us to identify its coefficients, and we get

ω 1 h(t)C 1 (E f (t), U ) dt = U h(f (x)) ∇f (x) dx,
and

ω 2 h(t)C 0 (E f (t), U ) dt = - 1 2 U h(f (x)) ∇f (x) Tr(Dν f (x)) -ν f (x) t Dν f (x)ν f (x) dx = - 1 2 U h(f (x)) Tr(D 2 f (x)) - ∇f (x) t D 2 f (x)∇f (x) ∇f (x) 2 dx.
We have thus recovered the formula for the curvature in dimension d = 2.

When the dimension is d = 3, considering a direct orthonormal basis of R 3 of the form (ν f (x), T 1 (x), T 2 (x)), developing the above determinant in this basis, we get

det(DF (x, s)) = ∇f (x) (1 -s Tr(Dν f (x)) + s 2 det(ν f (x), Dν f (x) t T 1 (x), Dν f (x) t T 2 (x))),
which is positive for s small enough. Therefore, computing the integral on s between 0 and r, and identifying the coefficients, we get

ω 1 h(t)C 2 (E f (t), U ) dt = U h(f (x)) ∇f (x) dx, ω 2 h(t)C 1 (E f (t), U ) dt = - 1 2 U h(f (x)) ∇f (x) Tr(Dν f (x)) dx = - 1 2 U h(f (x)) Tr(D 2 f (x)) - ∇f (x) t D 2 f (x)∇f (x) ∇f (x) 2 dx,
and

ω 3 h(t)C 0 (E f (t), U ) dt = 1 3 U h(f (x)) ∇f (x) det(ν f (x), Dν f (x) t T 1 (x), Dν f (x) t T 2 (x))) dx = 1 3 U h(f (x)) ∇f (x) t adj(D 2 f (x))∇f (x) ∇f (x) 3 dx,
where adj(D 2 f (x)) is the adjugate matrix of D 2 f (x) (that is the matrix of its cofactors). We have thus recovered the formulas for the mean curvature and the Gaussian curvature for implicit surfaces in R 3 .

Note that when x ∈ R s → L(h, x) is a Lipschitz function, assumption [START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF] implies that Y = (Y k ) k∈Z d is a stationary (BL, θ) dependent field. Then, by Theorem 1.12 p. 178 of [START_REF] Bulinski | Limit theorems for associated random fields and related systems[END_REF], the sequence (S n (h, x)) n will have a Gaussian limit with variance σ 2 as soon as (Var(S n (h, x))) n tends to σ 2 . However, for our specific application, we can not assume that the functional L is Lipschitz but we can assume that L may be approximated by a Lipschitz function according to the following hypothesis:

There exist α > 0 and C ≥ 1 such that ∀ε ∈ (0, 1], there is T ε : R s → [0, 1] such that

F h,ε : x ∈ R s → L(h, x)T ε (x) ∈ R is Lipschitz with Lip(F h,ε ) ≤ C h q ε -1 , E L(h, X(x)) 2 ≤ C h 2 q,0 and E L(h, X(x)) 2 (1 -T ε (X(x))) ≤ C h 2 q,0 ε 2α . (LA)
Proposition 3. Assuming that (25) and (LA) hold, if

i∈Z d sup y ∞≤1 |ρ(y + i)| < +∞,
then the random field Y is (BL, ψ, θ) dependent with

ψ(|I|, |J|, f, g) = Lip(f ) g ∞ |I| + Lip(g) f ∞ |J| + min(|I|, |J|)Lip(f )Lip(g) (27
)

and for C p (h) = 5C 3 p 2 h q (1 + h q ), θ r = C p (h) max      i∈Z d ;|i|≥r sup y ∞ ≤1 |ρ(y + i)|   α α+2 , i∈Z d ;|i|≥r sup y ∞≤1 |ρ(y + i)|    . (28)
Proof. Let r ∈ N and denote θr = i∈Z d ;|i|≥r sup y ∞≤1 |ρ(y + i)|. We set I, J ⊂ Z d with dist(I, J) ≥ r and |J| ≤ |I| such that for any y ∈ [-1, 1] d i∈I ρ(y + ij) ≤ θ(r), for all j ∈ J.

For f ∈ BL(p|I|) and ε ∈ (0, 1] to be fixed later, we write f • F h,ε with a slight abuse of notation to denote the function that is in BL(sp|I|) with Lip(f • F h,ε ) ≤ Lip(f )Lip(F h,ε ). Hence for g ∈ BL(p|J|), we get by ( 25)

|Cov(f • F h,ε (X x,I ), g • F h,ε (X x,J )| ≤ CLip(f )Lip(g)Lip(F h,ε ) 2 i∈I j∈J p l=1 p l =1 ρ(x l + i -(x l + j)) ≤ CLip(f )Lip(g)Lip(F h,ε ) 2 |J|p 2 θr . (29)
Now, let us remark that

|Cov(f (Y I ) -f • F h,ε (X x,I ), g(Y J ))| ≤ 2 g ∞ E(|f (Y I ) -f • F h,ε (X x,I )|) ≤ 2 g ∞ Lip(f )E ( L(h, X x,I ) -F h,ε (X x,I )) ∞ ) ≤ 2 g ∞ Lip(f ) k∈I p l=1 E (|L(h, X x l +k )|(1 -T ε (X x l +k ))) ≤ 2C 1/2 p|I| g ∞ Lip(f ) h q,0 ε α ,
by assumption (LA), using Cauchy-Schwarz inequality. Similarly, one has

|Cov(f • F h,ε (X x,I ), g(Y J ) -g • F h,ε (X x,J ))| ≤ 2C 1/2 p|J| f ∞ Lip(g) h q,0 ε α .
Since C ≥ 1, by definition of [START_REF] Legland | Computation of Minkowski measures on 2D and 3D binary images[END_REF], it follows that

|Cov (f (Y I ), g(Y J ))| ≤ C p (h) 5 ψ(|I|, |J|, f, g) θr ε -2 + 4ε α .
Now, when θr ∈ (0, 1] we can choose ε = ( θr )

1 α+2 to get |Cov (f (Y I ), g(Y J ))| ≤ C p (h)ψ(|I|, |J|, f, g) θ α α+2 r .
Otherwise, if θr = 0 the result comes from letting ε tend to 0 and if θr ≥ 1 by choosing ε = 1.

However, we can in fact be more precise in term of covariance control.

Proposition 4. Assuming that (25) and (LA) hold, one has for all

x, x ∈ R d |Cov (L(h, X(x)), L(h, X(x )))| ≤ 3C 3 h 2 q max ρ(x -x ) α α+2 , ρ(x -x ) .
Proof. As previously, let us choose ε ∈ (0, 1] to be fixed later and write for

x ∈ R d L(h, X(x)) = F h,ε (X(x)) + L(h, X(x)) (1 -T ε (X(x))) .
Noting that x may be written as x = x + k for some x ∈ [0, 1] d and k ∈ Z d , by [START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF] we simply get for

x ∈ R d , |Cov(F h,ε (X(x)), F h,ε (X(x ))| ≤ CLip(F h,ε ) 2 ρ(x -x ) ≤ C 3 h 2 q ε -2 ρ(x -x ).
Now, using Cauchy-Schwarz inequality, we obtain

|Cov(L(h, X(x)) (1 -T ε (X(x))) , L(h, X(x ))| ≤ Var(L(h, X(x)) (1 -T ε (X(x))) 1/2 Var(L(h, X(x )) 1/2
≤ C h 2 q,0 ε α , by (LA) using that the variance is smaller than the second order moment and the fact that T ε has values in [0, 1]. Moreover, we obtain similarly

|Cov(F h,ε (X(x)), L(h, X(x )) (1 -T ε (X(x ))))| ≤ Var(F h,ε (X(x)) 1/2 Var(L(h, X(x )) (1 -T ε (X(x ))) 1/2 ≤ C h 2 q,0 ε α , since E F h,ε (X(x)) 2 ≤ E L(h, X(x)) 2 . Hence, |Cov (L(h, X(x)), L(h, X(x )))| ≤ C 3 h 2 q ε -2 ρ(x -x ) + 2ε α .
Note that if ρ(x -x ) = 0 we get the upper-bound by letting ε tends to 0

. If ρ(x -x ) ∈ (0, 1] we choose ε = ρ(x -x ) 1
α+2 and otherwise we choose ε = 1 to get the result.

Then, under a natural decay assumption, we can define the asymptotic covariances of the sequence (S n (h, x)) n . In particular this ensures the so-called finite susceptibility assumption (see (5.3) p.90).

Corollary 3. Assuming that (25) and (LA) hold, we also assume that there exists β >

d(α+2) α such that, for all x ∈ R d , ρ(x) ≤ C(1 + |x|) -β . (30)
Then, there exists C α > 0 such that, for all Proof. Let x ∈ R d with x ≤ 1. First note that we may assume that ω X (x) ∈ (0, C]. Let us remark that, as previously, using the fact that for ε ∈ (0, 1]

x ∈ R d with x ≤ 1, E |L(h, X(x)) -L(h, X(0))| 2 1/2 ≤ C α h q ω X (x) α 1+α . (33
L(h, X(x)) = F h,ε (X(x)) + L(h, X(x))(1 -T ε (X(x))),
we have

E |L(h, X(x)) -L(h, X(0))| 2 1/2 ≤ Lip(F h,ε )E [ X(x) -X(0) ∞ ] + 2C 1/2 h q,0 ε α ≤ C h q ε -1 ω X (x) + 2ε α . Taking ε = ω X (x) C 1 1+α
we obtain [START_REF] Spivak | A comprehensive introduction to differential geometry[END_REF].

Assuming moreover [START_REF] Rataj | Curvature measures of singular sets[END_REF] for

β > d(α+2) α , for K ≥ 1, E |B(h, x) -B(h, 0)| 2 = 2 k∈Z d Cov(L(h, X(k)), L(h, X(0)) -L(h, X(x)), ≤ 2C 1/2 h q,0 E |L(h, X(x)) -L(h, X(0))| 2 1/2 (2K + 3) d + 6C 3 h 2 q |k|>K+1 [max ρ(k + x) α α+2 , ρ(k + x) + max ρ(k) α α+2 , ρ(k) ],
using Cauchy-Schwarz inequality and Proposition 4. Hence, by [START_REF] Rataj | Curvature measures of singular sets[END_REF] and [START_REF] Spivak | A comprehensive introduction to differential geometry[END_REF], we can find Cα,β,d > 0 such that

E |B(h, x) -B(h, 0)| 2 ≤ Cα,β,d h 2 q K d ω X (x) α 1+α + K -( βα α+2 -d) . Then, choosing K = ω X (x) -α+2 αβ α 1+α , there exists C α,β,d > 0 such that E |B(h, x) -B(h, 0)| 2 ≤ C α,β,d h 2 q ω X (x) α 1+α (1-d(α+2) αβ
) .

In view of [START_REF] Shashkin | Quasi-associatedness of a Gaussian system of random vectors[END_REF] and by stationarity of (B(h, x))

n , h = L d (U n ) 1 L d (U n ) Un L(h, X(x))dx -E(L(h, X(0))) d -→ n→+∞ [0,1] d B(h, x)dx,
where (B(h, x)) x∈[0,1] d is the continuous Gaussian random field with covariance function given by [START_REF] Schneider | Stochastic and integral geometry[END_REF]. The random variable [0,1] d B(h, x)dx is centered Gaussian with variance given by

R d Cov (L(h, X(z)), L(h, X(0))) dz. (34)
Hence, there exists B a generalized Gaussian random process on (S , B(S )), with characteristic functional given by

L B (h) := E e i B,h = E e i [0,1] d B(h,x)dx = exp - 1 2 R d Cov (L(h, X(z)), L(h, X(0))) dz such that T n distrib.
-→ n→+∞ B, where the convergence holds in distribution with respect to the strong topology of S (R).

Proof. In view of Proposition 5, by continuity of (B(h, x)) x∈[0,1] d we have the convergence of the Riemann sums

[0,1] d B(h, x)dx = lim m→+∞ 1 m d l∈(Z∩[0,m)) d B h, l m .
Hence the random variable [0,1] d B(h, x)dx is centered Gaussian with variance given by

[0,1] d [0,1] d Cov (B(h, x), B(h, y)) dxdy = [0,1] d [0,1] d k∈Z d Cov (L(h, X(x + k)), L(h, X(y))) dxdy = R d Cov (L(h, X(z)), L(h, X(0))) dz,
by change of variables and Fubini's theorem. Note also that for

x ∈ [0, 1] d , m > 0 and l ∈ (Z ∩ [0, m)) d we have S n h, x + l m -S n h, l m = 1 n d/2 k ∈nB∩Z d L h, X x + l m + k -L h, X l m + k .
Then, by stationarity,

Var S n h, x + l m -S n h, l m ≤ k ∞≤n |r(h, k, x)|, where r(h, k, x) = Cov (L(h, X(x)) -L(h, X(0)), L(h, X(x + k)) -L(h, X(k)))
satisfies, by Proposition 5 and Proposition 4 together with [START_REF] Rataj | Curvature measures of singular sets[END_REF],

|r(h, k, x)| ≤ C 2 α h 2 q ω X (x) 2α 1+α and |r(h, k, x)| ≤ 12C 3 h 2 q |k| -βα α+2 .
It follows that we can find C α,β,d > 0 such that

E S n h, x + l m -S n h, l m 2 = Var S n h, x + l m -S n h, l m ≤ C α,β,d h 2 q ω X (x) 2α 1+α (1-d(α+2) αβ ) .
Then,

[0,1] d S n (h, x)dx - 1 m d l∈(Z∩[0,m]) d S n h, l m L 2 (Ω) ≤ l∈(Z∩[0,m]) d [0,1/m] d S n (h, x + l m ) -S n h, l m dx L 2 (Ω)
Using Cauchy-Schwarz inequality,

E [0,1/m] d S n (h, x + l m ) -S n h, l m dx 2 ≤ 1 m d [0,1/m] d E S n (h, x + l m ) -S n h, l m 2 dx.
Then,

[0,1] d S n (h, x)dx - 1 m d l∈(Z∩[0,m]) d S n h, l m L 2 (Ω) ≤ C 1/2 α,β,d h q m d/2 [0,1/m] d ω X (x) 2α 1+α (1-d(α+2) αβ ) dx 1/2
.

In view of [START_REF] Shashkin | Quasi-associatedness of a Gaussian system of random vectors[END_REF], this implies that

[0,1] d S n (h, x)dx = lim m→+∞ 1 m d l∈Z∩[0,m] d S n h, l m ,
where the convergence holds in L 2 (Ω), uniformly in n. Since Cov (L(h, X(x + k)), L(h, X(y))) .

L d (U n ) 1 L d (U n ) Un L(h, X(x))dx -E(L(h, X(0))) = [0,1] d S n (h, x)dx,
It follows by Proposition 4 that, on the one hand we can use Fubini's theorem to obtain [START_REF] Thäle | 50 years sets with positive reach-a survey[END_REF], and on the other hand σ 2 (h) ≤ C h 2 q and therefore the characteristic functional L B : S(R) → C given by L B (h) = e -1 2 σ 2 (h) , is continuous at 0 and the conclusion comes from Theorem 2.3 of [START_REF] Biermé | Generalized random fields and Lévy's continuity theorem on the space of tempered distributions[END_REF].

A.2. Gaussian case and Proof of Theorem 2.1

Assuming that X is a Gaussian field satisfying the assumptions of Theorem 2.1, by the Cramer-Wold device, it is enough to prove that for (a j ) 0≤j≤d ∈ R d+1 one has

L d (U n ) d j=0 a j C j (E X (•), U n ) L d (U n ) -C X j , distrib.
-→ n→+∞ d j=0 a j B j , where (B j ) 0≤j≤d are the centered Gaussian generalized processes introduced in Theorem 2.1 and the convergence holds in distribution with respect to the strong topology of S (R). Therefore we let

T n = L d (U n ) d j=0 a j C j (E X (•), U n ) L d (U n ) -C X j
and B = d j=0 a j B j .

Hence we need to check that the assumptions of Theorem A.2 are satisfied for h ∈ S and L(h, X(x)) := d j=0 a j L j (h, X(x)), since

T n , h = L d (U n ) 1 L d (U n ) Un
L(h, X(x))dx -E (L(h, X(0))) .

Then we will identify the Gaussian limit. First, since X is a Gaussian field the quasi-association assumption ( 25) is satisfied (see [START_REF] Shashkin | Quasi-associatedness of a Gaussian system of random vectors[END_REF]). Second, for x = (u, v, w) ∈ R × R d × R d(d+1)/2 with a slight abuse of notation we write 2 with C that does not depend on M, δ, q = 2, and h q = h q,0 + h q,1 . Now P( ∇X(x) ≤ δ) ≤ Cδ Note also that for any p > 1, choosing p < d for d = 3, one has E(L(h, X(x)) 2p ) ≤ C p h 2p q,0 . Hence, by Hölder inequality, inequalities (LA) are satisfied for α ∈ (0, d 2 (1 -1/p)). Then, under Assumption [START_REF] Cheng | The mean Euler characteristic and excursion probability of Gaussian random fields with stationary increments[END_REF] with β-d > 12, when d = 2, since (A 1 ) is satisfied for all δ > 0, and 2d(1+4/d) = 12, one can find p > 1 large enough such that 2d(1 + 2 (1 -1/p) we can check both (A 1 ), (A 2 ) and [START_REF] Rataj | Curvature measures of singular sets[END_REF]. Finally, note that for x ≤ 1 one has

E [X ij (x) -X ij (0)] 2 = 2(-1) ε(i)+ε(j) ∂ 2ε(i)+2ε(j) ijij ρ(0) -∂ 2ε(i)+2ε(j) ijij ρ(x) ,
such that (32) is satisfied by (H).

It follows that all assumptions and therefore conclusions of Theorem 2.1 hold true.

To conclude it only remains to identify the asymptotic covariance. By linearity of B and L for Estimating the anisotropy from one level. We consider here the same two Gaussian random fields as in Figure 3. First line: estimation of γ Per from one level set of one sample (of size 1000 × 1000), as a function of the quantile, in the isotropic case (left figure) and in the anisotropic one (right figure). Second line: same experiment but for the estimation of γ TC . Third line : same experiment but for the estimation of the anisotropy ratio R, that is 1 on the left and 0.92 on the right. The confidence intervals have been obtained with 100 samples. Notice how the estimations of γ TC and of R degenerate around the quantile q = 0.5 since it corresponds to the effective level t = 0, and the division by t is therefore unstable.

  Fig 1.Two excursion sets of Gaussian random fields. The left one is isotropic and the right one is not. They have the same area, the same perimeter but different total curvature. Their geometry can be summarized as a point that is on the almond curve in the isotropic case and that is inside the almond domain in the anisotropic case.

Fig 2 .

 2 Fig 2. On the left, plot of the function α → R. On the right, the almond curves of anisotropy: the blue plain curve is the curve C 1 corresponding to isotropy, and the black dotted curve is a curve C R , with here R 0.92.

Fig 3 . 1 + γ 2 x 2 2 )

 312 Fig 3.First line: left, a sample of a centered isotropic Gaussian field X, size 1000 × 1000 pixels, with variance σ 2 = 1 and covariance of the form ρ(x 1 , x 2 ) = exp(-γ 1 x 2 1 + γ 2 x 2 2) with γ 1 = γ 2 = 0.005. Middle: the almond curve and some points from the isotropic field (100 samples). Right: corresponding estimations of the anisotropy ratio R (= 1 here) as a function of the quantiles. Second line: same experiments but with now γ 1 = 0.005 and γ 2 = 0.002, that corresponds to an anisotropic field with R 0.92.

Fig 4 .

 4 Fig 4. Densities C X 0 (first line), C X 1 (second line) and C X 2 (third line), as a function of the level or the quantile in the isotropic case (left part of the figure) and in the anisotropic one (right part of the figure). The Gaussian fields are the same as the ones of Figure 3, here with one sample of size 1000 × 1000 (n = 1000).

Fig 5 . 3 and α 2 := γ 2 γ 3 .

 5323 Fig 5. The anisotropy ratios R (on the left) and R G (on the right) as functions of α 1 := γ 1 γ 3 and α 2 := γ 2 γ 3 .

Fig 6 .

 6 Fig 6. Lipschitz-Killing curvatures densities of the excursion sets of an isotropic Gaussian volume. Here we have taken 20 samples of a 3D volume of size 200 3 voxels, of a Gaussian field with mean m = 0, variance σ 2 = 1 and covariance of the form exp(-γ 1 x 2 1 -γ 2 x 2 2 -γ 3 x23 ) with γ 1 = 0.01, γ 2 = 0.02 and γ 3 = 0.05. An example of such a Gaussian volume is shown on the left of the first line, with one excursion set on the right. Using the Matlab toolbox ImMinkowski[START_REF] Legland | Computation of Minkowski measures on 2D and 3D binary images[END_REF], we have estimated the volume densities (second line, left), the surface area densities (second line, right), the total mean curvature densities (third line, left) and the total Gaussian curvature densities (third line, right). The densities are plotted else as functions of the level t or as functions of the quantile q. The stars are the values for the 20 different samples.

Fig 7 .

 7 Fig 7. Estimation of R and R G from the Lipschitz-Killing curvatures densities of an excursion set. The samples are the same as the ones of Figure 6. Notice that, because of the term 1/ t, the estimation of R is inaccurate for values around t = 0, whereas for R G the inaccuracy of estimation occurs around t = ±1.

) 2 α 1+α 1 -

 21 Assuming moreover[START_REF] Rataj | Curvature measures of singular sets[END_REF] for β > d(α+2) α , one can choose a version of (B(h, x)) x∈R d , given in Corollary 3, that is γ-Hölder continuous on any compact sets of R d for all γ < δ d(α+2) αβ .

  the result will simply follows from the fact that for m > 0, by Corollary 4Note that the variance of the Gaussian variable [0,1] d B(h, x)dx is explicitly given byσ 2 (h) := [0,1] 2 [0,1] 2Cov(B(h, x), B(h, y))dxdy, with Cov(B(h, x), B(h, y)) = k∈Z d

F 0

 0 (x) = v , F 1 (x) = -Tr(w) + v t wv v 2 and F 2 (x) = -v t adj(w)v v 3 ,and remark that x → L d (h, x) is bounded Lipschitz and clearly satisfies (LA). However for 0≤ j ≤ d -1 we need to conveniently approximate x → L j (h, x) = h(u)F d-1-j (x). For δ > 0 we can find a C 1 function T δ : R d → [0, 1] such that T δ (v) = 1 for v > δ and ∇T δ ≤ δ -1 . Moreover for M > 0 we can find a C 1 function T M : R d × R d(d+1)/2 → [0, 1] such that T M (v, w) = 1 for (v, w) < M and T M (v, w) = 0 for (v, w) > 2M with ∇T M ∞ ≤ 1. Then x → L(h, x)T M (v, w)T δ (v) is C 1 with gradient bounded by C h q M δ

  d and one can find c d > 1 such that for M large enough P( (∇X(x), D 2 X(x)) ≥ M ) ≤ e -c d M . For ε ∈ (0, 1], we set δ = d 2c d ε 1/2 | log(ε/2)| and M = d 2c d | log(ε/2)| such that δ M = ε and for all p ≥ 1, we have E((1 -T M (∇X(x), D 2 X(x))T δ (∇X(x))) p ) ≤ P( ∇X(x) ≤ δ) + P( (∇X(x), D 2 X(x)) ≥ M ) ≤ Cε d/2 (1 + | log(ε/2)| d ).

4 d( 1 1 1 4 d 1 1- 1

 411411 -1/p) ) < β -d and take p > p with α ∈ (0, d 2 (1-1/p )) and then α α+2 (β-d) > 2d so that (A 2 ) and together with (30) are also satisfied. When d = 3, since 2d 4 d -1/d = 12, one can find δ < 1 and p < d such that d(1+ 1 δ ) /p < β -d. Then taking α = d

Fig 8 .

 8 Fig 8. Checking the central limit behavior. First line: estimation of R in the isotropic case (left) and in the anisotropic one (right, here R 0.92). The estimation is done for three different quantiles (corresponding to the three colors) and on different domain size T (x-axis of the plot). The confidence intervals have been obtained with 100 samples. Second line: standard deviation of the estimated value of R in the isotropic case (left) and in the anisotropic one (right), for three different quantiles as a function of the inverse domain size 1/T . This plot shows experimentally that the empirical standard deviation behaves like 1/T , with a slope that depends on the quantile.

Fig 9 .

 9 Fig 9.Estimating the anisotropy from one level. We consider here the same two Gaussian random fields as in Figure3. First line: estimation of γ Per from one level set of one sample (of size 1000 × 1000), as a function of the quantile, in the isotropic case (left figure) and in the anisotropic one (right figure). Second line: same experiment but for the estimation of γ TC . Third line : same experiment but for the estimation of the anisotropy ratio R, that is 1 on the left and 0.92 on the right. The confidence intervals have been obtained with 100 samples. Notice how the estimations of γ TC and of R degenerate around the quantile q = 0.5 since it corresponds to the effective level t = 0, and the division by t is therefore unstable.

Fig 10 . 1 -γ 2 x 2 2 -γ 3 x 2 3 )

 10123 Fig 10. Lipschitz-Killing curvature densities of the excursion sets of an isotropic Gaussian volume. Here we have taken 20 samples of a 3D volume of size 200 3 voxels, of a Gaussian field with mean m = 0, variance σ 2 = 1 and covariance of the form exp(-γ 1 x 2 1 -γ 2 x 2 2 -γ 3 x 2 3) with here γ 1 = γ 2 = γ 3 = 0.01. An example of such a Gaussian volume is shown on the left of the first line, with one excursion set on the right. Using the Matlab toolbox ImMinkowski[START_REF] Legland | Computation of Minkowski measures on 2D and 3D binary images[END_REF], we have estimated the volume densities (second line, left), the surface area densities (second line, right), the total mean curvature densities (third line, left) and the total Gaussian curvature densities (third line, right). The densities are plot else as functions of the level t or as functions of the quantile q. The stars are the values for the 20 different samples.

Fig 11 .

 11 Fig 11. First line: the samples of the previous figure are plotted as points on the almond and heart curve of isotropy. Second line: estimation of R and R G from the Lipschitz-Killing curvature densities of the excursion sets of the samples.

  Then, using the definition of the F k in (2), we have for a.e. t ∈ R

∂E f (t)∩U

  measurable function, as well as its first and second order partial derivatives. By stationarity, we deduce from (10) that U |F k (X(x))|1 I ∇X(x)) >0 dx < +∞ a.s. and (8) is satisfied for k ∈ {0, 1, d -1}. Hence, there exists a negligible set N of A such that for j∈ {0, d -2, d -1}, C j (E X(ω) (•), U ) ∈ L 1 (R) ⊂ S for ω /∈ N and we can set C j (E X(ω) (•), U ) = 0 for ω ∈ N . We

	can therefore define LC	X(ω) j	(h, U ), for h ∈ S, satisfying (9) for ω / ∈ N and LC X(ω) j

  x∈R d , we can choose a version such that (B(h, x)) x∈R d is γ-Hölder continuous on any compact set of R d for γ < δ

	2	α 1+α 1 -d(α+2) αβ	(see Proposition 9 in
	[6] for instance).
	Theorem A.2. Under the assumptions of Corollary 4, assuming moreover (32), one has for
	h ∈ S,
	T

available at https://www.mathworks.com/matlabcentral/fileexchange/33690-imminkowski

Acknowledgments

This work has been partly supported by the project ANR MISTIC (ANR-19-CE40-0005). The authors would like to thank Paul Doukhan for pointing them the appropriate reference on λdependent fields.

Appendix A: Central Limit Theorem

A.1. General results

We assume that X is a C 2 stationary second order field (not necessarily Gaussian) and write as before X = (X, ∇X, D 2 X) the vector-valued field with values in R s = R × R d × R d(d+1)/2 , following [START_REF] Estrade | A central limit theorem for the Euler characteristic of a Gaussian excursion set[END_REF]. We recall that since X is C 2 stationary second order, the covariance function x → ρ(x) = Cov(X(x), X(0)) is C 4 and for 0 ≤ i, j, k, l ≤ d Cov(X ij (x), X kl (0)) = (-1) ε(k)+ε(l) ∂ ε(i)+ε(j)+ε(k)+ε(l) ijkl ρ(x), where ε(i) = 1 if i = 0 and X 0j (x) = X j0 (x) = X j (x) while X 00 (x) = X(x). Hence we introduce

In order to prove a general central limit theorem, we will work in this section under a strong quasi-association assumption, namely, we assume that there exists C ≥ 1 such that, for all finite subsets I and J of Z d , for all p ≥ 1 and all x = (x 1 , . . . , x p ) ∈ [0, 1] dp , for all Lipschitz functions f : R s|I|p → R and g : R s|J|p → R, we have

where X x,I = (X(x l + i), i ∈ I, 1 ≤ l ≤ p). For h ∈ S and q ∈ N, we write h q = h q,0 + h q,1 . Note that h is a Lipschitz bounded function with Lip(h) = h 0,1 ≤ h q,1 and that h 0,0 ≤ h q,0 . We consider a sequence (T n ) n of generalized random processes given, for n ∈ N, U n = (0, n) d and h ∈ S, by

The quasi-association assumption [START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF] will imply some dependence properties on the stationary vector-valued random field Y = (Y k ) k∈Z d .

In the following, we will use the notation BL(q) for q ≥ 1 to denote the set of functions from R q to R that are Lipschitz and bounded.

Definition A.1 (Definition 5.14 p.94 of [START_REF] Bulinski | Limit theorems for associated random fields and related systems[END_REF]). The random field Y is called (BL, ψ, θ) dependent if there exists a non-increasing sequence θ = (θ r ) r∈N with lim r→+∞ θ r = 0 such that for any disjoint finite sets I, J ⊂ Z d with dist(I, J) ≥ r, and any bounded Lipschitz functions f ∈ BL(p|I|), g ∈ BL(p|J|), one has

The case where ψ(|I|, |J|, f, g) = Lip(f )Lip(g) min(|I|, |J|) is simply called (BL, θ) dependence.

Then, there exists C α,β,d > 0 such that, for all x, y ∈ R d ,

Hence, there exists a real-valued stationary Gaussian random field (B(h, x)) x∈R d with covariance function given by

Proof. Note that by Proposition 4, since β α α+2 > d, one has for all x, y ∈ R d ,

and one can choose

. Moreover by stationarity of (L(h, X(x)))

Hence one can define the stationary Gaussian field (B(h, x)) x∈R d with covariance function given by ( 31) since this function is of non-negative type as a limit of covariance functions.

Note that when x → L(h, x) is Lipschitz with E(L(h, X(x)) 2 ) < +∞, assumption (LA) is satisfied for all α > 0 since one can choose T ε = 1. Then, as soon as there exists β > d such that (30) holds true, we obtain the asymptotic normality of the R p -valued vector (S n (h, x)) n for any p ≥ 1 and x = (x 1 , . . . , x p ) ∈ ([0, 1] d ) p , by Corollary 1.13 p.180 of [START_REF] Bulinski | Limit theorems for associated random fields and related systems[END_REF], using the (BL, θ) dependence of Y . In the general case, noting that the ψ-dependency coefficient in [START_REF] Legland | Computation of Minkowski measures on 2D and 3D binary images[END_REF] satisfies

the field Y is also λ-dependent as defined in [START_REF] Doukhan | Weak dependence, models and some applications[END_REF]. By [START_REF] Rataj | Curvature measures of singular sets[END_REF] we obtain that λ Y (r)(= θ r ) = O r→+∞ (r -λ ) with λ = α α+2 (β -d) and θ r given by [START_REF] Osserman | The isoperimetric inequality[END_REF]. Then, using Cramer-Wold device and Theorem 2 p.219 in [START_REF] Doukhan | Weak dependence, models and some applications[END_REF] we obtain the asymptotic normality under stronger assumptions (A 1 ) and (A 2 ) as stated in the following corollary.

Corollary 4. Assuming that (25), (LA) and (30) hold, assuming moreover that

This result will allow us to get a central limit theorem for the sequence of generalized processes given by [START_REF] Kratz | Central limit theorems for level functionals of stationary Gaussian processes and fields[END_REF], by approximating integrals by Riemann sums, under an appropriate continuity assumption.

Proposition 5. Assuming that (25) and (LA) hold, if there exists δ > 0 such that for all

h, h ∈ S we obtain that

where we use the fact that, by stationarity and change of variables, one has Cov L j (h, X(z)), L k ( h, X(0)) dz.

Appendix B: Additional numerical experiments

We provide in this section several additional figures to illustrate, through numerical experiments, the statistical estimation of the anisotropy from the Lipschitz-Killing curvatures of an excursion set.

B.1. In the 2D case

On Figures 8 and9 we show error bars on the estimation of the anisotropy from one excursion set of a 2D stationary Gaussian field, and we illustrate the central limit behavior by plotting the standard deviation of the estimation of the anisotropy ratio R as a function of the observation window size. See the caption of the figures for detailed comments on these numerical experiments.

B.2. In the 3D case

On Figures 10 and11, we show some additional 3D experiments in the isotropic case. As in the anisotropic case, the estimation of R = R G = 1 here is inaccurate when the effective level is close to 0 (for R) or close to ±1 for R G . Now, these two ratios could be combined to accurately estimate the isotropy whatever the effective level.