Towards contrast-agnostic soft segmentation of the spinal cord.

Introduction

Spinal cord segmentation is clinically relevant, notably to compute cross-sectional area (CSA) for the diagnosis and monitoring of atrophy in multiple sclerosis (MS) ( Barkhof et al., 1997 ), spinal cord injury (SCI) ( Trolle et al., 2023 ), and in characterizing spinal cord compression ( Martin et al., 2018 ). While several approaches for semi-automatic and automatic segmentation of the spinal cord have been introduced ( De Leener et al., 2016 ; Gros et al., 2019 ; Horsfield et al., 2010 ), they all suffer from the same limitation: the output segmentation depends on the MRI contrast and acquisition parameters of the input image ( Gros et al., 2019 ). For instance, the CSA measured from T2-weighted (T2w) images is approximately 8% higher than that measured from T1-weighted (T1w) images ( Cohen-Adad et al., 2021a ; Kim et al., 2015 ). This contrast-dependency is partly due to the varying appearance of the boundary between the spinal cord and the cerebrospinal fluid because of differences in MR properties (e.g., relaxation times, spin density, flow). Different acquisition parameters and pulse sequences produce different contrast and sharpness of the spinal cord boundary, which consequently affect the output of the segmentation methods, whether they are manual, semi-automatic or automatic. The contrast-sensitive CSA also adds variability in multi-site studies, thereby reducing the sensitivity to detect subtle atrophies ( Bautin and Cohen-Adad, 2021 ).

One way to mitigate the impact of MRI contrast on metrics derived from the segmentation is to compute the CSA on various contrasts and estimate a scaling factor based on the contrast type as done by Cohen-Adad et al. (2021a) for T1w and T2w contrasts and MRI vendors. However, the scaling factors themselves are highly dependent on the MRI vendor, pulse sequence and imaging parameters, thus limiting their application to other studies.

Recent work addressed the contrast dependency of automatic segmentation in terms of model performance. Gros et al. (2019) trained separate deep learning models for each contrast. However, the ground truth (GT) masks used for training were generated using a combination of automatic PropSeg ( De Leener et al., 2014 ) and manual corrections, which were done separately for each contrast. As a result, the GT masks for each contrast were already biased, resulting in models that robustly segmented the spinal cord, but produced different CSA across contrasts (e.g., higher T2w CSA than T1w CSA). SynthSeg ( Billot et al., 2020a , b , 2023 ), a deep learning-based method primarily used for the segmentation of brain MRI scans of any contrast or resolution, leveraged GT label maps during training to synthetically generate brain images of various contrasts. While SynthSeg provides segmentations that are inherently agnostic to the input image contrast, it relies on a domain randomization strategy, where parameters such as orientation, resolution, and contrast are randomly sampled from a uniform distribution to synthetically generate the training scans. However, this requires a large number of GT segmentations, which are difficult to obtain for spinal cord scans that often include various structures such as the vertebrae, the spinal canal, intervertebral discs, nerve rootlets, surrounding muscles, the lungs and the heart.

Generalization to unseen domains is a paramount objective for deep learning algorithms. Domains can be defined as sets of images acquired from different sites and scanners, images consisting of contrasts other than those in training, or even images containing pathologies (i.e., lesions) when trained on healthy images. Domain generalization methods in the literature treat this as a domain shift problem at the fundamental level, where each contrast, for example, is seen as a different but related domain with minor differences in their marginal distributions ( Guan and Liu, 2022 ). Such methods propose to use domain adaptation techniques to transfer the differences between the source and target domains, either by mapping both domains to a shared latent space ( Kamnitsas et al., 2017 ; Dou et al., 2019b ; Ganin et al., 2016 ) or by generatively adapting the source to target domains by image-to-image translation methods ( Hoffman et al., 2018 ; Chen et al., 2019 ). Expanding further into the concept of learning domain-invariant features, regularization, as a means of creating a representative feature space consisting of various domains has also been explored ( Dou et al., 2019a ; Li et al., 2020 ; Zhang et al., 2022 ). Other related works include meta-learning for adapting model for few-shot weakly-supervised segmentation tasks ( Oliveira et al., 2022 ) and adversarial training for increasing the diversity of the training data ( Zhang et al., 2023 ). While unsupervised methods in domain adaptation alleviate the need for labeled training data, such methods still need re-training on each subsequent target domain ( Bateson et al., 2022 ), which is impractical.

Data augmentation-based domain generalization methods aim to model the domain shifts via a series of transformations applied to the input images at the source domains during training. For instance, Zhao et al. (2019) proposed a learning model for spatial transformations to synthesize additional labeled examples for one-shot segmentation in brain MRI scans. Ouyang et al. (2023) used causality-driven data augmentation specifically targeting domain shifts and acquisition shifts, while Zixian et al. (2023) proposed a location-scale augmentation using Bezier transformations, both in the context of single-source domain generalization. Zhang et al. (2020) showed that simply relying on sequential stack of data augmentation transforms based on image quality, appearance and spatial configuration, results in good generalization to unseen domains.

In addition to the contrast-dependent issues discussed above, one of the main limitations of traditional segmentation methods is that they produce binary (hard) segmentation masks, which do not account for partial volume effects ( Billot et al., 2020b ; Chaves et al., 2021 ). Partial volume effect is characterized by mixing of signals from different tissues within the same voxel, resulting in averaged intensities which are not representative of any of the underlying tissues. Binary masks do not provide calibrated output probabilities for the partial volume information of the tissue. With soft labels, the segmentation is encoded with continuous values between 0 and 1 and can therefore encode partial volume information, while resulting in better generalization ( Gros et al., 2021 ), faster learning ( Müller et al., 2019 ), and increase the precision of voxel-based morphometry or CSA measurements ( Lemay et al., 2022 ).

Contributions

In this work, we present a convolutional neural network (CNN) model for the automatic soft segmentation of spinal cord across various contrasts. Our model reduces the variability in CSA across contrasts and generalizes to spinal cord images of unseen contrasts and pathologies.

Our original contributions are as follows:

We introduce a new pipeline for generating a unique, soft GT that represents the segmentations across various MRI contrasts.

Contrary to Gros et al. (2021) where the softness was obtained implicitly after data augmentation, we propose to apply the data-augmentation transforms directly on the soft GT masks and train a contrast-agnostic SoftSeg model for spinal cord segmentation.

We show that the proposed model reduces variability morphometric measures (i.e., produces stable soft segmentations across contrasts) and shows significant improvement over prior work using contrast-specific models and domain-generalization methods.

The model is open-source and the code for pre-processing/training/inference can be found in the following GitHub release. 2 It is also integrated into the Spinal Cord Toolbox and available in v6.2 3 and higher. The rest of the paper is structured as follows: in Section 2 , we describe the training dataset, preprocessing pipeline, the training and evaluation protocols. In Section 3 , we show the results from the various validation experiments and comparisons with baselines and state-of-the-art methods. Different features and perspectives of the proposed spinal cord segmentation model are discussed in Section 4 , followed by the conclusion in Section 5 .

Materials and methods

Dataset

We used the Spine Generic Public Database 4 (Multi-Subject) ( Cohen-Adad et al., 2021a ) consisting of 267 healthy participants scanned across multiple MRI vendors (Siemens, GE and Philips) and scanner models. Each participant has a 3D T1-weighted MPRAGE (T1w) at 1 mm isotropic resolution, 3D T2w at 0.8 mm isotropic resolution, 2D T2*w axial at 0.5 × 0.5 × 3 mm (multi-echo GRE), 3D axial gradien-techo with (MT-on) and without (GRE-T1w, with a shorter TR and a higher flip angle compared to the MT-on scan) magnetization transfer pulse at 0.9 × 0.9 × 5 mm, and an axial diffusion-weighted scan motion corrected and averaged across diffusion directions at 0.9 × 0.9 × 5 mm. This multi-contrast dataset was chosen because it is publicly available, and it includes a large variety of MR contrasts that are popular in the MR community. Participants with missing contrasts or excessive artifacts were excluded from our experiments ( n = 24 out of 267).

The final dataset included 243 participants with 6 contrasts each, resulting in 1458 3D volumes in total. These were split according to 60/20/20 train/validation/test splits, resulting in 145 participants (870 volumes) for training, 49 participants (294 volumes) for validation and 49 participants (294 volumes) for testing.

Data preprocessing for ground truth generation

To eliminate the differences in CSA within the GT across contrasts, we used a unique segmentation averaged over all contrasts as the GT for training. Our objective here was to obtain the soft segmentation resulting from each contrast-specific hard segmentation. Fig. 1 shows an overview of the procedure for generating the GT using SCT ( De Leener et al., 2017 ). The GT soft segmentations are generated by averaging 6 different contrasts (T1w, T2w, T2*w, MT-on, GRE-T1w and DWI). For each participant and contrast, the spinal cord was segmented using SCT’s sct_deepseg_sc to generate a binary segmentation. Manual corrections were made when significant segmentation errors (i.e., leaking and under-segmentation) were observed in SCT’s quality control report. Since sct_deepsg_sc (DeepSeg2D) is considered the state-of-the-art for spinal cord segmentation and creating GT from scratch is highly time consuming, we obtained an initial batch of segmentations followed by manual quality control.

All images and binary segmentations were registered to the T2w image space as it has the highest resolution (0.8 mm isotropic). The registration was done using SCT’s sct_register_multimodal center-of-mass algorithm. It consists of a slice-by-slice alignment of the center of mass of the input and target segmentations (rotation and translation in x and y directions). The registration was performed in 10 iterations with a gradient step of 0.5. The segmentations of all 6 contrasts were then averaged within the T2w space to obtain a unique average soft segmentation (ranging from 0 to 1). The average of the segmentations was weighted according to the field-of-view (FoV) of each contrast. More precisely, we created a mask of the FoV of each contrast by dilating the spinal cord segmentation in the axial plane to get the complete superior–inferior coverage. Then, we registered these FoV masks to the common T2w space (see Fig. 1 , step 4), such that the contrasts with overlapping FoVs are weighted more in the unified soft GT. Then, the averaged soft segmentation was brought back to each contrast’s native space by applying the inverse warping field and using linear interpolation. This step was important because having the GT in the native space of each contrast eliminates biases due to various resolutions and FoVs during training. The original images along with the soft GT masks (both in their native space), were then used during training.

The vertebral levels were automatically labeled using SCT’s sct_ label_vertebrae command on the T1w and T2w images. Quality control was done using the sct_qc command and when necessary, labels were manually created at the posterior tip of each intervertebral disc. Since the contrasts T2*w, MT-on, GRE-T1w and DWI are axial acquisitions with thick slices, the manual or automatic labeling of the discs is not reliable. To generate the vertebral labels for those contrasts, we warped the T2w intervertebral discs to each contrast’s native space using the generated warping field in the previous step. Finally, we computed spinal cord CSA on the soft average segmentation GT and the binary segmentations averaged over the C2-C3 vertebral levels using SCT’s sct_process_segmentation .

Training protocol

In this section, we describe the (online) preprocessing, data augmentation, our proposed model and the training strategy used for the contrast agnostic segmentation of the spinal cord.

Preprocessing

All data were resampled to 1mm isotropic resolution and re-oriented to right-posterior-inferior (RPI) before training. The images and the GT labels were resampled using spline interpolation and linear interpolation, respectively. The min, median, and max shapes of all the images in the training set after resampling were 45 × 31 × 70, 192 × 230 × 106, and 230 × 320 × 320, respectively. As a consequence of having images with different orientations (3D, axial) and FoVs (cervical, cervico-thoracic, thoracic), we found center cropping to be extremely useful. Notably, the images were heavily cropped in the RL and A-P directions to keep the spinal cord at focus while the S-I direction was left uncropped. The final patch size for center-cropping was set to 64 × 192 × 320.

Data augmentation

Given the heterogeneity across contrasts, heavy data augmentation was crucial for the performance of the model. All data augmentation transforms are random, applied with a pre-defined probability and called in the following order: affine transformation with spline and linear interpolation for images and labels respectively (p = 0.9) and the rotation, scaling and translation parameters ranging between [−20,20], [−0.2,0.2], and [−0.1,0.1], respectively, elastic deformation (p = 0.5) by sampling a grid of random offsets within [25,35] and Gaussian smoothing the grid with the standard deviation (STD) between [3.5,5.5], simulation of low resolution (p = 0.25) with a downsampling and upsampling factors sampled uniformly from [0.5,1.0], gamma correction (p = 0.5) with magnitude between [0.5,3.0], where 1.0 gives the original image and smaller/larger value makes image lighter/darker, respectively, bias field adjustment (p = 0.3) with the range of random coefficients between [0.0,0.5], Gaussian noise addition (p = 0.1) with mean 0.0 and the STD spread uniformly between [0.0,0.1], Gaussian smoothing (p = 0.3) with the STD of the smoothing kernel ranging from [0.0,2.0] for all axes, intensity scaling (p = 0.15) by multiplying in the range [−0.25,1.0], random mirroring (p = 0.3) (across all axes). Lastly, all images were normalized (independently) using z -score normalization by subtracting the mean intensity and dividing by standard deviation intensity. These augmentation transforms are readily implemented in MONAI ( Jorge Cardoso et al., 2022 ).

Model architecture

Given the popularity of the nnUNet ( Isensee et al., 2021 ), we used the same architectural template found in nnUNet’s 3d_fullres model. 5 Each layer in the encoder and decoder contains two blocks, each consisting of a convolutional layer, instance normalization ( Ulyanov et al., 2016 ) and leakyReLU non-linearity ( Maas et al., 2013 ). Strided convolutions are used for downsampling while transposed convolutions are used for upsampling. Additionally, the network is trained with deep supervision ( Dou et al., 2017 ), where auxiliary losses from the feature maps at each upsampling resolution are added to the final loss. This allows for the gradients to be injected deeper into the network, thus facilitating the training of all layers. The encoder made up of 5 layers, starting with 32 feature maps at the initial layer and ending with 320 feature maps at the bottleneck (i.e. 32 → 64 → 128 → 256 → 320).

Unlike nnUNet, which uses softmax activation on the logits, we followed the SoftSeg approach ( Gros et al., 2021 ) and used normalized ReLU (NormReLU) as the final activation function. This choice is made from the observation that activation functions like sigmoid and softmax have a polarizing effect that undesirably shorten the range of soft values that carry valuable partial volume information at the boundaries. NormReLU simply normalizes the output of the ReLU activation using the maximum value, which is given by:

(1)
NormReLU ( x ) = ReLU ( x ) max ReLU ( x ) if max ReLU ( x ) ≠ 0 0 otherwise

This offers the advantage of preserving the useful properties of the ReLU activation function while ensuring that the predictions are normalized within the range of 0 and 1. A graphical representation of the model architecture is shown in Fig. 2 .

Loss function

An issue with the commonly used DiceLoss ( Milletari et al., 2016 ) is that it yields segmentation masks with sharp edges ( Deng et al., 2018 ). It has also been shown the optimizing with soft Dice leads to volumetric biases (due to under-/over-segmentation) with high inherent uncertainty ( Bertels et al., 2020 ). More specifically, for our contrast agnostic segmentation problem, DiceLoss does not drive the model towards optimizing for accurate segmentations at the boundary between the spinal cord and the cerebrospinal fluid, once it obtains a good enough segmentation of the spinal cord. This was also observed by Jia et al. (2019) in the context of left atrium segmentation in cardiac images. Therefore, the inability of DiceLoss to adapt its behavior upon reaching closer to convergence (i.e., failing to distinguish between the spinal cord/cerebrospinal fluid boundary), is the primary factor contributing towards the differences in CSA across contrasts. Although not as severe as just the DiceLoss, the same applies to other loss functions such as the Focal ( Lin et al., 2017 ), Tversky ( Salehi et al., 2017 ), and Dice-Cross Entropy (DiceCE) losses.

As suggested in Gros et al. (2021) , we considered the segmentation task as a pixel-wise regression problem (instead of classification) and trained our model with adaptive wing loss ( Wang et al., 2019 ). An immediate advantage is that a regression objective produces outputs with proper calibration while allowing soft outputs lying in [0 − 1]. Originally proposed in the context of alignment of facial landmarks via heatmap regression, there are two reasons why adaptive wing loss is a suitable candidate for obtaining soft segmentations. First , in heatmap regression, the model regresses against the GT heatmap generated by plotting Gaussian distributions centered at each facial landmark. The mode of the Gaussian (i.e. the landmark) and the pixels in its immediate neighborhood are considered as foreground, while the rest is background. In our case, this presents a similar class imbalance problem where the pixels at the spinal cord/cerebrospinal fluid boundary are outnumbered by pixels at the core of the spinal cord. Second , loss functions assigning equal weights to all pixels during training (such as DiceLoss) do not result in accurate predictions at the boundaries. Moreover, pertaining to medical image segmentation, adaptive logarithmic losses have been shown to converge faster and mitigate class imbalance ( Kaul et al., 2021 ; Gros et al., 2021 ).

The loss function is defined as follows:

(2)
AWing ( y , y ^ ) = ω ln 1 + ( y − y ^ ) ϵ ( α − y ) if | ( y − y ^ ) | < θ A | ( y − y ^ ) | − C otherwise ,

where y and y ^ correspond to GT and the predicted labels, and ω , ϵ , θ , α are the hyperparameters. Briefly, the piece-wise loss function has non-linear and linear parts. The former ensures that error between the prediction and GT smaller than θ have a larger influence in the loss function (via larger gradients during backpropagation), while the latter makes the loss function behave like the mean-squared error loss with equal weights to all voxels. The definitions of the adaptive factor A , the constant term C and the hyperparameters can be found in Section 4.2 of Wang et al. (2019) . The following hyperparameter values ω = 8.0 , ϵ = 1.0 , θ = 0.5 , α = 2.1 were set for training. We also experimented with ω = 12.0 , ϵ = 0.5 (i.e. larger ω and smaller ϵ ) as suggested in their original work, but did not observe substantial improvement in performance.

Hyperparameters & training details

We used the Adam optimizer ( Kingma and Ba, 2014 ) with a learning rate of 0.001 and a cosine annealing scheduler. The model was trained for a maximum of 200 epochs, and the batch size was set to 2. The patch size for training and sliding window inference was set to 64 × 192 × 320, same as the center cropping size. All the models were trained using the MONAI ( Jorge Cardoso et al., 2022 ) and PyTorch Lightning 6 frameworks on a single 48 GB NVIDIA A6000 GPU.

Evaluation protocol

In this section, we describe the evaluation protocol to assess the model’s performance.

Evaluation metrics

To quantitatively evaluate the segmentation accuracy, we computed the Dice coefficient (on the binarized predictions thresholded at 0.5), average surface distance (ASD), and relative volume error (RVE) for each contrast across all test participants. To assess the variability of CSA across contrasts, we computed CSA averaged over C2-C3 vertebral levels of the cervical spinal cord on all the test set predictions for each evaluated model. The following metrics were used for quantitative evaluation:

STD CSA : The standard deviation (STD) of CSA across contrasts for each participant to assess CSA variability,

Absolute CSA Error : The absolute error between the CSA of GT segmentation and the prediction for each participant.

Mathematically, the absolute CSA error ϵ is given by:

(3)
ϵ = y CSA − y ^ CSA ,

where y CSA corresponds to the CSA of the GT segmentation mask and y ^ CSA to the CSA of the prediction averaged at C2-C3 vertebral levels. The metrics are computed in mm 2 and the lower the STD and absolute CSA errors the better is the model, with the underlying assumption being that one participant should have the same spinal cord CSA across contrasts.

Baselines

We evaluated our model against 3 baselines, each with a different training strategy as described below.

Soft vs. Hard ground truth:

To assess the impact of the type of GT mask used for training and its consequential effect on CSA variability, we trained two models: one with the soft GT generated using the procedure described in Section 2.2 , and the second model with the contrast-specific hard GT generated using sct_deepseg_sc and manually corrected as required.

Single model for all contrasts vs. Contrast-specific models:

Using the soft GT masks, we trained one model per contrast (i.e. 6 models) and compared it against a single model trained on all 6 contrasts. This experiment is useful in understanding whether a single model, exposed to all contrasts during training, is capable of mitigating the CSA bias.

DiceCE loss vs. Adaptive wing loss:

As mentioned in Section 2.3.4 , optimizing only for the Dice coefficient is insufficient for accurate segmentations at the spinal cord/cerebrospinal fluid boundary. For empirically validating this hypothesis, we treated the loss function as a hyperparameter and trained two models with: (i) DiceCross entropy loss, and (ii) adaptive wing loss (the proposed model) while using both the soft and hard GT segmentations.

Among all baselines, the performance of the models were evaluated in terms of the STD CSA and the absolute CSA errors across all contrasts for each participant in the test set. Except for the ablation comparing the loss functions, all other models were trained with the adaptive wing loss.

Comparison with the state of the art (SOTA)

We also compared our model’s performance with a few SOTA methods, adapting it for spinal cord segmentation wherever necessary.

PropSeg:

PropSeg ( De Leener et al., 2014 ) is based on the iterative propagation of deformable models for spinal cord segmentation. The algorithm consists of three steps: (i) an initialization step for detecting and orientating the position of the spinal cord using a circular Hough transform, (ii) a propagation step that initializes a deformable model for its propagation along the spinal cord, and (iii), a refinement step for robust and accurate segmentation of the spinal cord.

DeepSeg:

DeepSeg ( Gros et al., 2019 ), implemented in SCT as sct_deepseg_sc , features a two-stage process: (i), the spinal cord centerline is detected using a 2D CNN with dilated convolutions, and (ii), the cord is segmented along the centerline using a 2D or 3D CNN with standard convolutions. This model was trained on ‘real world’ retrospective data from 30 sites including both healthy participants and pathological patients. Images were acquired from various vendors (Siemens, GE, Philips) and included 4 contrasts (T1w, T2w, T2*w and DWI) with a variety of image resolutions and FoVs (axial and sagittal). Because of its robustness to multi-site data, DeepSeg is an appropriate benchmark method.

nnUNet:

The nnUNet framework ( Isensee et al., 2021 ) is the SOTA in various segmentation tasks across several challenges. We used the latest version of nnUNet (i.e. nnUNetv2 ) and train both 2D and 3D variants with the default, self-configured parameters on a single fold for 1000 epochs using all contrasts together and soft GT segmentations binarized using a threshold of 0.5. This was done because nnUNet does not yet support training with soft GT labels.

SoftSeg:

SoftSeg ( Gros et al., 2021 ) showed that by skipping the binarization step after data augmentation, one can obtain the soft labels ‘for free’ and training on these soft GT results in better generalization and calibrated models. Contrary to our approach of creating soft labels by averaging the segmentations of multiple input contrasts and applying the data augmentation transforms directly on the soft labels, SoftSeg started with hard labels and trained on the soft labels obtained implicitly after data augmentation.

BigAug:

BigAug ( Zhang et al., 2020 ), is a data augmentation-based domain generalization approach, that applied a series of 9 stacked augmentation transforms based on image quality, appearance and spatial configuration to model domain shifts. While Zhang et al. (2020) reported generalization across sites/scanners only within a single contrast (T2w), we adapted their method to compare generalization across different contrasts. Specifically, BigAug was trained on a collection of all 6 contrasts with hard GT labels using Dice loss and evaluated on the basis of STD CSA, absolute CSA error and generalization to unseen contrasts.

SynthSeg:

SynthSeg ( Billot et al., 2023 ) is the SOTA method for contrast-agnostic segmentation of brain MRI scans. Starting with the original ground-truth segmentations of T1w scans, we used KMeans clustering to generate the labels for all the regions outside the spinal cord. Each T1w scan in the training set was automatically clustered into 3–10 regions and the SynthSeg model was re-trained. The output segmentations from SynthSeg were not comparable to the rest of the methods. The model performed well on T1w contrast while failing to segment the cord properly on the remaining contrasts. Hence we report its results in Section 4 of the supplementary material along with our re-implementation details.

Generalization to unseen data

As described in Section 2.1 , the Spine Generic Public Database ( Cohen-Adad et al., 2021a ) consists of healthy participants only. To evaluate our model’s ability to generalize to real world clinical data, we tested our model on three datasets of patients presenting various spinal cord pathologies, contrasts and/or on FoVs unseen during training. All patients provided written informed consent following Institutional Review Board approval and the Declaration of Helsinki.

Traumatic Spinal Cord Injury ( sci-t2w ):

This dataset consists of axial thoraco-lumbar T2w images of 80 patients with chronic traumatic spinal cord injury from the University of Colorado Anschutz Medical Campus. Acquisition was performed using MRI systems from 2 vendors (Siemens: n = 16, GE: n = 63) with 2 different field strengths (3T: n = 17, 1.5T: n = 62) and image resolutions ranging between {0.31 − 0.78} × {0.31 − 0.78} × {3 − 6} mm 3 . The challenge for the model is to be able to segment the spinal cord, in the presence of spinal cord compression, broken vertebrae and hyperintense lesions (likely edema).

Multiple Sclerosis ( ms-mp2rage ):

This dataset consists of sagittal MP2RAGE “UNI” images (1 × 1 × 1 mm resolution) of 103 healthy controls and 180 multiple sclerosis patients with visible lesions from the University of Basel acquired on a Siemens MRI scanner. The challenge for the segmentation model here is that the MP2RAGE contrast is unseen during training, and that the hypointense lesions can lead to under-segmentation mainly due to the similar signal intensity as the surrounding cerebrospinal fluid.

Cervical Radiculopathy ( radiculopathy-epi ):

This dataset consists of resting state axial gradient-echo echo-planar-imaging (GREEPI) images (0.89 × 0.89 × 5 mm resolution) of 24 participants with cervical radiculopathy and 28 age- and sex-matched healthy controls from Stanford University acquired on a GE MRI scanner. This dataset was acquired in the context of a resting state functional MRI experiment, and consists of 245 volumes that were motion corrected and averaged. Cervical radiculopathy is characterized by degenerative changes to the cervical spine, which can compress the spinal nerve roots and compromise the normal anatomy of the spinal cord, and the T2*w EPI images that can lead to strong image distortions and signal dropout making the segmentation difficult.

Lastly, we also tested our model on the FLAIR contrast. A few qualitative segmentations are shown in the supplementary material Section 2 .

Ablations on the number of contrasts

As six contrasts could be considered more than what are typically acquired in MRI examinations, we performed two more experiments ablating the number of contrasts in the preprocessing and training stages, and evaluated its downstream effect on the reduction of the CSA variability across all 6 contrasts. Starting with n = 2 contrasts (T1w and T2w), we followed the same preprocessing pipeline described in Section 2.2 and trained a model on the soft masks generated by averaging the T1w and T2w contrasts together. The same experiment was repeated for n = 4 contrasts (T1w, T2w, DWI, T2*w). The results are reported in the supplementary material Section 3 .

Results

In this section, we present the results from our proposed contrast-agnostic spinal cord segmentation model ( Section 3.1 ) and evaluate them against the baselines ( Section 3.2 ) and the existing SOTA methods ( Section 3.3 ). Then, we show the generalization capabilities of our model on unseen, out-of-distribution data ( Section 3.4 ). Lastly, in Section 3.5 , we compare the CPU inference times between various methods.

In all the plots in the following sections, the proposed model is denoted by soft_all , meaning that the model was trained with a soft GT averaged from the individual segmentations of each of the 6 contrasts and adaptive wing loss was used as the loss function.

Contrast-agnostic spinal cord segmentation

Table 1 shows the quantitative results for the proposed contrast-agnostic spinal cord segmentation model soft_all . For each contrast, we present the mean ± standard deviation across test participants for Dice coefficients, relative volume errors (in %), and average surface distances. While the Dice coefficients are consistent across all contrasts, we note a slight under-segmentation in the case of MT-on and DWI contrasts (reflected by the negative RVE) and an over-segmentation for the T1w, T2w, T2*w and GRE-T1w contrasts.

Fig. 3 shows the violin plot with absolute CSA error between the predictions and the GT across 6 contrasts (the lower the better). The mean CSA error is less than 2 mm 2 across all contrasts, which is encouraging given that 2 mm 2 represents only two pixels at an axial resolution of 1 × 1 mm.

Fig. 4 shows the comparison of CSA across contrasts between two models trained with soft GT (top panel) and hard GT (bottom panel). Training with soft GT resulted in similar CSA across all contrasts, while the hard GT training resulted in drastically uneven distributions of CSA. Note that training with contrast-specific hard GT masks with adaptive wing loss is precisely the training strategy used in SoftSeg ( Gros et al., 2021 ), hence we denote this model as hard_all_SoftSeg . We conducted a one-way paired ANOVA on CSA across contrasts for both soft_all and hard_all_SoftSeg models. The MRI contrast had a significant effect on CSA for both methods ( p < 0.05). A followup posthoc analysis (two-sided Bonferroni-corrected non-parametric Wilcoxon signed-rank test) revealed that for soft_all, T2w/T2*w, T2*w/T1w, T1w/GRE-T1w, T2w/GRE-T1w, T2w/MT-on, T2w/DWI, T1w/MT-on, and T1w/DWI pairs of contrast showed a significant difference between CSA ( p < 0.05). While for hard_all_SoftSeg , all pairs of contrasts show significantly different CSA values ( p < 0.05) except for the T2w/MT-on pair.

Despite significant paired differences across contrasts for both soft _all and hard_all_SoftSeg models, the variability of CSA across contrasts did indeed reduce significantly, as described in the next section.

Comparison with baselines

Fig. 5 compares the performance of the soft_all model with the baselines in terms of the STD of CSA across contrasts. The STD is computed across all 6 contrasts for each participant and each test participant is represented by an individual point in the scatter plot. Starting with the GT violin plots on the left (gray), we observe that the root cause of CSA variability can be mitigated using soft average segmentations as the GT during training. This is also supported by the fact that the hard_all_SoftSeg model, trained on binary GT, results in higher STD when compared to its soft counterparts. Within the models trained using soft GT, the performance of a single model trained on all contrasts ( soft_all ) is similar to the average of 6 models trained individually on each contrast ( soft_per_contrast ). Lastly, on comparing the effect of training with DiceCE loss ( soft_all_diceCE_loss ) and adaptive wing loss ( soft_all ), we observe significantly lower CSA STD across contrasts when using the regression-based adaptive wing loss (p < 0.001).

While the previous figure showed the STD of CSA across contrasts for each of the test participant, Fig. 6 compares the absolute CSA error between the prediction and the GT for our model and the baselines. The points in the scatter plots represent each test image and the mean CSA error is given on top of the violin plots. The superior performance of soft_all suggests that a combination of soft segmentations GT along with adaptive wing loss is crucial for mitigating CSA variability. When comparing the CSA errors between soft_per_contrast and soft_all models, we observe significantly lower CSA errors (p < 0.001) with the latter as also depicted in the violin plot containing a high density of scatter points between 0 − 2 mm 2 range. A few more plots comparing the variability in absolute CSA errors across each contrast between the baselines are shown in Section 1.2 of the supplementary material .

Comparison with the state of the art

Given the results from the comparison with baselines in the previous section, we considered soft_all as the best model for subsequent comparisons with the existing SOTA methods. Recall that soft_all denotes the model that has been trained with soft GT averaged from the segmentations of each of the 6 contrasts and adaptive wing loss was used as the loss function. Table 2 shows the quantitative results for the proposed segmentation model ( soft_all ), and the existing SOTA methods. The soft_all model outperforms the SOTA methods in terms on Dice coefficient and performs slightly worse in terms of RVE compared to DeepSeg2D and average surface distance (ASD) compared to nnUNet2D, respectively. DeepSeg2D and both nnUNet models show under-segmentation (reflected by the negative RVE value) while PropSeg, DeepSeg3D, BigAug and SoftSeg models present oversegmentation. Also note the relatively low Dice STD for the soft_all model, suggesting higher robustness (i.e., other models fail more often in presence of difficult images).

Fig. 7 compares the STD of CSA across contrasts for our best model and the existing methods. DeepSeg2D has the highest STD across contrasts, followed by hard_all_SoftSeg and hard_all_BigAug .

Interestingly, despite the hard requirement of having binarized GT and training with DiceCE loss, both nnUNet models achieved similar STD across contrasts with the 3D model showing lower STD across contrasts. Overall, both nnUNet 2D and 3D models showed higher STD when compared to our best model soft_all . In Supplementary Figure S8 , we show that for nnUNet3D, the relatively lower STD values do not necessarily correspond to lower absolute CSA errors across individual contrasts (as our model achieves in Fig. 3 ).

Fig. 8 shows a comparison between our best model ( soft_all ) and the other methods in terms of the absolute CSA error. Similar to the trend observed in Fig. 7 , we noted that soft_all achieves the lowest mean absolute CSA error with 1.64 ± 1.42 mm 2 across all test images. A closer look at the absolute CSA errors per contrast for each method is shown in Section 1.3 of the supplementary material . In supplementary Figures S9 and S10 , we show the comparison with all the methods including DeepSeg 3D and PropSeg to ovoid overcrowding Figs. 7 and 8 .

To understand the impact of data augmentation and loss functions on training with soft masks, we compared our models trained using DiceCE and adaptive wing loss with nnUNet in Fig. 9 . Keeping the loss function fixed (i.e. DiceCE), we can observe that data augmentation transforms in nnUNet3D play a stronger role as it achieves lower STD and absolute CSA errors compared to the soft_all_diceCE_loss model. However, as shown by soft_all , switching to the adaptive wing loss irrespective of data augmentation transforms leads to further reduction in the CSA variability across contrasts.

In Fig. 10 , we plotted the level of agreement between the CSA estimated by our model ( soft_all ) and the SOTA methods on T1w and T2w contrasts. The models trained with individual hard masks (namely, deepseg2D, hard_all_SoftSeg , and hard_all_BigAug ) showed large discrepancies between the estimated CSA. Interestingly, the models trained with soft masks obtained by averaging all contrasts moved closer to the diagonal line representing perfect alignment between T1w/T2w CSA. Within soft_all and nnUNet3D , our model achieves better alignment between the two contrasts, thus confirming that it reduces discrepancies between these two popular contrasts in spinal cord imaging. The correlation plots for the remaining pairs of contrasts for soft_all are shown in Figure S1 in the supplementary material .

Generalization to unseen data

Fig. 11 shows the predictions for 8 representative patients with spinal cord injury from the sci-t2w dataset. Despite the presence of spinal cord lesions, we notice that soft_all and nnUNet models were able to correctly segment the spinal cord, while DeepSeg2D T2w and soft_per_contrast T2w models under-segmented the spinal cord (except one over-segmentation pointed by the red arrow). The hard_all_SoftSeg, hard_all_BigAug, hard_all_diceCE and soft_all_diceCE models showed under-segmentation of the hyperintense lesions in the spinal cord as indicated by the red arrows.

Fig. 12 shows the predictions for 6 representative patients with multiple sclerosis and 2 healthy participants from the ms-mp2rage dataset. Despite the presence of spinal cord MS lesions, we observe that the soft_all model was able to correctly segment the spinal cord, while DeepSeg2D T1w model under-segmented the spinal cord typically at the lesion location. The soft_per_contrast T1w and nnUNet models were unable to properly segment the spinal cord in the presence of hypointense lesions. hard_all_BigAug showed some under-segmentation of the location of the hypointense lesions. hard_all_SoftSeg, hard_all_diceCE and soft_all_ diceCE showed similar performance and were able to properly segment the spinal cord.

Fig. 13 shows the predictions for 4 representative patients with cervical radiculopathy and 4 healthy controls from the radiculopathy -epi dataset. The soft_all model was able to correctly segment the spinal cord even with the poor image quality of the gradient echo EPI. In contrast, the soft_per_contrast T2*w was unable to segment the spinal cord in almost all cases. The DeepSeg2D T2*w and nnUNet performed slightly worst than soft_all , notably in slices affected by signal drop out (e.g., HC013), and nnUNet, hard_all_diceCE , hard_all_SoftSeg had a tendency to leak into the cerebrospinal fluid (red arrows). The hard_all_BigAug model had trouble with getting the shape of the spinal cord, mainly in the presence of signal drop-out, shown by the red arrow. Representative examples of generalization to the FLAIR contrast can be found in the supplementary material Figure S11 .

Table 3 presents the quantitative metrics for the models’ performance on all three external datasets. We do not report the Dice coefficients for DeepSeg2D since the GT masks were generated using this model with manual corrections on the slices presenting over/under-segmentations, hence biasing the scores. We observed that the nnUNet3D performed well on the T2w and GRE-EPI contrasts while performing poorly on the MP2RAGE “UNI” contrast (Dice = 0.24 ± 0.25) due to large under-segmentation overall (RVE = −82.51 ± 19.27%) while both models using DiceCE ( hard_all_DiceCE and soft_all_ DiceCE ) performed well on MP2RAGE “UNI” contrast, as we can also observe in Fig. 12 . Both BigAug ( Zhang et al., 2020 ) and SoftSeg ( Gros et al., 2021 ) models outperformed nnUNet3D only for the MP2RAGE contrast. The soft_per_contrast model performed poorly on GRE-EPI (Dice = 0.29 ± 0.18) but performed well on the MP2RAGE “UNI” and T2w data. In terms of RVE, we see that nnUNet 3D consistently shows under-segmentation on all datasets. The soft_all model outperforms all the tested methods for radiculopathyepi datasets in terms of Dice coefficients and performs similar to nnUNet on the sci-t2w dataset, and similar to SoftSeg on the msmp2rage dataset.

Inference time

Fig. 14 compares the average inference time across participants for the DeepSeg2D, nnUNet3D and soft_all methods. The inference is done on a CPU (Intel Xeon E7–4850 @ 2.10 GHz) and the time is shown in seconds. soft_all takes up to 2 min per prediction on average irrespective of the contrast, whereas nnUNet 3D’s inference times are highly variable and longer. For example, nnUNet3D takes about 3000 s (50 mins) for segmenting a T1w image (not shown in the plot for clarity), while soft_all requires only about 2 min. Unsurprisingly, the inference time largely depends on the size of the input volume. In addition to obtaining good segmentations, the average inference time per image is an important factor for consideration before deploying the model in real-world clinical settings. Models such as nnUNet3D requiring long inference times are impractical when used on large cohorts.

Discussion

We presented an automatic method for the contrast-agnostic soft segmentation of the spinal cord. Starting with the creation of GT masks, we proposed a new framework for generating a unique (soft) GT that represents the segmentations across various MR contrasts. Using these soft GT masks along with the corresponding images of all contrasts as inputs, we trained a U-Net model with aggressive data augmentation, regression-based adaptive wing loss, and NormReLU as the final activation function for spinal cord segmentation. We evaluated our model against 3 categories of baseline models to assess the impact of soft vs. hard GT masks, training a single model with all contrasts vs. training one model per contrast, and the impact of the loss functions on the subsequent morphometric measures (i.e., CSA). We compared our model’s performance with the SOTA methods for spinal cord segmentation, namely, PropSeg, DeepSeg, SoftSeg, BigAug and nnUNet. To demonstrate our model’s domain generalization capabilities on unseen contrasts and images with pathology (lesions), we presented a qualitative comparison against spinal cord segmentations of the SOTA methods. Lastly, we provided a graphical illustration of the average inference times of the existing methods and highlighted that our model takes only a fraction of the time per image irrespective of the contrast, while obtaining better segmentations that reduce the bias/variability in morphometric measures.

Preprocessing for soft GT

Our preprocessing framework to generates a unique soft segmentation from individual hard segmentations. The procedure required aligning all the images to the T2w image space, including resampling and reorientation. The soft segmentation was then brought back in each contrast’s native space. Then, in order to have a fixed patch size for training with all contrasts, we performed an up-sampling to 1 mm isotropic on all images and labels as a preprocessing step during data augmentation. One could question the reason for this up-down-up sampling in the preprocessing workflow, a few points are discussed below:

Comparison with baselines : In order to provide a fair comparison with the baselines, we ensured that all methods used the ground truth defined in the native space of each contrast.

Computation of CSA and disc labeling : Continuing the previous point and considering the evaluation of CSA variability across contrasts, it was important to compute the CSA in the native space, notably because of the highly variable spatial resolution across contrasts, which limited our ability to ensure that vertebral coverage would be the same when computed in the native image vs. in the reference image (T2w 0.8 mm isotropic).

Patch size : To ensure a uniform patch size during training for all the 6 contrasts, resampling the images to a unique, common, resolution was necessary. After experimenting with various spatial resolutions (trade-off between computational resource and precision), we chose 1 mm as the target resolution.

Introduce realistic variability : One of the advantages of training the model in the native space is that, for a given subject, the spinal cords are not perfectly aligned across contrasts.

Variability of CSA across contrasts

When comparing the performance on all 6 MR contrasts ( Fig. 3 ), DWI and MT-on contrasts yielded the highest absolute CSA error with 2.37 ± 1.09 mm 2 and 2.38 ± 1.77 mm 2 , respectively. This can be attributed to the fact that DWI and MT-on contrasts have less well defined boundaries between the spinal cord and the cerebrospinal fluid due to the coarser resolution of the images (0.9 × 0.9 × 5 mm). Furthermore, the presence of susceptibility related artifacts on DWI and MT-on (MT-on is based on GRE readout and hence suffers from signal dropout), and/or the ghosting and motion artifacts that particularly affect MT-on data ( Cohen-Adad et al., 2021b ) could explain the larger errors.

The agreement between T1w and T2w CSA ( Fig. 10 ) leads to a linear equation given by 0.95 x + 3.89 with an r 2 = 0.95, which is very close to the (dashed) identity line. In Cohen-Adad et al. (2021b) where sct_deeepseg_sc was used for generating spinal cord segmentations, the authors reported poor agreement between CSA computed on T1w and T2w images, regardless of the MR vendor.

Furthermore, training an individual model for each contrast ( soft _per_contrast ) yielded similar performance (albeit with a slightly higher error for each contrast) compared to training a single model for all contrasts ( soft_all ), as seen in Fig. 5 . Training an individual model for each contrast does not help in mitigating the CSA variability across contrasts as each model is optimizing for spinal cord segmentations asynchronously, thus leading to different CSA for a given participant despite using soft GT segmentations. On the other hand, a single model trained on all contrasts together is exposed to the wide heterogeneity in the images, thus leading to better estimation of the CSA and better generalization (as shown by our results on ms-mp2rage and radiculopathy-epi datasets in Figs. 12 and 13 ). Moreover, from a deployment standpoint, packaging and distributing a single model is more convenient and intuitive for researchers to use.

Overall, the proposed contrast-agnostic soft_all model outperforms the baselines and existing state-of-the-art spinal cord segmentation methods while minimizing contrast-specific CSA biases.

Effects of ground truth masks and loss functions

Results from Figs. 4 , 5 , and 6 demonstrated that the types of GT masks and loss functions used play a crucial role in the (downstream) computation of CSA. Using the unique soft GT generated by averaging the segmentations across 6 contrasts compared to traditional hard GT leads to lower CSA variability on the predictions. For more details, see Fig. 4 (qualitative assessment) and Fig. 5 (statistical assessment). Notably, the bias inherent to the individual GT for each contrast is propagated less when using the unique soft GT averaged from different contrasts.

For soft_all , the CSA did not differ significantly between T1w and T2w contrasts, while it was significantly different between the T1w/T2w and the other contrasts (T2*w, GRE-T1w, MT-on, DWI). Interestingly, T1w and T2w images have similar isotropic resolutions (respectively 1 × 1 × 1 mm and 0.8 × 0.8 × 0.8 mm), whereas the other contrasts feature highly anisotropic axial acquisitions (> 3 mm slice thickness). It is therefore possible that excessive partial volume effect along the superior–inferior axis created a bias in the CSA estimation. Another possible explanation of the discrepancy between isotropic and anisotropic scans, is the uncertainty in the estimation of the C2-C3 vertebral labels across contrasts. Since it was not possible to directly label the vertebral levels on the anisotropic scans (because discs are poorly visible), we used the disc labels created from the T2w images and applied the warping field to the labels to the target contrast. This likely resulted in slightly higher CSA STD across contrasts.

When the Dice or DiceCE loss functions ( Milletari et al., 2016 ; Isensee et al., 2021 ) were used in combination with hard GT masks, our results suggested that using Dice metric in the training objective is not sufficient for achieving accurate segmentations at the spinal cord/cerebrospinal fluid boundary. In fact, the model trained with Dice loss ( hard_all_BigAug ) showed subtle under-segmentations upon quality control as supported by larger absolute CSA errors in Fig. 6 . Instead, using adaptive wing loss that switches to the logarithmic (nonlinear) part of the loss when the error between the prediction and the GT are small, helps the model in refining the segmentations at the boundaries of the spinal cord, thus mitigating the CSA bias across contrasts. Similar observations have been reported about the effectiveness of regression-based ( Gros et al., 2021 ) and logarithmic ( Kaul et al., 2021 ) loss functions.

The benefits of using soft GT and adaptive wing loss in our model soft_all can be seen in Figs. 7 and 8 . PropSeg, DeepSeg, SoftSeg and BigAug methods, which used hard GT that are inherently biased due to the procedure of their GT generation, resulted in higher STD across contrasts per participant. As nnUNet does not support soft training yet, using soft segmentations averaged from all contrasts (but binarized at 0.5 threshold), still resulted in slightly higher CSA variability. Furthermore, DeepSeg used Dice loss and nnUNet used DiceCE loss during training, thus explaining the larger errors. Thus, the subtle yet important difference of training on implicitly obtained soft masks via data augmentation vs. applying the augmentation transforms directly on the soft masks generated from multiple contrasts has a significant downstream impact on the reduction of CSA variability across contrasts.

Generalization to unseen data

The proposed soft_all model demonstrates remarkable generalization on the unseen MP2RAGE “UNI” ( ms-mp2rage ) and resting state GRE-EPI ( radiculopathy-epi ) contrasts. Despite being trained only on healthy participants, it performed well on patients with MS lesions, traumatic spinal cord injury and cervical radiculopathy. We noticed that the tested models performed similarly on a contrast that was included during training (T2w) as reflected by the Dice coefficients of Table 3 and segmentations in Fig. 11 , even in the presence of traumatic spinal cord injury. Surprisingly, there is a marked difference between the qualitative segmentations of ( soft_all ) and nnUNet3D ( Figs. 12 and 13 ) on the ms-mp2rage and radiculopathy-epi datasets, respectively. With the former dataset, nnUNet3D performs poorly by completely missing the spinal cord in the presence of multiple sclerosis lesions, whereas for the latter, we observed cases with oversegmentation of the cord leaking into the cerebrospinal fluid. For both these contrasts, soft_all obtains accurate segmentations of the spinal cord under the presence of lesions and along the spinal cord boundaries.

The difference in the segmentations between our model and nnUNet is likely due to our improved training strategy involving cropping along the center of the spinal cord, regression-based adaptive wing loss, and most importantly, training directly on the soft GT masks (unlike in nnUNet where the soft GT were binarized). The localization of the spinal cord, mainly through cropping, has been a recurrent prerequisite step in the literature ( Gros et al., 2019 ; Benjdira et al., 2020 ; Lemay et al., 2021 ), suggesting its necessity for obtaining robust segmentations.

Furthermore, our model soft_all is able to better delineate the shape of the spinal cord even in the presence of lesions compared to models trained specifically on one contrast ( soft_per_contrast T2w and DeepSeg2D ). This enhanced performance can likely be attributed to the model’s comprehensive exposure to diverse contrasts of the spinal cord, cerebrospinal fluid, gray matter and white matter, that are featured in the Spine Generic ( Cohen-Adad et al., 2021a ) dataset.

Limitations & future work

Application to thoracic and lumbar levels

The proposed model is trained on a dataset of healthy participants containing cervical and upper-thoracic spinal cords only. Future research will add images with lumbar cords to further improve the generalizability of the model towards different FoVs.

Center cropping

The center cropping step in the online preprocessing of the images during training and inference assumes that spinal cord is at the center of the image. While this is a reasonable assumption, there might be some edge cases containing lumbar spines or participants with scoliosis on which the automatic predictions might fail. In order to mitigate this issue, the SCT function that runs inference with the proposed model provides a flag to change the default cropping (allowing researchers to adjust the cropping sizes based on their images).

Binarization threshold

At prediction time, the soft output was binarized for comparison with other methods. That binarization was done using a 0.5 threshold. As discussed in Gros et al. (2021) , that threshold could potentially be optimized to further reduce the variability of CSA across contrasts. However, this would imply arbitrarily categorizing images into a given contrast, which defeats the purpose of the current contrast-agnostic method, wherein the model can be used as is regardless of the acquisition parameters. Moreover, with MRI acquisitions, it is difficult to define the contrasts accurately as, for instance, some combinations of parameters could lead to more/less T2w or more/less magnetization transfer saturation depending on the offset of the MT pulse and/or the presence of saturation bands.

Validation in pathologies

As mentioned in the introduction, one of the advantages of soft segmentation is the ability to encode volumetric measures with finer precision compared to binary segmentation. Considering that changes in the spinal cord happen at a very slow rate at early stages of MS, we expect that soft segmentation of the spinal cord will help detect subtle atrophies. We expect that the soft segmentation of the spinal cord will increase the precision of spinal cord CSA measurements, and thus lead to lower arm-size in trials ( Bautin and Cohen-Adad, 2021 ; Cawley et al., 2018 ) . This will especially be of interest in the context of cross-sectional studies where protocols can vary. In patients with degenerative cervical myelopathy, which is characterized by a progressive compression of the spinal cord, being able to precisely segment the spinal cord could also lead to better prognostication and therapeutic strategies ( Martin et al., 2018 ).

Continual model refinement through active learning

The improvement in the mitigation of biases in morphometric measures between MRI contrasts holds exciting prospects for future work. Thanks to the remarkable generalization of the proposed model to unseen contrasts, datasets containing other contrasts (e.g., phase-sensitive inversion recovery, short tau inversion recovery, susceptibility-weighted imaging, MP2RAGE) can be added to the existing cohort to further improve the generalizability of the model. Moreover, enriching the model with images from patients with spinal pathologies will likely improve zero-shot generalization on participants with lesions and/or spinal cord compressions. For these advancements, human-in-the-loop active learning ( Budd et al., 2021 ) involving an initial batch of segmentations from our model followed by manual corrections of under-/over-segmentations and then re-training the model on the larger datasets until fully automatic predictions are obtained is an attractive strategy. This approach for gradually aggregating large-scale datasets while improving the model simultaneously is similar to the recently proposed Segment Anything ( Kirillov et al., 2023 ) model, thus paving way towards a foundational model for contrast-agnostic spinal cord segmentation.

Conclusion

We presented an automatic method for the soft segmentation of the spinal cord across various MRI contrasts. Using a novel framework for generating soft GT masks that leverages information from all contrasts along with aggressive data augmentation and a regression loss function, our model produces soft segmentations that encode partial volume information at the spinal cord and cerebrospinal fluid boundary and produces segmentations that are stable across contrasts. More importantly, the proposed model reduces variability in CSA across all contrasts and generalizes well to unseen contrasts, pathologies, and MRI vendors. Overall, the proposed method could potentially increase the sensitivity of CSA as a biomarker to detect subtle spinal cord atrophy ( Bautin and Cohen-Adad, 2021 ).

Supplementary Material