dataesr/scientific-paragraphs-categorization huggingface dataset:
The Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework  (Trouvé, 1998; Beg et al. , 2005; Younes, 2010)  provides a natural setting for quantifying deformations between shapes or images. This framework provides diffeomorphic transformations which preserve the topology and also provides a metric between shapes. The LDDMM framework is also a natural setting for estimating templates from a population of shapes, because such templates can be defined as means in the induced shape space. Various methods have been proposed to estimate templates of a given population using the LDDMM framework  (Vaillant et al. , 2004; Glaunès and Joshi, 2006; Durrleman et al. , 2008; Ma et al. , 2008). All methods are computationally expensive due the complexity of the deformation model. This is a limitation for the study of large databases.
In this paper, we present a fast approach for template-based statistical analysis of large datasets in the LDDMM setting, and apply it to a population of 1,000 hippocampal shapes. The template estimation is based on iterative diffeomorphic centroid approaches, which two of them (IC1 and IC2) were introduced at the Geometric Science of Information conference GSI13  (Cury et al. , 2013) , while the PW approach was introduced in a chapter of the Geometric Theory of Information book  (Cury et al. , 2014b). The advantage of such method is the possibility of parallelising the process, easily adding new subjects and saving computation time without losing accuracy in the template estimation. Furthermore, the iterative nature of the methodology is useful when data are being continuously collected or preprocessed.
Here, we propose to use these centroid estimations directly for template-based statistical shape analysis. The analysis is done on the tangent space to the template shape, either directly through Kernel Principal Component Analysis (Kernel PCA  Schölkopf et al. , 1997)  or to approximate distances between subjects. We perform a thorough evaluation of the approach using three datasets: one synthetic dataset and two real datasets composed of 50 and 1,000 subjects respectively. In particular, we study extensively the impact of different centroids on statistical analysis, and compare the results to those obtained using a standard variational template method. We will also use the large database to predict, using the shape parameters extracted from a centroid estimation of the population, some anatomical variations of the hippocampus in the normal population, called Incomplete Hippocampal Inversions and present in 17% of the normal population  (Cury et al.
Remarks on the IC Algorithms These three methods depend on the ordering of subjects. In a previous work  (Cury et al. , 2014b) , we showed empirically that different orderings result in very similar final centroids. Here we focus on the use of such centroid for statistical shape analysis. Regarding memory performances of the algorithm, computing a mean shape with the algorithm IC1 on meshes composed by 500 vertices, requires 1. 2 MB. For the algorithm IC2, which combines surfaces at each iterations, requires more memory, for 50 meshes IC2 needs 20 MB. We used the real dataset RD50 (see section 5. 1 for further details) to estimate memory performances.
A second strategy consists in analysing the set of pairwise distances between subjects. Then, the distance matrix can be entered into analysis methods such as Isomap  Tenenbaum et al. (2000) , Locally Linear Embedding  Roweis and Saul (2000) ,  and Yang X. et al. (2011)  or spectral clustering algorithms  (Von Luxburg, 2007). Here, we tested two approaches: (i) the analysis of initial momentum vectors using a Kernel Principal Component Analysis for the first strategy; (ii) the approximation of pairwise distance matrices for the second strategy. These tests allow us both to validate the different iterative centroid methods and to show the feasibility of such analysis on large databases.
Principal Component Analysis on Initial Momentum Vectors The Principal Component Analysis (PCA) on initial momentum vectors from the template to the subjects of the population is an adaptation of PCA in which Euclidean scalar products between observations are replaced by scalar products using a kernel. Here the kernel is K V the kernel of the R. K. H. S V. This adaptation can be seen as a Kernel PCA  (Schölkopf et al. , 1997). PCA on initial momentum vectors has previously been used in morphometric studies in the LDDMM setting  (Vaillant et al. , 2004; Durrleman et al. , 2009)  and it is sometimes referred to tangent PCA. We briefly recall that, in standard PCA, the principal components of a dataset of N observations a i ∈ R P with i ∈ {1,. , N} are defined by the eigenvectors of the covariance matrix C with entries: C(i, j) = 1 N -1 (a i -ā) t (a j -ā) (21) with a i given as a column vector, ā = 1 N N i=1 a i , and x t denotes the transposition of a vector x.
Distance Matrix Approximation Various methods such as Isomap  (Tenenbaum et al. , 2000)  or Locally Linear Embedding  (Roweis and Saul, 2000; Yang X. et al. , 2011)  use as input a matrix of pairwise distances between subjects. In the LDDMM setting, it can be computed using diffeomorphic distances: ρ(S i , S j ) = D(id, ϕ ij ). However, for large datasets, computing all pairwise deformation distance is computationally very expensive, as it involves O(N 2 ) matchings. An alternative is to approximate the pairwise distance between two subjects through their matching from the centroid or template. This approach has been introduced in Yang X. F. et al. (2011). Here we use a first order approximation to estimate the diffeomorphic distance between two subjects: ρ(S i , S j ) = α j (0) -α i (0) , K V (x(0))(α j (0) -α i (0) , ( 25 ) with x(0) the vertices of the estimated centroid or template and α i (0) is the vector of initial momentum vectors computed by matching the template to S i.
EXPERIMENTS AND RESULTS In this section, we evaluate the use of iterative centroids for statistical shape analysis. Specifically, we investigate the centring of the centroids within the population, their impact on population analysis based on Kernel PCA and on the computation of distance matrices. For our experiments, we used three different datasets: two real datasets and a synthetic one. In all datasets shapes are hippocampi. The hippocampus is an anatomical structure of the temporal lobe of the brain, involved in different memory processes.
Data The two real datasets are from the European database IMAGEN  (Schumann et al. , 2010)  1 composed of young healthy subjects. We segmented the hippocampi from T1-weighted Magnetic Resonance Images (MRI) of subjects using the SACHA software  (Chupin et al. , 2009)  (see Figure  2 ). The synthetic dataset was built using deformations of a single hippocampal shape of the IMAGEN database.
The Synthetic Dataset SD SD is composed of synthetic deformations of a single shape S 0 , designed such that this single shape becomes the exact center of the population. We will thus be able to compare the computed centroids to this exact center. We generated 50 subjects for this synthetic dataset from S 0 , along geodesics in different directions. We randomly chose two orthogonal momentum vectors β 1 and β 2 in R 3×n. We then computed momentum vectors α i ,i ∈ {1,. , 25} of the form k i 1 β 1 + k i 2 β 2 + k i 3 β 3 with (k i 1 , k i 2 , k i 3 ) ∈ R 3 , ∀i ∈ {1,. , 25},k i j ∼ N (0, σ j ) with σ 1 > σ 2 ≫ σ 3 and β 3 a randomly selected momentum vector, adding some noise to the generated 2D space. We computed momentum vectors α j ,j ∈ {26,. , 50} such as α j = -α j-25. We generated the 50 subjects of the population by computing geodesic shootings of S 0 using the initial momentum vectors α i ,i ∈ {1,.
The Real Dataset RD50 RD50 is composed of 50 left hippocampi from the IMAGEN database. We applied the following preprocessing steps to each individual MRI. First, the MRI was linearly registered toward the MNI152 atlas, using the FLIRT procedure  (Jenkinson et al. , 2002)  of the FSL software 2. The computed linear transformation was then applied to the binary mask of the hippocampal segmentation. A mesh of this segmentation was then computed 1 http://www. imagen-europe. com/ 2 http://fsl. fmrib. ox. ac. uk/fsl/fslwiki/FslOverview from the binary mask using the BrainVISA software 3. All meshes were then aligned using rigid transformations to one subject of the population. For this rigid registration, we used a similarity term based on measures (as in  Glaunès et al. , 2004). All meshes were decimated in order to keep a reasonable number of vertices: meshes have on average 500 vertices.
The Real Database RD1000 RD1000 is composed of 1,000 left hippocampi from the IMAGEN database. We applied the same preprocessing steps to the MRI data as for the dataset RD50. This dataset has also a score of Incomplete Hippocampal Inversion (IHI)  (Cury et al. , 2015) , which is an anatomical variant of the hippocampus, present in 17% of the normal population.
Experiments For the datasets SD and RD50 (which both contain 50 subjects), we compared the results of the three different iterative centroid algorithms (IC1, IC2, and PW). We also investigated the possibility of computing variational templates, initialized by the centroids, based on the approach presented in section 3. 5. We could thus compare the results obtained when using the centroid directly to those obtained when using the most expensive (in term of computation time) template estimation. We thus computed 6 different centres: IC1, IC2, PW and the corresponding variational templates T(IC1), T(IC2), T(PW). For the synthetic dataset SD, we could also compare those 6 estimated centres to the exact centre of the population. For the real dataset RD1000 (with 1,000 subjects), we only computed the iterative centroid IC1.
For all computed centres and all datasets, we investigated: (1) the computation time; (2) whether the centres are close to a critical point of the Fréchet functional on the manifold discretised by the population; (3) the impact of the estimated centres on the results of Kernel PCA; (4) their impacts on approximated distance matrices. Computation times are given for the different template estimation methods. For the variational template initialized with an estimated centroid, the computation time depends on the number of iteration needed to converge. We fixed this number to 5 iter (see section 3. 5), which is usually, for the kind of shape we are using here, the number of iteration needed by the variational template to converge. For the variational template with standard initialisation (all objects of the population are considered, see section 3.
Now a null ratio R would mean that N i=1 v i (•, 0) = 0, which is precisely the equation satisfied by critical points of the Fréchet functional. In practice, this ratio is expected to be close to zero but not null due to matchings inaccuracies. We compared the results of Kernel PCA computed from these different centres by comparing the principal modes and the cumulative explained variance for different number of dimensions. Finally, we compared the approximated distance matrices to the direct distance matrix. For the RD1000 dataset, we will try to predict an anatomical variant of the normal population, the Incomplete Hippocampal Inversion (IHI), present in only 17% of the normal population.
"Centring" of the Estimated Centres Since in practice a computed centre is never at the exact centre, and its estimation may vary accordingly to the discretisation of the underlying shape space, we decided to generate another 49 synthetic populations as detailed in section 5. 1, so we have 50 different discretisations of the shape space. For each of these populations, we computed the 3 centroids and the 3 variational templates initialized with these centroids. We calculated the ratio R described in the previous section for each estimated centre. Table  1  presents the mean and standard deviation values of the ratio R for each centroid and template, computed over these 50 populations. In a pure Riemannian setting (i. e. , disregarding the fact that matchings are not exact), a zero ratio would mean that we are at a critical point of the Fréchet functional, and under some reasonable assumptions on the curvature of the shape space in the neighborhood of the dataset (which we cannot check however), it would mean that we are at the Fréchet mean.
By construction, the ratio computed from the exact centre using the initial momentum vectors α i used for the construction of subjects S i (as presented in section 5. 1) is zero. Ratios R are close to zero for all centroids and variational templates, indicating that they are close to the exact centre. Furthermore, the value of R may be partly due to the nonexactitude of the matchings between the estimated centres and the subjects. To become aware of this non-exactitude, we matched the exact centre toward all subjects of the SD dataset. The resulting ratio is R = 0. 05. This is of the same order of magnitude as the ratios obtained in Table  1 , indicating that the estimated centres are indeed very close to the exact centre.
PCA on Initial Momentum Vectors We performed a PCA computed with the initial momentum vectors (see section 4 for details) from our different estimated centres (3 centroids, 3 variational templates and the exact centre). We computed the cumulative explained variance for different number of dimensions of the PCA. Results are presented in Table  2. The cumulative explained variances are very similar for the different centres for any number of dimensions. We wanted to take advantage of the construction of this synthethic dataset to answer the question: Do the principal components explain the same deformations? The SD dataset allows to visualize the principal component relatively to the real position of the generator vectors and the population itself. Such visualization is not possible for real dataset since shape spaces are not generated by only 2 vectors.
For this synthetic dataset, we can project principal components on the 2D space spanned by β 1 and β 2 as described in the previous paragraph. This projection allows displaying in the same 2D space subjects in their native space, and principal axes computed from the different Kernel PCAs. To visualize the first component (respectively the second one), we shot from the associated centre in the direction k * m 1 (resp. m 2 ) with k ∈ [-2 √ λ 1 ; +2 √ λ 1 ] (resp. √ λ 2 ). Results are presented in Figure  3. The deformations captured by the 2 principal axes are extremely similar for all centres. The principal axes for the 7 centres, have all the same position within the 2D shape space. So for a similar amount of explained variance, the axes describe the same deformation. Overall, for this synthetic dataset, the 6 estimated centres give very similar PCA results.
Distance Matrices We then studied the impact of different centres on the approximated distance matrices. We computed the seven approximated distance matrices corresponding to the seven centres, and the direct pairwise distance matrix computed by matching all subjects to each other. Computation of the direct distance matrix took 1,000 min (17 h) for this synthetic dataset of 50 subjects. In the following, we denote as aM(C) the approximated distance matrix computed from the centre C. To quantify the difference between these matrices, we used the following error e: e(M 1 , M 2 ) = 1 N 2 N i,j=1 M 1 (i, j) -M 2 (i, j) max(M 1 (i, j), M 2 (i, j)) (27) with M 1 and M 2 two distance matrices. Results are reported in Table  3. For visualization of the errors against the pairwise computed distance matrices, we also computed the error between the direct distance matrix, by computing pairwise deformations (23h h of computation per population), for 3 populations randomly selected.
Centring of the Centres As for the synthetic dataset, we assessed the centring of these six different centres. To that purpose, we first computed the ratio R of equation (  26 ), for the centres estimated via the centroids methods, IC1 has a R = 0. 25, for IC2 the ratio is R = 0. 33 and for PW it is R = 0. 32. For centres estimated via the variational templates initialized by those centroids, the ratio for T(IC1) is R = 0. 21, for T(IC2) is R = 0. 31 and for T(PW) is R = 0. 26. The ratios are higher than for the synthetic dataset indicating that centres are less centred. This was predictable since the population is not built from one surface via geodesic shootings as the synthetic dataset. In order to better understand these values, FIGURE 3 | Synthetic dataset SD. Illustration of the two principal components of the 6 centres projected into the 2D space of the SD dataset.
PCA on Initial Momentum Vectors As for the synthetic dataset, we performed six PCAs from the estimated centres. Figure  5  and Table  4  show the proportion of cumulative explained variance for different number of modes. We can note that for any given number of modes, all PCAs result have really similar proportions of explained variance. The highest differences in cumulative explained variance, at a given number of component, between all 6 centres is 0. 034 ( >4% of the total variance) and it is between IC2 and T(PW) for the 20th mode of variation (see Table  4 ). To reach 30% of total variance, a KPCA computed from T(C) for any C ∈ {IC1, IC2, PW} needs 3 components while C needs 4 components. To reach 50% KPCA computed from all 6 centres need 7 components. Finally to reach 90% of total variance, T(C) needs 21 components while C needs 22.
Distance Matrices As for the synthetic dataset, we then studied the impact of these different centres on the approximated distance matrices. A direct distance matrix was also computed (around 90 h of computation time). We compared the approximated distance matrices of the different centres to: (i) the approximated matrix computed with IC1; (ii) the direct distance matrix. We computed the errors e(M 1 , M 2 ) defined in equation 27. Results are presented in Table  5. Errors are small and with the same order of magnitude. Figure  6 -left shows scatter plots between the direct distance matrix and the six approximated distance matrices. Interestingly, we can note that the results are similar to those obtained by Yang  X. F. et al. (2011) , Figure  2 ). Figure  6 -right shows scatter plots between the approximated distance matrix from IC1 and the five others approximated distance matrices.
Real Dataset RD1000 Results on the real dataset RD50 and the synthetic SD showed that results were highly similar for the 6 different centres. In light of these results and because of the large size of the real dataset RD1000, we only computed IC1 for this last dataset. The computation time was about 832 min (13. 8 h) for the computation of the centroid using the algorithm IC1, and 12. 6 h for matching the centroid to the population. The ratio R of equation 26 computed from the IC1 centroid was 0. 1011, indicating that the centroid is well centred within the population. We then performed a Kernel PCA on the initial momentum vectors from this IC1 centroid to the 1,000 shapes of the population. The proportions of cumulative explained variance from this centroid are 0. 07 for the 1st mode, 0.
12 for the 2nd mode, 0. 48 for the 10th mode, 0. 71 for the 20th mode, 0. 85 for the 30th mode, 0. 93 for the 40th mode, 0. 97 for the 50th mode and 1. 0 from the 100th mode. In addition, we explored the evolution of the cumulative explained variance when considering varying numbers of subjects in the analysis. Results are displayed in Figure  7. We can first note that about 50 dimensions are sufficient to describe the variability of our population of hippocampal shapes from healthy young subjects. Moreover, for large number of subjects, this dimensionality seems to be stable. When considering increasing number of subjects in the analysis, the dimension increases and converges around 50. Finally, we computed the approximated distance matrix. Its histogram is shown in Figure  8. It can be interesting to note that, as for RD50, the average pairwise distance between the subject is around 12, which means nothing by itself, but the points cloud on Figure  6 -left and the histogram on Figure  8 , show no pairwise distances below 6, while the minimal pairwise distance for the SD50 dataset -generated by a 2D spaceis zero.
We now apply our approach to predict incomplete hippocampal inversions (IHI) from hippocampal shape parameters. Specifically, we predict the visual IHI score, which corresponds to the sum of the individual criteria as defined in  Cury et al. (2015). Only 17% of the population have an IHI score higher than 3. 5. We studied whether it is possible to predict the IHI score using statistical shape analysis on the RD1000 dataset composed of 1,000 healthy subjects (left hippocampus). The deformation parameters characterizing the shapes of the population are the eigenvectors computed from the centroid IC1, and they are the independent variables we will use to predict the IHI scores. As we saw in the previous step, 40 eigenvectors are enough to explain 93% of the total anatomical variability of the population.
We use the centred and normalized principal eigenvectors X 1,i ,. , X 40,i computed from the RD1000 database with i ∈ {1,. , 1000} to predict the IHI score Y. We simply used a multiple linear regression model  (Hastie et al. , 2009)  which is written as f (X) = β 0 + 40 i=1 X i β i where β 0 , β 1 ,. β 40 are the regression coefficients to estimate. The standard method to estimate the regression coefficients is the least squares estimation method in which the coefficients β i minimize the residual sum of squares RSS(β) = N j=1 (y jβ 0 -p i=1 x ji β i ) 2 , which leads  to the estimated β (with matrix notations) β = (X T X) -1 X T Y. For each number of dimensions p ∈ {1,. , 40} we validated the quality of the computed model with the adjusted coefficient of determination R 2 adj , which expresses the part of explained variance of the model with respect to the total variance: R 2 adj = 1 - SSE/(N -p) SST/(N -1) (28) with SSE = N i (y i -(X T 1.
DISCUSSION AND CONCLUSION In this paper, we proposed a method for template-based shape analysis using diffeomorphic centroids. This approach leads to a reasonable computation time making it applicable to large datasets. It was thoroughly evaluated on different datasets including a large population of 1,000 subjects. Codes can be found here: https://github. com/cclairec/Iterative_Centroid and synthetic data can be found here: https://doi. org/10. 5281/zenodo. 1420395. The results demonstrate that the method adequately captures the variability of the population of hippocampal shapes with a reasonable number of dimensions. In particular, Kernel PCA analysis showed that the large population of left hippocampi of young healthy subjects can be explained, for the metric we used, by a relatively small number of variables (around 50).
Moreover, when a large enough number of subjects was considered, the number of dimensions was independent of the number of subjects. The comparisons performed on the two small datasets show that the different centroids or variational templates lead to very similar results. This can be explained by the fact that in all cases the analysis is performed on the tangent space to the template, which correctly approximates the population in the shape space. Moreover, we showed that the different estimated centres are all close to the Frechet mean of the population. While all centres (centroids or variational templates) yield comparable results, they have different computation times. IC1 and PW centroids are the fastest approaches and can save between 70 and 90% of computation time over the variational template.