Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields

: In the present paper we study three geometrical characteristics for the excursion sets of a two-dimensional stationary isotropic random ﬁeld. First, we show that these characteristics can be estimated without bias if the considered ﬁeld satisﬁes a kinematic formula, this is for instance the case of ﬁelds given by a function of smooth Gaussian ﬁelds or of some shot noise ﬁelds. By using the proposed estimators of these geometric characteristics, we describe some inference procedures for the estimation of the parameters of the ﬁeld. An extensive simulation study illustrates the performances of each estimator. Then, we use the Euler characteristic estimator to build a test to determine whether a given ﬁeld is Gaussian or not, when compared to various alternatives. The test is based on a sparse information, i.e. , the excursion sets for two diﬀerent levels of the ﬁeld to be tested. Finally, the proposed test is adapted to an applied case, synthesized 2D digital mammograms.


Introduction
In this paper we are interested in using simple geometrical objects to build inference and testing procedures to recover global characteristics of a twodimensional stationary isotropic random field by using a sparse information.Assume that we only have access to the excursion set of one realization of a field over one (or two) level, namely we observe a black-and-white image depicting where the field is above or below some level u (see, e.g., the second and third images of the first row of Figures 1,3,5,7 and 8).From this image, it is possible to compute several Lipschitz-Killing (LK) curvatures of the given excursion set.In our bivariate framework, there are three LK curvatures available.Roughly speaking, they are the following quantities: surface area, half perimeter and Euler characteristic.Each of these three functionals carries a different information on the geometry or topology of the excursion set.Intuitively, the surface area is related to its occupation density, the perimeter to its regularity and the Euler characteristic to its connectivity.In the present paper we have multiple objectives.Firstly, we aim at proposing a unified definition of these quantities for a large class of stationary and isotropic random fields that can be either smooth or present discontinuities.We jointly use these quantities to build unbiased estimators for the expected values of the LK curvatures of the excursion sets by solving a triangular system.In the whole paper, the use of the term unbiased means that our estimators take into account the observation bias due to the fact that the excursion set is observed on some finite (large) window and not on the whole plane.Moreover, we construct a test procedure to determine whether the excursion sets of a given field come from the realization of a Gaussian field.
The Lipschitz-Killing curvatures, depending on the scientific domain under consideration, are also called intrinsic volumes or Minkowski functionals, they can also be assimilated to genus types.Considering real-life data as a realization of a random field and analyzing them with the help of LK curvatures is an effective approach that has already been successfully applied in various disciplines in the past decades.For instance, in cosmology, the question of Gaussianity and anisotropy of the Cosmic Microwave Background (CMB) radiation has been tackled in a huge amount of publications (see [10] for a recent overview) and part of them are explicitly based on Minkowski functionals (see, e.g., [31] or [17]).The methodology of these studies is often the following: choose a specific parametric model for the CMB signal and fit the parameters so that the observed Minkowski functional matches the theoretical one.A similar methodology is also applied for analyzing the distribution of galaxies (see, for instance [18]), where again the best-fitting model is the target.Another domain where the LK curvatures are used as methodological tools is in brain imaging (see [2], Section 5, and references therein).In this context, the main purpose is to find the locations of high brain activity assuming the brain is a manifold with very complex geometry.The tools that are involved concern both the LK curvatures of high level excursions of the signal and the LK curvatures of the brain itself.
As seen in the aforementioned domains of application, the field of interest can be either a smooth (e.g.Gaussian or transformations of Gaussian fields) or a discontinuous object (e.g.shot noise fields).A shot noise field is a multidimensional extension of one dimensional models such as compound Poisson processes.It is obtained by throwing random shapes at random locations and counting how many shapes cover each point of the space with eventually random weights.Those models are very popular since explicit computations can (sometimes) be performed and simulations are rather simple.They are natural alternatives to Gaussian fields both for theoretical and applicative purposes.This underlines the need for a framework that remains robust if the nature of the field changes.
That is why we focus our study on providing tools and methods that can be applied to both Gaussian or transformations of Gaussian and shot noise fields.
The same objective can be found in [6] and [7], which have also inspired the present study.Note that all the random fields that we will consider are defined on the continuous space R 2 .There is clearly an abundant literature on parameter estimation and tests for multidimensional stationary random fields.Inference methods and tests usually rely on the estimation of the covariance function and/or on the estimation of the finite dimensional distributions of the field (see [14], [28] or [27] among others).This usually requires the observation of the entire field and/or independent copies of the field.In the present paper, we have a completely different approach based on the observation of the geometry of a single excursion set above a fixed level, considering the random field as a random surface.Let us mention some famous precursors, for instance [25], [37] or [24], who mainly worked on the modeling of sea waves considering them as realizations of a two-dimensional random field.Let us also mention two major books on this topics, [1] and [3].
In the mathematical domain, the LK curvatures are the subject of several works studying their probabilistic and statistical properties.We mainly refer to [32] and [35] for a precise definition of these geometrical and topological objects.They have been applied to the excursion sets of random fields in many situations.For instance, in [9], [6] and [5], the length of the level sets (i.e. the perimeter of the excursion sets) is taken into account, while several limit theorems are obtained for the area in [8] and [33].In [12] and [16], the Euler characteristic of excursion sets is used to test Gaussianity or isotropy.In the second reference, looking simultaneously at the intersections of excursion sets with rectangles and segments is a way to get additional information.Although the studies often focus on one of the LK curvatures, considering all the LK curvatures of the excursion sets at the same time is the purpose of recent papers such as [19], [20], [26] or [7].With the same objective -joint study of all Minkowski functionals-, but in another context -excursion sets of random fields defined on a discrete space-, one can also quote the recent papers [29] and [13].The present paper lies within the same scope since we provide three unbiased estimators, one for each LK curvature, valid for a wide class of stationary isotropic random fields.The estimator of the LK curvature corresponding to the Euler characteristic allows us to construct a test to determine if two observed excursion sets result from a Gaussian field or not.The strength of this test is that it is independent of a specific choice of the correlation function of the considered random field, in particular it does not depend on its second spectral moment.Since this function is unknown in many applications, it is an important and very useful property of the proposed test.In this sense, it can be considered as more efficient than the Gaussianity test introduced in [12].
To summarize, the originality of the present study is threefold.First, we use a very sparse information on the field, namely excursion sets above one or two levels, to recover global characteristics of the considered random field.Our second asset is to gather in the same study smooth Gaussian type fields and shot noise fields, which are archetypes of non continuous fields.The third key point of our paper is the joint use of several LK curvatures for a given excursion set.
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 Lipschitz-Killing densities for the excursion set of a two-dimensional random field.Section 2 is devoted to the study of unbiased estimators with edge correction of these LK densities from the observation of one excursion set.We examine a wide range of random fields, namely fields of Gaussian type and shot noise fields.Moreover, we show how the knowledge of the LK densities permits to recover and to infer on parameters of the considered fields.The problem of the consistency of the proposed estimators is studied, when it is possible.Furthermore, their performances are numerically analyzed.In Section 3 we build a test to detect whether a given field is Gaussian or not based on the knowledge of two excursion sets corresponding to two different levels.A variant of this test is put into practice in Section 4 on synthesized 2D digital mammograms provided by GE Healthcare France (department Mammography).Finally, Appendix section gathers the proofs of the technical results.

Lipschitz-Killing densities of excursion sets
There are three additive functionals, L j for j = 0, 1, 2, defined on subsets of Borelians in R 2 that are extensively used in the literature.Depending on the mathematical domain, they are called either intrinsic volumes, Minkowski functionals or Lipschitz-Killing curvatures.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.
All over the paper, we will denote by | • | the two-dimensional Lebesgue measure of any Borelian set in R 2 and by | • | 1 its one-dimensional Hausdorff measure.In particular, when T is a bounded rectangle in R 2 with non empty interior, where ∂T stands for the frontier of the set T .The rest of this section is dedicated to precisely define the Lipschitz-Killing (LK) curvatures of excursion sets as well as the associated LK densities.
Let T be a bounded rectangle in R 2 with non empty interior and u be a real number.For X a real-valued stationary random field defined on R 2 , we consider the excursion set within T above level u: We are interested in the mean of the three LK curvatures of the excursion set, assuming they are well defined.The case of the area L 2 is particularly easy to handle with.Actually, since X is stationary, Fubini's Theorem gives immediately This formula is valid without any further assumption on X, whereas the existence and the exact value of the two other LK curvatures is more involving.In particular it is often needed to make strong regularity assumptions on the field in order that one can consider the LK curvatures of its excursion sets.
In the following definitions, we introduce positive reach sets and curvature measures for those sets (see [35], based on the seminal work of Federer [15]).

Definition 1.1 (Positive reach sets). For a set
Then, the reach of A is defined as Intuitively, A is a positive reach set if one can roll a ball of positive radius along the exterior boundary of A keeping in touch with A. Let us remark that positive reach sets are necessarily closed sets.Moreover, convex sets in R 2 and compact C 2 submanifolds of R 2 have positive reach, see e.g., Proposition 14 in [35]).

Definition 1.2 (Curvature measures).
Let A be a positive reach set.Its curvature measures Φ i (A, •), for i = 0, 1, 2, are defined for any Borel set U ⊂ R 2 by where TC(∂A, U ) denotes the integral over U of the curvature along the positively oriented curve ∂A.
For a compact positive reach set A such that A ⊂ U , with U an open set, let us note that Φ 0 (A, U ) coincides with the Euler characteristic of A. For more details on the total curvature TC, we refer to Definition 9 and Theorem 31 in [35] or to Definition 1 and Theorem 1 in [7].
Since the curvature measures Φ i (A, •) are additive functionals with respect to A, it is natural to deal with unions of positive reach sets.Therefore, we introduce the next definition (see [35] again).

Definition 1.3 (UPR class). Let UPR be the class of locally finite unions of sets with positive reach.
In the sequel, we are interested in the curvature measures of the excursions above the level u of the random field X, namely we consider A = T ∩ E X (u).Note that the random set T ∩ E X (u) belongs to the UPR class a.s.when, for instance, the random field X is of class C 2 a.s., (actually, in this case the random set T ∩ E X (u) has positive reach, as E X (u) is a C 2 submanifold of R 2 and its intersection with the rectangle T provides compactness and positive reach property), or when E X (u) is locally given by a finite union of disks.Definition 1.4 (LK curvatures and associated densities for UPR excursion sets).Let X be a stationary random field defined on R 2 and let T be a bounded rectangle in R 2 with non empty interior.Assuming that T ∩ E X (u) is a UPR set, define the LK curvatures of the excursion set E X (u) within T by

Define the normalized LK curvatures by
and, assuming the limits exist, the associated LK densities are where lim T R 2 stands for the limit along any sequence of bounded rectangles that grows to R 2 .
As already noticed, the case i = 2 in the above definition boils down to

Kinematic formulas and inference
From now on, all the rectangles T in R 2 are bounded with non empty interior.Notation T R 2 stands for the limit along any sequence of bounded rectangles that grows to R 2 .Then, |∂T |1  |T | always goes to 0 as T R 2 .In this section, we build unbiased estimators of the LK densities in (3), for i ∈ {0, 1}.The problem is more involved than it appears.Actually, a naive approach is to consider (2), that can be easily computed from T ∩ E X (u), as an estimator of (3).But, the main difference between theses quantities comes from the boundary terms.Indeed, we can write where, loosely speaking, |T | −1 E Φ i (E X (u), T ) contains all the information on C * i (X, u) and boundary terms involving ∂T is blurring the estimation for fixed |T |.However, as |∂T | 1 /|T | → 0 when T R 2 , this bias term vanishes.Hereafter, we derive the exact relations between (2) and (3) which are valid for a wide range of fields, that are called standard in the sequel, and permit to build unbiased estimators of C * i (X, u), i ∈ {0, 1} with an edge correction.Starting from those estimators, we are able to infer on some parameters of the random field X.

Kinematic formulas
We use kinematic formulas to get unbiased estimators of C * i (X, u), i ∈ {0, 1} as it is usual to proceed in convex geometry.To this aim, we first introduce the definition of a standard random field.We borrow the adjective standard to [32] because the fields we call standard have excursion sets that are standard in the sense of Definition 9.2.1 in [32].Definition 2.1 (Standard random field).Let X be a stationary isotropic random field defined on R 2 .We say that By using (4), no edge correction is necessary for estimating C * 2 (X, u) whatever the field X is.It explains why we do not include a third constraint in the definition of a standard field.The first example of a standard field is given by the large class of fields of Gaussian type defined below.Definition 2.2 (Fields of Gaussian type).We call field of Gaussian type any random field X = F (G), where for some k ∈ N * , F : R k → R is a C 2 function and G = (G 1 , . . ., G k ) is a family of i.i.d.Gaussian random fields defined on R 2 that are C 3 , stationary, isotropic, centered, with unit variance, and such that Var G i (0) = λI 2 for some λ > 0, with G i denoting the gradient of G i for all i ∈ {1, . . ., k} and I 2 the 2 × 2 identity matrix.Proposition 2.3.Let X be a field of Gaussian type as in Definition 2.2, then X is standard at any level u ∈ R.
Proof.Following Definition 2.2, we write X = F (G) and we denote by λ the second spectral moment of the Gaussian vectorial field G.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 [1] or Theorem 4.8.1 in [2]), for u ∈ R and i = 0, 1, 2, where L j (T ), j = 0, 1, 2 are defined in (1), Lebesgue measure of the k-dimensional unit ball (w 0 = 1, w 1 = 2 and w 2 = π), and, following Formula (3.5.2) in [2], the coefficients M l (X, u), l = 0, 1, 2 are obtained having an expansion in ρ at order 2 of the probability that G(0) belongs to T ube(F, ρ) as ρ → 0 + .Namely, the expansion is given by Dividing Equation ( 5) by |T | and letting T grow to R 2 imply the existence of the LK densities and Injecting these in (5) and dividing by |T | we obtain that X is standard at level u.
Another way to get standard fields relies on convex geometry tools especially developed for Boolean fields with convex grains.Proposition 2.4.Let X be a stationary isotropic random field defined on R 2 such that for u ∈ R, its excursion set E X (u) satisfies: • a.s.T ∩ E X (u) is the union of a finite number N (T ∩ E X (u)) of compact convex sets with non empty interior, Then, X is standard at level u.
Proof.First let us note that by stationarity and isotropy of X we have that E X (u) is a stationary and isotropic random set and the assumptions imply that E X (u) is a standard random set in the sense of Definition 9.2.1 of [32].Following the notations of this book, C i (X, u, T ) = V i (T ∩ E X (u)) with V i the so-called intrinsic volume and ) with V i the intrinsic density.It holds according to Theorem 9.4.1 of [32], for i ∈ {0, 1, 2}, where, by Equation (5.5) of [32], c 0,2 0 = c 2,0 0 = c 2,1 1 = c 1,2 1 = 1 and c 1,1 0 = 2 π .We deduce that X is standard at level u.
When X is a random field that is standard at level u, we are now able to build unbiased estimators of C * i (X, u), for i = 0, 1, 2.
Proposition 2.5 (Unbiased estimators for LK densities).Let u ∈ R and let X be a stationary isotropic random field defined on R 2 that is standard at level u.
Assume we observe T ∩ E X (u) for T a rectangle in R 2 .The following quantities are unbiased estimator of We highlight that the term unbiased in Proposition 2.5 does not refer to the pixelization error that arise numerically due to the discretized representation of images.Here, the use of the term unbiased means that our estimation of the LK densities takes into account the observation bias due to the intersection of the excursion set with the observation window T .Other common edge correction methods include for instance the toroidal correction, where the edge on one side can be thought of as being wrapped around to the opposite edge.Numerically for large data sets, the minus-sampling methods are often used (see e.g., [34]; [11]), where the points near the edges are ignored for the estimation but are still considered as neighbours for interior points.However, it is statistically inefficient because it discards a substantial amount of data.
In the following two sections we analyse how the unbiased estimators in Proposition 2.5 for standard fields can be used for parameter inference purposes.Section 2.2 is devoted to fields of Gaussian type whereas Section 2.3 concerns shot-noise fields.

Gaussian field
We start by considering the simplest case of Gaussian type field, i.e. with notation of Definition 2.2, k = 1 and F : x → x, then X = G 1 := G.To compute the LK densities, we use expansion (6) to write out with ψ being the Gaussian tail distribution with zero mean and unit variance.Hence, one easily gets where λ denotes the second spectral moment of G as in Definition 2.2.
Parameter estimation.Suppose we observe T ∩E G (u) for a given level u and a given bounded rectangle T in R 2 .Then we can use the unbiased estimators of Proposition 2.5 to propose a consistent estimator for the second spectral moment of the Gaussian field G.Each of the quantities C i (G, u, T ) satisfies a Central Limit Theorem (CLT) with normalizing term equals to |T | 1/2 (see, e.g., [20] or [26]).Concerning various levels and/or various disjoint domains, a joint CLT for the Euler characteristic is proved in [12].To get the asymptotic normality of the estimator of λ, we moreover assume that the Gaussian field G satisfies the following condition.
(A) For any fixed x in R 2 , the covariance matrix of the random vector (G(x), G (x), G (x)) has full rank and r, the covariance function of G, is such that, Proposition 2.6 (CLT for the spectral moment estimator).Assume that G satisfies condition (A).Consider C 0,T (G, u) the estimator defined in (8) built on the observation T ∩ E G (u), u = 0 being fixed.By using (11), we define the estimator of λ as Then, it holds that where d − → stands for the convergence in distribution.
Remark that for u = 0, the Euler density in (11) always vanishes, independently of λ.The estimator is therefore not defined for u = 0 in Proposition 2.6.
Theorem 4 in [12] leads to I 0 (T ) ) with a finite asymptotic variance V (u) given by Equation ( 9) in [12].Furthermore, the CLT in [20] gives that 1 (G, u)] admits a Gaussian centered limit distribution and therefore I 1 (T ) Numerical illustration.Figure 1 displays a realization of a Gaussian random field, two excursion sets and an illustration for the performance of the three estimators x 2 and second spectral moment λ = 2κ 2 .Let us remark that the parameter κ is chosen such that |T |κ 2 remains bounded, which explains small values of κ and λ in Figure 1.The quantities 2 (G, u) are computed with the Matlab functions bweuler, bwperim and bwarea, respectively.When it is required to specify the connectivity, we average between the 4th and the 8th connectivity.Since C * 1 is defined as the average half perimeter, we divide by 2 the output derived from bwperim.
From a numerical point of view, bweuler and bwarea functions seem very precise contrary to the bwperim function which performs less well.It was expected due to the pixelisation effect.This behavior will be observed also in other random field cases (see, for instance, Figures 7 and 8 in some shot-noise cases).Figure 1 (center) illustrates that C /T 1 (G, u) (green dashed line) does not well approximate C * 1 (G, u) (blue plain line), especially for small levels u and that the correction induced by (9) (red stars) improves the approximation.In Figure 1 (left), we provide an analogous bias correction for the Euler characteristic by using (8).However in this case, the discrepancy is less evident than in the perimeter case.
In Figure 2 (left), we observe the unbiased estimator λ T (u) for different values of u.The asymptotic variance Σ(u) in Proposition 2.6 is empirically estimated on M = 100 sample simulations (Figure 2, right).This allows us to identify some choices of levels u where the variance is minimum.Furthermore, we remark that for small and large values of u, less statistics is available than for intermediate values of |u|.

Chi-square field
We consider a chi-square random field with k degrees of freedom Z k defined with with notations of Definition 2.2.In order to deal with a random field with zero mean and unit variance, we also introduce Z k as Hence, Z k is C 3 (R 2 ), stationary, isotropic, centered, with unit variance and Var Z k (0) = 2λI 2 .
We assume that Z k is observed on a rectangle T ⊂ R 2 through its excursion set above a fixed level u; ).Let us remark that the above excursion is a proper subset only for a level u such that u > − k/2.Moreover, the LK densities satisfy The C * j (Z k , •) can be computed using (7) in the framework of Gaussian type fields described previously.Indeed, we have (6) has to be written with where χ 2 k stands for a chi-square random variable with k degrees of freedom.We recover the following formulas that can be found, for instance, in Theorem 15.10.1 in [1] or in Section 5.2.1 in [2]: These formulas were originally established in [36].Figure 3 displays a realization of a normalized chi-square random field Z k with k = 2 degrees of freedom, two excursion sets and an illustration of the performance of the three estimators Parameter estimation.Suppose we observe T ∩ E Z (u), the excursion of a centered chi-square field Z with unit variance, unknown degree of freedom K and unknown Gaussian second spectral moment λ.Using that for a centered chisquare field with unit variance and k degrees of freedom C * 2 ( Z, u) only depends on k, we propose the following estimator of K: with C 2,T as in (10) and k max a large positive integer.This estimator can be plugged in C * 0 ( Z, u) to derive an estimator of the Gaussian second spectral moment of Z: with C 0,T given by (8), K(u) as in ( 14) and where M 2, Z (k, u) does not depend on λ and is derived from ( 13) An illustration of this inference procedure is provided in Figure 4 for a centered chi-square field with K = 2 degrees of freedom.

Student field
Following Definition 2.2, let k ≥ 3 be an integer and consider F : Student random field with k degrees of freedom.Assumption k ≥ 3 ensures that all the G i (t), for i = 1, . . ., k, cannot vanish at some point t ∈ R 2 with positive probability.Then, . Furthermore, T k and −T k have the same distribution.In order to deal with a centered and unit variance field, we introduce We again assume that T k is observed on a rectangle T ⊂ R 2 through its excursion set above a fixed level u, i.e., ).Let us remark that the symmetry property of the distribution of T k allows us to only consider levels u that are non-negative.Moreover and the LK densities are given in the following result.
Proposition 2.7 (Student LK densities).Let k ≥ 3 and T k be a centered Student random field with unit variance and k degrees of freedom.Then, it holds that The proof of Proposition 2.7 is postponed to Appendix (see Section A.1) as well as the proof of the following Lemma.
where λ is the spectral moment of the underlying Gaussian fields.Note that  Parameter estimation.Suppose we observe T ∩E T (u) the excursion at level u of T a centered Student field with unit variance, unknown degree of freedom K and unknown second spectral moment λ of the underlying Gaussian fields.As in Section 2.2.2, using that for a centered Student field with unit variance with k degrees of freedom C * 2 ( T k , u) only depends on k, we propose the following estimator of K: with C 2,T as in (10) and C * 2 ( T k , u) given by Proposition 2.7.Furthermore, the Gaussian second spectral moment of the considered Student field can be estimated by with C 0,T as in (8), K(u) as in (17) and where M 2, T k (k, u) does not depend on λ and is given by An illustration of the performance of this inference procedure is given in Figure 6 for a centered Student field with K = 4 degrees of freedom.Moreover by using Lemma 2.8 and the two previous estimators one can derive an estimator for λ Stu (k).

Shot-noise field
In this section, we consider a shot-noise field S Φ : R 2 → R defined as in [6] and [7].It is prescribed by where Φ is a stationary Poisson point process on R 2 × R * × R + with intensity measure ν Leb R 2 ⊗ dF B ⊗ dF R where ν > 0, F B is a probability measure on R * , F R is a probability measure on R + and D is the unit disc in R 2 .Under the condition with B and R two random variables with respective distributions F B and F R , the field S Φ is well-defined, stationary and isotropic.It is moreover integrable and, for any , ].We also introduce p = 2πE(R).
In the case B = 1 a.s., the marginals of S Φ are Poisson distributed with parameter νā.In the case where B is uniformly distributed in {−1, +1}, the field S Φ is symmetric and its marginals have a Skellam distribution with both parameters equal to νā/2 (see [30]), that is to say that they coincide with the difference of two independent Poisson random variables with parameter νā/2.In that case, the marginal distribution of S Φ is given by k ∈ Z → e −νā I |k| (νā), where I .denotes the modified Bessel function of the first kind, defined as follows: m!Γ(m+n+1) , for any n ∈ Z + and any x ∈ R. When B is positive a.s., we prove in the following result that the considered shot-noise field S Φ is standard as in Definition 2.1.Proposition 2.9.Let R max > 0. If B > 0 and R ≤ R max a.s., the random field S Φ is standard at any level u ∈ R.
Proof.Since B > 0 a.s., T ∩ E SΦ (u) = T , for any u ≤ 0. Hence we may assume that u > 0. Since R ≤ R max a.s. the number of discs intercepting T is bounded by N Φ (T ) := #{x i ; d(x i , T ) ≤ R max } that is a Poisson random variable of finite intensity.It follows that T ∩ E SΦ (u) must be given by the union of a finite number N (T ∩ E SΦ (u)) ≤ N Φ (T ) of some finite intersections of discs that are convex bodies.Since N Φ (T ) is a Poisson random variable, we clearly have E(2 N (T ∩E S Φ (u)) ) < +∞ so that Proposition 2.4 allows us to conclude that S Φ is standard at level u.
If B changes sign, excursion sets T ∩E SΦ (u) are no more a.s.closed, nor locally convex, nor necessarily of positive reach, however they are elementary sets in the sense of Definition 1 in [7], with boundaries included in the discontinuity set of S Φ so that one can compute length and total curvature.Therefore, we still define Φ i and C i as in Definitions 1.2 and 1.4.We state the following lemma, in the same vein as Theorem 6 in [7].
Lemma 2.10.Let B be a discrete random variable with values in Z\{0}.Then, for any level u, the LK densities C * 1 (S Φ , u) and C * 0 (S Φ , u) are well defined and their Fourier transforms are given, for t = 0, by where ϕ Z (t) stands for the characteristic function E[e itZ ] of a random variable Z, B + = max(B, 0) and B − = min(B, 0).
The proof of Lemma 2.10 is established in the Appendix Section A.3.Turning back to the special cases where B = 1, a.s. or B uniformly distributed in {−1, +1} and following the same inverse Fourier procedure as Theorem 6 in [7], we obtain the explicit formulas below.
where • denotes the greatest integer less than u.
Figure 7 (resp.Figure 8) displays two excursion sets and an illustration of the shape of C * i (S Φ , u), for i ∈ {0, 1, 2} when B = 1, a.s.(resp.when B is uniformly distributed in {−1, +1}).Remark that if B is uniformly distributed in {−1, +1} we do not know if S Φ is a standard random field in the sense of Definition 2.1.Therefore in Figure 8, we do not apply the unbiased formulas of Proposition 2.5 and we use C Let us quote general results of [21] that can be used to obtain CLT for the volume, the perimeter and the Euler characteristic of shot noise excursion sets (see especially Proposition 4.1 in [21]).In particular, in our setting, by Theorem 2.4 of [21], one can get for R ≤ R max and A similar question is the purpose of [19] where a joint CLT is established for all the intrinsic volumes of a Boolean model.Note that the occupied domain of a Boolean model is nothing but the excursion set above level 1 of a shot-noise field with B = 1, a.s.In our framework, we provide a non-asymptotic result to control the variance on a finite domain T and prove consistency of C i,T (S Φ , u), for i = 0, 1, 2. Theorem 2.12.If B takes values in Z \ {0} and R ≤ R max a.s., there exists a constant C, depending on ν, p, a and R max only, such that, for i = 0, 1, 2,

We recover that choosing T such that |T | → +∞ and |∂T |1
|T | → 0, leads to the consistency of C i,T (S Φ , u).The proof of Theorem 2.12 is postponed to Section A.4.To establish this result we first derive variance bounds for the weak versions of the perimeter and the total curvature (see Proposition A.1, Section A.2) which present interest on its own and is related to the framework of [21].
Parameter estimation.Suppose we observe T ∩ E SΦ (u) for S Φ a shot-noise with B = 1, a.s. and for u in R + \ Z + with unknown parameters ν, ā and p. Define the estimator of νā, for some positive A, ν a(u) = arg min with C 2,T as in (10) and νā → C * 2 (S Φ , u) given by Proposition 2.11.In the same spirit, an estimator of ν p is obtained using Proposition 2.11 with C 1,T as in (9) and ν a(u) as in (20).Then, we define the following estimator for ν where C 0,T (S Φ , u) is as in Equation ( 8), ν a(u) as in (20) and ν p(u) as in (21).Finally, the obtained ν(u) can be used to isolate a(u) and p(u) in Equations ( 20) and ( 21), respectively.An illustration of this inference procedure is provided in Figure 9 (first row).
As mentioned earlier, the unsatisfactory quality of the Matlab function bwperim for the estimation of C * 1 (S Φ , u) strongly impacts the final performance of the estimation of ν p(u) (see Figure 9, center panel).For this reason, we also present here the simplified case where R is constant.In this particular case νa = νp 2 4π .Then, we estimate ν a(u) as in (20) and ν(u) by using the simplified version of (22), i.e., An illustration of the performance of these inference procedures is provided in Figure 9 (second row).

Using LK densities for testing Gaussianity
Let X be a stationary, isotropic, centered with unit variance random field.We observe an excursion of X above some levels u on the domain T , namely T ∩ E X (u) and we aim at testing whether X is Gaussian.

Testing against a Student field
In this section we are interested in testing H0 : X is Gaussian versus where under H0 and H1 the fields are defined through Definition 2.2, where F : Under H0 and H1 the field X is symmetric and then centered.Suppose that for 1 ≤ u 1 < u 2 we observe T ∩ E X (u), for u ∈ {u 1 , u 2 }.Notice that the fact that under H0 and H1 the fields have unit variance enables the comparison of the LK densities of their excursions at the same meaningful levels.Assuming u 1 and u 2 are positive and greater than 1 is not necessary.Indeed the test can be easily modified for negative values of u 1 or u 2 by using the symmetry of X under H0 and H1.Furthermore, under H1 we do not assume that k is known.Unsurprisingly, the performance of the test will depend on the unknown k.Moreover, we underline that under H0 and H1 we do not impose any constraint on the shape of the covariance function nor on the spectral moment other than (A).In particular, the spectral moments of X under H0 and H1 can be different.
A test statistic.Our test statistics is built from the previously studied LK densities.As the perimeter C /T 1 is hard to evaluate in practice (see the earlier comment on the Matlab function bwperim), we do not use it to build a test.Moreover, Gaussian and Student distributions differ in the tails, which would lead to consider large values for u 2 .But for large levels we do not observe many excursion sets which deteriorates the estimation of C * 2 .Therefore to work with intermediate values of u 1 and u 2 , we consider C * 0 to build our test statistic: it is well evaluated in practice at these levels.Finally, in order to get statistic free in λ, we take the ratio of C * 0 between to different levels u 1 and u 2 .This is why we need to observe two distinct excursion sets.
Without loss of generality assume that u 2 = γ u 1 , for some γ > 1.For the following empirically accessible ratios, we derive from Equation (11) and from Proposition 2.7 that Quantities v(γ, ∞) and v(γ, k) in ( 24) and ( 25), for values of u 1 and γ not too large, are quite different.Furthermore, for all k it holds that v(γ, k) > v(γ, ∞) and k → v(γ, k) is decreasing.We build a non symmetric test where we reject H0 whenever the associated empirical ratio is too large compared to its expected behavior under H0.It is important to observe that v(γ, ∞) does not depend on the spectral moment of X nor on its covariance function.The choice of u 1 ≥ 1 and γ > 1 are left to the practitioner, or might be imposed by the data-set.We expect that the larger γ is and the smaller k is the better the test will be.However, in practice γ should not be too large as we need to estimate C * 0 (X, u 2 ), which requires to have sufficiently enough excursions above the level u 2 = γ u 1 .
Let T 1 and T 2 be two rectangles in R 2 such that dist(T 1 , T 2 ) > 0 and For any positive integer N , we define T with C 0,T defined in (8).
Test with asymptotic level α.We can establish the asymptotic normality of the statistics R γ,N .
Proof.As the sets T (N ) 1 and are such that dist(T ) → +∞ and | → +∞, from Proposition 5a in [12] it holds that where u → V (u) is defined in Equation ( 9) of [12].Moreover, one can prove that the limit of is non degenerate (see [20] or [26]).Then, the same decomposition as in the proof of Proposition 2.6 leads to

Testing against a power of a Gaussian field
Consider the alternative where G is a centered Gaussian random field with unit variance and σ 2 (η) := 2 1−η Γ(2η) Γ(η) , such that the field X has unit variance under H 1 (η).It is centered by construction and it is unnecessary to impose any assumption on the covariance function of G.It is straightforward to get C * 0 (X, u) = C * 0 (G, sign(u) |u| 1 η σ(η)) and the statistics in (25) for these alternatives H1 becomes, for u 2 = γu 1 and γ > 1, Note that under H0, η = 1 and v(γ, ∞) previously defined in (24) coincides with v(γ, 1).For γ > 1, the quantity in ( 29) is either smaller or larger than (24) depending on η.Symmetrizing the previous test, it follows that to test H0 : η = 1 against H1(η) : η = 1 we may consider the test with asymptotic level α In this case we have a better understanding of the behavior of the test statistic under H1.Under (A), applying (27) and with f η (u) := sign(u)|u| 1/η σ(η).If we suppose that X = sign(G) |G| η σ −1 (η) with η = 1 it holds that where from (30) we get I 1 bounded from below, which is not easy to establish, we would have P H1 ( φ , which would ensure the consistency of the test φ T (N ) .

Illustration on 2D digital mammograms
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 whole model is formulated as a spatial marked point process.The contour of each ellipsoid is blurred to render the model more realistic (for details see [22], Section 2.2 and Figure 1).Finally, considered mammograms images were simulated by x-ray projection.Evaluation provided in [22] has shown that simulated mammograms and digital breast tomosynthesis images are visually similar, according to medical experts.We can remark the discrepancy with the Gaussianity hypothesis for all images of group (C) and for all considered significant α-levels.The same consideration holds true for the last three images of group (A).Conversely, some simulated 2D digital mammograms seem to be not so far from the Gaussianity in term of the studied R γ,N ratio.Roughly speaking, images of group (B) are closer to Gaussinity than group (A), which is itself closer to Gaussianity than (C).This might be viewed in parallel with the intensities considered for each group: ν C < ν A < ν B .This chosen parameter setting may explain why our test rejects Gaussianity more easily in group C than in group B. The interested reader is referred for instance to the first four images of group (B); a very different behavior is realized by the last image of this group where the Gaussianity hypothesis is rejected for almost all considered levels u 2 .Finally, the robustness with respect to the chosen significant α-levels can be observed in Table 1.

Conclusions and discussion
We have presented new statistical tools for inferring parameters and testing Gaussianity when only a sparse observation of a 2D random field is available, namely only the excursion set(s) above one or two level(s) within a large window.These tools are based on the three LK curvatures of the excursion sets, which are loosely given by the Euler characteristic, the half perimeter and the area.
The idea of considering the statistical characteristics of the excursion sets has been originally developed for one-dimensional processes with the powerful theory of crossings.Let us comment how our two-dimensional results can be adapted to dimension one.First, we recall that only two LK curvatures are available for Borelian subsets of R: the one-dimensional Lebesgue measure (i.e. the length) and the Euler characteristic C 0 (i.e. the number of connected components).A similar statement as Proposition 2.5 provides an estimator of C * 0 (X, u) for any standard stationary process X.
In particular, when X is a Gaussian process defined on R, an adapted version of Proposition 2.6 allows to build a consistent estimator of the second spectral moment of X, which is asymptotically normal with known asymptotic variance.Indeed, the latter variance can be computed from the asymptotic variance V (u) appearing in [12] with an explicit formula in dimension one (see Proposition 8 in [12]).Let us again insist on the fact that the consistent estimator of the second spectral moment that we propose only relies on the observation of a single excursion of X.It would be interesting to compare it with more usual estimators based for instance on the estimation of the covariance function of X, which requires the observation of the whole process, or with the naive empirical estimator based on the number of crossings at different levels that is described in [23] and used in [12], Section 4.1.2.
In the same vein, for a chi-square process defined on R, the finiteness of the asymptotic variance of the Euler characteristic of excursion sets is proved in Proposition 7 in [12].Using this result, one could expect to obtain the consistency of the estimator K(u) of the degrees of freedom and the estimator λ T, K(u) (u) of the Gaussian second spectral moment, which are obtained from Equations ( 14) and ( 15) adapted to the univariate framework.
Turning back to the bivariate framework, let us discuss about the potential improvement provided by the joint observations of several LK densities.In Section 4, in order to analyse the simulated 2D digital mammograms data-set, we have used C * 0 (•, u).Let us now investigate the use of C * 2 (•, u) and the associated unbiased estimator in Proposition 2.5.In Figure 12 we display C 2,T (•, u) for each image in the 3 groups compared with the tail distribution with zero mean and unit variance ψ(u) (see Equation (11)).In order to improve the readability for extreme values of u, in Figure 12 we chose a logarithmic scale for the y−axis.(11), i.e., the Gaussian tail distribution with zero mean and unit variance.We take here a logarithmic scale for the y−axis.
A clear assessment is that the two LK densities C * 0 and C * 2 bring two different points of view on the 15 images.For instance, for extreme levels u, say above level 2.5, the C 2,T (•, u) of the mammograms clearly have a different behavior than the Gaussian one (see Figure 12).Moreover, the difference is greater for images of type (C) than for images of type (B), which is already greater than images of type (A).A similar behavior was not observed with the C 0,T (•, u).This point clearly deserves to be studied.As it involves larger levels, it is certainly related to extreme values theory.
Following the idea of taking advantage of the joint observation of several LK densities, our study can be continued with the development of efficient numerical tools adapted to different sorts of medical images.For instance, considering 2D x-ray bone images in order to detect osteoporosis, one can expect getting information on the bone density through C * 2 and on the bone connectivity through C * 0 .Both quantities really affect the mechanical bone resistance and should be both evaluated with the target of reducing the number of misclassified patients.Another direction that merits to be explored in the use of 3D images that are now quite common for medical diagnoses.Hence, one should improve our methods by including now four LK curvatures that are directly linked to the bone tissue itself instead of projected images.
In this article, we have not fully used the joint estimation of two (nor three) LK densities.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 [19] or in the case of pixelated images in [29] and binary images in [13].Those results do not apply in our context.If it exists, such a theorem would enable to build asymptotic confident regions that take into account the whole information of LK densities of the observed excursion set.
Finally, in the present study, we became aware that the Matlab function bwperim seems to not perform very well on the excursion sets of the twodimensional fields we considered, yielding an overestimation of C * 1 .The observed pixelization errors arise numerically due to the discretized representation of the smooth level set from a pixelated image, we mentioned this drawback in Section 2. It should be interesting to link our study with the numerous literature on pixelization effects.The already mentioned papers [29] and [13] should provide good tools for such an approach.

A.1. Proof of Proposition 2.7
Computation of the LK densities for T k .We use (7) after observing that Denote by F k the cumulative chi distribution (the square root of a chi-square random variable) with k degrees of freedom.Using (31) and u > 0, we have that where we have neglected the circular part of T ube(F, ρ).Then, using The latter formula is established for u > 0. The case u < 0 is derived using symmetry arguments.If u = 0, F −1 ([0, ∞)) = R + × R k and Equation ( 6) becomes Then, we derive the following formulas, valid for k ≥ 3, These quantities can be simplified using the following result E[|N (0, 1) . It follows from simple computations that Then, it is straightforward to derive the Student LK densities from (7).The above quantities do not depend on the parameter dimension of the Student field T k , which is equal to 2 in our bivariate context.
Computation of the second spectral moment of T k (Formula ( 16)).The Student field is isotropic, the second spectral moment is where for t = (T where 2 is a chi-square random field with k degrees of freedom independent of the Gaussian field G 0 .We have which is centered, as G 0 and ∂ 1 G 0 are independent, centered and independent of (Z k , ∂ 1 Z k ).Therefore, we get where we used that G 0 and ∂ 1 G 0 are independent, centered, with unit variance and independent of ( To complete the computation of λ Stu , we first have that, for k ≥ 3, E Second, using the definition of a chi-square distribution, it holds that . Taking expectation and using that ∂ 1 G i is independent of (∂ 1 G j , G j ) j =i , G i and centered, together with the definition of λ, we derive that where Z k−1 is a chi-square random field with k − 1 degrees of freedom and independent of G 1 .Then, where we made the change of variable z = x 2 y.We recognize the (k − 5)th moment of a N (0, (1 + y) −1 ) random variable in the last line.Let μ( ) := E[N (0, 1) ] with ∈ N, it holds that where B(x, y) denotes the beta function.It simplifies in Finally, for k ≥ 5,

A.2. An auxiliary result
In the following we state and prove an auxiliary result which is crucial in the derivation of Proposition 2.11 and Theorem 2.12.In the sequel we will make extensive use of notations and definitions introduced in [7], some of them will be redefined, for the others the reader can refer to [7].As in Section 2.3, we consider a shot-noise field S Φ prescribed by where Φ is a stationary Poisson point process on where τ x y = x + y.By Proposition 2 and Theorem 4 of [7], for H a primitive of h and H 1 the one dimensional Hausdorff measure (recall that H 1  where we use the fact that H 1 (τ x ∂D r ∩ τ x ∂D r ) = 0 a.s.for all different points of Φ and denote b + = max(b, 0) and b − = min(b, 0).Besides, we have r the curvature of the curve ∂D r at point z ∈ ∂D r , and with dist the geodesic distance on the circle.We can now state the following result that has interest on its own: it allows to obtain uniform bounds for the integral functionals, these bounds being usually hard to obtain when considering an excursion set at a fixed level.
First consider R SΦ (h, U ), we prove similarly

A.3. Proof of Lemma 2.10
From now on we assume that B is a discrete random variable with values in Z * .Then, the random field S Φ has values in Z, and it follows that for all u ∈ R, E SΦ (u) = E SΦ ( u ).For u ∈ Z, choosing h a continuous function with compact support in (0, 1) and such that R h(t)dt = 1 one has As a consequence C * j (S Φ , u) exists for j ∈ {0, 1} and is equal to Proof.Let us remark that the boundary contributions (see Figure 2 of [7]) may be controlled using the fact that a.s.

Fig 1 .
Fig 1. Gaussian random field as in Section 2.2.1 with covariance r(x) = e −κ 2 x 2 , for κ = 100/2 10 in a domain of size 2 10 × 2 10 pixels.First row: A realization of a Gaussian random field (left) and two excursion sets for u = 0 (center) and u = 1 (right).Second row: C 0,T (G, u) (left), C 1,T (G, u) (center) and C 2,T (G, u) (right) as a function of the level u.We display the averaged values on M = 100 sample simulations (red stars) and the associated empirical intervals (vertical red lines).Theoretical u → C * 0 (G, u), C * 1 (G, u) and C * 2 (G, u) in (11) are drawn in blue lines.We also present u → C /T 0 (G, u) and C /T 1 (G, u) in green dashed lines.These samples have been obtained with Matlab using circulant embedding matrix.

Fig 2 .
Fig 2. Gaussian random field with covariance r(x) = e −κ 2 x 2 , for κ = 100/2 10 in a domain of size 2 10 × 2 10 pixels.Estimate λ T (u) with associated confidence intervals for M = 100 sample simulations as prescribed by the CLT in (12), for different values of u (left).Theoretical value λ = 0.019 is represented by the horizontal line.The empirically estimated variance Σ(u) is displayed for different values of u in the right panel.

Fig 3 .
Fig 3. Chi-square field as in Section 2.2.2 with 2 degrees of freedom and λ = 0.019 in a domain of size 2 10 ×2 10 pixels.First row: A realization of a normalized chi-square random field (left) and two excursion sets for u = 0 (center) and u = 1 (right).Second row: C 0,T ( Z k , u) (left), C 1,T ( Z k , u) (center) and C 2,T ( Z k , u) (right) as a function of the level u.We display the averaged values on M = 100 sample simulations (red stars) and the associated empirical intervals (vertical red lines).Theoretical u → C * 0 ( Z k , u), C * 1 ( Z k , u) and C * 2 ( Z k , u) are draw in blue lines.

Lemma 2 . 8 .
If k ≥ 5, the second spectral moment λ Stu (k) of T k is finite and

Figure 5
Figure 5 displays a realization of a normalized Student random field T k with k = 4 degrees of freedom, two excursion sets and an illustration for the performance of the three estimators C 0,T ( T k , u), C 1,T ( T k , u) and C 2,T ( T k , u).

Fig 5 .
Fig 5. Student field as in Section 2.2.3 with 4 degrees of freedom and λ = 0.019 in a domain of size 2 10 × 2 10 pixels.First row: A realization of a normalized Student random field (left) and two excursion sets for u = 0 (center) and u = 1 (right).Second row: C 0,T ( T k , u) (left), C 1,T ( T k , u) (center) and C 2,T ( T k , u) (right) as a function of the level u.We display the averaged values on M = 100 sample simulations (red stars) and the associated empirical intervals (vertical red lines).Theoretical u → C * 0 ( T k , u), C * 1 ( T k , u) and C * 2 ( T k , u) are draw in blue lines.

Fig 12 .
Fig 12.C 2,T (•, u) in (10) for the 3 groups for different values u.The full black line represents ψ(u) in Equation(11), i.e., the Gaussian tail distribution with zero mean and unit variance.We take here a logarithmic scale for the y−axis.

2
and denote g m (x) = b1 rD (x) with D the unit disk.The function g m is an elementary function on R 2 (see Definition 3 in [7]) with discontinuity set given by S gm = R gm = ∂D r , where ∂D r is the circle of radius r.According to Theorem 4 of [7], a.s.for any U bounded open set, the shot noise field S Φ is an elementary function on U with discontinuity set given by U ∩ S SΦ with

Fig 14 .
Fig 14.Illustration with an elementary function given by f = 1 Dr 1 + 1 Dr 2 with 0 < r 1 < r 2 .Left: its discontinuity set S f = R f ∪ I f with I f given by the red crosses.Center and Right: E f (1) and E f (2): the boundaries are regular except at corner points given by I f .

Proposition A. 1 (
Moments of LP SΦ and LTC SΦ ).Assume that E(R 2 ) < +∞ and E(|B|) < +∞ and recall that p = 2πE(R) and a = πE(R 2 ).Then, for any U bounded open set and any continuous function h with primitive H, one has

Table 1
We consider 15 simulated 2D digital images generated by this texture model.The images were kindly provided by GE Healthcare France, department Mammography.From a clinical point of view, radiologists use the Breast Imaging Reporting and Data System (or BI-RADS) to classify breast density into four categories.They go from almost all fatty tissue to extremely dense tissue with very little fat.In this latter category, it can be hard to see small tumors in or around the dense tissue.The images we studied belong to first three morphologic situation groups: Number of p−values associated to the 1000 different values of u 2 ∈ [−3, 3] that are smaller than the significant α-levels.In bold text the numbers larger than α × 1000 for which H0 is rejected.
Description of data-set.
[7]robability measure on R + and D is the unit disc in R 2 , well defined under the condition E[|B|] E[R 2 ] < +∞, with B and R two random variables with respective distributions F B and F R .Since excursion sets of S Φ are not always UPR, we need to introduce the level perimeter LP and the level total curvature LTC integrals as introduced in[7].Using this setting, the excursions E SΦ (u) are elementary sets (see Definition 1[7]).Then, for a bounded open set U and h a continuous bounded function, one can define LP ). Var (LTCS Φ (h, U )) ≤ CH |U | (Rmaxν) 2 + ν + (νp) 2 (Rmaxνp) 2 + νa + 1 ,for a constant C H that only depends on H ∞ and such that C H ≤ 2 7 π H2 iv).