Elena Di

Elena Di Bernardino

Céline Duval

Statistics for Gaussian Random Fields with Unknown Location and Scale using Lipschitz-Killing Curvatures

Keywords: Geometric inference, Excursion sets, Image analysis AMS Classification: 62H11, 62M40, 60G60

HAL is

Introduction

Lipschitz-Killing (LK) curvatures are geometrical tools which permit to analyse d dimensional objects. Considering a black and white image in dimension d = 2, there are three LK curvatures: the surface area, the half perimeter and the Euler characteristic. Each of them brings a distinct information on the geometry of the black (resp. white) zone. The surface area is related to its occupation density, the perimeter to its regularity and the Euler characteristic to its connectivity. the image of one excursion set if the field is symmetric, which is the case of Gaussian fields, the expected mean of the Euler characteristic of this set will be null at level u = 0 and the expected mean of its area, will be 1/2 of the size of this image at this same level. More generally, the mean of the field indicates where the excursion set shows black and white zones in comparable proportions, whereas the variance conveys information on the range of levels where the excursion sets are non degenerate (all white or all black images). If one only observes one excursion set -and not the whole field X-it is often impossible to test whether the underlying field is centered or has unit variance.

It is impossible to estimate the mean and variance of X from the sole excursion set. Then, from a statistical viewpoint, it is unrealistic to impose that the unobserved field X is standard.

The knowledge of the mean and variance of X is also crucial when it comes to compare the LK curvatures of several excursion sets. To compare the geometries of two given excursion sets, this information turns out to be crucial. Indeed, if the two excursion sets correspond to fields with distant mean, observing them at the same level will generate different pictures, but it does not mean that the underlying fields are necessarily different, e.g. one can be the translation of the other. This nonstandard context has deserved only a limited attention in the literature, we can mention, for instance, [START_REF] David | Testing for a Signal with Unknown Location and Scale in a Stationary Gaussian Random Field[END_REF] based on the Hadwiger characteristic.

In this article, in the case where X is a stationary and isotropic Gaussian random field, we introduce two quantities that we call effective level and effective spectral moment of the field (see Definition 1.2) that contain simultaneously the knowledge of the variance and mean of the field for every possible observation levels. Comparing excursion sets at the same effective level is meaningful as illustrated in Figure 1. In the first row we observe three different digital mammograms at the same fixed level. The obtained excursion sets look completely different; from this perspective the mammograms are dissimilar. But on the second row, we change the perspective and look at these objects at the same effective level, compared to the first row, the proportions of black and white zones are visually alike.

We propose consistent and asymptotically normal estimators for these effective level and effective spectral moment when X is a stationary and isotropic Gaussian random field (see Propositions 2.1 and 2.2). We relate these quantities to the model parameters and we show how these estimators permit to recover some model parameters and to build testing procedures.

Comparing images of exclusion sets

To compare two different fields from the observation of one excursion set of each using LK curvature devices, observing them at the same level u is irrelevant if the underlying fields have distant mean and/or variance. It is meaningful to compare them only if they are observed at the same effective level and effective spectral moment. This is why it is important to develop a consistent test to decide whether two images of excursion sets have the same effective level or not. This test is a mandatory preliminary step before putting in comparison two images of excursion sets. Its construction makes a full use of the asymptotically normal estimator of effective level as well as the consistent estimator of its limit variance, introduced in the present work (see Corollary 4.1).

Asymptotic variance estimation Frequently, when calibrating testing procedures or studying the deviations of the estimators, the proper calibration depends on the unknown asymptotic variances. Therefore, having consistent estimators of the asymptotic variances of the empirical LK curvatures is essential for the statistical study of these objects. If estimation strategies are known and used in practice in a wide range of contexts, there are few theoretical results concerning the consistency of these procedures. Some results exist in the case of the area if X is a Gaussian random field and under conditions on the correlation function of X that are difficult to check (see [START_REF] Mattfeldt | Characterization of mammary gland tissue using joint estimators of Minkowski functionals[END_REF], [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF]).

We establish in Theorem 3.1 the consistency of the empirical variance of the area when X is a stationary and isotropic Gaussian field whose correlation function decays polynomially (see Assumption (A1) below). For that we use the decomposition of Itô-Wiener chaos of the LK curvatures jointly with the diagram formula (see [START_REF] Taqqu | Law of the iterated logarithm for sums of non-linear functions of Gaussian variables that exhibit a long range dependence[END_REF]). We also argue why this estimator should remain consistent to estimate the Euler characteristic variance (see Corollary 3.1). We point out that in this article, we are mainly interested in the Euler characteristic and area as they are numerically more stable to compute (see, e.g., the discussion about the well known numerical instability of the perimeter estimation in [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF]). Moreover, as we prove results on the asymptotic variance of the area, the numerical results in Section 4.1 rely only on this latter quantity.

Outline of the paper The paper is organized as follows. In the remaining of this section we define the three objects of interest, i.e., the LK densities for the excursion sets of a two-dimensional Gaussian random field which is non centered and whose variance is not necessarily one. These formulae depend on two identifiable parameters, called effective level and effective spectral moment. In Section 2 we estimate these parameters and establish asymptotic normality results. The associated asymptotic variances are estimated in Section 3, using a subwindow technique. The consistency of the proposed subwindow empirical variance estimators is studied in Section 3.1 as well as some statistical implications in Section 3.2. In Section 4 we build a test based on the effective level to compare two images of excursion sets. This test is put into practice in Section 4.1 on synthesized 2D digital mammograms provided by GE Healthcare France (department Mammography) and in Section 4.2 on real digital mammograms provided by the Mammographic Image Analysis Society (MIAS). The proofs of the technical results are gathered in Appendix A.

Definitions and preliminary notions

In the whole paper, | • | denotes equally the absolute value or the two dimensional Lebesgue measure and by | • | 1 its one-dimensional Hausdorff measure, • denotes the Euclidian norm. We consider X being a Gaussian field defined on R 2 satisfying the following hypothesis.

(A0) The Gaussian field X is stationary, isotropic with mean E[X(0)] = µ, variance V(X(0)) = σ 2 and V(X (0)) = λI 2 for λ > 0, the second spectral moment, σ > 0, µ ∈ R and I 2 the 2 × 2 identity matrix. Moreover, the trajectories of X are almost surely of class C 3 . Denote by r the covariance function of X, i.e., for any fixed t, r(t) = Cov(X(0), X(t)) and by ρ the correlation function ρ(t) = corr(X(0), X(t)). Under Assumption (A0) the field is suppose to be C 3 , then the covariance function is C 6 . In some cases, e.g., to get the asymptotic normality or consistency of our estimators, we impose additional assumptions on covariance and correlation functions of X gathered in the following.

(A1) |ρ(t)| ≤ (1 + t ) -γ , γ > 2.
(A2) For any fixed t in R 2 , the covariance matrix of the random vector (X(t), X (t), X (t)) has full rank and the covariance function r of X is such that,

R 2 r(s) ds > 0, M r (t) → 0 when t → +∞ and M r ∈ L 1 (R 2 ) , where M r (t) = max ∂ k ∂t k r(t) ; k = (i 1 , . . . , i ) ∈ {1, 2} , 0 ≤ ≤ 4 and ∂ k r ∂t k (t) = ∂ ∂x i 1 ...∂x i r(t).
Assumptions (A0), (A1) and (A2) are standard when studying limit laws of non linear functionals of stationary Gaussian random fields. Assumptions (A0) and (A1) are used when working with the area, whereas for the Euler characteristic we rely on (A0) and (A2).

Rectangles, denoted by T , in R 2 are bounded and with non empty interior. In the following notation T R 2 stands for the limit along any sequence of bounded rectangles that grows to R 2 . For that, let T be a bounded closed rectangle in R 2 with 0 ∈ T and set N > 0 and define

T (N ) := N t : t ∈ T
the image of a fixed rectangle T by the dilatation t → N t and then letting T

R 2 is equivalent to N → ∞. Then, T (N )

Lipschitz-Killing curvatures of a given excursion set

Let u ∈ R. For X a real-valued stationary random field defined on R 2 , we consider the excursion set within T above level u:

{t ∈ T : X(t) ≥ u} = T ∩ E X (u)
, where E X (u) := X -1 ([u, +∞)).

Definition 1.1 (LK curvatures of E X (u)). Let X be a Gaussian field satisfying Assumption (A0).

Define the following Lipschitz-Killing curvatures for the excursion set E X (u), for u ∈ R,

L 2 (X, u, T ) : = |T ∩ E X (u)| = T 1 {E X (u)} (t) dt, (1)
L 1 (X, u, T ) : = |∂(T ∩ E X (u))| 1 = 1 2 lim ε→0 T δ ]u-ε,u+ε[ (X(t)) ∇X(t) dt, (2)
L 0 (X, u, T ) : = connected components in T ∩ E X (u) -holes in T ∩ E X (u), (3)
= 1 |T | (T ∩ C extr. ) -(T ∩ C saddle ) , (4)
where C extr. = t :

X(t) ≥ u, ∇X(t) = 0, t is a local extremum , C saddle = t : X(t) ≥
u, ∇X(t) = 0, t is a saddle point , δ u the Dirac mass at u and ∇X is the gradient field.

These three additive functionals, L j for j = 0, 1, 2 in Definition 1.1, called Lipschitz-Killing curvatures, are also closely related to intrinsic volumes and Minkowski functionals. Roughly speaking, for A a Borelian set in R 2 , L 0 (A) stands for the Euler characteristic of A, L 1 (A) for the half perimeter of the boundary of A and L 2 (A) is equal to the area of A, i.e., the two-dimensional Lebesgue measure.

In particular, when T is a bounded rectangle in R 2 with non empty interior,

L 0 (T ) = 1, L 1 (T ) = 1 2 |∂T | 1 , L 2 (T ) = |T |. (5)
Using the same formalism as in [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF], the normalized LK curvatures are given by

C /T i (X, u) := L i (X, u, T ) |T | , for i = 0, 1, 2,
and the associated LK densities by

C * i (X, u) := lim T R 2 E[C /T i (X, u)], for i = 0, 1, 2. ( 6
)
In the case i = 2, one can easily write 1), ( 2) and (3) can be evaluated in a given image by using the Matlab functions bwarea, bwperim and bweuler respectively. The numerical evaluation of (4) can be obtained for instance by finding local maxima, local minima and saddle points with the Matlab function imextrema. Definition 1.2. Let X be a Gaussian random field satisfying Assumption (A0). We define the effective observation level (effective level in the sequel)

C * 2 (X, u) = E[C /T 2 (X, u)] = P(X(0) ≥ u). Computationally C /T 2 (X, u) C /T 1 (X, u) and C /T 0 (X, u) in (
s u := u -µ σ
and the effective second spectral moment (effective spectral moment in the sequel) of the field X

w := λ σ 2 .
Notice that, if the underlying Gaussian field X is centered with unit variance, the effective level simply coincides with the observation level u and the effective spectral moment with the spectral moment of X. Using s u and w as in Definition 1.2, we compute the expected value of these empirical LK curvatures of excursion sets of a Gaussian random field with unknown location and scale and derive the corresponding LK densities by using the following adapted Gaussian Kinematic Formula.

Proposition 1.1. Let X be a Gaussian random field satisfying Assumption (A0), denote ψ(x) = P(N (0, 1) ≥ x). Then, it holds that

E[C /T 0 (X, u)] = ψ (s u ) |T | + √ w 2π exp - 1 2 s 2 u |∂T | 1 2|T | + w (2π) 3/2 exp - 1 2 s 2 u s u , E[C /T 1 (X, u)] = ψ (s u ) |∂T | 1 2|T | + √ w 4 exp - 1 2 s 2 u , E[C /T 2 (X, u)] = ψ (s u ) .
Having T R 2 , it follows that the LK densities defined in (6) are given by

C * 0 (X, u) = w (2π) 3/2 exp{- 1 2 s 2 u }s u , C * 1 (X, u) = √ w 4 exp{- 1 2 s 2 u }, C * 2 (X, u) = ψ (s u ) . (7)
Proof of Proposition 1.1 is postponed to Section A.2. The LK densities given by (7) enable to identify s u and w.

Remark 1 (Unbiased estimators for LK densities). Let u ∈ R and X as in Assumption (A0).

Assume we observe T ∩ E X (u) for T a rectangle in R 2 . The following quantities are unbiased estimator of C * i (X, u) in (7) (see also [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF]):

C 2,T (X, u) = C /T 2 (X, u), (8) C 1,T (X, u) = C /T 1 (X, u) - |∂T | 1 2|T | C /T 2 (X, u), (9)
C 0,T (X, u) = C /T 0 (X, u) - |∂T | 1 π|T | C /T 1 (X, u) + 1 2π |∂T | 1 |T | 2 - 1 |T | C /T 2 (X, u). ( 10
)
2 Estimators of the effective level and effective spectral moment

We focus on stationary Gaussian random fields with unknown location and scale parameters. In this section we build consistent estimators for s u and w in Definition 1.2 relying on (8), (10) and Proposition 1.1. The importance of estimating these quantities lies in the fact that a comparison between two images of excursion sets of two fields (with unknown and eventually different location and scale parameters) can be proposed only if the considered observation levels make sense in terms of their respective mean and variance. These effective level and effective spectral moment can be used to determine if such a comparison is meaningful.

Definition 2.1 (Estimators for s u and w). Define the estimator of s u built on the observation

T ∩ E X (u), u being fixed, by arg min s∈R |ψ(s) -C 2,T (X, u)|, i.e. s u,T := ψ -1 ( C 2,T (X, u)). ( 11
)
Let T 1 and T 2 be two rectangles in R 2 such that dist(T 1 , T 2 ) > 0 and

|T 1 | = |T 2 | > 0.
Define the estimator of w, u being fixed and u = µ, by

w u,T 1 ,T 2 := C 0,T 1 (X, u)(2π) 3/2 s u,T 2 exp{-1 2 ( s u,T 2 ) 2 } . ( 12
)
Remark 2. Notice that in Equation ( 12) we use two disjoint rectangles with same Lebesgue measure. To get the asymptotic normality of this quantity, we rely on a joint central limit theorem established for ( C 0,T 1 (X, u), C 2,T 2 (X, u)). However, if T 1 = T 2 such a result is not known, it can only be achieved for domains such that asymptotically dist(T 1 , T 2 ) → ∞ (see the proof of Proposition 2.2). Indeed, we are not aware of any joint central limit theorem for the LK densities of excursion sets of stationary random fields on R 2 , except in the case of a Boolean model in [START_REF] Hug | Second-order properties and central limit theorems for geometric functionals of Boolean models[END_REF] or in the case of pixelated images in [START_REF] Reddy | Central limit theorem for exponentially quasi-local statistics of spin models on cayley graphs[END_REF] and binary images in [START_REF] Ebner | Goodness-of-fit tests for complete spatial randomness based on Minkowski functionals of binary images[END_REF].

Such a constraint on T 1 and T 2 is satisfied for instance taking

T 1 ∪ T 2 ⊂ T such that T 1 ∩ T 2 = ∅, |T 1 | = |T 2 |, |T 1 | > c|T |, 0 < c < 1/2 and having T ↑ R 2 .
From a practical point of view, one only has one (large) image, the inference procedure can be implemented by partitioning it in sub-windows T 1 and T 2 separated by a "large" band (for instance of width √ size of the image) in order to build the estimator w u,T 1 ,T 2 .

Remark 3 (Bias for s u,T and w u,T 1 ,T 2 ). Estimator s u,T in ( 11) is the quantile of the standard Gaussian distribution at random level 1 -C 2,T (X, u), where C 2,T (X, u) defined in (8) is unbiased.

It follows that s u,T is biased, as the function ψ -1 is convex in (0, 0.5) (resp. concave in (0.5, 1)),

then it holds E[ s u,T ] > s u if C /T 2 (X, u) in (0,0.5) (resp. E[ s u,T ] < s u if C /T
2 (X, u) in (0.5,1)). The bias of s u,T certainly generates a bias also in the estimation of w. This bias for s u,T can be assessed in Figure 2. Remark 4. Let X be a stationary random field, non necessarily Gaussian, with finite variance.

It holds E[C /T 2 (X, u)] = P(X(0) ≥ u) = P X(0)-µ σ ≥ s u =: F (s u ).
Then, given a family for the marginal distribution F of the field, known up two its first two moments, the effective level of the field can always be recovered using s u = F -1 (C /T 2 (X, u)). This is not derived from the kinematic formula contrary to the equation leading to the computing of w.

Proposition 2.1 (Asymptotic normality of ( s u 1 ,T , . . . , s um,T )). Let X be a Gaussian random field satisfying Assumptions (A0) and (A1). For a positive integer N , consider T (N ) = {N t : t ∈ T } and s u i ,T (N ) the estimator defined in (11) built on the observation T (N ) ∩ E X (u i ), where u 1 , . . . , u m are fixed. Then, |T (N ) |( s u 1 ,T (N ) -s u 1 , . . . , s um,T (N ) -s um ) converges in distribution to a centered Gaussian vector with covariance matrix (Σ 2 s,(u i ,u j ) ) 1≤i,j≤m , with

Σ 2 s,(u i ,u j ) = R 2 ρ(t) 0 h (u i ,u j ) (r)dr dt ∈ (0, +∞), where h (u i ,u j ) (r) = 1 √ 1 -r 2 exp 2rs u i s u j -r 2 s 2 u i -r 2 s 2 u j 2(1 -r 2 ) .
Proposition 2.1 is a consequence of Proposition A.2 and a multidimensional version of the delta method.

In the sequel we denote by σ 2 su := Σ 2 s,(u,u) i.e.,

σ 2 su = R 2 ρ(t) 0 1 √ 1 -r 2 exp s 2 u r 1 + r dr dt ∈ (0, +∞). ( 13
)
The function s 2 u → σ 2 su is an increasing function, it is minimal when s 2 u is minimal i.e. at u = µ. This leads to the following result.

Corollary 2.1 (Asymptotic variance for s u=µ,T ). Under assumptions of Proposition 2.1, the variance σ 2 su in (13) attains its minimum at u = µ, i.e., s u = 0, and this minimum is

σ 2 sµ = R 2 arcsin(ρ(t)) dt.
Corollary 2.1 will be useful in Sections 4.1 and 4.2, where we frequently choose the observation level u such that the associated effective level s u is null. This procedure will guarantee that the estimator has minimal variance. Corollary 2.1 is consistent with recent results on the so-called Berry's cancellation phenomenon, i.e., the fact that the variance of the LK curvatures at level s u reaches a minimum, because of the disappearance of a large number of chaotic terms (see [START_REF] Cammarota | A reduction principle for the critical values of random spherical harmonics[END_REF]).

Proposition 2.2 (Asymptotic normality of ( w u 1 , . . . , w um )). Let X be a Gaussian random field satisfying Assumptions (A0), (A1) and (A2). Let T 1 and T 2 be two rectangles in R 2 such that

dist(T 1 , T 2 ) > 0 and |T 1 | = |T 2 | > 0. For a positive integer N , let T (N ) i = {N t : t ∈ T i }, i ∈ {1, 2}. Consider w u i ,T (N ) 1 ,T (N ) 2
the estimator defined in (12) with u i = µ being fixed, for i = 1, . . . , m. Then,

|T (N ) 1 | w u 1 ,T (N ) 1 ,T (N ) 2 -w, . . . , w um,T (N ) 1 ,T (N ) 2
-w converges in distribution to a centered Gaussian vector with covariance matrix (Σ 2 w,(u i ,u j ) ) 1≤i,j≤m given by

Σ 2 w,(u i ,u j ) = (2π) 3 Σ 2 C * 0 ,(u i ,u j ) + Σ 2 s,(u i ,u j ) (s 2 u i -1)(s 2 u j -1)C * 0 (X, u i )C * 0 (X, u j ) s 2 u i s 2 u j exp{- s 2 u i 2 + s 2 u j 2 } < +∞, with Σ 2 C * 0 ,(u i ,u j )
as in Proposition A.1 and Σ 2 s,(u i ,u j ) as in Proposition 2.1.

Proof of Proposition 2.2 is postponed to Section A.2.

Remark 5. The fact that Σ 2 C * 0 ,(u,u) is nonzero for all levels u still is an open problem. In [START_REF] Estrade | A Central Limit Theorem for the Euler characteristic of a Gaussian excursion set[END_REF] the first Itô-Wiener chaos element of the series of Σ 2 C * 0 ,(u,u) has been calculated for all levels u. However, to guarantee that it is nonzero, higher order elements of this series need to be explored. Conversely, the first Itô-Wiener chaos element of the series of σ 2 C * 2 ,u (see, e.g., [START_REF] Müller | A central limit theorem for Lipschitz-Killing curvatures of Gaussian excursions[END_REF] and [START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF]) guarantees it is non degenerate for all u ∈ R.

Numerical illustrations

In the following we provide an illustration of the finite sample performance of the proposed estimator s u,T for several values of u (see Figure 2). We consider here M = 100 sample simulations of Gaussian random fields as in Assumption (A0) with µ = 12, σ 2 = 4 and covariance r(x) = e -κ 2 x 2 , for κ = 100/2 10 in domains of size 2 10 × 2 10 pixels.

As expected from Proposition 2.1, the quality of the consistency can be evaluated in the displayed boxplots. Furthermore the associated asymptotic empirical variance is analysed and compared with the theoretical one in the case of u → σ 2 su in (13) (see right panel in Figure 2). Furthermore a bias in the estimation of s u can be observed for values of level u far form µ (see Remark 3). Similar results have been obtained for the estimator of w u,T 1 ,T 2 but they are omitted here for the sake of brevity.

3 Subwindow empirical variance estimation

Consistent variance estimator

Having asymptotically normal estimators for s u or w is insufficient if we cannot estimate their limit variance. Indeed, this variance plays a central role for instance to calibrate testing procedures, such as a test to determine if two images of excursion sets are comparable (see Section 4). Due to the specificity of our observations, which are excursion sets, proving that a variance estimator is consistent turns out to be a delicate, while fundamental, problem Definition of the subwindow estimators We introduce a technique to estimate the asymptotic variances appearing in Propositions 2.1 and 2.2. Inspired by the cutting of T (N ) introduced in Pantle et al. ( 2010) (see also [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF], Section 5) we consider a classical empirical variance estimator. To establish its consistency, we decompose this estimator on domains that are infinitely distant, mimicking the classical context of independent and identically distributed random variables.

Consider the following cutting of T (N ) : set M N ∈ N such that M N → ∞ as N → ∞ and consider the grid pattern of T (N ) defined by (V (N,(i,j)) ) 1≤i,j≤M N -such that ∪ 1≤i,j≤M N V (N,(i,j)) = T (N ) , for (i, j) = (i , j ), V (N,(j,j)) ∩ V (N,(i ,j )) = ∅ and |V (N,(i,j)

) | := r N , ∀ 1 ≤ i, j ≤ M N where r N → ∞ as N → ∞. Observe that, whenever max{|i -i |, |j -j |} ≥ 2 then dist(V (N,(i,j)) , V (N,(i ,j )) ) → ∞
as N → ∞. Then, by Propositions A.1 and A.2 the estimators C /V (N,(i,j)) (X, u), for ∈ {0, 2}, are consistent and asymptotically independent. Therefore, we estimate the asymptotic variances appearing in Propositions A.1 and A.2 by

Σ 2 C * ,(u,v) = 1 M 2 N -1 M N i,j=1 ξ (i,j) N (u) ξ (i,j) N (v) - 1 M 2 N -1 M N i,j=1 ξ (i,j) N (u) 1 M 2 N -1 M N i,j=1 ξ (i,j) N (v) , ( 14
)
where we define ξ

(i,j) N (u) := C /V (N,(i,j)) (X, u), for all u ∈ R, ∈ {0, 2}, 1 ≤ i, j ≤ M N .
The key point is that these estimators are identically distributed, as the field X is stationary, and are asymptotically independent whenever max{|i -i |, |j -j |} ≥ 2 in view of Propositions A.1 and A.2 and the fact that dist(V (N,(i,j)) , V (N,(i ,j )) ) → ∞, as N → ∞.

Case of the area ( = 2 in Equation ( 14)) To the knowledge of the authors, the only result establishing the consistency of ( 14), is Theorem 6 of [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF] in the case of the area ( = 2). However, this results imposes a condition (39) on the fourth-order cumulant that is difficult to check and for which the example provided is "a random field X with finite dependence range", which means that ρ = ρ1 D for some finite domain D. In the following Theorem 3.1 we prove consistency of the estimator Σ 2 C * 2 ,(u,v) defined in ( 14) under Assumption (A1).

Theorem 3.1. Let X be a Gaussian random field satisfying Assumptions (A0) and (A1). Let

Σ 2 C * 2 ,(uq,u k ) be as in (14) with q, k ∈ {1, . . . , m}. Then, it holds that Σ 2 C * 2 ,(uq,u k ) P -→ N →∞ Σ 2 C * 2 ,(uq,u k ) .
Proof of Theorem 3.1. We show the consistency of estimators Σ 2 C * 2 ,(uq,u k ) in ( 14), for all q, k ∈ {1, . . . , m}. First, note that Σ 2 C * 2 ,(uq,u k ) is given by the difference of two terms, write

Σ 2 C * 2 ,(uq,u k ) := S 1 (u, v) -S 2 (u)S 2 (v). ( 15
)
By Proposition 2.1 and the additivity of T → C /T 2 (X, u), we get that for all u, S 2 (u)

P ----→ N →∞ C * 2 (X, u). (2i, 2j) (2i, 2j + 1) (2i, 2j + 2) (2i, 2j + 3) (2i + 1, 2j) (2i + 1, 2j + 1) (2i + 1, 2j + 2) (2i + 1, 2j + 3) (2i + 2, 2j) (2i + 2, 2j + 1) (2i + 2, 2j + 2) (2i + 2, 2j + 3) S 1,1 S 1,1 S 1,1 S 1,1 S 1,2 S 1,2 S 1,2 S 1,2 S 1,3 S 1,3 S 1,4 S 1,4
Figure 3: Decomposition of the sum S 1 in ( 15) in four sums.

The remaining of the proof consists is noticing that the sum S 1 can be rewritten as four sums each composed of identically distributed terms that are asymptotically independent. Decompose S 1 as follows (see Figure 3) and without loss of generality assume that M N is even, then

S 1 = 1 M 2 N -1 M N i,j=1 ξ (2i,2j) N (u q ) ξ (2i,2j) N (u k ) + 1 M 2 N -1 M N i,j=1 ξ (2i,2j+1) N (u q ) ξ (2i,2j+1) N (u k ) + 1 M 2 N -1 M N i,j=1 ξ (2i+1,2j) N (u q ) ξ (2i+1,2j) N (u k ) + 1 M 2 N -1 M N i,j=1 ξ (2i+1,2j+1) N (u q ) ξ (2i+1,2j+1) N (u k ) : = S 1,1 + S 1,2 + S 1,3 + S 1,4 . To show that S 1,p -1 4 E[ ξ (1,1) N (u) ξ (1,1) N (v)] P ----→ N →∞
0, 1 ≤ p ≤ 4, we show that the convergence holds in L 2 . We perform computations for S 1,1 , other sums are treated similarly. It holds, using the stationarity of the field, that

V(S 1,1 ) = 1 (M 2 N -1) 2 M N i,i ,j,j =1 Cov ξ (2i,2j) N (u q ) ξ (2i,2j) N (u k ), ξ (2i ,2j ) N (u q ) ξ (2i ,2j ) N (u k ) = 1 (M 2 N -1) 2 M N (i,j) =(i ,j )=1 Cov ξ (2i,2j) N (u q ) ξ (2i,2j) N (u k ), ξ (2i ,2j ) N (u q ) ξ (2i ,2j ) N (u k ) (16) + 1 (M 2 N -1) 2 M N i,j=1
Cov ξ

(2i,2j) N (u q ) ξ (2i,2j) N (u k ), ξ (2i,2j) N (u q ) ξ (2i,2j) N (u k ) . ( 17
) For all (i, j) ∈ {1, . . . , M N } 2 , if we set G = (X -µ)/σ it holds that ξ (i,j) N (u) = C /V (N,(i,j)) 2 (X, u) = C /V (N,(i,j)) 2 (G, s u ) = L 2 (G, s u , V (N,(i,j)) ) |V (N,(i,j)) | .
The remaining of the proof is a direct consequence of the following result.

Proposition 3.1. Let X be a Gaussian random field satisfying Assumptions (A0) and (A1) and set G = (X -µ)/σ. Let u 1 , u 2 , u 3 , u 4 be four fixed levels in R and associated effective levels

s u i = (u i -µ)/σ, for i ∈ {1, 2, 3, 4}. Let T and T be such that |T | = |T | and dist(T, T ) → ∞.
Then, it holds that,

E[L 2 (G, s u 1 , T )L 2 (G, s u 2 , T )L 2 (G, s u 3 , T )L 2 (G, s u 4 , T )] = ψ(s u 1 )ψ(s u 2 )ψ(s u 3 )ψ(s u 4 )|T | 4 + o(|T | 3 ).
Proof of Proposition 3.1 is postponed to Section A.3. Proposition 3.1 ensures that ( 16) goes to 0.

Let (i, j) = (i , j ) ∈ {1, . . . , M N } 4 and T = V (N,(i,j)) and T = V (N,(i ,j )) , for which for N large enough it holds dist(V (N,(i,j)) , V (N,(i ,j )) ) > 1 and

|V (N,(i,j)) | = |V (N,(i ,j )) | → ∞ as N → ∞, then we derive lim n→∞ Cov ξ (2i,2j) N (u q ) ξ (2i,2j) N (u k ), ξ (2i ,2j ) N (u q ) ξ (2i ,2j ) N (u k ) = 0.
This implies that ( 16) goes to 0 as N → ∞ using Cesàro Lemma. Moreover, Equation ( 17) vanishes

using that | ξ N | ≤ 1 almost surely, implying (17) is bounded by M 2 N /(M 2 N -1) 2 → 0, as N → ∞.
This implies that V(S 1,1 ) → 0 as N → ∞ using Cesàro Lemma. The same behaviour holds true for V(S 1,p ) for all p ∈ {1, 2, 3, 4}, implying the desired result in Theorem 3.1.

Discussion Theorem 3.1 establishes consistency of this subwindow variance estimator in the case of the area under Assumption (A1). This assumption seems weak as it is the same assumption as the one imposed to get the asymptotic normality of C 2 (X, u) (see Theorem 4 of [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF]).

A key element in the proof of Theorem 3.1 is Proposition 3.1 which is similar to the constraint (39) of [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF]. More precisely, in Proposition 3.1 we control covariance terms between distant domains. Proof of Proposition 3.1 relies on two ingredients: i. the Itô-Wiener chaos decomposition of L 2 (X, u, T ) and ii. the diagram formula (see [START_REF] Taqqu | Law of the iterated logarithm for sums of non-linear functions of Gaussian variables that exhibit a long range dependence[END_REF], Lemma 3.2) applied up to the dimension 4. This diagram formula in dimensions 3 and 4 leads to technical and lengthly computations where the constraint dist(V (N,(i,j)) , V (N,(i ,j )) ) → ∞ as N → ∞, and the cutting of T (N ) , plays a crucial part.

A key question from an applied point of view is the choice of M N and |V (N,(1,1) | (see Section 5 of [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF]). Theoretically, it is only required is that M N and |V (N,(1,1) | increase to infinity, but in practice, we need a cutting of T in infinitely many rectangles such that in each subrectangle the excursion set is sufficiently rich to properly estimate ξ (i,j) N . A coarse comparison of the different order of magnitudes required to get a normalized central limit theorem using Σ 2

C * 2 (u,v) is to consider M N |T | 1/4 and |V (N,(1,1) | |T |.
The result below is a direct consequence of Theorem 3.1: the consistency of the estimator for Σ 2 s,(u i ,u j ) , introduced in Proposition 2.1.

Corollary 3.1. Let X be a Gaussian random field satisfying Assumptions (A0) and (A1). Define, the estimator of Σ 2 s,(u i ,u k ) , where (i, k) ∈ {1, . . . , m} 2 , by

Σ 2 s,(u i ,u k ) = 2π exp 1 2 ( s 2 u i + s 2 u k ) Σ 2 C * 2 ,(u i ,u k ) , (18)
with s u := s u,T (N ) as in Proposition 2.1. Then, it holds that, for

(i, k) ∈ {1, . . . , m} 2 , Σ 2 s,(u i ,u k ) P -→ N →∞ Σ 2 s,(u i ,u k ) .
Case of the Euler characteristic ( = 0 in Equation ( 14)) In the case of the Euler characteristic ( = 0), the same proof should generalize under the following additional assumption.

(A3) Suppose that Assumption (A2) holds true with |M r (t

)| ≤ C M /(1 + t ) τ where C M is a positive constant and τ > 2.
Indeed there is a Itô-Wiener chaos decomposition of L 0 (X, u, T ) (see the Additional Supporting Information) with similar properties on the coefficients. However, computing the fourth moment with this chaos decomposition would imply to use the diagram formula in dimensions 18 and 24, which become extremely technical and burdensome. Therefore, we give in the Additional Supporting Information some ingredients to establish the following Conjecture 3.1 which is an equivalent of Proposition 3.1 in the case of the Euler characteristic having for all (i, j) ∈ {1, . . . , M N } 2 , ξ (i,j)

N (u) = C /V (N,(i,j)) 0 (X, u) = L 0 (X, u, V (N,(i,j)) ) |V (N,(i,j)) | .
The missing step in the proof to have the result is the computation of the diagram formula in dimensions 18 and 24.

Conjecture 3.1. Let X be a Gaussian random field satisfying Assumptions (A0) and (A3). Let T and T be such that |T | = |T | and dist(T, T ) → ∞. Let u 1 , u 2 , u 3 , u 4 be four fixed levels in R.

Then it holds that

E[L 0 (X, u 1 , T )L 0 (X, u 2 , T )L 0 (X, u 3 , T )L 0 (X, u 4 , T )] = w 4 (2π) 6 exp{- 1 2 (s 2 u 1 + s 2 u 2 + s 2 u 3 + s 2 u 4 )}s u 1 s u 2 s u 3 s u 4 |T | 4 + O(|T | 3 ),
where w denotes the effective spectral moment of X in Definition 1.2.

Mimicking the proof of Theorem 3.1, we derive that a consequence of Conjecture 3.1 is that

Σ 2 C * 0 ,(uq,u k ) P -→ N →∞ Σ 2 C * 0 ,(uq,u k ) . (Conjecture)
Then, using (Conjecture), one can conjecture the consistency of the estimator for Σ 2 w,(u i ,u j ) , introduced in Proposition 2.2. Conjecture 3.2. Let X be a Gaussian random field satisfying Assumptions (A0), (A1) and (A3).

Define the estimator of Σ 2 w,(u i ,u k ) , where

(i, k) ∈ {1, . . . , m} 2 , Σ 2 w,(u i ,u k ) = (2π) 3 Σ 2 C * 0 ,(u i ,u k ) + Σ 2 s,(u i ,u k ) ( s 2 u i -1)( s 2 u k -1) C 0,T (N ) (X, u i ) C 0,T (N ) (X, u k ) s 2 u i s 2 u k exp - s 2 u i 2 + s 2 u k 2 ,
with s u := s u,T (N ) as in Proposition 2.1, Σ 2 s,(u i ,u k ) as in (18) and C 0,T (N ) (X, u k ) as in (10). Then,

Σ 2 w,(u i ,u k ) P -→ N →∞ Σ 2 w,(u i ,u k ) .

Statistical implications and numerical illustrations

In this section we consider several useful applications of theoretical results provided in Sections 2 and 3.1. Indeed, the asymptotically Gaussian estimators of the unknown location and scale parameters and the subwindow empirical variance estimators are used in the following to build a test of Gaussianity using LK densities (see Section 3.2.1). A procedure to build an interval for the unknown location µ of the field is given in Section 3.2.2.

A Gaussianity test for non standard fields

Here we generalize the test in Section 3.1 of [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF] in cases where the field is not supposed to be centered and with unit variance. In this section we are interested in testing H 0 : X is a Gaussian field with unknown mean and unknown variance.

Consider the quantities

R(u 1 , u 2 ) := C * 0 (X, u 2 ) C * 0 (X, u 1 ) under H 0 = s u 2 s u 1 exp - 1 2 (s 2 u 2 -s 2 u 1 ) := R H 0 (u 1 , u 2 ).
This ratio is empirically accessible provided the field is observed at two distinct levels u 1 and u 2 ,

R T (u 1 , u 2 ) := C 0,T (X, u 2 ) C 0,T (X, u 1 ) . ( 19
)
We establish below a central limit theorem for R T (u 1 , u 2 ). Notice that one can readily establish

a CLT for |T | R T (u 1 , u 2 ) -R H 0 (u 1 , u 2 ) using,
R H 0 T (u 1 , u 2 ) := s u 2 s u 1 exp - 1 2 s 2 u 2 -s 2 u 1 . ( 20
)
For technical reasons, we estimate R T 1 (u 1 , u 2 ) and R H 0 T 2 (u 1 , u 2 ) on rectangles T 1 and T 2 that are asymptotically infinitely distant so that both excursion sets are asymptotically independent and lead to independent estimators (see also Remark 2). Proposition 3.2. Assume that X is a Gaussian field satisfying Assumptions (A0), (A1) and (A2).

Let u 1 , u 2 = µ being fixed. Let T 1 and T 2 be two rectangles in R 2 such that dist(T 1 , T 2 ) > 0 and

|T 1 | = |T 2 | > 0.
For any positive integer N , we define T

(N ) i = {N t : t ∈ T i }, for i = 1, 2. Then, under H 0 it holds that |T (N ) 1 | R T (N ) 1 (u 1 , u 2 ) -R H 0 T (N ) 2 (u 1 , u 2 ) d ----→ N →∞ N (0, σ 2 Ru 1 ,u 2 ), where σ 2 Ru 1 ,u 2 < ∞ and R T (N ) 1 (u 1 , u 2 ) (resp. R H 0 T (N ) 2 (u 1 , u 2 )) is defined as in (19) (resp. in (20)) built on the observation T (N ) 1 ∩ E X (u 1 ) and T (N ) 1 ∩ E X (u 2 ) (resp. T (N ) 2 ∩ E X (u 1 ) and T (N ) 2 ∩ E X (u 2 )).
Proof of Proposition 3.2 is postponed to Section A.2.

Subwindow empirical estimation of variance σ

2 Ru 1 ,u 2
We provide a normalized version of Proposition 3.2, we build a consistent estimator for σ 2

Ru 1 ,u 2 , i.e., σ 2 Ru 1 ,u 2 = σ 2 g(C * 0 ),(u 1 ,u 2 ) + σ 2 h(s),(u 1 ,u 2 ) := V( R T (N ) 1 (u 1 , u 2 )) |T (N ) 1 | + V( R H 0 T (N ) 2 (u 1 , u 2 )) |T (N ) 2 |. (21) Consistent estimators for V( R T (N ) (u 1 , u 2 )) and V( R H 0 T (N ) (u 1 , u 2
)) are obtained using the subwindows technique of Section 3.1. Define σ 2 h(s),(u 1 ,u 2 ) in ( 21) by

σ 2 h(s),(u 1 ,u 2 ) = e -( s 2 u 2 -s 2 u 1 ) σ 2 su 2 1-s 2 u 2 su 1 2 + 2 Σ 2 s,(u 1 ,u 2 ) 1-s 2 u 2 su 1 s u 2 - su 2 s 2 u 1 + σ 2 su 1 s u 2 - su 2 su 1 2 2 , with s u i := s 2 u i ,T (N ) 2
defined in (11), Σ 2 s,(u 1 ,u 2 ) and σ 2 su i as in Corollary 3.1. From Proposition 2.1 and Corollary 3.1, it holds that σ 2 h(s),(u 1 ,u 2 )

P -→ N →∞ σ 2 h(s),(u 1 ,u 2 ) . Similarly, define σ 2 g(C * 0 ),(u 1 ,u 2 ) σ 2 g(C * 0 ),(u1,u2) = σ 2 C * 0 ,u2 C 0,T (N ) 1 (X, u 1 ) 2 -2 C 0,T (N ) 1 (X, u 2 ) C 0,T (N ) 1 (X, u 1 ) 3 Σ 2 C * 0 ,(u1,u2) + C 0,T (N ) 1 (X, u 2 ) 2 C 0,T (N ) 1 (X, u 1 ) 4 σ 2 C * 0 ,u1 , with C 0,T (N ) 1 (X, u i ) as in (10) and Σ 2 C * 0 ,(u 1 ,u 2 )
as in (Conjecture). Take a confidence level α ∈ (0, 1) and set q 1-α 2 such that P(|N (0, 1

)| ≤ q 1-α 2 ) = 1 -α 2 .
We define the symmetric test φ T (N ) with asymptotic level α as

φ T (N ) 1 ,T (N ) 2 = 1    |T (N ) 1 | σ 2 Ru 1 ,u 2 R T (N ) 1 (u 1 ,u 2 )-R H 0 T (N ) 2 (u 1 ,u 2 ) ≥q 1-α 2    , (22)
where

σ 2 Ru 1 ,u 2
is the consistent variance estimator built before.

Comments The test statistic is well defined for all levels u = µ, similarly to [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF] where we could not use the level u = 0. Here µ is unknown and choosing u 1 close to µ will make the ratio R T (u 1 , u 2 ) unstable. A solution is to modify the test statistics and consider the difference R T (u 1 , u 2 ) = C 0,T (X, u 2 ) -C 0,T (X, u 1 ) instead. The test could be modified accordingly at the expense of more tedious calculations, in particular due to the fact that the unknown quantity w := λ σ 2 will appear in the test statistics R T (u 1 , u 2 ) under H 0 .

Furthermore, under H 0 we do not impose any constraint on the shape of the covariance function nor on the spectral moment other than Assumptions (A1) and (A2). However, contrary to the case where the field is centered with unit variance (see [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF]), the test statistic under H 0 depends on unknown quantities (µ, σ 2 ).

Establishing the consistency of the test, even for particular alternatives H 1 , is a delicate problem as it requires the knowledge of a limit theorem for C 0,T (X, u) as well as a good understanding of the value R(u 1 , u 2 ) under the alternative. One possible hypothesis H 1 is

H 1 : ∃ k ≥ 3, X is Student random field with k unknown degrees of freedom.
Under H 1 , we can derive the ratio

R(u 1 , u 2 ) := C * 0 (X, u 2 ) C * 0 (X, u 1 ) under H 1 = u 2 u 1 k + u 2 2 k + u 2 1 1-k 2 := R H 1 (u 1 , u 2 ).
As explained in Section 3.1 in [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF], if a joint central limit theorem for Lipschitz-Killing curvatures for Student fields was known this would entail the consistency of the test in ( 22), i.e., P H 1 (φ

T (N ) 1 ,T (N ) 2 = 1) -→ N →∞ 1.

Interval containing the unknown location µ

We propose the following construction for an interval containing the unknown location µ of the field based on the observations T ∩ E X (u 1 ), for J fixed levels u 1 < . . . < u J .

Procedure to build interval for the unknown location of the field Input :

Let J ∈ N be fixed, the field X is observed at levels u 1 < . . . < u J .

Estimation :

Let σ 2 su j := Σ 2 s,(u i ,u i ) , with j ∈ {1, . . . , J} and Σ 2 s,(u i ,u i ) as in (18). Define j = argmin j∈{1,...,J} σ 2 su j and

j ±1 =        j + 1 if σ 2 su j+1 ≤ σ 2 su j-1 j -1 if σ 2 su j-1 ≤ σ 2 su j+1 ,

Final output :

The following interval contains µ with large probability:

I µ := min{u j , u j ±1 }, max{u j , u j ±1 } with convention that if j = 0, I µ := [-u 1 , u 0 ] and if j = J, I µ := [u J , -u J-1 ]. Proposition 3.3. It holds that P µ ∈ min{u j , u j ±1 }, max{u j , u j ±1 } -→ T R 2 1.
Proof of Proposition 3.3 is postponed to Section A.2.

Test to compare two images of excursion sets

Let Y and Z be two Gaussian fields satisfying Assumptions (A0) and (A1) with possibly different mean, variance, spectral moment or correlation function. Suppose one has two images of the excursion sets of these fields observed at the respective levels u Y and u Z and wants to know if it is possible to compare the Lipschitz-Killing curvatures of theses two images, i.e. one want to test whether the effective level of Y (denoted by s u Y (Y )) is equal to the effective level of Z (denoted by s u Z (Z)):

H 0 : s u Y (Y ) = s u Z (Z) H 1 : s u Y (Y ) = s u Z (Z).
Introduce the quantity S(Y, Z)

:= s u Y -s u Z and S T (N ) (Y, Z) := s u Y ,T (N ) -s u Z ,T (N ) its empirical counterpart.
Notice that S(Y, Z) = 0 under the null hypothesis. We can now state the autonormalised central limit theorem for the test statistics S T (N ) (Y, Z).

Corollary 4.1. Let X be a Gaussian random field satisfying Assumptions (A0) and (A1). It holds

that 1 Σ Y,Z S T (N ) (Y, Z) d, H 0 ----→ N →∞ N (0, 1), where Σ Y,Z = σ 2 s,u Y + σ 2 s,u Z , with σ 2 s,u Y := Σ 2 s,(u Y ,u Y ) ( resp. σ 2 s,u Z := Σ 2 s,(u Z ,u Z )
), the consistent estimator of the variance in (18).

Corollary 4.1 is proved by using Proposition 2.1, Corollary 3.1 and Theorem 3.1. Notice that due to the independence of the two considered images, we get a simplified expression for the asymptotic variance.

Let q 1-α 2 be such that P(|N (0, 1)| ≤ q 1-α 2 ) = 1 -α 2 . Finally, we introduce the symmetric test φ T (N ) with asymptotic level α:

φ T (N ) = 1 1 Σ Y,Z S T (N ) (Y,Z) ≥ q 1-α 2 . ( 23
)
Clearly the test proposed in Equation ( 23) is consistent: under

H 1 : s u Y (Y ) = s u Z (Z), the test statistics S T (N ) (Y, Z) diverges.
In the following, we evaluate the performances of the proposed test φ T (N ) in ( 23) on simulated mammograms (from a digital texture model, see Section 4.1) and on real mammograms (from mini-MIAS database, see Section 4.2). In particular, we aim to test whether the effective level of an image is equal to the effective level of another one.

Comparing images of excursion sets: a synthetic mammograms study

The data-set In this section we consider images from a recent solid breast texture model inspired by the morphology of medium and small scale fibro-glandular and adipose tissue observed in clinical breast computed tomography (bCT) images (UC Davis database). Each adipose compartment is modeled as a union of overlapping ellipsoids and the underlying dynamic is dictated by a spatial marked point process. The contour of each ellipsoid is blurred to render the model more realistic The first image from each group is reported in Figure 4. As observed in Section 3 in [START_REF] Li | A novel 3D stochastic solid breast texture model for X-Ray breast imaging[END_REF], these simulated digital mammograms from groups (F), (FG) and (D) show a high visual realism compared to real images in these 3 different clinical situations (see also Figure 8).

Importance of the effective level In Figure 5 (left) we estimate the effective levels s u for different values of u. One can appreciate that the estimated effective levels are different for each group of mammograms. Moreover, it seems there exists an order between the three groups of mammograms in terms of effective levels. This means that if one want to compare images from different groups, one cannot use the same value of u to perform the comparison. Indeed, in Figure 5 (right-up) we display the excursion sets of the first image of each groups (denoted image 1.F, 1.FG and 1.D in Figure 4) for the same fixed level u = 2200. The resulting excursion sets look completely different from one group to the other. However this difference is not necessarily the result of a intrinsic difference between the images, but a problem in the calibration of the level u used for the comparison. On this example, as we observed in Figure 5 (left) a difference between the groups of effective levels, we calibrate a level u for each group such that they all have an effective level close to 0, i.e., adaptive levels u, such that | s u | < , for = 10 -2 . Notice that this choice guarantees a minimal variance for the estimated effectively level (see Corollary 2.1). Figure 5 We observe now that the excursion sets in Figure 5 (right-down) look visually similar. Moreover, at these levels u, we theoretically expect to get estimated Euler characteristic values close to zero (see Equation ( 7)). This behavior is illustrated in the boxplots gathered in Figure 6. Then, any geometric estimated quantity (as for instance the Lipschitz-Killing curvatures) on these excursion sets cannot be compared at the same levels but at the same effective levels.

Testing two images of excursion sets We now test for the 1.F, 1.FG, 1.D images in Figure 4:

H 0 : s u Y (Y ) = s u Z (Z) versus H 1 : s u Y (Y ) = s u Z (Z), (24)
for Y, Z ∈ {1.F,1.FG, 1.D}, where u Y and u Z are the adaptive levels previously defined and such that | s u | < 10 -2 , i.e., the associated C 0,T ( u) is zero (see right panel of Figure 5 and boxplots in Figure 6).

Indeed, for these levels u 1.F , u 1.F G , u 1.D , the p-values of the test ( 24), gathered in We now consider 1000 different values of u and for each u we perform the following test

: s u Y (Y ) = s u Z (Z)
H 0 : s u (Y ) = s u (Z) versus H 1 : s u (Y ) = s u (Z), (25)
for Y, Z images of this synthetic mammograms data-set. Notice that, contrary to the test (24), we did not previously choose adaptive levels. Then the test ( 25) is performed for possible different effective levels s u . In Table 2 we display the number of p-values associated to the 1000 values of u that are smaller than the significant level α = 0.2. In particular we perform an intra-class analysis (left panel in Table 2) and an inter-classes analysis (right panel in

Comparing images of excursion sets: a real digital mammograms study

The data-set The Mammographic Image Analysis Society (MIAS) is an organisation of UK research groups interested in the understanding of mammograms and has produced a database of real digital mammograms. The range of intensity in all images is represented from 0 to 255 and we consider images size of 250 × 250 pixels. We study 211 mammograms classified in terms of the character of background tissue: fatty tissue group (F) (66 images), fatty-glandular (FG) (67 images), and dense (D) (77 images). Mammographic images are available online: http://peipa.essex.ac.

uk/info/mias.html (see also [START_REF] Suckling | The Mammographic Image Analysis Society Digital Mammogram Database[END_REF]). One image of each group is displayed in Figure 8. Inference and testing for effective level For this real mammograms data-set, similarly to

Figure 7, we perform the test in ( 25) and we display in bold marked points the cases when the test rejects H 0 for a significant level α = 0.2 (see Figure 9 below). These points are drawn on the estimates s u , for the considered values of u. For the sake of brevity we only display here the inter-classes analysis. Notice that in the first panel of Figure 9 we choose the difficult comparison between the two closest images between groups F and FG. Remark that the H 0 hypothesis is accepted for almost all the levels u in the considered grid in the first and the last panels of Figure 9, due to the visible proximity of the effective levels s u . In the second, third and fourth panel, the H 0 hypothesis is rejected for all u and the associated excursion sets can not be compared without a preliminary image processing. The fifth panel in Figure 9 represents an interesting hybrid situation. Almost everywhere H 0 is rejected except for a small interval of u's values where H 0 is accepted. In this small range of u, the relative excursion sets can be considered visually similar.

Analogously to Table 2, we now perform the test ( 25 As one can expect, in Figure 10 we observe differences in terms of the medians and of the variances:

the intra-classes boxplots (first three boxplots) have smaller median and variance values with respect to the inter-classes ones (last three boxplots). As in Section 4.1, this study suggests that the effective level scaling procedure seems to be necessary especially in the comparison of two images of excursion sets belonging to two different background tissue groups.

Conclusions and discussion

In this paper, we have presented new statistics based on the average LK curvatures of the excursion set of a stationary non-standard isotropic Gaussian field X on R 2 , in particular on the Euler characteristic, the half perimeter and the area. These tools allow to build consistent inference procedures based only on a sparse observation of the Gaussian 2D random field with unknown location and scale parameters. A byproduct is the construction of a test to determine if two images of excursion sets can be compared. Here we discuss some potential improvements of the results

proposed in this work.

Firstly, notice that in Sections 3.2 and 4, our testing procedure results relied mostly on the area and not on the Euler characteristic devise. The reason being that for the Euler characteristic we did not fully establish the consistency of its asymptotic variance estimator, therefore limiting the possibility to properly calibrate statistical procedures. If Conjecture 3.1 was proved, this would enlarge the options for testing and inference. Besides, if a joint auto-normalized central limit theorem for

C /T 0 (X, u), C /T
2 (X, u) was available, this would imply results on the joint behavior of ( s u,T , w u,T ) in Definition 2.1 and the possibility to build consolidated tests using both quantities.

Secondly, inspired by the analysis of Sections 4.1 and 4.2, the effective level could be useful to build a "classification criterion" between the three groups (F), (FG) and (D). Indeed, the estimated s u seems globally able to distinguish images coming from different groups relying exclusively on the sparse information of the excursion set with a same level u.

To explore this idea, in Figure 11 we represent for the real digital data-set the adaptive level u (y axis) for all considered 211 images such that the estimated effective level s u is approximated equal to δ, for several values of δ and where the different groups are distinguished in notation. Visually, it appears that for small values of δ (i.e., for intermediate values of u) the groups might be classified from their effective level. Obviously, for large values of δ (i.e., for extreme values of u) the quality of the separation of groups is less good, as one can expect (see last panel of Figure 11). Moreover, one can expect that for large values of the level u, the obtained classes will be more distant also due to the large variance of the effective level estimator (see Corollary 2.1). However globally, this tool seems adequate to recognize the underlying background breast tissue. However an investigation on a rigorous classification procedure based on effective levels remains an open point and is left for future works.

Finally, in the present paper, the empirical estimator of the variance is obtained from the cutting of the domain T (N ) , as described in Section 3.1 and illustrated in Figure 3. While the theoretical constraints of the cutting procedure are clear, its practical and numerical implementation has to be treated carefully. Indeed we need simultaneously that M N → ∞ (i.e., a large number of subwindow domains) and that dist(V (N,(i,j)) , V (N,(i ,j )) ) → ∞ (i.e., the size of each subwindow domain is large).

As previously described in [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF], this is the crucial and well known compromise between variance and bias. Indeed, [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF] provide a numerical study in the very specific case of a spherical covariance model to find the best subwindow size in order to minimize the mean error for the variance estimator for three intermediate levels u. The appropriate size and the form (e.g., square, rectangular, . . . ) of the subwindow is a crucial issue in applications (see Section 4). Furthermore, the cutting of T (N ) has to guarantee that in each subwindow one can find observations of the excursion set at the chosen level u ∈ R. Then, the subwindow procedure is also implicitly related to the choice of the (intermediate, large or extremes) value of observation level u.

A Proofs

A.1 Preliminary results

In the following we prove two auxiliary consistency results for C 0,T (N ) (X, u) and C 2,T (N ) (X, u).

Proposition A.1. Let X be a Gaussian random field satisfying Assumptions (A0) and (A2).

Q (N ) i := |T (N ) i | ( C 0,T (N ) i (X, u i ) -C * 0 (X, u i ))
, for i ∈ {1, . . . , m}.

As N → +∞, Q (ii) Let T be a cube in R d , u 1 , . . . , u m be m levels in R. For any integer N > 0 and C 0,T (N ) (X, u i )

as in (10), let

S (N ) i := |T (N ) | ( C 0,T (N ) (X, u i ) -C * 0 (X, u i )), for i ∈ {1, . . . , m}.
As N → +∞, S Proof. Let i ∈ {1, . . . , m}. The ith coordinate of the considered multivariate central limit theorem can be written using the following decomposition

Q (N ) i = |T (N ) i |( C 0,T (N ) i (X, u i ) -C * 0 (X, u i )) = |T (N ) i | C /T (N ) i 0 (X, u i ) -E[C /T (N ) i 0 (X, u i )] - 1 π |T | C /T (N ) i 1 (X, u i ) -E[C /T (N ) i 1 (X, u i )] |∂T (N ) i | 1 |T (N ) i | + |T (N ) i | C /T (N ) i 2 (X, u i ) -E[C /T (N ) i 2 (X, u i )] 1 2π |∂T (N ) i | 2 1 |T (N ) i | 2 - 1 |T (N ) i | := I 0 (T (N ) i ) + I 1 (T (N ) i ) + I 2 (T (N ) i
).

For a Gaussian random field satisfying Assumptions (A0) and (A2), Theorem 1.1 in [START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF] and Theorem 2.1 in [START_REF] Müller | A central limit theorem for Lipschitz-Killing curvatures of Gaussian excursions[END_REF] give

|T (N ) i | C /T (N ) i 1 (X, u i )-E[C /T (N ) i 1 (X, u i )]
admits a Gaussian centered limit distribution and therefore I 1 (T

(N ) i ) P ----→ T R 2 0. The same discussion holds for I 2 (T (N ) i ). Let Z (N ) i = |T (N ) i |(C /T (N ) i 0 (X, u i ) -E[C /T (N ) i 0 (X, u i )]
) for i ∈ {1, . . . , m}.

One can adapt the proof of the CLT in Proposition 5 (a) in Di [START_REF] Di Bernardino | A test of Gaussianity based on the Euler characteristic of excursion sets[END_REF] to prove that for X satisfying Assumptions (A0) and (A2), X(0) ∼ N (µ, σ 2 ) and as N → +∞, (Z

(N ) 1 , . . . , Z (N ) m )
converges in distribution to a centered Gaussian vector with diagonal covariance matrix and finite elements. Finally, applying the multivariate Slutsky's theorem we get the item (i). The proof of (ii) comes down in a similar way using item (b) of Theorem 2.5 in [START_REF] Estrade | A Central Limit Theorem for the Euler characteristic of a Gaussian excursion set[END_REF].

Proposition A.2 (Asymptotic normality of ( C 2,T (N ) (X, u 1 ), . . . , C 2,T (N ) (X, u m ))). Let X be a Gaussian random field satisfying Assumptions (A0) and (A1). For a positive integer N , consider T (N ) = {N t : t ∈ T } and C 2,T (N ) (X, u) the estimator defined in (8) built on the observation

T (N ) ∩ E X (u i ),
where u 1 , . . . , u m are fixed. Then,

|T (N ) |( C 2,T (N ) (X, u 1 ) -C * 2 (X, u 1 ), . . . , C 2,T (N ) (X, u m ) -C * 2 (X, u m ))
converges in distribution to a centered Gaussian vector with m×m covariance matrix (Σ 2

C * 2 ,(u i ,u j ) ) 1≤i,j≤m given by Σ 2 C * 2 ,(u i ,u j ) = 1 2π R 2 ρ(t) 0 g (u i ,u j ) (r)dr dt ∈ (0, +∞) where g (u i ,u j ) (r) = 1 √ 1 -r 2 exp - (u i -µ) 2 -2r(u i -µ)(u j -µ) + (u j -µ) 2 2σ 2 (1 -r 2 ) .
The proof of Proposition A.2 comes down from Theorem 4 in [START_REF] Bulinski | Central limit theorems for the excursion set volumes of weakly dependent random fields[END_REF], together with [START_REF] Shashkin | Quasi-associatedness of a Gaussian system of random vectors[END_REF] ensuring Gaussian fields are quasi-associated.

Σ2

C * 0 = Σ 2 C * 0 ,(u i ,u j ) 1≤i,j≤m is defined in Proposition A.1, Σ 2 1 h = J 1 h Σ 2 s J t 1 h with Σ 2 s = Σ 2 s,(u i ,u j ) 1≤i,j≤m
defined in Proposition 2.1 and J 1 h the m × m Jacobian matrix

J 1 h = (2π) 3 2 diag (s 2 u 1 -1) s 2 u 1 exp s 2 u 1 2 , . . . , (s 2 um -1) s 2 um exp s 2 um 2 .
After computations we get AΣ 2

C * 0 A t = Σ 2 C * 0 ,(u i ,u j ) h(su i )h(su j ) 1≤i,j≤m and
BΣ 2 1 h B t = (2π) 3 C * 0 (X, u i )C * 0 (X, u j ) (s 2 u i -1)(s 2 u j -1) s 2 u i s 2 u j exp 1 2 (s 2 u i + s 2 u j ) Σ 2 s,(u i ,u j ) 1≤i,j≤m
.

Adding these expressions provides the desired formula.

Proof of Proposition 3.2

First, write

|T (N ) 1 | R T (N ) 1 (u 1 , u 2 ) -R H 0 T (N ) 2 (u 1 , u 2 ) = |T (N ) 1 | R T (N ) 1 (u 1 , u 2 ) -R H 0 (u 1 , u 2 ) (26) + |T (N ) 2 | R H 0 (u 1 , u 2 ) -R H 0 T (N ) 2 (u 1 , u 2 ) .
For the first term in (26) we use Proposition A.1 (ii) and the delta method with the function g : (x, y) → x y . We get after computations, using that R T

1 (u 1 , u 2 ) = g C 0,T 1 (X, u 2 ), C 0,T 1 (X, u 1 ) , |T (N ) 1 | g C 0,T (N ) 1 (X, u 2 ), C 0,T (N ) 1 (X, u 1 ) -g C * 0 (X, u 2 ), C * 0 (X, u 1 ) d, H 0 ----→ N →∞ N 0, σ 2 g(C * 0 ),(u 1 ,u 2 )
where

σ 2 g(C * 0 ),(u 1 ,u 2 ) = σ 2 C * 0 ,u 2 C * 0 (X, u 1 ) 2 -2 C * 0 (X, u 2 ) C * 0 (X, u 1 ) 3 Σ 2 C * 0 ,(u 1 ,u 2 ) + C * 0 (X, u 2 ) 2 C * 0 (X, u 1 ) 4 σ 2 C * 0 ,u 1 .
For the second term in (26) we use Proposition 2.1 and apply the delta method with the function h : (x, y) → x y exp{-1 2 (x 2 -y 2 )}. We get after computations, using that R H0

T (N ) 2 (u 1 , u 2 ) = h( s u 2 ,T (N ) 2 , s u 1 ,T (N ) 2 ), |T (N ) 2 | h(s u 2 , s u 1 ) -h( s u 2 ,T (N ) 2 , s u 1 ,T (N ) 2 ) d, H 0 ----→ N →∞ N 0, σ 2 h(s),(u 1 ,u 2 ) where σ 2 h(s),(u 1 ,u 2 ) = ∇h(s u 2 , s u 1 ) Σ 2 s,(u 1 ,u 2 ) ∇h(s u 2 , s u 1 ), Σ 2 s,(u 1 ,u 2 ) is defined as in Proposition 2.1, σ 2 h(s),(u1,u2) = e -(s 2 u 2 -s 2 u 1 ) σ 2 su 2 1 -s 2 u2 s u1 2 + 2Σ 2 s,(u1,u2) 1 -s 2 u2 s u1 - s u2 s 2 u1 + s u2 + σ 2 su 1 - s u2 s 2 u1 + s u2 2 .
Using dist(T

(N ) 1 , T (N ) 2 ) → ∞ as N → ∞ and σ 2 Ru 1 ,u 2 = σ 2 g(C *
Since J is a finite set and using Corollary 3.1, it immediately holds that

σ 2 su 1 , . . . , σ 2 su J P ----→ T R 2 σ 2 su 1 , . . . , σ 2 su J .
Denote by j * = argmin j∈{1,...,J} σ 2 su j

, it follows

P( j = j * ) = P j =j * { σ 2 su j ≤ σ 2 su j * } ≤ j =j * P( σ 2 su j -σ 2 su j + (σ 2 su j -σ 2 su j * ) ≤ σ 2 su j * -σ 2 su j * ) ≤ j =j * P(|σ 2 su j -σ 2 su j * | ≤ | σ 2 su j * -σ 2 su j * | + | σ 2 su j -σ 2 su j |) ≤ j =j * P( 1 2 |σ 2 su j -σ 2 su j * | ≤ | σ 2 su j * -σ 2 su j * |) + P( 1 2 |σ 2 su j -σ 2 su j * | ≤ | σ 2 su j -σ 2 su j |) ----→ T R 2 0,
where we used that J is finite, ∀j = j * , |σ 2 su j -σ 2 su j * | > 0 and Corollary 3.1. Finally, we derive that j P -→ j * , for T R 2 , together with the definition of j * and the fact that u → σ 2 su is decreasing on (-∞, µ) and increasing on (µ, ∞), the result follows.

A.3 Proof of Proposition 3.1

Preliminaries on the Itô-Wiener chaos decomposition for L 2 Let G := (X -µ)/σ be the centered and unit variance Gaussian random field associated to X, for all fixed levels u. Recall that s u = (u -µ)/σ, then L 2,T (X, u) = L 2,T (G, s u ), which is a square-integrable functional of the Gaussian field G. It admits an orthogonal decomposition into Itô-Wiener chaos in the L 2 sense (see, e.g., [START_REF] Marinucci | Stein-Malliavin approximations for nonlinear functionals of random eigenfunctions on S d[END_REF], [START_REF] Nourdin | Normal approximation with Malliavin calculus[END_REF]):

L 2 (G, s u , T ) = +∞ q=0 β q (s u ) q! T H q (G(t)) dt,
where H q is the q-th Hermite polynomial, i.e., for z ∈ R, H 0 (z) = 1 and

H q (z) = (-1) q exp{z 2 /2} d q dz q exp{-z 2 /2}, if q ≥ 1.
The above series converges in L 2 (P) and for Z ∼ N (0, 1), β q (s u ) := E 1 {Z≥su} H q (Z) . The chaotic coefficients (β q (s u )) q≥0 for L 2 (G, s u , T ) are given by: β 0 (s u ) = ψ(s u ) and

β q (s u ) = +∞ su ϕ(z) (-1) q ϕ(z) d q dz q ϕ(z)dz = (-1) q-1 d q-1 dz q-1 ϕ(s u ) = ϕ(s u )H q-1 (s u ),
for q ≥ 1, where ϕ is the probability density function of Z. Denoting H φ q the "physicist Hermite polynomials" we have H φ q (x) = 2 q/2 H q ( √ 2x) and it holds ∀x, |(2 q q! √ π) -1/2 H φ q (x)e -x 2 /2 | ≤ C ∞ /(q + 1) 1/12 (see [START_REF] Szegő | Orthogonal Polynomials[END_REF] for the the value of the constant C ∞ ). As, it holds ∀x ∈ R

H φ q (x)e -x 2 /2 = 2 q/2 H q ( √ 2x)e -x 2 /2 = e x 2 /2 2 q/2 H q ( √ 2x)e -( √ 2x) 2 /2 , it follows ∀x ∈ R |β q+1 ( √ 2x)| = |H q ( √ 2x)e -( √ 2x) 2 /2 | ≤ e -x 2 /2 q! √ πC ∞ /(q + 1) 1 12 ≤ C ∞ π 1 4 √ q! (q + 1) 1 12
.

We derive the following inequality, frequently used in the sequel

β 0 ∞ ≤ 1 and β q ∞ ≤ c β (q -1)! q 1 12 , q ≥ 1. ( 27
)
The mean of L 2 (G, s u , T ) is its projection onto the 0-th Itô-Wiener chaos

E [L 2 (G, s u , T )] = ψ(s u )|T |.
The remaining of the proof consists in controlling the order of the fourth moment of L 2 (G, s . , T ), computed using the chaos decomposition as follows

E[L 2 (G, s u 1 , T )L 2 (G, s u 2 , T )L 2 (G, s u 3 , T )L 2 (G, s u 4 , T )] = +∞ k 1 ,k 2 ,k 3 ,k 4 =0 β k 1 (s u 1 )β k 2 (s u 2 )β k 3 (s u 3 )β k 4 (s u 4 ) k 1 !k 2 !k 3 !k 4 ! × T T T T E[H k 1 (G(t 1 ))H k 2 (G(t 2 ))H k 3 (G(t 3 ))H k 4 (G(t 4 ))] dt 1 dt 2 dt 3 dt 4 .
Control of a fourth moment of L 2 (G, s . , T ) We sometimes denote T 1 = T 2 := T and T 3 =

T 4 := T and use that |T | = |T |. We write, E[L 2 (G, s u 1 , T )L 2 (G, s u 2 , T )L 2 (G, s u 3 , T )L 2 (G, s u 4 , T )] = V 0 + V 1 + V 2 + V 3 + V 4 , (28)
where V 0 : = |T | 4 β 0 (s u 1 )β 0 (s u 2 )β 0 (s u 3 )β 0 (s u 4 )

V 1 : = |T | 3 4 j=1 β 0 (s u 1 )β 0 (s u 2 )β 0 (s u 3 )β 0 (s u 4 ) β 0 (s u j )

+∞ k j =1 β k j (s u j ) k j ! T j E[H k j (G(t j ))] dt j V 2 : = |T | 2 j 1 =j 2 ∈{1,...,4} +∞ k i 1 ,k i 2 =1 i 1 ,i 2 =j 1 ,j 2
β k i 1 (s u i 1 )β k i 2 (s u i 2 )β 0 (s u j 1 )β 0 (s u j 2 ) β 0 (s u j 1 )β 0 (s u j 2 )

k i 1 !k i 2 ! T i 1 T i 2 E[H k i 1 (G(t i 1 ))H k i 2 (G(t i 2 ))] dt i 1 dt i 2 V 3 : = |T | 4 j=1 β 0 (s u j ) +∞ k i 1 ,k i 2 ,k i 3 =1 i 1 =j,i 2 =j,i 3 =j β k i 1 (s u i 1 )β k i 2 (s u i 2 )β k i 3 (s u i 3 ) k i 1 !k i 2 !k i 3 ! × T i 1 T i 2 T i 3 E[H k i 1 (G(t i 1 ))H k i 2 (G(t i 2 ))H k i 3 (G(t i 3 ))] dt i 1 dt i 2 dt i 3 V 4 : = +∞ k 1 ,k 2 ,k 3 ,k 4 =1 β k 1 (s u 1 )β k 2 (s u 2 )β k 3 (s u 3 )β k 4 (s u 4 ) k 1 !k 2 !k 3 !k 4 ! × T T T T E[H k 1 (G(
+∞ k=1 i 1 ,i 2 =j 1 ,j 2 β k (s u i 1 )β k (s u i 2 ) k! T i 1 T i 2 ρ(t i 1 -t i 2 ) k dt i 1 dt i 2 = |T | 2 j 1 =j 2 ∈I 2,1 + j 1 =j 2 ∈I 2,2
β 0 (s u j 1 )β 0 (s u j 2 )

+∞ k=1 i 1 ,i 2 =j 1 ,j 2 β k (s u i 1 )β k (s u i 2 ) k! T i 1 T i 2 ρ(t i 1 -t i 2 ) k dt i 1 dt i 2 =: V 2,1 + V 2,2 ,
where I 2,1 := {(1, 3), (1, 4), (2, 3), (2, 4)} and I 2,2 := {(1, 2), (3, 4)}. Using ( 27), (A1) together with the fact that on I 2,1 the integration is made on distinct rectangles T and T , we get

|V 2,1 | ≤ 4c 2 β |T | 2 +∞ k=1 (k -1)! k!k 1/6 T T |ρ(t -t )| k dt dt ≤ 4c 2 β |T | 2 +∞ k=1 1 k 1+1/6 T T 1 (1 + t -t ) γk dt dt ≤ 4c 2 β |T | 2 +∞ k=1 1 k 1+1/6 1 dist(T, T ) γk-2 T T 1 (1 + t -t ) 2 dt dt ≤ 4c 2 β |T | 2 +∞ k=1 1 dist(T, T ) γk-2 R 2 1 (1 + t -t ) 2 dt dt ≤ C 2,1 |T | 2 dist(T, T ) γ-2 ,
where we used Assumption (A1) and dist(T, T ) > 2, then C 2,1 is given by 8c 2 β R 2 (1+ t-t ) -2 dt dt .

For V 2,2 , using (27) and Assumption (A1), we have Control of V 3 . To control this term we rely on the diagram formula and the following decomposition

V 3 = |T | 4 j=1 β 0 (s u j ) +∞ k i 1 ,k i 2 ,k i 3 =1 i 1 =j,i 2 =j,i 3 =j β k i 1 (s u i 1 )β k i 2 (s u i 2 )β k i 3 (s u i 3 ) k i 1 !k i 2 !k i 3 ! × T i 1 T i 2 T i 3 E[H k i 1 (G(t i 1 ))H k i 2 (G(t i 2 ))H k i 3 (G(t i 3 ))] dt i 1 dt i 2 dt i 3 = |T | 4 j=1 β 0 (s u j ) +∞ k i 1 ,k i 2 ,k i 3 =1 i 1 =j,i 2 =j,i 3 =j k i 1 +k i 2 +k i 3 ≤2N T,T + +∞ k i 1 ,k i 2 ,k i 3 =1 i 1 =j,i 2 =j,i 3 =j k i 1 +k i 2 +k i 3 >2N T,T β k i 1 (s u i 1 )β k i 2 (s u i 2 )β k i 3 (s u i 3 ) k i 1 !k i 2 !k i 3 ! × T i 1 T i 2 T i 3 E[H k i 1 (G(t i 1 ))H k i 2 (G(t i 2 ))H k i 3 (G(t i 3 ))] dt i 1 dt i 2 dt i 3 = V 3,≤N T,T + V 3,>N T,T ,
for some positive integer N T,T depending on T and T and such that N T,T → ∞ as dist(T, T ) → ∞.

First, note that V 3,>N T,T = o(|T | 2 ) as dist(T, T ) → ∞: using that L 2 (G, s u , T ) ≤ |T | a.s. and that the Itô-Wiener chaos decomposition holds in the L 2 sense, we write

V 3,>N T,T ≤ C|T | 2 max u∈{u 1 ,u 2 ,u 3 ,u 4 } E q≥N T,T β q (s u ) q! T H q (G(t))dt 2 ,
which tends to 0 as the remainder of a convergent series. To control this term, we take advantage of the cutting described in Figure 3.

We now focus on V 3,≤N T,T , by symmetry we set k 4 = 0 and compute

T k 1 ,k 2 ,k 3 := E[H k 1 (G 1 )H k 2 (G 2 )H k 3 (G 3 )]
for k j ≥ 1, j ∈ {1, 2, 3} and (G 1 , G 2 , G 3 ) a standard Gaussian (see [START_REF] Taqqu | Law of the iterated logarithm for sums of non-linear functions of Gaussian variables that exhibit a long range dependence[END_REF] Definition 3.1), i.e.

a centered Gaussian vector with E[G 2 i ] = 1 and ρ i,j := E[G i G j ] such that |ρ i,j | ≤ 1, 1 ≤ i, j ≤ 3. The diagram formula (see [START_REF] Taqqu | Law of the iterated logarithm for sums of non-linear functions of Gaussian variables that exhibit a long range dependence[END_REF], Lemma 3.2) gives

T k 1 ,k 2 ,k 3 =          k 1 !k 2 !k 3 ! 2 q q! I(k 1 ,k 2 ,k 3 ) ρ i 1 ,j 1 . . . ρ iq,jq if k 1 + k 2 + k 3 = 2q, 1 ≤ k 1 , k 2 , k 3 ≤ q 0 otherwise (29)
where the set of indices I(k 1 , k 2 , k 3 ) is the set of all indices (i 1 , j 1 , . . . , i q , j q ) such that i) (i 1 , j 1 , . . . , i q , j q ) ∈ {1, 2, 3} 2q , ii) i 1 = j 1 , . . . , i q = j q , iii) there are k 1 indices 1, k 2 indices 2 and k 3 indices 3.

As ρ i,j = ρ j,i , after multiplying T k 1 ,k 2 ,k 3 by 2 q , I(k 1 , k 2 , k 3 ) simplifies in I(k 1 , k 2 , k 3 ) the set of all indices (i 1 , j 1 , . . . , i q , j q ) such that i+ii) ((i 1 , j 1 ), . . . , (i q , j q )) ∈ {(1, 2), (1, 3), (2, 3)} q , iii) there are q 1 times (1, 2), q 2 times (1, 3) and q 3 times (2, 3) with q 1 = k 1 + k 2 -q = q -k 3 , q 2 = k 1 + k 3 -q = q -k 2 and q 3 = k 2 + k 3 -q = q -k 1 .

The cardinality of I(k 1 , k 2 , k 3 ) is   q q 1 , q 2 , q 3   = q! q 1 !q 2 !q 3 ! . Set J 3 (q) = {(k 1 , k 2 , k 3 ) ∈ {1, . . . , q} 3 , k 1 +

k 2 + k 3 =
|β k 1 (s u * )β k 2 (s u * )β k 3 (s u * )| q 1 !q 2 !q 3 ! × T T T |ρ(t 1 -t 2 )| q-k 3 |ρ(t 1 -t 3 )| q-k 2 |ρ(t 2 -t 3 )| q-k 1 dt 1 dt 2 dt 3 ≤ 4c 3 β |T | N T,T q=4 (k 1 ,k 2 ,k 3 )∈J 3 (q) (k 1 -1)!(k 2 -1)!(k 3 -1)! (k 1 k 2 k 3 ) 1 12 (q -k 3 )!(q -k 2 )!(q -k 1 )! × T T T
dt 1 dt 2 dt 3 (1 + t 1 -t 2 ) γ(q-k 3 ) (1 + t 1 -t 3 ) γ(q-k 2 ) (1 + t 2 -t 3 ) γ(q-k 1 ) .

Set

A q (k 1 , k 2 , k 3 ) := (k 1 -1)!(k 2 -1)!(k 3 -1)! (k 1 k 2 k 3 ) 1 12 (q -k 1 )!(q -k 2 )!(q -k 3 )!

, k 1 + k 2 + k 3 = 2q.

Figure 1 :

 1 Figure 1: Excursion sets of 3 images of synthetic digital mammograms observed at the same level (first row) and at the same effective level (second row).

  is a Van Hove-growing sequence (VH-growing sequence, see Definition 6 in Bulinski et al. (2012)), i.e., |∂T (N ) | 1 /|T (N ) | → 0 as N → ∞, where ∂T stands for the frontier of the set T . In the sequel, we sometimes drop the dependency in N of the rectangle T to soften notation.

(Figure 4 :

 4 Figure 4: Synthetic digital mammograms study. Image 1.F from group (F) (left), 1.FG from group (FG) (center) and 1.D from group (D) (right). Image size: 251 × 251.

  (right-down) is obtained by considering u 1.F = 2133, u 1.F G = 2291, u 1.D = 2518.

Figure 5 :

 5 Figure 5: Synthetic digital mammograms study. Left: Estimated s u for several values of u for each image: group (F) in blue curves, (FG) in green curves and (D) in red ones. Right: Excursion sets for a fixed level u = 2200 (first row) and for the three adaptive levels u, such that for each u it holds that | s u | < , for = 10 -2 (second row).

Figure 6 :

 6 Figure 6: Synthetic digital mammograms study. Boxplots of the estimates C 0,T ( u) in (10), for each group of images (group (F) in left panel, (FG) in center and (D) in right one) with adaptive levels u, such that for each u it holds that | s u | < , for = 10 -2 .

  for excursion sets in Figure 5 (right, second row) with u 1.F = 2133, u 1.F G = 2291, u 1.D = 2518.

Figure 7 :

 7 Figure 7: Synthetic digital mammograms study. Estimation of s u for 1000 different values of u and couples of images: 2.F and 3.F (first panel); 1.F versus 5.D (second panel); 1.F and 3.FG (third panel). In bold marked points we represent the cases where the test (25) rejects H 0 for a significant level α = 0.2. Group (F) is displayed using blue curves, (FG) green curves and (D) red ones.

Figure 8 :

 8 Figure 8: Real digital mammograms study. Image from group fatty background tissue (F) (left), from group fatty-glandular (FG) (center) and from group dense (D) (right). Image size: 251 × 251.

Figure 9 :

 9 Figure 9: Real digital mammograms study. Estimation of s u for several levels u and a couple of images. From left to right: 1.F versus 1.FG, 1.F versus 1.D, 1.FG versus 1.D, 27.F versus 19.FG, 27.F versus 20.D, 19.FG versus 20.D. In bold marked points we represent the relative u values such that the test (25) rejects H 0 for a significant level α = 0.2. Group (F) is displayed using blue curves, (FG) green curves and (D) red ones.

Figure 10 :

 10 Figure 10 we display the number of p-values associated to intra and inter-classes analyses for 200 different values of u ∈ [100, 240] that are smaller than the significant level α = 0.2. The reference value is represented by the horizontal line at level α × 200 = 40, above which H 0 is rejected.

Figure 11 :

 11 Figure 11: Real digital mammograms study. Adaptive level u (i.e., level u such that s u is close to a given value δ, y axis) for all images (x axis), such that | s u -δ| < , for = 10 -2 and δ = 0 (first panel) δ = 0.5 (second panel), δ = 1 (third panel) and δ = 3 (fourth panel). Images from group fatty (F) are displayed by blue starts, fatty-glandular group (FG) by green crosses and dense group (D) by red points.

  (i) Let T 1 , . . . , T m be m cubes in R d such that |T 1 | = . . . = |T m | and dist(T i , T j ) > 0 for all i = j ∈ {1, . . . , m} 2 . Let u 1 , . . . , u m be m levels in R, for any integer N > 0 and C 0,T (N ) i (X, u i ) as in (10), let

  i < +∞, ∀i ∈ {1, . . . , m}.

  u i ,u j ) 1≤i,j≤m with Σ 2 C * 0 ,(u i ,u j ) < +∞, for i, j ∈ {1, . . . , m}.

  |V 2,2 | ≤ |T | 2 |β 0 (s u 1 )β 0 (s u 2 )| ∞ k=1 |β k (s u 3 )β k (s u 4 )| k! T T |ρ(t -s)| k dt ds + |β 0 (s u 3 )β 0 (s u 4 )| ∞ k=1 |β k (s u 1 )β k (s u 2 t -s ) γ ≤ C 2,2 |T | 2 , where C 2,2 = 2c 2 β R 2 (1 + t -t ) -2 dt dt k≥1 k -1-1 6 is a constant independent of T . It follows that V 2 = O(|T | 2 ).

  2q}, the set of indices where T k 1 ,k 2 ,k 3 is nonzero. Using (27), Assumption (A1) and the cardinality of I(k 1 , k 2 , k 3 ), it follows from (29) that,|V 3,≤N T,T | ≤ 4 |T | max u * ∈(u 1 ,u 2 ,u 3 ,u 4 ) |β 0 (s u * )| 4 N T,T q=4 (k 1 ,k 2 ,k 3 )∈J 3 (q)

  for instance, Proposition A.1, but conversely to Section 3.1 of[START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF], as R H 0 (u 1 , u 2 ) is now unknown, this CLT does not allow to determine a rejection level. To overcome that problem, define

Table 1 :

 1 Table 1 below, lead to accept the H 0 hypothesis. This preliminary adjustment of the image level is necessary to properly compare the excursion sets. In case only two excursion sets are available and not the whole images, this test (24) can be a preliminary prerequisite to determine whether or not it is legitimate to perform a comparison between these images. Synthetic digital mammograms study. p-values associated to the test H 0

	1.F versus 1.FG 1.F versus 1.D 1.FG versus 1.D
	0.9858	0.9511	0.9642

Table 2

 2 Finally, we provide a graphical illustration of the test in (25) for three couple of images (2.F and 3.F, first panel; 1.F versus 5.D, second panel; 1.F and 3.FG, third panel). In bold marked points we represent the cases when the test (25) rejects H 0 for at level α = 0.2. These points are drawn on the estimates s u , for 1000 values of u. Coherently, the test accepts H 0 (see Table2) in the first

	).

  Wschebor, M. (2006). Surfaces aléatoires: mesure géométrique des ensembles de niveau, volume 1147. Springer.

	Corresponding author's address: Elena Di Bernardino, Laboratoire J.A. Dieudonné, UMR
	CNRS 7351, Université Côte d'Azur, Parc Valrose, 06108 Nice, Cedex 2, France; elenadb@unice.fr

  t 1 ))H k 2 (G(t 2 ))H k 3 (G(t 3 ))H k 4 (G(t 4 ))] dt 1 dt 2 dt 3 dt 4 . It holds V 0 = ψ(s u 1 )ψ(s u 2 )ψ(s u 3 )ψ(s u 4 )|T | 4 and by orthogonality of H k 1 with H 0 = 1, ∀ k 1 ≥ 1, we get V 1 = 0.It remains to control the terms V 2 , V 3 and V 4 .

Control of V 2 . Using the fact that E[H k (G(t))H l (G(s))] = δ k,l k!ρ(t -s) k (see e.g. Equation (2.1)

in

[START_REF] Breuer | Central limit theorems for non-linear functionals of Gaussian fields[END_REF]

), we get

V 2 = |T | 2 j 1 =j 2 ∈{1,...,4}

The same atypical behavior for figure 5.F has been previously observed in[START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF], where the same data-set is studied.

),(u

,u

) + σ 2 h(s),(u 1 ,u 2 ) we get the result.

Acknowledgments:

The authors express their gratitude to two anonymous Referees for their valuable comments on this article. The authors are grateful to Anne Estrade for fruitful discussions.

The authors acknowledge Z. Li, (GE Healthcare France, department Mammography) for simulated 2D digital mammograms data-sets in Section 4.1 and the related discussion. Mammographic Image Analysis Society database can be found here http://peipa.essex.ac.uk/info/mias.html. This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002. This work has been partially supported by the project ANR MISTIC (ANR-19-CE40-0005).

Supporting information. Additional information for this article is available online.

values of u that are smaller than the significant level α = 0.2. The numbers larger than α × 1000 = 200 for which H 0 is rejected.

In the intra-class analysis only tests involving the image 5.F lead to a number of p-values slightly larger than α × 1000 = 200, for which H 0 is rejected 1 . In all the other cases, the test accepts the H 0 hypothesis (see left-side Table 2). Conversely in the inter-classes analysis all the obtained numbers of p-values are much higher than 200 (right-side table ). This means that, to properly compare excursion set of a given level u for images belonging to different classes, a previous effective level

A.2 Proofs of the obtained results

Proof of Proposition 1.1

The Gaussian kinematic formula provides the mean LK curvatures of excursion sets of X within a rectangle T (see, e.g., Theorem 15.9.5 in Adler andTaylor (2007) or Theorem 4.8.1, 4.3.1 in Adler and[START_REF] Adler | Topological complexity of smooth random functions[END_REF]), for u ∈ R and i = 0, 1, 2,

where L j (T ), j = 0, 1, 2 are defined in (5),

with ω k the Lebesgue measure of the k-dimensional unit ball (w 0 = 1, w 1 = 2 and w 2 = π), and, following Formula (3.5.2) in [START_REF] Adler | Topological complexity of smooth random functions[END_REF], the coefficients M l (X, u), l = 0, 1, 2 are obtained having an expansion in θ at order 2 of

where G(0) is a Gaussian random variable with zero mean and unit variance and ψ its tail distribution. This concludes the proof.

Proof of Proposition 2.2

Set h : s → s exp{-s 2 /2})/(2π) 3/2 and note that |T

-w) can be decomposed as follows

where s u,2 := s u,T (N )

2

. Since h is continuous, using Proposition 2.1 we get h(s

for all 1 ≤ i ≤ m. Then, by using Propositions 2.1 and A.1, the multivariate delta method and that dist(T

converges in distribution to a centered Gaussian vector with covariance matrix given by AΣ 2

, where . t denotes the matrix transposition, A and B are the diagonal matrices

Consider the decomposition,

First, we consider V 3,1 , i.e. the set {1 ≤ k 3 ≤ q/4 }. Straightforward computations and the fact

For the integral in (31), after two successive changes of variables we get, for all Q > 1 that

where the last inequality follows from an integration by part and C is a positive constant. Next, for all k 3 ≤ q/4 , it holds

where we used that k 3 ≤ q/4 and that → ( + q -k 3 -1)!(q --1)! is symmetric and maximal for = 0 (or = k 3 ). Injecting (32) and ( 33) in (31) and using that 1 ≤ k 3 ≤ q/4 , lead -for a positive constant C whose value may change from line to line-to

k is decreasing and bounded by

Second, we study the set { q/4 + 1 ≤ k 3 ≤ q}, similarly we get

, and the Stirling formula gives A q ( 2q/3 , 2q/3 , 2q/3 ) = O(q -1 2 q ). Moreover,

It follows that for dist(T, T ) > 2 4/γ , V 3,2 = o(|T |). Observe that the condition for convergence of this term depends on γ: the larger γ, the weaker is the condition on dist(T, T ) is. Gathering (34) and ( 35) in (30), we derive that

Control of V 4 . Similarly to V 3 , we rely on the diagram formula and the following decomposition

for some positive integer N T,T depending on T and T and such that N T,T → ∞ as dist(T, T ) → ∞.

Using similar arguments as for V 3,>N T,T , we derive that V 4,>N T,T = O(|T | 2 ) as dist(T, T ) → ∞.

We focus on V 4,≤N T,T and compute [START_REF] Taqqu | Law of the iterated logarithm for sums of non-linear functions of Gaussian variables that exhibit a long range dependence[END_REF] Definition 3.1) with

The diagram formula -taking into account that ρ i,j = ρ j,i , which comes down to multiplying the one given in [START_REF] Taqqu | Law of the iterated logarithm for sums of non-linear functions of Gaussian variables that exhibit a long range dependence[END_REF], Lemma 3.2, by 2 q (see also proof of Lemma 10.7 in Azaïs and Wschebor (2009)-gives

where the set of indices I(k 1 , k 2 , k 3 , k 4 ) is the set of all indices (i 1 , j 1 , . . . , i q , j q ) such that i+ii) ((i 1 , j 1 ), . . . , (i q , j q )) ∈ {(1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4)} q =: I pair q , iii) there are k 1 indices 1, k 2 indices 2, k 3 indices 3 and k 4 indices 4.

For q ∈ N, q ≥ 4 define the set of indices

Formulae ( 36) and ( 27) lead to, for a positive constant C whose value may change from line to line,

where q a for is the number of occurrences of the pair w and

Denote I pair (q) = q a , a ∈ I pair , a∈I pair q a = q, (i,j) q (i,j) 1 {i=l}∪{j=l} = k l , the latter can be rewritten as

1 12 q q (1,2) , q (1,3) , q (1,4) , q (2,3) , q (2,4) , q (3,4)

Similarly to V 3,≤N T,T , we decompose this majorant of V 4,≤N T,T in V 4,1 + V 4,2 according to {q (1,2) + q (3,4) ≤ 3q/4 } (which plays the same role as q -k 3 in V 3,≤N T,T ) and its complementary set.

First, we study the set {q (1,2) + q (3,4) ≤ 3q/4 }, where q (1,3) + q (1,4) + q (2,3) + q (2,4) ≥ q 4 , it follows

, where we used that k . ≥ 1, the cardinality of I pair (q) is bounded by q 5 and that the ratio in ( 38)

is maximal for q a = q/6 , ∀a ∈ I pair . Using that

e m provides the bound 3q 6 ! 2 q 6 ! 6 ≤ C3 q /q 2 and we get for dist(T,

Second, consider the set {q (1,2) + q (3,4) > 3q/4 },

Note that, if q (1,2) + q (3,4) > 3q/4 , it follows that max{q (1,2) , q (3,4) } > 3q/8 , then using (37) we derive that on this set (k 1 k 2 k 3 k 4 ) -1/12 ≤ Cq -1/6 , for some positive constant C. Next, we decompose V 4,2 according to the values of k = q (1,3) + q (1,4) + q (2,3) + q (2,4) ≤ q/4 and use Assumption (A1), ( 32) and (37) to get

where the constant C contains R 2 (1 + t -s ) -2 dtds. Observe that the function in q (1,2) in the last summand is symmetric with respect to (q -k)/2 and is maximal for

×

(q (1,3) + q (1,4) )!(q (2,3) + q (2,4) )!(q (1,3) + q (2,3) + q -k -2)!(q (1,4) + q (2,4) + q -k -2)! (q -k -1)!(q -k -1) .

Considering the function in q (1,3) , q (1,4) , q (2,3) , q (2,4) under the square root in the last display, it is easy -considering all possible cases-to derive that it is maximal at (k, 0, 0, 0) (the maximum is not unique). We obtain

q -1 6 q/4 k=0 1 (dist(T, T )) γk k! q (1,3) +q (1,4) + q (2,3) +q (2,4) =k k! q (1,3) !q (1,4) !q (2,3) !q (2,4) ! k!(q -2)!(q -k -2)! (q -k -1)!(q -k -1)

where we used the multinomial theorem. Therefore, we derive that

q -2-1 6 q/4 k=0 4 k (dist(T, T )) γk q k Finally, as for V 3,2 if we set η T,T := 4/(dist(T, T )) γ the function k → η k T,T q k is decreasing if q ≤ (η T,T ) -2 . Fix N T,T := (η T,T ) -2 , note that we have N T,T → ∞ as dist(T, T ) → ∞. It Gathering the control of four terms, V 0 to V 4 , leads to the announced result in Proposition 3.1.


Frédéric Michon

Loïc Forest

Elodie Collomb

Jacques Demongeot

Danielle Dhouailly

BMP2 and BMP7 play antagonistic roles in

Figure 1 .

 1 Figure 1. Michon et al.