SynthSeg: Segmentation of brain MRI scans of any contrast and resolution without retraining — TEI-tagged view

Benjamin BillotpersName, Douglas N. GrevepersName, Oula PuontipersName, Axel ThielscherpersName, Koen Van LeemputpersName, Bruce FischlpersName, Adrian V. DalcapersName, Juan Eugenio IglesiaspersName (2023date). “SynthSeg: Segmentation of brain MRI scans of any contrast and resolution without retraining.” Medical Image Analysistitle·j.

Object of study —

Research question —

Contribution —

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Abstract. Despite advances in data augmentation and transfer learning, convolutional neural networks (CNNsterm·dom) difficultly generalise to unseen domains. When segmenting brain scans, CNNsterm·dom are highly sensitive to changes in resolution and contrast: even within the same MRIterm·dom modality, performance can decrease across datasets. Here we introduce SynthSegterm·dom, the first segmentationterm·dom CNNterm·dom robust against changes in contrast and resolution. SynthSegterm·dom is trained with synthetic data sampled from a generative model conditioned on segmentationsterm·dom. Crucially, we adopt a domain randomisation strategy where we fully randomise the contrast and resolution of the synthetic training data. Consequently, SynthSegterm·dom can segment real scans from a wide range of target domains without retraining or fine-tuning, which enables straightforward analysis of huge amounts of heterogeneous clinical data. Because SynthSegterm·dom only requires segmentationsterm·dom to be trained (no images), it can learn from labels obtained by automated methods on diverse populations (e.g., ageing and diseased), thus achieving robustness to a wide range of morphological variability. We demonstrate SynthSegterm·dom on 5,000 scans of six modalities (including CT) and ten resolutions, where it exhibits unparallelled generalisation compared with supervised CNNsterm·dom, state-of-the-art domain adaptation, and Bayesian segmentationterm·dom. Finally, we demonstrate the generalisability of SynthSegterm·dom by applying it to cardiac MRIterm·dom and CT scans. 1. Introduction 1.1. Motivation Segmentationterm·dom of brain scans is of paramount importance in neuroimaging, as it enables volumetric and shape analyses ( Hynd et al., 1991 ). Although manual delineation is considered the gold standard in segmentationterm·dom, this procedure is tedious and costly, thus preventing the analysis of large datasets. Moreover, manual segmentationterm·dom of brain scans requires expertise in neuroanatomy, which, even if available, suffers from severe inter- and intra-rater variability issues ( Warfield et al., 2004 ). For these reasons, automated segmentationterm·dom methods have been proposed as a fast and reproducible alternative solution. Most recent automated segmentationterm·dom methods rely on convolutional neural networks (CNNsterm·dom) ( Ronneberger et al., 2015 ; Milletari et al., 2016 ; Kamnitsas et al., 2017b ). These are widespread in research, where the abundance of high quality scans (i.e., at high isotropic resolution and with good contrasts between tissues) enables CNNsterm·dom to obtain accurate 3D segmentationsterm·dom that can then be used in subsequent analyses such as connectivity study ( Müller et al., 2011 ). However, supervised CNNsterm·dom are far less employed in clinical settings, where physicians prefer 2D acquisitions with a sparse set of high-resolution slices, which enables faster inspection under time constraints. This leads to a huge variability in image orientation (axial, coronal, or sagittal), slice spacing, and in-plane resolution. Moreover, such 2D scans often use thick slices to increase the signal-to-noise ratio, thus introducing considerable partial voluming (PV). This effect arises when several tissue types are mixed within the same voxel, resulting in averaged intensities that are not necessarily representative of underlying tissues, often causing segmentationterm·dom methods to underperform ( Van LeemputpersName et al., 2003 ). Additionally, imaging protocols also span a huge diversity in sequences and modalities, each studying different tissue properties, and the resulting variations in intensity distributions drastically decrease the accuracy of supervised CNNsterm·dom ( Chen et al., 2019 ). Overall, the lack of tools that can cope with the large variability in MR data hinders the adoption of quantitative morphometry in the clinic. Moreover, it precludes the analysis of vast amounts of clinical scans, currently left unexplored in picture archiving and communication systems (PACS) in hospitals around the world. The ability to derive morphometric measurements from these scans would enable neuroimaging studies with sample sizes in the millions, and thus much higher statistical power than current research studies. Therefore, there is a clear need for a fast, accurate, and reproducible automated method, for segmentationterm·dom of brain scans of any contrast and resolution, and that can adapt to a wide range of populations. 1.2. Contributions In this article, we present SynthSegterm·dom, the first neural network to segment brain scans of a wide range of contrasts and resolutions, without having to be retrained or fine-tuned ( Fig. 1 ). Specifically, SynthSegterm·dom is trained with synthetic scans sampled on the fly from a generative model inspired by the Bayesian segmentationterm·dom framework, and is thus never exposed to real scans during training. Our main contribution is the adoption of a domain randomisation strategy ( Tobin et al., 2017 ), where all the parameters of the generative model (including orientation, contrast, resolution, artefacts) are fully randomised. This exposes the network to vastly different examples at each mini-batch, and thus forces it to learn domain-independent features. Moreover, we apply a random subset of common preprocessing operations to each example (e.g., skull stripping, bias field correction), such that SynthSegterm·dom can segment scans with or without preprocessing. With this domain randomisation strategy, our method only needs to be trained once. This is a considerable improvement over supervised CNNsterm·dom and domain adaptation strategies, which all need retraining or fine-tuning for each new contrast or resolution, thus hindering clinical applications. Moreover, training SynthSegterm·dom is greatly facilitated by the fact that it only requires a set of anatomical label maps to be trained (and no real images, since all training scans are synthetic). Furthermore, these maps can be obtained automatically (rather than manually), since the training scans are directly generated from their ground truths, and are thus perfectly aligned with them. This enables us to greatly improve the robustness of SynthSegterm·dom by including automated training maps from highly diverse populations. Overall, SynthSegterm·dom yields almost the accuracy of supervised CNNsterm·dom on their training domain, but unlike them, exhibits a remarkable generalisation ability. Indeed, SynthSegterm·dom consistently outperforms state-of-the-art domain adaptation strategies and Bayesian segmentationterm·dom on all tested datasets. Moreover, we demonstrate the generalisability of SynthSegterm·dom by obtaining state-of-the-art results in cross-modality cardiac segmentationterm·dom. This work extends our recent articles on contrast-adaptiveness ( BillotpersName et al., 2020a ) and PV simulation at a specific resolution ( BillotpersName et al., 2020b ; IglesiaspersName et al., 2021 ), by building, for the first time, robustness to both contrast and resolution without retraining. Our method is thoroughly evaluated in four new experiments. The code and trained model are available at https://github.com/BBillot/SynthSegterm·dom as well as in the widespread neuroimaging package FreeSurfer ( FischlpersName, 2012 ). 2. Related works Contrast-invariance in brain segmentationterm·dom has traditionally been addressed with Bayesian segmentationterm·dom. This technique is based on a generative model, which combines an anatomical prior (often a statistical atlas) and an intensity likelihood (typically a Gaussian Mixture Model, GMM). Scans are then segmented by “inverting” this model with Bayesian inference ( Wells et al., 1996 ; FischlpersName et al., 2002 ). Contrast-robustness is achieved by using an unsupervised likelihood model, with parameters estimated on each test scan ( Van LeemputpersName et al., 1999 ; Ashburner and Friston, 2005 ). However, Bayesian segmentationterm·dom requires approximately 15 min per scan ( PuontipersName et al., 2016 ), which precludes its use in time-sensitive settings. Additionally, its accuracy is limited at low resolution (LR) by PV effects ( Choi et al., 1991 ). Indeed, even if Bayesian methods can easily model PV ( Van LeemputpersName et al., 2003 ), inferring high resolution (HR) segmentationsterm·dom from LR scans quickly becomes intractable, as it requires marginalising over all possible configurations of HR labels within each LR supervoxel. While simplifications can be made ( Van LeemputpersName et al., 2003 ), PV-aware Bayesian segmentationterm·dom may still be infeasible in clinical settings. Supervised CNNsterm·dom prevail in recent medical image segmentationterm·dom ( Milletari et al., 2016 ; Kamnitsas et al., 2017b ), and are best represented by the UNet architecture ( Ronneberger et al., 2015 ). While these networks obtain fast and accurate results on their training domain, they do not generalise well to unseen contrasts ( Karani et al., 2018 ) and resolutions ( Ghafoorian et al., 2017 ), an issue known as the “domain-gap” problemterm·def ( Pan and Yang, 2010 ). Therefore, such networks need to be retrained for any new combination of contrast and resolution, often requiring new costly labelled data. This problem can partly be ameliorated by training on multi-modality scans with modality dropout ( Havaei et al., 2016 ), which results in a network able to individually segment each training modality, but that still cannot be applied to unseen domains. Data augmentation improves the robustness of CNNsterm·dom by applying simple spatial and intensity transforms to the training data ( Zhang et al., 2020 ). While such transforms often relies on handcrafted (and thus suboptimal) parameters, recent semi-supervised methods, such as adversarial augmentation, explicitly optimise the augmentation parameters during training ( Zhang et al., 2017 ; Chaitanya et al., 2019 ; Chen et al., 2022 ). Alternatively, contrastive learning methods have been proposed to leverage unsupervised data for improved generalisation ability ( Chaitanya et al., 2020 ; You et al., 2022b ). Overall, although these techniques improve generalisation in intra-modality applications ( Zhao et al., 2019 ), they generally remain insufficient in cross-modality settings ( Karani et al., 2018 ). Domain adaptation explicitly seeks to bridge a given domain gap between a source domain with labelled data, and a specific target domain without labels. A first solution is to map both domains to a common latent space, where a classifier can be trained ( Kamnitsas et al., 2017a ; Dou et al., 2019 ; Ganin et al., 2017 ; You et al., 2022a ). In comparison, generative adaptation methods seek to match the source images to the target domain with image-to-image translation methods ( Sandfort et al., 2019 ; Huo et al., 2019 ; Zhang et al., 2018 ). Since these approaches are complementary, recent methods propose to operate in both feature and image space, which leads to state-to-the-art results in cross-modality segmentationterm·dom ( Chen et al., 2019 ; Hoffman et al., 2018 ). In contrast, state-of-the-art results in intra-modality adaptation are obtained with test-time adaptation methods ( Karani et al., 2021 ; He et al., 2021 ), which rely on light fine-tuning at test-time. More generally, even though domain adaptation alleviates the need for supervision in the target domain, it still needs retraining for each new domain. Synthetic training data can be used to increase robustness by introducing surrogate domain variations, either generated with physics-based models ( Jog and FischlpersName, 2018 ), or adversarial generative networks ( Frid-Adar et al., 2018 ; Chartsias et al., 2018 ), possibly conditioned on label maps for improved semantic content ( Mahmood et al., 2020 ; Isola et al., 2017 ). These strategies enable to generate huge training datasets with perfect ground truth obtained by construction rather than human annotation ( Richter et al., 2016 ). However, although generated images may look remarkably realistic, they still suffer from a “reality gap” ( Jakobi et al., 1995 ). In addition, these methods still require retraining for every new domain, and thus do not solve the lack of generalisation of neural networks. To the best of our knowledge, no current learning method can segment medical scans of any contrast and/or resolution without retraining. Domain randomisation is a recent strategy that relies on physics-based generative models, which, unlike learning-based methods, offer full control over the generation process. Instead of handcrafting ( Jog and FischlpersName, 2018 ) or optimising ( Chen et al., 2022 ) this kind of generative model to match a specific domain, Domain Randomisation (DR) proposes to considerably enlarge the distribution of the synthetic data by fully randomising the generation parameters ( Tobin et al., 2017 ). This learning strategy is motivated by converging evidence that augmentation beyond realism leads to improved generalisation ( Bengio et al., 2011 ; Zhao et al., 2019 ). If pushed to the extreme, DR yields highly unrealistic samples, in which case real images are encompassed within the landscape of the synthetic training data ( Tremblay et al., 2018 ). As a result, this approach seeks to bridge all domain gaps in a given semantic space, rather than solving this problem for each domain gap separately. So far, DR has been used to control robotic arms ( Tobin et al., 2017 ), and for car detection in street views ( Tremblay et al., 2018 ). Here we combine DR with a generative model inspired by Bayesian segmentationterm·dom, in order to achieve, for the first time, segmentationterm·dom of brain MRIterm·dom scans of a wide range of contrasts and resolutions without retraining. 3. Methods 3.1. Generative model SynthSegterm·dom relies on a generative model from which we sample synthetic scans to train a segmentationterm·dom network ( BillotpersName et al., 2020a , b ; IglesiaspersName et al., 2021 ). Crucially, the training images are all generated on the fly with fully randomised parameters, such that the network is exposed to a different combination of contrast, resolution, morphology, artefacts, and noise at each mini-batch ( Fig. 2 ). Here we describe the generative model, which is illustrated in Fig. 3 and exemplified in Supplement 1 . 3.1.1. Label map selection and spatial augmentation The proposed generative model assumes the availability of N training label maps S n n = 1 N defined over discrete spatial coordinates ( x , y , z ) at high resolution r H R . Let all label maps take their values from a set of K labels: S n ( x , y , z ) ∈ { 1 , … , K } . We emphasise that these training label maps can be obtained manually, automatically (by segmenting brain scans with an automated method), or even can be a combination thereof – as long as they share the same labelling convention. The generative process starts by randomly selecting a segmentationterm·dom S i from the training dataset ( Fig. 3a ). In order to increase the variability of the available segmentationsterm·dom, S i is deformed with a random spatial transform ϕ , which is the composition of an affine and a non-linear transform. The affine transformation ϕ aff is the composition of three rotations θ x , θ y , θ z , three scalings s x , s y , s z , three shearings s h x , s h y , s h z , and three translations t x , t y , t z , whose parameters are sampled from uniform distributions: (1) θ x , θ y , θ z ∼ 𝒰 a r o t , b r o t , (2) s x , s y , s z ∼ 𝒰 a s c , b s c , (3) s h x , s h y , s h z ∼ 𝒰 a s h , b s h , (4) t x , t y , t z ∼ 𝒰 a t r , b t r , (5) ϕ aff = A f f θ x , θ y , θ z , s x , s y , s z , s h x , s h y , s h z , t x , t y , t z , where a r o t , b r o t , a s c , b s c , a s h , b s h , a t r , b t r are the predefined bounds of the uniform distributions, and A f f ( ⋅ ) refers to the composition of the aforementioned affine transforms. The non-linear component ϕ nonlin is a diffeomorphic transform obtained as follows. First, we sample a small vector field of size 10×10× 10 × 3 from a zero-mean Gaussian distribution of standard deviationterm·stat σ S V F drawn from 𝒰 0 , b nonlin . This field is then upsampled to full image size with trilinear interpolation to obtain a stationary velocity field (SVF). Finally, we integrate this SVF with a scale-and-square approach ( Arsigny et al., 2006 ) to yield a diffeomorphic deformation field that does not produce holes or foldings: (6) σ S V F ∼ 𝒰 0 , b nonlin , (7) SVF ′ ∼ 𝒩 10 × 10 × 10 × 3 0 , σ SVF , (8) SVF = Resample ( SVF ′ ; r H R , (9) ϕ nonlin = Integrate ( SVF ) . Finally, we obtain an augmented map L by applying ϕ to S i using nearest neighbour interpolation ( Fig. 3b ): (10) L = S i ∘ ϕ = S i ∘ ϕ a f f ∘ ϕ nonlin . 3.1.2. Initial HR synthetic image After deforming the input segmentationterm·dom, we generate an initial synthetic scan G at HR by sampling a GMM conditioned on L ( Fig. 3c ). For convenience, we regroup all the means and standard deviationsterm·stat of the GMM in M G = μ k 1 ≤ k ≤ K and Σ G = σ k 1 ≤ k ≤ K respectively. Crucially, in order to randomise the contrast of G , all the parameters in M G and Σ G are sampled at each mini-batch from uniform distributions of range a μ , b μ and a σ , b σ , respectively. We highlight that Σ G jointly models tissue heterogeneities as well as the thermal noise of the scanner. G is then formed by independently sampling at each location ( x , y , z ) the distribution indexed by L ( x , y , z ) : (11) μ k ∼ 𝒰 a μ , b μ , (12) σ k ∼ 𝒰 a σ , b σ , (13) G x , y , z ∼ 𝒩 ( μ L x , y , z , σ L x , y , z 2 ) . 3.1.3. Bias field and intensity augmentation We then simulate bias field artefacts to make SynthSegterm·dom robust to such effects. We sample a small volume of shape 4 3 from a zero-mean Gaussian distribution of random standard deviationterm·stat σ B . We then upsample this small volume to full image size, and take the voxel-wise exponential to obtain a smooth and non-negative field B . Finally, we multiply G by B to obtain a biased image G B ( Fig. 3d ), where the previous exponential ensures that division and multiplication by the same factor are equally likely ( Van LeemputpersName et al., 1999 ; Ashburner and Friston, 2005 ): (14) σ B ∼ 𝒰 0 , b B , (15) B ′ ∼ 𝒩 4 × 4 × 4 0 , σ B 2 , (16) B = U p s a m p l e B ′ , (17) G B x , y , z = G x , y , z × exp B x , y , z . Then, a final HR image I H R is produced by rescaling G B between 0 and 1, and applying a random Gamma transform (voxel-wise exponentiation) to further augment the intensity distribution of the synthetic scans. This transform enables us to skew the distribution while leaving intensities in the [0,1] interval. In practice, the exponent is sampled in the logarithmic domain from a zero-mean Gaussian distribution of standard deviationterm·stat σ γ . As a result, I H R is given by: (18) γ ∼ 𝒩 0 , σ γ 2 , (19) I H R ( x , y , z ) = ( G ( x , y , z ) − min x , y , z G max x , y , z G − min x , y , z G ) exp ( γ ) . 3.1.4. Simulation of resolution variability In order to make the network robust against changes in resolution, we now model differences in acquisition direction (i.e., axial, coronal, sagittal), slice spacing, and slice thickness. After randomly selecting a direction, the slice spacing r s p a c and slice thickness r t h i c k are respectively drawn from 𝒰 ( r H R , b r e s ) and 𝒰 r H R , r s p a c . Note that r t h i c k is bound by r s p a c as slices very rarely overlap in practice. Once all resolution parameters have been sampled, we first simulate slice thickness by blurring I H R into I σ with a Gaussian kernel that approximates the real slice excitation profile. Specifically, its standard deviationterm·stat σ t h i c k is designed to divide the power of the HR signal by 10 at the cut-off frequency ( BillotpersName et al., 2020b ). Moreover, σ t h i c k is multiplied by a random coefficient α to introduce small deviations from the nominal thickness, and to mitigate the Gaussian assumption. Slice spacing is then modelled by downsampling I σ to I L R at the prescribed low resolution r s p a c with trilinear interpolation ( Fig. 3e ) ( Van LeemputpersName et al., 2003 ). Finally, I L R is upsampled back to r H R (typically 1 mm), such that the CNNterm·dom is trained to produce crisp HR segmentationsterm·dom, regardless of the simulated resolution. This process can be summarised as: (20) r s p a c ∼ 𝒰 ( r H R , b r e s ) , (21) r t h i c k ∼ 𝒰 ( r H R , r s p a c ) , (22) α ∼ 𝒰 ( a α , b α ) , (23) σ t h i c k = 2 α log ( 10 ) ( 2 π ) − 1 r t h i c k / r H R , (24) I σ = I H R * 𝒩 ( 0 , σ t h i c k ) , (25) I L R = Resample ( I σ ; r s p a c ) , (26) I = Resample ( I L R ; r H R ) . 3.1.5. Model output and segmentationterm·dom target At each training step, our method produces two volumes: an image I sampled from the generative model, and its segmentationterm·dom target T . The latter is obtained by taking the deformed map L in ( 10 ), and resetting to background all the label values that we do not wish to segment (i.e., labels for the background structures, which are of no interest to segment). Thus, T has K ′ ≤ K labels ( Fig. 3f ). We emphasise that the central contribution of this work lies in the adopted domain randomisation strategy. The values of the hyperparameters controlling the uniform priors (listed in Supplement 2 ) are tuned using a validation set, and are the object of a sensitivity analysis in Section 5.2 . 3.2. Segmentationterm·dom network and learning Given the described generative model, a segmentationterm·dom network is trained by sampling pairs { I , T } on the fly. Here we employ a 3D UNet architecture ( Ronneberger et al., 2015 ) that we used in previous works with synthetic scans ( BillotpersName et al., 2020a ). Specifically, it consists of five levels, each separated by a batch normalisation layer ( Ioffe and Szegedy, 2015 ) along with a max-pooling (contracting path), or upsampling operation (expanding path). All levels comprise two convolution layers with 3 × 3 × 3 kernels. Every convolutional layer is associated with an Exponential Linear Unit activation ( Clevert et al., 2016 ), except for the last one, which uses a softmax. While the first layer counts 24 feature maps, this number is doubled after each max-pooling, and halved after each upsampling. Following the UNet architecture, we use skip connections across the contracting and expanding paths. Note that the network architecture is not a focus of this work: while we employ a UNet ( Ronneberger et al., 2015 ) (the most widespread network for medical images), it could in principle be replaced with any other segmentationterm·dom architecture. We use the soft Dice loss for training ( Milletari et al., 2016 ): (27) Loss ( Y , T ) = 1 − ∑ k = 1 K ′ 2 × ∑ x , y , z Y k ( x , y , z ) T k ( x , y , z ) ∑ x , y , z Y k ( x , y , z ) 2 + T k ( x , y , z ) 2 , where Y k is the soft prediction for label k ∈ 1 , … , K ′ , and T k is its associated ground truth in one-hot encoding. We use the Adam optimiser ( Kingma and Ba, 2017 ) for 300,000 steps with a learning rate of 10 −4 , and a batch size of 1. The network is trained twice, and the weights are saved every 10,000 steps. The retained model is then selected relatively to a validation set. In practice, the generative model and the segmentationterm·dom network are concatenated within a single model, which is entirely implemented on the GPU in Keras ( Chollet, 2015 ) with a Tensorflow backend ( Abadi et al., 2016 ). In total, training takes around seven days on a Nvidia Quadro RTX 6000 GPU. 3.3. Inference At test time, the input is resampled to r H R with trilinear interpolation (such that the output of the CNNterm·dom is at HR), and its intensities are rescaled between 0 and 1 with min–max normalisation (using the 1st and 99th percentiles). Preprocessed scans are then fed to the network to obtain soft predictions maps for each label. In practice, we also perform test-time augmentation ( Moshkov et al., 2020 ), which slightly improved results on the validation set. Specifically, we segment two versions of each test scan: the original one, and a right-left flipped version of it. The soft predictions of the flipped input are then flipped back to native space (while ensuring that right-left labels end up on the correct side), and averaged with the predictions of the original scan. Once test-time augmentation has been performed, final segmentationsterm·dom are obtained by keeping the biggest connected component for each label. On average, inference takes ten seconds on a Nvidia TitanXP GPU (12 GB), including preprocessing, prediction, and postprocessing. 4. General experimental setup 4.1. Brain scans and ground truths Our experiments employ eight datasets comprising 5000 scans of six different modalities and ten resolutions. The splits between training, validation, and testing are given in Table 1 . T1-39: 39 T1-weighted (T1) scans with manual labels for 30 structures ( FischlpersName et al., 2002 ). They were acquired with an MP-RAGE sequence at 1 mm isotropic resolution. HCP: 500 T1 scans of young subjects from the Human Connectome Project ( Van Essen et al., 2012 ), acquired at 0.7 mm resolution, and that we resample at 1 mm isotropic resolution. ADNI: 1500 T1 scans from the Alzheimer’s Disease Neuroimaging Initiative (ADNI. 2 ) All scans are acquired at 1 mm isotropic resolution from a wide array of scanners and protocols. In contrast to HCP, this dataset comprises ageing subjects, some diagnosed with mild cognitive impairment (MCI) or Alzheimer’s Disease (AD). As such, many subjects present strong atrophy patterns and white matter lesions. T1mix: 1000 T1 scans at 1 mm isotropic resolution from seven datasets: ABIDE ( Di Martino et al., 2014 ), ADHD200 ( The ADHD-200 Consortium, 2012 ), GSP ( Holmes et al., 2015 ), HABS ( Dagley et al., 2017 ), MCIC ( Gollub et al., 2013 ), OASIS ( Marcus et al., 2007 ), and PPMI ( Marek et al., 2011 ). We use this heterogeneous dataset to assess robustness against intra-modality contrast variations due to different acquisition protocols. FSM: 18 subjects with T1 and two other MRIterm·dom contrasts: T2-weighted (T2) and a sequence used for deep brain stimulation (DBS) ( IglesiaspersName et al., 2018 ). All scans are at 1 mm resolution. MSp: 8 subjects with T1 and proton density (PD) acquisitions at 1 mm isotropic resolution ( FischlpersName et al., 2004 ). These scans were skull stripped prior to availability, and are manually delineated for the same labels as T1-39. FLAIR: 2393 fluid-attenuated inversion recovery (FLAIR) scans at 1 × 1 × 5 mm axial resolution. These subjects are from another subset of the ADNI database, and hence also present morphological patterns related to ageing and AD. This dataset enables assessment on scans that are representative of clinical acquisitions with real-life slice selection profiles (as opposed to simulated LR, see below). Matching 1 mm T1 scans are also available, but they are not used for testing. CT: 6 computed tomography (CT) scans at 1 × 1 × 3 mm axial resolution ( West et al., 1997 ), with the aim of assessing SynthSegterm·dom on imaging modalities other than MRIterm·dom. As for the FLAIR dataset, matching 1 mm T1 scans are also available. In order to evaluate SynthSegterm·dom on more resolutions, we artificially downsample all modalities from the T1-39, FSM, and MSp datasets (all at 1 mm isotropic resolution) to nine different LR: 3, 5 and 7 mm spacing in axial, coronal, and sagittal directions. These simulations do not use real-life slice selection profiles, but are nonetheless very informative since they enable to study the segmentationterm·dom accuracy as a function of resolution. Except for T1-39 and MSp, which are available with manual labels, segmentationterm·dom ground truths are obtained by running FreeSurfer ( FischlpersName, 2012 ) on the T1 scans of each dataset, and undergo a thorough visual quality control to ensure anatomical correctness. FreeSurfer has been shown to be very robust across numerous independent T1 datasets and yields Dice scores in the range of 0.85–0.88 ( FischlpersName et al., 2002 ; Tae et al., 2008 ). Therefore, its use as silver standard enables reliable assessment of Dice below 0.85; any scores above that level are considered equally good. Crucially, using FreeSurfer segmentationsterm·dom enables us to evaluate SynthSegterm·dom on vast amounts of scans with very diverse contrasts and resolutions, which would have been infeasible with manual tracings only. 4.2. Training segmentationsterm·dom and population robustness As indicated in Table 1 , the training set for SynthSegterm·dom comprises 20 label maps from T1-39, 500 from HCP, and 500 from ADNI. Mixing these label maps considerably increases the morphological variety of the synthetic scans (far beyond the capacity of the proposed spatial augmentation alone), and thus enlarges the robustness of SynthSegterm·dom to a wide range of populations. We emphasise that using automated label maps for training is possible because synthetic images are by design perfectly aligned with their segmentationsterm·dom. We highlight that the training data does not include any real scan. Because SynthSegterm·dom requires modelling all tissue types in the images, we complement the training segmentationsterm·dom with extra-cerebral labels ( Supplement 3 ) obtained with a Bayesian segmentationterm·dom approach ( PuontipersName et al., 2020 ). Note that these new labels are dropped with 50% chances during generation, to make SynthSegterm·dom compatible with skull stripped images. Moreover, we randomly “paste” lesion labels from FreeSurfer with 50% probability, to build robustness against white matter lesions ( Supplement 4 ). Finally, we further increase the variability of the training data by randomly left/right flipping segmentationsterm·dom and cropping them to 160 3 volumes. 4.3. Competing methods We compare SynthSegterm·dom against five other approaches: T1 baseline ( Zhang et al., 2020 ): A supervised network trained on real T1 scans ( Table 1 ). This baseline seeks to assess the performance of supervised CNNsterm·dom on their source domain, as well as their generalisation to intra-modality (T1) contrast variations. For comparison purposes, we use the same UNet architecture and augmentation scheme (i.e., spatial deformation, intensity augmentation, bias field corruption) as for SynthSegterm·dom. nnUNet ( Isensee et al., 2021 ) 3 : A state-of-the-art supervised approach, very similar to the T1 baseline, except that the architecture, augmentation, pre- and postprocessing are automated with respect to the (real) T1 input data. Test-time adaptation (TTA) ( Karani et al., 2021 ) 4 : A state-of-the-art domain adaptation method relying on fine-tuning. Briefly, this strategy uses three CNNterm·dom modules: an image normaliser (five convolutional layers), a segmentationterm·dom UNet, and a denoising auto-encoder (DAE). At first, the normaliser and the UNet are jointly trained on supervised data of a source domain, while the DAE is trained separately to correct erroneous segmentationsterm·dom. At test time, the UNet and DAE are frozen, and the normaliser is fine-tuned on scans from different target domains by using the denoised predictions of the UNet as ground truth. SIFA ( Chen et al., 2019 ) 5 : A state-of-the-art unsupervised domain adaptation strategy, where image-to-image translation and segmentationterm·dom modules are jointly trained (with shared layers). SIFA seeks to align each target domain to the source data in both feature and image spaces for improved adaptation. SAMSEG ( PuontipersName et al., 2016 ): A state-of-the-art Bayesian segmentationterm·dom framework with unsupervised likelihood. As such, SAMSEG is contrast-adaptive, and can segment at any resolution, albeit not accounting for PV effects. SAMSEG does not need to be trained as it solves an optimisation problem for test each scan. Here we use the version distributed with FreeSurfer 7.0, which runs in approximately 15 min. The predictions of all methods are postprocessed as in Section 3.3 , except for nnUNet, which uses its own postprocessing. We use the default implementation for all competing methods, except for a few minor points that are listed in Supplement 5 . All learning-based methods are trained twice, and models are chosen relatively to the validation set. Segmentationsterm·dom are assessed by computing (hard) Dice scores and the 95th percentile of the surface distance (SD95, in millimetres). 5. Experiments and results Here we present four experiments that evaluates the accuracy and generalisation of SynthSegterm·dom. First, we compare it against all competing methods on every dataset. Then, we conduct an ablation study on the proposed method. The third experiment validates SynthSegterm·dom in a proof-of-concept neuroimaging group study. Finally, we demonstrate the generalisability of our method by extending it to cardiac MRIterm·dom and CT. 5.1. Robustness against contrast and resolution In this experiment, we assess the generalisation ability of SynthSegterm·dom by comparing it against all competing methods for every dataset at native resolution ( Fig. 4 , Table 2 ). Remarkably, despite SynthSegterm·dom has never been exposed to a real image during training, it reaches almost the same level of accuracy as supervised networks (T1 baseline and nnUNet) on their training domain (0.88 against 0.91 Dice scores on T1-39). Moreover, SynthSegterm·dom generalises better than supervised networks against intra-modality contrast variations, both in terms of mean (average difference of 2.5 Dice points with the T1 baseline for T1 datasets other than T1-39) and robustness (much higher lower-quartiles for SynthSegterm·dom). Crucially, the employed DR strategy yields very good generalisation, as SynthSegterm·dom sustains a remarkable accuracy across all tested contrasts and resolutions, which is infeasible with supervised networks alone. Indeed, SynthSegterm·dom outputs high-quality segmentationsterm·dom for all domains, even for FLAIR and CT scans at LR ( Fig. 6 ). Quantitatively, SynthSegterm·dom produces the best scores for all nine target domains, six of which with statistical significance for Dice and nine for SD95 ( Table 2 ). This flexibility is exemplified in Fig. 5 , where SynthSegterm·dom produces features that are almost identical for a 1 mm T1 and a 5 mm axial T2 of the same subject, the latter being effectively super-resolved to 1 mm. Although the tested domain adaptation approaches (TTA, SIFA) considerably increase the generalisation of supervised networks, they are still outperformed by SynthSegterm·dom for all contrasts and resolutions. This is a remarkable result since, as opposed to domain adaptation strategies, SynthSegterm·dom does not require any retraining. We note that fine-tuning the TTA framework makes it more robust than supervised methods for intra-modality applications (noticeably higher lower-quartiles), but its results can substantially fluctuate for larger domain gaps (e.g., on FSM-DBS, FLAIR, and CT). This is partly corrected by SIFA (improvement of 14.92 mm in SD95 for CT), which is better suited for larger domain gaps ( Karani et al., 2021 ), albeit some abrupt variations (e.g., MSp-PD). In comparison, SAMSEG yields much more constant results across MR contrasts at 1 mm resolution (average Dice score of 0.83). However, because it does not model PV, its accuracy greatly declines at low resolution: Dice scores decrease to 0.71 on the 3 mm CT dataset (5 points below SynthSegterm·dom), and to 0.64 on the 5 mm FLAIR dataset (14 points below SynthSegterm·dom). To further validate the flexibility of the proposed approach to different resolutions, we test SynthSegterm·dom on all artificially downsampled data ( Table 1 ), and we compare it against T1 baselines retrained at each resolution, as well as SAMSEG. The results show that SynthSegterm·dom maintains a very good accuracy for all tested resolutions ( Fig. 7 ). Despite the considerable loss of information at LR and heavy PV effects, SynthSegterm·dom only loses 3.8 Dice points between 1 mm and 7 mm slice spacing on average, mainly due to thin structures like the cortex ( Fig. 8 ). Meanwhile, SAMSEG is strongly affected by PV, and loses 7.6 Dice points across the same range. As before, the T1 baselines obtain remarkable results on scans similar to their training data, but generalise poorly to unseen domains (i.e., FSM-T1 and MSp-T1), where SynthSegterm·dom is clearly superior. Moreover, the gap between them on the training data progressively narrows with decreasing resolution, until it almost vanishes at 7 mm, thus making SynthSegterm·dom particularly useful for LR scans. 5.2. Ablations on DR and training label maps We now validate several aspects of our method, starting with the DR strategy. We first focus on the intensity profiles of the synthetic scans by training four variants: (i) SynthSeg-R, which is resolution-specific; (ii) SynthSeg-RC, which we retrain for every new combination of contrast and resolution by using domain-specific Gaussian priors for the GMM and resolution parameters ( BillotpersName et al., 2020b , 2021 ), (iii) a variant using slightly tighter GMM uniform priors ( μ ∈ [ 10,240 ] , σ ∈ [ 1,25 ] , instead of μ ∈ [ 0,255 ] , σ ∈ [ 0,35 ] ); and (iv) a variant with even tighter priors ( μ ∈ [ 50,200 ] , σ ∈ [ 1,15 ] ). SynthSeg-R and SynthSeg-RC assess the effect of constraining the synthetic intensity profiles to look more realistic, whereas the two last variants study the sensitivity of the chosen GMM uniform priors. Finally, we train three more networks by ablating the lesion simulation, bias field, and spatial augmentation, respectively. Fig. 9 shows that, crucially, narrowing the distributions of the generated scans in SynthSeg-R and SynthSeg-RC to simulate a specific contrast and/or resolution, leads to a consistent decrease in accuracy: despite retraining them on each target domain, they are on average lower than SynthSegterm·dom by 1.4 and 2.6 Dice points respectively. Interestingly, the variant with slightly tighter GMM priors obtains scores almost identical to the reference SynthSegterm·dom, whereas further restricting these priors (at the risk of excluding intensities encountered in real scans) leads to poorer performance (2.1 fewer Dice points on average). Finally, the bias and deformation ablations highlight the impact of those two augmentations (loss of 3.7 and 4.3 Dice points, respectively), whereas ablating the lesion simulation mainly affects the ADNI and FLAIR datasets, where the ageing subjects are more likely to present lesions (average loss of 3.9 Dice points). In a second set of experiments, we evaluate the effect of using different numbers of segmentationsterm·dom during training. Hence, we retrain SynthSegterm·dom on increasing numbers of label maps randomly selected from T1-39 ( N ∈ { 1,5 , 10,15,20 } , see Supplement 6 ). Moreover, we include the version of SynthSegterm·dom trained on all available maps, to quantify the effect of adding automated segmentationsterm·dom from diverse populations. All networks are evaluated on six representative datasets ( Fig. 10 ). The results reveal that using only one training map already attains decent Dice scores (between 0.68 and 0.80 for all datasets). As expected, the accuracy increases when adding more maps, and Dice scores plateau at N = 5 (except for MSp-PD, which levels off at N = 10 ). Interestingly, Fig. 10 also shows that SynthSegterm·dom requires fewer training examples than the T1 baseline to converge towards its maximum accuracy. Meanwhile, adding a large amount of training automated maps enables us to improve robustness to morphological variability, especially for the ADNI and FLAIR datasets with ageing and diseased subjects (Dice scores increase by 1.9 and 2.0 points respectively). To confirm this trend, we study the 3% of ADNI subjects with the largest ventricular volumes (relatively to the intracranial volume, ICV), whose morphology substantially deviates from the 20 manual training maps. For these ADNI cases (30 in total), the average Dice score increases by 4.7 Dice points for the network trained on all label maps compared with the one trained on manual maps only. This result further demonstrates the gain in robustness obtained by adding automated label maps to the training set of SynthSegterm·dom ( Fig. 11 ). 5.3. Alzheimer’s disease volumetric study In this experiment, we evaluate SynthSegterm·dom in a proof-of-concept volumetric group study, where we assess its ability to detect hippocampal atrophy related to AD ( Chupin et al., 2009 ). Specifically, we study whether SynthSegterm·dom can detect similar atrophy patterns for subjects who have been imaged with different protocols. As such, we run SynthSegterm·dom on a separate set of 100 ADNI subjects (50 controls, 50 AD), all with 1 mm isotropic T1 scans as well as FLAIR acquisitions at 5 mm axial resolution. We measure atrophy with effect sizes in predicted volumes between controls and diseased populations. Effect sizes are computed with Cohen’s d ( Cohen, 1988 ): (28) d = μ C − μ A D s , s = n C − 1 s C 2 + n A D − 1 s A D 2 n C + n A D − 2 , where μ C , s C 2 and μ A D , s A D 2 are the means and variances of the volumes for the two groups, and n C and n A D are their sizes. Hippocampal volumes are computed by summing the corresponding soft predictions, thus accounting for segmentationterm·dom uncertainties. All measured volumes are corrected for age, gender, and ICV (estimated with FreeSurfer) using a linear model. In addition to SynthSegterm·dom, we evaluate the performance of SAMSEG, and all Cohen’ d are compared to a silver standard obtained by running FreeSurfer on the T1 scans ( FischlpersName, 2012 ). The results, reported in Table 3 , reveal that both methods yield a Cohen’s d close to the ground truth for the HR T1 scans. We emphasise that, while segmenting the hippocampus in 1 mm T1 scans is of modest complexity, this task is much more difficult for 5 mm axial FLAIR scans, since the hippocampus only appears in two to three slices, and with heavy PV. As such, the accuracy of SAMSEG greatly degrades on FLAIR scans, where it obtains less than half the expected effect size. In contrast, SynthSegterm·dom sustains a high accuracy on the FLAIR scans, producing a Cohen’s d much closer to the reference value, which was obtained at HR. 5.4. Extension to cardiac segmentationterm·dom In this last experiment, we demonstrate the generalisability of SynthSegterm·dom by applying it to cardiac segmentationterm·dom. With this purpose, we employ two new datasets: MMWHS ( Zhuang et al., 2019 ), and LASC13 ( Tobon-Gomez et al., 2015 ). MMWHS includes 20 MRIterm·dom scans with in-plane resolutions from 0.78 to 1.21 mm, and slice spacings between 0.9 and 1.6 mm. MMWHS also contains 20 CT scans of non-overlapping subjects at high resolution (0.28–0.58 mm in-plane, 0.45–0.62 mm slice spacing). All these scans are available with manual labels for seven regions (see Table 4 ). On the other hand, LASC13 includes 10 MRIterm·dom heart scans at 1.25 × 1.25× 1.37 mm axial resolution, with manual labels for the left atrium only. We form the training set for SynthSegterm·dom by randomly drawing 13 label maps from MMWHS MRIterm·dom. For consistency, these training segmentationsterm·dom are all resampled at a common 1 mm isotropic resolution. Finally, the validation set consists of two more scans from MMWHS, while all the remaining scans are used for testing. Nevertheless, the training labels maps only model the target regions to segment, whereas SynthSegterm·dom requires labels for all the tissues present in the test images. Therefore, we enhance the training segmentationsterm·dom by subdividing all their labels (background and foreground) into finer subregions. This is achieved by clustering the intensities of the associated image with the Expectation Maximisation algorithm ( Dempster et al., 1977 ). First, each foreground label is divided into two regions to model blood pools. Then, the background region is split into a random number of N regions N ∈ [ 3,10 ] , which aim at representing the surrounding structures with different levels of granularity ( Supplement 6 ). All these label maps are precomputed to alleviate computing resources during training. We also emphasise that all sub-labels are merged back with their initial label for the loss computation during training. The network is trained twice as described in Section 3.2 with hyperparameters values obtained relatively to the validation set (see values in Supplement 8 ). Inference is then performed as in Section 3.3 , except for the test-time flipping augmentation that is now disabled. The results are reported in Table 4 , and show that SynthSegterm·dom segments all seven regions with very high precision (all Dice scores are above 0.8). Moreover, it maintains a very good accuracy across all tested datasets, with mean Dice scores of 0.84 and 0.88 for MMWHS MRIterm·dom and CT respectively. Interestingly, these scores are similar to the state-of-the-art results in cross-modality cardiac segmentationterm·dom obtained by Chen et al. (2019) (Dice score of 0.82), despite not being directly comparable due to differences in resolution (2 mm for Chen et al. (2019) , 1 mm for SynthSegterm·dom). Overall, segmenting all datasets at such a level of accuracy ( Fig. 12 ) is remarkable for SynthSegterm·dom, since, as opposed to Chen et al. (2019) , it is not retrained on any of them. 6. Discussion We have proposed a method for segmentationterm·dom of brain MRIterm·dom scans that is robust against changes in resolution and contrast (including CT) without retraining or fine-tuning. Our main contribution lies in the adopted domain randomisation strategy, where a segmentationterm·dom network is trained with synthetic scans of fully randomised contrast and resolution. By producing highly diverse samples that make no attempt at realism, this approach forces the network to learn domain-independent features. The impact of the DR strategy is demonstrated by the domain-constrained SynthSegterm·dom variants, for which training contrast and resolution-specific networks yields poorer performance ( Section 5.2 ). We believe this outcome is likely a combination of two phenomena. First, randomising the generation parameters enables us to mitigate the assumptions made when designing the model (e.g., Gaussian intensity distribution for each region, slice selection profile, etc.). Second, this result is consistent with converging evidence that augmenting the data beyond realism often leads to better generalisation ( Tobin et al., 2017 ; Chaitanya et al., 2019 ; Bengio et al., 2011 ). Additionally, SynthSegterm·dom enables to greatly alleviate the labelling labour for training purposes. First, it is only trained once and only requires a single set of anatomical segmentationsterm·dom (no real images), as opposed to supervised methods, which need paired images and labels for every new domain. Second, our results show that SynthSegterm·dom typically requires less training examples than supervised CNNsterm·dom to converge to its maximum performance ( Section 5.2 ). And third, parts of the training dataset can be acquired at almost no cost, by including label maps obtained by segmenting real brain scans with automated methods, and visually checking the results to ensure reasonable quality and anatomical plausibility. We highlight that while automated segmentationsterm·dom are generally not used for training (since they are prone to errors), this is made possible here by the fact that synthetic scans are, by design, perfectly aligned with their ground truths. We also emphasise that using automated segmentationsterm·dom to train SynthSegterm·dom is not only possible, but recommended, as the inclusion of such segmentationsterm·dom greatly improves robustness against highly different morphologies caused by anatomical variability (e.g., ageing subjects). Nevertheless, the employed Gaussian model imposes that the training label maps encompass tracings of all tissues present in the test scans. However, this is not a limitation in practice, since automated labels can be obtained for missing structures by simple intensity clustering. This strategy enabled us to obtain state-of-the-art results for cardiac segmentationterm·dom, where the original label maps did not describe the complex distribution of finer structures (blood pools in cardiac chambers) and surrounding tissues (vessels, bronchi, bones, etc.). Moreover, our results show that SynthSegterm·dom can handle deviations from the Gaussian model within a given structure if they are mild (like the thalamus in brain MRIterm·dom), or far away from the regions to segment (like the neck in brain MRIterm·dom). A limitation of this work is the high proportion of automated label maps used for evaluation. This choice was initially motivated by the wish to evaluate SynthSegterm·dom on a wide variety of contrasts and resolutions, which would have been infeasible with manual labels only. Nonetheless, we emphasise that a lot of testing datasets still use manual segmentationsterm·dom (T1-39, and MSp for brain segmentationterm·dom; MMWHS and LASC13 for the heart experiment), and that the remaining datasets have all undergone thorough visual quality control. Importantly, SynthSegterm·dom has shown the same remarkable generalisation ability when evaluated with manual or automated ground truths. Finally, the conclusions of this paper are further reinforced by the indirect evaluation performed in Section 5.3 , which demonstrates the accuracy and clinical utility of SynthSegterm·dom. Thanks to its unprecedented generalisation ability, SynthSegterm·dom yields direct applications in the analysis of clinical scans, for which no general segmentationterm·dom routines are available due to their highly variable acquisition procedures (sequence, resolution, hardware). Indeed, current methods deployed in the clinic include running FreeSurfer on companion 1 mm T1 scans and/or using such labels to train a supervised network (possibly with domain adaptation) to segment other sequences. However, these methods preclude the analysis of the majority of clinical datasets, where 1 mm T1 scans are rarely available. Moreover, training neural networks in the clinic is difficult in practice, since it requires corresponding expertise. In contrast, SynthSegterm·dom achieves comparable results to supervised CNNsterm·dom on their training domain (especially at LR), and can be deployed much more easily since it does not need to be retrained. 7. Conclusion In this article, we have presented SynthSegterm·dom, a learning strategy for segmentationterm·dom of brain MRIterm·dom and CT scans, where robustness against a wide range of contrasts and resolutions is achieved without any retraining or fine-tuning. First, we have demonstrated SynthSegterm·dom on 5000 scans spanning eight datasets, six modalities and 10 resolutions, where it maintains a uniform accuracy and almost attains the performance of supervised CNNsterm·dom on their training domain. SynthSegterm·dom obtains slightly better scores than state-of-the-art domain adaptation methods for small domain gaps, while considerably outperforming them for larger domain shifts. Additionally, the proposed method is consistently more accurate than Bayesian segmentationterm·dom, while being robust against PV effects and running much faster. SynthSegterm·dom can reliably be used in clinical neuroimaging studies, as it precisely detects AD atrophy patterns on HR and LR scans alike. Finally, by obtaining state-of-the-art results in cardiac cross-modality segmentationterm·dom, we have shown that SynthSegterm·dom has the potential to be applied to other medical imaging problems. While this article focuses on the use of domain randomisation to build robustness against changes in contrast and resolution, future work will seek to further improve the accuracy of the proposed method. As such, we will explore the use of adversarial networks to enhance the quality of the synthetic scans. Then, we plan to investigate the use of CNNsterm·dom to “denoise” output segmentationsterm·dom for improved robustness, and we will examine other architectures to replace the UNet employed in this work. Finally, while the ablation of the lesion simulation in Section 5.2 is a first evidence of the robustness of SynthSegterm·dom to the presence of lesions, future work will seek to precisely quantify the performance of SynthSegterm·dom when exposed to various types of lesions, tumours, and pathologies. The trained model is distributed with FreeSurfer. Relying on a single model will greatly facilitate the use of SynthSegterm·dom by researchers, since it eliminates the need for retraining, and thus the associated requirements in terms of hardware and deep learning expertise. By producing robust and reproducible segmentationsterm·dom of nearly any brain scan, SynthSegterm·dom will enable quantitative analyses of huge amounts of existing clinical data, which could greatly improve the characterisation and diagnosis of neurological disorders. Supplementary Material appendix A