The effect of discretization on the mean geometry of a 2D random field

. — The study of the geometry of excursion sets of 2D random ﬁelds is a question of interest from both the theoretical and the applied viewpoints. In this paper we are interested in the relationship between the perimeter (resp. the total curvature, related to the Euler characteristic by Gauss–Bonnet Theorem) of the excursion sets of a function and the ones of its discretization. Our approach is a weak framework in which we consider the functions that map the level of the excursion set to the perimeter (resp. the total curvature) of the excursion set. We will be also interested in a stochastic framework in which the sets are the excursion sets of 2D random ﬁelds. We show in particular that, under some stationarity and isotropy conditions on the random ﬁeld, in expectation, the perimeter is always biased (with a 4 /π factor), whereas the total curvature is not. We illustrate all our results on diﬀerent examples of random ﬁelds. Résumé. — L’étude de la géométrie des ensembles d’excursion des champs aléatoires 2D est une question importante tant d’un point de vue théorique qu’appliqué. Dans cet article nous nous intéressons à la relation qu’il existe entre le périmètre (resp. la courbure totale, liée à la caractéristique d’Euler par le théorème de Gauss–Bonnet) des ensembles d’excursion d’une fonction et de sa discrétisée. Nous utilisons une formulation faible de cette quantité vue comme une fonction qui à un niveau lui associe le périmètre (resp. la courbure totale) de l’excursion correspondante. Nous nous intéressons également à un cadre stochastique où les fonctions sont remplacées par des champs aléatoires. Nous montrons en particulier que, sous des hypothèses de stationarité et d’isotropie sur le champ aléatoire, en moyenne, le périmètre est toujours biaisé (avec un facteur 4 /π ) contrairement à la courbure totale. Nous illustrons nos résultats sur diﬀérents exemples de champs aléatoires.


Introduction
Understanding the geometry of excursion sets of random fields is a question that receives much attention from both the theoretical and the applied point of view (see [Adl00] for instance). This is partly due to numerous applications in image processing [Ser82,Wor96] for pattern detection, segmentation or image model understanding. Moreover, important strong results have been already obtained especially for smooth Gaussian and related fields [AT07]. This allows to consider some geometrical characteristics of a given image considered as the realization of a random field, related to Minkowski functionals in convex geometry [SKM87] or Lipschitz-Killing curvatures in differential geometry [Thä08]. Roughly speaking, the considered quantities are the surface area, the perimeter and the Euler characteristic, i.e. the number of connected components minus the number of holes (also related to the total curvature), of a black-and-white image obtained by thresholding a gray-level image at some fixed level, corresponding to an excursion set. There exists an abundant literature studying these geometrical features, let us cite for instance [AW09,DBEL17,EL16,KV18,LR19]. Most of these mentioned results rely on strong assumptions on the smoothness of the underlying random fields.
But when making numerical computations in applications, we rarely have access to functions defined on a continuous domain U , we rather have access to the function taken at points on a discrete grid. The main example is the one of digital images that are made of pixels, where the excursion sets are obtained through discrete sets.
The link between the discrete geometry of a set and its "true" underlying continuous geometry has of course been already studied a lot in different fields: for instance in discrete geometry [KSS06,Sva14,Sva15], in systematic sampling [GJ87], in digital topology [Gra71,Pra07] or in mathematical morphology [Mat75,Ser82]. This list is far from being exhaustive.
This discretization procedure also induces a switch of functional framework since piecewise constant functions instead of smooth ones have to be considered. For the perimeter, the nice functional framework of functions of bounded variation [AFP00] allows to unify both approaches by considering perimeter as a function of the level and adopting a weak formulation [BD16].
In our previous paper [BD20], we have introduced functionals that allow us to give (weak) formulas not only for the perimeter but also for the total curvature (related to the Euler Characteristic, by Gauss-Bonnet Theorem) of the excursion sets of a function defined on an open set of R 2 . More precisely, the framework is the following.
Let U = (0, T ) 2 with T > 0, be a square domain of R 2 . Let f be a real-valued function defined on R 2 , and such that for almost every t, the boundary of the excursion set above level t in U is a piecewise C 2 curve that has finite length and finite total curvature. For t ∈ R, we denote the excursion set of f above the level t by Under suitable assumptions on f , we define the level perimeter integral (LP) and the level total curvature integral (LTC) of f , as the functional defined for every h ∈ C b (R), the space of bounded continuous functions on R, by where, denoting by H 1 the 1-dimensional Hausdorff measure, we have and TC is the total curvature of a curve. It is defined, for any piecewise C 2 oriented curve Γ by where κ Γ is the signed curvature of Γ defined at regular points, and α i are the turning angles at the singular points (corners) of Γ. Thanks to the Gauss-Bonnet theorem, the total curvature of the positively oriented curve ∂E f (t) is closely related to the Euler characteristic of E f (t) (see [DC76,p. 274] for instance). Considering for h the constant function equal to 1, we will simply denote LP f (U ) := LP f (1, U ) and LTC f (U ) := LTC f (1, U ).
By the coarea formula ( [AFP00] or [EG92]), LP f (U ) is equal to the total variation of f in U .
To have all three Minkowski functionals (or Lipschitz Killing curvatures), we could also define the level area functional as, where h now needs also to be integrable, h ∈ L 1 (R), and L(E) denotes the Lebesgue measure (area) of a set E. Now, this level area can be written as where H is any primitive of h. Here the integral LA that was defined on the levels t ∈ R has been rewritten as an integral on the domain U . This can also be done for LP and LTC. More precisely, for the level perimeter, when f ∈ C 1 (R), we obtain in [BD16] the following formula, for h ∈ C b (R), and in particular that is the coarea formula. For the level total curvature, when f ∈ C 2 (R), we obtained in [BD20], for h ∈ C b (R), where if u and v are two vectors of R 2 , the notation D 2 f (x).(u, v) stands for u t D 2 f (x)v where here D 2 f (x) is seen as a 2 × 2 symmetric matrix. We also obtained explicit formulas when f is no more smooth but piecewise constant on nice sets (it is then called an elementary function) in [BD20, Equation (17) for LP and (18) for LTC]. We investigate in this paper how this point of view can be adapted to functions that are piecewise constant on a regular tiling, where the geometry of the tiling will also play an important role. Despite the fact that such functions are no more elementary in the sense of our previous paper, this is a natural framework for numerical computations as soon as one has to consider discretization of functions. Hence we will consider here two situations. The first one where we have a tiling of the plane with regular hexagons. The second one, is a more realistic case, where we have a tiling with squares (pixels). Assuming some regularity of f (C 1 or Lipschitz on R 2 for instance), we can use an approximation inequality such as the one of Proposition 4.1, to show that the level area of a discretized version f ε of f converges to the level area of f , when ε goes to 0. In this paper, we will focus on what happens for the level perimeter and the level total curvature of a discretized version f ε of f . The geometry of the tiling is important, and in the case of pixels, the connectivity is not well defined since both 4-and 8-connectivity can be considered. The two cases will be studied. Now, the specificity of our approach here is that we follow our "functional" point of view (through LP and LTC), but also our random field approach, replacing the deterministic function f by a random one X and considering the expectation of LP or LTC. This allows us to provide explicit mean formulas in particular when the random field X is stationary and isotropic.
The paper is organized as follows. In Section 2 we give formulas for the level perimeter integral and the level curvature integral of discrete deterministic functions defined on else an hexagonal or a square tiling. Then, in Section 3, we derive expressions for the expectation of these integral in the case of discrete random fields. More precisely, we are interested in white noise and in positively correlated Gaussian random fields. Now, another way to obtain discrete functions is to discretize a smooth function (or random field). This is what we do in Section 4, and we give the limits as the tile size goes to 0, showing that the level curvature integral behaves well, whereas the level perimeter integral has a bias that we quantify. We illustrate this with some numerical experiments. In the Appendix, we have postponed some technical proofs and also we propose an unbiased way to compute the level perimeter integral.

Geometry of discrete functions
2.1. The hexagonal tiling case On the right: the domains U (black square), U ε (red rectangle) and U ε (blue square).
We first introduce some notations for the tiling with hexagons. For θ ∈ R, we will denote by e θ the unit vector of coordinates (cos θ, sin θ). Let ε > 0 and let us consider a regular tiling with hexagons of "size" ε where the set of the centers of the hexagons is given by The distance between the centers of two neighbouring hexagons is √ 3ε, the side length of the hexagons is ε and the area of each hexagon is 3 √ 3 2 ε 2 . The vertices of the hexagons are the set of points V ε given by V ε = C ε + εe π 6 +n π 3 ; 0 n 5 . On Figure 2.1, we show such a tiling with regular hexagons. The points of C ε are plotted with black stars and the points of V ε are the vertices of the hexagons marked by small red circles. For z ∈ C ε we will denote by D(z, ε) the (open) hexagon of center z and size ε. Notice that the distance between a vertex v ∈ V ε and the centers of its three neighbouring hexagons is equal to ε, that is also the side length of the hexagons.
Finally we will denote by E ε the set of edges. Each edge is a segment of length ε between two neighbouring vertices of V ε and we will sometimes identify an edge w ∈ E ε with its middle point. The set of edges is the union of three sets, depending on the orientation of the edge, and that are denoted by E π/2 ε , E π/6 ε and E −π/6 ε . In order to remove boundary effects, when considering a square domain U = (0, T ) 2 , we will consider the enlarged domain U ε = (− ε 2 , T + ε 2 ) × (− ε 2 , T + ε 2 ) as well as the restricted domain This will ensure that no edge (seen as an open segment) of the tiling in U ε intersects ∂U ε , and that each midpoint w ∈ E ε ∩ U ε is the middle of two centers in C ε ∩ U ε (see Figure 2.1 right).
To give some order of magnitudes, notice that the cardinality of the different sets of points are These equivalents also hold when we consider U ε or U ε in place of U . We denote by PC Hex ε (U ε ) the set of piecewise constant functions on the hexagonal tiling in U ε . A function f ∈ PC Hex ε (U ε ) can be identified with the finite set of values {f (y)} y ∈ Cε ∩ U ε . To have a function that is defined everywhere, we adopt the convention that the value of f on an edge is equal to the mean value of its two neighbouring centers, and the value at a vertex is the mean value of its three neighbouring centers. For f ∈ PC Hex ε (U ε ), we denote for each vertex v ∈ V ε , the three ordered neighbouring values at v by . And for each w ∈ E ε , we denote by f + (w) and f − (w), respectively the maximum and the minimum of the two values of f on the two sides of w.
. The function f has a finite total variation in U ε and for h ∈ C b (R) and H a primitive of h, the level perimeter integral of f satisfies Moreover, the function f is of finite level total curvature integral and the level total curvature integral of f satisfies In particular, Proof. -Let us start with the level perimeter integral. Since f ∈ PC Hex ε (U ε ), for any t ∈ R, the excursion set E f (t) ∩ U ε is a union of hexagons (or parts of hexagons on the boundary), and an edge w ∈ E ε is part of the boundary of E f (t) in U ε if and only if f − (w) < t f + (w). We also recall that all edges in E ε have the same length, that is equal to ε and that an edge in U ε is entirely contained in U ε . Moreover, since U ε is bounded, t → Per(E f (t), U ε ) is piecewise constant with compact support and therefore we have, for h ∈ C b (R), denoting H a primitive of h,  For the level total curvature integral the computations are similar. The boundary of an excursion set E f (t) in U ε is a curve that is piecewise linear since it is made of edges in E ε . Its curvature at regular points is then 0, and it has only corner points at vertices v ∈ V ε , where the turning angle is else π Figure 2.2). Therefore The centers (set C ε ) of the squares are the black stars, and the vertices (set V ε ) are the points marked by a red circle. On the right: the domain U (black square) and U ε (red square).

The square tiling case
We now consider the case of a tiling with squares. This is the case used in practice for digital images since they are defined on (square) pixels (contraction of picture elements). Let ε > 0 and let us consider a regular tiling with squares of "size" (side length) ε where the set of the centers of the squares is given by The side length of the squares is ε and the area of each square is ε 2 . The vertices of the squares are the set of points V ε given by We will denote by E ε the set of edges. Each w ∈ E ε is a segment of length ε that is else horizontal or vertical. For z ∈ C ε we will denote by D(z, ε) the (open) square of center z and size ε. Finally, notice that the distance between a vertex v ∈ V ε and the centers of its four neighbouring squares is equal to ε √ 2/2. When considering a square domain U = (0, T ) 2 and ε > 0, we will define here the restricted domain U ε = (0, ε T ε ) 2 . For the enlarged domain U ε , since we already have that each midpoint w ∈ E ε ∩ U ε is the middle of two centers in C ε ∩ U (see Figure 2.3 right), we can simply set U ε = U .
Let us notice that here we have When dealing with a tiling with squares, the definition of connectivity is not unique. Indeed we can say that two squares are neighbours if they have a common edge (this is the 4-connectivity), or only as soon as they have a common corner (this is the 8-connectivity). Now, in fact, these two connectivities are "complementary". Indeed, if we want a discrete version of the Jordan curve theorem to hold, we have to state it in the following way ( [Ros79]) : the complement of a 4-connected simple closed discrete curve (sequence of squares) is made of exactly two 8-connected components.
We will denote by PC Sq ε (U ) the set of functions f defined on U that are piecewise constant on the tiling with regular squares of size ε > 0. Such a function can be simply identified to the finite set of values {f (z)} z ∈ Cε ∩ U . The value of f along an edge is taken as being the mean value of its two neighbouring centers, while its value at a vertex is given by the mean value of its four neighbouring centers. Since we have to consider two different total curvatures according to the choice of connectivity, we write TC 4 (∂E f (t) ∩ U ε ) and TC 8 (∂E f (t) ∩ U ε ) such that, for h a bounded continuous function, We denote for each vertex v ∈ V ε , the four ordered neighbouring values at v by . And for each w ∈ E ε , we denote by f + (w) and f − (w), respectively the maximum and the minimum of the two values of f on the two sides of w.
. The function f has a finite total variation in U ε and for h ∈ C b (R) and H a primitive of h, the level perimeter integral of f satisfies Moreover, the function f is of finite level total curvature integral and the level total curvature integrals of f satisfy and where c(v) = cross denotes the event that the configuration at v is "a cross " (meaning that f (1) and f (2) are achieved at two "opposite" squares (see Figure 2.5)).
Proof. -Let us start with the level perimeter integral. Since f ∈ PC Sq ε (U ), for any t ∈ R, the excursion set E f (t) ∩ U ε is a union of squares, and an edge w ∈ E ε is part We also recall that all edges in E ε have the same length, that is equal to ε.
and the configuration at v is not a cross (since in that case the set {f t} is made of two adjacent squares, see the right-most figure).
For the level total curvature integral the computations are also similar to the ones in the case of hexagons. However, we have to consider the two different types of connectivity. The boundary of an excursion set E f (t) in U ε is a curve that is piecewise linear since it is made of edges in E ε . Its curvature at regular points is then 0, and it has only corner points at vertices v ∈ V ε , where the turning angle β is (see Figure 2.4): and the configuration at v is a cross (left figure), the turning angle at a vertex v is π = π/2 + π/2 (in 4−connectivity, since it is equivalent to the "zoom" presented in the middle figure) or −π = −π/2−π/2 (in 8−connectivity, see the "zoom" on the right figure).
, then β = 0 if the configuration at v is not a "cross" (see Figure 2.4), whereas if the configuration at v is a cross (see Figure 2.5), then β = π in 4−connectivity and β = −π in 8−connectivity. Therefore, For LTC 8 f (h, U ε ) the computation is the same, except that the +π in front of the second sum is changed into −π.
A convenient way to get rid of the connectivity ambiguity is to consider a kind of "6-connectivity" by setting Then, the "cross" configuration doesn't appear anymore in the formula, since, using the above results, we simply have Remark 2.3. -Let us note that these formulas are of course linked with numerical computations of discrete topology. Actually, when considering a set E ⊂ U we can choose f ∈ PC Sq ε (U ) corresponding to its discretization of size ε by taking f (z) = 1 when z ∈ C ε ∩ E and f (z) = 0 otherwise. Now, since the values {f (z)} z ∈ Cε ∩ U are in {0, 1} and those of f in {0, 1/4, 1/2, 3/4, 1}, one can take h ∈ C b (R) a non-negative function with support in (3/4, 1) such that R h = 1 3/4 h = 1. On the one hand, we clearly have On the other hand, choosing Among them we can distinguish three configurations. The first one when f (2) (v) = 0 and f (3) (v) = 1, only contributes to the second sum for cross events with +1; the other ones contribute only to the first sum with +1 when Hence it is enough to count the number of such configurations. By the Gauss-Bonnet theorem, since the Euler characteristic corresponds to the total curvature divided by 2π, this coincides with the algorithms proposed for computing the Euler characteristic of discrete sets as for example the function bweuler in Matlab [Gra71,Pra07] with respect to the two different connectivities.

The mean geometry of discrete random fields
In this section we introduce (Ω, A, P) a complete probability space and replace Then LP and LTC are now real random variables and we will focus on their mean values given by expectations when they can be defined.

Perimeter and total curvature of a white noise
In this first part we investigate the case of a white noise obtained choosing {X(z)} z ∈ Cε ∩ U ε independent identically distributed real random variables of common distribution function F . We will note X Hex ∈ PC Hex ε (U ε ) and X Sq ∈ PC Sq ε (U ε ) according to the tiling when considering LTC.
Proposition 3.1. -Assume that the X(z), z ∈ C ε ∩ U ε are independent identically distributed on R with distribution function F . Then, for h ∈ C b (R) ∩ L 1 (R), for both the hexagonal and the square tiling case, LP and LTC have finite expectation and we have In the hexagonal case, we have while in the square case we have where we have the sign + for LTC 4 and the sign − for LTC 8 .
Proof. -Note that choosing h ∈ C b (R) ∩ L 1 (R) ensures that we can choose a bounded primitive function H in such a way that H(X(z)) are bounded random variables and therefore they all have finite expectation. It ensures that LP and LTC have finite expectation as finite sums of such variables. Now, since the X(z), z ∈ C ε , are independent identically distributed on R with distribution function F , we have for any w ∈ E ε , (X − (w), X + (w)) d = (min(X 1 , X 2 ), max(X 1 , X 2 )) where X 1 and X 2 are independent and follow the distribution F . Therefore, Hence we have to compute where we use the notations of [AN01], meaning that X 1, n . . . X n, n are the ordered observations of X 1 , . . . , X n , for n 2. We will denote by F k, n the distribution function of X k, n . In this setting, we have that for 1 k n, where F k, n (t) = P(X k, n t). Now, we can write by Fubini Theorem, since h ∈ L 1 (R), and this completes the proof of the formula for the level total perimeter.

TOME 4 (2021)
For the level total curvature, the computation is very similar, except that for the hexagonal tiling we have now three independent random variables X 1 , X 2 and X 3 of the same distribution F , and we consider their max, min and median. We have Now, as above we can write and this completes the proof of the formula for the level total curvature integral. Finally for the square tiling we have now four independent random variables X 1 , X 2 , X 3 and X 4 of the same law F to order. We have Now, for any vertex v we also have 1 I c(v) = cross d = 1 I cross 1 I X 2, 4 t < X 3, 4 with E 1 I cross 1 I X 2, 4 t < X 3, 4 X 2, 4 , X 3, 4 = 1 3 1 I X 2, 4 t < X 3, 4 , since there are 2 configurations over 6 possible ones to get a cross. Thus Finally, we get Notice that, as a consequence, we get which is, up to a constant factor that depends on the geometry of the tiling (angles between the edges and number of vertices), the same as E(LTC X Hex (h, U ε )). Let us also remark that choosing a distribution F such that F (1 − F ) ∈ L 1 (R) we can deduce that, for almost every t ∈ R, and similarly for the mean values of total curvatures. We insist on the fact that this holds for almost every t. Actually, considering a Bernoulli noise of parameter p ∈ (0, 1), the distribution function F p has two jumps at t = 0 and t = 1 and to compute the mean perimeter of the excursion set at these jumps values. We illustrate these results in the case of tiling with squares on Figures 3.1 and 3.2. Here we consider the square domain U = (0, 1) 2 and ε = 1/200. The random field is a white noise with uniform distribution on [0, 1] of size 200 × 200 pixels, i.e. here F is continuous with This example leads us to two remarks. The first one is that the empirical curves on one large sample are very close to the theoretical mean values, suggesting that the variances of LP X and LTC X are very small. Computing these variances is doable in theory and it would be an interesting direction for future investigations. The second remark is the question of knowing if there is a relationship between the . On the right, the total curvature: empirical values (plotted with stars) and theoretical mean total curvatures given by t → 2π values where t → E(TC d (∂E X (t, U ε )))) crosses 0 and percolation thresholds. Indeed in the hexagonal case, the percolation threshold is p c = 0.5 and this is also the value t c at which t → 2π And in the square case, the percolation threshold is p c 0.593 and t c = 1 2 + 1− √ 5 2 0.618 is the positive zero of t → E(TC 4 (∂E X (t), U ε )), hence it seems that |p c − t c | is "small", and studying this fact to know if it can be generalized would be interesting.

Perimeter and total curvature of positively correlated Gaussian fields
It is more difficult to get explicit computations when considering non-independent random variables without adding assumptions on their distribution. In this section we consider the discretization of a standard centered Gaussian stationary field X = (X(x)) x ∈ R 2 , that is also positively correlated, with covariance function ρ, meaning that Cov(X(x), X(y)) = ρ(x − y) 0 with ρ(0) = 1. Note that the case where the variance of X(x), given by σ 2 := ρ(0), is not equal to 1 can easily be deduced from this one considering X/σ. For ε > 0 we consider the discretization of X given by the set of the values (X(z)) z ∈ Cε . The main quantities of interest will be Note that the behavior of β θ (ε) is linked with the regularity of the field. Actually, mean square regularity is related to sample paths continuity for Gaussian fields (see [AT07] for instance). Adding stationarity, one can deduce directional regularity from the behavior of β θ (ε) when ε tends to zero. For instance, when there exists α ∈ (0, 1] and λ 2α (θ) > 0 such that ε −2α β θ (ε) → λ 2α (θ) as ε goes to 0, one can find a modification of X such that, for x ∈ R 2 , t ∈ R → X(x + te θ ) is almost surely α -Hölder continuous for any α < α (we refer the interested reader to [Bie19, Part 2]). Let us also emphasize that when X is a.s. C 1 one has ε −2 β θ (ε) → λ 2 (θ), where λ 2 (θ) = Var(∂ θ X(0)), with ∂ θ X the partial derivative of X in the direction e θ .
When the field X is isotropic this value does not depend on θ and the common value denoted as λ 2 is usually called second spectral moment.
In the following Theorem 3.2 we focus on asymptotics for mean LP and LTC obtained for the discretization of X on a tiling as ε goes to zero. Our results are mainly based on ordered statistics of order 2 for LP, 3 for LTC in the hexagonal tiling case and 4 for LTC in the square tiling case. Even with a Gaussian distribution, there are few results available in our dependent setting and we need to impose extra assumptions on the dependency given by the covariance function. In particular we are working with positively correlated variables meaning that ρ is a non-negative function. Moreover, we will need the following assumptions: (A1) there exists α ∈ (0, 1] and real numbers λ 2α (θ) 0 such that for any edge orientation e θ of the tiling. (A2) Assumption (A1) holds and ρ(εe θ ) = ρ(εe π/2 ) for any edge orientation e θ of the tiling, hence we write λ 2α the common value of λ 2α (θ). (A3) Assumption (A2) holds for the square tiling and ρ(εe π/2 )−ρ(ε(e 0 +e π/2 )) > 0 and 1 − 2ρ(εe π/2 ) + ρ(ε(e 0 + e π/2 )) 0 with Theorem 3.2. -We consider the discretization X ε of a centered stationary standard Gaussian and positively correlated random field X.
Finally, under (A3), then where we recall that LTC 6 := 1 2 (LTC 4 + LTC 8 ). The proof of this theorem is technical and it is postponed to Appendix A.1 Remark 3.3. -When X is assumed to be a.s. C 3 and isotropic, we have for all t ∈ R, E (Per(E X (t), U )) = 2L(U )C * 1 (X, t) and E (TC(∂E X (t) ∩ U )) = 2πL(U )C * 0 (X, t), where C * 1 and C * 0 are the Lipschitz-Killing curvatures densities (see [BDBDE19]), given by where λ 2 denotes the spectral moment of X corresponding to Var(∂ 1 X(0)) = Var(∂ 2 X(0)) by isotropy. Since the field X will satisfy (A1) with α = 1 and λ 2α (θ) = λ 2 for any orientation θ by isotropy, as soon as ρ is non-negative in order to ensure the positive dependence assumption. Hence, since it also satisfies (A2) by isotropy, by Theorem 3.2 we obtain It follows that we have a weak-convergence Now, if we assume moreover that ρ(x) = ρ( x 2 ) with ρ a non-negative function that is C 2 on a neighbourhood of 0 and such that ρ (0) < 0 and ρ (0) > 0 we easily check, using Taylor formula, the additional assumptions for the square tiling and also obtain the weak-convergence ) .

Discretization of smooth functions
In this section, we will study the limits of LP fε and LTC fε as the size ε of the tiling goes to 0, when f is a smooth function. In particular, we would like to know if the limits coincide with LP f and LTC f . For this aim, let U = (0, T ) 2 with T > 0 and we recall first the main formulas obtained in our previous paper [BD20] for a smooth (C 2 ) function f on R 2 , given by (1 .1) and (1.3) When h is also assumed to be C 1 on R, denoting by H a primitive of h, a simple computation leads to, for any vector e ∈ R 2 , where ·, · is the Euclidean scalar product on R 2 . Using this formula with e = ∇f (x) ⊥ |∇f (x)| , we get, Hence we are reduced back to the case h = 1. In a similar way, decomposing a bounded continuous function h into h = h 1 − h 2 , where h 1 = h + 2 h ∞ and h 2 = 2 h ∞ are both positive continuous and bounded functions, we get that Observe that we also have LP Since we can proceed similarly for the level total curvature we will focus on the case where h = 1 in the following.
Usually, to infer an approximation error between the integral of a function and its discretized version, one uses an approximation inequality like the Koksma-Hlawka inequality [PS15]. Now here we will need a similar result that is given by the following proposition.
Proposition 4.1 (Approximation Inequality). -Let W be a rectangular domain in R 2 . Let g be a bounded, Lipschitz function defined on R 2 . Let us consider a regular tiling with a shape H ε (that can be an hexagon, a square or a rhomb) of "size" ε. Let a ε = L(H ε ) be the area of H ε and let d ε be the diameter of H ε (that is the maximal distance between two points of H ε ). Let C ε be the set of centers of the tiles. Then Let A ⊂ W be an open or closed subset of W . Then the cardinality of C ε ∩ A is bounded: where B(0, d ε ) is the ball of center 0 and radius d ε , and ⊕ denotes the Minkowski sum, defined by A ⊕ B := {x + y ; x ∈ A and y ∈ B}.
Proof. -Let us start with the first part of the proposition. We notice that a ε y ∈ Cε ∩ W g(y) = y ∈ Cε ∩ W Hε g(y) dz.

ANNALES HENRI LEBESGUE
Bounding L(W ε ) by L(W ) + d ε H 1 (∂W ) + 4d 2 ε , we have the result. For the second part of the proposition, we first notice that A ⊕ B(0, d ε ), we have the announced inequality.

Notations
When f is a C 2 function on U ⊂ R 2 , we will use the notations for the partial derivatives of f at point x ∈ U .

Limits as the hexagon's size goes to 0
Let U = (0, T ) 2 be a fixed domain. Let ε 0 > 0, and let us consider a tiling with regular hexagons of size ε ∈ (0, ε 0 ]. Let f be a C 2 function defined on U ε 0 . We then consider a discretized version where the D(z, ε) are the hexagonal tiles, and the boundary conditions are defined as in Section 2.1. The formulas for LP fε (U ε ) and LTC fε (U ε ) were given in Proposition 2.1. We are interested in their limits as ε goes to 0, and the links with LP f (U ) (Equation (1.2)) and LTC f (U ) (Equation (1.4)).
Theorem 4.2. -Let f be a function defined on U and assume that f is C 2 on For ε ∈ (0, ε 0 ], let f ε be the discretized version of f on U ε . Then, C being a numerical constant independent of everything.
Then, there exists a constant C Hex Proof. -We detail here the result concerning the level perimeter integral of f ε , as ε goes to 0. We assume that f ε is the discretized version of a C 2 function f defined on U ε 0 with U ε ⊂ U ⊂ U ε ⊂ U ε 0 for ε ε 0 . We have by Proposition 2.1 that Let w ∈ E ε be an edge, that is the boundary between two neighbouring hexagons, and let z w be the center of the right-most hexagon (i.e. among the two hexagon centers, z w is the one that has the largest first coordinate). The center of the other hexagon is then z w + ε √ 3e ⊥ w , where e w is the unit length vector oriented as the edge w, and e ⊥ w is its π 2 -rotation. See also Figure 4.1 left. Then where |r 1 (z w , ε)| 3 2 ε 2 D 2 f ∞ . Now, each center z ∈ C ε is the z w of three different w, with respective normal orientation e ⊥ w equal to e 2π/3 , e π = −e 0 and e 4π/3 = −e π/3 . Therefore, we can rewrite
The proof for LTC also relies on Taylor formulas but now of order 2 instead of 1 and needs a clever grouping of vertices (see Figure 4.1 right). The details are postponed to Appendix A.2.

Limit as the square's size ε goes to 0
Again, let U = (0, T ) 2 be a fixed domain. Let ε 0 > 0, and let us now consider a tiling with squares of size ε ∈ (0, ε 0 ]. Let f be a C 2 function defined on U ε 0 . We then consider a discretized version where the D(z, ε) are the square tiles, and the boundary conditions are defined as in Section 2.2. The formulas for LP fε (U ε ) and LTC fε (U ε ) were given in Proposition 2.2. We are interested in their limits as ε goes to 0, and the links with LP f (U ) (Equation (1.2)) and LTC f (U ) (Equation (1.4)).
Theorem 4.3. -Let f be a function defined on U and assume that f is C 2 on U ε 0 with ∇f ∞ := max U ε 0 ∇f < +∞ and D 2 f ∞ := max U ε 0 D 2 f < +∞. For ε ∈ (0, ε 0 ], let f ε be the square discretized version of f on U ε . Then, C being a numerical constant independent of everything.
ANNALES HENRI LEBESGUE Proof. -Let us consider the level perimeter integral of f ε , as ε goes to 0, with f ε the square discretized version of a C 2 function f defined on U ε 0 . The computations here will be very similar to the ones in the hexagonal case. We have, by Proposition 2.2, that Let w ∈ E ε be an edge, that is the boundary between two neighbouring squares, and let z w be the center of the left-most (if the edge is vertical), or bottom-most (if the edge is horizontal) square. The center of the other square is then z w + εe ⊥ w . See Figure 4.2 left. Then where |r 1 (z w , ε)| ε 2 D 2 f ∞ , by Taylor formula. Now, each center z ∈ C ε is the z w of two different w, with respective normal orientation e ⊥ w equal to e 0 and e π/2 . Therefore, we can rewrite with |r 2 (ε)| Cε (L(U ) D 2 f ∞ + H 1 (∂U ) ∇f ∞ ), using the fact that |C ε ∩ U | CL(U )ε −2 . Then, since the area of each square D(z, ε) is equal to a ε = ε 2 , by Proposition 4.1, we finally get . This ends the proof for the level perimeter integral.
The proof for LTC is similar but more technical, and it is postponed to Appendix A.3.

Discretizing a smooth random field
In this section, we will see what happens to the mean level perimeter integral and to the mean level total curvature integral of a discretized smooth stationary random field, in both the hexagonal tiling and the square tiling cases. Roughly speaking, we will see that the perimeter is always biased, whereas the total curvature is not, under an additional isotropy assumption.
In all this section, as previously, we consider a fixed domain U = (0, T ) 2 .
Proposition 4.4. -Let X be a stationary C 2 random field on R 2 such that have finite expectations for some ε 0 > 0. Then LP X (U ), LP Hex X (U ) and LP Sq X (U ) are in L 1 (Ω).
Let us consider the discretization X ε ∈ PC Hex ε (U ε ) as in (4.1), respectively X ε ∈ PC Sq ε (U ε ) as in (4.4). Then LP Xε (U ε ) converges to LP Hex X (U ), respectively to LP Sq X (U ), in L 1 (Ω), as ε goes to 0. Moreover, Under the additional assumption that X is isotropic we have To give some hints on the numerical values: 2 1.27. This shows that, whatever the field, whatever the smallness of the hexagons, there is always a bias when approximating the level perimeter integral of the field X by the one of its discretized version on an hexagonal tiling. There is also a bias on a square tiling, except if the smooth field X has a gradient that is everywhere aligned with e 0 or e π 2 . The strongest bias is obtained when the gradient is everywhere aligned with the diagonal directions e π 4 or e − π 4 . These last remarks are consequences of the proofs below.
Proof. -Under our assumptions it is clear that LP Xε (U ε ), LP X (U ), LP Hex X (U ) and LP Sq X (U ) are in L 1 (Ω). Moreover we also have that C Hex LP (X, U ) and C Sq LP (X, U ) are in L 1 (Ω) so that the convergence results hold taking expectation from Theorems 4.2 and 4.3. According to the beginning of Section 4, by Fubini theorem and the stationarity of X, we have Now, a simple computation shows that for any θ ∈ R, we have In the case of a tiling with squares, since for any θ ∈ R we have we obtain When the smooth stationary random field X is moreover isotropic, ∇X is rotationally invariant and, according to [BB99,Proposition 4.10], its gradient direction ∇X(x)/ ∇X(x) is independent from ∇X(x) and uniform on S 1 . Thus we get, for any θ ∈ [0, 2π), This shows that in the isotropic case and since 4 π > 1, there is always a bias.
Remark 4.5. -Assuming moreover that E( ∇X 2 ∞ ) < +∞, for any h : R → R bounded C 1 function with derivative h ∈ C b (R), and denoting by H a primitive of h, then the random field H • X will satisfy the assumptions of Proposition 4.4. By linearity this allows us to state the convergence results for LP Xε (h, U ε ) and obtain in the isotropic case the weak convergence as remarked in the Gaussian setting of Section 3.2. This is illustrated on Figure 3.3 and it explains why computing perimeters from discrete images is not easy. However, solutions exist to obtain non-biased estimates of the level perimeter integrals from pixelated images, and we propose such a solution in the Appendix B. The idea behind the unbiased estimation of the perimeter given in the Appendix B in the square tilling framework is to linearly interpolate the function inside each dual square and approximate the boundary of each level set by a polygonal line where now segments are not only horizontal or vertical (as for the discretized function). We show in the Appendix why this linear interpolate provides unbiased estimates of the level perimeter integral. See also Figure 4.3, where we used a non-Gaussian smooth isotropic shot noise field as considered in [BD20].
For the level total curvature, things are different: there is no bias. Intuitively this can be explained by the fact that the total curvature is related to the Euler characteristic that counts the number of connected components and the number of holes, and these numbers remain (almost) the same when the function is discretized on very small hexagons or squares.

ANNALES HENRI LEBESGUE
But for j = 1, 2, by Markov inequality. Since lim ε→0 P(|∂ j X(0)| < ε 1/2 ) = P(∂ j X(0) = 0) = 0 by assumption, we can conclude that L(U ε (X, U )) converges to 0 in L 1 (Ω) and thus in probability. Hence D 2 X ∞ L(U ε (X, U )) converges to 0 in probability and since the variables { D 2 X ∞ L(U ε (X, U )); ε ∈ (0, ε 0 ]} are uniformly integrable (because they are uniformly bounded by L(U ) D 2 X ∞ ) we also have that D 2 X ∞ L(U ε (X, U )) converges to 0 in L 1 (Ω). According to Theorem 4.3 this implies that LTC d Xε (U ε ) converges to LTC Sq X (U ) in L 1 (Ω). Moreover by stationarity we obtain Now let us assume also that X is isotropic and remark that our assumption implies that ∇X(0) = 0 a.s. Hence let us define Θ as the argument of the gradient ∇X(0) and write where g is the π periodic function piecewise C 1 defined by g(θ) = 1 if θ ∈ (0, π/2), g(θ) = −1 if θ ∈ (π/2, π) and g(θ) = 0 if θ ∈ {0, π 2 }. Then, using the Fourier series of g and the isotropy of X, we can show (see the Appendix A.4) that Now let us consider the hexagonal tiling case. The first part follows similarly using Theorem 4.2, once we have remarked that Then, by stationarity we obtain Assuming moreover that X is isotropic, again our assumption implies that ∇X(0) = 0 a.s. As previously we define Θ as the argument of the gradient ∇X(0) and write where g is now the π-periodic function define on [−π/6, 5π/6] by g = 1 I (π/6, 5π/6) − 21 I (−π/6, π/6) . Here again, using the Fourier series of g and the isotropy of X (see the technical details in the Appendix A.4), we can show that Let us remark that assuming moreover that E( ∇X 3 ∞ ) < +∞ and E( D 2 X 2 ∞ ) < +∞, for any h : R → R bounded C 2 function with h and h bounded, denoting by H a primitive of h, the random field H • X will also satisfies assumptions of Proposition 4.6 as soon as P(h(X(0)) = 0) = 0. By linearity this allows us to state the convergence results for LTC d Xε (h, U ε ) and obtain in the isotropic case the weak-convergence (according to the set of admissible test functions h) as remarked in the Gaussian setting. It explains that there is no bias on the level total curvature when discretizing a smooth isotropic stationary random field. This is illustrated on Figure 3.3 for a Gaussian field and on Figure 4.3 for a smooth isotropic shot noise field. This last figure also shows the robustness of TC with respect to the scales of resolution.
Remark 4.7. -Let us notice that by the formula for LTC Xε in Proposition 2.1, we have in the hexagonal tiling case, We see that it involves X (2) (v), that is the median value around v. In the square tiling case, there are two median values given by X (2) (v) and X (3) (v). Hence we have here shown an interesting link between median value and curvature since, as ε goes to 0, the limit of LTC Xε (U ε ) is, in expectation, the curvature of the smooth function X. This link was already well-known in the field of mathematical image processing where the median filter, a commonly used filtering method for images, converges (when iterated) to the so-called mean curvature motion (see [Cao03] for instance). , a sample on (0, 1) 2 of a smooth shot noise random field X with Gaussian kernel [BD20] on a digital image of size 4000 × 4000 pixels, i.e. field X ε with here ε = 1/4000. On the right, same field X but now discretized on a 62 × 62 pixels grid. It corresponds to a "scale" s = 6, since 62 = 4000/2 s with s = 6. Bottom line: on the left, the perimeter of the excursion sets of X ε as a function of the level t (in abscissa) -only values of t multiples of .5 have been used, hence the plot is a polygonal curve -for different scales s (different colors). The plain curve is the perimeter as defined for discretized fields, the dashed curve is the unbiased perimeter computed as in Appendix B, and the dotted curve is 4/π times the unbiased perimeter. It fits quite well the plain curve, illustrating Proposition 4.4. On the bottom right, the total curvature TC 6 of the excursion sets of X ε as a function of the level t (in abscissa) -again, only values of t multiples of .5 have been used, hence the plot is a polygonal curvefor different scales s (different colors). As intuitively expected, and except for the coarsest scale s = 6, whatever the size of the discretization, the values of the total curvature remain almost the same.
First assume that we are in the hexagonal tiling case. Then we write 3 √ 3 L(U ) for 1 i 3, and by (A1), It follows that On the other hand, for the square tiling case, It follows that Now, let us consider the level total curvature where we assume moreover that ρ(εe θ ) = ρ(εe π/2 ) for any edge orientation and denote λ 2α the common value in view of (A2).
Now, since O ε (U, f ) is defined as the set of x ∈ U such that ||∂ 1 f (x)| − | √ 3∂ 2 f (x)|| < 6ε D 2 f ∞ , by Taylor formula we have that, if x ∈ O ε (U, f ) and y ∈ B(0, 3ε) then, and this concludes the proof of Theorem 4.2.
Therefore, the Fejer sum σ N (g) = 1 N N −1 n=0 S n (g) also converges pointwise towards g with |σ N (g)(θ)| g ∞ = 1 for all N 1 so that by Lebesgue theorem we have