Cortical Surface Parcellation using Spherical Convolutional Neural Networks

TL;DR

This paper introduces a spherical convolutional neural network approach for cortical surface parcellation that significantly reduces processing time and improves accuracy over traditional methods by utilizing deformation-based data augmentation.

Contribution

The authors propose a novel method that uses deformation fields to generate training data for spherical CNNs, enabling fast and accurate cortical parcellation without extensive surface registration.

Findings

Outperforms traditional multi-atlas methods

Achieves full cortical parcellation in less than a minute

Validated on 427 adult brains with 49 labels

Abstract

We present cortical surface parcellation using spherical deep convolutional neural networks. Traditional multi-atlas cortical surface parcellation requires inter-subject surface registration using geometric features with high processing time on a single subject (2-3 hours). Moreover, even optimal surface registration does not necessarily produce optimal cortical parcellation as parcel boundaries are not fully matched to the geometric features. In this context, a choice of training features is important for accurate cortical parcellation. To utilize the networks efficiently, we propose cortical parcellation-specific input data from an irregular and complicated structure of cortical surfaces. To this end, we align ground-truth cortical parcel boundaries and use their resulting deformation fields to generate new pairs of deformed geometric features and parcellation maps. To extend the capability of the networks, we then smoothly morph cortical geometric features and parcellation maps using the intermediate deformation fields. We validate our method on 427 adult brains for 49 labels. The experimental results show that our method out-performs traditional multi-atlas and naive spherical U-Net approaches, while achieving full cortical parcellation in less than a minute.

Equations

$$\hat{\mathbf{x}}=M(\mathbf{x}).$$

$$\partial \Omega =\left\{\mathbf{s}\in \Omega |F(\mathbf{s})\ne F(\mathbf{x}):\forall \mathbf{x}\in \mathbf{N}(\mathbf{s})\right\},$$

$$\Vert \nabla T(\mathbf{x})\Vert =1,$$

$$T(\mathbf{x})=0,$$

$$\bigcup _{l=0}^{10}\left\{\left[iH(M_l(\mathbf{x})),SD(M_l(\mathbf{x})),H(M_l(\mathbf{x})),F(M_l(\mathbf{x}))\right]\right\}.$$

Full text

1 Introduction

Regional-based morphological analysis is a widely adapted approach in neurodevelopmental studies. For valid regional analysis, cortical surfaces need to be consistently subdivided into multi-regions based on cortical parcellation protocols in anatomical or functional fashion [5, 10, 16]. However, consistent labeling of cortical regions is challenging due to the complicated cortical folds and inter-subject variability. Typically, manual labeling is tedious and time-consuming, and there exists labeling inconsistency even across experts. In contrast, a multi-atlas cortical parcellation approach [7] expedites the labeling task with algorithmic consistency. It generally tends to provide better performance as the number of atlases increases. Unfortunately, inter-subject registration is unavoidable in this approach to align multiple atlases to a target subject with significant computational demands proportional to the number of atlases.

With an increasing quantity of imaging data, convolutional neural networks (CNNs) are readily available to handle image segmentation problems on a structured grid. Yet, traditional CNNs architectures are still immature in handling non-uniform data with high complexity. This is mainly because spatial coherence incorporated with existing deep architectures is optimized on standard Euclidean image grids in addition to large memory requirement. In this regard, spherical CNNs recently emerge with efficient operations on a spherical domain. Cohen et al. and Esteves et al. [3, 6] proposed spherical CNNs architectures to achieve computational efficiency as well as numerical accuracy. Although they work effectively on classification or regression tasks, semantic segmentation tasks were not fully addressed. Later, general semantic segmentation in a spherical domain was well discussed in [9].

Cortical surface mesh is of high complexity that still hampers practical use of existing CNNs due to their limited scalability on large size mesh. A few recent pioneering studies led into drawing the attention of CNNs to surface parcellation, unlike well-developed volumetric segmentation. Cucurull et al. [4] targeted cortical parcellation on only a few ROIs due to memory capacity. Gopinath et al. [8] proposed better capability with their graph CNNs for full cortical parcellation on adult brains with comparable results to a traditional approach [7]. The equal importance of training features is also emphasized in recent studies with the central theme being the specific design of the features for accurate cortical parcellation. For example, Gopinath et al. [8] utilized spectral features for better cortical alignments. Wu et al. [15] proposed geometry-aware spherical features to use a standard image CNNs architecture.

In this paper, we propose a novel cortical parcellation approach using a deep spherical U-Net [9] that can naturally encode relatively large surface mesh. In particular, we focus on parcellation-specific inputs and their augmentation for efficient utilization of the architecture and accurate parcellation results. Specifically, we compute deformation fields to generate deformed geometric features that best fit ground-truth parcel boundaries using a spherical surface registration method [12]. Since the networks lack generalization of input features, we further propose data augmentation driven by intermediate deformation fields rather than dipole moment variation that overcomes only rotational invariance. This can thus offer a rich set of plausible training samples by leveraging geometric features and their deformation. The key contributions include (i) novel features optimized over cortical parcel boundaries and (ii) data augmentation driven by their intermediate deformation fields. Figure 1 shows an overview of the proposed method.

2 Methods

2.1 Objective

We denote the $i$th cortical label by $z_{i}\in\mathbb{Z}^{+}$. Given $N$ cortical labels $\mathbf{L}=\{z_{1},\cdots,z_{N}\}$ and a cortical surface $\Omega\in\mathbb{R}^{3}$, our objective is to estimate a mapping $F:\Omega\rightarrow\mathbf{L}$ to determine a label for each cortical location.

2.2 Parcel Boundary Alignment

2.2.1 Deformation field.

Let $M:\mathbb{S}^{2}\rightarrow\mathbb{S}^{2}$ denote a continuous spherical deformation field. Given $\mathbf{x}\in\Omega$ and its corresponding location $\hat{\mathbf{x}}$, the deformation field $M$ holds

[TABLE]

To estimate $M$, we use a spherical surface registration method [12] that reconstructs $M$ by a linear combination of spherical harmonics coefficients; i.e., $M$ is a function of spherical harmonics degree $l$. For convenience, let $M_{l}(\cdot)$ denote a deformation field at degree $l$ in the remainder of the paper. An advantage of this method is to easily generate incremental deformation fields by adding basis functions due to orthonormality of spherical harmonics bases that smoothly morph subjects to a target template (e.g., $M_{0}$ = rigid body alignment, $M_{10}$ = more local non-rigid deformation). Hence, once the deformation fields are computed with a high degree, the intermediate fields can be reconstructed without recomputing low degree again. We can then use all intermediate deformation fields for data augmentation by adding basis functions later.

2.2.2 Boundary map.

Optimal geometric alignment does not necessarily provide optimal cortical parcellation despite their high correlation (see precentral gyrus in Fig. 2 for example). Also, it is important in training to have well-shaped features. Therefore, we compute deformation fields that align parcel boundaries for more accurate prediction. To compute such deformation fields, we need two steps: (1) boundary extraction and (2) the extracted boundaries as a continuous function. Given ground-truth parcel labels $F$, we can obtain boundaries by finding points:

[TABLE]

where $\mathbf{N}(\cdot)$ is a set of neighboring vertices on $\Omega$. Now, we need to represent boundaries as a continuous function on $\Omega$ to allow derivatives required for the objective function in [12]. The idea is to compute the geodesic distance between $\partial\Omega$. Let $T(\mathbf{x}):\Omega\rightarrow\mathbb{R}^{+}$ denote the minimum travel-time $\partial\Omega$ to $\forall\mathbf{x}\in\Omega$. The travel-time $T(\mathbf{x})$ holds the following eikonal equation with a unit propagation speed:

[TABLE]

The solution is thus equivalent to the geodesic distance from $\partial\Omega$. In this work, we use the ordered upwind method [14]. The distance map $T$ is of different scale for each region across subjects. For better surface registration, we further normalize $T$ with respect to the maximum distance per parcel, similar to the distance map normalization in [13], which provides consistent measurements across parcellation maps.

2.2.3 Deformed data.

For input features for training, we use standard cortical geometric features: mean curvature ($iH(\mathbf{x})\in\mathbb{R}$) from inflated surface (for global cortical folding agreement), sulcal depth ($SD(\mathbf{x})\in\mathbb{R}$) and mean curvature ($H(\mathbf{x})\in\mathbb{R}$) from cortical surface (for local cortical folding agreement). To create a template, we co-register training samples in an iterative averaging manner [11]. Here, we compute a distance map of the mode (most frequent) cortical labels across the training set after their registration to the template using the three geometric features. We then register the normalized distance map $T$ to the template distance map at $l=10$, which produces deformation fields $M_{10}$. Note that we found no noticeable improvement of the boundary alignment after $l$ becomes greater than 10 in practice. Finally, the deformed data in training are given by a tuple $P_{10}(\mathbf{x})=\left[iH(M_{10}(\mathbf{x})),SD(M_{10}(\mathbf{x})),H(M_{10}(\mathbf{x})),F(M_{10}(\mathbf{x}))\right]$.

2.3 Data Augmentation

The proposed feature deformation is latent. It is valid only if an unseen surface has similar geometric patterns to the deformed features. However, it is unlikely to happen unless a fairly large number of training data are given, which suggests data augmentation to predict unseen data better. Thus, our goal is to generate intermediate deformed features between subjects and the target template. In this way, we can include smooth deformation trajectories as additional plausible training samples. Specifically, we create all intermediate samples as follows:

[TABLE]

To exploit more samples, we also compute deformation fields that align the three geometric features to the template in a similar manner. Figure 2 illustrates an example of deformed features along their deformation trajectory.

2.4 Deep Learning Architecture

We adapt a state-of-the-art spherical U-net architecture designed for segmentation tasks [9] that can be naturally extended to cortical spherical parametrization. In this method, the convolutional kernels are predefined as differential operators for the 1st and 2nd derivatives, which yields fast convolution as well as superior performance over existing spherical networks in their benchmarks. In our framework, three geometric features with their augmentation are provided to input channels and $N$ labels (after the deformation) to output channels. In training, we incrementally reconstruct deformation fields from [math] to $10$, which generates $11\times 2$ times of the original size of the training set (deformation driven by parcel boundary and geometric feature). At the end of the testing stage, we refine predicted parcellation maps with a standard graph cut technique [2] to remove potential isolated regions and to create smooth parcel boundaries.

3 Experimental Setup

We used T1-weighted scans on healthy adults ($n=427$) from 23 to 34 years old, acquired from a Phillips 3T scanner. The cortical surfaces and their spherical mapping were reconstructed via a standard FreeSurfer pipeline with a large number of vertices ($\approx 160k$). We used only left hemispheres. The BrainCOLOR protocol [10] ($N=49$ ROIs) was used for labeling with manual correction.

We trained the spherical U-Net on NVIDIA Titan Xp with a batch size of 4 at level 5 of the icosahedral subdivision due to memory capacity. We used the cross-entropy loss, and a total of 5,205,008 parameters were optimized by the Adam optimizer. The initial learning rate was set to $0.01$ with a step decay of 0.9 per 20 epochs. We randomly divided our data into three sets: training (80%), validation (10%), and testing (10%). Thus, $385\times 11\times 2=8,470$ training samples were used in our framework after data augmentation. The optimal weights with the lowest validation loss were chosen up to 100 epochs, and each epoch took about 41 minutes for training of the 8,470 data. For a fair comparison, we applied the same graph-cut technique [2] on all the baseline methods. To avoid potential errors introduced by misalignment, we also used the aligned features rigidly to the template, i.e., $P_{0}(\mathbf{x})=\left[iH(M_{0}(\mathbf{x})),SD(M_{0}(\mathbf{x})),H(M_{0}(\mathbf{x})),F(M_{0}(\mathbf{x}))\right]$).

4 Results

For proof of concept, we trained a spherical U-Net model [9] with the proposed deformed features driven by only $M_{10}$. From the testing set, we then provided the deformed geometric features $P_{10}$ driven by their optimal boundary alignments. The Dice overlap was $88.53\pm 1.05\%$. This indicates that prediction is quite accurate if boundary-driven geometric features $P_{10}$ are provided, which is a strong assumption in practice. We observed low Dice overlap of $78.24\pm 4.48\%$ when we fed the rigid features $P_{0}$ from the same testing set to the networks, which is expected as the networks lack generalization. After the proposed data augmentation, we observed Dice overlap of $86.59\pm 1.53\%$ closer to that with the deformed features driven optimal boundary alignment.

In comparison, we performed surface parcellation using multi-atlas and spherical U-Net [9] with $P_{0}$. In multi-atlas, we propagated labels from all training samples to a single subject after surface registration [12], and their final labels were determined by majority voting. Such a large number of atlases ($=385$) generally results in accurate parcellation due to less bias to atlas selection with computational demands (about a day: registration for 3-5 minutes per atlas). Also, the spherical U-Net was trained with $P_{0}$. We note that the spherical U-Net with $P_{0}$ is presented in this paper first time for evaluation.

The Dice overlap was $82.73\pm 1.86\%$ and $85.23\pm 1.57\%$ for multi-atlas and spherical U-Net approaches, respectively. Of these approaches, ours achieved the highest Dice overlap with statistical significance in paired $t$-tests ($p<0.05$). Note that both spherical U-Net and our approach used exactly the same input features $P_{0}$ and no deformed features were provided (i.e., no registration step involved), which hence yields very fast cortical parcellation ($<$ a minute). Figure 3 shows an example of resulting cortical parcellation maps for the three approaches. We further performed paired $t$-tests to observe ROI-wise improvement on individual parcels. The test statistics revealed that our approach significantly improved parcellation accuracy after false discovery rate (FDR) [1] for multi-comparison correction ($q=0.05$). Our approach outperforms multi-atlas (46 regions) and spherical U-Net (24 regions) as shown in Fig. 4. It is noteworthy that no regions were found with significantly reduced Dice overlap.

5 Conclusion

We presented a cortical parcellation method using spherical U-Net with novel features optimized over cortical parcellation boundaries. To enhance the capability of the spherical U-Net, we also incorporated intermediate deformed features along trajectories of the deformation fields. In the experiments, the proposed method achieved qualitatively and quantitatively better performance. Furthermore, full cortical parcellation was obtained in less than a minute.

Acknowledgments.

This work was supported in part by the National Institutes of Health under Grants R01EB017230, R01MH102266, R01NS097783, R01MH102272, and R01MH098098, in part by the National Science Foundation under Grant CAREER IIS 1452485, in part by the VISE/VICTR under Grant VR3029, and in part by NVIDIA Corporation under GPU Grant Program.