Anatomical Priors in Convolutional Networks for Unsupervised Biomedical Segmentation

TL;DR

This paper introduces a generative probabilistic model using anatomical priors in convolutional neural networks to enable fast, unsupervised biomedical image segmentation without paired training data, demonstrated on brain MRI datasets.

Contribution

It presents a novel approach that incorporates anatomical priors into CNNs for unsupervised segmentation, applicable across different datasets and imaging modalities.

Findings

Enables fast unsupervised segmentation with anatomical priors.

Effective on large multi-study brain MRI dataset.

Facilitates segmentation in data-scarce clinical scenarios.

Abstract

We consider the problem of segmenting a biomedical image into anatomical regions of interest. We specifically address the frequent scenario where we have no paired training data that contains images and their manual segmentations. Instead, we employ unpaired segmentation images to build an anatomical prior. Critically these segmentations can be derived from imaging data from a different dataset and imaging modality than the current task. We introduce a generative probabilistic model that employs the learned prior through a convolutional neural network to compute segmentations in an unsupervised setting. We conducted an empirical analysis of the proposed approach in the context of structural brain MRI segmentation, using a multi-study dataset of more than 14,000 scans. Our results show that an anatomical prior can enable fast unsupervised segmentation which is typically not possible using standard convolutional networks. The integration of anatomical priors can facilitate CNN-based anatomical segmentation in a range of novel clinical problems, where few or no annotations are available and thus standard networks are not trainable. The code is freely available at http://github.com/adalca/neuron.

Equations

$$p(\mathbf{z})=\mathcal{N}(\mathbf{z};0,\text{mathbbm}1),$$

$$p_{{\mathbf{\theta }}_{s|z}}(\mathbf{s}|\mathbf{z})$$

$$p_{{\mathbf{\theta }}_{x|s}}(\mathbf{x}|\mathbf{s})=\prod _j\prod _l\mathcal{N}(\mathbf{x}[j];{\mathbf{\mu }}_l,{\sigma }_l{)}^{\delta (\mathbf{s}[j]=l)},$$

$${\hat{\mathbf{s}}}_i$$

$$=\arg \underset{{\mathbf{s}}_i}{\max }\log p({\mathbf{s}}_i,{\mathbf{x}}_i;\mathbf{\theta })$$

$$\left[q_{\mathbf{\varphi }}(\mathbf{z}|\mathbf{x},\mathbf{s})||p_{\mathbf{\theta }}(\mathbf{z}|\mathbf{x},\mathbf{s})\right]$$

$$={\mathrm{I}\mathrm{E}}_q\left[\log q_{\mathbf{\varphi }}(\mathbf{z}|\mathbf{x},\mathbf{s})-\log p_{\mathbf{\theta }}(\mathbf{z}|\mathbf{x},\mathbf{s})\right]$$

$$={\mathrm{I}\mathrm{E}}_q\left[\log q_{\mathbf{\varphi }}(\mathbf{z}|\mathbf{x},\mathbf{s})-\log p_{\mathbf{\theta }}(\mathbf{x},\mathbf{s},\mathbf{z})\right]+\log p_{\mathbf{\theta }}(\mathbf{x},\mathbf{s}).$$

$$\log p(\mathbf{x},\mathbf{s})$$

$$+{\mathrm{I}\mathrm{E}}_q\left[\log p_{\mathbf{\theta }}(\mathbf{x},\mathbf{s},\mathbf{z})-\log q_{\mathbf{\varphi }}(\mathbf{z}|\mathbf{x},\mathbf{s})\right].$$

$${\mathcal{V}}_{\text{model}}($$

$$={\mathrm{I}\mathrm{E}}_q\left[\log p_{\mathbf{\theta }}(\mathbf{x},\mathbf{s},\mathbf{z})-\log q_{\mathbf{\varphi }}(\mathbf{z}|\mathbf{x},\mathbf{s})\right]$$

$$={\mathrm{I}\mathrm{E}}_q\left[\log p_{\mathbf{\theta }}(\mathbf{x},\mathbf{s}|\mathbf{z})\right]-\text{KL}\left[q_{\mathbf{\varphi }}(\mathbf{z}|\mathbf{x},\mathbf{s})||p(\mathbf{z})\right].$$

$$q_{\mathbf{\varphi }}(\mathbf{z}|\mathbf{x},\mathbf{s})$$

$$=\mathcal{N}(\mathbf{z};{\mathbf{\mu }}_{z|x},{\mathbf{\Sigma }}_{z|x}).$$

$$q_{\mathbf{\psi }}(\mathbf{z}|\mathbf{s})$$

$${\mathcal{V}}_{\text{prior}}($$

$$={\mathrm{I}\mathrm{E}}_q\left[\log p_{\mathbf{\theta }}(\mathbf{s},\mathbf{z})-\log q_{\mathbf{\varphi }}(\mathbf{z}|\mathbf{s})\right]$$

$$={\mathrm{I}\mathrm{E}}_q\left[\log p_{\mathbf{\theta }}(\mathbf{s}|\mathbf{z})\right]-\text{KL}\left[q_{\mathbf{\varphi }}(\mathbf{z}|\mathbf{s})||p(\mathbf{z})\right].$$

$${\mathcal{L}}_{\text{prior}}($$

$$=\text{KL}\left[\log q_{\mathbf{\psi }}(\mathbf{z}|{\mathbf{s}}_i)||\log p(\mathbf{z})\right]-\log p_{{\mathbf{\theta }}_{s|z}}({\mathbf{s}}_i|{\mathbf{z}}_k)$$

$$=\frac{1}{2}\sum _j\left(1+\log {\mathbf{\Sigma }}_{z|s_i}[j]-{\mathbf{\mu }}_{z|s_i}^2[j]-{\mathbf{\Sigma }}_{z|s_i}[j]\right)$$

$$-\sum _j{\mathbf{s}}_i[j]\log f({\mathbf{z}}_k;{\mathbf{\theta }}_{s|z})[j].$$

$$\mathcal{L}$$

$$=\text{KL}\left[q_{\mathbf{\varphi }}(\mathbf{z}|{\mathbf{x}}_i)||p(\mathbf{z})\right]-\log p_{{\mathbf{\theta }}_{s|z}}({\mathbf{x}}_i,{\mathbf{s}}_i|{\mathbf{z}}_k).$$

$$=\text{KL}\left[q_{\mathbf{\varphi }}(\mathbf{z}|{\mathbf{x}}_i)||p(\mathbf{z})\right]-\log p_{{\mathbf{\theta }}_{s|z}}({\mathbf{s}}_i|{\mathbf{z}}_k)$$

$$-\log p_{{\mathbf{\theta }}_{x|s}}({\mathbf{x}}_i|{\mathbf{s}}_i).$$

$$=\frac{1}{2}\sum _j\left(1+\log {\mathbf{\Sigma }}_{z|x_i}[j]-{\mathbf{\mu }}_{z|x_i}^2[j]-{\mathbf{\Sigma }}_{z|x_i}[j]\right)$$

$$-\sum _j{\mathbf{s}}_i[j]\log f(\mathbf{z};{\mathbf{\theta }}_{s|z})[j]$$

$$+\sum _j\sum _l\frac{\delta ({\mathbf{s}}_i[j]=l)}{2{\sigma }_l^2}({\mathbf{x}}_i-{\mathbf{\mu }}_l),$$

$${\mathrm{I}\mathrm{E}}_q$$

$$={\mathrm{I}\mathrm{E}}_q\left[\log {\int }_{\mathbf{s}}{p}_{\mathbf{\theta }}(\mathbf{x}|\mathbf{s})p_{{\mathbf{\theta }}_{s|z}}(\mathbf{s}|\mathbf{z})d\mathbf{s}\right]$$

$$\ge {\mathrm{I}\mathrm{E}}_q\left[{\int }_{\mathbf{s}}{p}_{{\mathbf{\theta }}_{s|z}}(\mathbf{s}|\mathbf{z})\log p_{\mathbf{\theta }}(\mathbf{x}|\mathbf{s})d\mathbf{s}\right]$$

$$={\mathrm{I}\mathrm{E}}_{q,p(\mathbf{s}|\mathbf{z})}\left[\log p_{\mathbf{\theta }}(\mathbf{x}|\mathbf{s})\right]$$

$${\mathcal{L}}_{\text{model}}($$

$$=\frac{1}{2}\sum _j\left(1+\log {\mathbf{\Sigma }}_{z|x_i}[j]-{\mathbf{\mu }}_{z|x_i}^2[j]-{\mathbf{\Sigma }}_{z|x_i}[j]\right)$$

$$+\sum _j\sum _l\frac{f_{j,l}({\mathbf{z}}_k|{\theta }_{s|z})}{2{\sigma }_l^2}({\mathbf{x}}_i-{\mathbf{\mu }}_l),$$

$$H(\mathbf{s}[j])$$

$$=-\sum _{l}p(\delta (\mathbf{s}[j]=l)|{\mathbf{x}}_i)\log (p(\delta (\mathbf{s}[j]=l)|{\mathbf{x}}_i)).$$

$${\hat{\mathbf{s}}}_p={\mathbf{f}}_L\circ \dots \circ {\mathbf{f}}_1(\mathbf{x}),$$

$$\text{Dice}(\hat{\mathbf{s}},{\mathbf{s}}_t)=\frac{2{\sum }_j\delta (\hat{\mathbf{s}}[j]=l)\delta ({\mathbf{s}}_t[j]=l)}{{\sum }_j\delta (\hat{\mathbf{s}}[j]=l)+{\sum }_j\delta ({\mathbf{s}}_t[j]=l)}.$$

Full text

1 Introduction

Biomedical image segmentation plays a crucial role in many applications, such as population analysis, disease progression modelling, or treatment planning. Convolutional neural networks (CNNs), a class of deep learning methods, have been employed to derive powerful biomedical segmentation algorithms, showing promise of overcoming limitations in previous methods [3, 4, 29, 34]. However, CNN-based approaches most often depend on (large-scale) training data, particularly in the form of image scans paired with segmentations. These annotations are often costly and challenging to obtain because they require the tedious effort of a trained expert, taking several expert hours per scan.

1.1 Contributions

To our knowledge, there has not been a theoretically rigorous effort to integrate rich probabilistic anatomical priors with a CNN-based segmentation model in a computationally effective manner. We introduce a generative model for biomedical segmentation that employs a deep anatomical prior. We describe a principled derivation that follows directly from our generative model. We demonstrate that this yields intuitive cost functions and simpler models. We use an auto-encoding variational CNN to characterize the anatomical prior, and an encoder-decoder CNN to provide fast segmentation of medical images in unsupervised settings.

We demonstrate the method in an unsupervised biomedical image segmentation setting where paired annotations are not available. Our proposed strategy is general and computationally efficient, provides a natural framework for sampling possible subject-specific segmentations of a scan, and provides uncertainty estimates for these segmentations.

2 Related Work

2.1 Segmentation Convolutional Neural Networks

CNN-based segmentation approaches generally rely on fully convolutional architectures applied to image data. They extract hierarchical and multi-resolution features that are in turn combined to compute a semantic segmentation [23, 29, 31, 34].

A popular discriminative segmentation architecture, U-net [29], involves a convolutional encoder or downsampling network, followed by a convolutional decoder or upsampling network, and skip-connections between layers. The encoder captures relevant features of the input image at different resolutions. The decoder then synthesizes a high-resolution segmentation, using the skip connections to achieve voxel-level precision. While the exact architecture of these networks, such as the number of layers and levels, size of convolution kernels, or application of batch normalization vary, they typically involve millions of parameters and necessitate large datasets and data augmentation techniques to train.

CNN-based segmentation models have two major shortcomings: the dependency on annotated data, limiting their use in unsupervised settings; and their lack of anatomical knowledge. The latter limits the network’s ability to be faithful to known anatomical shapes during segmentation.

In our work, we use CNN architectures to learn anatomical priors and segment medical images. The prior eliminates the burden of providing paired example segmentations.

2.2 Priors for Convolutional Neural Networks

A clinical expert performing manual delineation relies on spatial coordinates and prior knowledge about anatomy, and may use a template of the structures to constrain the task. This process draws on the anatomical similarity across patient scans. This is in stark contrast with typical computer vision problems that have led to many popular CNN architectures, where object location, shape, and appearance can be unpredictable.

Convolutional methods are often limited in incorporating domain expertise. For example, U-Net [29] and its derivatives produce segmentation algorithms that do not exploit location information or other explicit anatomical priors. A CNN might have difficulty differentiating two distinct objects that are consistently in two specific parts of the scan, if they have the same intensity and context (as in bilateral structures in two hemispheres)111Assuming that the field of view of the network is constrained to not include the other object’s vicinity. While increasingly more complex networks that extend receptive fields may tease out object differences in supervised settings, the problem would be trivial if we consider anatomical knowledge like spatial location. Furthermore, in these modalities, image contrast can be weak or noisy in certain regions resulting in uncertainty of the segmentations. An anatomical prior can resolve these ambiguities, while making the segmentation task easier.

A popular strategy to explicitly employ prior structure in CNNs for biomedical image segmentation is to use a conditional random field (CRF) as a post-processing step [13, 30, 34]. However, CRFs only capture local constraints, and add to the computational burden. Location information has been included as a feature in patch-based CNN segmentation networks [34]. While this addition carries prior location information, it is network-specific, increases the parameter burden on the network, and does not capture shape information.

Recent methods have employed shape priors for neural network solutions in supervised problems [26, 28]. In particular, they often design a series of networks that learn representations of images and segmentations in a supervised setting. They propose ad-hoc cost functions that encourage the computed segmentations to be similar to both the learned shape and the ground truth. These methods attempt to correct segmentations produced by standard CNNs by adding a prior constraint.

Convolutional image generative models, such as generative adversarial nets, have grown in popularity. They have recently been applied to biomedical image segmentations [16, 24] in a supervised setting where standard loss functions are combined with adversarial losses. A series of recent papers in the computer vision community removes the requirement for paired data by introducing a cycle dependency [37]. However, these methods are less applicable in medical image segmentation with many anatomical labels, as an image signal can pass through the rich networks at low cost, leading to a perfect cycle loss, circumventing the required constraints [37].

Variational Bayes auto-encoders have been used for various tasks to learn probabilistic generative models, and often use convolutional networks [18]. Our method builds on these models to combine anatomical priors with image generation.

2.3 Classical Generative Models

Encoding and exploiting prior knowledge is common in generative models. Our inspiration comes from classical atlas-based probabilistic segmentation methods that estimate the maximum a posteriori (MAP) probability based on a generative model involving a prior probability and likelihood  [5, 10, 17, 27, 32, 33, 35].

The prior term captures knowledge of underlying anatomy and usually involves a probabilistic atlas and a spatial deformation that models geometric variation. The spatial deformation can be explicitly solved using a registration algorithm or accounted for in a unified segmentation framework [1].

The likelihood models the physical process that yields medical image intensities, sometimes called the appearance model, conditioned on the latent anatomy. These appearance models are often simpler, relying on additive and/or multiplicative Gaussian or Rician noise models [35]. Model parameters are most often estimated using training data, such as annotated image pairs, for example using maximum likelihood.

Given a new image, most popular segmentation algorithms use numerical non-convex optimization and can take several hours per image on a modern CPU.

In our model, we draw on ideas from classical model-based biomedical segmentation algorithms, convolutional neural networks (CNNs) used in semantic segmentation, and recent developments in variational Bayes approximations using neural networks. In our experiments, we consider the segmentation of structural brain MRI scans into cortical and subcortical regions of interest (ROIs). Our results show that the proposed anatomical prior enables rapid unsupervised segmentation. While complex, specialized tools exist for segmenting some specific scan modalities or particular diseases, they do not generalize to other modalities and can take hours to process one scan. Our goal is to provide a first general approach to biomedical image segmentation in an unsupervised setting.

3 Generative Model

We let $\boldsymbol{x}$ be an (MR) 3D volume, and assume it is generated from a 3D anatomical segmentation map $\boldsymbol{s}$. We will use $x[j]$ and $s[j]$ to denote the image intensity and label at voxel $j$, respectively.

We use a generative model to describe the spatial distribution, shape, and appearance of anatomical structures. Figure 1 provides a graphical representation.

The prior captures our knowledge about spatial distributions and shape of anatomy. We let $\boldsymbol{z}$ be a latent variable representing an embedding of these shapes, and model the prior probability of this embedding as normal with mean $\boldsymbol{0}$ and an identity covariance matrix:

[TABLE]

where $\mathcal{N}(\cdot;\boldsymbol{\mu},\boldsymbol{\Sigma})$ is the normal distribution parametrized by mean $\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$.

We let $\boldsymbol{s}$ be drawn from a categorical prior distribution determined by the low-dimensional embedding $\boldsymbol{z}$ via $p_{\boldsymbol{\theta}_{s|z}}(\boldsymbol{s}|\boldsymbol{z})$:

[TABLE]

where $f_{j,l}(\cdot;\boldsymbol{\theta}_{s|z})$ is the probability of label $l$ at voxel $j$.

Finally, given the label map $\boldsymbol{s}$, the intensity observations are generated via $p_{\boldsymbol{\theta}_{x|s}}(\boldsymbol{x}|\boldsymbol{s})$, sampled at each voxel from a normal distribution:

[TABLE]

where $\boldsymbol{\theta}_{z|s}=\{\boldsymbol{\mu}_{l},\sigma_{l}\}$, and $\delta(\boldsymbol{s}[j]=l)$ is the indicator function that evaluates to 1 if $\boldsymbol{s}[j]=l$ and 0 otherwise. The joint likelihood is therefore $p_{\boldsymbol{\theta}}(\boldsymbol{x},\boldsymbol{s}|\boldsymbol{z})=p_{\boldsymbol{\theta}_{x|s}}(\boldsymbol{x}|\boldsymbol{s})p_{\boldsymbol{\theta}_{s|z}}(\boldsymbol{s}|\boldsymbol{z})$, where $\boldsymbol{\theta}=\{\boldsymbol{\theta}_{s|z},\boldsymbol{\theta}_{x|s}\}$. Intuitively, the embedding $\boldsymbol{z}$ determines the possible anatomical shapes in $\boldsymbol{s}$, which in turn determine the possible observed images $\boldsymbol{x}$.

We describe the learning procedure in the next section. Given learned parameters, to obtain the segmentation $\boldsymbol{s}_{i}$ given a new image $\boldsymbol{x}_{i}$, we perform MAP estimation:

[TABLE]

4 Learning

In this section, we describe a learning strategy that uses convolutional neural networks to estimate anatomical representations and optimize posterior segmentation distributions. This procedure is applicable to broad modelling choices for the probability distributions described above. We also discuss a separate learning procedure for the anatomical prior, uncertainty estimation, and implementation.

Without assuming voxel independence of the segmentation map given an image, estimating the posterior probability $p_{\boldsymbol{\theta}}(\boldsymbol{s}|\boldsymbol{x})$ is intractable since it involves integrating over the latent variable $\boldsymbol{z}$. Estimating $p_{\boldsymbol{\theta}}(\boldsymbol{z}|\boldsymbol{x},\boldsymbol{s})$ is similarly intractable, making the Expectation Maximization algorithm not pertinent.

We first introduce an encoding probability $q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x},\boldsymbol{s})$ as an approximation to the intractable $p_{\boldsymbol{\theta}}(\boldsymbol{z}|\boldsymbol{x},\boldsymbol{s})$, similar to [18]. Consider the KL divergence between the approximate distribution $q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x},\boldsymbol{s})$ and the true posterior $p_{\boldsymbol{\theta}}(\boldsymbol{z}|\boldsymbol{x},\boldsymbol{s})$:

[TABLE]

Rearranging terms, we obtain

[TABLE]

Since the KL divergence of the approximate and true posterior of $\boldsymbol{z}$ is non-negative, the second term is referred to as the variational lower bound of the model evidence or joint probability. For a given approximate distribution $q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x},\boldsymbol{s})$, we can estimate $\boldsymbol{\theta}$ by optimizing the lower bound:

[TABLE]

We model the approximating posterior $q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x},\boldsymbol{s})$ as a normal that depends on the image only:

[TABLE]

where $\boldsymbol{\Sigma}_{z|x}$ is diagonal.

We estimate the parameters of the approximating distribution using convolutional neural networks. We design an encoding convolutional neural network $\text{enc}_{\boldsymbol{\phi}}(\boldsymbol{x})$ that takes as input $\boldsymbol{x}$ and outputs the parameters of the approximating posterior distribution $\boldsymbol{\mu}_{z|x}(\boldsymbol{x})$, and $\boldsymbol{\Sigma}_{z|x}(\boldsymbol{x})$. This network learns how to embed an entire (MR) image into the most likely low-dimensional anatomical embedding $\boldsymbol{z}$ and its variance.

Conditioned on $\boldsymbol{z}$, the probability of the segmentation can be computed with a decoder network $\text{dec}_{\boldsymbol{\theta}_{s|z}}(\boldsymbol{z})$ that takes $\boldsymbol{z}$ as input and outputs the parameters $f(\boldsymbol{z};\boldsymbol{\theta}_{s|z})$ of the segmentation categorical distribution $p_{\boldsymbol{\theta}_{s|z}}(\boldsymbol{s}|\boldsymbol{z})$. The parameters $\boldsymbol{\theta}_{s|z}$ of this decoder can be learned using a separate set of segmentations, as described below.

The final part of the generative model, the appearance or likelihood model, can also be learned with a neural network that takes a segmentation probability map as input and computes the parameters $\boldsymbol{\mu}_{l}$. We separately estimate $\sigma_{l}$, assuming additive zero mean Gaussian noise in an image, using a difference of Laplacian filters [15].

4.1 Auto-Encoding Anatomical Prior

In this work, we learn a prior independently from an unpaired segmentation dataset. This enables the flexibility of having an external description of the anatomy that need not be available in the current data. Unfortunately, as before, estimating the probability distribution $p(\boldsymbol{s})$ is intractable. Following a derivation similar to the previous section and to the auto-encoding variational Bayes framework, we introduce an approximation $q_{\boldsymbol{\psi}}(\boldsymbol{z}|\boldsymbol{s})$ to the posterior $p(\boldsymbol{z}|\boldsymbol{s})$ as a normal distribution:

[TABLE]

where $\boldsymbol{\Sigma}_{z|s}$ is diagonal, leading to the following lower bound:

[TABLE]

This can be optimized using a Stochastic Gradient Variational Bayes (SGVB) estimator that uses mini-batches. The reparametrization trick enables us to sample $\boldsymbol{z}_{k}\sim q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{s})$, leading to an approximation of the expectation ${\rm I\kern-3.00003ptE}\left[\cdot\right]$ [18]. The loss $\mathcal{L}_{i}(\boldsymbol{\theta},\boldsymbol{\phi};\boldsymbol{s}_{i})$ for each data point $\boldsymbol{s}_{i}$ and sample $\boldsymbol{z}_{k}\sim q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{s}_{i})$ is

[TABLE]

We design an encoding network $\text{enc}_{\boldsymbol{\psi}}(\boldsymbol{s})$ that takes a segmentation map as input and outputs the parameters $\boldsymbol{\mu}_{z|s}$ and $\boldsymbol{\Sigma}_{z|s}$. Importantly, we learn the parameters of the encoding network given only a set of segmentations $\{\boldsymbol{s}_{i}\}$, which can be derived from other imaging modalities and/or datasets. The segmentation prior therefore does not require paired training data in the traditional sense. For example, we can use a prior computed using publicly available annotated datasets such as [19] in a problem that involves a different imaging modality than in the current task.

4.2 Unsupervised Learning

We assume we have learned a segmentation prior using the Auto-Encoding Anatomical Prior described in the previous section. In particular, we will utilize the decoder component of the prior model, namely $p_{\boldsymbol{\theta}_{s|z}}(\boldsymbol{s}|\boldsymbol{z})$.

If we had annotated pairs $\{\boldsymbol{x}_{i},\boldsymbol{s}_{i}\}$, we could jointly learn model parameters $\boldsymbol{\theta}_{x|s}$, and variational parameters $\boldsymbol{\phi}$ by optimizing the evidence lower bound objective (7), similar to the previous section. For each sample $\{\boldsymbol{x}_{i},\boldsymbol{s}_{i}\}$ and sample $\boldsymbol{z}_{k}\sim q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x}_{i},\boldsymbol{s}_{i})$, the loss function would be

[TABLE]

resulting in terms of KL divergence, segmentation map categorical cross-entropy, and intensity-based mean squared error, respectively. During training, these terms would ensure that the probability $q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x}_{i})$ stays close to the standard normal, while explaining the segmentations, and that the model parameters $\boldsymbol{\theta}_{x|s}=\{\boldsymbol{\mu}_{l},\boldsymbol{\sigma}_{l}\}$ capture the relationship between the segmentations and the images.

However, in this paper we tackle the unsupervised setting, where annotated pairs $\{\boldsymbol{x}_{i},\boldsymbol{s}_{i}\}$ are not available, and we only have the images $\{\boldsymbol{x}_{i}\}$. Therefore, we cannot compute the categorical cross entropy term in (12). Instead, we marginalize over the segmentation $\boldsymbol{s}$ in the second term of the variational lower bound (7):

[TABLE]

where we used Jensen’s inequality. We therefore arrive at the following upper bound of the loss function:

[TABLE]

where we used the factorization of $p_{\theta{s|z}}(\boldsymbol{s}|\boldsymbol{z})$ over voxels from (2), and sample $\boldsymbol{z}_{k}\sim q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x}_{i})$.

4.3 Inference and uncertainty

Given a new image $\boldsymbol{x}$, we approximate $\widehat{\boldsymbol{s}}=\arg\max_{\boldsymbol{s}}p_{\boldsymbol{\theta}}(\boldsymbol{s}|\boldsymbol{x})$ by first obtaining $\boldsymbol{\mu}_{z}$ using the encoder $\text{enc}_{\boldsymbol{\theta}_{z|x}}(\boldsymbol{x})$, and taking the maximum segmentation at each voxel $\widehat{\boldsymbol{s}}=\arg\max_{\boldsymbol{s}}\text{dec}_{\boldsymbol{\theta}_{s|z}}(\boldsymbol{\mu}_{z})$. The operations are fast, since both are feed-forward neural networks.

This model also enables sampling segmentations conditioned on a particular image and enables estimation of uncertainty. Given an input image $\boldsymbol{x}_{i}$, we can create samples $\boldsymbol{z}_{k}\sim q_{\boldsymbol{\phi}}(\boldsymbol{z}|\boldsymbol{x})$ and $\boldsymbol{s}_{k}\sim p_{\boldsymbol{\theta}_{s|z}}(\boldsymbol{s}|\boldsymbol{z}_{k})$, simulating different plausible segmentations for a given subject. We can estimate the uncertainty of our segmentation given a new image $\boldsymbol{x}_{i}$ using

[TABLE]

4.4 Implementation

A CNN can be seen as a hierarchical function, a set of concatenated functions, or layers. For example, CNNs often map some input image $\boldsymbol{x}$ to an output probability $\hat{\boldsymbol{s}}_{p}$:

[TABLE]

where $\circ$ denotes concatenation, $\boldsymbol{f}_{i}$ is often some nonlinear function such as a rectified linear unit or ReLU or max-pool [12] applied to (linear) convolutions of the output of the previous layer $\boldsymbol{f}_{i-1}$ (with $\boldsymbol{f}_{0}=\boldsymbol{x}$).

Although we operate on 3D images, we use a 2D architecture in our experiments. We experimented with 3D architectures as well, but found little gain while facing significant challenges related to limitations between GPU memory, batch size, the number of features, and the number of labels possible. Each encoder consists of five downsampling levels of one convolution layer each, with 3x3 convolution kernels with elu activations, and 32 features for each kernel. The final layer is dense, with $1000$-long encoding of the means and standard deviations representations.

The decoder is a mirror of this design, but upsamples instead of downsampling and ends with a sigmoid activation. In addition, we use a final layer that implements a pixel-wise spatially-varying voxel-wise (location) prior $p_{\text{loc}}(\boldsymbol{s})$, which is multiplied with dec$(\boldsymbol{z})$ (in practice, we add the logarithms). As is common in the atlas-based segmentation literature, the prior $p_{\text{loc}}(\boldsymbol{s})$ was computed as the frequency of labels in the held out prior dataset, in affine-normalized coordinate system. This layer discourages any extreme decodings of $\boldsymbol{z}$ but does not capture shape properties, which is encoded in dec$(\boldsymbol{z})$.

We implement the normal probability $p(\boldsymbol{x}|\boldsymbol{s})$ with a single linear layer. We also find it useful to pre-train the image encoder using an image variational auto-encoder similar to the segmentation one. The encoder weights are used as initialization only. During training, we used the Adadelta optimizer [36].

For the latent encoding layers representing $\boldsymbol{\mu}_{z}$ and $\boldsymbol{\Sigma}_{z}$, we introduce an activation function that discourages the sample activations from being too large, helping limit numerical issues stemming in sampling from these layers during the reparametrization trick. We use concepts from the softsign and tanh activations to define our function as $\text{act}(\boldsymbol{x})=\text{softsign}_{\alpha}(\boldsymbol{x})\log(2+\alpha*|\boldsymbol{x}|)$.

5 Experiments

We demonstrate our model on two datasets. For the first dataset, we obtain ground truth segmentations using a specialized algorithm with intense computational requirements, combined with manual work and QC [9]. We use a subset of these segmentations to learn the prior probability parameters. We treat the rest of the dataset as unsupervised, where we only use the ground truth segmentations as validation. For the second dataset we do not have ground truth, offering a realistic scenario. Figure 3 shows example images from the two datasets, highlighting the difference, and the difficulty of the task.

5.1 Data

T1w scan dataset

We gathered a large-scale multi-site, multi-study dataset of more than 14,000 T1-weighted brain MRI scans from eight publicly available datasets: including data from ADNI [25], OASIS [20], ABIDE [7], ADHD200 [22], MCIC [11], PPMI [21], HABS [6], and Harvard GSP [14]. Subject age ranges, health states, and acquisition details vary with each dataset, but all scans were resampled to a 256x256x256 grid with 1mm isotropic voxels, and all images cropped to 160x192x224 to eliminate entirely-background voxels.

We carry out standard pre-processing steps, including affine spatial normalization using FreeSurfer for each scan [9]. All MRIs were also segmented with FreeSurfer - a task that takes several CPU hours per scan. We also applied quality control (QC) using visual inspection to catch gross errors in segmentation results.

We partitioned the data into a prior training subset of 5,000 images, where we only used the annotations. The rest of the data was treated as an unannotated dataset, where QCed segmentations were only used for validation.

While developing the network architectures we partitioned the rest of the data into training, validation and test sets. Once the architecture was fixed, we reported results on the test dataset by training and evaluating the model in an unsupervised fashion.

T2-FLAIR scan dataset

We also gathered a dataset of more than 3800 T2-FLAIR scans, a significantly different MR modality, from the ADNI cohort. These scans exhibit significantly different tissue properties compared to the T1w images, lower acquisition quality, and exhibit $5$mm slice spacing (Figure 3). They provide a good test of our hypothesis that priors learned are useful for segmenting image data with different tissue properties. To our knowledge there is no automatic method to obtain detailed anatomical segmentations for these images. We affinely align these images to the same space as the T1 images using mutual information based affine registration with ANTs [2]. We perform brain extraction using an in-house developed neural network-based algorithm that uses a UNet architecture and extensive data augmentation.

In the set of annotations that we used to train the prior, we avoided including any annotations coming from ADNI subjects whose T2-FLAIR scans are in this dataset.

5.2 Evaluation

We evaluate our results both visually and quantitatively. For the T1w dataset, we use a volume overlap measure, Dice, to quantify the automatic segmentation results [8]:

[TABLE]

where $\widehat{\boldsymbol{s}}$ is the predicted segmentation map, and $\boldsymbol{s}_{t}$ indicates the ground truth (FreeSurfer) label at each location. A Dice score of 1 indicates perfect segmentation.

We experimented segmenting in the unsupervised setting with standard UNet architectures, using the image MSE and mutual information loss functions. Because of the many structures that share similar intensities, these architectures are not able to produce sensible segmentations that resemble the correct segmentations, and we omit them from these results. Classical unsupervised methods that include sophisticated prior anatomical information take a significant amount time to run, and for T1w we regard FreeSurfer results as an optimistic bound for the T1w data. However, as these methods tend to be focused on specific modalities, there is no annotation tool for cortical and subcortical regions in T2-FLAIR. We evaluate the T2-FLAIR segmentation visually in Figure 5.

5.3 Results

At test time, a new subject only needs to be affinely registered to a template, after which the proposed CNN model evaluates a segmentation estimate. The entire process takes less than a few seconds on an NVidia Titan X GPU.

Fig. 4 shows a series of example segmentations for the T1w dataset demonstrating that our method is able to estimate approximate anatomical structures, reproducing the general location as well as the shape of structures. Fine details, such as details of cortex folding, is not easily captured by the prior encoding, leading to smooth segmentation predictions. Fig. 6 illustrates the average Dice measure across several anatomical regions for T1 scans. We focus on the most prevalent (larger) structures, which can also be evaluated in detail in the visualizations of Figure 4.

Fig. 5 shows results of our algorithm on T2-FLAIR scans. Even with the significantly lower image quality and different tissue contrasts, our algorithm is able to use the prior information to predict plausible segmentations, even given challenging images in this unsupervised scenario.

5.4 Discussion

Our method is able to detect anatomical structures that are guided by image contrast while respecting anatomical shapes according to the prior. Rapid, zero-shot segmentation is a challenging task, and has not been tackled by previous methods. As such, the absence of prior results makes it difficult to fully interpret current results. The detailed FreeSurfer results are an upper bound, which any model is unlikely to achieve in the unsupervised setting. We omit showing results from lower bound (simplistic) baselines, such as the unsupervised U-Net model described above, since these models yielded nonsensical segmentations. To the best of our knowledge, our results are the first for zero-shot neural-network based segmentation of brain structures.

6 Conclusion

In this paper, we introduced a generative probabilistic model that employs a prior model learned through a convolutional neural network to compute segmentations in an unsupervised setting. We can interpret the anatomical prior as encouraging the neural network to predicting segmentation maps that come from a known distribution characterized by $\boldsymbol{z}$ while simultaneously producing images that agree with the observed scan. We demonstrate that our model enables segmentation using convolutional networks leading to rapid inference in a setting where segmentation is traditionally not possible, or takes hours to obtain for a single scan. The integration of priors promises to facilitate accurate anatomical segmentation in a variety of novel clinical problems with limited dataset availability.