Uncertainty-Aware Cortical Surface Reconstruction and Spherical Registration for Heterogeneous Brain MRI
Authors: Jesper Duemose Nielsen, Karthik Gopinath, Andrew Hoopes, Adrian Dalca, Axel Thielscher, Juan Eugenio Iglesias, Oula Puonti
Organizations: Danish Research Centre for Magnetic Resonance, Department of Radiology and Nuclear Medicine, Copenhagen University Hospital – Amager and Hvidovre, Copenhagen, Denmark · Department of Health Technology, Technical University of Denmark, Kgs. Lyngby, Denmark · Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charlestown, MA, USA · Computer Science and Artificial Intelligence Lab, Massachusetts Institute of Technology, Cambridge, MA, USA · Department of Radiology, Harvard Medical School, Boston, MA, USA · Hawkes Institute, University College London, London, UK
Surface-based cortical analysis is useful in many applications as it is a natural way to model the highly folded structure of the cerebral cortex, allowing quantities such as curvature and thickness to be easily calculated. Existing deep learning (DL) methods have focused on 1x1x1 mm3 T1-weighted (T1w) scans, thus replicating classical pipelines. Here we present a neural network that is trained exclusively on synthetic scans with random contrast and resolution. Consequently, it can be applied to almost any kind of magnetic resonance imaging (MRI) scan. In addition, our model provides vertex-wise uncertainty estimates of the predicted surface as well as a spherical registration to fsaverage space. Our model performs on par with other state-of-the-art DL methods on 1x1x1 mm3 T1w scans and outperforms the only other method capable of estimating cortical surfaces from low resolution, non-T1w scans. The predicted uncertainty, moreover, correlates with overall scan quality on low resolution scans and highlights areas of poor surface reconstruction on 1x1x1 mm3 T1w scans, suggesting that it can be used for quality control. Finally, the predicted spherical registrations prove comparable to other DL and classical registration methods. These qualities, when combined with the fast inference time, will provide the ability to analyze large, existing databases of clinical MRI scans, with the potential for finding disease-specific patterns that would otherwise be infeasible in a research setting.
Figures & tables
t1w
synth
synth
p
1mm
1mm
random
Select input image
x
x
x
T1w
x
Labeling
x
x
Sample intensities
x
x
Gamma transform
0.33
x
x
Table 1: Data augmentation steps for input scan when training the different models. FOV adjustment ensures that the data matrix is valid for processing with the CNN . p denotes the probability that an operation is applied (blank indicates that it is always applied). Indentation denotes a sub-step and x indicates that a step applies to this pipeline (if no sub-steps are marked, all sub-steps apply).
Figure 1: Examples of input data used to train the different models. The input template (blue) as well as the target WM (yellow) and GM (red) surfaces are overlaid on the input scans.
Figure 2: Overview of network training and inference. The input scan (here we show a single slice of a T1w image) is just an example; for the synth models this could be any contrast and for the synth-random model the resolution could be different. f is the image feature extractor (a UNet in our case) which extracts the features F . gi are the GCN blocks which predict the vectors used to deform the input template surface. Note that V(7) from the WM pipeline is the final WM surface prediction and serves as the input for the GM pipeline. The predicted uncertainty is used to color the cortical surfaces as they evolve; the sphere is simply colored by vertices index to show its structure (we do not color by, for example, convexity or curvature since this is not used during training). Lx refers to the losses described in table 2 (the WM - GM loss is not shown).
Task
Loss
Weight
0
20000
40000
WM
chamfer
1.0
0.0
0.0
−L
0.0
0.5
0.5
Taubin
40.0
40.0
40.0
edge
5.0
2.0
1.0
quality
2.5
1.0
0.5
Table 2: Overview of the loss weights at each stage of the training process. The first number in the “weight” columns denotes the gradient step after which those weights apply. (w) indicates that it is weighted by the uncertainty estimate.
Dataset
Subj.
Males
Diag.
Age
Scans
Train
ABIDE Di Martino et al. (2014)
1112
1001
539
7−64
807
ADHD200 Milham et al. (2012)
947
575
362
7−26
820
ADNI3 Weiner et al. (2017)
331
170
199
61−97
284
AIBL Ellis et al. (2009)
668
294
668
55−96
627
Buckner40 Fischl et al. (2002)
40
19
0
18−93
38
Table 3: Overview of training and test datasets. Datasets in the “train” section were used for training, validation, and testing (splits 0.8, 0.1, 0.1). Datasets in the “generalization” section were used to test on unseen subjects from unseen datasets (generalizability). The number in the “Diag.” (diagnosis) column refer to the number of subjects with any kind of diagnosis. “Scans” refers to the number of scans that passed quality control.
Figure 3: Illustration of the problem of correlating closest distance with uncertainty when the surface reconstruction error is large. The blue surface is reconstructed using FreeSurfer on a corresponding 1\times1\times1mm T1w scan (not shown) whereas the red and yellow surface is estimated from a 1\times1\times5mm FLAIR scans (shown). Red and yellow indicate low and high uncertainty, respectively. The three encircled parts show examples of patches where the (closest) distance to the FreeSurfer surface is small even though the surface estimate is not correct thus warranting a high uncertainty. Spherical registration aligns the folding patterns of the surfaces and by resampling one surface onto the other, we hope to achieve a more appropriate alignment of the surfaces from which to calculate the distances.
Figure 4: Comparison of different surface reconstruction methods showing which cases are handled by the different methods. The ex vivo FLASH scan is from Edlow et al. (2019) . Its original resolution is 0.1\times0.1\times0.1mm but downsampled to 1\times1\times1mm . N/A denotes that the model is not applicable to the given contrast-resolution combination.
White matter
Gray matter
WM-GM
ASSD
HD90
SIF
σ(E)
TQ
ASSD
HD90
SIF
σ(E)
TQ
IF
Test ( n=490 )
Ours
T1w-1mm
0.143
0.346
0.001
0.143
0.925
0.197
0.480
0.004
0.202
0.889
0.005
synth-1mm
0.232
0.501
0.001
0.145
0.929
0.305
0.685
0.001
0.212
0.894
0.026
synth-random
0.282
0.622
0.001
0.142
0.930
0.366
0.816
0.001
0.221
0.892
0.016
Table 4: Evaluation on 1\times1\times1mm T1w scans (except where noted) from the test and generalization sets. Metrics ASSD , HD90 , and σ(E) (standard deviation of edge lengths) are in mm; TQ (triangle quality; from eq. 16 ) is unitless; SIF and IF are in %. Please note that surfaces produced by RAC are postprocessed with FreeSurfer to, for example, remove SIF . The best value within each dataset is highlighted. n denotes the combined size of the datasets in each section.
Figure 5: Relationship between average uncertainty and surface reconstruction error for the generalization set using predictions from the T1w-1mm model. Pearson correlation cofficients are r=0.16 ( WM ) and r=0.46 ( GM ).
Figure 6: Relationship between average uncertainty and surface reconstruction error for the low resolution generalization set using predictions from the synth-random model. Pearson correlation cofficients are r=0.67 ( WM ) and r=0.69 ( GM ).
Figure 7: Correlation between uncertainty estimate and smoothed surface errors within subjects. The bar in the violin plot indicates the mean of the distribution.
Figure 8: Average uncertainty and error over subjects from the generalization set mapped to fsaverage . Surface predictions performed using the T1w-1mm model on 1\times1\times1mm T1w scans. The units are mm .
Figure 9: Average uncertainty and error over subjects from the ADNI -GO2 test set mapped to fsaverage . Surface predictions performed using the synth-random model on 1\times1\times5mm FLAIR . The units are mm .
Figure 10: Boxplots of accuracy scores across all regions of the Desikan-Killiany atlas Desikan et al. (2006) . (left) Spherical registration of surfaces predicted on 1\times1\times1mm T1w scans from the generalization set using the model T1w -1mm. (right) Spherical registration of surfaces predicted on 1\times1\times1mm T1w scans from the ADNI GO-2 dataset using the model synth-random. The reference against which the accuracy score is compute is the default spherical registration pipeline from FreeSurfer. The method denoted “map” uses coordinates obtained by spherical registration of our template to fsaverage .
Figure 11: WM surface reconstructions (FreeSurfer and our synth-random model) on the subject from ISBI2015 with the highest reconstruction error. This illustrates that FreeSurfer (blue) is sensitive to WM lesions (hypointensities) whereas the synth-random model (red-yellow where red and yellow corresponds to low and high uncertainty as estimated by the model, respectively) is less so. Additionally, we see that the estimated uncertainty highlights these difficult or problematic areas (even if the surface prediction might be correct).
Figure 12: Surface reconstructions on 1\times1\times1mm T1w scans with ( 12(a) ) no added noise, ( 12(b) ) low amount of local noise, or ( 12(c) ) high amount of local noise. Black (yellow/white) indicates low (high) values of uncertainty. Blue is the surface reconstruction from ( 12(a) ) with no overlay. Noise is sampled from a uniform distribution and added in a window of size 14\times14\times14voxels .
Appendix figures & tables8 assets
Supplementary material from the paper’s appendix.
Appendix
White matter
Gray matter
Contrast
T1w
synth
T1w
synth
Resolution
1 mm
1 mm
rand
1 mm
1 mm
rand
Model
DS
Ours
V2C
Ours
Ours
DS
Ours
V2C
Ours
Ours
ABIDE
0.153
0.150
0.211
0.247
0.291
0.210
0.239
0.328
0.353
0.409
ADHD200
0.144
0.138
0.187
0.213
0.259
0.206
0.232
0.329
0.332
0.389
ADNI3
0.197
0.179
0.198
0.273
0.344
0.171
0.186
0.200
0.293
0.367
Appendix
Table 1: Evaluation on 1\times1\times1mm T1w scans from the test set (top) and generalization set (bottom) using ASSD in mm. DS (DeepSurfer), V2C (V2C-Flow). Contrast and resolution refers to the scans on which the models were trained.
Figure 1: White matter surface. The view shows a sulcus from below.
Figure 2: Gray matter surface. The view shows a sulcus from below.
Figure 3: Relationship between average uncertainty and surface reconstruction error for the test set using predictions from the T1w-1mm model. Pearson correlation coefficients are r=0.22 ( WM ) and r=0.74 ( GM ).
Figure 4: Average uncertainty and error over subjects from the test set mapped to fsaverage . Surface predictions performed using the T1w-1mm model on 1\times1\times1mm T1w scans. The units are mm .
CAT
JOSA
ours (est)
ours (map)
Region
bankssts
0.84
0.82
0.83
0.39
caudalanteriorcingulate
0.86
0.82
0.85
0.36
caudalmiddlefrontal
0.86
0.84
0.86
0.65
cuneus
0.87
0.75
0.82
0.70
entorhinal
0.73
0.65
0.73
0.32
Appendix
Table 2: Classification accuracy per cortical region of the Desikan-Killiany atlas Desikan et al. (2006) . The results are computed using on 1\times1\times1mm T1w scans from the generalization set (AOMIC-ID1000, COBRE, ISBI2015, MCIC). We use our model trained on 1\times1\times1mm T1w scans. The reference, against which the accuracy is computed, is the default spherical registration pipeline from FreeSurfer.
CAT
JOSA
ours (est)
ours (map)
Region
bankssts
0.81
0.82
0.79
0.57
caudalanteriorcingulate
0.84
0.83
0.87
0.49
caudalmiddlefrontal
0.87
0.85
0.86
0.57
cuneus
0.87
0.75
0.82
0.70
entorhinal
0.67
0.64
0.72
0.47
Appendix
Table 3: Classification accuracy per cortical region of the Desikan-Killiany atlas Desikan et al. (2006) . The results are computed using 1\times1\times5mm FLAIR scans from the ADNI GO-2 dataset. We use our model trained on synthetic scans (random contrast, random resolution). The reference, against which the accuracy is computed, is the default spherical registration pipeline from FreeSurfer.
Figure 5: ( 5(a) ) Original white matter surface from FreeSurfer’s recon-all pipeline. ( 5(b) ) White matter surface resampled to fsaverage . Sheared triangles are clearly visible.
Reconstructing cortical WM and pial surfaces from structural magnetic resonance imaging (MRI) is a prerequisite for surface-based neuroanatomical analysis, yet remains challenging because the cortex is thin and tightly folded. Reconstruction methods can produce geometric artifacts such as mesh self-intersections and collisions between cortical surfaces, and although recent deep learning methods have reduced reconstruction time from hours to minutes, these artifacts persist. We propose SimCortex v2, a deep learning framework for simultaneous reconstruction of the left and right WM and pial surfaces from T1-weighted MRI. SimCortex v2 estimates topologically correct initial surfaces from a volumetric segmentation and refines all four jointly using multi-scale stationary velocity fields predicted by a ribbon-conditioned, U-Net-like network. We evaluated SimCortex v2 on 560 cases from 14 cohorts, thirteen of them unseen during training, spanning ages 6-89, healthy and clinical populations, and scanners from three vendors. SimCortex v2 matched the surface-distance accuracy of the strongest baseline (average symmetric surface distance 0.253 mm) while showing no detected inter-surface collision in 92.14% of cases and the lowest self-intersection fraction (0.044%) among learning-based methods, whereas every baseline produced at least one collision in every case. Source code, configuration files, pretrained weights, preprocessed data, and the exact evaluation splits are publicly released.
Deep learning (DL) has transformed neuroimaging by delivering state-of-the-art performance with reduced computation times. Yet, the numerical uncertainty inherent to DL training remains largely underexplored despite its potential to significantly impact the reliability of model outcomes. We show that training the FastSurfer segmentation model introduces substantial numerical uncertainty that exceeds its non-DL counterpart (FreeSurfer 7.3.2) in cortical regions, potentially impacting downstream clinical results. We also characterize this training-time uncertainty using random seed perturbations and demonstrate that seed-induced variability is structurally comparable to numerical variability. We then show that seed variability can be leveraged as a data augmentation technique through ensembling to improve downstream brain age regression performance. These findings position numerical uncertainty during DL training as a substantive factor in neuroimaging reliability, with measurable consequences for downstream tasks, and demonstrate that it can simultaneously be harnessed as a data augmentation technique.
Department of Computer Science Concordia University Montreal, QC H3G 1M8, Canada · Department of Computer Science Concordia University2026 Montreal, QC H3G 1M8, Canada · Krembil Centre for Neuroinformatics Centre for Addiction and Mental Health Toronto, ON M5T 0S8, Canada +1
Cortical surface registration is often driven by local geometric descriptors (e.g., sulcal depth and curvature). While this approach achieves geometric correspondence, it neglects the long-range wiring constraints imposed by white-matter anatomy. Diffusion MRI tractography offers these crucial constraints; however, prior connectivity-informed pipelines typically align precomputed connectivity matrices, making the optimization highly sensitive to connectivity estimation and its resolution. In this paper, we introduce a novel connectivity-based surface registration method that aligns cortical surfaces by operating directly on white-matter fiber-tract endpoints. We model tract endpoints as a point cloud on the product manifold Ω×Ω, where Ω represents the spherical domain of the inflated cortical hemispheres. Our alignment method iteratively (i) computes a small diffeomorphic warp for Ω by minimizing connectivity mismatch, and (ii) updates the endpoints based on this warp. The method relies on a geometric framework that ensures output warps are diffeomorphisms and has a final goal that optimizes the matching of well-known fiber bundles. Experiments on Human Connectome Project (HCP) data demonstrate improved tract-level correspondence, achieving higher connectivity-level overlap coefficients on major fiber bundles and stronger robustness across grid resolutions for Ω compared to state-of-the-art methods such as ENCORE and MSMAll.
Yang Xiang, Martin Cole, Zhengwu Zhang
Department of Statistics and Operations Research, The University of North Carolina at Chapel Hill · Department of Psychiatry, University of Rochester