domain: sdv.mhep stat math

text:
Sieun Lee

Nicolas Charon

Benjamin Charlier

Karteek Popuri

Evgeniy Lebed

Marinko V Sarunic

Alain Trouvé

Mirza Faisal Beg

Atlas-based Shape Analysis and Classification of Retinal Optical Coherence Tomography Images using the Functional Shape (fshape) Framework

Keywords: Atlas-based shape analysis, Atlas generation, Computational

We propose a novel approach for quantitative shape variability analysis in retinal optical coherence tomography images using the functional shape (fshape) framework. The fshape framework uses surface geometry together with functional measures, such as retinal layer thickness defined on the layer surface, for registration across anatomical shapes. This is used to generate a population mean template of the geometry-function measures from each individual. Shape variability across multiple retinas can be measured by the geometrical deformation and functional residual between the template and each of the observations. To demonstrate the clinical relevance and application of the framework, we generated atlases of the inner layer surface and layer thickness of the Retinal Nerve Fiber Layer (RNFL) of glaucomatous and normal subjects, visualizing detailed spatial pattern of RNFL loss in glaucoma. Additionally, a regularized linear discriminant analysis classifier was used to automatically classify glaucoma, glaucoma-suspect, and control cases based on RNFL fshape metrics.

anatomy, Classification, Optical coherence tomography

Introduction

In the past two decades, optical coherence tomography (OCT) has been widely adopted in ophthalmology for noninvasive, in-vivo, micrometer-resolution imaging of the anterior and posterior segments of the eye, bringing new insights in optic diseases along with more informed clinical decision making. Since the seminal paper in 1991 [START_REF] Huang | Optical coherence tomography[END_REF], significant progresses in OCT imaging have made possible highly detailed 3D images. Fully utilizing the information in such volumetric data, especially in large population studies, requires robust, quantitative measures of variability in the images. The conventional approach is to define an anatomical parameter, measure it in each image, and compare the values. Examples of OCT image parameters include optic cup and disk measurements and thickness of retinal layers [START_REF] Medeiros | Evaluation of retinal nerve fiber layer, optic nerve head, and macular thickness measurements for glaucoma detection using optical coherence tomography[END_REF][START_REF] González-García | Reproducibility of RTVue retinal nerve fiber layer thickness and optic disc measurements and agreement with Stratus optical coherence tomography measurements[END_REF]. However, the spatial and anatomical correspondence across measurements from multiple images is often limited. This can be due to ambiguity in the parameter definition, such as in the optic cup-to-disc ratio, where the boundaries of the cup and disc often do not correspond to a single anatomical structure [START_REF] Chauhan | From clinical examination of the optic disc to clinical assessment of the optic nerve head: A paradigm change[END_REF][START_REF] Young | Comparison of the clinical disc margin seen in stereo disc photographs with neural canal opening seen in optical coherence tomography images[END_REF].

Poor intra-and inter-subject correspondence can be also due to a lack of clear spatial references, as in the case of retinal layer thickness maps.

In the latter case, comparing or computing a group average from multiple OCT images requires artificially defined regional correspondence. This often relies on gross measures of orientation and distance from an anatomical landmark. In the common sectoral retinal layer thickness analysis of OCT scans, sectors are delimited by superior, inferior, temporal, and nasal orientations, and distance from the foveal pit or optic disc. Such an approach uses limited information, and leaves questions such as how comparable measurements from two eye would be if one eye's optic disc was larger or more elliptical than the other.

Averaging over a region may mitigate some of the variability due to unknown mismatch, but comes at the cost of spatial sensitivity of the measurement.

One way to address the issue is to define parameters that are more anatomically consistent. Several groups [START_REF] Strouthidis | Comparison of clinical and spectral domain optical coherence tomography optic disc margin anatomy[END_REF][START_REF] Chauhan | From clinical examination of the optic disc to clinical assessment of the optic nerve head: A paradigm change[END_REF][START_REF] Young | Comparison of the clinical disc margin seen in stereo disc photographs with neural canal opening seen in optical coherence tomography images[END_REF], suggested Bruch's Membrane Opening (BMO) based parameters as more reliable structural measurements than the cup-to-disk ratio. Peripapillary sector shapes were adjusted to better align to individual BMO dimensions [START_REF] Lee | Optic nerve head and peripapillary morphometrics in myopic glaucoma[END_REF], and foveal-BMO axis has been shown to be a better reference for angular orientation than image frame in its anatomical and physiological justification and measurement repeatability [START_REF] Chauhan | From clinical examination of the optic disc to clinical assessment of the optic nerve head: A paradigm change[END_REF][START_REF] He | Anatomic vs. acquired image frame discordance in spectral domain optical coherence tomography minimum rim measurements[END_REF].

Registration of OCT images fundamentally addresses the correspondence issue; however, most of its application has been in averaging multiple images for noise reduction and motion correction [START_REF] Jørgensen | Enhancing the signal-to-noise ratio in ophthalmic optical coherence tomography by image registration-method and clinical examples[END_REF][START_REF] Young | Real-time high-speed volumetric imaging using compressive sampling optical coherence tomography[END_REF], or rigid alignment of time-course images [START_REF] Niemeijer | Registration of 3D spectral OCT volumes using 3D SIFT feature point matching[END_REF]. Three previous studies have specifically examined nonrigid registration of OCT data from multiple subjects in the context of shape variability analysis. In Gibson et al. [START_REF] Gibson | Optic nerve head registration via hemispherical surface and volume registration[END_REF], segmented optic cups were registered to a single template optic cup, first by rigid and nonrigid intensity-based volumetric registration followed by spherical mapping and demons algorithm [START_REF] Yeo | Spherical demons: Fast diffeomorphic landmark-free surface registration[END_REF].

In Chen et al. [START_REF] Chen | Analysis of macular OCT images using deformable registration[END_REF], registration of macular OCT scans used rigid alignment of foveae, affine registration of A-scans to match the inner and outer retinal boundaries, and smooth deformation of A-scans using radial basis functions for refined alignment of the retinal layers. [START_REF] Lee | Exact surface registration of retinal surfaces from 3-d optical coherence tomography images[END_REF] extended upon [START_REF] Gibson | Optic nerve head registration via hemispherical surface and volume registration[END_REF] and represented segmented retinal surfaces utilizing the framework of mathematical currents. Two surfaces were brought into proximity by minimizing a functional of reproducing kernel Hilbert space (RKHS) norm-based energy and a dissimilarity term. This was followed by spherical demons registration to establish point-to-point correspondence. In above works, target subjects were registered to one subject as the template, which leaves the problem of the bias in template choice.

In this report, we propose to exploit a novel approach to analyzing variability in retinal OCT images using the functional shape or fshape framework. Fshapes and some of its abstract mathematical properties were introduced in [START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF]. The fshape framework considers a geometrical surface, such as a retinal layer surface, and functions defined on the surface, such as retinal layer thickness, together as a single object. The "distance", i.e. the shape difference between two retinas is measured by both geometry and one or more physiological or morphological signals defined on the geometrical surface. This is a conceptual departure from previous works in which the goal is to first establish anatomical correspondence between two eyes and then to compare the functional values at resultant corresponding locations, and the first work on atlas-based morphological analysis of the retina. There are several advantages to the fshape approach. First, the mathematical abstraction allows variability analysis, or comparison, beyond the frame of the anatomy; given any one or more surfaces and functional values, the fshape framework can measure the inter-subject or inter-time points difference. This allows for combination of any number of features -for example, Inner Limiting Membrane (ILM) geometry and total retinal thickness, Bruch's membrane (BM) geometry with RNFL thickness, or ILM and BM surface geometry -and to investigate individually and jointly which features are more or less significant in differentiating between a disease and control group. The core of the fshape framework is generation of a mean template of multiple fshapes. An initial template or hypertemplate is taken as a simple model of a prototype fshape. This template is evolved through an optimization process that simultaneously minimizes i) geometric-functional distance from the current template to the observations, and ii) dissimilarity between the transformed mean template and the observations. This approach eliminates the need to choose one of the existing data as a template, mitigating the issue of template selection and bias. In addition, the fshape framework registration does not rely on specific anatomical or image features based on prior knowledge. Given decent quality segmentation and measurements, the generality and versatility of the algorithm allows it to be applied broadly and robustly.

The fshape framework combines a distance of fshape transformations within the fshape bundle of a template and dissimilarity measures between arbitrary fshapes based on an extension of varifold spaces and norms. By optimization of these metrics, a mean atlas can be created, and serve as the reference to all fshapes from the dataset. In addition, the algorithm simultaneously provides the variations in shape and signal of this template to each observation. These outputs then constitute the basis of inter-subject analysis and group classification through a wide choice of possible statistical analysis tools.

This paper proposes a very first set of applications of the fshape methodology to the OCT dataset while extending the scope of the framework to statistical learning and classification. Namely, after a short application-oriented presentation of fshapes, we propose a formal Bayesian derivation of the template estimation problem studied in [START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF] and augment the current approach with a statistical analysis module for the estimated geometric and functional features. To demonstrate the clinical relevance and application of the algorithm, we generated a mean atlas and performed automated classification of glaucoma eyes and healthy control eyes based on both geometry and thickness of posterior RNFL surface and RNFL thickness. We also present group averages and t-test of the RNFL thickness over the mean surface between the healthy and glaucomatous groups to emphasize how the fshape approach can reveal morphological and functional patterns in cross-sectional or longitudinal data.

Methods

Image acquisition and processing

The OCT images in this study were acquired at the Eye Care Centre at Vancouver General Hospital in Vancouver, British Columbia, Canada, using a custom prototype swept-source OCT machine developed at Simon Fraser University (SFU) with a 1060-nm wavelength light source. Three-dimensional volumetric images centered at optic nerve head were acquired over a 5-8 mm 2 region with 2.8 mm in depth. The image resolution was approximately 6 µm in the axial direction and 12-21 µm in the lateral direction, depending on the subject eye's axial length. Each volume consisted of 1024 x 400 x 400 voxels, with 1024 voxels in axial direction. The images were corrected for axial motion using cross-correlation of adjacent frames. A bounded variation regularization method was used to reduce the effect of speckles and enhance retinal boundaries for segmentation and visualization. Retinal layers were segmented automatically in 3D using a graph-cut based surface segmentation algorithm implemented in MATLAB [START_REF] Li | Optimal surface segmentation in volumetric images -a graph-theoretic approach[END_REF][START_REF] Garvin | Intraretinal layer segmentation of macular optical coherence tomography images using optimal 3-D graph search[END_REF][START_REF] Lee | Comparative analysis of repeatability of manual and automated choroidal thickness measurements in nonneovascular age-related macular degeneration[END_REF]. In this study four boundaries were segmented for two layers: RNFL and choroid. Left (OS) eyes were flipped so that all eyes were in the right (OD) eye orientation. Layer thickness was measured at each point of the layer's posterior surface as the closest distance to the anterior surface. In addition to the automated segmentation, BMO in each eye was segmented manually, and a best-fit ellipse was generated by Principal Component Analysis (PCA). The BMO ellipse was used to crop the surfaces near the BMO where the layers terminate. Fig. 1 shows an example optic nerve head (ONH) scan, segmentation, and the resulting posterior RNFL surface with RNFL thickness mapped on the surface.

Atlas estimation

The central goal of the fshape framework is to recover the inter-subject variability both in the geometry of the retina surfaces as well as in their signals (i.e the thickness maps). Following the standard process in computational anatomy [START_REF] Joshi | Unbiased diffeomorphic atlas construction for computational anatomy[END_REF][START_REF] Ma | A Bayesian generative model for surface template estimation[END_REF][START_REF] Zhang | Bayesian Estimation of Regularization and Atlas Building in Diffeomorphic Image Registration[END_REF], the primary step is to estimate an atlas from the population, which serves as a template object with the variation of the subjects in the population. The obvious difficulty in this situation is that both geometrical supports and functional maps on the supports vary concurrently. The notion of functional shape or fshape [START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF] has the crucial advantage of treating both geometry and function together while providing important flexibility for the atlas generation and registration algorithms.

The following sections give a compact and application-oriented exposition of the theoretical presentation presented in [START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF] and emphasizes the Bayesian interpretation of atlas estimation in the fshapes context.

Background: fshapes and fshape spaces

In the fshape framework, a geometrical structure and its associated scalar field are considered as a single object and processed jointly. Thus, in general, an fshape consists of a pair (X, f ) where X is the geometrical support, i.e a surface in the 3D ambient space, and f is a function defined on this surface; in our particular application, the pair of a retinal layer surface and a retinal layer thickness mapped on the surface.

Transformation of an fshape (X, f ) can be modeled by combining a deformation φ of the geometrical support and a signal change ζ. In the simplest setting considered in [START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF], φ is a diffeomorphism of R 3 and ζ a residual function on X. The combination (φ, ζ) then 'acts' on the fshape (X, f ) as

(φ, ζ) • (X, f ) = (φ(X), (f + ζ) • φ -1 ) (1)
meaning that the surface is transported by deformation φ while signal is modified by adding the residual ζ and then mapped onto the deformed surface φ(X).

Quantifying such transformations is done first by introducing a model of deformation group, our reference model in this paper is the widely studied Large Diffeomorphic Metric Mapping (LDDMM) framework of [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF], in which diffeomorphisms are constructed as the flow of time-dependent velocity fields v ∈ L 2 ([0, 1], V ) with V a reproducing kernel Hilbert space (RKHS) of smooth velocity fields of R 3 . On the other hand, we shall consider L 2 signals on surfaces, i.e f, ζ ∈ L 2 (X) with the corresponding surface 2-norm on the shape X. With φ = φ v 1 the flow of v up to time t = 1, the energy of the fshape transformation (φ, ζ) of (X, f ) considered in [START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF] is given by:

E X (v, ζ) 2 = 1 2γ 2 V 1 0 v t 2 V dt + 1 2γ 2 f ζ 2 L 2 (X) (2)
where σ V , σ f are weighting parameters between the geometric and functional energies.

The energy in (2) gives a way to measure a notion of distance between two fshapes but only in the situation where these can be exactly mapped on each other in the transformation model introduced above. In practice, this is unrealistic or not desirable since datasets present inter-subject variability caused by noisy irregularity in the shape and signals, variations that are not well represented by the smooth geometric-functional transformations of the previous section.

In registration problems, it is thus common to introduce additional dissimilarity (or data fidelity) terms to the deformation cost, which may be interpreted equivalently as a noise model on the observations as we shall see in the next section. In the situation of functional shapes, unlike usual images, the absence of any explicit correspondence between their geometrical supports does not enable direct comparison of their signals.

A possibility to overcome this issue was examined thoroughly in (Charon and Trouvé, 2013;[START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF] and consisted in extending previous works on curves and surfaces [START_REF] Glaunès | Diffeomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching[END_REF]Charon and Trouvé, 2013). In the fvarifold framework that we very briefly sum up, an fshape (X, f ) is represented as a distribution on the product space R 3 × P (R 3 ) × R, with P (R 3 ) being the projective space of all lines in R 3 . This distribution, written µ (X,f ) is formally a sum of Diracs involving the position of the shape points x ∈ X, the attached (unoriented) direction of the normal vector ← → n (x) ∈ P (R 3 ) and signal function f (x).

These distributions are compared via metrics derived from tensor product of positive kernels of the form k g ⊗k n ⊗k f on R 3 ×P (R 3 )×R. The induced Hilbert (pseudo-)metric between any two fshapes (X, f ) and (Y, g) writes explicitly:

µ (X,f ) ,µ (Y,g) W * = X×Y k g (x, y)k n ( ← → n X (x), ← → n Y (y))k f (f (x), g(y))dσ(x)dσ(y) (3)
where dσ denotes the surface area measures of X and Y respectively. Then

A . = µ (X,f ) -µ (Y,g) 2 W * = µ (X,f ) -µ (Y,g) , µ (X,f ) -µ (Y,g) W *
can be taken as a suitable dissimilarity term between (X, f ) and (Y, g) since it enforces proximity in the respective geometries and signal functions concurrently.

Heuristics of atlas estimation

Let us now consider a population of N observations (X i , f i ) i=1,...,N where the X i 's are the retinal surfaces and f i the corresponding thickness maps. We introduce a forward generative model in the line of [START_REF] Durrleman | A forward model to build unbiased atlases from curves and surfaces[END_REF] for which observations are noisy geometric-functional transformations of a common unknown template fshape (X, f ) plus additional noise terms:

(X i , f i ) = (φ i , ζ i ) • (X, f ) + ε i , for all i = 1, . . . , N . (4)
Above, φ i is the flow of a vector field

v i ∈ L 2 ([0, 1], V ) and (v i , ζ i ) are regarded
as hidden latent variables of the transformations from template to subjects, ε i 's are noise variables.

Considering i.i.d. variables ε i , we may define a noise model on fshapes based on fvarifold metrics:

p(ε i ) = p (X i , f i )|(X, f, v i , ζ i ) ∝ e - µ (φ i ,ζ i )•(X,f ) -µ (X i ,f i ) 2 W * 2γ 2 W
which is a Gaussian model with respect to the metric . W * . Note that this is only formal for the infinite dimensional space of fvarifolds but can be given a rigorous sense if restricted to a predefined discrete grid, similarly to [START_REF] Gori | Bayesian Atlas Estimation for the Variability Analysis of Shape Complexes[END_REF].

As for the latent variables (v i , ζ i ), we take the following prior deriving from the energy (2):

p((v i , ζ i )|(X, f )) ∝ e -E X (v i ,ζ i ) 2
which is essentially assuming independent (formal) Gaussian distribution on v and ζ in their respective metric spaces.

Finally, we also model the template (X, f ) as a random variable itself. Inspired from the hypertemplate model for shape atlases of [START_REF] Ma | A Bayesian generative model for surface template estimation[END_REF], we represent (X, f ) as a transformation of a given hypertemplate fshape (X

0 , f 0 ), i.e (X, f ) = (φ 0 , ζ 0 ) • (X 0 , ζ 0 ) for a deformation φ 0 and a residual ζ 0 ∈ L 2 (X 0 ).
As previously, the prior on the template is:

p((v 0 , ζ 0 )) ∝ e -E X 0 (v 0 ,ζ 0 ) 2 .
With an hypertemplate (X 0 , f 0 ) fixed by the user, estimating the template then amounts to computing the maximum a posteriori (MAP) estimate of (v 0 , ζ 0 ) knowing the observations (X i , f i ). With Bayes rules this leads to minimizing:

- N i=1 log p(X i , f i |v 0 , ζ 0 , v i , ζ i )p(v i , ζ i ) -log p(v 0 , ζ 0 ) .
The first term involves the integral with respect to the probability distribution of the latent variables (v i , ζ i ). As there is no closed form expression of this integral, we use the standard Fast Approximation with Modes and replace it by

max v i ,ζ i p(X i , f i |v 0 , ζ 0 , v i , ζ i )p(v i , ζ i
) leading eventually to the following variational problem:

v 0 * , ζ 0 * , (φ i * ) i , (ζ i * ) i ) = arginf v 0 ,ζ 0 ,(v i ,ζ i )i 1 2γ 2 V 1 0 v 0 t 2 V dt + 1 2γ 2 f ζ 0 2 L 2 (X 0 ) + N i=1 1 2γ 2 V 1 0 v i t 2 V dt + 1 2γ 2 f ζ i 2 L 2 (X 0 ) + 1 2γ 2 W µ (φ i ,ζ i )•(X,f ) -µ (X i ,f i ) 2 W * (5)
where we remind that for i = 0, . . . , N deformations φ i are the flows of the v i and X = φ 0 (X 0 ).

200

The previous paragraphs underline the Bayesian interpretation behind variational problem (5), for which the existence of solutions under certain conditions was addressed in [START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF]. Note that the variances γ 2 V , γ 2 f , γ 2 W act as weighting coefficients between the different terms. In the rest of the paper, we will consider these as parameters of the algorithm. However, we wish to 205 point out that with a proper discretization of the previous probability densities, this could be improved by also estimating the optimal variances, extending approaches developed for images in [START_REF] Zhang | Bayesian Estimation of Regularization and Atlas Building in Diffeomorphic Image Registration[END_REF] and surfaces in Ma et al.

(2010); [START_REF] Gori | Bayesian Atlas Estimation for the Variability Analysis of Shape Complexes[END_REF].

Discrete fshapes 210

In the discrete setting, a functional surface is a textured mesh given as a P × 3 matrix x of coordinates of the P vertices x k in R 3 , a P × 1 column vector

f of the P values f k . = f (x k ) ∈ R
of the signal at the vertices and a T × 3 connectivity matrix C (a triangulation for simplicity). The L 2 norm of f on the surface is here approximated with P0 finite elements

f 2 L 2 (X) ≈ T =1 f 2 |T | (6)
where for the th triangle, |T | is the area of the triangle and f is the average of the three signal values at its vertices.

Optimal deformation fields minimize the metric

1 0 v t 2
V dt for given final time conditions and are therefore geodesics in the context of diffeomorphism groups. It has been shown [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Arguillere | Shape deformation analysis from the optimal control viewpoint[END_REF] that such geodesic flows are governed by Hamiltonian equations. In the present case of discrete set of particles, geodesics are parametrized by initial momenta p = (p k ) ∈ R P ×3 attached to every vertex x k and the following shooting equations:

   ẋk (t) = v t (x k (t)) = P l=1 K V (x l (t), x k (t))p l (t) ṗk (t) = - P l=1 p k (t) • p l (t)∂ 1 K V (x l (t), x k (t)) (7)
given the initial conditions x(0) = x and p(0) = p and where K V is the 3 × 3 matrix-valued kernel function associated to the RKHS V .

We thus parametrize geometric-functional transformations of the discrete fshape (x, f , C) by the momenta vectors p and a residual discrete signal ζ :

(φ, ζ) • (x, f ) = (x k (1)) 1≤k≤P , (f k + ζ k ) 1≤k≤P . (8)
As for the fvarifold-norm fidelity term of (3), we use a similar discrete approximation where each face triangle T is approximated by a single Dirac with position x the vertices' barycenter, f the mean signal value, and ← → n the unit unoriented normal vector to the triangle. With two discrete fshapes (x, f , C 1 ) and (y, g, C 2 ),

µ (x,f ) , µ (y,g) W * = T1 k=1 T2 =1 |T x k ||T y |k g (x k , y )k n ( ← → n k x , ← → n y )k f (f k , g ) . (9)

Numerical scheme for atlas estimation

With the above discrete model and a given discrete hypertemplate (x 0 , f 0 , C 0 ) the atlas estimation problem of (5) reduces to the minimization of the function:

J (p 0 , ζ 0 , (p i ) i , (ζ i ) i ) . = 1 2γ 2 V (p 0 ) T K V (x 0 , x 0 )p 0 + 1 2γ 2 f ζ 0 2 L 2 + N i=1 1 2γ 2 V (p i ) * K V (x, x)p i + 1 2γ 2 f ζ i 2 L 2 + 1 2γ 2 W µ (x i , f i ) -µ (x i ,f i ) 2 W * (10)
where • L 2 is given by ( 6) and the deformed fshapes (x, f ) .

= (φ 0 , ζ 0 )•(x 0 , f 0 ), (x i , f i ) . = (φ i , ζ i ) • (x, f
) are obtained by the shooting equations ( 7). We will

write (p 0 * , ζ 0 * , (p i * ) i , (ζ i * ) i
) for a minimizer of J . The overall principle is illustrated by Figure 2. This is now a finite yet very high-dimensional optimization problem since each variable is roughly of dimension P , the number of hypertemplate vertices (which is typically of the order of P ≈ 5000 in the simulations of this paper) and it is also non-convex. We thus implement a Polak-Ribiere conjugate gradient descent where all variables are simultaneously updated at each step. The algorithm is also coupled with basic line search for each variable to increase the speed of convergence. The gradients with respect to signal variables ζ 0 and (ζ i ) are easily computable from the expression in ( 6) and ( 9). Gradients with respect to initial momenta variables p 0 , (p i ) are computed based on forwardbackward shooting procedure as in usual geometric registration, using Euler

f 3 * ) = (φ 3 * , ζ 3 * ) • (x * , f * ) ∈ F
midpoint scheme for the time-discretization of Hamiltonian and adjoint Hamiltonian systems. We refer to [START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF] for the detailed expressions.

The computational burden of the algorithm concentrates in repeated computations of sums of kernels and kernel derivatives in the shooting equations and evaluation of the functional, these specific parts are implemented in CUDA to achieve competitive speed and precision.

Apart from the construction of a discrete fshape hypertemplate, the rest of the estimation is fully automatic but depends on a few parameters including the weighting coefficients γ V , γ f , γ W between the different energy terms, but also the parameters associated to the kernels that define the metrics on velocity fields and on fvarifolds. In the applications of this paper, we use isotropic Gaussian kernels for k g and k f in (3), which are determined by two variance parameters σ g and σ f . Heuristically, σ g and σ f measure the sensitivity scale of the fidelity term respectively in the spatial and signal domain. While too big values would typically result in a very imprecise overlap of fshapes, too small values on the other hand could either induce over-fitting or result in no evolution at all in the case where fshapes are initially too far apart relative to these scales. We address this issue by completing the approach of [START_REF] Charlier | The Fshape Framework for the Variability Analysis of Functional Shapes[END_REF] with a coarseto-fine multiscale scheme that allows for successively decreasing σ g and σ f in the minimization. A last important parameter is the scale of the vector kernel defining the space of velocity fields V . Instead of using a single Gaussian kernel, we also allow for multiscale deformations by incorporating sums of Gaussians with different scales, as proposed in [START_REF] Bruveris | Mixture of Kernels and Iterated Semidirect Product of Diffeomorphisms Groups[END_REF].

Variability Analysis and Classification

Augmenting the current fshape framework with a classification module is a natural extension with several merits. On experimental data with confirmed diagnosis, automated classification can demonstrate the mean template-based variability of retinas has indeed anatomical and clinical relevance. Furthermore, a classification module can determine whether a particular anatomical feature, thickness or topology of a layer, is significantly correlated with a disease. This can lead to better understanding of the structural manifestation of the disease in terms of cause or effect and directly assist in clinical decision making.

The previous section demonstrated how, given a set of observations, in our case retinal surfaces, the fshape approach can generate an estimated mean atlas of the dataset. Such a group average retina in itself is useful for purposes such as qualitative group comparisons. To quantify more accurately group differences and perform statistical classification, it is necessary to rely on the latent variables (φ i , ζ i ) that characterize deviations in shape and thickness from the template for each subject. These are automatically estimated in the previous framework and, as we saw, are optimal or 'geodesic' transformation paths for the energy of (2) which can serve as metrics for inter-eye variability. Statistical analysis on shape spaces typically exploits the linearization on diffeomorphism groups given by the initial velocity fields (or momenta) and extract main shape variance components via dimensionality-reduction methods like PCA [START_REF] Vaillant | Statistics of diffeomorphisms via tangent space representation[END_REF]. We easily extend that principle to fshapes by ap-pying PCA to the initial momenta (p i * ) i and functional residuals (ζ i * ) i of Section 2.2.4. However, for classification tasks, PCA directions need not be discriminative with regard to the differences across populations (e.g controls vs glaucoma in our case). Supervised methods such as a Linear Discriminant Analysis (LDA)

classifier [START_REF] Hastie | The Elements of Statistical Learning[END_REF] are generally more adequate for pathology detection [START_REF] Durrleman | Deformetrics : morphometry of shape complexes with space deformations[END_REF]. The method we propose in the present work is a variation of LDA tuned to the specific structure of variables issued by fshape atlases, which deals with a relatively low number of high-dimensional samples in particular Hilbert spaces. The general setting is described below and the results of the classification experiments using this method are presented in section 3.3.

Regularized LDA on Hilbert space valued dataset

LDA is a classic linear classification method that may be viewed as a weighted principal component analysis (see chapter 4 of [START_REF] Hastie | The Elements of Statistical Learning[END_REF] for an introduction). We discuss hereafter the general case in which data points

x 1 , • • • , x N
belong to a separable Hilbert space H of a possibly infinite dimension. In the context of fshapes and the applications of this paper (cf section 3.3), the x i 's are either functional residuals for which H = L 2 (X) the space of L 2 functions on the template (or the discrete equivalent with the metric given by ( 6)), and initial momenta in which case H = V * .

Within-and between-class scatter operators

We assume that the N observations at hand are divided into K classes

C 1 , • • • , C k such that {1, • • • , N } = K k=1 C k . The mean of the k-th class is xk = 1 N k i∈C k x i where N k = Card(C k ). Let us define the within-class scatter operator S w : H → H S w (•) = 1 N K k=1 i∈C k x i -xk , • H (x i -xk )
and the between-class scatter operator

S b : H → H S b (•) = 1 N K k=1 N k xk -x, • H (x k -x) .
In the case where H is of finite dimension, the matrix associated to the operator S w (or S b ) is the standard within-class (or between-class) scatter matrix.

Finite dimensional representation

We denote H 0 = Span (x i , . . . , x N ) ⊂ H and q = dim(H 0 ). In our framework, H is high (or infinite) dimensional, and, in practice, q equals to N .

We then denote f H0 : H 0 → H 0 the restriction to H 0 of any linear mapping

f : H → H such that f (H 0 ) ⊂ H 0 .
In particular, we may consider the restrictions S w,H0 and S b,H0 of the scatter operators S w and S b , respectively.

Let us now choose an arbitrary isometric linear mapping L : H 0 → R q in order to get a finite representation of the data

xi = Lx i ∈ R q .
By definition we have xi R q = x i H . From now on, we will work with this new representation of the data and we may consider their corresponding scatter operators Sw , Sb : R q → R q defined as, Sw = LS w,H0 L † and Sb = LS b,H0 L † , where L † : R q → H 0 is given by L † a, x H = a, Lx R q for any x ∈ H 0 and a ∈ R q .

When the x i 's form a basis of H 0 , an effective way to build an isometric mapping L is to introduce the mapping γ(x) = ( x, x i H ) N i=1 for any x ∈ H 0 and the Gram matrix

G = [ x i , x j H ] N i,j=1 ∈ R N ×N .
It is then easy to check that the linear mapping Lx = G -1 /2 γ(x) is isometric.

Discriminant axes in the finite dimensional space

In our framework, the within-class scatter operator Sw may not be invertible in general. We consider the following regularization of the within-class scatter operator

Sε w = Sw + εId R q , ( 11
)
where Id R q is the identity matrix and ε > 0 is a regularization parameter which has to be calibrated by the user as described below.

The discriminant spaces of the LDA are given by the eigendecomposition of Ãε = ( Sε w ) -1 Sb .

The number of non-vanishing eigenvalues is limited by the rank of Sb which is less than K -1. In general, we have the K -1 corresponding unit eigenvectors ũ1 , . . . , ũK-1 ∈ R q that will be used to derive the discriminant axes in H.

Classification with regularized LDA

The rationale behind the isometric dimension reduction is that the unit vectors

u = L † ũ ∈ H 0 , = 1, . . . , K -1
are the eigenvectors corresponding to the non-vanishing eigenvalues of A ε H0 = (S w,H0 + εId H0 ) -1 S b,H0 , as in [START_REF] Friedman | Regularized discriminant analysis[END_REF]. This classic trick avoids numerical issues as the matrix inversion is performed in the small dimensional space R q .

Then, given a new observation y ∈ H we use directly the discriminant axes u 1 , . . . , u K-1 to define a classification rule. For instance, when K = 2 (a twoclass classifier) we have u 1 ∝ (S ε w ) -1 (x 1 -x2 ) where the ∝ symbol means "collinear to". The discriminant rule is then a threshold on y → u 1 , y -x H .

Calibration of the regularization parameter

The regularization parameter ε > 0 can be optimized with a leave-p-out cross-validation (CV) procedure. Note that we do not need to compute the Gram matrix G appearing in the definition of the isometric mapping L at each stage of the CV. This is very helpful since the operation may be costly depending on the magnitude of N . Instead, the N × N "full" Gram matrix G f with all the observations is computed once, and the (N -p) × (N -p) Gram matrices G are obtained by deleting the p rows and columns of G f corresponding to the p left-out observations.

Experimental Results

Peripapillary OCT images from 53 eyes with confirmed diagnosis from 10 controls, 10 bilateral glaucoma patients, and 7 unilateral glaucoma patients were included in the experiment. Written consent forms were obtained from all participants and ethics review was approved by the Office of Research Ethics at Simon Fraser University (SFU) and the Research Ethics Board of the University of British Columbia (UBC). In addition to OCT imaging, all participants were subject to a battery of standard tests, including dilated stereoscopic examination of the optic nerve, stereo disc photography analysis and visual field abnormality check, to ensure there was no other pathology present. The acquired images were smoothed, segmented, and measured for retinal layer thickness as described in Section 2.1.

Mean template generation

The mean templates of retinal nerve fiber layer (RNFL) posterior surfaces and associated RNFL thickness maps were generated with the atlas estimation algorithm described above. The surfaces were rigidly registered prior to the mean template generation, as shown in in Fig. 3 (a). The physical dimensions of the images vary as the imaging field of view changes depending on the axial length of the eye. We also used a multiscale approach for the different parameters: the deformation kernel is a mixture of Gaussian of scales 2.4, 1.2, 0.6 and 0.3 mm while the algorithm is run successively with three sets of scale parameters σ g = 0.8, 0.4, 0.2 and σ f = 0.3, 0.2, 0.1mm for the fidelity term as described in Section 2.2.4. Fig. 3 (e), (f), and (g) show typical sectoral averaging of the RNFL thickness in all eyes, bilaterally normal eyes, and glaucomatous eyes groups. The sectorization process is described in detail in [START_REF] Lee | Optic nerve head and peripapillary morphometrics in myopic glaucoma[END_REF]. Briefly, sectors are drawn in each image by the distance from the BMO (boundaries here are at 0.25, 0.5, 1.0, 1.25 mm from the BMO) and angular sections of superior, inferior, nasal, temporal, and in-between regions determined by the relative position from the reference line extending horizontally from the BMO centroid to the right frame border of the image. Unlike the fshape metrics, sectoral average thickness is measured in the sectors that are defined individually in each eye, and a group average is given by averaging the values in the corresponding sec- Fig. 4 shows the observed RNFL surfaces (top row) and their approximations (bottom row) from the common mean template of 53 RNFLS in Fig. 3 (b). As described by equation ( 8), the i th approximation is a deformed version of the template, such that xi

* . = φ i * (x * ), f i * . = (f * + ζ i * ) • (φ i * ) -1 .
The salient shape features are reproduced, in particular the RNFL thickness pattern and the BMO location and size. The level of detail in approximation can be tuned by user specified parameters including the kernel sizes for geometric The log of the p-value presented on the entire surface.

deformation and dissimilarity metric, and runtime parameters in the gradient descent optimization. Although sharper reconstructions can be achieved by adjusting the parameter γ W , we note that this will not necessarily induce a better discriminative power due to the risk of over-fitting.

T-test and z-score map between healthy and glaucomatous RNFL thickness residuals

The mean templates of the normal (c) and glaucomatous (d) RNFLs in Fig. 3 qualitatively show the difference between the two groups. In order to identify the spatial locations of the thickness difference, a vertex-wise t-test was performed using the functional (thickness) residual ζ i , the vertex-wise thickness offset between the template and its approximation of the i th observation. The t-test compared the residuals of 26 normal RNFLs and 27 glaucoma RNFLs at each vertex on the mean template of all RNFLs in Fig. 3 (b).

The result is shown in Fig. 5. On the left, the voxels with p-values less than 0.05 are marked in red and indicate the regions with statistically significant RNFL thickness difference between the normal and glaucomatous RNFLs.

Most of the peripapillary region shows significance, except at the image boundaries, where the fshape correspondence may be less reliable due to the different imaging field of view size, and in the region immediately temporal to BMO. with previous studies [START_REF] Morrison | Glaucoma: Science and Practice[END_REF][START_REF] Kanamori | Evaluation of the glaucomatous damage on retinal nerve fiber layer thickness measured by optical coherence tomography[END_REF] that the inferior peripapillary region is the most distinguishing between normal and glaucomatous eyes, especially in the early stage of the disease. Unlike the previous studies, however, our analysis shows full spatial detail without averaging in sectors, and reveals a clear pattern of statistical significance that resembles the characteristic thickness pattern of a healthy RNFL. The result suggests that the most significant amount of glaucomatous thinning, or the earliest of the thinning, may occur along the RNFL ridges where the RNFL is naturally the thickest. The map shows distinctions between the healthy and early glaucoma groups, and between the early glaucoma and moderate to sever glaucoma groups, along with the individual variability in the relationship between the degrees of RNFL loss and vision loss among the glaucomatous eyes. It is noteworthy that even within the group diagnosed as healthy, the eyes with lower VFMD show generally low ζ values.

Classification

With the methods detailed in Section 2.3, classification experiments were performed using the LDA classifier on the momenta p i , the vertex-wise geometric deformation of the template to the i th observation, and functional residuals ζ i , the vertex-wise thickness offset between the template and the approximation of the i th retina. The classification experiments are presented as a demonstration of the discriminative power of the fshape metrics and how it captures the anatomical variability in RNFLs.

For each experiment, the LDA classifier was trained by leave-one-out crossvalidation, with the best-performing regularization parameter ε (see section 2.4.5) selected by a nested leave-one-out procedure from values between 0.001 to 1.

Healthy vs. Glaucoma

Age-matched (59.6 ± 6.7) normal (N = 10) and glaucomatous (N = 18)

RNFLs were classified with the result in Table 1. The result indicates that the fshape metrics of the peripapillary RNFL posterior surface and RNFL thickness can predict the clinical diagnosis of glaucoma with high accuracy, and confirms the connection between vision loss that bases glaucoma diagnosis, and characteristic morphological changes in the RNFL.

Healthy vs. Suspect

Glaucoma is generally bilateral, affecting both eyes of the patient, and often asymmetric, such that the affected fellow eyes exhibit different degrees of severity. In unilateral primary glaucoma in which only one eye is diagnosed with glaucoma, the healthy fellow eye is at a greater risk of developing glaucoma in the future than healthy eyes of bilaterally normal subjects [START_REF] Kass | Prognostic factors in glaucomatous visual field loss[END_REF][START_REF] Susanna | The visual prognosis of the fellow eye in uniocular chronic open-angle glaucoma[END_REF]. Based on this, we labeled 7 healthy fellow eyes from unilateral glaucoma cases as suspect, and attempted to detect these from healthy eyes of bilaterally nonglaucomatous subjects. In LDA classification, an intuitive way to understand the result is to visualize the classifier, or the direction that yields the maximum between-class variance.

The classifier for the Healthy vs. Glaucoma dataset is visualized in Fig. 7 on the mean template RNFL, in functional residual (RNFL thickness), and x, y, and z coordinates of initial momenta representing the template deformation.

The relative magnitude can be interpreted as the degree of contribution to the classification, and it can be seen that RNFL thickness was a more decisive factor in the classification of healthy and glaucomatous eyes than posterior RNFL surface geometry. The map in Fig. 7 A, which is the spatial pattern of classification contribution of RNFL thickness residual, is similar to the left image of Fig. 5, which is the statistical significance map of the group difference between the two classes. This confirms that the classification was the most influenced by the regions where there are the most significant difference between the healthy and glaucomatous eyes.

Classification with RNFL thickness sectoral averages

The above classifications for Healthy vs. Glaucoma and Healthy vs. Suspect were repeated with the same eyes with the RNFL thickness sectoral averages as described in Section 3.1. In this case, correspondence is assumed between the corresponding sectors across the subjects. Comparing the results in Table 2 with that in Table 1 of the fshape metrics, the classifier performance is significantly reduced. The higher classification success rate between the healthy and glaucomatous groups than between the healthy and suspect groups is consistent with the fshape metrics result in 1 and can be attributed to the large RNFL thickness difference between the two groups, visualized in the sectoral group averages in Fig. 3. The thinning in the glaucomatous RNFLs is distinguishable even after sectoral averaging. However, any difference between healthy and suspect RNFLs is likely more subtle, as the eyes in the suspect groups are nonglaucomatous with normal visual function and without clinically observable structural degradation. Such fine distinctions may be smoothed out by the sectoral averaging. The result of this experiment, along with the results presented above, suggests that spatially more detailed comparison and analysis are possible with the fshape metrics than the conventional sectorization approach.

Finer sectorization for more localized comparison will be limited by decreasing confidence in the correspondence between the same sectors from different eyes.

Comparison with standard surface registration

To illustrate the influence of functions on the estimated deformation and residuals, we performed a simple matching experiment between two surfaces, by a standard surface LDDMM [START_REF] Lee | Exact surface registration of retinal surfaces from 3-d optical coherence tomography images[END_REF] and by the fshape approach.

In the former, the residuals were retrieved by closest point projection from the target to the mapped surface. As shown in second column of Figure 8 ssLDDMM benefits from the pre-alignment that occurs in the imaging system by using the chin and forehead restraint providing a vertical orientation of the head for retinal imaging. This normalizes the pose for the superior-inferior crescent shaped nerve fiber thickness profile as it is consistently located across individuals albeit the individual variability is still present. Hence, before ssLD-DMM was utilized here, the pose of the source and target were already a-priori aligned, which aligns the function signal and benefits ssLDDMM. Without this benefit, if the original pose were to be arbitrary, and absent surface geometry features for registration, as is the case here, the ssLDDMM will return close to identity warps. These warps which would align the geometry but fail to align the function signal on the geometry as this aspect is not taken into account in ssLDDMM registration. Hence, the competing state-of-art of ssLDDMM will fail as compared to the fshape approach. RNFL within 0.25 mm from BMO was excluded, because the layer boundary in the region is often ambiguous. The input surfaces were then concentric with respect to BMO and centrally overlapping, but the outer edges were mismatched due to the varying image sizes. Whereas BMO is an anatomical structure, the image boundary is artificial and depends on the subject's axial length. We were wary that the estimated momenta at the boundary region would be affected by this mismatch. This could also have confounded the classification as the metrics at the boundary region have less spatial correspondence across eyes. However, the regions near the outer boundary also made the least contribution to the classification, as it can be seen in Fig. 5 and7. One way to possibly mitigate the boundary effect is to crop the images, for example, to the largest commonly overlapping region. This however can waste a significant amount of data, and the cropping must be done in a way that does not introduce additional bias.

Choosing a good initial template is important in the computation and quality of the final mean template. The initial template must be topologically equivalent to the observations. An implicit assumption is that the observations are topologically similar. Retinal layer surfaces in OCT images can be modeled relatively simply, as a planar surface in the macular region, and as a planar surface with a hole in the center in the peripapillary region. The initial template can incorporate the general dimensions of the observations. A future work would address the question of the sensitivity of the method to the choice of the initial template. This will be an important issue with more complex structures than retinal surfaces, like lamina cribrosa.

Future work will include creation of population normative atlases and statistics using larger cohorts, and a comprehensive investigation of the optic nerve and macular morphology in a larger data using the fshape and classification modules. The goal is to better detect and understand the shape changes or differences correlated with diseases such as glaucoma. In this report only RNFL was examined, but any other retinal layers can be similarly analyzed, in any combination of layers, surface topology, and signal types. Shape and signal from multiple eyes will be made comparable by a common atlas, and significance of a metric, region, or anatomical structure can be tested by statistical analysis including the classification modules above. In this work the mean template generation and classification process yielded convincing results despite of the mismatch between the RNFL sizes and relative BMO locations. Another direction of future research is to apply the fshape framework to structures like lamina cribrosa, which is not only more topologically complex but also much less consistently visible across eyes than retinal layers. The combined framework may be used to extract common shape information from multiple lamina cribrosa images with varying and limited visibility. Lastly, it will be of value and interest to comparatively investigate different classification methods and optimization techniques for the retinal fshape metrics.

Figure 1 :

 1 Figure 1: (a) 3D visualization of the ONH scan, (b) segmentation of ILM (magenta), posterior RNFL boundary (blue), BM (green), and posterior choroidal boundary (cyan) of a smoothed volume, (c) a posterior RNFL boundary surface color mapped with RNFL thickness.

Figure 2 :

 2 Figure 2: Illustration of the fshape atlas estimation procedure. the estimated mean template (x * , f * ) is a deformation of an initial hypertemplate (x 0 , f 0 ). The hypertemplate transformed to each target belongs to the orbit F of the hypertemplate for example (x 3 * , f 3 * ) = (φ 3 * , ζ 3 * ) • (x * , f * ) ∈ F

Fig. 3

 3 Fig. 3 (b), (c), and (d) respectively show the mean templates generated from all eyes (N = 53), bilaterally normal eyes (N = 20), and glaucomatous eyes (N = 26). Qualitatively, the normal mean template in Fig. 3 (c) displays the

Figure 3 :

 3 Figure 3: (a) Aligned RNFL surfaces with RNFL thickness mapping, (b) mean template of all RNFLs, (c) mean template of normal RNFLs only, and (d) mean template of glaucomatous RNFLs only. Note the low estimated RNFL thickness of the mean glaucomatous template as compared to that of the mean normal template. (e)-(g) show sectoral thickness averages of all RNFLs, normal RNFLs, and glaucomatous RNFLs, respectively. The sectors in (e)-(g) are not in the physical space as the mean templates in (b)-(d), and plotted for visualization with each sector value given by the average RNFL thickness in the corresponding sector across the eyes in the group.

Figure 4 :

 4 Figure 4: Observed RNFLs (top row) and their reconstructions from a common mean template (bottom row). Note that the reconstruction agrees with the pattern of the original RNFL thickness with an overall smooth and noise-reduced profile.

Figure 5 :

 5 Figure 5: Left: Voxel-wise t-test significance map of retinal nerve fiber layer (RNFL) thickness between normal (N = 26) and glaucomatous (N = 27) eyes. The number of vertices is 5355, and the red region indicates Bonferroni corrected p < 0.05/5355 i. e. p < 9.34 × 10 -6 . Right:

Figure 6 :

 6 Figure 6: Pointwise Z-score map of retinal nerve fiber layer (RNFL) fshape residuals for agematched (59.6 ± 6.7) healthy eyes and glaucomatous eyes in early, moderate, and severe stages of the disease, with visual field mean deviation -a measure of glaucomatous vision loss.

Fig. 6

 6 Fig. 6 plots the pointwise Z-score map of ζ for age-matched (59.6 ± 6.7) healthy eyes and glaucomatous eyes in early, moderate, and severe stages of the disease. Visual field mean deviation (VFMD), a measure of glaucomatous vision loss, is also plotted for each eye. The z-score was computed at each point of the mean template by subtracting the mean ζ of the healthy samples and dividing

Figure 7 :

 7 Figure 7: LDA classifier for Healthy vs. Glaucoma data in A) functional residual, B) initial momenta in x-direction, C) initial momenta in y-direction, and D) initial momenta in zdirection.

Accu

  , the ssLDDMM algorithm returns a geometric map close to identity between the given source and target, whereas the fshapes returns a geometric map further from the identity indicating that it was driven by the function signal for creating overlap. This has a direct impact on the error remaining after mapping as well (as shown in the third column). The geometry-based transformation of source signal to target derived from ssLDDMM shows high values of remaining error after mapping -the error profile shows the presence of two unregistered function signals as the transformed source function is subtracted from the target signal. In comparison, by incorporating the function signal in the fshapes registration, the function profiles of the source and target are matched leading to a lower error profile obtained by fshape approach. A typical effect on the atlas estimation is shown in Figure 9 in which mean templates and group difference significance maps for the older normals and older glaucoma groups are compared. The mean function estimated using fshapes is sharper and retains more high frequency features of the population since the functions are aligned, whereas the ssLDDMM mean is more blurry and retains less of the function profiles since the geometry-only registration does not incorporate the alignment of functional signals. An important point to note is that the geometry-only-based

Figure 8 :

 8 Figure 8: The first column shows a source RNFL thickness function and a target RNFL thickness function. Top row shows results with standard-surface LDDMM (ssLDDMM) and bottom row shows results using the fshape algorithm. The ssLDDMM algorithm returns a geometric map close to identity to map the given source and target, and the remainder after transforming the functions from this map show considerable residue that this mapping is unable to register. The bottom row show the same experiments with fshapes, where the residue left after mapping the source to the target function is smaller and more diffuse since the overall functions have been matched.

565 4 .

 4 Discussion and future workWe presented a novel application of the fshape framework for variability analysis of shape and associated signals in retinal optical coherence tomography (OCT) images, comparing healthy, glaucomatous, and suspect peripapillary

Figure 9 :

 9 Figure 9: Top row shows results with standard-surface LDDMM (ssLDDMM) and bottom rowshows results using the fshape algorithm. The first column is the group average of healthy individuals, the second is group average of individuals with Glaucoma and the third column shows t-statistic from group difference using these two methods. Note that the ssLDDMM results are more blurry due to mis-registration of the functional signal, and the group differences are missing the superior part of the crescent shape of pathological changes as well as overall being more diffuse and less strong everywhere.

Table 1 :

 1 The classification result of the bilaterally healthy eyes (N = 19, mean age: 43.4 ± 14.8) and suspect eyes (N = 7, mean age: 57.1 ± 12.4) are shown in Table1.That the suspect eyes are distinguished from the bilaterally healthy eyes is noteworthy, considering that the suspect eyes are pre-diagnosis and without functional loss or other conventional clinical features of glaucoma. As expected, the classification rates are lower than those between confirmed glaucomatous eyes and bilaterally healthy eyes in the previous experiment, but still relatively high. The result suggests the fshape metrics may capture some morphological changes in the RNFL that precedes vision loss in glaucoma. Accuracies, sensitivities, and specificities of the classification of healthy, glaucomatous, and suspect RNFLs based on the fshape metrics of RNFL posterior surface geometry and RNFL thickness. Both Healthy vs. Glaucoma and Healthy vs. Suspect show high classification success rates.

		Accu. (%) Sens. (%) Spec. (%)
	Healthy vs. Glaucoma	92.9	94.4	90.0
	Healthy vs. Suspect	88.5	71.4	94.7

Table 2 :

 2 Accuracies, sensitivities, and specificities of the classification of healthy, glaucomatous, and suspect RNFLs based on RNFL thickness sectoral averages. In comparison with the results in Table1with RNFL fshape metrics, the classification performance is lower.

	3.3.4. RNFL thickness vs. RNFL surface geometry
	In order to compare the discriminating power of RNFL thickness against
	RNFL posterior surface geometry, and to confirm what was shown in Fig. 7,
	we repeated the Healthy vs. Glaucomatous classification with RNFL functional
	residual (RNFL thickness) and initial momenta (RNFL posterior surface ge-
	ometry) separately. The result is summarized in Table 3. As expected from
	Fig. 7, the classification result with RNFL thickness only is comparable to that
	with both RNFL thickness and geometry, and superior to the result with RNFL
	geometry only.

Table 3 :

 3 Accuracies, sensitivities, and specificities of the classification of healthy and glaucomatous RNFLs by RNFL thickness only (functional residual, top row) and RNFL posterior surface geometry only (geometrical momenta, bottom row) (BMO-based rigid pre-registration)Prior to the fshape template generation, the RNFL surfaces were rigidly registered by the BMO centroids. The fshape geometrical deformation then contains information for both RNFL surface topology, and the distance between RNFL and BMO. The experiment was repeated with rigid pre-registration by aligning the central opening of the RNFL surface instead of the BMO centroid, and the result is summarized in Table4. The functional-based classification is comparable; however, the geometry-based classification is worse in that most of the glaucoma cases are misclassified. This suggests that the factors that deter-

	mine the distance between the RNFL and BMO, such as the optic canal skew
	or post-RNFL retinal thickness, may be affected by glaucomatous structural
	change, more so than RNFL posterior surface topology alone.
		Accu. (%) Sens. (%) Spec. (%)
	Functional residual	96.4	94.4	100.0
	Geometrical momenta	53.6	33.3	90.0

Table 4 :

 4 Accuracies, sensitivities, and specificities of the classification of healthy and glauco-

	matous RNFLs by RNFL thickness only (functional residual, top row) and RNFL posterior
	surface geometry only (geometrical momenta, bottom row) (RNFL opening-based rigid pre-
	registration)

Acknowledgment This work was supported by Canadian Institute of Health Research(CIHR), Natural Science and Engineering Research Council of Canada (NSERC), Michael Smith Foundation for Health Research (MSFHR), Alzheimer Society Research Program (ASRP), and Brain Canada.


Sandrine Rosin-Paumier

Farimah Masrouri

Evaluating the impact of thermal variations on the penetration test parameters

INTRODUCTION

Thermo-active geostructures consist in adding to the classical sustainable function of geostructure a second function of thermal exchange with the ground. These systems use the subsurface as a heat source or as a cooling source. A liquid circulates through a ground circulation closed-loop system integrated to the geostructure. Geostructures may be piles, barrettes, diaphragm walls, tunnel linings, basement slabs or walls [START_REF] Fromentin | Pieux échangeurs : Conception et règles de pré-dimensionnement[END_REF][START_REF] Laloui | Comportement d'un pieu bi-fonction, fondation et échangeur de chaleur[END_REF][START_REF] Brandl | Energy foundations and other thermo-active ground structures[END_REF]. The standard use of thermoactive foundations modifies the current ground temperature (12°C) in a scale of 1 to 30°C [START_REF] Peron | A method for the geotechnical design of heat exchanger piles, Proceedings of the Geo-Frontiers Conference[END_REF]. Geo-exchange systems in ground structures require consideration of the complex interaction between temperature change and induced stresses or deformations, which may affect building performance [START_REF] Mccartney | Application of Geoexchange Experience to Geothermal Foundations[END_REF]. Foundation piles work thanks to the combination of two effects: the cone tip resistance and the sleeve friction along the pile. The sleeve friction can be associated with the shear resistance. The shear resis-tance is the relative resistance of the ground that goes against the sliding caused by a load. The impact of temperature on shear resistance is complex. Some authors such as [START_REF] Mitchell | Shearing resistance of soils as a rate process[END_REF], [START_REF] Sherif | Temperature effects on the unconfined shear strength of saturated, cohesive soil. Effects of Temperature and Heat on Engineering Behavior of Soils[END_REF], [START_REF] De Bruyn | The influence of temperature on mechanical characteristics of Boom clay: The results of an initial laboratory program[END_REF], [START_REF] Lingnau | Effects of temperature on strength and compressibility of sand-bentonite buffer[END_REF] and Wiebe et al. (1998) showed that ground heating leads to a decreasing of shear resistance whereas [START_REF] Noble | Effect of temperature on strength behavior of cohesive soil[END_REF], [START_REF] Houston | Thermo-mechanical behavior of seafloor sediments[END_REF], [START_REF] Tanaka | Stressstrain behaviour of reconstituted illitic clay at different temperatures[END_REF] and Cekerevac and Laloui (2004) noticed the opposite results. In fact, heating and loading cause opposite effects. On the one hand, the heating causes decreasing of the viscosity of the interstitial water. In consequence, the resistance between aggregate decreases as well as the shear resistance. On the other hand, the heating of a material under a mechanical stress causes a densification of the ground which results in to increase contacts between aggregates. In consequence, the resistance between aggregate increases as well as the shear resistance. So the two kinds of conclusions found in published works are due to a combination of these two effects in different soils. Considering that thermo-active foundations may be built in a large extent of soil nature and given the difficulty in forecasting how a temperature variation Evaluating the impact of thermal variations on the penetration test parameters.

G. Franchomme, S. Rosin & F. Masrouri LEMTA (CNRS, UMR 7563) 2 avenue de la Forêt de Haye,BP 160,France ABSTRACT: Geostructures or thermo-active foundations consist in establishing a thermal exchange with the ground directly through building foundations or other structures buried in the soil. Geostructures may extract or inject energy from or into the ground to produce heating or cooling in a building. Their cost is currently about 30% lower than the usual geothermal solutions using vertical probes. However, adapting foundation piles in whole or in part to energy piles brings questions on the impact of the temperature variations on the geotechnical parameters and load-bearing capacity of the soil. The bearing capacity calculations of deep foundations are currently based on penetrometer testing results. In this study, a high resolution mini-penetration test, adapted from current static penetration test, was developed to characterize laboratory samples. A kaolinite/sand mixture was prepared to test a material in an intermediate state between a cohesive and a non-cohesive soil. Samples were performed at various water contents (9 to 19%) and dry densities (1.58 to 1.91 Mg/m 3 ) and submitted to linear and cyclic thermal stresses (1 to 50°C). Results showed a high repeatability of these tests. Two parameters were measured: the cone tip resistance (q c ) and the sleeve friction (fs). These parameters highly depend on the water content of the material. We showed that the temperature variations with constant water content did not affect the measured shear stress parameters. The proposed penetration tests can be carried out as a quick and easy-to-use test to compare the impact of thermal stresses on the soil parameters.

can affect the shear resistance of a material, it seems necessary to develop a quick and easy-to-use test to compare the impact of a thermal stress on the soil parameters. In this study, a high resolution minipenetration test was developed. The method supplies parameters close to the penetrometer test currently used to calculate the bearing capacity of piles. The aim is to define quickly the impact of a thermal stress on the shear resistance of different materials. In this study, the mini-penetration test was first presented then tested on a kaolinite and fine sand mixture. Several samples were prepared at the proctor optimum and and submitted to temperature from 1 to 50°C, which included the maximal variation of the soil temperature for a current use of thermal piles. Then, the main results were summarized. Finally, conclusions were made on the impact of the tested temperatures on the chosen material and the possible prospects of the use of the mini-penetration test.

MATERIAL AND METHODS

Preparation and basic characterization of the samples

The tested material was a mixture of 50% of kaolinite Argirec B24 from Charente basin and 50% of Hostun sand HN 38 (France), prepared in laboratory.

The materials were dried separately at least 24 hours in an oven (60°C) then mechanically mixed in desired proportion. Water was added to obtain different initial water contents comprised between 9 and 19%.

Then the mixture was packed in hermetic bags and left for homogenization at least 24 hours (Table 1). The liquid limit of the material (L L ) was 26.9%, the plastic limit (P L ) was 15.4% and the plasticity index (PI) was 11.5. The particle size distribution was determined using a Laser Fritsch analysette 2© (Figure 1). The material was compacted in 6 layers in a mold of 152 mm in diameter and 152 mm in height. The applied energy was the same than that used in the Proctor method. The compaction in 6 layers permitted to ensure the homogeneity of the material in the mold (Table 1). The obtained compaction curve was superimposed to the Standard Proctor curve (Figure 2). The optimum water content was 13% and the maximum dry density was 1.88 Mg/m 3 (Figure 2).

Thermal cure

A thermal curing at 1, 30 or 50 °C was applied to a sample equipped with a thermal sensor to determine the minimal duration of thermal curing (Table 2). The initial temperature of the sample was 20°C. 113, AFNOR, 1996). Our test consisted of pushing a probe with a cone having a 60° apex continuously at 30 mm/s into each sample (Figure 3).Two stainless steel probes, respectively named P1 and P2 were used (Figure 3). The length (250 mm) as well as the angle (60°) and the maximal external diameter (12 mm) of the cones are similar. Diameters of the probes differ from 10 mm for P1 to 12 mm for P2.

The maximum penetration length of each probe is 100 mm. The cone tip resistance was measured with P1, whereas using P2, a combination of the cone tip resistance and the sleeve friction was obtained.

At the equilibrium temperature, 3 mini-penetration tests were carried out in each mold. The position of each penetration test is fixed to avoid edge effects (Figure 4). Three penetration velocities were tested: 1 mm/min, 30 mm/min and 90 mm/min. The velocity of the probe determined the test duration. This duration should be minimized to avoid the changes of water content of the samples. Preliminary tests showed noimpact of probe penetration velocity on the measured parameters. For that reason an intermediate speed of 30 mm/min.

Figure 4: Position of penetration tests in the mold.

RESULTS

Mini-penetration curves typology

The penetration curves obtained by P1 and P2 probes were respectively named E and E* tests. The variation of the total strength versus depth is shown on Figure 5. In the first part of the curves (a), the apex of the probe penetrated into the soil, this part is similar for two probes as the angle and the diameter of the cone apex are identical. Part b is peculiar to the use of P1. As the friction sleeve was not measured, the strength variation was constant for a homogeneous soil. Part c is peculiar to P2. The strength increased as the probe enters in contact with the soil. The cone tip resistance (q c ) was calculated by averaging the strength values obtained on the linear part (b) of P1 test. Equations 1 and 2 allowed the transformation of the measured stresses in a resistance.

q c = F c / (πR 1 ²) (1) fs = (Ft -Fc)/( 2πR 2 (L-Lp)) (2)
where Fc is the cone tip strength ; R 1 is the radius of the cone apex, R 2 is the probe radius; fs is the sleeve friction; L is the penetration length and Lp is the height of the cone apex.

Coupled influence of water content and dry density at 20°C

This part presents the results of the tests carried out at 20°C on the samples described in Table 1. The samples E3 and E8 were prepared at the same initial dry density (1.82 Mg/m 3 ). The difference of 750 N between the obtained strengths could be only due to the differences in water contents: for E3, w = 11.9 % and for E8, w= 15.1%. Generally it could be noticed that the variation of the dry density has a lower impact on these results.

Cone tip resistance

The Figure 7 superimposes the cone tip resistance at 90 mm penetration depth to the compaction curve.

The lower the water content, the higher the tip resistance. On the other hand the decrease of the cone tip resistance is sharper at the dry side (w<14.5%, equation 3) than the wet side (w>14.5%, equation 4). The equations of these two curves were defined in order to be used as a reference for the following tests. q c = -2028.2w + 29770

(3) q c = 49.589w² -1875.2w + 17919 (4)

Sleeve friction

For each sample, a test with P1 and a test with P2 were performed. The difference between the two curves gave the sleeve friction (Figure 8). The comparison of sleeve friction at different couples of water content and dry density let us draw the same conclusion than for the cone tip resistance. The decrease of the water content of the samples increased the sleeve friction (Figure 9). The equations of the two curves for sleeve friction are showed in the equation 5 for water contents from 12% to 14.5% and in the equation 6 for water contents varying from 14.5% to 19. fs = -74.353w + 1113.9

(5) fs = 1. 5969w² -59.743w + 571.35 (6) The results of penetration tests carried out on different samples (E11 to E13) after a thermal curing of 1, 30, and 50°C were compared to the results of the same test at 20°C (E6) (See Table 1 and3). The cone tip resistance and the sleeve friction did not greatly change with the temperature variations between 1 to 30°C. However an increase of the temperature to 50°C, increased dramatically these two parameters (Figure 10). This phenomenon could be due to the fact that when heating the samples small water content changes occurred from 13.2% to 12.5% despite the complete covering of the mold during the thermal curing. As already shown in Figure 7 and Figure 9 (equation 3 and 5), this decrease in water content could increase dramatically both the cone tip resistance and sleeve friction.

Test E4 performed at 20°C on a sample with similar water content and dry density than E13, showed the same results (Table 3).

CONCLUSION

The purpose of the study was to develop a high resolution mini-penetration test to propose an easy-touse method to define quickly the influence of different temperatures on cone tip resistance (q c ) and sleeve friction (fs) of a soil. The results showed a good repeatability and reliability of the test. The study was divided into several phases: (1) definition of an adapted compaction protocol, (2) definition of the thermal cure protocol, (3) definition of penetration test conditions (depth, velocity, …). The effects of water content and dry density and thermal cure on the penetration test parameters were studied. We showed that the decrease of the water content of the samples increased the cone tip resistance and the sleeve friction. The variation of the dry density of the sample has a lower impact than the variation of its water content on the penetration parameters. The decrease of the cone tip resistance and the sleeve friction is nevertheless sharper at the dry side of the compaction curve (w<14.5%, equation 3) than at the wet side (w>14.5%, equation 4). The whole results let us conclude on the small influence of thermal cure between 1 and 30°C on penetration parameters for the mixture kaolinite/Hostun sand. Thermal curing at higher temperature (50°C) lead to a decrease of the sample water content and, as a result, an increase in cone tip resistance and sleeve friction.

Figure 1 :

 1 Figure 1 : Particle size distribution of kaolinitesand mixture.

Figure 2 :

 2 Figure 2 : Comparison between the Proctor curve and the specific compaction curve of the material. Sr : Saturation rate.

Figure 3 :

 3 Figure 3: a. Mini-penetration test with tractioncompression machine and probes for minipenetration test (b. P1; c. P2).

Figure 5 :

 5 Figure 5: Total strength versus depth with P1 (E) and P2 (E*) probes.

Figure 6

 6 Figure 6 shows the evolution of P1 penetration strength for tests carried out in samples with different values of water content and dry density. The b parts of the curves are constant indicating a good homogeneity of the samples. The maximal strength obtained during the tests increased with the decrease of the water content of the samples.The samples E3 and E8 were prepared at the same initial dry density (1.82 Mg/m 3 ). The difference of 750 N between the obtained strengths could be only due to the differences in water contents: for E3, w = 11.9 % and for E8, w= 15.1%. Generally it could be noticed that the variation of the dry density has a lower impact on these results.

Figure 6 :

 6 Figure 6: Strength measured along the penetration of P1 in samples at different water contents and dry densities.

Figure 7 :

 7 Figure 7: Comparison of the cone tip resistance at 90 mm penetration depth with the compaction curve

Figure 8 :

 8 Figure 8: Sleeve friction calculated from the penetration tests with P1 and P2 at different water contents and dry densities

Table 1 :

 1

	E1	9.3	1.60	20	9.3
	E2	11.3	1.74	20	11.3
	E3	11.9	1.82	20	11.9
	E4,E*4	12.3	1.87	20	12.3
	E5, E*5	13.0	1.88	20	13.0
	E6	13.2	1.89	20	13.2
	E7, E*7	14.2	1.85	20	14.2
	E8, E*8	15.1	1.82	20	15.1
	E9, E*9	16.9	1.73	20	16.9
	E10, E*10 18.5	1.66	20	18.5
	E11, E*11 13.5	1.90	1	13.4
	E12, E*12 13.5	1.90	30	13.2
	E13	13.4	1.87	50	12.5

Water content determines after the compaction of the sample (w i ) and after the penetration test (w f ), dry density (ρ d ) and temperature (T) of each sample.

Test w i (%) ρ d (Mg/m 3 ) T (°C) w f (%)

Table 2 :

 2 Parameters of the thermal cure protocol

	Climatic chamber	1	30	50
	temperature (°C)
	Duration of the	9h20	4h20 12h50
	thermal cure (h)
	Sample temperature (°C) 1.5	29.3 48.8
	2.3 Mini penetration test
	The mini-penetration test is a laboratory test adapted
	from the normalized static penetration test (NF P94-

Table 3 :

 3 Measured and theoretical cone tip resistance (q c ) and sleeve friction (f s ) at 90 mm penetration depth

	Sample number	E11 E6	E12 E13
	Temperature (°C)	1	20	30	50
	q c (KPa)	2323 2504 2336 4075
	f s (KPa)	105 131* 106 317
	*Calculated with equation 5