Cortical Surface Parcellation using Spherical Convolutional Neural Networks

Cortical Surface P arcellation using Spherical Con v olutional Neural Net w orks Prasanna P arv athaneni 1 ? , Sh unxing Bao 1 ? , Vish wesh Nath 1 , Neil D. W o odward 2 , Daniel O. Claassen 3 , Carissa J. Cascio 2 , Da vid H. Zald 4 , Y uank ai Huo 1 , Bennett A. Landman 1 , Ilw o o Lyu 1 1 Electrical Engineering and Computer Science, V anderbilt Universit y , TN, USA ilwoo.lyu@vanderbilt.edu 2 Psyc hiatry & Behavioral Sciences, V anderbilt Universit y Medical Cen ter, TN, USA 3 Neurology , V anderbilt Universit y Medical Center, TN, USA 4 Psyc hology , V anderbilt Univ ersity , TN, USA Abstract. W e present cortical surface parcellation using spherical deep con volutional neural netw orks. T raditional multi-atlas cortical surface parcellation requires in ter-sub ject surface registration using geometric features with high processing time on a single sub ject (2-3 hours). More- o ver, even optimal surface registration do es not necessarily pro duce op- timal cortical parcellation as parcel b oundaries are not fully matched to the geometric features. In this context, a choice of training features is imp ortan t for accurate cortical parcellation. T o utilize the net w orks efficien tly , w e propose cortical parcellation-sp ecific input data from an irregular and complicated structure of cortical surfaces. T o this end, we align ground-truth cortical parcel boundaries and use their resulting de- formation fields to generate new pairs of deformed geometric features and parcellation maps. T o extend the capabilit y of the net works, w e then smo othly morph cortical geometric features and parcellation maps using the in termediate deformation fields. W e v alidate our method on 427 adult brains for 49 labels. The exp erimen tal results show that our metho d out- p erforms traditional multi-atlas and naive spherical U-Net approac hes, while ac hieving full cortical parcellation in less than a minute. Keyw ords: cortical surface parcellation, spherical deformation, spheri- cal U-Net, surface registration 1 In tro duction Regional-based morphological analysis is a widely adapted approac h in neuro de- v elopmental studies. F or v alid regional analysis, cortical surfaces need to b e consisten tly sub divided into multi-regions based on cortical parcellation proto- cols in anatomical or functional fashion [ 5 , 10 , 16 ]. How ever, consistent lab eling of cortical regions is challenging due to the complicated cortical folds and inter- sub ject v ariabilit y . Typically , manual lab eling is tedious and time-consuming, ? P . Parv athaneni and S. Bao contributed equally to this work. 2 P . Parv athaneni et al. and there exists lab eling inconsistency ev en across exp erts. In contrast, a multi- atlas cortical parcellation approach [ 7 ] exp edites the lab eling task with algorith- mic consistency . It generally tends to provide b etter performance as the num b er of atlases increases. Unfortunately , inter-sub ject registration is unav oidable in this approac h to align multiple atlases to a target sub ject with significant com- putational demands prop ortional to the num ber of atlases. With an increasing quantit y of imaging data, conv olutional neural netw orks (CNNs) are readily av ailable to handle image segmentation problems on a struc- tured grid. Y et, traditional CNNs arc hitectures are still immature in handling non-uniform data with high complexity . This is mainly because spatial coherence incorp orated with existing deep arc hitectures is optimized on standard Euclidean image grids in addition to large memory requiremen t. In this regard, spherical CNNs recently emerge with efficient op erations on a spherical domain. Cohen et al. and Esteves et al. [ 3 , 6 ] proposed spherical CNNs arc hitectures to ac hieve computational efficiency as w ell as n umerical accuracy . Although they w ork ef- fectiv ely on classification or regression tasks, semantic segmentation tasks were not fully addressed. Later, general semantic segmentation in a spherical domain w as well discussed in [ 9 ]. Cortical surface mesh is of high complexity that still hamp ers practical use of existing CNNs due to their limited scalabilit y on large size mesh. A few recen t pioneering studies led in to dra wing the attention of CNNs to surface parcella- tion, unlik e well-dev elop ed volumetric segmentation. Cucurull et al. [ 4 ] targeted cortical parcellation on only a few ROIs due to memory capacity . Gopinath et al. [ 8 ] prop osed b etter capabilit y with their graph CNNs for full cortical parcel- lation on adult brains with comparable results to a traditional approach [ 7 ]. The equal imp ortance of training features is also emphasized in recent studies with the cen tral theme b eing the sp ecific design of the features for accurate cortical parcellation. F or example, Gopinath et al. [ 8 ] utilized spectral features for b etter cortical alignmen ts. W u et al. [ 15 ] prop osed geometry-aw are spherical features to use a standard image CNNs arc hitecture. In this pap er, we prop ose a nov el cortical parcellation approach using a deep spherical U-Net [ 9 ] that can naturally enco de relatively large surface mesh. In particular, w e focus on parcellation-specific inputs and their augmentation for efficien t utilization of the architecture and accurate parcellation results. Sp ecif- ically , w e compute deformation fields to generate deformed geometric features that b est fit ground-truth parcel b oundaries using a spherical surface registra- tion metho d [ 12 ]. Since the net works lac k generalization of input features, we further prop ose data augmentation driven by intermediate deformation fields rather than dip ole moment v ariation that o vercomes only rotational inv ariance. This can thus offer a ric h set of plausible training samples by lev eraging geo- metric features and their deformation. The k ey con tributions include (i) nov el features optimized ov er cortical parcel b oundaries and (ii) data augmentation driv en b y their in termediate deformation fields. Figure 1 shows an ov erview of the prop osed metho d. Cortical Surface Pa rcellation using Spherical Conv olutional Neural Netw orks 3 L0 L5 L5 L4 L4 L3 L3 L2 L2 L1 L1 D ow nsa m pl e Res b l o ck M e s hC onv T Res b l o ck Ge om e t r i c F e a t ur e s T r a i ni ng Da t a l =0 l =1 l =2 l =3 l =4 l =5 l =6 l =7 l =8 l =9 l =1 0 G e om e t ry Bounda ry T ra i ni ng S a m pl e s T emp late Int e rm e di a t e D e form a t i on H SD iH T emp late S phe r i c a l U - Ne t L a be l P r e di c t i on Fig. 1. An o v erview of the proposed method. Three geometric features ( iH , S D , H ) are used for training the spherical U-Net to predict 49 cortical parcellation lab els. F or eac h geometric prop ert y , intermediate deformation fields draw a total of 11+11 resp ectiv e samples after b oundary and geometric alignment for data augmentation. The cortical parcellation is then p erformed using the original geometric features of testing sub jects. 2 Metho ds 2.1 Ob jective W e denote the i th cortical lab el by z i ∈ Z + . Giv en N cortical lab els L = { z 1 , · · · , z N } and a cortical surface Ω ∈ R 3 , our ob jective is to estimate a map- ping F : Ω → L to determine a lab el for each cortical lo cation. 2.2 P arcel Boundary Alignment Deformation field. Let M : S 2 → S 2 denote a contin uous spherical deforma- tion field. Giv en x ∈ Ω and its corresp onding lo cation ˆ x , the deformation field M holds ˆ x = M ( x ) . (1) T o estimate M , we use a spherical surface registration metho d [ 12 ] that recon- structs M b y a linear combination of spherical harmonics co efficien ts; i.e., M is a function of spherical harmonics degree l . F or con venience, let M l ( · ) denote a deformation field at degree l in the remainder of the pap er. An adv antage of this metho d is to easily generate incremen tal deformation fields by adding ba- sis functions due to orthonormality of spherical harmonics bases that smo othly morph sub jects to a target template (e.g., M 0 = rigid b ody alignment, M 10 = more lo cal non-rigid deformation). Hence, once the deformation fields are com- puted with a high degree, the intermediate fields can b e reconstructed without recomputing lo w degree again. W e can then use all in termediate deformation fields for data augmen tation by adding basis functions later. 4 P . Parv athaneni et al. Boundary map. Optimal geometric alignment do es not necessarily provide op- timal cortical parcellation despite their high correlation (see precen tral gyrus in Fig. 2 for example). Also, it is important in training to hav e w ell-shaped features. Therefore, we compute deformation fields that align parcel b oundaries for more accurate prediction. T o compute such deformation fields, we need tw o steps: (1) b oundary extraction and (2) the extracted b oundaries as a contin uous function. Giv en ground-truth parcel lab els F , we can obtain b oundaries b y finding p oin ts: ∂ Ω = { s ∈ Ω | F ( s ) 6 = F ( x ) : ∀ x ∈ N ( s ) } , (2) where N ( · ) is a set of neighboring vertices on Ω . No w, we need to represent b oundaries as a contin uous function on Ω to allow deriv atives required for the ob jective function in [ 12 ]. The idea is to compute the geo desic distance b et ween ∂ Ω . Let T ( x ) : Ω → R + denote the minim um trav el-time ∂ Ω to ∀ x ∈ Ω . The tra vel-time T ( x ) holds the follo wing eikonal equation with a unit propagation sp eed: k∇ T ( x ) k = 1 , x ∈ Ω , T ( x ) = 0 , x ∈ ∂ Ω . (3) The solution is thus equiv alent to the geo desic distance from ∂ Ω . In this work, w e use the ordered upwind metho d [ 14 ]. The distance map T is of different scale for eac h region across sub jects. F or b etter surface registration, w e further normalize T with resp ect to the maximum distance per parcel, similar to the distance map normalization in [ 13 ], which provides consistent measurements across parcellation maps. Deformed data. F or input features for training, we use standard cortical ge- ometric features: mean curv ature ( iH ( x ) ∈ R ) from inflated surface (for global cortical folding agreemen t), sulcal depth ( S D ( x ) ∈ R ) and mean curv ature ( H ( x ) ∈ R ) from cortical surface (for lo cal cortical folding agreement). T o create a template, we co-register training samples in an iterative av eraging manner [ 11 ]. Here, w e compute a distance map of the mo de (most frequen t) cortical labels across the training set after their registration to the template using the three geometric features. W e then register the normalized distance map T to the tem- plate distance map at l = 10, which pro duces deformation fields M 10 . Note that w e found no noticeable improv ement of the b oundary alignment after l b ecomes greater than 10 in practice. Finally , the deformed data in training are given by a tuple P 10 ( x ) = [ iH ( M 10 ( x )) , S D ( M 10 ( x )) , H ( M 10 ( x )) , F ( M 10 ( x ))]. 2.3 Data Augmentation The prop osed feature deformation is latent. It is v alid only if an unseen surface has similar geometric patterns to the deformed features. How ever, it is unlikely to happ en unless a fairly large n um b er of training data are giv en, which suggests data augmentation to predict unseen data b etter. Thus, our goal is to generate in termediate deformed features b et w een sub jects and the target template. In Cortical Surface Pa rcellation using Spherical Conv olutional Neural Netw orks 5 Boundary P arcellation F Boundaries ∂ Ω Distance Map T Alignmen t H F eature l = 0 (rigid) l = 5 l = 10 T emplate Fig. 2. Boundary extraction and alignment. ( 1st r ow ) F or inputs for training, parcel b oundaries are obtained from ground-truth lab els (Eq. ( 2 )). The b oundaries are used to generate distance map T by solving an eik onal equation, and ( 2nd r ow ) smooth tra jectory of its deformation to a template is represented by increasing spherical har- monics degree l . ( 3r d r ow ) The features for training are accordingly deformed by the deformation fields obtained by the b oundary alignment. Note that these boundaries are quite well matched to those of the template, whereas their corresp onding deformation on mean curv ature H does not fully agree with that of the template ( yel low circles). this wa y , we can include smo oth deformation tra jectories as additional plausible training samples. Sp ecifically , we create all in termediate samples as follows: 10 [ l =0 { [ iH ( M l ( x )) , S D ( M l ( x )) , H ( M l ( x )) , F ( M l ( x ))] } . (4) T o exploit more samples, w e 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 tra jectory . 2.4 Deep Learning Arc hitecture W e adapt a state-of-the-art spherical U-net architecture designed for segmenta- tion tasks [ 9 ] that can b e naturally extended to cortical spherical parametriza- tion. In this metho d, the conv olutional kernels are predefined as differential op erators for the 1st and 2nd deriv atives, which yields fast conv olution as well as superior p erformance o v er existing spherical net works in their b enc hmarks. In our framework, three geometric features with their augmentation are provided to input c hannels and N lab els (after the deformation) to output channels. In 6 P . Parv athaneni et al. training, we incrementally reconstruct deformation fields from 0 to 10, which generates 11 × 2 times of the original size of the training set (deformation driven b y parcel b oundary and geometric feature). At the end of the testing stage, we refine predicted parcellation maps with a standard graph cut technique [ 2 ] to remo ve p oten tial isolated regions and to create smo oth parcel b oundaries. 3 Exp erimen tal Setup W e used T 1 -w eighted 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 F reeSurfer pip eline with a large n umber of vertices ( ≈ 160 k ). W e used only left hemispheres. The BrainCOLOR proto col [ 10 ] ( N = 49 ROIs) was used for lab eling with manual correction. W e 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 . W e used the cross-en tropy loss, and a total of 5,205,008 parameters w ere optimized by the Adam optimizer. The initial learning rate was set to 0 . 01 with a step deca y of 0.9 per 20 ep ochs. W e randomly divided our data into three sets: training (80%), v alidation (10%), and testing (10%). Thus, 385 × 11 × 2 = 8 , 470 training samples w ere used in our framework after data augmentation. The optimal weigh ts with the lo west v alidation loss w ere chosen up to 100 epo c hs, and eac h ep och to ok ab out 41 minutes for training of the 8,470 data. F or a fair comparison, w e applied the same graph-cut technique [ 2 ] on all the baseline metho ds. T o av oid p oten tial errors in tro duced b y misalignment, we also used the aligned features rigidly to the template, i.e., P 0 ( x ) = [ iH ( M 0 ( x )) , S D ( M 0 ( x )) , H ( M 0 ( x )) , F ( M 0 ( x ))]). 4 Results F or pro of of concept, we trained a spherical U-Net mo del [ 9 ] with the prop osed deformed features driv en by only M 10 . F rom the testing set, w e then provided the deformed geometric features P 10 driv en by their optimal b oundary alignments. The Dice o verlap was 88 . 53 ± 1 . 05%. This indicates that prediction is quite accurate if b oundary-driv en geometric features P 10 are pro vided, which is a strong assumption in practice. W e observed low Dice ov erlap of 78 . 24 ± 4 . 48% when we fed the rigid features P 0 from the same testing set to the netw orks, whic h is exp ected as the net works lack generalization. After the prop osed data augmen tation, w e observ ed Dice ov erlap of 86 . 59 ± 1 . 53% closer to that with the deformed features driv en optimal b oundary alignment. In comparison, w e p erformed surface parcellation using m ulti-atlas and spher- ical U-Net [ 9 ] with P 0 . In m ulti-atlas, we propagated lab els from all training sam- ples to a single sub ject after surface registration [ 12 ], and their final lab els were determined b y ma jority v oting. Such a large n um b er of atlases (= 385) generally results in accurate parcellation due to less bias to atlas selection with computa- tional demands (ab out a da y: registration for 3-5 minutes p er atlas). Also, the Cortical Surface Pa rcellation using Spherical Conv olutional Neural Netw orks 7 Ground-truth Multi-atlas Spherical U-Net Boundary F eature Mean Dice 82 . 73 ± 1 . 86% 85 . 23 ± 1 . 57% 86 . 59 ± 1 . 53 % Min/Max Dice 75.99%/85.54% 80.21%/88.25% 81.13% / 88.78% Fig. 3. Qualitative comparison: ground-truth, multi-atlas, spherical U-Net, and spher- ical U-Net with the proposed features. Our approach shows b etter p erformance than the other metho ds. The arrows highligh t the mismatching regions to the ground-truth. ACgG AIns AOrG AnG Calc CO Cun Ent FO FRP FuG GRe IOG ITG LiG LOrG MCgG MFC MFG MOG MOrG MPoG MPrG MSFG MTG 60 70 80 90 100 Dice (%) OCP OFuG OpIFG OrIFG PCgG PCu PHG PIns PO PoG POrG PP PrG PT SCA SFG SMC SMG SOG SPL STG TMP TrIFG TTG 60 70 80 90 100 Dice (%) Multi-atlas Spherical U-Net Boundary Feature Fig. 4. Dice o verlap of 49 regions on the left hemisphere. P aired t -tests rev eal improv ed regions with statistical significance after the FDR correction ( q = 0 . 05). 46 and 24 out of 49 regions are improv ed against multi-atlas and spherical U-Net approac hes, resp ectiv ely . The color in the x -axis lab els indicates the improv ed regions: multi-atlas ( blue ), b oth approaches ( gr e en ), and no improv ement ( black ). spherical U-Net was trained with P 0 . W e note that the spherical U-Net with P 0 is presen ted in this pap er first time for ev aluation. The Dice ov erlap w as 82 . 73 ± 1 . 86% and 85 . 23 ± 1 . 57% for m ulti-atlas and spherical U-Net approac hes, resp ectiv ely . Of these approac hes, ours ac hieved the highest Dice ov erlap with statistical significance in paired t -tests ( p < 0 . 05). Note that b oth spherical U-Net and our approac h used exactly the same input features P 0 and no deformed features were provided (i.e., no registration step in volv ed), whic h hence yields very fast cortical parcellation ( < a min ute). Figure 3 shows an example of resulting cortical parcellation maps for the three approaches. W e further p erformed paired t -tests to observe R OI-wise improv ement on individual parcels. The test statistics revealed that our approach significantly improv ed parcellation accuracy after false disco v ery rate (FDR) [ 1 ] for m ulti-comparison correction ( q = 0 . 05). Our approach outp erforms m ulti-atlas (46 regions) and spherical U-Net (24 regions) as sho wn in Fig. 4 . It is noteworth y that no regions w ere found with significantly reduced Dice ov erlap. 8 P . Parv athaneni et al. 5 Conclusion W e presen ted a cortical parcellation method using spherical U-Net with no v el features optimized ov er cortical parcellation b oundaries. T o enhance the capabil- it y of the spherical U-Net, we also incorp orated intermediate deformed features along tra jectories of the deformation fields. In the exp erimen ts, the prop osed metho d achiev ed qualitatively and quantitativ ely b etter p erformance. F urther- more, full cortical parcellation w as obtained in less than a minute. Ac knowledgmen ts. This w ork was supp orted in part by the National In- stitutes of Health under Gran ts R01EB017230, R01MH102266, R01NS097783, R01MH102272, and R01MH098098, in part by the National Science F oundation under Grant CAREER I IS 1452485, in part by the VISE/VICTR under Grant VR3029, and in part b y NVIDIA Corp oration under GPU Grant Program. References 1. Benjamini, Y., Ho c hberg, Y.: Controlling the false discov ery rate: a practical and p o w erful approac h to multiple testing. Journal of the Ro yal statistical so ciet y: series B (Metho dological) 57 (1), 289–300 (1995) 2. Bo yko v, Y., Kolmogoro v, V.: An experimental comparison of min-cut/max-flo w al- gorithms for energy minimization in vision. IEEE T ransactions on P attern Analysis & Mac hine Intelligence 26 (9), 1124–1137 (2004) 3. Cohen, T.S., Geiger, M., K¨ ohler, J., W elling, M.: Spherical cnns. arXiv preprin t arXiv:1801.10130 (2018) 4. Cucurull, G., W agst yl, K., Casanov a, A., V eli ˇ cko vi´ c, P ., Jak obsen, E., Drozdzal, M., Romero, A., Ev ans, A., Bengio, Y.: Conv olutional neural net works for mesh- based parcellation of the cerebral cortex. In: International Conference on Medical Imaging with Deep Learning (2018) 5. Desik an, R.S., S´ egonne, F., Fischl, B., Quinn, B.T., Dick erson, B.C., Black er, D., Buc kner, R.L., Dale, A.M., Maguire, R.P ., Hyman, B.T., et al.: An automated lab eling system for sub dividing the human cerebral cortex on mri scans into gyral based regions of interest. Neuroimage 31 (3), 968–980 (2006) 6. Estev es, C., Allen-Blanchette, C., Mak adia, A., Daniilidis, K.: Learning so(3) equiv- arian t representations with spherical cnns. In: Europ ean Conference on Computer Vision (2018) 7. Fisc hl, B., V an Der Kouw e, A., Destrieux, C., Halgren, E., S´ egonne, F., Salat, D.H., Busa, E., Seidman, L.J., Goldstein, J., Kennedy , D., et al.: Automatically parcellating the h uman cerebral cortex. Cerebral cortex 14 (1), 11–22 (2004) 8. Gopinath, K., Desrosiers, C., Lombaert, H.: Graph conv olutions on sp ectral em- b eddings: Learning of cortical surface data. arXiv preprin t arXiv:1803.10336 (2018) 9. Jiang, C., Huang, J., Kashinath, K., Marcus, P ., Niessner, M., et al.: Spherical cnns on unstructured grids. arXiv preprint arXiv:1901.02039 (2019) 10. Klein, A., Dal Can ton, T., Ghosh, S.S., Landman, B., Lee, J., W orth, A.: Op en lab els: online feedback for a public resource of manually lab eled brain images. In: 16th Ann ual Meeting for the Organization of Human Brain Mapping (2010) Cortical Surface Pa rcellation using Spherical Conv olutional Neural Netw orks 9 11. Lyttelton, O., Bouc her, M., Robbins, S., Ev ans, A.: An un biased iterative group registration template for cortical surface analysis. Neuroimage 34 (4), 1535–1544 (2007) 12. Lyu, I., Kang, H., W oo dw ard, N.D., St yner, M.A., Landman, B.A.: Hierarc hical spherical deformation for cortical surface registration. Medical image analysis 57 , 72–88 (2019) 13. Lyu, I., Kim, S.H., Girault, J.B., Gilmore, J.H., Styner, M.A.: A cortical shap e- adaptiv e approach to lo cal gyrification index. Medical image analysis 48 , 244–258 (2018) 14. Sethian, J.A., Vladimirsky , A.: Ordered up wind metho ds for static hamilton–jacobi equations: Theory and algorithms. SIAM Journal on Numerical Analysis 41 (1), 325–363 (2003) 15. W u, Z., Li, G., W ang, L., Shi, F., Lin, W., Gilmore, J.H., Shen, D.: Registration- free infant cortical surface parcellation using deep conv olutional neural netw orks. In: In ternational Conference on Medical Image Computing and Computer-Assisted In terven tion. pp. 672–680. Springer (2018) 16. Y eo, B.T., Krienen, F.M., Sepulcre, J., Sabuncu, M.R., Lashk ari, D., Hollinshead, M., Roffman, J.L., Smoller, J.W., Z¨ ollei, L., Polimeni, J.R., et al.: The organization of the human cerebral cortex estimated by intrinsic functional connectivity . Journal of neuroph ysiology 106 (3), 1125 (2011)