Mariem Abaach
email: mariem.abaach@u-paris.fr

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

Elena Di

Elena Di Bernardino
email: elena.dibernardino@unice.fr

Testing marginal symmetry of digital noise images through the perimeter of excursion sets

Keywords: Primary 60G60, 60D05 ; Secondary 62F03, 60G10 Binary image, Crossings, Excursion sets, Test of symmetry, Threshold procedure

Introduction

The stochastic modeling of images by spatial random fields allows to set up a convenient statistical framework for different issues in image processing for image denoising, pattern detection, segmentation or classification [START_REF] Pratt | Digital Image Processing[END_REF][START_REF] Jähne | Digital Image Processing 6th Edition[END_REF][START_REF] Desolneux | Stochastic methods for image analysis[END_REF]. A particular interest directed towards the study of geometric attributes and features of objects has grown these last years. Roughly speaking the considered quantities are the surface area, perimeter and Euler characteristic, which is a topological invariant that is equal to the number of connected components minus the number of holes, in a black and white binary image obtained by thresholding a gray image at some fixed level. These quantities are related to intrinsic volumes and Minkowski functionals particularly studied in the integral and stochastic geometry [START_REF] Chiu | Stochastic Geometry and Its Applications[END_REF][START_REF] Schneider | Stochastic and integral geometry[END_REF] or Lipschitz-Killing (LK) curvatures of excursion sets (see [START_REF] Thäle | 50 years sets with positive reach -a survey[END_REF] for a formal introduction to the subject). Those functionals are robust and efficient shape descriptors that have been largely applied to a variety of domains: cosmology (e.g. the morphological analysis of cosmic microwave background [START_REF] Schmalzing | Minkowski functionals used in the morphological analysis of cosmic microwave background anisotropy maps[END_REF][START_REF] Marinucci | Testing for non-Gaussianity on cosmic microwave background radiation: a review[END_REF]), astrophysics (e.g. modeling galaxy formation [START_REF] Gott | Genus Topology of Structure in the Sloan Digital Sky Survey: Model Testing[END_REF]), military (e.g. mine field detection, [START_REF] Lake | Identifying minefields in clutter via collinearity and regularity detection[END_REF]) medicine (e.g. brain imaging, [START_REF] Flandin | Topological inference[END_REF], study of synthesized 2D digital mammograms, [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF][START_REF] Bernardino | Statistics for Gaussian random fields with unknown location and scale using Lipschitz-Killing curvatures[END_REF]).

Very important results have been obtained in the framework of smooth stationary random fields, especially for Gaussian related fields (see [START_REF] Adler | The Geometry of Random Field[END_REF]). In this framework, the theoretical means of LK curvatures are explicitly computed with respect to the parameters of the fields [START_REF] Wschebor | Surfaces aléatoires: mesure géométrique des ensembles de niveau[END_REF][START_REF] Worsley | Local maxima and the expected Euler characteristic of excursion sets of χ 2 , F and t fields[END_REF][START_REF] Azaïs | Level sets and extrema of random processes and fields[END_REF] and can be efficiently estimated from images [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF]. The consistent estimation of related variances had been recently studied in [START_REF] Bernardino | Statistics for Gaussian random fields with unknown location and scale using Lipschitz-Killing curvatures[END_REF]. Moreover, Central Limit Theorems (CLT) have been proven for Gaussian fields [START_REF] Estrade | A central limit theorem for the Euler characteristic of a Gaussian excursion set[END_REF][START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF] that paved the road for statistical hypothesis tests. Several extensions for non-Gaussian fields, namely shot noise random fields are also available for means of LK curvatures [START_REF] Lachièze-Rey | Bicovariograms and Euler characteristic of random fields excursions[END_REF] and associated CLT's [START_REF] Lachièze-Rey | Normal convergence of nonlocalised geometric functionals and shot-noise excursions[END_REF]. However the framework of smooth fields, indexed by continuous space variables, is not appropriate with the discrete framework induced by digital images. Let us remark that links with discrete setting have been recently explored in [START_REF] Biermé | The effect of discretization on the mean geometry of a 2D random field[END_REF].

In this paper we consider this discrete framework for digital images for which white noise is well defined. More precisely, we assume that pixels values are given by (X i,j ) i,j≥1 a sequence of independent and identically distributed (i.i.d.) variables within an observation window S. We have in mind the residuals of a denoising procedure or a linear regression, or the difference between two images as for brain activities study. Therefore, we are interested to test the natural symmetry hypothesis, that is to know whether X i,j is drawn from a symmetrical distribution. This assumption of marginal distribution symmetry is called here the null hypothesis H 0 and includes Gaussian or Student distributions. Being able to formally test for symmetry hypotheses is an important topic in many fields since this basic assumption contains important information regarding the underlying model, that would allow to validate it before further investigation. This problem has been largely considered in the literature for i.i.d. sample, especially due to its importance for time series analysis. Different approaches to construct a symmetry test include the use of empirical distribution for known center of symmetry [START_REF] Butler | A test for symmetry using the sample distribution function[END_REF], the characteristic symmetry function and study of its behavior to know if the distribution is symmetric or not [START_REF] Csörgő | Testing for Symmetry[END_REF] with unknown center. Other tests rely on the study of the skewness of the distribution function as initiated by [START_REF] Gupta | An asymptotically nonparametric test of symmetry[END_REF]; see [START_REF] Mira | Distribution-free test for symmetry based on Bonferroni's measure[END_REF] for a free-distribution test based on Bonferroni's measure. For an overview and comparison of the existing methods we refer to [START_REF] Quessy | On Consistent Nonparametric Statistical Tests of Symmetry Hypotheses[END_REF] and [START_REF] Ivanovic | Comparison of symmetry tests against some skew-symmetric alternatives in i.i.d. and non-i.i.d. setting[END_REF].

In our original approach we attempt to use our digital image setting and the geometrical features of its excursion sets. More precisely, for a given threshold t ∈ R we consider the binary image given by the i.i.d. Bernoulli variables (1 {Xi,j ≥t} ) i,j≥1 , coloring black pixels for values equal to one and considering the area of black components as well as its perimeter inside a given domain. This allows us to build an unbiased estimator of the perimeter, without border effects. This framework is close from the one used by [START_REF] Ebner | Goodness-of-fit tests for complete spatial randomness based on Minkowski functionals of binary images[END_REF] to set up a goodness of fit test for complete spatial randomness. In the aforementioned paper the authors consider the counting process obtained from a homogeneous Poisson point process such that the observed values (X i,j ) are i.i.d. Poisson random variables. They obtain explicit expectation and covariance matrix for the three Lipschitz-Killing curvatures at a fixed level as well as CLT to build their test. In the present paper, we focus on the perimeter and obtain a multivariate CLT with respect to several threshold levels with an explicit covariance matrix, under a general i.i.d. assumption, in a dense tiling framework. The choice of our unbiased perimeter estimator permits us to preserve the symmetry of the distribution, in contrast with perimeter with border. This lead us to construct the statistical test of symmetry, by considering the behavior of the ratio between two specific thresholds. We also explore the asymptotic in large threshold limit regime and obtain the equivalence between the asymptotic behavior of the expectation value and the variance establishing a joint auto-normalized Central Limit Theorem for dense tiling and large threshold limit regime in a white noise framework, similarly to [START_REF] Dalmao | Asymptotic normality of high levellarge time crossings of a Gaussian process[END_REF] for crossings in dimension 1.

The paper is organized as follows. In Section 2 we introduce our mathematical framework for binary images and we proceed to the construction of an unbiased estimator for the associated perimeter. Then, in Section 3 we study the first and second moments of the perimeter, which allows us to prove asymptotic normality results with an explicit covariance matrix, for which we consider several consistent estimators. The construction of our consistent symmetry test via empirical accessible test statistics is proposed in Section 4. Numerical evaluations are presented in Section 5. The case of large threshold values is investigated in Section 6, that yields to a new a central limit theorem and test statistics illustrated with numerical comparisons. Finally, we postpone the technical proofs to Section 7, and additional numerical evaluations to the Appendix.

Mathematical framework

Construction of the binary image

Square tiling Let m be an integer with m ≥ 2, without loss of generality, we consider our observation window as the unit square S = [0, 1] 2 and we divide our window into m 2 pairwise disjoint squares. This provides a regular tiling of S with squares of "size" (side length) 1/m, i.e.,

C (m) i,j := [ i -1 m , i m ] × [ j -1 m , j
m ], for i, j ∈ {1, . . . , m}.

The C

(m) i,j will be referred to as cells. We denote by E m the set of edges in S = (0, 1) 2 , each w ∈ E m is a segment of length 1/m. C m i,j . Let (X i,j ) 1≤i,j≤m be a sequence of random variables defined over the same probability space (Ω, A, P), that are independent and identically distributed. Each pixel value X i,j is associated to the cell C (m) i,j . In Figure 2, X i,j ∼ U(0, 1) and m = 20. Let us consider a threshold parameter t ∈ R. In order to create the associated binary image, we introduce a random m × m-matrix Z(t) = (Z i,j (t)) 1≤i,j≤m , where

Z (m)
i,j (t) := 1 {Xi,j ≥t} , for i, j ∈ {1, . . . , m}.

Each cell C

(m) i,j is associated to black or white according to whether Z i,j (t) = 0 or Z i,j (t) = 1. Then, Z i,j (t) follows a binomial distribution of parameter (1, p t ), where

p t := P(Z i,j (t) = 1) = P(X i,j ≥ t) = 1 -F (t -),
with F the cumulative distribution function associated to X i,j .

Perimeter of a binary image

Let Z = Z(t) be the binary image Z at that given threshold t ∈ R. Following the approach presented in [START_REF] Biermé | The effect of discretization on the mean geometry of a 2D random field[END_REF], for each edge w ∈ E m , we aim to know if w contributes to the perimeter of the black component of Z. Let w be a horizontal edge of the form w = l-1 m , l m × { k-1 m }, w belongs to both cells C l,k-1 and C l,k , which means that w is a contribution to the perimeter if Z l-1,k = Z l,k . Following this consideration, one can consider

f (t) 1 (l, k) = 1 {(Z l,k-1 (t) = 0 ∩ Z l,k (t) = 1) ∪ (Z l,k-1 (t) = 1 ∩ Z l,k (t) = 0)}
to count the horizontal contributions and

f (t) 2 (k, l) = 1 {(Z k-1,l (t) = 0 ∩ Z k,l (t) = 1) ∪ (Z k-1,l (t) = 1 ∩ Z k,l (t) = 0)}
to count the vertical ones. Definition 2.1 (Perimeter and scaled perimeter of a binary image). We denote by P (1) m (t) = m l=1 m k=2 f (t) 1 (l, k) the sum of all horizontal contributions and P (2)

m (t) = m l=1 m k=2 f (t)
2 (k, l) the sum of the vertical ones, then, the perimeter and scaled perimeter are given by P m (t) := P (1) m (t) + P (2) m (t), Pm (t) := 1 m 2 (P (1) m (t) + P (2) m (t)) .

An equivalent way to compute the perimeter is by considering, for each edge w, the maximal and minimal values on the two sides of w. Then, the perimeter is equivalently given by

P m (t) := w∈Em (f (t) + (w) -f (t) -(w)), (2)
where f (t) + (w) = max(Z l,k-1 (t), Z l,k (t)) and f (t) -(w) = min(Z l,k-1 (t), Z l,k (t)) for w the common edge between C l,k-1 and C l,k . The interested reader is referred to [START_REF] Biermé | The effect of discretization on the mean geometry of a 2D random field[END_REF]. The only edges w ∈ E m that contribute to the computation are those for which f (t) + (w) = 1 and f (t) -(w) = 0 which are exactly the edges belonging to the perimeter intersected with the interior of our observation window.

Algorithm 1 Perimeter of a binary image computation using Equation (2)

Initialization Binary image Z of shape (M, N ),

Do L 1 = sum((Z[0 : M -1, 0 : N ] == 1) × (Z[1 : M, 0 : N ] == 0)) L 2 = sum((Z[0 : M, 0 : N -1] == 1) × (Z[0 : M, 1 : N ] == 0)) L 3 = sum((Z[0 : M -1, 0 : N ] == 0) × (Z[1 : M, 0 : N ] == 1)) L 4 = sum((Z[0 : M, 0 : N -1] == 0) × (Z[0 : M, 1 : N ] == 1)) P 1 = L 1 + L 3 P 2 = L 2 + L 4 Return P 1 + P 2
Remark 1. Notice that the proposed method in Definition 2.1 (see Equation (2)) does not take into account the contribution of the edges that belongs to the frontier of the observation window S. Indeed, in Algorithm 1, the border cells can only contribute to the perimeter with one or two edges contrary to the other cells that can have up to four edges contributing to the perimeter. In this sense this method is unbiased in comparison with [START_REF] Ebner | Goodness-of-fit tests for complete spatial randomness based on Minkowski functionals of binary images[END_REF]'s technique. In the aforementioned paper, the authors introduce a white artificial frontier around the image and then, compute the perimeter of the extended new image (white frontier included), making all cells contributing to the perimeter have a four edge contribution.

Statistics of the perimeter of a binary image

Moments and asymptotic normality

In this section we investigate the first and second moments of the perimeter. This preliminary study will be useful to state our multivariate Central Limit Theorem (see Theorem 3.3 below). Proposition 3.1 (First moment of the perimeter). The expected value of the scaled perimeter in (1) is given by

E( Pm (t)) = 4p t (1 -p t ) 1 - 1 m := µ P (p t , m). (3)
The proof of Proposition 3.1 is postponed to Section 7.

Remark 2. Assume that X i,j is drawn from a symmetrical continuous distribution of axis θ ∈ R, i.e p θ-t = 1 -p θ+t , ∀t ∈ R. Then, µ P (p θ-t , m) = µ P (p θ+t , m) and therefore, t → µ P (p t , m) has an axis of symmetry on θ. This implies that the ratio ( Pm (θ-t))/( Pm (θ+t)) would be distributed around 1 in the case of symmetry. This consideration will be crucial for the proposed symmetry testing procedure in Section 4.

We now study the second moment of the perimeter of the considered binary image.

Proposition 3.2 (Covariance between perimeters at two thresholds). Let t, s ∈ R be two given thresholds. The covariance between the scaled perimeter taken at levels t and s is given by

Cov Pm (t), Pm (s) = 4 m 4 p max(t,s) (1 -p min(t,s) ) (4m 2 -7m + 2) -2p min(t,s) (1 -p max(t,s) )(7m 2 -13m + 4) + 2(p min(t,s) -p max(t,s) )(3m 2 -6m + 2) .
Taking t = s, we obtain the following variance formula

Var( Pm (t)) = 4 m 4 p t (1 -p t ) (4m 2 -7m + 2) -2p t (1 -p t )(7m 2 -13m + 4) := σ 2 P (t, m). (4)
The proof of Proposition 3.2 is postponed to Section 7. If we add in Proposition 3.2 the bias induced by the boundary of the window S to the variance computation, we obtain the same results as those provided in [START_REF] Ebner | Goodness-of-fit tests for complete spatial randomness based on Minkowski functionals of binary images[END_REF] (see Theorem 4.1). Beyond Propositions 3.1 and 3.2, we can prove the multivariate asymptotic normality as m → ∞ of our geometrical feature for a given vector of thresholds.

Theorem 3.3 (Multivariate CLT for r-thresholds). Let r be a positive integer, m ≥ 2 and t 1 , . . . , t r ∈ R, then, m           Pm (t 1 ) Pm (t 2 ) . . . Pm (t r )      -      E( Pm (t 1 )) E( Pm (t 2 )) . . . E( Pm (t r ))           d ----→ m→∞ N (0, Σ r ) ,
where d -→ holds for the convergence in distribution and N (0, Σ r ) for the rdimensional centered Gaussian distribution with covariance matrix Σ r given by Σ r (i, j) := 4p max(ti,tj ) (1 -p min(ti,tj ) ) 4 -14p min(ti,tj ) (1 -p max(ti,tj ) ) + 6(p min(ti,tj ) -p max(ti,tj ) ) .

(

The proof of Theorem 3.3 is postponed to Section 7. Furthermore, some numerical illustrations of Theorem 3.3 in the case r = 2 can be found in Appendix A.

Finally, as can be seen in ( 5), the last theoretical result that will be needed for the construction of an empirical accessible test statistic is an estimation of p t . This will be discussed in the next section.

3.2. Two consistent estimators for p t .

In the following, we will present two estimators and study their consistency properties. In this section, for the sake of clarity, we denote p(t) := p t , ∀ t ∈ R.

One way to estimate p(t) is by considering the area of the excursion set. Another idea for the estimation of p(t) is built on the relationship that exists between the expectation value of the perimeter at a level t and p(t) (see Equation ( 3)):

E Pm (t) (1 -1 m ) = 4p(t)(1 -p(t)). (6)
Definition 3.4 (Area and perimeter based estimator for p(t)). Given a threshold t ∈ R, we define

p A m (t) = 1 m 2 m i,j=1 Z i,j (t).
Besides, we partition the image S in m 2 /4 sub-images S i 2 of size (2 × 2). We denote by Pi 2 (t) the value of the scaled perimeter for each sub-image. Let S m (t) := 4

m 2 i
Pi 2 (t) and the continuous mapping g : [0, 1] → R defined by

x → g(x) =      1 2 1 - √ 1 -2x if x < 1 2 , 1 2 else.
Then, for t > θ with θ the median value of the distribution, we define the perimeter based estimator of p(t) as

p P m (t) := g(S m (t)).

Proposition 3.5 (Asymptotic normality for the proposed estimators for p(t)).

Let p A m (t) and p P m (t) as in Definition 3.4. Then, it holds that

p A m (t) a.s ----→ m→∞ p(t) and m ( p A m (t) -p(t)) d ----→ m→∞ N (0, p(t)(1 -p(t))).
Similarly, it holds

p P m (t) a.s ----→ m→∞ p(t)
and

m ( p P m (t) -p(t)) d ----→ m→∞ N   0, p(t)(1 -p(t)) 1 -3p(t)(1 -p(t)) 4(1 -2p(t)) 2   .
The proof of Proposition 3.5 is postponed to Section 7.

Remark 3. It is also possible to estimate p(t) for t < θ by considering the con-

tinuous mapping h : [0, 1] → R, x → h(x) =      1 2 1 + √ 1 -2x if x < 1 2 , 1 2
else, and applying the same procedure.

Perimeter based symmetry test

As previously discussed in Remark 2, we now aim to build a test using the specific structure of the mean perimeter in the case of symmetry.

Let I = x ∈ R, 0 < p(x) < 1 2 and I θ = I -θ, it follows that ∀ t ∈ I θ , 0 < p(θ + t) < 1 2 .
A test statistic Following the idea in [START_REF] Biermé | Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields[END_REF]

Pm (t 2 ) - µ Pm (t 1 ) µ Pm (t 2 ) d ----→ m→∞ N 0, σ 2 (t 1 , t 2 ) , with σ 2 (t 1 , t 2 ) = p t1 (1 -p t1 ) (p t2 (1 -p t2 )) 3 p t2 (1 -p t2 ) -2p t2 (1 -p t1 ) -3p t2 (p t1 -p t2 ) (1 -p t1 ) + p t1 (1 -p t1 ) . Let t ∈ I θ and R m,θ (t) := Pm (θ -t) Pm (θ + t) . Then, under H 0 (t), it holds that m (R m,θ (t) -1) d,H0 ----→ m→∞ N 0, σ 2 (θ + t) , (7)
with σ 2 (θ + t) = (2p θ+t -1)(3p θ+t -2) p θ+t (1 -p θ+t ) 2 .
We postpone the proof of Proposition 4.1 in Section 7.

Proposed test with asymptotic level α

Let α ∈ (0, 1) and q 1-α/2 such that P(N (0, 1) ≤ q 1-α/2 ) = 1 -α/2. Then, using Proposition 4.1, we define the statistic of the test with asymptotic level α,

φ P m (σ) = 1 | m σ(θ+t) (R m,θ (t) -1)| ≥ q 1-α/2 . ( 8
)
Firstly, we study the consistency of the proposed test statistic (see Proposition 4.2 below). Secondly, we focus on some estimation procedures for σ(θ + t) (see Section 4.2).

Let t ∈ I θ and H 1 (t) be the alternative hypothesis

H 1 (t) : p θ-t = 1 -p θ+t .
Proposition 4.2 (Consistency of the proposed symmetry test). For t ∈ I θ , it holds that P H1(t) (φ P m (σ) = 1) → 1 for m → ∞, with φ P m (σ) as in Equation (8). We postpone the proof of Proposition 4.2 in Section 7.

Empirical accessible test statistic

Notice that, if we take a consistent empirical estimator σ m (θ + t) of the variance of the considered ratio at level θ+t, then, we can consider the empirical accessible test statistic

φ P m ( σ) = 1 | m σm(θ+t) (R m,θ (t) -1)| ≥ q 1-α/2 . (9)
Using the estimators of p(t) presented in Proposition 3.5, we can now construct estimators of the variance of the ratio to built φ P m ( σ) in ( 9). Proposition 4.3 (Plug-in variance estimators). Let p m (θ + t) be a consistent estimator of p(θ+t) (i.e p m (θ + t) → p(θ + t), a.s.) and σ(θ+t) as in Proposition 4.1. We define

σ m (θ + t) : = (2 p m (θ + t) -1)(3 p m (θ + t) -2) p m (θ + t)(1 -p m (θ + t)) 2 . ( 10
)
Then, σ m (θ + t) a.s.

----→ m→∞ σ(θ + t) and

m σ m (θ + t) (R m,θ (t) -1) d,H0 ----→ m→∞ N (0, 1). ( 11
)
The proof of Proposition 4.3 is postponed to Section 7.

Corollary 4.4. Using the same framework as in Proposition 4.3, we denote by σ P m (θ + t) and σ A m (θ + t) the plug-in estimators based on Equation (10) using respectively p P m (θ + t) and p A m (θ + t) as in Proposition 3.5. Then, σ P m (θ + t) and σ A m (θ + t) are two consistent estimators of σ(θ + t) and asymptotic normality in [START_REF] Desolneux | Stochastic methods for image analysis[END_REF] holds for these two estimators.

Thus, from Corollary 4.4 and using Equation ( 9), we get two empirical accessible test statistics

φ P m ( σ P ) = 1 | m σ P m (θ+t) (R m,θ (t) -1)| ≥ q 1-α/2 (12) and φ P m ( σ A ) = 1 | m σ A m (θ+t) (R m,θ (t) -1)| ≥ q 1-α/2 . ( 13
)

Comparison with Kolmogorov-Smirnov based estimator

Remark that under H 0 (t), the cumulative distribution function F satisfies F (θ -t) = 1 -F (θ + t). Then, considering the empirical cumulative

distribution function t → F m (t) := 1
m 2 i,j
m 1 -p A m (θ -t) p A m (θ + t) -1 d,H0 ----→ m→∞ N 0, (σ k-s ) 2 (θ + t) , (14)
with (σ k-s ) 2 (θ + t) = 2 p(θ+t) . The proof of Proposition 4.5 is postponed to Section 7. Let us denote

σ k-s m (θ + t) = 2 p A m (θ + t) and R k-s m,θ (t) := 1 -p A m (θ -t) p A m (θ + t) .
One can provide an alternative Kolmogorov-Smirnov based test statistic with asymptotic level α

φ k-s m ( σ k-s ) = 1 | m σ k-s m (θ+t) (R k-s m,θ (t) -1)| ≥ q 1-α/2 . ( 15
) Remark 4. In particular for s ∈ I we can prove that σ 2 (s) < σ 2 k-s (s)
, where σ 2 (s) is given by ( 7) and σ 2 k-s (s) by ( 14), as shown in Figure 9 in Appendix B. Then, the perimeter based test has a smaller variance than the alternative Kolmogorov-Smirnov one.

Some numerical studies for the proposed symmetry test

In this section, we report simulation results for our symmetry test, for samples of 500 images of size m × m, for m = 100, 512, 1024 and different choices of threshold level. Firstly, we display the results for three symmetrical distributions (Gaussian, Student and Uniform distributions, see Figure 4) and two asymmetrical ones (Exponential and Skew-normal distributions, see Figure 5).

For a better approximation of θ, we use the mean estimator provided by the NumPy library, making use of the fact that for a symmetrical distribution it is equivalent to use the mean or median to characterize the symmetry of a distribu-tion [i.e. a random variable X is symmetrically distributed around θ, the center of symmetry, if its distribution function

F satisfies F (θ -x) = 1 -F (θ + x)].
As we can appreciate in Figure 4, the test is able to successfully accept the H 0 hypothesis. However, the choice of the threshold is crucial in regard to the quality of the test. Unsurprisingly, for extreme thresholds, the test is less precise. We will focus on this important point in Section 6 below. Notice that the performances of our test statistic φ P m with the two proposed estimators of the variance σ A and σ P (in red crosses and blue diamonds, respectively) seem to be globally similar. Furthermore, these satisfactory results are comparable to those obtained with the Kolmogorov-Smirnov based estimator in (15) (black stars). m = 1) for different thresholds t, for 500 Montecarlo simulations for m = 512 with a skew-normal marginal distribution with different choices of parameter a = 0.1, 0.5, 1, 3 (from first to fourth panel) and Exp(1) distribution (fifth panel). We consider different test statistics: φ P m ( σ P ) as in (12) (blue diamonds), φ P m ( σ A ) as in (13) (red crosses) and φ k-s m ( σ k-s ) as in (15) (black stars).

Skew-normal and Exponential marginal distributions

Figure 5 numerically describes the good performance of our test φ P m to distinguish H 1 (t) from H 0 (t), for t ∈ R + , in two asymmetric distributional cases: Skew-normal distribution with probability density function given by f (x) = 2φ(x)Φ(ax), for x ∈ R, with a ≥ 0, φ(•) (resp. Φ(•)) the standard normal probability density function (resp. standard normal cumulative distribution function) and Exponential distribution with parameter 1. Unsurprising, in the skew-normal case, we observe that when a is larger the power of the test increases, indeed the smaller the parameter a, the less asymmetric the distribution will be. This behavior can be easily observed for the four chosen levels of skewness parameter a = 0.1, 0.5, 1, 3 in Figure 5 (from first to fourth panel).

Perimeter based symmetry test for large threshold and dense tiling

Following the considerations in Section 5, we aim to improve the quality of the proposed test for large thresholds. As in [START_REF] Dalmao | Asymptotic normality of high levellarge time crossings of a Gaussian process[END_REF], we firstly study the convergence of the ratio between the first two moments of the unscaled perimeter P m (see Equations ( 3) and ( 4)). In Lemma 6.1 below, the expectation and the variance are proven to have the same order of magnitude for large threshold and dense tiling.

Lemma 6.1 (First two moments ratio for large threshold). Let P m (t) as in (1) and µ P (p t , m) and σ 2 P (p t , m) the associated mean and variance. Let (t m ) be a sequence of positive real numbers such that p tm → 0 as m → ∞, then,

σ 2 P (p tm , m) µ P (p tm , m) ----→ m→∞ 4,
and furthermore, for m large enough,

0.5 ≤ σ 2 P (p tm , m) µ P (p tm , m) ≤ 4.
The proof of Lemma 6.1 is postponed to Section 7. By using Lemma 6.1 we can now formulate a bivariate Central Limit Theorem for our geometrical feature for m → ∞ and large thresholds.

Proposition 6.2 (Bivariate Central Limit Theorem for large threshold and dense tiling). Let 0 < γ < 2 and t (m,γ) and s (m,γ) two points of continuity for the cumulative distribution function F , such that

t (m,γ) := F -1 ( 1 m γ ) and s (m,γ) := F -1 (1 -1 m γ ). Then, it holds that m 2 E( Pm (t (m,γ) )) 4           Pm (t (m,γ) ) E( Pm (t (m,γ) )) Pm (s (m,γ) ) E( Pm (s (m,γ) ))      - 1 1      d ----→ m→∞ N 0 0 , 1 0 0 1 . ( 16
)
The proof of Proposition 6.2 is postponed to Section 7. Note that we have written E Pm (t (m,γ) ) and E Pm (s (m,γ) ) in ( 16) for visual symmetry purposes although those two quantities are equal (see Equation ( 3)). Finally, we can easily derive the following result.

Corollary 6.3. Under assumptions of Proposition 6.2, it holds that

m 2 E( P(t (m,γ) )) 8 P(t (m,γ) ) P(s (m,γ) ) -1 d ----→ m→∞ N (0, 1). ( 17
)
Furthermore,

m 2 Pm (t (m,γ) ) 8 P(t (m,γ) ) P(s (m,γ) ) -1 d ----→ m→∞ N (0, 1). ( 18
)
Equation ( 17) is obtained via the Delta method applied to Proposition 6.2. Finally, Equation ( 18) is given by using that, for m → ∞, Pm (t (m,γ) )/E( Pm (t (m,γ) ))

P -→ 1 and by applying the Slutsky Theorem.

Using Corollary 6.3, one can build a modified version of the empirical accessible test statistic φ P m in ( 9) with asymptotic level α, adapted to large thresholds. Under H 0 (t) hypothesis, let θ be the median value of the distribution and u (m,γ) ∈ R + so that θ +u (m,γ) := F -1 (1-1 m γ ), and θ -u (m,γ) := F -1 ( 1 m γ ). Coming back to the ratio R m,θ introduced in Section 4, we have

R m,θ (u (m,γ) ) := Pm (θ -u (m,γ) ) Pm (θ + u (m,γ) ) .
Then, from Corollary 6.3, it holds that

v (m,γ) (R m,θ (u (m,γ) ) -1) d,H0 ----→ m→∞ N (0, 1) with v (m,γ) = m 2 P(θ -u (m,γ) ) 8 .
This asymptotic allows us to define a new symmetry test statistic by prescribing

φ P m ( v (m,γ) ) = 1 v (m,γ) |R m,θ (u (m,γ) )-1|≥q 1-α/2 . ( 19
)
Conversely to Equation ( 9), the main statistical interest of the test in ( 19) is to avoid the estimation of the variance. Furthermore, as one can remark for instance in the numerical study in Figure 6, global performances of the test in [START_REF] Heinrich | Stable limit theorems for sums of multiply indexed mdependent random variables[END_REF] are similar to those obtained in the same distributional framework in Figure 4. This means that we can obtain similar p-values with less effort in terms of moment estimations. Furthermore, notice that, since γ = -ln(pu m,γ ) ln(m)

and γ ∈ (0, 2), the range of admissible values for u m,γ from Proposition 6.2 is growing with m.

In order to better appreciate the statistical equivalent performance of the four proposed test statistics, we gathered in Table 1 averaged values of the empirical P H0(t) (φ P m = 1) and associated standard deviations for different thresholds and several distributional models. We underline in bold some values for which the last proposed test in [START_REF] Heinrich | Stable limit theorems for sums of multiply indexed mdependent random variables[END_REF] seems slightly outperform the other ones especially for large threshold values (see Table 1). Thresholds tm,γ 0.7 0.78 0.85 0.93 b P H 0 (φ P m ( σ P ) = 1) 0.002 (0.045) 0.0 (0.0) 0.004 (0.063) 0.094 (0.292) b P H 0 (φ A m ( σ A ) = 1) 0.002 (0.045) 0.0 (0.0) 0.004 (0.063) 0.094 (0.292) bblP H 0 (φ k-s m ( σ k-s ) = 1) 0.00 (0.0) 0.0 (0.0) 0.004 (0.063) 0.096 (0.295) P H 0 (φ P m ( v) = 1) 0.00 (0.0) 0.0 (0.0) 0.00 (0.0) 0.064 (0.245)

Table 1

Under H 0 hypothesis. Average of the empirical values of P H 0 (t) (φ P m = 1) and associated standard deviations for different thresholds tm,γ , for 500 Montecarlo simulations, for m = 1024. We consider different test statistics: φ P m ( σ P ) as in [START_REF] Bernardino | Statistics for Gaussian random fields with unknown location and scale using Lipschitz-Killing curvatures[END_REF], φ P m ( σ A ) as in (13), φ k-s m ( σ k-s ) as in [START_REF] Estrade | A central limit theorem for the Euler characteristic of a Gaussian excursion set[END_REF] and finally φ P m ( v (m,γ) ) as in [START_REF] Heinrich | Stable limit theorems for sums of multiply indexed mdependent random variables[END_REF].

Proofs

Proof of Proposition 3.1. One can start by observing that both f (t) 1 (l, k) and f (t) 2 (k, l) follow a Bernoulli distribution with parameter 2p t (1 -p t ). Thus, using Equation ( 1), one can write

E(P m (t)) = E(P (1) m (t) + P (2) m (t)) = 2 m l=1 m k=2 E(f (t) 1 (l, k)) = 4m(m -1)p t (1 -p t ).
Thus we obtain the result for E Pm dividing by m 2 .

Proof of Proposition 3.2. Without loss of generality, one can assume that t ≤ s.

We use here Definition 1 of the perimeter. Let us start by making use of the fact that (P (1) m (s), P (1) m (t)) d = (P (2) m (s), P (2) m (t)). Thus, Cov (P m (s), P m (t)) = 2Cov (P (1) m (s), P (1) m (t)) + 2Cov (P (2) m (s), P (1) m (t)) .

Firstly,

Cov (P (1) m (s), P (1) m (t)) = Cov   m l=1 m k=2 f (s) 1 (l, k), m i=1 m j=2 f (t) 1 (i, j)   = m l=1 m k=2 m j=2 Cov f (s) 1 (l, k), f (t) 1 (i, j) ,
where we have used that for ∀k, l,

, j) = 0. In the above sum, there will be only three cases for which the covariance function is different from 0: j = k, j = k + 1 and j = k -1.

1. If j = k, E(f (s) 1 (l, k)f (t) 1 (l, k)) = P f (s) 1 (l, k) = 1 ∩ f (t) 1 (l, k) = 1 = P ((X l,k-1 < s ∩ X l,k ≥ s) ∪ (X l,k-1 ≥ s ∩ X l,k < s) ∩ ((X l,k-1 < t ∩ X l,k ≥ t) ∪ (X l,k-1 ≥ t ∩ X l,k < t)) = P X l,k-1 < s ∩ X l,k ≥ s ∩ X l,k-1 < t ∩ X l,k ≥ t ∪ X l,k-1 < s ∩ X l,k ≥ s ∩ X l,k < t =∅ ∩ X l,k-1 ≥ t ∪ X l,k-1 ≥ s ∩ X l,k-1 < t =∅ ∩ X l,k < s ∩ X l,k ≥ t ∪ (X l,k-1 ≥ s ∩ X l,k < s ∩ X l,k-1 ≥ t ∩ X l,k < t) = P X l,k-1 < t ∩ X l,k ≥ s + P X l,k-1 ≥ s ∩ X l,k < t = p s (1 -p t ) + (1 -p t )p s = 2p s (1 -p t ). Thus, Cov f (s) 1 (l, k), f (t) 1 (l, k) = 2p s (1 -p t ) -2p s (1 -p s )2p t (1 -p t ) = 2p s (1 -p t )(1 -2p t + 2p t p s ). 2. If j = k + 1, E(f (s) 1 (l, k)f (t) 1 (l, k + 1)) = P f (s) 1 (l, k) = 1 ∩ f (t) 1 (l, k + 1) = 1 = P ((X l,k-1 < s ∩ X l,k ≥ s) ∪ (X l,k-1 ≥ s ∩ X l,k < s) ∩ ((X l,k < t ∩ X l,k+1 ≥ t) ∪ (X l,k ≥ t ∩ X l,k+1 < t))) = P (X l,k-1 < s ∩ X l,k ≥ s ∩ X l,k < t =∅ ∩ X l,k+1 ≥ t) ∪ (X l,k-1 < s ∩ X l,k ≥ s ∩ X l,k ≥ t ∩ X l,k+1 < t) ∪ (X l,k-1 ≥ s ∩ X l,k < s ∩ X l,k < t ∩ X l,k+1 ≥ t) ∪ (X l,k-1 ≥ s ∩ X l,k < s ∩ X l,k ≥ t ∩ X l,k+1 < t) = P(X l,k-1 < s)P(X l,k ≥ s)P(X l,k+1 < t) + P(X l,k-1 ≥ s)P(X l,k < t)P(X l,k+1 ≥ t) + P(X l,k-1 ≥ s)P(X l,k ≥ t ∩ X l,k < s)P(X l,k+1 < t) = (1 -p s )p s (1 -p t ) + (1 -p t )p t p s + p s (p t -p s )(1 -p t ) = p s (1 -p t )(1 + 2p t -2p s ). Thus, Cov f (s) 1 (l, k), f (t) 1 (l, k + 1) = p s (1 -p t )(1 + 2p t -2p s ) -2p s (1 -p s )2p s (1 -p t ) = p s (1 -p t )(1 -2p t -2p s + 4p t p s ). 3. Similarly, if j = k -1, Cov(f (s) 1 (l, k), f (t) 1 (l, k -1)) = Cov(f (s) 1 (l, k), f (t) 1 (l, k + 1)) = p s (1 -p t )(1 -2p t -2p s + 4p t p s ).
Thus, Cov P (1) m (s), P (1)

m (t) = m l=1 m k=2 Cov f (s) 1 (l, k), f (t) 1 (l, k) + m l=1 m k=3 Cov f (s) 1 (l, k), f (t) 1 (l, k + 1) + m l=1 m k=3 Cov f (s) 1 (l, k), f (t) 1 (l, k -1) = 2m(m -1)p s (1 -p t )(1 -2p t (1 -p s )) + 2m(m -2)p s (1 -p t )(1 + 2p t -2p s -4p t (1 -p s )).
Now, let us compute the inter-level covariance between the horizontal contributions to the perimeter and the vertical ones:

Cov (P (1) m (s), P

2 (j, i) .

To this aim, notice that the covariance between cells is non equal to zero only if the cells are neighbors, that implies that

for i / ∈ {k -1, k} or for l / ∈ {j -1, j} Cov(f 1 (l, k), f 2 (j, i)) = 0.
Thus,

Cov(P 1 m (s), P 2 m (t)) = m k=2 m j=2
Cov(f

(s) 1 (j -1, k), f (t) 2 (j, k -1)) + Cov(f (s) 1 (j, k), f (t) 2 (j, k -1)) + Cov(f (s) 1 (j -1, k), f (t) 2 (j -1, k)) + Cov(f (s) 1 (j, k), f (t) 2 (j, k)) .
Due to the i.i.d setting (see Figure 7), one can remark that the four types of covariance in the above sum are equal and are equal to Cov(f

(s) 1 (l, k), f (t) 2 (l, k + 1)). ∼ Fig 7.
Representation of the configuration between (f 1 (l, k), f 2 (l + 1, k -1)) (left panel) and (f 1 (l, k), f 1 (l, k + 1)) (right panel), the black dot represents the cell (k, l).

Thus,

Cov(P 1 m (s), P 2 m (t)) = 4(m -1) 2 Cov(f (s) 1 (l, k), f (t) 2 (l, k + 1)) = 4(m -1) 2 p s (1 -p t )(1 + 2p t -2p s -4p t (1 -p s )).
Putting all elements together, we get the covariance function of the unscaled perimeter: Cov (P m (s), P m (t)) = 4 p s (1 -p t )((4m 2 -7m + 2)

-2p t (1 -p s )(7m 2 -13m + 4) + 2(p t -p s )(3m 2 -6m + 2)).
Hence, the result.

Proof of Theorem 3.3. In order to prove this result, we will use the well known Cramèr-Wold method. For the sake of completeness we recall it below.

Theorem 7.1 (Cramèr-Wold method, see, e.g., [START_REF] Billingsley | Probability and Measure[END_REF]

, p. 383). Let X m = (X m 1 , • • • , X m r ) and X = (X 1 , • • • , X r ) be random vectors of dimension r ≥ 1. Then, X m con- verges in distribution to X if and only if r i=1 a i X m i d -----→ m→+∞ r i=1 a i X i ,
for each (a 1 , • • • , a r ) ∈ R r , i.e., if every fixed linear combination of the coordinates of X m converges in distribution to the correspondent linear combination of coordinates of X.

We first need to introduce some notations. Let V m = {1, . . . , m} × {2, . . . , m} a 1 , . . . , a r ∈ R, t 1 , . . . , t r ∈ R, (l, k) ∈ V m and i ∈ {1, . . . , r}. Let P m (t i ) as in Equation ( 1). We introduce

Y l,k (t i ) = 1 (Z (t i ) l,k-1 = 0 ∩ Z (t i ) l,k = 1) ∪ (Z (t i ) l,k-1 = 1 ∩ Z (t i ) l,k = 0) + 1 (Z (t i ) k-1,l = 0 ∩ Z (t i ) k,l = 1) ∪ (Z (t i ) k-1,l = 1 ∩ Z (t i ) k,l = 0) . Note that 0 ≤ Y l,k (t i ) ≤ 2 a.s. and |Y l,k (t i ) -E(Y l,k (t i ))| ≤ 2, thus | r i=1 a i (Y k,l (t i ) -E(Y l,k (t i )))| ≤ 2 r i=1 |a i |. We note W m l,k = 1 m r i=1 a i Y l,k (t i ) -E(Y l,k (t i ))) , (20)
and finally,

S m = m l=1 m k=2 W m l,k = a 1 m (P m (t 1 ) -E(P m (t 1 ))) + . . . + a r m (P m (t r ) -E(P m (t r ))).
Using the Cramèr-Wold method, we need to prove that

S m d ----→ m→∞ N   0, r i=1 a 2 i Σ r (i, i) + 0≤i<j≤r 2a i a j Σ r (i, j)   ,
with Σ r as in [START_REF] Billingsley | Probability and Measure[END_REF]. To this purpose we use the following weaker version of Theorem 1 of [START_REF] Heinrich | Stable limit theorems for sums of multiply indexed mdependent random variables[END_REF] for bounded variables (see p.200 in [START_REF] Heinrich | Stable limit theorems for sums of multiply indexed mdependent random variables[END_REF] for details).

Proposition 7.2. Let V m = {1, . . . , m} × {2, . . . , m} and {X m z , z ∈ V m } be a sequence of 1-dependent random field, and S m = z∈Vm X m z . If the following three conditions are satisfied:

(i) ∀ > 0, z∈Vm P(|X m z | ≥ ) ----→ m→∞ 0, (ii) z∈Vm E((X m z ) 2 ) ≤ C < +∞, (iii) E(S m ) ----→ m→∞ 0 and Var(S m ) ----→ m→∞ σ 2 with σ 2 ∈ (0, ∞), then, S m d ----→ m→∞ N (0, σ 2 ).
We start by showing that {W m z , z ∈ V m } satisfies (i). For > 0, we get that,

(l,k)∈Vm P W m l,k ≥ ≤ (l,k)∈Vm E( W m l,k 3 ) 3 ≤ (2 r i=1 |a i |) 3 × m(m -1) (m ) 3 ----→ m→∞ 0.
To prove that {W m z , z ∈ V m } satisfies (ii), we write

(l,k)∈Vm E W m l,k 2 ≤ m(m -1) 2 m r i=1 |a i | 2 ≤ 2 r i=1 |a i | 2 .
For (iii) we observe that Var(S m ) = Var a1 m P m (t 1 ) + . . . + ar m P m (t r ) thus, as m → ∞, Var a1 m P m (t 1 ) + . . . + ar m P m (t r ) → r i=1 a 2 i Σ r (i, i) + 0≤i<j≤r 2a i a j Σ r (i, j), by Equation [START_REF] Billingsley | Probability and Measure[END_REF]. Using the auxiliary Proposition 7.2 proves Theorem 3.3.

Proof of Proposition 3.5. Let θ be the median value of the underlying distribution. We first start with solving the quadratic Equation ( 6) for m = 2 under the condition that t > θ (i.e p(t) < 1 2 ). We get that p(t) =

m 2 4 S m (t) -E P2 (t) d ----→ m→∞ , N (0, σ 2 P (t, 2)), (21)
with σ 2 P (t, 2) = Var P2 (t) = p t (1 -p t ) (1 -3p t (1 -p t )). We apply the Delta method to Equation (21) using the function g as defined in Proposition 3.4. Note that g is differentiable on (0, 1 2 ) and that for t

= 1 2 , E P2 (t) < 1 2 Thus m 2 4 g(S m (t)) -g(E P2 (t) ) d ----→ m→∞ g (E P2 (t) N (0, σ 2 P (t, 2)), with g (x) = 1 2 √ 1 -2x . Since 1 -2E P2 (t) = 1 -2p t , we get m ( p P m (t) -p(t)) d ----→ m→∞ N   0, p(t)(1 -p(t)) 1 -3p(t)(1 -p(t)) 4(1 -2p(t)) 2   . + P H1(t) m σ R m,θ (t) - µ Pm (θ -t) µ Pm (θ + t) - m σ δ ≤ -q 1-α/2 σ σ . Then P H1(t) (φ P m (σ) ----→ m→∞ 1,
hence, we get the result.

Proof of Corollary 4.3. Considering the continuous mapping

ψ : x → (2x -1)(3x -2) x(1 -x) 2 -1 2
and the fact that p m (θ +t) is a robust estimator of p(θ + t), we get that ψ( p m (θ + t))

a.s.

----→ m→∞ ψ(p(θ + t)) where

ψ(p(θ + t)) = 1 σ(θ + t)
. Thus, applying the Slutsky Theorem to (7), we get

mψ( p m (θ + t))(R m,θ (t) -1) d,H0 ----→ m→∞ ψ(p(θ + t))N (0, σ 2 (θ + t)).
Hence the result.

Proof of Proposition 4.5. For the sake of simplicity, let us note θ -t = u and

θ + t = s, for t ∈ (0, ∞) and µ = 1 -p u 1 -p s . We recall that, ∀x ∈ R F m (x -) = 1 m 2 m i=1 m j=1
1 {Xi,j <x} .

Furthermore, from the bivariate Central Limit Theorem,

m F m (u -) F m (s -) - F (u -) F (s -) d ----→ m→∞ N 0 0 , ( 1
-p u )p u (1 -p u )p s (1 -p u )p s (1 -p s )p s .
Applying the Delta method using the function g : where the second part of the inequality is due to [START_REF] Lachièze-Rey | Bicovariograms and Euler characteristic of random fields excursions[END_REF] and [START_REF] Lachièze-Rey | Normal convergence of nonlocalised geometric functionals and shot-noise excursions[END_REF].

x y → x 1 -y , it holds that m F m (u -) 1 -F m (s -) -1 d ----→ m→∞ N 0, ∇g(µ) T (1 -p u )p u (1 -p u )p s (1 -p u )p s (1 -p s )p s ∇g(µ) with ∇g(µ) = 1 ps 1-pu (ps) 2 . Since under H 0 (t), p u = 1 -p s , then, ∇g(µ) T (1 -p u )p u (1 -p u )p
Proof of Proposition 6.2. Let us first notice that E (P(t m,γ )) = E(P(s m,γ )), since p tm,γ = 1 -1 m γ = 1 -p sm,γ . In order to prove the result, we will establish that ∀ a 1 , a 2 ∈ R To prove this statement, we will follow the same idea as in the proof of Theorem 3.3, introducing the following notations, similar to (20), 2 ), we will check that the family { W m z , z ∈ V m } satisfies the three conditions of Proposition 7.2.

To show that { W m z , z ∈ V m } satisfies (i), we proceed as in the proof of Theorem 3.3. For > 0 and choosing p >

Fig 1 .

 1 Fig 1. Tilling with squares for m = 5 and associated vertices (red circles).

Fig 2 .

 2 Fig 2. Left and center panels: Image of size (20 × 20) realization of a Uniform white noise model. Right panel: Obtained binary image for t = 0.5.

Fig 3 .

 3 Fig 3. Computation of the perimeter of a binary image with m = 3. Here P 3 = 5.

Fig 4 .

 4 Fig 4. Under H 0 hypothesis. Average of the empirical values of P H 0 (t) (φ P m = 1) for different thresholds t, for 500 Montecarlo simulations and for m = 100 (left panels), 512 (center panels) and 1024 (right panels). We consider Gaussian marginal distribution N (0, 1) (first row), Student marginal distribution with ν = 3 degrees of freedom (second row), Uniform marginal distribution U (0, 1) (third row) and three different test statistics: φ P m ( σ P ) as in (12) (blue diamonds), φ P m ( σ A ) as in (13) (red crosses) and φ k-s m ( σ k-s ) as in (15) (black stars).

Fig 5 .

 5 Fig 5.Under H 1 hypothesis. Average of the empirical values of P H 1 (t) (φ P m = 1) for different thresholds t, for 500 Montecarlo simulations for m = 512 with a skew-normal marginal distribution with different choices of parameter a = 0.1, 0.5, 1, 3 (from first to fourth panel) and Exp(1) distribution (fifth panel). We consider different test statistics: φ P m ( σ P ) as in (12) (blue diamonds), φ P m ( σ A ) as in (13) (red crosses) and φ k-s m ( σ k-s ) as in (15) (black stars).

Fig 6 .

 6 Fig 6. Under H 0 hypothesis. Average of the empirical values of P H 0 (t) (φ P m = 1) for different thresholds tm,γ , for 500 Montecarlo simulations, for m = 100 (left panels), 512 (center panels) and 1024 (right panels). We consider a Gaussian marginal distribution N (0, 1) and different test statistics: φ P m ( σ A ) as in (13) (red stars), φ k-s m ( σ k-s ) as in (15) (black stars) and φ P m ( v (m,γ) ) as in (19) (green diamonds).

  t m,γ )) (a 1 (P (t m,γ ) -E (P (t m,γ ))) + a 2 (P (s m,γ ) -E (P (s m,γ ))))

  tm,γ)) a1 (Pm(tm,γ) -E(Pm(tm,γ)))+a2 (Pm(sm,γ) -E (Pm(sm,γ)))and finally, S m =

4 ×,≤

 4 4m(m -1) × 1 m γ (1 -1 m γ ) p Cm 2-1 2 p(2-γ) , for some positive constant C. As p > 4 2 -γ , Cm 2-1 2 p(2-γ) ----→ m→∞ 0.To show that { W m z , z ∈ V m } satisfies (ii),

  we propose a method to test the symmetry of the marginal distribution marginal of the field at equidistant levels from the median. Let us consider the null hypothesis for t ∈ I

θ H 0 (t) : p θ-t = 1 -p θ+t , where θ is the median value of the distribution. Let us first present a normality asymptotic result in the general case and a second result under the null hypothesis.

Proposition 4.1 (Ratio between perimeters at two thresholds). Let t 1 ∈ R and t 2 ∈ I with t 1 < t 2 , and Pm (•) be the scaled perimeter in Definition 2.1 computed for each threshold. Then, it holds that, m Pm (t 1 )

1

  Xi,j ≤t , which is a robust estimator of F , one can build an alternative symmetry test. Let us note that F m (t -) = 1 -p A m (t) with p A m (t) as in Definition 3.4. Proposition 4.5. Let t ∈ (0, ∞). Under H 0 (t), it holds that

  are i.i.d, and, by using the Law of Large numbers, S m (t) :=4

	1 2	1 -1 -2E P2 (t) .
	m 2 estimator of E P2 (t) . Applying the Central Limit Theorem, we get i	Pi 2 (t) is an

Since we are in an i.i.d. framework, the variables Pi 2 (t) for each sub-image S i

  s (1 -p u )p s (1 -p s )p s ∇g(µ) = 1Proof of Lemma 6.1. We recall from Equation (3) that we have, µ P (p t , m) = 4p t (1 -p t )m(m -1), and from Equation (4)σ 2 P (p t , m) = 4p t (1 -p t )((4m 2 -7m + 2) -2p t (1 -p t )(7m 2 -13m + 4)).Then, the ratio is equal to P (p t , m) µ P (p t , m) tends to 4 = lim m→∞ r m .For the second item, we can note that p t (1 -p t ) reaches its maximal value at p t = 1/2. It follows that

		σ 2 P (p t , m) µ P (p t , m)	=	(4 -7 m + 2 m 2 ) -2p t (1 -p t )(7 -13 m + 4 m 2 ) (1 -1 m )	.	(23)
	Let r m :=	(4 -7 m + 2 m 2 ) m 1 -1 , then,
	σ 2 P (p t , m) µ P (p t , m) ≤ |14p Hence, uniformly in m, -r m ≤ 2p t (1 -p t )(7m -13 + 4 m ) (m -1)
							σ 2 P (p t , m) µ P (p t , m)	---→ pt→0	r m .	(24)
	Thus, as m → ∞,	σ 2
		σ 2 P (p t , m) µ P (p t , m)	≥	(4m 2 -7m + 2) -1 2 (7m 2 -13m + 4) m(m -1)	=	1 2	.
	Hence, we have proven that for m large enough,
							1 2	≤	σ 2 P (p t , m) µ P (p t , m)	≤ 4,
								ps	1 ps	p u p s p 2 s p 2 s p u p s	1 ps 1 ps	=	2 p s	.
	Then,
			m			F m (u -) 1 -F m (s -)	-1	d,H0 ----→ m→∞	N 0,	2 p s	.
	Hence the result.

t (1 -p t )| .

  1 4E (P(t m,γ ))Var (a1 Y l,k (t m,γ ) + a 2 Y l,k (s m,γ )) Var (Y 1,2 (t m,γ )) + a 2 2 Var (Y 1,2 (s m,γ )) .Let us recall that Y 1,2 (t) is a sum of two Bernoulli random variables with parameter 2p t (1 -p t ). Hence, Var (Y 1,2 (t m,γ )) ≤ 2p tm,γ (1 -p tm,γ ). Since p tm,γ = 1 -p sm,γ , we have P m (t m,γ ) + a 2 P m (s m,γ ) Var(P(t m,γ )) + 2a 1 a 2 Cov(P(t m,γ ), P(s m,γ )) + a 2 2 Var(P(t m,γ )) 4p tm,γ (1 -p tm,γ ) Cov(P(t m,γ ), P(s m,γ )). Note that s m,γ > t m,γ , thus, by Proposition 3.2, Cov (P(t m,γ ), P(s m,γ)) = 4 p sm,γ (1 -p tm,γ )((4m 2 -7m + 2) + 2(p tm,γ -p sm,γ )(3m 2 -6m + 2)).Given that p tm,γ = 1 -p sm,γ , we haveCov (P(tm,γ), P(sm,γ)) = 4 (1 -pt m,γ ) 2 ((4m 2 -7m + 2) -2p 2 tm,γ (7m 2 -13m + 4) + 2(2pt m,γ -1)(3m 2 -6m + 2)).Recall that, p tm,γ → 1, thus, Cov (P(t m,γ ), P(s m,γ )) m(m -1)p tm,γ (1 -p tm,γ ) ----→

			m→∞	0. Hence, from
	Equation (25), we get that Var (S m ) ----→ m→∞	a 2 1 +a 2 2 , since	Var (P(t m,γ ) E (P(t m,γ )	----→ m→∞	4
	by Lemma 6.1.
	≤ 1 (l,k)∈Vm 2m(m -1) 4E (P (t m,γ )) a 2 E W m l,k 2 ≤ m(m -1) 4E (P (tm,γ)) (a 2 1 + a 2 2 ) × 8pt m,γ (1 -pt m,γ ) ≤ a 2 1 + a 2 2 .
	Thus, the sum is bounded. Finally, for (iii ),
	Var(S m ) = Var a 1 = 1 4E (P(t m,γ )) 1 4E (P(t m,γ )) a 2 1 = (a 2 1 + a 2 2 ) 4E(P(t m,γ )) 4m(m -1) Var(P(t m,γ )) + 1	2a 1 a 2
					(25)

Acknowledgments: This work has been partially supported by the project ANR MISTIC (ANR-19-CE40-0005) and by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.

Let us denote σ := σ(θ

the standard deviation of the ratio under H 0 (t) and σ the standard deviation in [START_REF] Kratz | Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields[END_REF]. Let φ P m (σ) as in Equation [START_REF] Csörgő | Testing for Symmetry[END_REF]. Under the alternative hypothesis H 1 (t), p θ-t = 1 -p θ+t . Besides, p θ-t = p θ+t .

Let define

. By hypothesis δ = 0. Recall that

Then, we get

In order to illustrate the multivariate CLT of Theorem 3.3, we follow the example given by [START_REF] Di Bernardino | Cross-Correlations and Joint Gaussianity in Multivariate Level Crossing Models[END_REF]. We start by simulating a sample of 900 images of size 1000 × 1000 pixels generated from a Uniform distribution and compute the perimeter without the border using the Algorithm 1 associated to the binary image for a threshold t = 0.5. We create the histogram of the variable:

σ 2 (t) [START_REF] Marinucci | Testing for non-Gaussianity on cosmic microwave background radiation: a review[END_REF] and compare it to the plot of Gaussian density. To illustrate the bivariate CLT, let s = 0.6 and t = 0.2, we use the Python function stats.normaltest from scipy's library to do a normal test based on [START_REF] Pearson | Tests for departure from normality. Empirical results for the distributions of b 2 and √ b 1[END_REF]'s method. This function tests whether a sample differs from a normal distribution. Thus, we plot the p-value of the test conducted on the projection cos(θ)A m (t) + sin(θ)A m (s) for θ ∈ [0, 2π). This method mimics the Cramèr-Wold Theorem 7.1. We also use the Mahalanobis distance to show the joint Gaussianity. This illustration is based on the fact that if X ∼ N d (µ, Σ) then, the Mahalanobis distance D 2 with entries

d -distribution with d degrees of freedom. Using the theoretical value of Σ and numerically estimating the sample average µ, we compute the value D 2 i and compare it with a χ 2 d -distribution, using the QQ-plot method. In Figure 8, we gathered the obtained results from these numerical studies which confirm the joint normal behaviour of (A m (t), A m (s)). We propose two validation methods of the multivariate Gaussianity of the bivariate vector (Am(s), Am(t)) with s = 0.6 and t = 0.2. (b.) p-value for cos(θ)Am(t) + sin(θ)Am(s) for θ ∈ [0, 2π), computed using the function stats.normaltest from scipy. (c.) Estimated probability density of D 2 (blue histogram) and theoretical Chi-squared density (red curve). QQ-plot between the theoretically predicted χ 2 2 quantiles and the empirical quantiles of D 2 .

Appendix B: Comparison between σ 2 (s) and σ 2 k-s (s) variances

Figure 9 below provides a comparison between the two variance values σ 2 (s) in ( 7) (red full line) and (σ k-s ) 2 (s) in ( 14) (black dashed line), for p s ∈ (0, 0.5).

As notice in Remark 4, Figure 9 shows that the variance of the ratio computed using the perimeter is lower than the variance of the ratio computed using a Kolmogorov Smirnov method.

Appendix C: Optimization problem for the variance function

For statistical purposes, it is important to identify the points where the variance function in Equation (4) reaches its maximum and minimum values.

Let us consider the following optimization problem :

By using XCAS, we get that the set of extrema of this variance function is equal to By doing an analytical study, we get the maximum value of the variance is reached on α 1 and α 3 and the minimum value for α 2 (see Figure 10 for an illustration in the Uniform white noise model).


Louise Guillon
email: louise.guillon@cea.fr

Bastien Cagna

Benoit Dufumier

Joël Chavas

Denis Rivière

Jean-François Mangin

Detection of abnormal folding patterns with unsupervised deep generative models

Keywords: Variational autoencoder, Brain architecture, Cortical folding, Anomaly benchmark, Anomaly detection

Although the main structures of cortical folding are present in each human brain, the folding pattern is unique to each individual. Because of this large normal variability, the identification of abnormal patterns associated to developmental disorders is a complex open challenge. In this paper, we tackle this problem as an anomaly detection task and explore the potential of deep generative models using benchmarks made up of synthetic anomalies. To focus learning on the folding geometry, brain MRI are preprocessed first to deal only with a skeleton-based negative cast of the cortex. A variational auto-encoder is trained to get a representation of the regional variability of the folding pattern of the general population. Then several synthetic benchmark datasets of abnormalities are designed. The latent space expressivity is assessed through classification experiments between control's and abnormal's latent codes. Finally, the properties encoded in the latent space are analyzed through perturbation of specific latent dimensions and observation of the resulting modification of the reconstructed images. The results have shown that the latent representation is rich enough to distinguish subtle differences like asymmetries between the right and left hemispheres.

Introduction

The cortex folds in utero to form numerous furrows called sulci, which delimit circumvolutions. Cortical folding is related to cortical architecture (architectony and connectivity) [START_REF] Fischl | Cortical Folding Patterns and Predicting Cytoarchitecture[END_REF]7] and can be impacted by developmental issues that lead to brain disorders [START_REF] Walsh | Genetic Malformations of the Human Cerebral Cortex[END_REF][START_REF] Guerrini | Abnormal development of the human cerebral cortex: genetics, functional consequences and treatment options[END_REF]. The identification of folding patterns acting as markers of developmental brain diseases would be a major breakthrough facilitating early diagnosis. However, although cortical morphology embeds a topography of the sulci sufficiently consistent across subjects to enable the design of automatic recognition tools [START_REF] Borne | Automatic labeling of cortical sulci using patch-or CNN-based segmentation techniques combined with bottom-up geometric constraints[END_REF], the shapes of the sulci present a high diversity, which hinders the modelling of the inter-individual variability necessary to define abnormalities [START_REF] Ono | Atlas of the cerebral sulci[END_REF]. Hence, the diversity of the folding pattern is often put aside and canceled out using spatial normalisation, namely warping all brains toward a template space.

Associations between folding patterns and developmental disorders have already been described. For instance, a very rare pattern called the Power Button Sign (PBS) has been linked to the epileptogenic zone of patients suffering from drug-resistant type 2 focal cortical dysplasia [START_REF] Mellerio | The Power Button Sign: A Newly Described Central Sulcal Pattern on Surface Rendering MR Images of Type 2 Focal Cortical Dysplasia[END_REF]. Similarly, it was demonstrated that the paracingulate sulcus morphology is correlated to hallucinations in patients suffering from schizophrenia [START_REF] Australian | Paracingulate sulcus morphology is associated with hallucinations in the human brain[END_REF]. Sulci shape deviations have also been observed in autism spectrum disorder (ASD) [START_REF] Auzias | Atypical sulcal anatomy in young children with autism spectrum disorder[END_REF][START_REF] Hotier | Social cognition in autism is associated with the neurodevelopment of the posterior superior temporal sulcus[END_REF]. Once abnormal folding patterns linked to a pathology have been identified, automatic detection techniques can be developed using supervised learning. For instance, the PBS can be detected with a supervised classifier [START_REF] Borne | Automatic recognition of specific local cortical folding patterns[END_REF]. However, the upstream process of identifying and defining new patterns of interest is tedious and difficult as each individual has a unique cortical folding geometry and spotting a recurrent abnormality throughout a set of patients is very complex. An unsupervised tool designed to uncover cortical folding abnormalities and potential biomarkers would be an important lever to harvest the potential meaning of unusual folding patterns.

In this paper we propose a dedicated framework based on deep learning and we test its potential through the detection of synthetic unusual folding patterns. The automatic detection of abnormal folding patterns is a challenging task that has not yet been addressed in the field of neuroimaging. In this work, we relate this objective with the general field of anomaly detection. Anomaly and novelty detection aims at identifying samples that do not fit the normal data distribution [START_REF] Pang | Deep Learning for Anomaly Detection: A Review[END_REF]. A few years ago, anomaly detection methods evolved towards deep learning approaches and specifically unsupervised deep learning due to the ability to detect potentially unseen events. Auto-encoder (AE) based methods have been particularly studied as they infer a latent space of interest with much fewer dimensions than the input space, enforcing to learn only the most common features of the training data. There exists a broad range of AEbased models. An extensive review on the detection of epilepsy lesions in MRI patches can be found in [START_REF] Alaverdyan | Unsupervised representation learning for anomaly detection on neuroimaging[END_REF]. Similarly, different methods dedicated to medical images have been compared in [START_REF] Baur | Autoencoders for Unsupervised Anomaly Segmentation in Brain MR Images: A Comparative Study[END_REF] leading to qualify the variational AE (VAE) architecture as the most efficient. Generative adversarial networks (GAN) have also been used in order to identify biomarkers in optical coherence tomography scans, reaching good performances [START_REF] Schlegl | Unsupervised Anomaly Detection with Generative Adversarial Networks to Guide Marker Discovery[END_REF]21]. Based on these initial results, a first framework was proposed for anomaly detection in computed tomography scans of 3D brains where anomalies consisted in labeled traumatic brain injuries [START_REF] Simarro | Unsupervised 3D Brain Anomaly Detection[END_REF]. β-VAE have also been successfully used to model the inter-individual variability in the mouse brain [START_REF] Liu | A Generative Modeling Approach for Interpreting Population-Level Variability in Brain Structure[END_REF]. More recently, very promising self-supervised methods have been applied to anomaly detection problems in medical images [START_REF] Bozorgtabar | SALAD: Self-supervised Aggregation Learning for Anomaly Detection on X-Rays[END_REF], but these methods lack the generative aspect provided by GAN and VAE, which is crucial in terms of explainability. All these works assessed images containing known lesions. Our aim however is to discover still unknown patterns linked to diseases, which leads to challenging evaluation issues. Therefore, this paper is focused on dedicated synthetic benchmark datasets.

In this paper, we propose the first technique aiming at bringing to light unusual and potentially abnormal folding patterns. For this purpose, we propose first a dedicated preprocessing leading to focus the learning on the cortical folding geometry of a specific region of interest (ROI). Then, like in [START_REF] Liu | A Generative Modeling Approach for Interpreting Population-Level Variability in Brain Structure[END_REF], a β -V AE is trained on a set of control data sampling the general population to get a latent representation of the folding pattern distribution in this ROI. We also create several benchmark datasets simulating unusual regional folding patterns to assess the ability of our model to detect them. Finally, we analyse the latent space capacity to separate regional patterns from the two hemispheres.

Methods

Focusing on Folding Information

Brain MRIs contain diverse information that are not all relevant to study folding patterns. Our method therefore includes a crucial first step of data pre-processing based on the BrainVisa/Morphologist pipeline (http://brainvisa.info) [START_REF] Mangin | Object-Based Morphometry of the Cerebral Cortex[END_REF]. This pipeline combines several steps such as bias correction, grey-white segmentation, and skeletonization to obtain a negative cast of the folding. Morphologist's skeletons were used as input of our learning model. They were first defined in [START_REF] Mangin | From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations[END_REF] and are obtained by skeletonization of the grey matter and cerebrospinal fluid union while preserving the topology. The result is 3D volumes with three values: inside of the brain, sulci skeletons and outside of the brain. The use of these simple images rather than raw MRIs puts the focus of learning on the folding geometry and discard a major confound related to the width of the sulci, which increases with local atrophies induced by aging or degenerative pathologies.

Generating Synthetic Brain Anomalies

One of our biggest challenges is the lack of consensual datasets of abnormalities to assess the approach. The examples mentioned in the introduction are either especially challenging in terms of shape (the PBS in epilepsy) or inter-subject local variability (Broca's area or Superior Temporal Sulcus (STS) branches in autism), or correspond to a stratification of the population into several frequent patterns, which is not in the scope of anomaly detection (paracingulate sulcus). Therefore, in this paper, we focus on a 3D ROI of 23x37x36 voxels with 2mm isotropic resolution, localized in only one of these challenging areas, the STS branches. This ROI has been defined in each subject using affine normalisation to the classical MNI reference space. The localization of the ROI in the MNI space has been learned from the open access training dataset with annotated sulci of the Morphologist pipeline [START_REF] Borne | Automatic labeling of cortical sulci using patch-or CNN-based segmentation techniques combined with bottom-up geometric constraints[END_REF]. We have designed several dedicated synthetic anomaly benchmark datasets in order to be able to evaluate the performances of our model (Fig. 1). Deletion: Our first benchmark dataset consists of skeletons in which we have randomly deleted one piece of sulcus, which is chosen among topologically elementary parts called simple surfaces and proposed by the Morphologist pipeline [START_REF] Mangin | From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations[END_REF]. To be deleted, a simple surface must be completely within the ROI and made up of more than 1000 voxels, which corresponds to about 17% of the average number of skeleton voxels in the ROI. This arbitrary threshold aims at performing modification of the geometry beyond the normal anatomical variability observed in the population. Random: A second benchmark dataset is composed of random ROIs of the same dimensions and overlapping the skeleton but localized in different positions in the cortex. This benchmark is expected to be very easily spotted as abnormal since the images are highly different. It is used to ensure that the model is able to identify inputs that are far away from the normal distribution and that what the model has learned is not only non-region specific features such as voxel proportion and sulci continuity. Asymmetry: The last benchmark dataset corresponds to the same ROI but defined in the other hemisphere and flipped. The flip is defined from the interhemispheric plane of the MNI space after affine spatial normalisation. This benchmark has a biological interest as hemispheric asymmetry is still an intense field of research.

Learning a representation of the normal variability

An effective way to model population variability is through AE-based networks. These architectures learn to project input data onto a lower dimensional manifold, also called latent space and to reconstruct from this space the input image. Simple AE are known to have some drawbacks and particularly the lack of regularization of the latent space. To overcome this issue, VAE model was introduced [START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF], and later β -V AE was proposed [START_REF] Higgins | beta-VAE: Learning Basic Visual Concepts with a Constrained Variational Framework[END_REF]. Like classical AE, β -V AE are composed of two parts: an encoder and a decoder but add a variational objective. Contrary to simple AE, an input from image space X , is encoded as a distribution in a latent space Z comprising L dimensions, leading to a twofold objective. First, the minimization of the reconstruction error of the input image. Second, the matching of the encoded distribution to a prior distribution, usually a Gaussian, which is done thanks to Kullback-Leibler divergence and enables to regularize the latent space. VAE is a β -V AE with KL divergence weighted at 1. Thus, β -V AE encoder θ and decoder φ are trained maximising the following objective:

L(θ, φ; x, z, β) = E q φ (z|x) [log p θ (x|z)] -βD KL (q φ (z|x)||p(z)) (1)
where p(z) is the prior distribution, a reduced centered Gaussian distribution in our work that is approximated with the posterior distribution q φ (z|x). Tuning β parameter enables to improve latent factors disentanglement [START_REF] Higgins | beta-VAE: Learning Basic Visual Concepts with a Constrained Variational Framework[END_REF].

Analysing the Latent Space. The analysis of the latent space to understand the meaning of the encoded features is capital to assess the potential of our model to highlight unusual folding patterns. As such, we first trained a β -V AE on normal data only, for our model to learn to encode normal variability. Next, normal and benchmark data unseen during training are projected to the latent space to perform this analysis. The resulting latent codes are used to train a gradient boosting algorithm to classify normal versus synthetic abnormal samples. These three classifiers are used first to ensure that the latent representations are able to capture some relevant information regarding the folding patterns. Then, we can focus our analysis on the features contributing the most to the success of the classification using the generative power of the β -V AE, like in [START_REF] Liu | A Generative Modeling Approach for Interpreting Population-Level Variability in Brain Structure[END_REF]. We travel throughout the latent space modifying only one of these features and observe the generated folding patterns.

Results

Datasets and Implementation

To learn the inter-individual variability of control subjects, the HCP database was used3 . MRIs were obtained with a single Siemens Skyra Connectom scanner with a resolution of 0.7mmx0.7mmx0.7mm. In our work, we studied only the right hemisphere of 997 right-handed subjects with high quality result of the Morphologist pipeline. 547 subjects were used for the training and 150 for the validation of the β -V AE. The remaining 300 subjects were used to train classifiers, half of them being used to create synthetic abnormal patterns following each of the benchmark experiment. Two third of these 300 subjects were used to estimate the β -V AE hyper-parameters using a grid-search driven by a classifier and the last third was used to explore the latent space organization.

The skeletons were spatially normalized with an affine transformation to the standard MNI space and were down-sampled to a voxel size of 2mm, which is sufficient to preserve the folding geometry. Voxels were set to 0 for inside the brain (26% of voxels), 1 for outside the brain (65%) and 2 for the sulci skeleton (9%). As mentioned above, the input to the β -V AE was a 3D ROI, whose bounding box in MNI space was learned from the BrainVISA open access training set of 64 annoted brains, in order to include the two posterior branches of the STS. This 3D ROI is made up by 23x37x36 voxels extended to 40x40x40 using 1-padding.

The complete pipeline is shown in figure 2. Our β -V AE was composed of fully convolutional encoder and decoder with symmetrical architectures comprising three convolutional blocks and 2 fully connected layers. We did a gridsearch (L=8-100, β=1-20, ranges are based on previous works [START_REF] Liu | A Generative Modeling Approach for Interpreting Population-Level Variability in Brain Structure[END_REF] and reconstruction abilities), where hyperparameters were chosen according to the classification performances of the Deletion classifier applied to 100 controls and 100 synthetic samples using a 5-fold stratified cross-validation. We selected L=100, β=2 and a learning rate of 2e-4. Training lasted for 400 epochs on a Nvidia Quadro RTX 8000 GPU and was completed in roughly 1 hour. Fig. 2. Whole pipeline. First, bounding boxes of sulci of interest are defined. The HCP database is processed with the Morphologist pipeline and cropped thanks to the bounding boxes. Crops are then downsampled before feeding the β -V AE.

Analysing Learned Folding Variability

First, we visually evaluated the reconstruction ability of our β-VAE. 2D slices of several reconstructed inputs are presented in figure 3. For the control set, deletion and asymmetry benchmarks, the reconstructions are approximate but retain the general geometry. The main sulci included in the ROI can be identified. In return, the random crops cannot be reconstructed by the decoder, which outputs an image looking like a disturbed configuration of STS branches. This suggests that the model has really learned the distribution of the specific geometry of this cortical region rather than a generic distribution covering any skeleton configuration. To evaluate our model latent space, three gradient boosting classifiers were trained using 50 controls and 50 synthetic samples from the test set, using a 5-fold stratified cross-validation. ROC curves are presented in figure 4.A). For each configuration, the AUC score is above chance. Highest performances are obtained on benchmark Random (AUC = 0.98), which was expected. This first result comforts the fact that our model is able to capture very obvious abnormalities. On benchmark Asymmetry, very good scores are also obtained (AUC = 0.85). Though for an inexperienced eye the difference between right hemisphere and flipped left hemisphere can seem subtle, this region is known to be asymmetrical in terms of length and tilt of the main sulci [START_REF] Rubens | Asymmetry of the lateral (sylvian) fissures in man[END_REF]. It confirms that our model is capable of representing specific anatomical structures included in the ROI. Finally, when deleting one large simple surface, results are slightly above chance (AUC=0.69) which indicates a potential of the latent space to detect such anomalies. Nevertheless further work is required using a larger test set to overcome potential limitations of this classifier experiment. Using gradient boosting classifiers gives insights on the most decisive dimensions of the latent space, which depend on the benchmark. Figure 4.B) shows a visualization of the datasets using the two most important latent features or tSNE 2D manifold projection. For Random and Asymmetry benchmarks, two groups can be clearly identified. Surprisingly, in the tSNE visualization, the random crops are surrounded by controls, which is counter-intuitive relative to the Gaussian prior and does not fit with the plot performed from the two most discrimative dimensions. Further experiments with unbalanced datasets including only a small ratio of abnormality will help us clarify this point. The Suppression benchmark is clearly the most difficult classification experiment with the current latent space.

The last experiment consists in sampling latent vectors to travel along the most important dimensions according to the classifiers. the Sylvian fissure (SF), seems to shorten from the left to the right of the dimension. This observation is consistent with Asymmetry benchmark's subjects distribution according to the two most important features on figure 4.B). Indeed, control subjects, i.e. right hemispheres correspond to higher values of the 88th dimension which is also suggested by the SF shortening on figure 4.C). This evolution is interesting as previous works demonstrated that the SF was shorter in the right hemisphere. Thus, this dimension could encode the length of the SF which is an important marker of laterality [START_REF] Idowu | Morphometry, Asymmetry and Variations of the Sylvian Fissure and Sulci Bordering and Within the Pars Triangularis and Pars Operculum: An Autopsy Study[END_REF].

Discussion and Conclusion

In this paper we developed a framework that shall lead to the discovery of abnormal folding patterns beyond reach for human cognition because of the high inter-subject normal variability. Our main contribution is the design of synthetic benchmarks used to decipher the organization of the latent space used to model this variability. We have shown that the regional specificity of the folding pattern can be learned and used to detect some deviations from the norm. Our methods achieved to detect obvious and more subtle deviations (respectively with Random and Asymmetry benchmarks), but detection is harder for more complex ones such as Deletion benchmark. However, we don't seek to detect benchmark subjects but rather to detect abnormal patterns. In the future we plan to perform further experiments to get more hints about the nature of the representation. Although some interesting features can be observed when varying latent vectors dimensions' values, it is very difficult to apprehend all possible shifts looking at 2D sections because we deal with 3D images. Additionally, in some cases, we observe that the generated images are not realistic, preventing the interpretation. Generating synthetic anomaly benchmarks is a very useful initial step but induces biases that must be acknowledged. In our case, with the Deletion benchmark we made the hypothesis that anomalies could be linked to "missing" simple surfaces. This work constitutes a first step of proof of concept and enables to control the complexity of abnormalities, however an important limitation is the high dependency of the results on simulated data as our benchmarks are made of synthetic anomalies. Nevertheless we stress out that the model has learned only control data, thus is totally unsupervised in regards to anomaly data, and only the evaluation results depend on simulated data. In future works, the ultimate benchmarks will have to be built from large datasets of neurodevelopmental disorders and will aim at discovering actual abnormal folding patterns.

Fig. 1 .

 1 Fig. 1. Generated anomaly benchmarks. Benchmark Deletion: Original crop and its modified version. Benchmark Random: two examples of random crops. Benchmark Asymmetry: crop of right hemisphere and crop of left hemisphere flipped.

Fig. 3 .

 3 Fig. 3. Reconstructions of test inputs. 2D sections from the 3D ROI presented in sagittal view at depth 18. First row: original images, second row: model outputs.

Figure 4 .

 4 C) shows the reconstruction provided by the decoder when following the dimension corresponding to the best discriminator of Asymmetry benchmark. All other dimensions were fixed at their mean value. As expected, the outputs look like approximated skeletons of the ROI. Subtle changes consistent during the travel can be observed: i.e. the upper part of the sulcus presented in pink on figure 4.C), called

Fig. 4 .

 4 Fig. 4. Analyses of latent space. A) ROC curves of GradientBoosting classifiers for classifications between controls and benchmarks. B) Distribution of data points using classifiers most important features and t-SNE. C) Travelling through the latent space.

https://www.humanconnectome.org/

Acknowledgments

This project has received funding from the FRMDIC20161236445, the ANR-19-CE45-0022-01 IFOPASUBA, the ANR-14-CE30-0014-02 APEX the ANR-20-CHIA-0027-01 FOLDDICO. Data were provided in part by the Human Connectome Project funded by the NIH. This work was performed using HPC resources from GENCI-IDRIS (Grant 2020-AD011011929).