Learning joint lesion and tissue segmentation from task-specific hetero-modal datasets

TL;DR

This paper introduces a novel deep learning approach for joint brain tissue and lesion segmentation from heterogeneous, partially annotated MRI datasets, effectively handling missing modalities and large anatomical variations due to pathology.

Contribution

It proposes a variational formulation and empirical risk decomposition enabling joint segmentation from task-specific hetero-modal datasets with missing modalities.

Findings

Achieves comparable performance to task-specific models.

Effectively handles missing imaging modalities.

Supports segmentation despite large anatomical changes.

Abstract

Brain tissue segmentation from multimodal MRI is a key building block of many neuroscience analysis pipelines. It could also play an important role in many clinical imaging scenarios. Established tissue segmentation approaches have however not been developed to cope with large anatomical changes resulting from pathology. The effect of the presence of brain lesions, for example, on their performance is thus currently uncontrolled and practically unpredictable. Contrastingly, with the advent of deep neural networks (DNNs), segmentation of brain lesions has matured significantly and is achieving performance levels making it of interest for clinical use. However, few existing approaches allow for jointly segmenting normal tissue and brain lesions. Developing a DNN for such joint task is currently hampered by the fact that annotated datasets typically address only one specific task and rely on a task-specific hetero-modal imaging protocol. In this work, we propose a novel approach to build a joint tissue and lesion segmentation model from task-specific hetero-modal and partially annotated datasets. Starting from a variational formulation of the joint problem, we show how the expected risk can be decomposed and optimised empirically. We exploit an upper-bound of the risk to deal with missing imaging modalities. For each task, our approach reaches comparable performance than task-specific and fully-supervised models.

Equations

$${\theta }^{\ast }={\text{argmin}}_{\theta }{\mathbb{E}}_{(X,Y)\sim \mathcal{D}}\left[\mathcal{L}\left(h_{\theta }(X),Y\right)\right]$$

$$\tag{\textbf{H}_{1}}\mathcal{L}(h_{\theta}(X),Y)=\mathcal{L}^{t}(h_{\theta}(X),Y^{t})+\mathcal{L}^{l}(h_{\theta}(X),Y^{l})$$

$${\theta }^{\ast }={\text{argmin}}_{\theta }\underset{{\mathcal{R}}^t(\theta )}{\underbrace{{\mathbb{E}}_{(X,Y)\sim \mathcal{D}}[{\mathcal{L}}^t(h_{\theta }(X),Y^t)]}}+\underset{{\mathcal{R}}^l(\theta )}{\underbrace{{\mathbb{E}}_{(X,Y)\sim \mathcal{D}}[{\mathcal{L}}^l(h_{\theta }(X),Y^l)]}}$$

$$\tag{\textbf{H}_{2}}\mathcal{D}\triangleq\mathcal{D}_{lesion}.$$

$$\tag{\textbf{H}_{3}}\forall(a,b,c)\in\mathcal{Y}^{3}:\ \mathcal{L}^{t}(a,c)\leq\mathcal{L}^{t}(a,b)+\mathcal{L}^{t}(b,c)$$

$${\mathcal{L}}^t(h_{\theta }(X),Y^t)\le {\mathcal{L}}^t(h_{\theta }(X),h_{\theta }(p(X)))+{\mathcal{L}}^t(h_{\theta }(p(X)),Y^t)$$

$${\mathcal{R}}^t(\theta )\le {\mathbb{E}}_{(X,Y)\sim {\mathcal{D}}_{lesion}}[{\mathcal{L}}^t(h_{\theta }(X),h_{\theta }(p(X)))]+{\mathbb{E}}_{(X,Y)\sim {\mathcal{D}}_{lesion}}[{\mathcal{L}}^t(h_{\theta }(p(X)),Y^t)]$$

$$\tag{\textbf{H}_{4}}\forall i\in\{1...N\}P_{\mathcal{D}_{lesion}}(x_{i},y_{i}|y_{i}\in\mathcal{C}_{tissue})=P_{\mathcal{D}_{control}}(x_{i},y_{i}|y_{i}\in\mathcal{C}_{tissue})$$

$${\mathcal{R}}^t(\theta )\le \underset{{\mathcal{R}}_1^t(\theta )}{\underbrace{{\mathbb{E}}_{(X,Y)\sim {\mathcal{D}}_{lesion}}[{\mathcal{L}}^t(h_{\theta }(X),h_{\theta }(p(X)))]}}+\underset{{\mathcal{R}}_2^t(\theta )}{\underbrace{{\mathbb{E}}_{(X,Y)\sim {\mathcal{D}}_{control}}[{\mathcal{L}}^t(h_{\theta }(p(X)),Y^t)]}}$$

$$\mathcal{L}(h_{\theta }(\text{x}),y)=1-\sum _{c\in C}{\omega }_c\frac{{\sum }_{j=1}^N{g}_{j,c}{p}_{j,c}}{{\sum }_{j=1}^N{g}_{j,c}^2+p_{j,c}^2-g_{j,c}{p}_{j,c}}\text{suchas}\sum _{c\in C}{\omega }_c=1$$

Full text

1 Introduction

Traditional approaches for tissue segmentation used in brain segmentation software packages such as FSL (Jenkinson et al., 2012), SPM (Ashburner and Friston, 2000) or NiftySeg (Cardoso et al., 2015) are based on subject-specific optimisation. FSL and SPM fit a Gaussian Mixture Model to the MR intensities using either a Markov Random Field (MRF) or tissue prior probability maps as regularisation. Alternatively, multi-atlas methods rely on label propagation and fusion from multiple fully-annotated images, i.e. atlases, to the target image (Iglesias and Sabuncu, 2015). These methods typically require extensive pre-processing, e.g. skull stripping, correction of bias field or registration. They are also often time-consuming, and are inherently only adapted for brains devoid of large anatomical changes induced by pathology. Indeed, it has been showed that the presence of lesions distorts the registration output (Sdika and Pelletier, 2009). Similarly, lesions introduce a bias in the MRF. This leads to a performance degradation in presence of lesions for brain volumes measurement Battaglini et al. (2012) and any subsequent analysis.

While quantitative analysis is expected to play a key role in improving the diagnosis and follow-up evaluations of patients with brain lesions, current tools mostly focus on the lesions themselves. Existing quantitative neuroimaging approaches allow the extraction of imaging biomarkers such as the largest diameter, volume, and count of the lesions. Thus, automatic segmentation of the lesions promises to speed up and improve the clinical decision-making process but more refined analysis would be feasible from tissue classification and region parcellation. As such, although very few works have addressed this problem yet, a joint model for lesion and tissue segmentation is expected to bring significant clinical impact.

Deep Neural Networks (DNNs) became the state-of-the-art for most of the segmentation tasks and one would now expect to train a joint lesion and tissue segmentation algorithm. Yet, DNNs require a large amount of annotated data to be successful. Existing annotated databases are usually task-specific, i.e. providing either scans with brain tissue annotations for patients/controls devoid of large pathology-induced anatomical changes, or lesion scans with only lesion annotations. For this reason, the imaging protocol used for the acquisition also typically differs from one dataset to another. Indeed, tissue segmentation is usually performed on T1 scans, unlike lesion segmentation which normally also encompasses Flair Barkhof and Scheltens (2002). Similarly, the resolution and contrast among databases may also vary. Given the large amount of resources, time and expertise required to annotate medical images, given the varying imaging requirement to support each individual task and given the availability of task-specific databases, it is unlikely that large databases for every joint problem, such as lesion and tissue segmentation, will become available for research purposes. There is thus a need to exploit task-specific databases to address joint problems. Learning a joint model from task-specific hetero-modal datasets is nonetheless challenging. This problems lies at the intersection of Multi-Task Learning, Domain Adaptation and Weakly Supervised Learning with idiosyncrasies making individual methods from these underpinning fields insufficient to address it completely.

Multi-Task Learning (MTL) aims to perform several tasks simultaneously by extracting some form of common knowledge or representation and introducing a task-specific back-end. When relying on DNN for MTL, usually the first layers of the network are shared, while the top layers are task-specific. The global loss function is often a sum of task-specific loss functions with manually tuned weights. Recently, Kendall and Gal (2017) proposed a Bayesian parameter-free method to estimate the MTL loss weights and Bragman et al. (2018) extended it to spatially adaptive task weighting and applied it to medical imaging. In addition to arguably subtle differences between MTL and joint learning discussed in more depth later, MTL approaches do not provide any mechanism for dealing with hetero-modal datasets and changes in imaging characteristics across task-specific databases.

Domain Adaptation (DA) is a solution for dealing with heterogeneous datasets. The main idea is to create a common feature space for the two sets of scans. Csurka (2017) proposed an extensive comparison of these methods in deep learning. Learning from hetero-modal datasets could be consider as a particular case of DA. Havaei et al. (2016) proposed a network architecture for dealing with missing modalities. However, DA methods focus on solving a single task and rely on either fully-supervised approaches or unsupervised adaptation as done by Kamnitsas et al. (2017).

Weakly-supervised Learning (WSL) deals with missing, inaccurate, or inexact annotations. Our problem is a particular case of learning with missing labels since each dataset provide a set of labels and the two sets are disjoint. Li and Hoiem (2017) proposed a method to learn a new task from a model trained on another task. This method combines DA through transfer learning and MTL. At the end, two models are created: one for the first task and one for the second one. Kim et al. (2018) extent this approach by using a knowledge distillation loss in order to create a unique joint model. This aims to alternatively learn one task without forgetting the other one. The WSL problem was thus decomposed into a MTL problem with similar limitations for our specific use case.

The contributions of this work are four-fold. First we propose a joint model that performs tissue and lesion segmentation as a unique joint task and thus exploits the interdependence between lesion and tissue segmentation tasks. Starting from a variational formulation of the joint problem, we exploit the disjointness of the label sets to propose a practical decomposition of the joint loss. Secondly, we introduce feature channel averaging across modalities to adapt existing networks for our hetero-modal problem. Thirdly, we develop a new method to minimise the expected risk under the constraint of missing modalities. Relying on reasonable assumptions, we show that the expected risk can be further decomposed and minimised via a tractable upper bound. To our knowledge, no such optimisation method for missing modalities in deep learning has been published before. Finally, we evaluate our framework for white matter lesions and tissue segmentation. We demonstrate that our joint approach can achieve, for each individual task, similar performance compared to a task-specific baseline. Albeit relying on different annotation protocols, results using a small fully-annotated joint dataset demonstrate efficient generalisability.

2 Tissue and lesion segmentation as a single task

In order to develop a joint model, we propose a mathematical variational formulation of the problem and a method to optimise it empirically.

2.1 Formal problem statement

Let $\textbf{x}=(x^{1},..,x^{M})\in\mathcal{X}=\mathbb{R}^{N\times M}$ be a vectorized multimodal image and $y\in\mathcal{Y}=\{0,..,C\}^{N}$ its associated segmentation map. $N$, $M$ and $C$ are respectively the number of voxels, modalities and classes. Our goal is to determine a predictive function $h_{\theta}:\mathcal{X}\mapsto\mathcal{Y}$ that minimises the discrepancy between the ground truth label vector $y$ and the prediction $h_{\theta}(\textbf{x})$. Let $\mathcal{L}$ be a loss function that computes this discrepancy. Following the formalism used by Bottou et al. (2018), given a probability distribution $\mathcal{D}$ over $(\mathcal{X},\mathcal{Y})$ and random variables $(X,Y)$ under this distribution, we want to find $\theta^{*}$ such that:

[TABLE]

Let $\mathcal{C}_{t}$, $\mathcal{C}_{l}$ and [math] be respectively the tissue classes, the lesion classes and the background class. Since $\mathcal{C}_{t}$ and $\mathcal{C}_{l}$ are disjoint, the segmentation map $y$ can be decomposed into two segmentation maps $y_{i}=y^{l}_{i}+y^{t}_{i}$ with $y^{t}_{i}\in\mathcal{C}_{t}\cup\{0\}$, $y^{l}_{i}\in\mathcal{C}_{l}\cup\{0\}$, as shown in \figurereffig:decomposition.

Let’s assume that the loss function $\mathcal{L}$ can also be decomposed into a tissue loss function $\mathcal{L}^{t}$ and a lesion loss function $\mathcal{L}^{l}$. This is common for multi-class segmentation loss functions in particular for those with one-versus-all strategies (e.g. Dice loss, Jaccard loss):

[TABLE]

Then, \equationrefeq:0 can be rewritten as:

[TABLE]

2.2 On the distribution \texorpdfstring$\mathcal{D}$D in the context of heterogeneous databases

As we expect different distributions across heterogeneous databases, two probability distributions of $(X,Y)$ over $(\mathcal{X},\mathcal{Y})$ can be distinguished: 1/ under $\mathcal{D}_{control}$, $(X,Y)$ corresponds to a multimodal scan and segmentation map of a patient without lesions. Note that although we use the term control for convenience, we expect to observe pathology with ”diffuse” anatomical impact, e.g. from dementia. 2/ under $\mathcal{D}_{lesion}$, $(X,Y)$ corresponds to a multimodal scan and segmentation map of a patient with lesions.

Since traditional methods are not adapted in the presence of lesions, the most important and challenging distribution $\mathcal{D}$ to address is the one for patients with lesions, $\mathcal{D}_{lesion}$. In the remainder of this work we thus assume that:

[TABLE]

2.3 Hetero-modal network architecture

In order to learn from hetero-modal datasets, we need a network architecture that allows for missing modalities. We proposed an architecture inspired by HeMIS (Havaei et al., 2016) and HighResNet (Li et al., 2017) shown in \figurereffig:network. Features of each modality are first extracted separately and are then averaged. The spatial resolution of the input and the output are the same. Dilated convolutions and residual connections are used to capture information at multiple scales and avoid the problem of vanishing gradients. This network with weights $\theta$ is used to capture the predictive function $h_{\theta}$.

2.4 Upper-bound for the tissue expected risk \texorpdfstring$\mathcal{R}^{t}$

Although thanks to its hetero-modal architecture, $h_{\theta}$ may now handle inputs with varying number of modalities, the current decomposition of \eqrefeq:0 assumes that all the modalities of $X$ are available for evaluating the loss. In our scenario, we have only access to T1 control scans with tissue annotations or T1 and Flair scans with only lesion annotations. Consequently, as we do not have any T1 and Flair images with tissue annotations, and as evaluating a loss with missing modalities would lead to a bias, estimating $\mathcal{R}^{t}$ is not straightforward.

In this section we propose an upper-bound of $\mathcal{R}^{t}$ using T1 control images with tissue annotations and outputs from the network. Let’s assume that the loss function $\mathcal{L}^{t}$ satisfies the triangle inequality (e.g. Jaccard loss):

[TABLE]

Let $p$ denote the projection of x (will all the modalities) to the T1 modality, $p:\textbf{x}=(x^{\textrm{T1}},x^{\textrm{Flair}})\mapsto x^{\textrm{T1}}$. Under \eqrefeq:H3, $\mathcal{L}^{t}$ satisfies the following inequality:

[TABLE]

In combination with \eqrefeq:H2, this leads to:

[TABLE]

The decomposition in \eqrefeq:upbnd_init requires comparison of inference done from T1 inputs, i.e. $h_{\theta}(p(X))$ with ground truth tissue maps $Y^{t}$. While this provides a step towards a practical evaluation of $\mathcal{R}^{t}$, we still face the challenge of not having tissue annotations $Y^{t}$ under $\mathcal{D}_{lesion}$. Let us further assume that the restriction of the distributions $\mathcal{D}_{lesion}$ and $\mathcal{D}_{control}$ to the parts of the brain not affected by lesions are the same, i.e.:

[TABLE]

By combining \eqrefeq:H3 and \eqrefeq:H4, an upper bound of $\mathcal{R}^{t}$ can be provided as:

[TABLE]

As observed in \eqrefeq:upper, the upper-bound is the sum of the expected loss between the T1 scan outputs and the labels and the expected loss between the outputs using either one or two modalities as input. We emphasise that, to the best of our knowledge, this second loss term does not appear in existing heteromodal approaches such as HeMIS (Havaei et al., 2016).

2.5 Empirical estimation of the decomposed loss

As is the norm in data-driven learning, we do not have access to the true joint probabilities $\mathcal{D}_{control}$ or $\mathcal{D}_{lesion}$. The common method is to estimate the expected risk using training samples. In our case, we have two hetero-modal training samples $\mathcal{S}_{control}$ and $\mathcal{S}_{lesion}$ with respectively tissue and lesion annotations. We can estimate the expected risks $\mathcal{R}^{l}(\theta)$, $\mathcal{R}^{t}_{1}(\theta)$, $\mathcal{R}^{t}_{2}(\theta)$ by respectively using lesion segmentation outputs of lesion T1+Flair scans, tissue segmentation outputs from T1 and T1+Flair scans and tissue segmentation outputs of control T1 scans. \figurereffig:optimisation illustrates the complete training procedure.

3 Experiments

While focusing on the white matter lesion and tissue segmentation problem, our goal in the following experiments is to predict six tissue classes (white matter, gray matter, basal ganglia, ventricles, cerebellum, brainstem), the white matter lesions and the background.

3.1 Data

To demonstrate the feasibility of our joint learning approach, we used three sets of data.

Lesion data $S_{lesion}$:

The White Matter Hyperintensities (WMH) database consists of 60 sets of brain MR images (T1 and Flair, $M=2$) with manual annotations of WMH (\urlhttp://wmh.isi.uu.nl/). The data comes from three different institutes.

Tissue data $S_{control}$:

Neuromorphometrics provided 32 T1 scans ($M^{{}^{\prime}}=1)$ for MICCAI 2012 with manual annotations of 155 structures of the brain from which we deduct the six tissue classes. In order to have balance training datasets for the two types of segmentation, and similar to Li et al. (2017), we added 28 T1 control scans from the ADNI2 dataset with bronze standard parcellation of the brain structures computed with the accurate but time-consuming algorithm of Cardoso et al. (2015).

Fully annotated data:

MRBrainS18 (\urlhttp://mrbrains18.isi.uu.nl/) is composed of 30 sets of brain MR images with tissue and lesions manual annotations. Only 7 MR images are publicly available. We used this data only for evaluation and not for training. To be consistent with the lesion data, the cerebrospinal fluid is considered as background.

To satisfy the assumption \eqrefeq:H4, we resampled the data to $1\times 1\times 3$ mm3, used a histogram-based scale Milletari et al. (2016) and a zero-mean unit-variance normalization.

3.2 Choice of the loss

We used the probabilistic version of the Jaccard loss for $\mathcal{L}$:

[TABLE]

\eqref

eq:H1 is satisfied because of the one-versus-all strategy, i.e. sum over the classes of a class-specific loss. In order to give the same weight to the lesion segmentation and the tissue segmentation, we choose for any tissue class $c$, $w_{c}=\frac{1}{16}$ and for the lesion class $l$, $w_{l}=\frac{1}{2}$. While the triangle inequality holds for the Jaccard distance Kosub (2018), the proof that its probabilistic version also satisfies it, i.e. \eqrefeq:H3, is left for future work.

3.3 Implementation details

We implemented our network in NiftyNet, a Tensorflow-based open-source platform for deep learning in medical imaging Gibson et al. (2018). Convolutional layers are initialised such as He et al. (2015). The scaling and shifting parameters in the batch normalisation layers were initialised to 1 and 0 respectively. As suggested by Ulyanov et al. (2016), we used instance normalization for inference. We used the Adam optimisation method (Kingma and Ba, 2014). The learning rate $l_{R}$, $\beta_{1}$, $\beta_{2}$ were respectively set up to 0.005, 0.9 and 0.999. At each training iteration, we feed the network with one image from the tissue dataset and one from the lesion dataset. $120\times 120\times 40$ sub-volumes were randomly sampled from the training data using an uniform sampling for the tissue data and a weighted sampling based on dilated lesions maps for the lesion data. The models were trained until we observed a plateau in performance on the validation set. We experimentally found that the inter-modality loss has to be skipped for the first (5000) iterations. We randomly spitted the data into $70\%$ for training, $10\%$ for validation and $20\%$ for testing for each of the 4 folds.

3.4 Results for the joint learning model

Joint learning versus single task learning

First, we compare individual models to the joint model using our approach. The lesion segmentation (W) model was trained on WMH dataset with the lesion annotations, the tissue segmentation (N) on Neuromorphometrics dataset with the tissue annotations, and our method (W+N) on WMH and Neuromorphometric datasets with their respective set of annotations. The similarity between the prediction and the ground truth is computed using the Dice Similarity Coefficient (DSC) for each class. \tablereftab:scores and \figurereffig:results show the results of these models on test images. The joint model and single task models achieve comparable performance. This suggest that learning from hetero-modal datasets via our method does not degrade the task-specific performance. Moreover, we observe in \figurereffig:results that the tissue knowledge learnt from T1 scans has been well generalised to multi-modal scans.

Joint model versus fully-supervised model

In this section, we compare our method (W+N) to the fully-supervised (M) model trained on MRBrainS18 using both tissue and lesion annotations. We evaluated the performance on the three datasets. On the one hand, we submitted our models to the online challenge MRBrainS18. One of the major benefits of evaluating our method on a challenge is to directly benchmark our method with existing methods, in particular with traditional methods such as SPM (Ashburner and Friston, 2000). On the other hand, we compared the performance on the tissue and lesion datasets using either all the scans (M) or the testing split (W+N). The DSC was computed for each class and \tablereftab:scores show the results. First our model outperforms SPM on 6 of the 7 classes. Secondly, the two models achieve very similar performance on lesion segmentation. Concerning the tissue segmentation, as expected, each of the network outperforms on its training datasets. However, the fully supervised model doesn’t generalise to Neuromorphometrics dataset. In contrast, the differences are smaller for the tissue segmentation classes on MRBrainS18. Especially, \figurereffig:differences shows differences in the annotation protocol between MRBrainS18 and Neuromorphometrics for the white matter, brainstem and cerebellum and how it affects the predictions.

4 Conclusion

We propose a joint model learned from hetero-modal datasets with disjoint heterogeneous annotations. Our approach is mathematically grounded, conceptually simple, new and relies on reasonable assumptions. We validated our approach by comparing our joint model with single-task learning models. We show that similar performance can be achieved for the tissue segmentation and lesion segmentation in comparison to task-specific baselines. Moreover, our model achieves comparable performance to a model trained on a small fully-annotated joint dataset. Our work shows that the knowledge learnt from one modality is preserved when more modalities are used as input. In the future, we will evaluate our approach on datasets with annotations protocols showing less variability. Furthermore, exploitation of recent techniques for domain adaptation could help us bridge the gap and improve the performance by helping to better satisfy some of our assumptions. Finally, we plan to integrate uncertainty measures in our framework as a future work. As one of the first work to methodologically address the problem of joint learning from hetero-modal datasets, we believe that our approach will help DNN make further impact in clinical scenarios.

\midlacknowledgments

This work was supported by the Wellcome Trust [203148/Z/16/Z] and the Engineering and Physical Sciences Research Council (EPSRC) [NS/A000049/1]. TV is supported by a Medtronic / Royal Academy of Engineering Research Chair [$RCSRF1819\backslash 7\backslash 34$].