Construction of a fetal spatio-temporal cortical surface atlas from in utero MRI: Application of spectral surface matching

In this study, we construct a spatio-temporal surface atlas of the developing cerebral cortex, which is an important tool for analysing and understanding normal and abnormal cortical development. In utero Magnetic Resonance Imaging (MRI) of 80 healthy fetuses was performed, with a gestational age range of 21.7 to 38.9 weeks. Topologically correct cortical surface models were extracted from reconstructed 3D MRI volumes. Accurate correspondences were obtained by applying a joint spectral analysis to cortices for sets of subjects close to a speci ﬁ c age. Sulcal alignment was found to be accurate in comparison to spherical demons, a state of the art registration technique for aligning 2D cortical representations (average Fréchet distance ≈ 0.4 mm at 30 weeks). We construct consistent, unbiased average cortical surface templates, for each week of gestation, from age-matched groups of surfaces by applying kernel regression in the spectral domain. These were found to accurately capture the average cortical shape of individuals within the cohort, suggesting a good alignment ofcorticalgeometry.Each spectralembedding anditscorrespondingcorticalsurfacetemplateprovideadualref-erence space where cortical geometry is aligned and a vertex-wise morphometric analysis can be undertaken. © 2015 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The cerebral cortex, a thin layer of neural tissue lining the cerebrum, plays a key role in high level brain functions such as memory, sensory perception, language and consciousness (Kandel, 2013).Cortical folding is an important neurodevelopmental process that occurs largely in the third trimester of pregnancy, where the surface of the cerebral cortex transforms rapidly from a smooth sheet into a highly convoluted arrangement of gyri and sulci (Garel et al., 2001).The increase in surface area allows more neurons to be packed into the skull's limited volume and is fundamental to increases in cognitive ability seen in humans and other intelligent species (Lui et al., 2011).It is thought that abnormal cortical folding could lead to neurocognitive deficiencies.For example, Dubois et al. (2008) found that fetuses with intra-uterine growth restriction, which is a risk factor for developing attention deficit hyperactivity disorder and schizophrenia (Geva et al., 2006), had measurable structural abnormalities at birth that were associated with a lower neurobehavioural development score at term equivalent age.
Our previous work focussed on quantifying cortical development, at a global and lobar scale, by computing a single scalar value, summarising the average extent of folding over a given region (Wright et al., 2014).
However, disturbances in development may be more localised and anomalies may remain undetected when averaging over a relatively large region.We therefore seek to characterise folding at any given point on the cortex.
An average model or atlas can facilitate this by providing a densely sampled reference space where the anatomies of a population can be co-aligned and morphometric analysis performed.Magnetic Resonance Imaging (MRI) based atlases have been used previously to align and compare neuroanatomy between cohorts on a voxel-wise basis using voxel-based morphometry (Ashburner & Friston, 2000).This methodology can also be applied to surfaces, undertaking morphometric analysis on a vertex-wise basis using mesh models (Chung et al., 2008).The aim in this paper is to develop an atlas construction framework to generate a set of accurate average cortical surface templates, which will give a longitudinal reference space for surface morphometry.This will enable us to quantify normal and abnormal folding patterns of the cerebral cortex in future work as well as provide a reference space for subsequent work on connectivity analysis in the developing brain.
Constructing a surface atlas requires finding smooth and accurate correspondences between cortical surfaces from different subjects at different ages, which can be a challenging and computationally intensive task.Previous registration based methods for aligning adult cortical surfaces have required some form of preprocessing to produce a 2D surface parameterisation.This is advantageous as a smooth correspondence between cortical surfaces can be estimated more robustly, while reducing the optimisation problem to two dimensions and decreasing the computation required.For example, Van Essen and Drury (1997) flattened the cortical surface by making a series of surface cuts to reduce distortion before warping one flat cortical map to another.Popular methods such as Freesurfer (Fischl et al., 1999), spherical demons (Yeo et al., 2010) and multi-modal surface matching (Robinson et al., 2014), adopt a 2D spherical coordinate system.Through a process of inflation, a surface can be transformed into a sphere while minimising geometric distortion.Surfaces can then be treated as 2D spherical images whose intensity is given by scalar information associated with each vertex (e.g.sulcal depth) to guide the registration.
Spectral graph theory provides a fast alternative to these methods using the eigenvalues and eigenvectors of the graph Laplacian matrix (Chung, 1997).The eigenvalues of the Laplacian may be interpreted as the natural oscillating frequencies, or resonances, of a physical shape.
Each of these resonances has a distinct pattern of displacement, a socalled normal mode of the shape, which is given by the eigenvectors.Equivalent locations on similarly shaped surfaces will have a similar relative amplitude and phase, at any given resonant frequency, which allows surface points to be matched in the spectral domain.Accuracy with spectral methods can be limited due to macroscopic variations in brain shape such as expansion or compression that can cause significant changes in a cortical surface's spectral representation.Previous attempts to address this problem include registering the derived spectral representations (Jain & Zhang, 2006;Mateus et al., 2008;Lombaert et al., 2011Lombaert et al., , 2013a)).
More recently, Lombaert et al. (2013b) proposed a novel spectralbased method to acquire a mapping between surfaces, that is comparable in accuracy to Freesurfer (Fischl et al., 1999) and spherical demons (Yeo et al., 2010), while maintaining a considerable speed advantage.Moreover, atlases generated are more consistent when varying the initial reference surface for aligning a group of surfaces.The standard deviation of atlas boundaries was found to be significantly lower for their spectral method (0.0014 mm ± 0.0009) compared to spherical demons (0.28 mm ± 0.08).In this method, a preliminary point-wise correspondence between two surface graphs is established from their spectral representations and additional information such as sulcal depth.A novel dual-layered graph is then constructed, by linking the corresponding vertices across the two graphs.The spectral decomposition of this duallayered graph provides an orthonormal basis where accurate correspondences between the two cortical surfaces can be estimated.Orasanu et al. (2014) used this methodology to analyse cortical development in neonates born very pre-term (around 26 weeks).Cortical surfaces for five subjects were modelled at two time-points, 30 weeks and 40 weeks equivalent gestational ages.For each individual, cortical correspondences were then estimated between the two time-points, enabling longitudinal changes in cortical geometry to be measured at matched locations.
Given the considerable diversity of cortical shape during gestation, a single anatomical atlas may not be appropriate.Indeed, a sulcus on the surface of an older more neurologically developed brain may not yet exist on a younger brain and thus establishing a meaningful correspondence may not be possible.Without any detectable surface sulcation, from some shape descriptor (e.g.mean curvature), correspondences can be determined by the spectral properties of the surfaces alone.Unfortunately, variations in cerebral shape can alter the spectral vibration modes limiting matching accuracy.Furthermore, macroscopic changes in cerebral shape during gestation alter vibration modes over time, therefore accuracy will also fall as the discrepancy in age increases.A spatio-temporal approach (Serag et al., 2012) can be used to avoid these difficulties, constructing an average anatomical template at several time points, within a specific time interval, where anatomy is similar.
To date, there have only been a handful of studies which have constructed atlases of the developing cortex.Hill et al. (2010) constructed a single cortical surface template from neonates born at term to compare hemispheric asymmetries in infants and adults.A landmark based registration algorithm (Van Essen, 2005) was used to align spherical representations of the cortex.Li et al. (2014) recently constructed an infant spatio-temporal surface atlas, generating atlas templates for 1, 3, 6, 9, 12, 18 and 24 months of age.This work also aligned cortices on the sphere, using spherical demons (Yeo et al., 2010).
A significant portion of the cortical folding process is completed before term, including the formation of primary and secondary sulci (Garel et al., 2001) and disturbances in development during this period may have a profound effect.Preterm neonatal imaging therefore provides an opportunity to visualise an important part of the cortical folding process.However, the complications associated with premature birth mean that this cohort is not ideally suited for characterising normal development (Ajayi-Obe et al., 2000;Kapellou et al., 2006).For this reason, in utero MR imaging may ultimately provide improved visualisation and quantification of the developing cortex over this critical period.Previously, fetal motion has greatly hindered the acquisition of coherent MR images in utero.However, recent advances in image reconstruction (Rousseau et al., 2006;Jiang et al., 2007;Rousseau et al., 2010;Kim et al., 2010;Gholipour et al., 2010;Kuklisova-Murgasova et al., 2012) have allowed the acquisition of high resolution 3D volumes that are free from motion artefacts.
Two studies have constructed average surface templates of the fetal cortex from in utero MRI.Habas et al. (2012) described an imageanalysis framework to investigate early folding patterns revealing asymmetric hemispheric development in a cohort of fetuses, from 20 to 28 weeks gestational age.MR image volumes were co-aligned at each week of gestation using image-based registration and an average cortical surface was extracted from the combined subject tissue maps.3D image-based registration is, however, ill-suited for cortical alignment at older gestational ages, where the geometry of the cortex is significantly more complex.
More recently, Clouchoux et al. (2012) constructed a spatio-temporal surface atlas consisting of four templates, based on twelve in utero MR images of healthy fetuses (25.2 to 35.1 weeks gestational age), enabling a quantitative analysis of cortical development.Surfaces were reconstructed using a deformable model, CLASP (Kim et al., 2005), and correspondences were found on the unit sphere using an iterative group registration scheme (Lyttelton et al., 2007).This methodology was also used to compare development in normal fetuses and fetuses with a severe form of congenital heart disease (Clouchoux et al., 2013), where a delay in cortical development was observed.
Most of the atlas construction techniques discussed have utilised a spherical coordinate system to register cortical surfaces.In this study, we present an alternative spectral method for constructing a spatiotemporal surface atlas of the developing fetal cortex.We generate average surface templates for each week of gestation from 23 to 37 weeks (15 templates in total), from a cohort of 80 healthy fetal subjects.To our knowledge there are no publicly available surface atlases for the developing fetal cortex over this age range built from such a large population, therefore we will make our atlas available to assist future brain development studies.
In the next section, we outline our spatio-temporal framework and template construction methodology.We describe a natural extension to the spectral-based method of Lombaert et al. (2013b) to embed a group of cortical surfaces simultaneously, thus eliminating bias due to the choice of an initial reference surface for alignment.We then demonstrate the use of kernels (Nadaraya, 1964), to compute a spatially and temporally weighted estimate of the average surface boundary, and show how this is readily converted to a mesh representation (Kazhdan & Hoppe, 2013).
We evaluate our methodology by examining sulcal alignment, regional overlaps, and embedding regularity in comparison with spherical demons and conventional 3D surface registration algorithms.We also demonstrate the consistency of the average templates by generating multiple instances from disjoint subsamples of our dataset of varying sizes.The implications of this work for characterising normal and abnormal cortical growth are then discussed alongside the possible directions for future work.

Data and surface modelling
Eighty healthy control subjects with a gestational age (GA) range of 21.7 to 38.9 weeks (mean 29.6 ± 4.6) were used for the study.In utero MR imaging was performed on a 1.5 Philips Intera scanner using a 32 channel coil.Individual slices were acquired using a singe-shot fast spin echo sequence which were then reconstructed into a 3D volume using a robust slice-to-volume registration method (Jiang et al., 2007).Reconstructed volumes had an isotropic voxel resolution of 1 mm.For more details on data preprocessing please see Wright et al. (2014).
As with our previous work (Wright et al., 2014), for the cortical modelling, we focussed on extracting the cortical grey matter (GM)/white matter (WM) surface as this is easier to delineate.MR images were first segmented into basic tissue types using an expectation maximisation (EM) algorithm (Makropoulos et al., 2014), including cortical GM, WM, cerebrospinal fluid (CSF), the deep grey nuclei, lateral ventricles, brain stem and cerebellum.The posteriors obtained from the EM segmentation were then combined to form a probability map for matter lying inside the cortical GM/WM boundary.An explicit surface representation was extracted by deforming a mesh with spherical topology towards the implicit surface, using a similar approach to Kim et al. (2005).The deformation maintains the topology of the mesh and fits the mesh directly to the probability map rather than the image intensities (see Appendix A for more details).All meshes were tested and confirmed to be topologically correct with no selfintersections.To estimate the surface extraction error, a distance map was computed for a segmented boundary (Maurer et al., 2003), and distance values were linearly interpolated at the mesh vertex locations.The average distance error for a vertex was 0.106 mm (±0.006) (see Fig. 1 for surface contour examples).

Framework overview
Given the temporal diversity of anatomy within the cohort, a spatiotemporal framework was adopted, building an average template of cortical shape for each week of gestation across the cohort (Fig. 2).In particular, the approach avoids the difficult task of aligning very disparate cortical anatomies and trying to find accurate correspondences between them.Each template was built from a group of cortical hemisphere surfaces (N = 19.6 ± 3.1) within a narrow time-window of ± 1 week, where cortical features are similar.Both hemisphere surfaces were used to generate the atlas templates (mirroring the right hemisphere), effectively doubling the sample size, allowing a more reliable estimate of the population average surface boundary.By choosing overlapping windows, correspondences between successive embeddings can still be established (Temporal correspondence section).For each template, the contributing surfaces were first affine-aligned by registering all pairs of surfaces, and deforming each surface using the computed average affine transformation.This improved matching accuracy by removing differences in global shape and size.Spectral matching was then performed between all pairs of surfaces, yielding initial vertex-wise correspondences.Surfaces were connected using these preliminary correspondences and a joint spectral analysis gave an embedding in which all the surfaces were co-aligned.The average surface position was then estimated using kernel regression at each point within the embedding, allowing a mean surface to be reconstructed.The following describes these steps in further detail.

Spectral matching
Spectral matching is a basic concept that has been used by several authors (Jain & Zhang, 2006;Mateus et al., 2008;Lombaert et al., 2011Lombaert et al., , 2013a,b;,b;Orasanu et al., 2014) to find correspondences between two surfaces from their spectra.To recap, we first define a surface model S as a pair S = {V, E}.V denotes a set of vertices V = (x 1 ,…, x |V| ) with position vectors x i , while E denotes an adjacency matrix, where e ij = 1, if vertices i and j form an edge, or 0 otherwise.From this we define the |V| × |V| weighted adjacency matrix W, is then given by summing all edge weights for each node, d i = ∑ j w ij .The |V| × |V| general graph Laplacian L = D − W can now be constructed from W and D.
Solving the eigenvector problem L = UΛU −1 gives the diagonal eigenvalue matrix, Λ = diag(λ 0 , λ 1 , …, λ |V | ) and the corresponding eigenvector matrix U = (u 0 , u 1 , …, u |V | ), where each eigenvector Each of the eigenvectors or eigenmodes, u i , can be viewed as a pattern of vibration for the surface model S (Fig. 3), associated with the resonant frequency λ i .The first eigenvector is not a vibration mode as λ 0 = 0 and is thus discarded.These eigenmodes give a useful parameterisation which can be used to match points with another similarly shaped surface.For correct matching, eigenmodes may need to be corrected for their sign ambiguity, multiplicity, and perturbation in isometry (Lombaert et al., 2011(Lombaert et al., , 2013a)).
To match points between surfaces the concept of a spectral representation is useful (Fig. 4).The spectral representation of S is denoted as a kdimensional embedding where the ith vertex has the spectral coordinates We define a correspondence map c which gives the point-wise correspondences between two surfaces S I and S J .For each vertex x i ∈ S I with spectral coordinates u i ′, the corresponding vertex y j ∈ S J with spectral coordinates v j ′ is determined by the shortest Euclidean distance in the spectral domain, To achieve a more accurate match, surface descriptors such as mean curvature or the surface normal direction may be added to the embedding as extra dimensions (Lombaert et al., 2011).This ensures that a sulcus (or gyrus) from one surface matches with a nearby sulcus (or gyrus) from another surface.The values for each descriptor may be scaled in order to weigh their influence when finding correspondences.
We have found, however, that this nearest neighbour approach to define an initial correspondence map c between two surfaces Fig. 2. Framework.An average cortical surface template is constructed for each week of gestation, with all subjects within a week of the target age contributing to the output.A spectral analysis yields spatial correspondences between a group of cortical surfaces, allowing the average surface position to be computed and a template surface constructed.Note the red/blue surface colour mapping depicts the mean curvature of the surfaces.
(Spectral matching Section) lacks regularisation and the matched locations may be incoherent, i.e. two vertices that are close on one surface may match with vertices that are distant on another surface.This can cause irregularities in the embedded surfaces (Embedding regularity section), particularly when incorporating scalar information into the embedding, which is not determined by surface position (e.g.mean curvature), or where scalar data may be noisy.
To address this, we used edge-based smoothing (Field, 1988) in the spectral domain to regularise the connections between two surfaces S I and S J given an initial correspondence map c = c IJ (Algorithm 1).For each vertex x i ∈ S I , with corresponding vertex y c(i) ∈ S J and corresponding spectral coordinate values v′ c(i) , we first obtain the smooth spectral coordinate v 0 i .This is achieved by taking a weighted average of the neighbouring spectral coordinates v′ c(j) , based on the edge set A more regular correspondence map can then be defined Algorithm 1. Regularise correspondences

Joint spectral analysis
In Lombaert et al. (2013b), a novel dual-layered surface graph is built from two meshes interconnected by a preliminary correspondence map generated by spectral matching (Spectral matching section).A joint spectral analysis of the connected surfaces produces a set of shared eigenmodes as opposed to two independent sets of eigenmodes.This produces a shared parameterisation of the surfaces.
We extend this approach by building a combined graph for k subjects with surface graphs S 1 ,…, S k , allowing a joint spectral analysis of all surfaces simultaneously.A multi-layered surface graph ) is a concatenation of all mesh vertices and E c is the combined connectivity matrix.E c preserves the intrasurface vertex connectivity defined in surface graphs S 1 ,…, S k as well as incorporating additional inter-surface connections.For each pair of surfaces S I and S J , an edge is formed between each vertex x i ∈ S I and its corresponding vertex y c(i) ∈ S J , where c = c IJ is a preliminary correspondence map between S I and S J .As surfaces are linked pairwise, the number of inter-surface edges may be far greater than the number of intra-surface edges.To compensate for this we introduce a parameter ϕ to control the relative weighting of intra-and inter-mesh edges.Intra-mesh edges are weighted as before, by the inverse squared Euclidean distance, w = ‖x i − x j ‖ −2 .For inter-surface connections we choose the weighting w . Note ϕ = 1 gives an equal weighting between intra-surface and inter-surface edges (see Appendix B for further details).
The spectral decomposition of the graph Laplacian built from the multi-layered surface graph S c , gives a set of shared eigenmodes U c and an orthonormal basis where surfaces S 1 ,…, S k are co-aligned (Fig. 5).

Kernel regression
The coordinate system provided by the shared eigenmodes (obtained from a joint spectral analysis of a group of surfaces) is referred to as the spectral domain or spectral embedding (Fig. 5).This coordinate system allows us to estimate the values of a scalar variable associated with the surfaces, at matched locations on the cortex.For example, each embedded vertex may have associated shape descriptors (e.g.mean curvature, sulcal depth) or regional anatomical labels (e.g.Brodmann areas) as well as spatial location.
Kernel regression (Nadaraya, 1964) is a non-parametric technique that can be used to estimate the value of one of these scalar variables y, at any given spectral location, by averaging the values of y associated with nearby vertices, using a kernel based weighting function w.Let y i denote the expected value of y at spectral location u i ′ and N(i) a set of vertices close to u i ′.Then y i is given by a weighted average of y for all vertices in N(i), y i ¼ ∑ j∈N i ð Þ y j w i; j ð Þ=∑ j∈N i ð Þ w i; j ð Þ.The vertices of each surface are evenly sampled in the spatial domain, however this is not the case in the spectral domain.Therefore, N(i) was defined as the k nearest vertices to u i ′ in order to achieve uniform spatial smoothness.
The weighting function w i; j ð Þ ¼ w u u 0 i ; u 0 j w t t j À Á =C j ð Þ was chosen which has three terms: a spatial weighting w u , a temporal weighting w t and a count function C. C j ð Þ simply counts the number of vertices in N(i) that originate from the same surface mesh as j.This removes the influence of different vertex densities for each surface at different locations within the embedding.A Gaussian kernel was used for the spatial weighting w u (u i ′, ) and the temporal weighting w t (t j ) = exp (−(t′ − t j 2 /2σ t 2 ), where t′ is a target GA.The temporal smoothness and spatial smoothness of y i are controlled by σ t and σ u respectively.σ u , can be chosen adaptively, so that the average distance of the k vertices is equal to the full width at half maximum of the kernel, . Note, the parameter k now effectively controls the spatial smoothness.Temporal smoothness, σ t , can be fixed for all time-points or adjusted within regions where the temporal sampling is lower (Serag et al., 2012).This will sacrifice temporal sharpness, when estimating y i , in order to smooth out variations in cortical anatomy by incorporating more subjects.

Surface reconstruction
For each embedded vertex i (Fig. 5), the average cortical surface position, x i , can be estimated using kernel regression, given its spectral location u i ′ and a target gestational age t′.This is achieved by treating each spatial dimension as a separate scalar field.The estimated location x i , for each vertex, gives a smooth point sampling of the average surface (Fig. 6).Reconstructing surfaces from point data is a well-studied problem driven by the need for measurement and visualisation of real world objects that can be geometrically sampled using, for example, 3D range scanning technology.We used the popular screened Poisson surface reconstruction algorithm (Kazhdan & Hoppe, 2013), to reconstruct an average surface template from x (Fig. 6).This technique also requires an estimate of the average surface normal direction, N i , for each point sample of the average surface boundary, x i .For further details on how this may be obtained using kernel regression see Appendix C.

Temporal correspondence
Establishing temporal correspondences between embeddings is a core aim of our work as it permits a longitudinal analysis of cortical growth.As a spectral analysis was performed across several overlapping time windows (Fig. 2), for each resulting embedding, a number of the contributing surfaces will also contribute to a neighbouring embedding.For the surfaces that are not shared, we can estimate the position of each of their vertices in the neighbouring embedding.This is achieved by treating the neighbouring coordinates as a labelling associated with each surface vertex within the current embedding.The neighbouring coordinates can then be estimated at unlabelled vertex locations using kernel regression (Kernel regression section), by restricting the k nearest vertices to the set of vertices already labelled.
This procedure is the equivalent of multi-atlas label propagation (Rohlfing et al., 2003;Heckemann et al., 2006), adapted for a set of embedding coordinates.This can also be iterated to propagate a set of coordinates to every surface within the cohort (Fig. 7).This ultimately   enables all surfaces to be embedded in a single coordinate system and for correspondences to be found between subjects of different ages even if these differ significantly in their cortical surface geometry (Fig. 8b).We further demonstrate the propagation of nineteen anatomical labels derived from a neonatal atlas set (Gousias et al., 2012).The MRI volumes for fetuses whose age was greater than 36 weeks GA were first parcellated by registering the neonatal atlas set to each volume and fusing the atlas labels.Each surface voxel was then assigned the anatomical label of its nearest voxel from the corresponding MR volume.This formed a surface atlas set for stepwise propagation to all other cortical surfaces (Fig. 8a).

Model parameters
Graph Laplacian edge weighting: ϕ ¼ 0:1 The parameter ϕ was introduced to control the weighting of intermesh edges relative to intra-mesh edges, when constructing the graph Laplacian for a multi-layered surface graph (ϕ = 1, gives an equal total weighting to both types of edges).The value of phi was chosen empirically to achieve a compromise between embedding regularity and sulcal alignment accuracy, determined using the validation methods discussed later in this section, and also by visual inspection of the 3D representations of the final embeddings.This value could be tuned further to increase sulcal alignment however, increasing the value of ϕ substantially can have a detrimental effect on embedding regularity.Kernel parameters: The kernel parameters were chosen by observing the smoothness of the final atlas templates.As the subjects in our particular dataset were fairly evenly distributed with age (Wright et al., 2014) a constant temporal smoothness value (σ t = 1) was chosen.Both parameters are dependant on the size of the dataset.Additionally the value chosen for k will depend on the required resolution of the surface meshes.

Sulcal alignment accuracy
To determine the precision of the sulcal alignment we manually delineated a common anatomical feature for each subject; the central sulcus.A "polyline" (a line defined by a successive sequence of vertices) was extracted by traversing each surface along the centre of the sulcus.These lines were then mapped pairwise between subjects within each embedding, using the established spectral correspondence (Fig. 9).The similarity between the central sulci after alignment was then assessed using two established measures, the discrete Fréchet distance (F d ) (Alt & Godau,   1995) and the average Fréchet distance (F a ) (Brakatsoulas et al., 2005).
The central sulcus was selected as it is a well-defined feature that can be easily delineated.It is the one of the first sulci to form, and is therefore detectable at an early GA (25 weeks) (Garel et al., 2001) and was observed in most of the subjects within the cohort (N 75%).
The Fréchet distance was chosen as it takes into account the location and ordering of the points along a curve, unlike the popular Hausdorff distance, and is therefore more sensitive to discontinuities in mapped polylines.Furthermore, the Fréchet distance has a well defined average (Brakatsoulas et al., 2005) while the Hausdorff distance is a maximum distance measure, increasing its sensitivity to outliers.The Fréchet distance is often referred to as the "dog walker's distance" as it can be described intuitively using the scenario of an owner walking its dog.It is given by the minimum leash length needed to connect a dog and its owner as they traverse two separate paths.The owner and dog may proceed at any speed and even stop, but may not backtrack.Formally, let S be a metric space, then a curve A in S is a continuous mapping from the unit interval A: [0, 1] → S. A reparameterisation α is a continuous monotonic mapping α: [0, 1] → [0, 1].If A and B are two curves in S, then their Fréchet distance F(A, B) is defined as the infimum, over all reparameterisations α and β, of the maximum distance in S between A(α(t)) and B(β(t)), for all t ∈ [0, 1].
Here d is a distance function of S.
The sulcal alignment error was compared with that given by spherical demons (Yeo et al., 2010) and two further baseline methods which register surfaces using transformation in 3D: iterative closest point (ICP) registration using affine transformations and ICP using free-form deformations (with a control point spacing of 2.5 mm) (Rueckert et al., 1999).
Alignment error was first compared within each embedding, i.e. between subjects of a similar age.The error between spectral matching and spherical demons was similar, with both outperforming the baseline methods, especially at older ages.The spectral alignment accuracy improved from 25 weeks onwards with the lowest mean error, F a ≈ 0.4 mm, at 30 weeks (Fig. 10).This trend was replicated with spherical demons and may be explained by the central sulcus becoming more distinct as it transforms from a shallow depression into a deep sulcus, which is more easily matched.After this the mean matching accuracy dropped for both methods as the cortical complexity increased, with the largest average error seen around 36 weeks (F a ≈ 1 mm).
The sulcal mapping accuracy was also compared between subjects of different ages, by first embedding both of their surfaces in the same coordinate system (Temporal correspondence section).Accuracy was then compared when mapping sulcal delineations for younger subjects (≈27 weeks) to progressively older subjects (Fig. 11a) and also when mapping sulci for older subjects (≈37 weeks) to progressively younger ones (Fig. 11b).Again performance of the two methods was similar and superior to the baseline methods.

Regional overlaps
To assess the alignment accuracy over the whole brain surface, we also evaluated the overlap between different regions of the brain.Each surface was parcellated using a publicly available neonatal MR atlas database (Gousias et al., 2012) into 6 regions (Fig. 12).MR atlases were first registered to each subject's MRI volume using a free-form deformation model (Rueckert et al., 1999), and the atlas labels were fused to create a volumetric parcellation for the subject.This parcellation was then refined using the method of Makropoulos et al. (2014).Each surface vertex was then assigned the anatomical label of its nearest voxel from the corresponding MR volume.Labels were propagated between each pair of surfaces within an embedding using the spectral (or spherical) correspondences, and their overlap measured using the Dice coefficient (Sørensen, 1948).
The results were again similar between the spectral method and spherical demons, over all six regions (Table 1).As expected, higher Dice scores were observed in the lobar regions that have a larger surface area.The lowest Dice scores were observed for the insula, which is likely due to the absence of a clearly defined inferior anatomical boundary.

Embedding regularity
The nth eigenvector of the Laplacian has the property of being smooth and monotonous over a surface between at most n poles (Lombaert et al., 2013b).This gives a diffeomorphic mapping between a surface and its spectral representation.The local structure of a surface graph is preserved in the embedding, i.e. an embedded vertex is central to its connected vertices.We conducted a simple experiment to ascertain whether the local structure of each surface graph is preserved when a group of surfaces is embedded using our methodology.
For each vertex i from a surface graph with spatial location x i and embedded position u i ′, we use edge-based smoothing (Field, 1988) to find the spectral centroid, ū i ′, of its neighbouring vertices j, i.e. ū i ′ = ∑ j ∈ N(i) w j u j ′, where w j ∝ 1/‖x i − x j ‖ and ∑ j ∈ N(i) w j = 1.If the local structure is preserved u i ′ and ū i ′ should be close.We find the closest vertex to ū i ′ in the spectral domain, with spatial position x i , and measure the distance x i −x i k k .This allows us to quantify, for each vertex, how distant its embedded position is from those of its mesh neighbours.
This experiment was repeated for spectral matching with and without regularisation of the surface links (Spectral matching section, final paragraph) and also for spherical demons replacing the spectral embedding coordinates with the method's associated spherical coordinates.For each graph the percentage of voxels where x i −x i k k¼ 0 was computed, i.e. the percentage of voxels where the local structure of the graph was preserved.The average percentage was high for spherical demons (99.57%) and both variations of spectral matching (with regularisation, 99.91%, and without, 98.83%) (Table 2).For points where x i −x i k kN0, the mean distance was also computed and this was found to be, on average, around the mesh resolution for all methods (mesh edge length ≈ 0.5 mm).The average maximum distance was   relatively high for spectral matching with irregular surface links (4.15 mm), however this decreased substantially with regularisation (0.724 mm) in line with spherical demons (0.949 mm).These results show that the local structure is largely preserved for all methods and that most irregularities are small.However, there were a small number of large irregularities with spectral matching without regularisation, particularly at older gestational ages as the cortical surfaces become more complex (Fig. 13).

Average cortical surfaces
Average cortical surfaces were generated for each week of gestation from the spectral embeddings as well as via spherical demons (Fig. 14).Generating an average surface from a spectral embedding was discussed in Surface reconstruction section, and for spherical demons we follow a similar approach.A regularly sampled spherical mesh gives the correct topology and provides the vertex connectivity for an average surface.This was readily generated by subdividing a platonic solid mesh and projecting the vertices back onto the sphere.Kernel regression (Kernel regression section) can then be used to find an average surface position at each vertex.Transforming the spherical mesh back to the spatial domain, using the average surface position at each vertex, gives an average cortical surface.
The average surfaces for both methods accurately captured the sulcal development of the individuals within the cohort, signifying a good alignment between surfaces.Both techniques produced similar average cortical surfaces, with spherical demons producing slightly sharper surfaces, suggesting a finer sulcal alignment.However, this may also be explained by the slight difference in the way surfaces are reconstructed across the methods.The timing of the sulcation observed in the average surfaces matched previous observations from MR images (Garel et al., 2001)  formed and the sylvian, calcarine, parieto-occipital and hipocampic fissures as well as the cingulate are all well defined on the average surfaces.
The first sulcus to form is the central sulcus, which is detectable from 24 weeks and very prominent after 27 weeks.By 29 weeks the precentral and postcentral sulci, frontal sulci, superior temporal sulcus and intraparietal sulcus are all clearly visible on the cortical surfaces.Visually the atlas templates are similar to those previously constructed by Habas et al. (2012) and Clouchoux et al. (2013).

Variability of average cortical surfaces
To estimate the variability of generated average surface templates, given a sample size N, we compute two distinct templates from two disjoint subsets, both of size N, of the cortical surface models within a temporal window, and then measure the average distance between them.This was performed for a range of sample sizes (2-11) and 4 target ages (24, 28, 32 and 36 weeks).For each sample size and target age, 10 iterations were performed using randomly generated disjoint pairs of subsets from the data, the average of these 10 iterations is shown as a point in Fig. 15.
At 36 weeks, where anatomical variation is the largest (for the target ages examined), the average surface generated was consistent to around 0.54 mm ± 0.04, for N = 9.The relationship between sample size and the average distance between generated average surfaces is well-modelled by a power law curve (R 2 = 0.992 ± 0.003), for sample sizes up to 11. Extrapolating this relationship to sample sizes above 11, we predict the consistency of the average surfaces generated from our complete dataset (Table 3).At 36 weeks, for our full sample size (N = 18) we therefore estimate the atlas templates generated to be consistent to around 0.40 mm, with smaller values for ages less than 36 weeks.

Conclusions
In this work, we have constructed a spatio-temporal surface atlas of the developing fetal brain.Establishing a set of preliminary correspondences between cortical surfaces allowed a group of surface graphs to be linked.This permitted an analysis of their graph Laplacian spectrum as a single entity, yielding an orthonormal basis for modes of an entire group of surfaces.This gave accurate and unbiased correspondences between surfaces, in a fraction of the time required by spherical registration methods.
The sulcal alignment was found to be significantly better than 3D deformation methods and comparable to spherical demons, a state of the art 2D surface registration technique.The regularity of the embedded surfaces was also examined, which is necessary for a diffeomorphic mapping across surfaces, and this was also found to be comparable with that provided by spherical demons.
We have constructed average cortical surface templates for each week of gestation by averaging the surface position over the spectral embeddings.These accurately capture the global shape and the timing of sulcal development over the cortex.Furthermore, the prominence of visible sulcation suggested a good alignment between embedded surfaces.We have also shown the average surfaces generated by our method are reproducible, with an average boundary variability of around 0.54 mm, for disjoint subsamples of our dataset (N = 9) for a target age of 36 weeks, where sulcal variation is high.
A spatio-temporal atlas is an important tool for analysing normal and abnormal folding patterns.This work enables a vertex-wise statistical analysis, of morphology across an embedding, which will in turn allow us to identify abnormalities associated with poor neurodevelopmental outcome in the future.
Groups of particular interest include fetal subjects with mild ventricular enlargement, one of the most common congenital anomalies, with a prevalence of around 1%.In most cases, prognosis is good (Gaglioti et al., 2005), however it is associated with a greater risk of cognitive, language, and behavioural impairments in infancy.It is thought that disturbances in cortical development may be indicative of poor outcome in these babies (Kyriakopoulou et al., 2014), and in future work we plan to investigate this hypothesis.

Fig. 1 .
Fig. 1.Cortical surface extraction.Cortical surface models for a selection of gestational ages are overlaid on their corresponding MR image volumes.The isosurface of the sub-cortical tissue probabilities is shown as a yellow contour while the cross-section of the cortical surface model is overlaid as a red contour.

Fig. 3 .
Fig. 3. Surface eigenmodes.The 1st three vibration modes of two surfaces are shown.Each mode depicts a pattern of displacement for a resonant frequency of the surface between poles, colour mapped to blue and red respectively.Figure adapted from Lombaert et al. (2013b).

Fig. 4 .
Fig. 4. Spectral matching.Vertex-wise correspondences between two surfaces (a) are given by the shortest Euclidean distance in the spectral domain (b).Colour mapping is given by the first three spectral coordinates which are mapped to RGB channels respectively.Figure adapted from Lombaert et al. (2013b).

Fig. 5 .
Fig. 5. 3D spectral embedding of a group of surfaces.Embedded vertices of 18 cortical surfaces are shown (29-31 weeks GA).Colour mapping depicts the mean curvature at each point on the original surfaces (red, convex; blue, concave).
(a) point sampling of the average cortical boundary (b) reconstructed surface

Fig. 6 .
Fig. 6.Average cortical surface template for a 30 week fetus.Kernel regression in the spectral domain is used to find the average surface position.This gives a point sampling of the average cortical surface (a).Poisson surface reconstruction is then used to extract a closed surface (b), from (a).

Fig. 7 .
Fig. 7.Stepwise embedding coordinate propagation.This diagram illustrates the propagation of a set of embedding coordinates (t′ = 37) to the cortical surfaces of younger subjects.Surface colour mapping shows the first 3 embedding coordinates encoded as RGB values.These coordinates can be viewed as a surface labelling, which may be propagated to neighbouring surfaces via surfaces that are shared across neighbouring embeddings.At each iteration, labelled subjects are treated as atlases, whose labels are propagated to any unlabelled surfaces within the same embedding using kernel regression.This process can be repeated iteratively until all surfaces are labelled, yielding a shared parameterisation for all surfaces.The dashed lines enclose the surfaces used to create each of the embeddings.

Fig. 9 .
Fig. 9. Aligned central sulci.An example of a manual delineation of the central sulcus (depicted in red), along with 38 central sulcus delineations mapped from other subjects, using the established spectral correspondence (shown in white).

Fig. 10 .
Fig.10.Intra-embedding sulcal mapping error.For each embedding, the sulcal delineations were mapped pairwise for all surfaces.The average alignment error is shown for four methods, quantified by the Fréchet distance, F d (top), and the average Fréchet distance, F a (bottom).Error bars show the standard deviation of the alignment error.

Fig. 11 .Fig. 12 .
Fig. 11.Inter-embedding sulcal mapping error.Central sulcus delineations were mapped between surfaces that contributed to different embeddings by first establishing temporal correspondences (Temporal correspondence section).The average alignment error, quantified by the average Fréchet distance, F a , is shown for four methods (a) when mapping delineations for younger source subjects (≈27 weeks) to progressively older target subjects and (b) when mapping delineations for older source subjects (≈37 weeks) to progressively younger target subjects.Error bars show the standard deviation of the alignment error.

Fig. 13 .
Fig.13.Embedding regularity.Each line signifies a vertex where the local structure of the surface was not preserved in the embedding.The length of the line is proportional to the distance that the embedded vertex is from the centroid of its neighbouring vertices, when both are mapped back to the spatial domain (the exact distance in mm is given by the colour mapping).Note the irregularities seen when the initial surface links are not regularised.

Fig. 14 .
Fig.14.Average cortical surface templates.Cortical surface templates were constructed for each week of gestation, for both spectral embeddings and spherical demons.Both methods produced visually similar templates, capturing an estimate of the average growth for the cohort.

Fig. 15 .
Fig.15.Effect of sample size on the variability of generated atlas templates.The average distance between distinct templates generated from disjoint subsets of the cortices was computed to give an estimate of the variability of generated surface templates.The points plotted show an average of 10 iterations for a particular sample size and age, with the error bars showing the standard deviation of the 10 iterations.The lines show a power law fit of the data for each target age.

Table 1
Cortical region overlap.Region labels were propagated between pairs of cortical surfaces using the established spectral or spherical correspondences, before measuring their overlap using the Dice coefficient.Average Dice scores are summarised.

Table 2
Embedding regularity.Percentage of mesh vertices where structure is preserved (see text).
. At 23 weeks, the global shape of the brain has already

Table 3
Atlas template variability prediction.By extrapolating the power law fit of the variability measured for sample sizes b11 (see Fig.15), it is possible to estimate the variability of surface templates generated using our full dataset size.