Irène Kaltenmark

Christine Deruelle

Lucile Brun

Julien Lefèvre

Olivier Coulon

Guillaume Auzias

Cortical inter-subject correspondences with optimal group-wise parcellation and sulcal pits labeling

Keywords:

analyze and to visualize the inter-subject variability. We implemented these methodological advances and we applied them on two different populations of healthy subjects. The first database of 137 adults allows us to compare our method to the previous one and the second database of 209 children illustrates the robustness of our method on heterogeneous data.

Introduction

The concept of sulcal roots emerged 20 years ago in [Rgis 95, Rgis 05] as indivisible atomic folding entities located in the deeper parts of sulci and stable across subjects. In [Lohm 00, Yang 08, Li 09], the authors proposed a framework to decompose the folds of the cortical surface into several parcels called sulcal basins using a watershed approach. Sulcal pits were then formalized in [Lohm 08] as the deepest points of each sulcal basin, to give a concrete representation of the abstract concept of sulcal roots.

A major interest of sulcal pits is to define a finite set of robust landmarks on the cortical surface allowing local comparison across individuals. They have often been used to describe adult cortical morphometry [Im 10, Auzi 15, Take 17, Le G 18] or atypical sulcal pattern in pathologies [Im 16]. The sulcal pits are also of interest for tracking the growth and the evolution of the geometry of the cortical surface during brain maturation [Lefv 09]. The longitudinal study of [Meng 14,Meng 18] on the first two years of life highlighted a stable spatial distribution of pits in the major deepest folds already present at term birth. Authors also reported a regionally heterogeneous increase of sulcal pits depth, especially in the high-order association cortex. In [Auzi 15], authors showed that the notion of sulcal pit should not be reduced to the deepest areas of the cortical surface. The first study on young children covering both deep and shallow folds showed that the number of shallow pits increases with age between 18 months and 10 years [Brun 16]. However, the precise location and time table of these changes remain unexplored. A thorough review on sulcal pits studies has recently been presented in [Im 18].

In this work, we highlight the methodological challenges that we are facing when working with sulcal pits covering the whole brain. The identification of homologous pits across subjects is made difficult by the large inter-subject variability. We therefore introduce a combination of new techniques to address this issue, paving the way for the fine-grained characterization of the spatio-temporal patterns of sulcal pits, including in shallower folds.

Context

Sulcal pits extraction with a watershed by flooding algorithm

A decade ago, [Im 10] proposed a surface-based approach to extract sulcal pits from an individual cortical surface. A preliminary step consists in the estimation of a sulcal depth map. A watershed by flooding algorithm [Rett 02, Yang 08] then orders the vertices of the cortical surface by their depth. The deepest vertex defines the first sulcal pit, the initial vertex of a sulcal basin. If the next vertex in the list is the neighbor of the previously identified sulcal basin, it is added to this sulcal basin. If all of its neighbors are unlabeled, it defines a new sulcal pit as a seed vertex for a new sulcal basin. This algorithm simulates water that would rise and fill the folds of the cortical surface with respect to the level sets of the depth map. By construction, the sulcal pits are the local maxima of this map. This watershed is applied to a smoothed depth map that still contains spurious extrema related to anatomically irrelevant variations in the folds of the cortical surface. During the flooding, neighboring basins are thus merged when one of the basin contains a noisy pit. The shallowest basin is merged with the deepest one according to three geometric criteria: the ridge height (R), the basin area (A), and the geodesic distance (D) between the two pits (see [Im 10] for details). [Auzi 15] adapted this method by delaying the filtering with respect to the basins' areas at the end of the flooding procedure.

Regarding the estimation of the sulcal depth, the approach used in [Im 10] (see also [Meng 14]) consists in computing the Euclidian distance to the closest point on the external hull of the brain. an estimation of the depth of the folds based on the concept of bending distance for curves. This method does not require the extraction of the brain hull and the resulting map is independent of the brain size which is convenient for parameter settings [Auzi 15].

Inter-subjects labeling of the sulcal pits

Once sulcal pits have been extracted from each individual cortical surface, we then need to identify homologous pits across subjects. The aim here is to associate a label that is common to all subjects to each sulcal pit, which defines a nomenclature of sulcal pits and a corresponding atlas of sulcal basins. In the ideal case of a perfect inter-individual match, each of the resulting atlas label would correspond to one and only one sulcal pit for each subject, hereby defining a consistent correspondence of pits across individuals. Because the pit extraction process is subject to noise and because of inter-subject cortical variability, we expect to find subjects for whom a given sulcal pit is missing, or for whom several pits end up being associated with a single atlas basin.

The method to label sulcal pits across subjects initially proposed in [Im 10] consists in a clustering algorithm by another watershed flooding on the density map of sulcal pits after registration of the subjects' cortical surfaces. The labeling is then straightforward: each pit inherits the label of the cluster that contains it. In that work, the authors manually selected the clusters located in the deep major sulci and excluded others. [Auzi 15] yet demonstrated that while deeper pits are less variable across subjects than superficial ones, some shallower sulcal pits also define robust landmarks of the cortical surface. The authors proposed to extend the method from [Im 10] in order to integrate the shallow pits in the labeling system and avoid the manual selection of deeper folds. Importantly, the one-to-one correspondence of sulcal pits across subjects was used as a measure of the sulcal pits extraction quality. This measure was then instrumental to determine optimal settings of the sulcal pits extraction algorithm. Between these two approaches, the number of sulcal pits by subject roughly increased from 50 to 100 per hemisphere. Accordingly, the number of group-level pit labels in the respective atlases increased from 48 in the left and 47 in the right hemisphere in [Im 10], to 104 and 114 respectively in [Auzi 15]. These numbers illustrate the tremendous increase in the complexity of inter-subjects labeling task when including shallower pits corresponding to more variable folds.

If a subject has more than one pit present in a cluster, there is no consensus on the automatic selection of the unique pit identified to the cluster. In [Im 10], authors use the distribution of sulcal pits in the cluster to select the pit that is closest to the densest point. In [Auzi 15, Le G 18], the deepest pit is preserved for analysis. In [Le G 18], authors investigated the distribution of the sulcal pits depth in each cluster prior to the selection of the deepest pit of each subject. In some deep clusters, they observed two underlying Gaussian distributions. A prominent Gaussian corresponds to the deep sulcal pit of the area and a minor one, which presents a higher variance and a shallower mean depth, corresponds either to a bump in the sulcus shape that leads to the identification of an additional shallower pit (see also some minor patterns in [Meng 18]) or to a noise artifact of the depth potential map. After selection of the deepest sulcal pits, the minor Gaussian is only partially reduced which illustrates how shallow sulcal pits can be wrongly identified with a deep sulcal pit when this last one is missing in a subject. Authors proposed therefore to also add a threshold on the depth to remove these outliers. However, after thorough observation of individual data, we noticed that these outliers could also result from incorrectly registered sulcal pits of neighboring clusters. Due to the spatial inter-individual variability, sulcal pits can be located at the boundary of some clusters, in which case their position are often unreliable. In our application detailed later, we count 8.2% of pits that are projected on the boundary of a cluster and 15.3% are very close to the boundary (in a two-ring neighborhood), questioning the robustness of their affectation since the labeling method relies strictly on the pit's position. Yet, the most natural way to manually label a sulcal pit is to identify the part of the sulcus where the pit is located. This means to observe the local and geometric environment of the sulcal pit. To automate this process, the idea developed in this paper is to exploit the sulcal basins generated by the individual watershed algorithm (used to extract sulcal pits on baseline cortical surfaces). This algorithm produces indeed a complete parcellation of each individual cortical surfaces. Likewise, the group-wise watershed algorithm on the density map of sulcal pits produces an atlas parcellation of the atlas surface. The clusters of this parcellation resembles the individual sulcal basins after registration (see Fig. 9 in [Auzi 15] or Fig. 6 in [Im 10]). We therefore propose to automate the geometric comparisons between individual sulcal basins and clusters of sulcal pits to improve the group-wise identification. Regarding more specifically the group-wise watershed, the optimization of its parameters on a heterogeneous pediatric database has proven to be difficult. The algorithm is subject to the same type of parameters A, R, and D than the individual watershed algorithm. No matter the choice of the parameters, we always obtain a small number of poorly robust clusters (this will be illustrated at the end of this paper). Moreover, the new geometric viewpoint that we proposed above raises a conceptual problem with the method itself. Indeed, the clustering is grounded on the density map of sulcal pits. This map only integrates the spatial distribution of the pits on a common space and does not exploit their associated sulcal basin. This limits the geometric similarities between the clusters and the individual sulcal basins. A thorough observation of the resulting atlases illustrates the inherent limitations of the current approach. The left panel in Fig. 2 is a zoom on the central sulcus (CS) of the left hemisphere. The median blue cluster (1) of the CS is significantly smaller than usual individual sulcal basins at this place (see e.g. Fig. 3 in [Meng 18]). The large red basin (2) on top of the gyrus actually represents a small shallow but robust sulcal basin whose shape is more similar to the white cluster (3) on the opposite gyrus. Moreover, some clusters present a tail as pointed by the arrows in Fig. 2. This type of shape is however never observed on actual individual sulcal basins. Consequently, when a sulcal pit is located in such a tail, it usually inherits a wrong label. These observations led us to propose a new method to build an atlas that takes into account the shape of individual sulcal basins.

Contributions

Introductory remark: we study in this paper the sulcal organization of individual cortical surfaces after surface registration (e.g. with FreeSurfer) on an average mesh. This organization consists of a parcellation in sulcal basins with a unique sulcal pit in each basin. Any mention of an individual feature, e.g. sulcal pit or sulcal basin, henceforth refers to the registered feature on the atlas mesh.

The core of this paper is the introduction of a new pipeline to generate a group-wise atlas of the sulcal basins and a corresponding labeling of individual sulcal. The first process, presented in Section 2, builds a group-wise parcellation of the cortical surface. Each region, that we call atlas

Geometric atlas

Algorithm 3: Geometrical pairing between individual sulcal basins and atlas basins

Varifold labeling

1) Atlas map of sulcal basins 2) Individual labeling of sulcal pits and sulcal basins with respect to the atlas

Final output

Old method

Clustering on sulcal pits with a watershed by flooding algorithm on the density map of sulcal pits Straightforward labeling: each pit inherits the label of its parent cluster basin, integrates the geometry of homologous individual sulcal basins. This new type of atlas is then the key element of the second process, presented in Section 3. It consists in an automatic labeling of sulcal pits by comparison between atlas and individual sulcal basins, with a similarity measure on surfaces. In the experiments, we use a measure induced by varifold metrics that are highly efficient as they are invariant to the parametrization of surfaces. A novelty is also to integrate in the labeling system the concept of unlabeled sulcal pit, called isolated pits, in order to guarantee for each subject the uniqueness of the sulcal pit identified by a label.

The pipeline is applied to two different populations of healthy subjects. The first database consists of 137 adults from [Auzi 15] to show the improvement compared to previous approach and the application to a second population of 209 children demonstrates the efficiency of the approach in a different context. Results are displayed and discussed in Section 4.

Note at last that according to our best knowledge, the variability of these sulcal basins has never been investigated. We highlight their robustness through density maps of homologous sulcal basins. These maps also define a new qualitative measure to evaluate the generated atlases.

Geometric atlas 2.1 Heuristic

Under the assumption that homologous sulcal pits across subjects are contained in similar sulcal basins, the aim of the method presented hereafter is to take into account the shape of individual sulcal basins to generate a group-wise parcellation of atlas basins on an average mesh, that will serve to define the sulcal pits labeling. More precisely, atlas basins will locally play the role of geometric representatives of individual sulcal basins (in terms of location and area).

The main ideas of the method are as follows (each step will be detailed in the following sections). Given a smoothed spatial distribution of all the sulcal pits of a population on an average mesh as found in [Auzi 15] or [Im 10], we extract vertices with high concentration of pits. These vertices are called seeds. A cluster is initialized around each seed as a small neighborhood of vertices of the atlas mesh. Clusters will then grow vertex by vertex until they form a complete parcellation of the mesh. By construction, as soon as a cluster is created, it already contains sulcal pits from a set of individuals. For each one of these pits, there is an associated individual sulcal basin whose shape should approach the final shape of the cluster. These sulcal basins will therefore drive the growth of the cluster. To integrate this geometric information, the cluster will start to grow inside the intersection of all its associated basins. Once this area is filled, we progressively label the neighboring vertices that belong to most of the basins associated to the cluster. For this purpose, we consider the density map of these basins. By construction, the cluster is located around the maxima of this local map. It will then grow by filling the level sets of this map: the vertex, on the external boundary of the cluster, that maximizes this density map is the next vertex that should be added to the cluster. At last, new sulcal basins are gradually associated to the cluster: each time a new vertex is assigned to a cluster, the pits located on this vertex are associated to the cluster, with their respective basin. As the number of associated basins increases, the density map becomes more robust. Therefore, this map is updated as soon as a new basin is associated to the cluster.

This new clustering method consists therefore in a watershed by flooding algorithm where the usual depth map is replaced by local density maps of sulcal basins. As illlustrated in Fig. 4, these local maps can be read as local depth maps. This approach has two particularities:

1. These density maps rely on the identification of the homologous sulcal basins across the population. They will be estimated with few subjects at the first step of the algorithm and progressively updated as the clusters will grow.

2. Some of these maps have a strictly positive value on the same set of vertices which induces a boundary uncertainty (see Fig. 4). Section 2.2 describes the initialization of the clusters. Section 2.3 introduces the local density maps of sulcal basins. Section 2.4 presents the main algorithm that provides a complete parcellation of clusters and Section 2.5 details the final step to remove noisy clusters. The remaining clusters are then finally called atlas basins as they represent average sulcal basins.

Mesh

Notation

Throughout this paper, we will denote B A (k) the k-th cluster of the atlas map. As these clusters will progressively evolve, they will finally be called atlas basins at their final state. Likewise, we will denote B i (k) the k-th sulcal basin of the i-th subject of a given population. At last, x ∈ X will denote a vertex of the average cortical mesh.

Initialization: seeds of the clusters

The projection of smoothed individual sulcal pits on the atlas mesh allows to sum across subjects hereby resulting in a density map of sulcal pits showing the consistency of localization of the pits across individuals. The sum is normalized by the number of subjects (see Fig. 5). Vertices with locally maximal density are extracted and called seeds of the clusters. For each seed s k , we initialize a new cluster B A (k) (to read k-th Atlas Basin) with the vertex of the seed and its two-ring neighborhood (see Fig. 5). On the individual cortical surface, sulcal pits are located in distinct sulcal basins and are by construction well separated. Two clusters that are not disjoint at this stage contain thus homologous sulcal pits across the population. In this case, we delete the cluster with the lowest peak of pits density. This ensures that all clusters are disjoint. Our algorithm will then add each remaining vertex to one of these initial clusters starting by those located at the boundary of these clusters (Algorithm 1 Step 1). Before the presentation of this algorithm, let us introduce the notion of influence.

Influence maps

For any cluster B A (k), we denote Subj k the list of individual sulcal basins associated to the cluster joint with their respective parent subject. More precisely, (i, j) ∈ Subj k says that the j-th sulcal basin of the i-th subject is associated to B A (k). This list will allow us to control that a subject cannot have two distinct sulcal basins associated to the same cluster. We can now introduce the notion of influence map which is the local density map of sulcal basins associated to a cluster. The influence of a cluster B A (k) on a vertex x is the percentage of individual sulcal basins associated to B A (k) that contains x. For any vertex x,

I k (x) = 100 #Subj k (i,j)∈Subj k 1 1 Bi(j) (x) , (1)
where B i (j) denotes the j-th sulcal basin of the i-th subject, #Subj k is the number of pair (i, j) in Subj k , i.e. the number of associated basins, and 1 1 Bi(j) is the indicator function of the basin B i (j). For example, I k (x) = 100% means that x belongs to all the sulcal basins associated to the cluster B A (k). Fig. 6 illustrates the construction of an influence map and how this map integrates both the average location and shape of homologous sulcal basins. In a perfect situation with no variability across subject, the influence map of a cluster would be an indicator function of all the identical sulcal basins: it would be equal to 100% in the intersection of the basins, and 0% outside. As we want to define clusters that retrieve the shape of their associated basins, the initial cluster implicitly represented in the left panel of Fig. 6 should extend in order to roughly match the red area of the right panel. Clusters will thus grow where their influence is high. For more illustrations, the right panel in Fig. 9 displays the influence maps of six initial clusters delimited by the blacks areas. Since we have as many influence maps I k : X → [0, 100] as clusters, we need to adapt the standard watershed by flooding algorithm. A vertex x can indeed be under the influence of multiple clusters, i.e. we can have I k (x) > 0 and I k (x) > 0 with k = k (see the boundary uncertainty in Fig. 4). To find the most pertinent map for each unlabeled vertex x ∈ X, we identify the cluster B A (k x ) with the maximal influence on x:

k x = arg max k∈[1:n] I k (x) , (2)
where k ∈ [1 : n] indexes the n clusters. This is to say that x should be added to the cluster B A (k x ).

Adaptive watershed algorithm

Algorithm 1 presents the core of our adaptive watershed algorithm. We will discuss hereafter the details of the different steps. The output of this algorithm is a parcellation of the average mesh. As we will see, a particularity of this algorithm is to exploit the one-to-one correspondences assumption, meaning that two sulcal pits of the same subject cannot be associated to the same cluster. For the sake of clarity, we will present in Annex a complete version of the algorithm, Algorithm 4, which additionally ensures that clusters remain connected.

Algorithm 1: Adaptive watershed algorithm (simplified version) Input : Subjects' sulcal pits and sulcal basins registered on an atlas mesh.

Output: The atlas basins, which form a parcellation of the atlas mesh. A labeling system with respect to this parcellation (given by the sets (Subj k ) k )).

Initialization:

• Compute the density map of pits.

• Extract the seeds, which are the maxima of the density map.

• Initialize the clusters by the 2-ring neighborhood of the seeds.

• Initialize the respective lists Subj k .

Watershed:

1. Compute the list Q ⊂ X of unlabeled vertices at the boundary of the clusters. If Q is empty, the atlas is complete. Otherwise:

2. Using the lists (Subj k ) k , compute the influence map of each cluster and extract the vertices under maximal influence:

arg max x max k∈[1:n] I k (x)
3. Select among them a vertex x with minimal conflict C kx (x).

4.

x inherits the label k x where k x denotes the cluster with maximal influence on x.

B A (k x ) ← B A (k x ) ∪ {x}.
5. Update Subj kx : pits on the vertex x are associated to the cluster B A (k x ).

6. Return to step (1).

Learning process: updating of the influence maps

Let us recall that the influence map of a cluster B A (k) is the normalized sum of the individual sulcal basins associated to this cluster by the set Subj k (see Fig. 6). The influence maps rely thus on the sets Subj k . The clusters' initialization to a 2-ring neighborhood of the seeds is intended to have enough sulcal pits to generate a robust first estimation of the influence maps of each cluster. However, during its growth, a cluster will progressively contains more sulcal pits. An increasing number of individual pairs of sulcal pit and sulcal basin can thus be associated to this cluster by the set Subj k . The input of these new subjects should strengthen the influence map of the cluster.

The simplest approach to define Subj k is to list of all sulcal pits whose projected position is inside B A (k). However, in order to improve the robustness of the influence maps, we exploit the one-to-one assumption across subject and we define a geometrical criterion to assign a new pit to a cluster. For each new vertex x added to the cluster B A (k x ) (Step 4 of the algorithm), we update the list Subj kx as follows (Step 5). We retrieve all subjects with a sulcal pit on x as a list of pairs (i, j) such that the i-th subject has the pit from its j-th basin at the new vertex x. Then each pair (i, j) is added to Subj kx if and only if it satisfies both following conditions:

1. subject i was not already indexed by Subj kx , 2. the sulcal basin j of subject i contains the seed of the cluster B A (k x ).

Figure 7: Two scenarios to illustrate the influence map updating criteria in two ambiguous situations. We display for few subjects associated to a growing cluster (not displayed) their respective basins (black dashed lines) and their pits (black dots). The seed of the cluster is marked by a +. We then illustrate the effect of the two criteria with a new subject in red. Left scenario: assume that the red pit on the left has already been associated to the cluster. When the cluster grows on the right and meet the second red pit, this last one is not assigned to the cluster as it does not satisfied (1). This prevents to add a second pit from that subject and its basin to the influence map of the cluster which would drive the cluster to keep extending on the right. Right scenario: the cluster does not yet contain the red pit on the left. If it first meets the red pit on the right, this last one will not be assigned to the cluster since it does not satisfied the geometrical criterion (2). If the cluster then grows on the left, it will still have the possibility to retrieve the left red pit and its basin.

The first condition ensures that two sulcal pits of the same subject cannot be associated to the same cluster. The second condition formalizes the influence of individual basins on the corresponding atlas basin only if their geometry are similar. See Fig. 7 for schematic illustration. Fig. 8 illustrates the growth of a cluster and the update of its associated sulcal pits and sulcal basins on real data.

Through this continuous update of the influence maps, the algorithm is progressively learning about the population (see the evolution of the influence map in Fig. 8). In order to take advantage of this learning process, we define a priority order between the unlabeled vertices by selecting those under maximal influence (Algorithm 1 Step 2) :

x = arg max x I kx (x) , (3)
where k x indexes the cluster with maximal influence on x (see eq. ( 2)). This selection automatically handles the simultaneous growth of all clusters instead of flooding one by one each cluster according to their respective influence map. These black clusters will then preferentially grow where they have a high influence.

Conflict of influence

To further exploit the learning process of the algorithm, we introduce the notion of conflict that allows to delay the most difficult labeling choices. Since eq. (3) usually have multiple solutions, we will refine the selection. Due to the variability of sulcal basins across individuals, the influence maps of neighboring clusters may overlap. This phenomenon highlights the uncertainty of the boundary of these clusters (see Fig. 4 and9). The conflict of influence will quantify this uncertainty. Consider x ∈ X an unlabeled vertex and B A (k x ) the cluster with the maximal influence on x. The conflict with B A (k x ) results from the influence of other basins and is defined by

C kx : X → R+, C kx (x) = k =kx I k (x) 2 . ( 4
)
This conflict measures the uncertainty of the labeling by k x at any vertex (see Fig. 10). Therefore, between two vertices under maximal influence (i.e. solutions of eq. ( 3)), we will first label the vertex with the lowest conflict (Algorithm 1 Step 3).

Figure 10: Left: influence map of an initial cluster. Right: cluster in red and conflict map of this cluster in blue. Note that the conflict is only relevant in the neighborhood of the cluster. Our approach consists in penalizing the cluster's growth in the blue area.

Example 2.1. Assume that the cluster B A (k) has the maximal influence on the unlabeled vertices and denote m this maximum. We have to choose a vertex at the boundary of B A (k) in the level line I -1 k (m). Assume for example that this level line intersect the boundary of the cluster on its right side and on its left side. If there is a close other cluster on the right side, the conflict should be high in this area. If the vertex that minimizes the conflict is on the left side of the cluster, the cluster will grow on this left side first. The notion of conflict thus drives the cluster to extend in the less uncertain area (see Fig. 11).

Remark 2.2. Equation (4) considers the squares of conflicting influence to penalize vertices with fewer conflicting clusters. Consider a vertex x with one conflicting cluster for which 40% of subjects contain x and a vertex y with two conflicting clusters for which 20% of subjects contain y, then C(y) = 20 2 + 20 2 < 40 2 = C(x). The labeling of x is thus considered more uncertain than the labeling of y. Hence, if x and y are under the same influence, y will be labeled before x.

Isolated sulcal pits

The robustness criteria presented in Section 2.4.1 can prevent a sulcal pit to be associated to the cluster that contains the pit's vertex. When a pair (i, j) does not satisfy one of these two criteria, the vertex is added to the cluster but we investigate if the pit and its basin B i (j) can be assigned to a neighboring cluster. If another cluster B A (k ) covers more than 50% of the sulcal basin B i (j) and if the two robustness criteria are satisfied with respect to Subj k and s k , then (i, j) is assigned to this last cluster. In Fig. 7, this means that the right red basin in both scenarios could be assigned to a neighboring cluster.

When a pit is rejected by one of the two conditions for all the potential clusters, the pit remains unlabeled. We call these pits isolated pits. They can be individual noisy pits or part of a minor local pattern in the population, e.g. shallow pits on some gyri.

Basins' filtering

The previous algorithm never merges two clusters. Yet, in practice, the algorithm provides about 170 clusters when subjects have about 100 sulcal pits, per hemisphere (see Fig. 5). This phenomenon results from the large inter-subject variability and the inability of the surface registration to remove this variability, hence leading to an irregular density map of pits that usually induces too many seeds (local maxima of density). In the next section, we detail how to delete a cluster and we then present second algorithm to select and delete expendable clusters.

Figure 11: Three growing clusters and their respective influence maps (their values are displayed by the gradient of intensity in the faded areas). The seed s k of the cluster B A (k) is marked by a +. The brown pixel, on the left, represents an unlabeled vertex under maximal influence: it will be the next pixel added to the red cluster. The black pixel in the middle, representing another unlabeled vertex, is under the influence of both red and green basins. If the red basin has the maximal influence, the influence of the green basin is then seen as the conflict on this vertex. Even if the influences on the brown pixel and the black pixel are equal, the brown pixel has no conflict and will therefore be labeled in priority.

Cluster deletion

With a standard watershed algorithm, the basins' filtering is achieved by a straightforward merging between neighboring basins. This solution is yet not always the most pertinent one. In some situations, we would rather split a basin between several neighbors. With our adaptive watershed algorithm, the influence maps play the role of the standard depth map. In order to delete a cluster B A (k D ), we will reset all the clusters to the initial 2-ring neighborhood of their respective seed after excluding the seed s k D so that the influence map I k D does no longer exist and then apply the Watershed function of Algorithm 1. Since I k D was playing the role of the local depth in the area of the cluster B A (k D ), this means that we have flattened the depth in this area so that vertices are now attracted by the neighboring clusters. See an example of cluster deletion in Fig. 12. As we will see hereafter, this process will be applied several times to quantify the effects of a cluster deletion (factorial complexity with respect to the initial number of clusters). As it is numerically expensive, we propose a faster version in which we only reset the labels of the vertices in the cluster to remove.

Cluster selection

To produce the final atlas basins, several clusters must be deleted iteratively. The method proposed hereafter relies on two key aspects: an automatic selection of a cluster to delete and a stop criterion. According to the hypothesis that there exist one-to-one inter-subject correspondences for most of sulcal pits, most clusters should have a number of associated sulcal pits close to the number of subjects. A basic approach would thus be to delete the clusters with the smallest number of sulcal pits, for example, with less than 25% of the population size. However, consider two clusters that devise what should be a single atlas basin. In the most critical situation, the two clusters may each label 50% of the population and we could thus not merge them with the previous heuristic.

When a cluster is deleted, the total number of individual sulcal pits associated to a cluster • Reset all the vertices and reinitialize all the clusters, but B A (k D ), with the 2-ring neighborhood of their seed.

• Delete the seed s k D .

• Run the Watershed function of Algorithm 1.

}

Main part of the algorithm

Compute N 1 .

Use Deletion complete version to delete one by one the clusters with a N 1 score lower than 10% and compute N 1 after each cluster deletion.

While min k∈[1:n] N 1 (k) < p, do: 1. For k 0 ∈ [1 : n], do: i) Simulate the deletion of the cluster B A (k 0 ) with Deletion fast version.
ii) Compute the temporary N 1 score, denoted N k0 1 , associated to the resulting parcellation without the cluster

B A (k 0 ). iii) Compute Tot(k 0 ) = k∈[1:n],k =k0 N k0 1 (k). 2. Select k D = arg max k∈[1:n]
Tot(k).

3. Delete the cluster B A (k D ) with Deletion complete version.

Re-index the remaining clusters

, n = n -1, N 1 = N k D 1 .
of the new parcellation tends to decrease. This means that some sulcal pits become isolated pits (see Section 2.4.3). The choice of the cluster to delete should aim at minimizing this number of lost pits, i.e. the number of additional isolated pits after a cluster deletion. For this purpose, let us recall the N 1 score, introduced in [Auzi 15]. This measure is defined on each cluster B A (k) as the percentage of subjects that have one and only one pit in the cluster:

N 1 (k) = #Subj k n S
where n S is the total number of subjects. The number of additional isolated pits induced by a cluster deletion is the difference between the sum of the N 1 scores before and after the cluster deletion. For example, consider again the example of two clusters that divide what should be a single atlas basin. Assume that both clusters have a N 1 score of 50%. In a perfect situation, the deletion of one of the two clusters should generate 0 additional isolated pits. All the pits associated to the deleted cluster would then be associated to the remaining cluster. The N 1 score of this last one would increase from 50% to 100%.

In Algorithm 2, we will estimate this number of additional isolated pits for every single cluster deletion possible and select a posteriori the optimal cluster to delete. In practice, minimizing the number of additional isolated pits is equivalent to maximizing the sum of the a posteriori N 1 scores across every possible single cluster deletion. In Algorithm 2, the deletion simulations are achieved in Step (1) and the optimal configuration is identified in Step (2). For these simulations, we use the fast cluster deletion. Once we have estimated the optimal configuration, i.e. the optimal cluster to delete, we compute in Step (3) the new atlas basins with the complete version of our deletion process.

Finally, we define the following stop criterion. The cluster deletion stops when the minimal N 1 score of the current clusters is above a given threshold p.

Remark 2.3 (Algorithm details).

To fasten the algorithm, we first delete the clusters with a N 1 score lower than 10%. We also only consider in Step (1) the clusters with a N 1 score lower than 70% since we do not want to delete well filled clusters. In practice and in order to produce a more robust stop criterion, we consider the average N 1 score of the 5 less reproducible clusters. The threshold for this average is set to p = 25%. The stop criterion of the while loop is:

mean(N 1 (k i )) i=1..5 < p , where for any k / ∈ (k i ) i , N 1 (k) > max i N 1 (k i ) .
Remark 2.4 (Runtimes). For one hemisphere mesh of 41000 vertices and a population of 137 subjects, Algorithm 1 takes about 30mn to run. The filtering Algorithm 2 is the most time-consuming as we simulate the deletion of each cluster at each iteration (factorial process with respect to the number of initial clusters). The trick to use a fast version of the cluster deletion function reduces the run time from 4-5 days to 2-3 hours.

Parameters of the algorithm

The robustness criteria introduced, in Section 2.4.1, in order to validate the identification of individual sulcal pits and sulcal basins, are implicit parameters of this algorithm. Although experiments showed us the need to ground the flooding algorithm on most robust basins, the choice of the criteria can be investigated in future work, especially the geometric criterion.

The main parameter of the algorithm is the threshold to stop the cluster filtering: the minimal N 1 score of the parcellation. This threshold is difficult to define as a small variation could lead to the deletion of an additional cluster and conversely. Moreover, this N 1 score is grounded on the labeling system of the clustering algorithm that is not intended to optimize the labeling. As we will see in Section 5, the identification of homologous sulcal pits in the less reproducible clusters (e.g. shallow pits on gyri) is underestimated with this labeling as a consequence of the robustness criteria. Therefore, we fixed a low threshold of 25% that may induce a higher lower bound after the labeling process presented in the next section. However, note that these contentious clusters have a low N 1 score and may often not be enough reproducible for subsequent statistical analysis (e.g. sulcal depth analysis). This leads yet to note that in order to compare atlases (for example between left and right hemispheres), one should rather compare the distributions of the N 1 score than the numbers of atlas basins. In Section 4.3, we will discuss and illustrate a new approach to compare two atlases.

At last, the previous method was grounded on a clustering on the density map of sulcal pits. The shape of this map strongly depends on the smoothing process of the sulcal pits. Our algorithm only extracts the extrema of this map. It is therefore significantly more robust to the smoothing parameters.

Individual pits labeling with varifold comparison

In Section 2, we presented a new automatic method to generate an atlas parcellation of a cortical mesh with respect to a registered population of cortical surfaces. Each region of this atlas, called atlas basin, integrates the average geometry of homologous individual sulcal basins across the population. We can now exploit these geometric similarities in order to optimally assign a subject pit to an atlas basin. For this purpose, we compare individual sulcal basins to atlas basins with a similarity measure M on surfaces. The optimal pairing of these entities will induce the labeling of the sulcal pits. We used an oriented varifold metric which compares the positions, areas, and orientations of the surfaces [Char 13, Kalt 17]. The varifold representation does not depend on the shape parametrization which allows to compare surfaces with a different number of vertices.

Remark 3.1.

1. Although the previous method also provides a labeling system of the individual sulcal pits, the identification is limited to the most robust sulcal pits to improve the quality of the final atlas. Difficult labeling choices are precluded under the notion of isolated pits. Given an atlas of sulcal basins, this second method aims at minimizing the number of such unlabeled pits.

2. This labeling method is independent of the atlas construction, in the sense that given a reference atlas, we can label new subjects that were not included for the atlas construction. This can be convenient for small databases whose size could question the robustness of a tailor-made atlas.

The approach relies on the reproducibility of the sulcal basins. Observation of the data showed that this well-founded hypothesis is yet patently not satisfied in one marginal situation: a subject misses pits in a cortical fold and this fold is represented by a unique large basin (see Fig. 13). This basin covers thus one or more atlas basins. We will therefore preprocess the identification of these basins. Given a registered subject, the identification with the atlas is achieved in two steps presented in Algorithm 3. We first detail these two steps. We denote A the function that returns the area of a surface.

• Step 1: This step handles the case of large individual basins that should not be compared to atlas basins as surfaces, as explained above. For any individual sulcal basin B, we compare its area with the area of the atlas basin B that contains the sulcal pit of B. If A(B) is twice larger than A(B), B is considered as a large basin and its pit inherits the label of B (up to a necessary condition introduced hereafter). During this step, a small number of individual basins are thus associated to atlas basins according to the position of their sulcal pits.

•

Step 2: This is the core of the algorithm. We will estimate for each atlas basin (that have not already been assigned in Step 1) the best geometric match among the individual sulcal basins (that have not already been assigned in Step 1). The most reproducible atlas basins should more likely find a match than small and shallow atlas basins. To estimate the robustness of these atlas basins, we used the maximal sulcal pits density in each atlas basins. Therefore, atlas basins are ordered in a list denoted Q by decreasing peak of sulcal pits density. Each atlas basin of this list is compared to all the unmatched subject basins in order to find the best candidate with respect to the similarity measure M. Once a subject basin is paired to an atlas basin, it can no longer be candidate for a new match.

As the number of pits varies across individuals, we do not expect to retrieve perfect correspondence between individual and atlas basins. This means that when an atlas basin has no anatomical homologous subject basin, its best geometric match is irrelevant. To prevent undesirable pairings, we added a necessary condition grounded on the overlap between two queried basins. This condition is formalized by inequalities of the type A(B ∩ B) > P %A(B) meaning that the intersection of two basins B and B fills at least P % of B, with P a given threshold.

Each labeled individual basin contains a unique sulcal pit that straightforwardly inherits the same label. A novelty of this approach, regarding the methods in the literature, is that a sulcal pit can be assigned to an atlas basin that does not contain it as illustrated in Fig. 14 (third example). However, the overlap criterion always ensures that this atlas basin does overlap subsequently with the pit's basin. Fig. 14 illustrates two other valid matchings and two rejected matchings.

At the end of Step 2, when all atlas basins have been queried for an optimal match, some individual sulcal basins remain unlabeled as well as their sulcal pits. By analogy with the method presented in Section 2, we call them again the isolated pits (see 2.4.3).

Remark 3.2 (Algorithm details). The overlap criterion allows to only compute the varifold metric between admissible pairs of basins. We also added in Step 1 (a) the automatic labeling of individual sulcal basin that strongly overlap with the atlas basin containing the sulcal pit. In fine, Algorithm 1 takes about 2 hours to run by hemisphere for 137 subjects.

Figure 14: Attempts of matching between a subject basin and an atlas cluster. The third matching illustrates the flexibility of our method: the sulcal pit is outside its associated atlas basin.

Algorithm 3: Pits labeling with varifold comparison

Input : Sulcal pits and sulcal basins of an individual registered on an atlas mesh. Atlas parcellation of this mesh in sulcal basins. Output: A set of pairings between the individual sulcal basins and the atlas basins.

Step 1: Order subject basins by decreasing area in a queue Q. For each basin B ∈ Q, consider the atlas basin B that contains the subject pit. Step 2: For each atlas basin B ∈ Q, ordered unassigned subject basins in a new list Q by decreasing distance to B with respect to the similarity measure M. Select the first subject basin B ∈ Q that satisfies the overlap criterion: either

A(B ∩ B) > 50%A(B) or A(B ∩ B) > 50%A(B).

The distance in

Step 2 is the norm of the difference between two basins represented by oriented varifolds. This norm is induced by the following scalar product:

µ X , µ Y = F X i=1 F Y j=1 k pos (x i , y j )k or ( t X i , t Y j )r X i r Y j (5)
where X and Y are the triangulated meshes of two basins B and B, µ X and µ Y are their respective varifold representation that is, e.g. for X, the distribution of the tangent vectors t X i , on the center points (x i ) i of the F X faces of X, and weighted by the areas (r X i ) i of these faces. In the experiments, we will use two Gaussian kernels

k pos (x, y) = e -|x-y| 2 σ 2 , k or ( t i , t j ) = e 2( t i • t j ) σ 2 s , ( 6
)
where the main scale parameter is set to σ = 15mm and the other one, on the angular sensitivity, is set to σ s = 0.5.

Experiments and results

In order to evaluate our pipeline (summarized in Fig. 3), we first used the adult database studied in [Auzi 15]. We will compare their results with ours. We will then present a second application to a children database with a large range of ages and we will show that we provide reliable atlases and inter-subject labeling on these two very different populations, with the same parameters. More precisely, we computed for each population and each hemisphere a specific atlas with the first part of our pipeline. We then labeled each subject according to their respective atlas with the varifold labeling method. This labeling will be evaluated at the group level through the N 1 score (the percentage of subjects that have a pit associated to the basin) and the density maps of homologous sulcal basins. We recall that subjects have no more than one sulcal pits associated to an atlas basin. Few sulcal pits remain unlabeled (the isolated pits).

Data acquisition and preprocessing

OASIS database: This data was used in [Auzi 15] and consists in 137 right-handed young adult healthy subjects, aged from 18 to 34, selected from the Open Access Series of Imaging Studies (OASIS) database (www.oasis-brains.org) [ Marc 07]. For each subject, three to four individual T1-weighted magnetization prepared rapid gradient echo (MP-RAGE) scans were obtained on a 1.5 T Vision system (Siemens, Erlangen, Germany) with the following protocol: in-plane resolution = 256 × 256 (1 mm × 1 mm), slice thickness= 1.25 mm, TR= 9.7 ms, TE = 4 ms, flip angle = 10u, TI = 20 ms, TD = 200 ms. Images were motion corrected and averaged to create a single image with a high contrast-to-noise ratio [ Marc 07]. Anatomical MR images were processed using Freesurfer v5.1.0, in order to extract the inner cortical surface mesh, surface area, intracranial volume and obtain spherical interindividual correspondences [Dale 99, Fisc 12]. CMIND database: 155 healthy children, between 9 months and 18-years-old, were selected form the Cincinnati MR Imaging of Neurodevelopment (C-MIND) project (04/30/2015 release). Among them 55 subjects were scanned longitudinally with two or three timesteps. For each timestep, one or two T1w images were acquired on a Philips 3T MRI scanner with the following protocol: matrix size = 256 × 224, spatial resolution = 1 mm × 1 mm × 1 mm, TR/TE = 8.1/3.7 ms, flip angle = 8 • , TI = 939 ms. These images were denoised using SPM8 (SANLM algorithm) and a grey/white matter segmentation was performed using FreeSurfer 5.3.0. For 10 of these subjects, manual correction using the control points of FreeSurfer was applied.

The extraction of sulcal pits was performed using the procedure of [Auzi 15], designed to yield reproducible sulcal pits especially in children populations (see Section 1.1.1). Correspondences between cortical meshes were then obtained by spherical alignment constrained by convexity information using Freesurfer [Fisc 12]. We used a single average mesh estimated with the OASIS database. To remove the non-cortex part of this mesh (cingular mask), we summed all the individual registered cingular masks and we extracted the vertices that belong to more than 70% of the population. Each individual map showing the localization of the sulcal pits was then projected onto a template surface, after smoothing with an FWHM of 5 mm, maintaining a peak value of 1 [Auzi 15].

Results on OASIS

On average, subjects have 88.3 ± 4.7 sulcal pits on the left hemisphere and 89.5 ± 5.1 on the right hemisphere. We identified 92 clusters on the left hemisphere and 90 clusters on the right hemisphere. Note that the numbers of clusters (or atlas basins) well correspond to the average number of sulcal pits. The old atlas has 104 clusters on the left hemisphere and 114 on the right hemisphere. With the old pipeline, a subject can have more than one sulcal pit in the same cluster. In order to compare the sulcal pit's labeling of the two pipelines, in case of multi-labeling with the method of [Auzi 15], we selected, for any following analysis, the closest pit to the peak of pits' density in the cluster. The old pipeline identified for each subject about 72.8 sulcal pits on the left hemisphere and 76.7 on the right hemisphere. The new pipeline identified respectively 72.7 and 72.9 sulcal pits. Despite the large reduction of the number of clusters, the identification of the sulcal pits is thus quantitatively similar (with 24 additional clusters, the old pipeline only labels 4 additional sulcal pits on the right hemisphere). See Tab. 2.

The percentages of sulcal pits associated to each atlas basin across the population, i.e. the N 1 scores, are displayed in Fig. 15. In overall, we retrieved mostly the same basins than [Auzi 15] but most basins with a low N 1 score no longer appear (blue basins). In order to make a fair comparison of the N 1 scores on the most reliable clusters, we excluded the less reproducible clusters of the old method. For the sake of simplicity, we aligned the number of clusters of each method, i.e. we selected for the old pipeline the 92 and 90 best N 1 scores on the left and right hemispheres. We then computed the average N 1 scores of the two methods. Results are reported in Tab. 2. Despite this selection, the new pipeline still presents the best scores and with a smaller variance. The specificity of our method is to generate atlas basins that are geometrically good representatives of the individual sulcal basins. Regarding this goal, we can observe that the boundaries of our atlas basins are very well aligned with the geometry of the atlas mesh which suggests that we indeed retrieved the average boundaries of sulcal basins. Conversely, purple arrows in Fig. 15 highlight atlas basins that cover gyri. Typically, the boundary of a basin, that should be aligned with the crest line, is overextended in the adjacent sulcus.

Figure 15: (OASIS) Final N 1 scores, percentage of subjects that have a pit associated to a given cluster, after the complete pipeline. First row: old method. Second row: new method.

The influence maps introduced in this paper are a powerful tool to analyze the geometry of sulcal basins at the group level. Let us recall that for any atlas basin, its influence map is the density map of its associated individual sulcal basins (see Fig. 6). In order to evaluate our complete pipeline, we estimated these maps with respect to the labeling of the varifold method presented in Section 3. As we will see, these maps give us a qualitative evaluation of the consistency of the atlas basins and an evaluation of the labeling.

In Fig. 16 (a), we illustrate the flexibility of the varifold matching. As the labeling is mostly driven by the shape of sulcal basins, sulcal pits associated to a cluster may be located outside this cluster. In most cases, these outliers are close to the boundary of the cluster, but in some highly variable regions like the temporal lobe, their distribution can be more spread out. We can then either individually verify the labeling or we can directly consider the influence map of this atlas basin as displayed in the background of Fig. 16 (a). This density map reaches 100% in the middle In order to visualize all the influence maps (at least partially), we display on a single mesh in Fig. 17 the restriction of each map to its respective atlas basin and we compare them to the influence map generated by the previous atlas of [Auzi 15]. We reported the purple arrows of Fig. 15. These arrows now all point blue areas where the density vanishes. This confirms that the respective clusters are too large and that the associated sulcal basins do not cover these areas. For a finer analysis, we denoted two clusters in the superior temporal sulcus (Fig. 15, top left, arrows 1 and 2). The density in Cluster 2 remains very high on the left boundary. This indicates that most of the associated basins overextend the cluster to the left. This is also confirmed by the low density in Cluster 1 on its right boundary. Most of the associated basins do not reach this boundary that should be slightly shifted towards the left as it is corrected in our new atlas. The misplacement of this boundary seems to induce labeling errors. Indeed, the right end of the pits distribution in Fig. 16 (a) belongs to the blue area on the right side of Cluster 1. These pits should yet be assigned to Cluster 2. This explains why the N 1 score of Cluster 2 is smaller than its equivalent cluster in the new atlas and conversely why the N 1 score of Cluster 1 is higher than its equivalent in the new atlas (see Fig. 15).

With the new pipeline, most of the clusters are significantly filled in red which indicates that they are geometrically good representatives of the individual basins. The peaks of density in each atlas basin is in overall significantly higher than with the old method. This means that the overlap between homologous sulcal basins is more important which suggests that the old pipeline provides less robust sets of homologous sulcal basins. A small decrease of the density at proximity of the boundaries ensures that these boundaries are well defined (thin yellow auras on the boundaries of each atlas basins). The blue-green areas illustrate regions where a minority of subjects present a different basin pattern. Most likely for these subjects, an additional cluster is present in these areas. We recall at last that these influence maps are normalized. Consequently, some clusters may have a low N 1 score and a high influence map. This means that they are good representatives of the subset of the population that have a basin at this place. Let us note that the notion of variability maps introduced in [Krug 18] is in the same spirit that our fusion of density maps and could also be investigated to evaluate the quality of an atlas.

To go further on the comparison of the two methods, we investigated the depth of sulcal pits, estimated with the depth potential map. In Fig. 18, we display for each pipeline the average depth in deep atlas basins. The new atlas tends more to the blue and is more homogeneous. This suggests that we successfully excluded noisy pits. More generally, Fig. 19 displays the standard deviation (STD) of the sulcal pits' depth. With the exception of shallowest basins where the variability can result from a large set of factors, the STD is globally smaller with the new atlas.

Figure 18: (OASIS) Mean of sulcal pits' depth of the deepest clusters (with a threshold at -1). Old pipeline (first row) and new pipeline (second row). Note that, on the whole cortical surface, the depth potential map (DPF) takes values between -2 and 0.3 with negative values for the deepest points.

We conclude this study on the OASIS database with the isolated pits. They are the unidentified sulcal pits at the end of the complete pipeline. They should result from noise or minor sulcal patterns. We computed the density map of these isolated pits to verify that their final number remains small with respect to the total number of sulcal pits and that their distribution is spread. We do observe some low peaks of density that indicate the location of some minor patterns in the distribution of these pits at the population level (see Fig. 20).

Results on CMIND

This database presents 97 ± 6.9 sulcal pits by subject on the left hemisphere and 99.2 ± 6.8 on the right hemisphere. We identified 100 atlas basins on the left hemisphere and 102 atlas basins on the right hemisphere. Regarding the labeling, our pipeline identified for each subject about 80.4 Right column: Density map of all sulcal pits. The (common) colormap has been saturated above 25% to emphasize the density of isolated pits. Each map is normalized by the total number of subjects.

± 4.3 sulcal pits on the left hemisphere and 79.8 ± 3.7 on the right hemisphere. The average N 1 score is 78.0 ± 18.1 % on the left hemisphere and 80.6 ± 16.8 % on the right hemisphere. See Tab. 3.

We recall that we used the same parameters for both populations (OASIS and CMIND) and both hemispheres. Despite the important differences regarding the average number of sulcal pits by subject between the two databases, our pipeline provided in each case a number of atlas basins equivalent to the number of sulcal pits. Moreover, the average and variance of the N 1 scores, that can be seen as global performance indicators, are also very similar.

The N 1 scores of the atlas basins are presented in Fig. 21. As for the adult database, only few atlas basins have a low N 1 (light blue). Once again, the boundaries of the atlas basins are very well aligned with the geometry of the atlas mesh which suggests that we indeed retrieved the average boundaries of sulcal basins. This is confirmed by the restriction of the influence maps to each atlas basin presented in Fig. 21 second row. Most of the basins are again significantly filled in red which indicates that they are geometrically good representatives of the individual basins. This also validates the labeling at the group level.

Atlas comparison

We end this paper with a small discussion on atlas comparison from a methodological viewpoint. We do not intend to perform a neurological comparison between the OASIS and CMIND atlases as this goes beyond the scope of this paper. Note however that these two databases illustrate that sulcal pits can depend on the acquisition and preprocessing of cortical meshes. Indeed, adult subjects of the first database have about 10% less sulcal pits than the younger subjects of the second database (age range from 0 to 18). This observation is not consistent with the fact that the number of sulcal pits in this last database increases with age and this emphasizes that comprehensive analyses of sulcal pits should be grounded on tailored atlases of sulcal pits for any new database. This being said, we will see that the density maps of sulcal basins are an efficient tool to investigate the differences between two atlases. The comparison of two atlases cannot be reduced to the number of clusters. The presence of clusters with a low N 1 scores is not robust. Indeed, they strongly depend on the threshold parameter p of Algorithm 2 (see Annex). To illustrate this point, we selected three examples where the CMIND atlas presents an additional cluster. In Fig. 22, 23, and 24, the red cross-hair points this additional cluster and the corresponding area on the OASIS atlas. On the first example, we can observe that the density maps of the two neighboring clusters in the OASIS atlas is particularly low (40%) at the position of the cross. This means that 60% of associated sulcal basins of these clusters do not reach this cross-hair. We can deduce that there is a similar minor pattern in the OASIS database of an additional cluster in this area. Most likely, this minor pattern was closely under the threshold value. On the second example, the density remains high (60%) which suggests a significant difference that should be deepen. For these two arguable examples, the N 1 scores of the additional cluster are low (51% and 49%). On the last example (Fig. 24), the density at the location of the missing cluster is very high (76%). This means that 76% of the subjects in the OASIS database have a single large basin in this area that fills the top of this fold. Therefore, if we ignore the variability of sulcal basins that remains after the registration and that could partially explain the density decrease, there are at most 24% of the subjects with an additional sulcal pit around the cross-hair. On the contrary, in the CMIND atlas, the N 1 score of this additional cluster is very high (78%). This confirms an important difference between the two atlases regarding this cluster.

Conclusion

In [Auzi 15], authors showed that the notion of sulcal pit should not be reduced to the deepest areas of the cortical surface. This approach doubles the number of these landmarks. The resulting highly complex organization of the sulcal pits required the new methods that we presented in this paper.

The first process, presented in Section 2, builds an atlas of the sulcal organization in a groupwise parcellation of the cortical surface. This parcellation results from a derivation of the watershed by flooding algorithm where density maps of sulcal basins play the role of local depth maps. We then addressed the inter-subject labeling problem with a geometrical viewpoint and independently from the atlas problem. Besides the synergy between the two methods presented in Sections 2 and 3, the second process can be used with any given atlas allowing to label additional subjects or small databases. We successfully applied the complete pipeline to a highly heterogeneous children database. The quantitative performances were similar to those obtained on an adult database.

Besides the geometric viewpoint, a novelty of both processes is to take advantage of the oneto-one correspondences assumption between sulcal pits across subjects. As a consequence, we introduced the notion of isolated sulcal pits. These are the sulcal pits that cannot be identified by the atlas as they are not reproducible enough across the population.

Our new atlases are the first insight in the literature on the average shape of sulcal basins after registration. The N 1 score and the merged density maps of sulcal basins validated the global assumption of one-to-one inter-subject correspondences of sulcal pits and sulcal basins on a major part of the brain. These density maps of sulcal basins, namely the influence maps, are a highly efficient tool to evaluate the labeling at the group level. They also seem to be a key to understand where the inter-individual variability is locally more complex. In future work, they could be used to validate and to illustrate patterns of sulcal basins in addition to the pattern analysis of sulcal pits as [Meng 18] 's work.

• OASIS: Cross-Sectional: Principal Investigators: D. Marcus, R, Buckner, J, Csernansky J. Morris; P50 AG05681, P01 AG03991, P01 AG026276, R01 AG021910, P20 MH071616, U24 RR021382

• the C-MIND Data Repository created by the C-MIND study of Normal Brain Development. This is a multisite, longitudinal study of typically developing children from ages newborn through young adulthood conducted by Cincinnati Children's Hospital Medical Center and UCLA and supported by the National Institute of Child Health and Human Development (Contract #s HHSN275200900018C, release 010715). A listing of the participating sites and a complete listing of the study investigators can be found at https://research.cchmc.org/cmind. This manuscript reflects the views of the authors and may not reflect the opinions or views of the NIH.

Annex

Varifold labeling

In order to evaluate the improvement of the labeling with the second part of the pipeline (the labeling with varifolds), we compared the N 1 score between the first labeling implicitly produced during the atlas construction and the final labeling. The N 1 increase is in average on both hemispheres 14.0 % for CMIND. However, this increase is not uniform and can be significantly higher in less robust basins as illustrated in Fig. 25. This figure also highlights a weakness of the geometric criterion regarding the sulcal basins whose pits are close to the boundary. In this case, the seed of the cluster is also close to the boundary and some individual basins that are slightly smaller than the cluster may not contain this seed. Although they may still significantly overlap with the cluster, the geometric criterion would reject them. This partially explains the important increase of sulcal basins associated to the three clusters with a red marker in Fig. 25.

Figure 25: (CMIND) Increase of the N 1 score between the constrained labeling of the method to estimate the atlas and the optimized labeling with varifolds.

Threshold of Algorithm 2

Let us recall that to finalize the atlas basins, we iteratively delete the less robust basins. In Remark 2.3, we explained that to produce a robust stop criterion, we consider at each iteration the 5 minimal N 1 scores of the clusters and we stop the procedure if the average of these scores is above a threshold set to p = 25%. Fig. 26 displays this average minimal score at each iteration. The choice of the threshold remains open to discussion as the slopes of these curves are quite low. However, the affected clusters by this threshold are not robust clusters and should not significantly impact the final labeling system of the sulcal pits. An interesting fact with the current threshold p = 25% is that the final numbers of clusters match the respective average numbers of sulcal pits in the four hemispheres of the OASIS and CMIND populations. We investigated few other criteria grounded on the maxima of density in the distribution of isolated pits, the ratio between the area of the cluster and the average area of its associated individual basins, the average conflict on a cluster, etc.

[

  Auzi 15] uses instead the Depth Potential Function (DPF) introduced in [Bouc 09]. The DPF is

Figure 1 :

 1 Figure 1: Standard watershed algorithm: merging of the yellow basin with the blue basin during flooding. The merging is due to the small ridge height illustrated in the middle by the red arrow.

Figure 2 :

 2 Figure 2: Geometric issues with the clusters resulting from a watershed by flooding algorithm on the distribution of sulcal pits.

  Construction of atlas basins with a watershed by flooding on multiple density maps of sulcal basins Algorithm 2: Filtering of spurious atlas basins

Figure 3 :

 3 Figure 3: Overview of the complete pipeline referring to Algorithm 1, 2 and 3.

Figure 4 :

 4 Figure 4: Density maps of sulcal basins playing the role of local depth maps for a watershed by flooding algorithm.

Figure 5 :

 5 Figure 5: From left to right: Map showing the density of sulcal pits. Clusters initialized with a small neighborhood around the seeds before and after basic filtering. Red arrows highlight clusters that overlap.

Figure 6 :

 6 Figure 6: From left to right: Concentration of sulcal pits in an initial cluster. Extraction of individual sulcal basins associated to the cluster's pits. Influence map of the cluster.

Figure 8 :

 8 Figure 8: Growth of a cluster during the adaptive watershed algorithm. Second row: Pits' concentration. Third row: Influence map of the cluster.

Figure 9 :

 9 Figure9: Six influence maps before (left) and after normalization (right) with their respective cluster in black (the maps are superimposed on a single mesh but their support can overlap). These black clusters will then preferentially grow where they have a high influence.

Algorithm 2 :

 2 Cluster's filtering Input : List of n clusters, their seed and their associated individual sulcal pits (the sets (Subj k ) k )). Threshold p. Output: The final atlas basins, which form a parcellation of the atlas mesh. Auxiliary functions: Cluster deletion Deletion fast version { Input : Index k D of the cluster to delete. • Reset the vertices in the cluster B A (k D ) and delete the seed s k D . • Run the Watershed function of Algorithm 1. } Deletion complete version { Input : Index k D of the cluster to delete.

Figure 12 :

 12 Figure 12: Example of one cluster deletion computed with the function Deletion complete version. From left to right: clusters before and after the deletion. The algorithm does not induce a fusion between two basins. The deleted cluster is here divided among the others. With Deletion fast version, only vertices in the purple cluster could inherit a new label. Note here that the right part of the boundary between the orange and the yellow clusters has slightly moved which would not have been possible with fast version.

Figure 13 :

 13 Figure 13: Example of a subject with an unusual large basin in the superior temporal sulcus and two others on the intern face.

  (a) If A(B ∩ B) > 80%A(B), the pit is assigned to B and B is removed of Q. (b) If the area of B is twice larger than the area of B, and if A(B ∩ B) > 50%A(B), the pit is assigned to B and B is removed of Q.

Figure 16 2 Figure 17

 16217 Figure 16: (OASIS) (a) Distribution of sulcal pits for a label of the temporal sulcus on top of the density map of the respective sulcal basins after the complete pipeline. (b) Density map of the individual sulcal basins for the label corresponding to the middle atlas basin of the central sulcus.

Figure 19 :

 19 Figure 19: (OASIS) Standard deviation of sulcal pits' depth for each cluster. Old pipeline (first row) and new pipeline (second row).

Figure 20 :

 20 Figure 20: (OASIS) Left column: Density map of isolated sulcal pits after the complete pipeline.Right column: Density map of all sulcal pits. The (common) colormap has been saturated above 25% to emphasize the density of isolated pits. Each map is normalized by the total number of subjects.

Figure 21 :

 21 Figure 21: (CMIND) Final N 1 score (first row) and fusion of the local influence maps (second row) after the complete pipeline for the right and left hemispheres.

Figure 22 :

 22 Figure22:The N 1 score of the additional atlas basin in CMIND is 51%.OASIS CMIND

Figure 26 :

 26 Figure26: Average minimal N 1 after each iteration. Horizontal axis is the number of current clusters at each iteration. The horizontal line illustrates the threshold. Consider for example the left hemisphere in OASIS (top plot, red marks). At the beginning of the process, the atlas contains 158 clusters for which the mean of the 5 minimal N 1 scores equals about 10% (left side of the plot). The cluster deletion stops when this mean exceeds 25% and the atlas contain at this stage 92 clusters.

Table 1 :

 1 Number of atlas basins, number of labeled sulcal pits and comparison with the old pipeline.

			New pipeline	Old pipeline
			# Atlas # Labeled # Pits # Labeled
	Hemisphere # Sulcal pits	basins	pits	clusters	pits
	Left	88.3 ± 4.7	92	72.7	104	72.8
	Right	89.5 ± 5.1	90	72.9	114	76.7

Table 2 :

 2 Average N 1 score.

		New pipeline	Old pipeline
			Filtered	# Excluded
	Hemisphere	Average N 1	average N 1	clusters
	Left	79.0 ± 17.4 % 76.8 ± 21.8 %	12
	Right	81.0 ± 18.0 % 78.3 ± 18.6 %	24

Table 3 :

 3 Number of atlas basins, number of labeled sulcal pits and average N 1 score. Hemisphere # Sulcal pits # Atlas basins # Labeled pits Average N 1

	Left	97 ± 6.9	103	80.4 ± 4.3	78.0 ± 18.1 %
	Right	99.2 ± 6.8	99	79.8 ± 3.7	80.6 ± 16.8 %

Acknowledgments

This work was funded by the Fondation de France (Grant #2015-00059546). Data used in the preparation of this article were obtained from

Connectedness of the clusters

Finally, let us return on Algorithm 1. In order to ensure that every cluster B A (k) remains connected, we consider a penalty term that modifies the influence

where d(x, B A (k)) equals 0 if x is an unlabeled vertex at the boundary of B A (k) and +∞ otherwise.

In practice, we introduce a set of connectedness penalty functions CP k : Algorithm 4: Adaptive watershed algorithm with a connectedness constraint Input : Subjects' sulcal pits and sulcal basins registered on an atlas mesh.

Output: The atlas basins, which form a parcellation of the atlas mesh. A labeling system with respect to this parcellation (given by the sets (Subj k ) k ).

Initialization:

• Compute the density map of pits.

• Extract the seeds, which are the maxima of the density map.

• Initialize the clusters by the 2-ring neighborhood of the seeds.

• Initialize the respective lists Subj k .

Watershed: Initialize all the CP k to 0.

1. Compute the list Q ⊂ X of unlabeled vertices at the boundary of the clusters. If Q is empty, the atlas is complete. Otherwise:

2. Using the lists (Subj k ) k , compute the influence map of each cluster and extract the vertices under maximal influence:

3. Select among them a vertex x with minimal conflict C(x).

If x belongs to the neighborhood of B

Update Subj kx , reset all the CP k to 0 and return to step (1).

5. If CP kx (x) = 0, then update it to CP kx (x) = I kx (x) and return to step (1).

6. If CP kx (x) = 0, x inherits the most frequent label in its 1-ring neighborhood. Update Subj kx , reset all the CP k to 0 and return to step (1).


Valérie Laperche
email: v.laperche@brgm.fr

Gabriel Billon

Bruno Lemiere

Ludovic Lesven

Pierre Jean Superville

Marija Margus

Claire Alary

Olivier Bour

Florian Liénard

Pascal Auger

Bruno Lemière

Mathieu Henry

Environmental monitoring at a sediment source site to qualify for reuse

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction:

The Interreg FWVL VALSE project aims to promote sustainable dredged sediment reuse [1]. Reuse options and supporting research were investigated during an earlier project [2]. In order to provide waterways operators and related industries, with case studies, full scale pilot tests are planned at Belgian and French sites. A former sediment disposal site of the French Waterways (VNF) at Saint Omer, Northern France, is considered as a potential source of mature sediment (14 to 40 years old) for reuse applications (landscaping, cycle paths, dikes). Sediments are now comparable to soil in terms of water content and mechanical behaviour. Sediment environmental characteristics are assessed with laboratory methods, and its variability with on-site methods [3]. Methods: The site survey was completed using a portable XRF (pXRF) to test the presence potential inorganic contaminants in stored sediments and their spatial variability. A new field colorimetric method was used for Acid Volatile Sulfides. Existing boreholes allowed a shallow groundwater survey for inorganic contaminants (water sampling, field electrochemical techniques). This survey is aimed at testing field measurements for effective site monitoring. Sediments were analysed similarly to raw soil with Olympus Delta Premium and Niton XL3t980 pXRFs at < 15 m intervals. No dehydration by filter press was needed to reach standard water contents (25-35%). In this range, measurements carried out on wet pressed pellets correlate well with laboratory analyses [3]. Results: The resulting map (Fig. 1) suggests that sediments are generally homogeneous for inorganic contaminants, with Pb-Zn-Cr levels similar to target site soils. One exception is the northwesternmost storage cell, with moderately higher values. In the sediments, the S-AVS contents measured in the field are quite low, even under the limit of detection (< 50 ppm) in most cases. Electrochemical measurements (EM) were able to detect Zn after a short acid attack in the field. Samples were analysed in the laboratory with further optimisation of EM, Cu and Pb were detected as well.

In porewater, EM on solid electrodes indicate the presence of metals at low concentrations, most of them under the detection limit (< 0.5 ppb). Further optimisation is required due to the complex matrix. Discussion: Site-specific threshold values for pXRF measurements can be derived from regulatory limits or risk-based targets after calibration. Measurements will be used during works to identify potentially offspecification sediments in quasi real time. pXRF measurements will be used at each reuse application site, prior to and after works, in order to monitor application impacts. The low AVS contents indicate that the sediment has been re-oxydized, and would promote a progressive change in metal speciation with time. On-site electrochemical analyses for contaminants in groundwater offer flexibility and rapidity when monitoring works impact, and provide an improved response for public awareness.

Fig. 1 :

 1 Fig. 1: Map of Pb measurements by pXRF (mg/kg).

  Acknowledgements: This research is funded by InterReg FWVL, project VALSE. References: [1] Alary, C. et al., (2017) SedNet conference; [2] Alary, C. et al., (2015) SedNet conference; [3] Lemière, B. et al., (2014) GEEA 22:222-233.