A multi-contrast MRI approach to thalamus segmentation
Veronica Corona Jan Lellmann Peter Nestor Carola-Bibiane Schoenlieb Julio Acosta-Cabronero







 A multi-contrast MRI approach to thalamus segmentation
 
Veronica Corona^a*, Jan Lellmann^b, Peter Nestor^c,d, Carola-Bibiane Schönlieb^a,1, Julio Acosta-Cabronero^e,f,1



^a Department of Applied Mathematics and Theoretical Physics, University of Cambridge, United Kingdom


^b Institute of Mathematics and Image Computing, University of Lübeck, Germany

^c Queensland Brain Institute, University of Queensland, Brisbane, Australia

^d Mater Hospital, South Brisbane, Queensland, Australia

^e Wellcome Centre for Human Neuroimaging, UCL Institute of Neurology, University College London, United Kingdom

^f German Center for Neurodegenerative Diseases (DZNE), Magdeburg, Germany


^* Corresponding author. Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, Wilberforce Rd, Cambridge, CB3 0WA, United Kingdom. E-mail address: vc324@cam.ac.uk 

^1 Shared senior role 

Keywords: thalamus segmentation, thalamic nuclei, magnetic resonance imaging, multi-contrast segmentation

Thalamic alterations are relevant to many neurological disorders including Alzheimer's disease, Parkinson's disease and multiple sclerosis. Routine interventions to improve symptom severity in movement disorders, for example, often consist of surgery or deep brain stimulation to diencephalic nuclei. Therefore, accurate delineation of grey matter thalamic subregions is of the upmost clinical importance. MRI is highly appropriate for structural segmentation as it provides different views of the anatomy from a single scanning session. Though with several contrasts potentially available, it is also of increasing importance to develop new image segmentation techniques that can operate multi-spectrally. We hereby propose a new segmentation method for use with multi-modality data, which we evaluated for automated segmentation of major thalamic subnuclear groups using T_1-, T_2^*-weighted and quantitative susceptibility mapping (QSM) information. The proposed method consists of four steps: highly iterative image co-registration, manual segmentation on the average training-data template, supervised learning for pattern recognition, and a final convex optimisation step imposing further spatial constraints to refine the solution. This led to solutions in greater agreement with manual segmentation than the standard Morel atlas based approach. Furthermore, we show that the multi-contrast approach boosts segmentation performances. We then investigated whether prior knowledge using the training-template contours could further improve convex segmentation accuracy and robustness, which led to highly precise multi-contrast segmentations in single subjects. This approach can be extended to most 3D imaging data types and any region of interest discernible in single scans or multi-subject templates. 





§ INTRODUCTION



The thalamus is composed of a complex set of sub-nuclei. It is considered the central relay station for sensory and motor information as nearly all sensory and motor signals are sent to the thalamus prior to reaching the cortex. It is also thought to have an integrative role as thalamic structures receive, process, sort and send information between specific subcortical and cortical areas, and might be involved in regulation of sleep and wakefulness, memory, emotion, consciousness, awareness, attention, ocular motility, learning and motor control processes <cit.>. 

The thalamus is composed of several major substructures. The internal medullary lamina is a thin sheet of white matter that runs longitudinally through the thalamus separating it into medial and lateral regions. In the anterior part, the internal medullary lamina forks to isolate the anterior thalamic nucleus; thus, thalamic nuclei can be broken down into four regions based on their position relative to the lamina, i.e., anterior, medial, lateral and posterior subnuclear groups <cit.>. 
Lesions to thalamic nuclei and their connections to the cortex can result in a wide range of neurological deficits. Thalamic alterations have been identified in several neurodegenerative diseases including Alzheimer's disease, Parkinson's disease, Huntington's disease and multiple sclerosis, the majority of which present evidence of atrophy in one or more substructures of the thalamus <cit.>. Neurological patients also often undergo brain surgery and deep brain stimulation targeting thalamic subnuclei, thus accurate and reliable localisation of such structures are key both for research and for delivering effective clinical treatments <cit.>.
New developments in imaging techniques, including 3-7 Tesla MRI, provide greater contrast and higher spatial specificity to study the thalamus. Therefore, new strategies need to be investigated for clinical and research applications, which could potentially lead to suitable tools for predicting cognitive impairment and for monitoring disease progression in neurological patients <cit.>.
To date, several methods have been proposed to delineate subthalamic regions with MRI, a few of which used diffusion MRI. For example, Behrens et al. <cit.> described a diffusion tensor imaging (DTI) based segmentation procedure based on coarse tractography patterns from the thalamus to the cortex, and Wiegell et al. <cit.> developed a k-means clustering algorithm to detect groups of coherent DTI-based fiber orientation. The use of the mean-shift algorithm <cit.> has also been proposed, whereby regional clusters and shapes are inferred from the local modes of a density estimator computed with a multivariate kernel <cit.>. Furthermore, Jonasson et al. <cit.> proposed a level-set method whereby a region-based force (defined using a diffusion similarity index between the most representive tensor of each level set and its neighbors) drives a set of coupled level-set functions each representing a segmented region. High angular resolution diffusion images (HARDI) have also been investigated for segmenting the thalamus. Grassi et al. <cit.> proposed a k-means clustering approach whereby a specific number of initialised centroids are updated based on a weighted sum of the Euclidean distance of voxels to centroids and Frobenius distance of their orientation distribution function. Notably, however, all diffusion MRI based methods are hampered by low spatial resolution. In an attempt to overcome this limitation, Deoni et al. <cit.> explored with some success the use of high-resolution quantitative MRI, namely T_1 and T_2 mapping, with a modified k-means clustering approach that combined T_1/T_2 information and center-of-mass distances to Morel atlas segmentations <cit.>. Further using anatomical images, Magon et al. <cit.> developed a method to segment thalamic subnuclei employing a voxel-wise majority vote after co-registration to multiple atlases.
Past efforts also focused on the MRI acquisition. Bender et al. <cit.>, for example, proposed an inversion time optimisation strategy to enhance the T_1-weighted contrast between gray and white matter using the 3D magnetization-prepared rapid acquisition of gradient echo (MPRAGE) sequence. Tourdias et al. <cit.> subsequently optimised MPRAGE for 7T MRI and proposed imaging at the white matter null regime both for enhancing the contrast between the thalamus and surrounding tissues and for depicting several subnuclear groups. 
Thalamus segmentation with quantitative susceptibility mapping (QSM) - a relatively new quantitative MRI contrast - has also gained increasing interest in recent times. Deistung et al. <cit.> illustrated that high-resolution QSM is a superior contrast to depict thalamic substructures than T_2^*-weighting, the local field or R_2^* maps. Therefore, considering QSM's ability to provide quantifiable information about iron content <cit.>, that iron accumulation has been associated with several neurological disorders <cit.> and that thalamic lesions are not uncommon in such disorders <cit.>, it is highly plausible that enabling reliable segmentation of thalamic substructures could have a major impact on the study of ageing and disease.
Traditionally, however, the anatomy of the thalamus has been inferred from post-mortem tissue examinations. The most widely used histological atlas was developed by Morel et al. <cit.> using an iterative approach for reconstructing the mean model from six series of maps derived from different stacks of histologically processed brains. The model, thus, is an average template incorporating topological and geometric features from only a few individuals. Morel's and other similar proprietary atlases are widely used for guiding MRI-based segmentations in neurosurgical planning, although notably, the direct superposition onto brain scans is often not fit for precision measurements, a situation often aggravated by age-related differences <cit.>. Therefore, the development of image-guided segmentation approaches are highly relevant in this context. 




  
Our contribution. This work proposes a new multi-contrast segmentation algorithm, and its optimisation, to exploit the full potential of T_1-, T_2^*-weighted and QSM contrast differences in thalamic subregions. 
We show that using multi-contrast information maximizes segmentation performance, by exploiting structures that become visible and enhanced in specific MR imaging protocols. In the proposed method, regions of interest defined in template space are learnt and then approximated in single subjects with spatial constraints that ensure robustness. Our multi-contrast segmentation framework can be extended to any data types and regions of interest.





§ METHODS


The proposed semi-automated method consists of four steps: spatial normalisation, manual (reference) segmentation, pattern recognition and a final refinement step using a convex formulation.
Details on study subjects, MRI acquisition and pre-processing are given below. For now, we will assume all subject data has been spatially co-registered to a common reference space, from which multi-subject templates (one for each contrast) have been computed. We will also assume hereafter (specific details given below) that region(s) of interest has(have) been manually traced (at least once) with the aid of such templates. We then consider the following multi-class labelling problem.



 §.§ Classification



For each voxel in the image volume domain Ω⊂ℝ^3, Ω={1, …, n_1}×{1, …, n_2}×{1, …, n_3}, we assign one of ℓ class labels, with each class referring to a segmented region. Let X={x_i, i=1,…, n}, where n=n_1 n_2  n_3, be the vectorised volume in template space. For each x_i, we have c image intensities or MRI parameter values, f_1(x_i), …, f_c(x_i), one from each imaging contrast available; in this study, T_1-, T_2^*-weighted signals and QSM. We then identify a set of possible class labels, {0, 1, …, ℓ-1}; in this particular context, we set 0 to be the background region, 1 the lateral thalamic subnuclear group, 2 the medial group and 3 the posterior group. The manual segmentation in template space is required to define the label set for the volume X as Y={y_i, i=1,…, n}, where y_i ∈{0,1,2,3}.



  
Feature space. In the reference coordinate system we then build the feature space: Ψ={ψ_ij, i=1, …,n;  j=1,…, m}, assigning x_i,  i=1,…,n, i.e., an m-dimensional feature vector, to each voxel. Features describe objects, in our case voxel information reflecting thalamic tissue properties. In this work we set out to develop a multi-spectral approach to exploit features from several contrast types, whereby the key features are intensity/MRI parameter values: f_k(x_i), k=1,2,3 from T_1, T_2^*-weighted MRI and QSM, which return different contrast characteristics for tissues with different local concentrations of water, iron, myelin, etc. For each contrast, we also selected additional features which are the result of an empirical study of the feature space. These are:  mean, μ, and standard deviation, σ, across the 26-neighbourhood, and intensity/MRI parameter values for the six closest 3D neighbours in each contrast, leading to a feature space of m=27 dimensions. All features were then scaled by their normalised variance (i.e., with mean shifted to the origin and total variance for all features scaled to 1). It should be noted that this feature space was optimised through an investigation of classification accuracy versus feature space dimensionality on a data subset. This might differ for other data types and/or target regions.




  
Classifiers. Let us consider the feature space Ψ and the label set Y for n voxels. Each template voxel is therefore described by the pair (ψ_i, y_i), where ψ_i is the m-dimensional feature vector of voxel x_i and y_i is its label. We define the labeled training dataset as 𝒯={(ψ_1,y_1),...,(ψ_n,y_n)}. We set out to solve a classification problem based on supervised learning, in which we train a classifier to derive a decision mapping for new observations. Initially, we explored the performance of several classification methods including Support Vector Machine, Random Forest, Naive-Bayes, k-Nearest Neighbours (k-NN) and Parzen classifiers using a data subset. In this preliminary study, we obtained greater accuracy with two classifiers: k-NN (k=3) and Parzen classifiers, which are described in <cit.>. As pointed out in the context of feature selection, the optimal choice of training classifier may also vary according to data type and/or target region. 







 §.§ Convex segmentation


Classification routines yield a posterior probability distribution p̂(u|f) ∈ℝ^n ×ℓ for each class, that is the probability for voxel x to be assigned class u(x)=l given the measured data f(x). From this, winner-takes-all segmentation can be derived selecting the class with the highest probability value in each voxel. This, however, often results in scattered clusters of misclassified voxels that break the smoothness and continuity of segmented regions. Hereby, therefore, we introduce an additional convex optimisation step to further improve the spatial cohesiveness of tissue segments.
More precisely, we consider a labelling function u: Ω→ℝ^ℓ that represents the unique assignment of a label to each voxel x in the image domain Ω. Because this is a hard combinatorial problem, we use a convex relaxation to facilitate the optimisation, see <cit.> for an overview. 
The notion of labelling function is relaxed to u taking values in the convex set defined by the unit simplex  Δ_ℓ := { u ∈ℝ^n →ℝ^ℓ |  u ≥ 0, ∑_i=1^ℓ u_i =1 }. Then, by choosing J convex, we solve the following convex segmentation (CS) problem:


    min_u:Ω→Δ^ℓ∑_x ∈Ω - log p̂(u(x) | f(x)) _data term+ λTV(u)_regulariser,

whereby the data term is the negative logarithm of the posterior probability distribution computed by the classifier, and the regulariser is the total variation (TV) of the labelling function u defined as the L^1-norm of a discrete finite-difference approximation of the two-dimensional gradient of u. The TV regulariser on the relaxed u is the convex equivalent to the length penalty on the original hardcoded labelling function and, as such, it can be thought of as introducing a penalty for long or irregular interfaces between different classes. The parameter λ > 0 balances the data term and the regulariser in the minimisation. We solve (<ref>) using the fast primal-dual algorithm described in <cit.>.





 §.§ Convex segmentation with additional priors


Individual datasets are overall inferior to group-wise templates in terms of signal- and contrast-to-noise ratio. With a view then to ensure segmentation robustness at the single-subject level, we extended the forward model by the introduction of a priori information on the manual segmentation of the training template. We enabled the weighting of posterior probabilities, p̂(u | f), according to template-based constraint as follows:


    min_u:Ω→Δ^ℓ∑_x ∈Ω - log( (1-w)   p̂(u(x)|f(x)) + w   n ) +λTV(u)

where w ∈ [0,1] is a normalised weight determining the level of prior information constraining the data term, n:Ω→{0,1}^ℓ is a labelled mask of thalamic subregions, and λ >0 is the regularisation parameter. 





 §.§ Study subjects



Training and test datasets consisted of N=43 (age: 59± 7, [50-69] years old, 19 female/24 male) and N=116 (age: 54± 19, [20-79] years old, 56 female/60 male) healthy subjects, respectively. The latter was an aging cohort previously investigated with QSM <cit.>. All elderly subjects (age>50 y.o) were free of neurological or major psychiatric illness and performed normally on cognitive screening: mini-mental state examination (MMSE>27) <cit.>. 





 §.§ MRI scanning protocol


The imaging protocol, QSM reconstruction and spatial normalization methods (briefly summarised below) are essentially identical to those in a previous aging study <cit.>. 

All participants were scanned on a Siemens Verio 3 Tesla MRI system with a 32-channel head coil (Siemens Medical Systems, Erlangen, Germany) under the following imaging protocol: 
T_2-weighed fast spin echo images were acquired and visually inspected to ensure vascular pathology was not significant in any subject as for standard clinical routine. Scan parameters were as follows: flip angle (α)/ echo time (TE)/ receiver bandwidth (RB)/ turbo factor/ number of echo trains per slice/ echo spacing/ repetition time (TR) = 150^∘/ 96 ms/ 220 Hz per pixel/ 18/ 13/ 9.64 ms/ 8160 ms; matrix, 320×320 (in-plane resolution: 0.7×0.7  mm^2); 45 axial slices for whole-brain coverage (thickness: 3 mm; gap: 0.9 mm); two-fold parallel acelleration was enabled for phase encoding giving a total scan time of 1:56 minutes.

T_1-weighed data were acquired using a 3D MPRAGE sequence <cit.> with the following scan parameters: inversion time/ α/ TE/ RB/ echo spacing/ TR = 1100 ms/ 7^∘/ 4.37 ms/ 140 Hz per pixel/ 11.1 ms/ 2500 ms; 256×256×192 matrix dimensions (straight-sagittal orientation), 1×1×1 mm^3 voxel size, two-fold parallel acceleration and further 7/8 partial Fourier undersampling for phase encoding. The total scan time was 5:08 minutes.

T_2^*-/susceptibility-weighted data were obtained from a fully flow-compensated, 3D spoiled gradient-echo sequence. Scan parameters were: α/ TE/ RB/ TR = 17^∘/ 20 ms/ 100 Hz per pixel/ 28 ms; matrix, 256×224×80 with voxel resolution of 1×1×2 mm^3; and two-fold parallel acceleration for phase encoding. The total scan time was 5:32 minutes. All susceptibility maps were inspected to exclude subjects with severe calcifications or extensive haemorrhages.




 §.§ MRI data pre-processing




  §.§.§ QSM reconstruction


QSM is a relatively new contrast modulated by the local content of chemical species that have different magnetic susceptibilities than soft tissue and water <cit.>. Myelin phospholipids and calcifications, for example, are more diamagnetic than water; whereas, iron - the most abundant transition metal in the human brain and the dominant source of QSM contrast - is greatly paramagnetic <cit.>. Specific details on susceptibility reconstruction can be found elsewhere <cit.>.




  §.§.§ Spatial standardisation


Radio-frequency (RF) bias corrected <cit.> T_2^*-weighted magnitude images were affine co-registered to their corresponding bias-corrected MPRAGE volume using ANTs (<http://stnava.github.io/ANTs/>) <cit.>. Subsequently, all bias-corrected anatomical T_1-weigthed MPRAGE images were used to generate a study-wise space using a previously described ANTs routine <cit.>. Finally, all T_2^*-weighted images and susceptibility maps were normalised to the same coordinate system through the warp composition of the above transformations and high-order interpolation.



  §.§.§ Manual (reference) segmentation


Three templates were subsequently obtained from averaging T_1-, T_2^*-weighted and QSM data across subjects in the study-wise space (<ref>). This was performed separately for training and testing data.

Then, three major thalamic subregions, namely lateral, medial and posterior nuclear groups, were manually traced as illustrated in <ref>. The manual annotations from the training-average template were utilised in the training phase of the algorithm as ground truth. In addition, thalamic subregions from the average test template (N=116) and for N=6 individual test datasets were delineated for algorithm validation (see next subsection).





 §.§ Performance evaluation
 

Segmentation performance was assessed through visualisation of the confusion matrix (incorporating exact error distributions). For simplicity, however, in this study we report two representative summary measures: the global classification error rate (i.e., the overall proportion of erroneously classified voxels) and the true positive (TP) rate for non-background (i.e., subnuclear group) regions. At the outset, general performance was evaluated on the high-contrast N=116 template dataset, which included a comparison with standard Morel atlas based segmentation. Subsequently, error measures were also computed for individual (noisier) test data.



 §.§ Methods summary


In Algorithm <ref>, we summarise the proposed methodology for multi-contrast MRI segmentation. The first stage of the algorithm trains a classifier for use in stage two. Given then a "new" multi-contrast MRI dataset to be segmented, all contrast images must be first realigned to a common space, then the "new subject" segmentation pipeline can be applied as follows:


  * supervised learning (testing), given the trained classifier, its mapping is applied to independent test data yielding class labels and posterior probabilities.

  * multi-class convex segmentation, where posterior probabilities are used in the data term of the convex optimisation formulation defined in <ref>.


The supervised learning and convex segmentation steps of the algorithm were implemented in MATLAB R2017b (The Mathworks Inc., Natick, MA, USA) and are available at
 <https://github.com/veronicacorona/multicontrastSegmentation.git>.



All MRI datasets in this study were spatially standardised via nonlinear co-registration to a common coordinate system. This enabled custom training from a single set of regional contours in template space. Future applications of this algorithm could alternatively consider using manual tracings from each individual training dataset. The only requirement is that all contrasts for a given subject must share a common frame of reference.




§ RESULTS
 

In what follows we present our numerical results obtained independently for the N=116 test dataset described in <ref>.



 §.§ Qualitative assessment


In this implementation, classifiers were set out to assign four posterior probabilities per voxel, i.e., those of belonging to background, lateral, medial and posterior subregions of the thalamus. <ref>-<ref> shows posterior probability maps using the Parzen classifier on N=116 average-template data. The figure indicates that accurate classification of specific subnuclear groups and the background region is feasible; supporting, thus, the choice of feature space and classifier. Overall, the best performing algorithms in our preliminary assessment were k-NN and Parzen classifiers. For k-NN, the optimal number of nearest neighbours, k, was k=3. For the Parzen classifier, the empirically optimal parameter h, i.e., the width of the Gaussian smoothing kernel, was h=0.1668. <ref> further confirmed that the Parzen classifier output is overall in agreement with a priori knowledge on the regional distribution of subnuclear groups. However, as predicted, winner-takes-all local classification introduced undesirable regional discontinuities. This was substantially mitigated through the additional CS step as shown in <ref>.




























































 §.§ General performance evaluation




  §.§.§ Convex segmentation validation


Indeed, as shown in <ref>, the introduction of convex segmentation systematically improved classification performance, which, in turn, also confirmed that posterior probability maps from both k-NN and Parzen classifiers are suitable pre-conditioners for the CS formulation in (<ref>).




  §.§.§ Algorithm comparison


<ref> illustrates segmentation results for all methods herein evaluated. Outputs from the proposed multi-contrast method were in greater agreement with the manual ground truth than atlas-based Morel segmentation, which is solely based on template co-registration.

It is worth noting that in this particular implementation the background region outsizes (approximately 4:1) the extent of putative thalamic regions. Therefore, segmentation results are reported in <ref> both as global classification errors and true positive rates; the latter computed for non-background regions only. Such an evaluation confirmed the proposed algorithm outperforms Morel atlas segmentation on all performance metrics. Pre-conditioning with 3-NN and Parzen based posterior probabilities minimised classification error and true positive rates, respectively.




  §.§.§ Regularisation parameter selection for convex segmentation


The 3-NN and Parzen based segmentation results in <ref> and <ref> were obtained through solving the convex optimisation problem defined in (<ref>), which has a regularisation multiplier, λ, that requires optimisation for optimal solution smoothness. In this study, λ was optimised empirically on a small subset: for 3-NN, we chose λ=1, and for Parzen λ=5. We then confirmed the validity of this choice calculating overall classification errors (on the N=116 template dataset) for a range of regularisation parameters. Results from this validation experiment are summarised in <ref>.






  §.§.§ Algorithm performance as a function of input data


A unique aspect of the proposed algorithm is that it can integrate information from several MRI contrasts capturing simultaneously different views of the anatomy. In this study, we hypothesised that T_1-, T_2^*-weighting and QSM all provide differentially relevant information about internal thalamic boundaries. In order to substantiate this claim, we estimated algorithm performance for all the available combinations, i.e., one, two or three data types, using the same 27-dimensional feature space that was previously optimised. CS errors are shown in <ref> for pre-conditioning with both 3-NN and Parzen classifiers. Interestingly, using single contrasts alone as input data led to sytematically greater error rates than when using QSM in combination with other contrast types. Confirming our hypothesis, the best segmentation performance was obtained when using all MRI information. Although some differences were observed, overall both pre-conditioning approaches, i.e., 3-NN and Parzen classification, yielded relatively similar error rates throughout.




























































































































































































 §.§ Convex segmentation with additional priors for increased performance in single subjects




We also confirmed that constraining the data term for fidelity with training-average tissue priors is feasible and desirable to improve accuracy and robustness in single-subject thalamic subnuclei segmentation. The consistency weight, w in (<ref>), represents a trade-off between the calculated posterior probabilities and template-based priors. On a N=6 test subset, we explored how classification error varies as a function of w. This is illustrated in <ref>, which indicated that 3-NN is generally preferred (to Parzen based training) in this context; 3-NN pre-conditioning resulted in greater accuracy with optimal performance for a critical w=0.4. 

















































































From a qualitative standpoint, <ref> illustrates the improvement in single-subject segmentation when weighting the fidelity term by the training-data based prior. Weighted segmentations with both supervised learning approaches (i.e., 3-NN and Parzen) converged to solutions that were overall in close agreement with the manual ground truth. Extended results for the same subject using 3-NN are shown in <ref>. 









§ DISCUSSION

In this study, we present a data-driven method to segment several internal thalamic boundaries using multi-contrast MRI data. We had three imaging contrasts available to drive the segmentation procedure: T_1, T_2^*-weighting and QSM. We confirmed that using all information maximised segmentation performance, and also found evidence suggesting that QSM was the most informative contrast type for this thalamic problem. Different data types or new implementations for other anatomical regions will require new calibration work.
This work was motivated by the observation that study-wise MRI templates obtained using highly iterative non-linear coregistration routines are showing superb anatomical detail over and above what can be inferred from individual datasets. It is therefore unsurprising these are being used to trace regions of interest that are not available from automated segmentation tools <cit.>. Although this is an effective strategy, it assumes both that there are no co-registration errors in the calculation of the study-wise template space and that the manual reference is an exact definition of the region of interest, which are somewhat inaccurate assumptions. In this work we broke away from this idealisation and propose to correct these errors with two additional steps: one of pattern recognition, followed by convex segmentation promoting (from a Bayesian standpoint) segmentation boundaries that are short, continuous and regular while also consistent with contrast variations in single subjects. Furthermore, to capitalize on the facts (i) that multiple MRI contrasts are typically acquired in the same imaging session, and (ii) that different MRI contrast types could act in concert to help resolve tissue boundaries, the algorithm was implemented in multivariate form. In turn, this new method yielded regional boundaries that were in good agreement with manually traced contours. This was in stark contrast with the output from Morel atlas based segmentation of the same subnuclear groups, confirming that data-driven approaches (such as that which is hereby proposed) signify an improvement (with respect to co-registration based atlas labelling methods) in terms of consistency with manual segmentation.
It is also worth noting that posterior probability maps from individual datasets can be noisy. In this study, with a view to make the convex segmentation algorithm more robust in this regime, we introduced an additional data fidelity weight in (<ref>) to enable additional prior knowledge to be transferred from the training reference to single-subject segmentations. Such an approach led to significant improvements for both classifiers, although we noted optimal performance (i.e., lower errors with respect to the gold standard) specifically for 3-NN based modelling and inclusion of 40 % prior knowledge. Intuitively, w-dependent errors reflect the complex interaction between co-registration performance, accuracy on training-template manual delineation and the algorithm's ability to identify biologically meaningful boundaries between tissue types. In other words, the finding that segmentation errors were systematically minimised by w<1 confirmed that the proposed algorithm effectively corrects for co-registration and/or manual initialisation errors. We cannot guarantee, however, that the proposed implementation (i.e., 3-NN classification with λ=1 and CS with w=0.4) will be optimal for other regions and/or data types. This warrants further investigation.
An additionally important consideration for early adopters of this method is that posterior probability maps can only be obtained from models trained on separate data. In this study, we had sufficient power to split the dataset into training and test subsets. However, future studies wanting to implement this algorithm with limited available data may consider e.g. an algorithmic extension for synthetic data augmentation.
In conclusion, this work presented a highly versatile multi-contrast segmentation framework and its successful application to identify thalamic substructures. Future work is warranted to extend this method for segmentation of other structures. In addition to developing appropriate forward pipelines for bootstrapping training data augmentation, further improvements might be obtained using e.g. deep learning, which may eliminate the need for additional spatial constraints.


 §.§ Data statement

The dataset used in this work and the proposed supervised learning and convex segmentation implementations are available at <https://github.com/veronicacorona/multicontrastSegmentation.git>. 



§ ACKNOWLEDGEMENTS

The Wellcome Centre for Human Neuroimaging is supported by core funding from the Wellcome (203147/Z/16/Z).
VC acknowledges the financial support of the Cambridge Cancer Centre and the Cancer Research UK. CBS acknowledges support from Leverhulme Trust project "Breaking the non-convexity barrier", EPSRC grant "EP/M00483X/1", EPSRC centre "EP/N014588/1", the Cantab Capital Institute for the Mathematics of Information, and from CHiPS and NoMADS (Horizon 2020 RISE project grant). Moreover, CBS is thankful for support by the Alan Turing Institute.

tocsectionBibliography
unsrt
  














































