Mixture Modeling of Global Shape Priors and Autoencoding Local Intensity Priors for Left Atrium Segmentation

TL;DR

This paper introduces a novel segmentation method combining mixture modeling of global shape priors with autoencoding local intensity priors, improving accuracy in complex MRI images of the left atrium.

Contribution

It proposes a mixture of Gaussians for shape priors learned in high-dimensional space, integrated with deep autoencoders for intensity modeling, within a MAP framework for improved segmentation.

Findings

Effective segmentation of left atrium in MRI images.

Handles nonlinear and multimodal shape variations.

Outperforms traditional methods on benchmark datasets.

Abstract

Difficult image segmentation problems, for instance left atrium MRI, can be addressed by incorporating shape priors to find solutions that are consistent with known objects. Nonetheless, a single multivariate Gaussian is not an adequate model in cases with significant nonlinear shape variation or where the prior distribution is multimodal. Nonparametric density estimation is more general, but has a ravenous appetite for training samples and poses serious challenges in optimization, especially in high dimensional spaces. Here, we propose a maximum-a-posteriori formulation that relies on a generative image model by incorporating both local intensity and global shape priors. We use deep autoencoders to capture the complex intensity distribution while avoiding the careful selection of hand-crafted features. We formulate the shape prior as a mixture of Gaussians and learn the corresponding parameters in a high-dimensional shape space rather than pre-projecting onto a low-dimensional subspace. In segmentation, we treat the identity of the mixture component as a latent variable and marginalize it within a generalized expectation-maximization framework. We present a conditional maximization-based scheme that alternates between a closed-form solution for component-specific shape parameters that provides a global update-based optimization strategy, and an intensity-based energy minimization that translates the global notion of a nonlinear shape prior into a set of local penalties. We demonstrate our approach on the left atrial segmentation from gadolinium-enhanced MRI, which is useful in quantifying the atrial geometry in patients with atrial fibrillation.

Tables

Table 1. Table 1: Dice scores for test results.

Autoencoder	yes	yes	no	no
Mixture Model	yes	no	yes	no
Median Dice Coefficient	0.743	0.747	0.729	0.739

Equations

$$\log p(\mathcal{S}|\mathcal{I})=\log p(\mathcal{S},\mathcal{Z}=z|\mathcal{I})-\log p(\mathcal{Z}=z|\mathcal{S},\mathcal{I}),$$

$$p(\mathcal{Z}|{\mathcal{S}}^{(t)},\mathcal{I})$$

$$\log p(\mathcal{S}|\mathcal{I})$$

$${\mathcal{S}}^{(t+1)}$$

$${\mathbb{S}}_k=\left\{{\mathcal{S}}_k^{\beta }={\mu }_k+{\Psi }_k{\vec{\beta }}_k|{\mathcal{S}}_k^{\beta }\in :{\mathbb{R}}^{M\times 3},{\vec{\beta }}_k\in {\mathbb{R}}^{d_k}\right\}$$

Full text

1 Introduction

Automatic image segmentation is an important enabling technology in most medical imaging applications that involve soft tissue imaging, e.g., neurology, cardiology, and oncology. In particular, the preoperative anatomical representation of the left atrium (LA) is important for ablation guidance, fibrosis quantification, and biophyiscal modeling in artial fibrillation patients [1, 2, 3, 4]. Difficult segmentation problems typically have two aspects. First, image features, such as intensity and texture are noisy and unreliable. Second, anatomical boundaries are ill-defined, often in the presence of low contrast, background clutter, and partial volumes, while irregular shapes with high variability limit the ability to find invariant features. Nonetheless, if one operates in a context with expectations of particular classes of anatomies (e.g., LA), such challenges can be addressed by means of shape prior information to guide and constrain the segmentation process; motivating a Bayesian formulation. In LA segmentation context, shape-driven methods were found to be the most appropriate in addressing inherent challenges [5]; thin myocardial wall, surrounding anatomical structures with similar image intensities, and topological variants pertaining to pulmonary viens arrangments [6]. In this paper, we propose a Bayesian surface-based segmentation framework that is based on a mixture-based global shape prior along with a feature-based local intensity prior.

Traditionally, when dealing with low-quality image segmentation, statistical shape information has proved to be helpful in delineating correct object boundaries (e.g., [7, 8]). For example, active shape models [9] and their variants incorporate over-restrictive shape constraints in the form of statistical shape models (SSM) by limiting the solution to some low-dimensional linear subspace defined via training shape exemplars. However, these approaches are limited in their ability to accommodate shapes that are not represented in a low-dimensional description; a typical situation arises in applications with small and large-scale shape variability (e.g., [6]). Further, these methods only handle unimodal Gaussian-like shape densities. When it comes to modeling complex shape distributions, a single Gaussian can not adequately model cases of nonlinear shape variation where the probability distribution is multimodal [10].

On the other end of the spectrum, multiatlas-based segmentation (MAS) approaches require a large database of atlases to capture wide range of shape variation where label maps are propagated to the testing image through registration. In essence, MAS can be viewed as a nonparametric estimate of the prior probabilities that converges to the true densities with the number of atlases [11]. To reduce the computational burden introduced by registration, atlas selection is usually performed to exclude irrelevant atlases that might misguide the segmentation process [12]. Generally, kernel-based methods (e.g., [8, 13]) are popular strategies for dealing with complex distributions in explicitly or implicitly represented training data. Under mild assumptions, they converge to the true distribution in the limit of infinite sample size [14]. For example, Cremers et al. [8] used kernel density estimate (KDE) to derive non-linear shape priors which are based on a shape distance between implicitly embedded training shapes. Nonetheless, such methods are prone to over-fitting due to small sample size in high-dimensional shape spaces, limiting the generality of the resulting models to fit unseen examples. Specifically, in multivariate density estimation, KDE requires larger kernel width to accommodate more exemplars as the dimension of a variable increases. This will eventually result in model under-fit due to high bias [15]. Further, for kernel based methods, the dimension of the parameter space is proportional to the training sample size.

A finite mixture of Gaussians can approximate a full kernel density, to capture nonlinearities in the distribution or subpopulations in the underlying shape subspace [16] where the expectation-maximization (EM) framework is usually used to find the maximum likelihood estimate of mixture parameters. However, modeling shape distribution after being projected onto low-dimensional subspace (e.g., [14, 16]) will often collapse or mix the subpopulations, which would derail learning of the mixture structure of the underlying shape space [17]. Nonlinearity of shape statistics can also be modeled by lifting training shapes to a higher, probably infinite, dimensional feature (a.k.a. kernel) space where the shape distribution is assumed to be Gaussian distributed [10], yet this approach results in an infinite-dimensional optimization scheme while sacrificing the efficiency of optimizing in low-dimensional subspaces [14]. Further, one can settle for only an approximate solution for the reverse mapping from feature space to shape space.

The proposed bayesian formulation relies on a maximum-a-posteriori estimation from a generative statistical model that incorporates global, nonlinear shape priors modeled as a mixture of Gaussian components, and local, nonlinear intensity prior as automatically learned image features via deep autoencoders. The Gaussian mixture model (GMM) takes into account the linear subspace spanned by each mixture component to avoid the classical problem of model over fitting in high-dimensional spaces [18], and can adequetely model non-linear shape variations. To use these shape priors in segmentation, we treat the identity of the mixture component as a latent variable while marginalizing it within a generalized expectation-maximization framework (GEM) [19]. We present a conditional maximization-based scheme that alternates between a closed-form solution for component-specific shape parameters that provides a global update-based optimization strategy, and an intensity-based energy minimization that translates the global notion of a nonlinear shape prior into a set of local penalties. Preliminary results show improved accuracy from traditional shape-based approaches.

2 Methods

Given an image $\mathcal{I}\in\mathbb{R}^{D}$, where $D$ is the total number of voxels and $\mathcal{I}(\mathbf{x})$ is the intensity value of the voxel located at $\mathbf{x}\in\mathbb{R}^{3}$, the Bayesian formulation of the segmentation problem amounts to finding the optimal surface $\mathcal{S}^{*}$ that maximizes the log-posterior probability $p(\mathcal{S}|\mathcal{I})$, where a shape’s surface is represented by a dense set of geometrically consistent $M-$ points (landmarks) $\{\mathbf{x}^{\mathcal{S}}_{i}\}_{i=1}^{M}$. In order to obtain segmentations that preserve the global shape characteristics of the shape population of interest, the segmentation process is influenced by the prior, in the form of the shape probability distribution. We model such a prior distribution as a finite mixture of Gaussians, parameterized by the mixture component identity $\mathcal{Z}$, which is treated as a latent variable, and the model parameters $\{\Theta_{z}\}_{z=1}^{K}$ for $K-$components. Hence, given a latent variable $\mathcal{Z}=z$, the log-posterior can be written as111For notational simplicity, we will refer to $p(\mathcal{Z}=z)=p(\mathcal{Z})$ hereafter.,

[TABLE]

where $p(\mathcal{Z}|\mathcal{S},\mathcal{I})$ defines a conditional probability of the latent variable $\mathcal{Z}$ given the input image $\mathcal{I}$ and the estimated surface $\mathcal{S}$.

2.1 Generalized EM for shape-driven segmentation

The use of log-posterior allows us to marginalize the over the latent variable $\mathcal{Z}$. Using generalized expectation-maximization (GEM) [19, p. 318], we iterate to find the mode of the marginal posterior $p(\mathcal{S}|\mathcal{I})$ by averaging over the latent variable $\mathcal{Z}$. GEM starts with an initial estimate $\mathcal{S}^{o}$ of the surface and iteratively refines it by marginalizing over $\mathcal{Z}$. Within each iteration, the posterior $p(\mathcal{S}|\mathcal{I})$ is guaranteed to increase, i.e., converge to a local maxima. Given the current surface guess $\mathcal{S}^{(t)}$, GEM consists of two steps: (1) E-step, which computes the latent conditional probability distribution $p(\mathcal{Z}|\mathcal{S}^{(t)},\mathcal{I})~{}\forall\mathcal{Z}\in\{1,..,K\}$ and (2) M-step, which finds a new surface $\mathcal{S}^{(t+1)}$ such that $p(\mathcal{S}^{(t+1)}|\mathcal{I})\geq p(\mathcal{S}^{(t)}|\mathcal{I})$. The details of the algorithm are as follows.

E-step: The conditional probability of the latent variable $\mathcal{Z}$ given the input image and its labeling function can be computed as,

[TABLE]

where the intensity model $p(\mathcal{I}|\mathcal{S}^{(t)},\mathcal{Z})$ in the vicinity of the current surface guess $\mathcal{S}^{(t)}$ is assumed to be conditionally independent of the shape generating component, leaving the E-step to be fully shape-driven. Note the effect of the prior probability $p(\mathcal{Z})$ on the expectation step, where the conditional component distribution would favor mixture components with higher support/proportion, the M-step afterwards would pull the new surface towards the most probable mixture component(s) which generated the shape to be segmented as being compatible with the local intensity and global shape priors.

M-step: Taking the expectations on both sides of (1) while treating $\mathcal{Z}$ as a random variable with the distribution $p(\mathcal{Z}|\mathcal{S}^{(t)},\mathcal{I})$ yields,

[TABLE]

where the left side of (1) does not depend on $\mathcal{Z}$. The second term in (3) is maximized when $\mathcal{S}=\mathcal{S}^{(t)}$ (the key result of EM [19]). Hence, it is sufficient for the new surface estimate $\mathcal{S}^{(t+1)}$ to maximize the first term of (3), which is the expected complete-data log-likelihood.

[TABLE]

The intensity term $p(\mathcal{I}|\mathcal{S},\mathcal{Z})$ gives rise to an intensity model in the vicinity of an estimated surface that resembles the $\mathcal{Z}-$th component. Intensity features (learned via a deep autoencoder [20]) $\mathcal{F}_{i}^{\mathcal{I}}$ are assumed to follow a normal distribution in the learned feature space with parameters $\Theta_{i}^{\mathcal{F}}=\{\mu_{i}^{\mathcal{F}},\Sigma_{i}^{\mathcal{F}}\}$ that are estimated from the training data, where $\mu_{i}^{\mathcal{F}}$ is the mean intensity feature at the $i-$th point and $\Sigma_{i}^{\mathcal{F}}$ is the corresponding covariance matrix. Notice that $\mathcal{S}^{(t+1)}$ increases the first term of (3) since it is obtained based on its maximization. In addition, $\mathcal{S}^{(t+1)}$ decreases the second term of (3), because such a term is maximized only when $\mathcal{S}=\mathcal{S}^{(t)}$. Hence, if $\mathcal{S}^{(t+1)}=\mathcal{S}^{(t)}$, then the GEM method has converged to a local maxima of the posterior $p(\mathcal{S}|\mathcal{I})$.

2.2 Mixture modeling of global shape priors

The combination of high dimensional shape space and a small number of example shapes usually makes modeling of shape priors subject to curse-of-dimensionality, which hinders the effectiveness of mixture learning and density estimation in high dimensional spaces. To avoid over fitting, we can parameterize the covariance matrix of a component, in a manner similar to [18], through its eigen decomposition. Hence we can learn the Gaussian mixture model taking into consideration the dominant linear subspace spanned by each mixture component, where the noise variance is assumed to be isotropic and contained in a subspace orthogonal to the component’s subspace.

We first model a dense set of $M$ homologous landmarks via particle based modeling. Assume that the anatomy-specific shape space of all shapes defined using the $M-$landmarks representation is comprised of $K-$mixture components. Each component is parameterized by $\Theta_{k}=\{\pi_{k},\mu_{k},\Lambda_{k},$ $\Psi_{k},d_{k},\sigma_{k}\}$ and learned via high-dimensional EM [18] where $\pi_{k}=p(\mathcal{Z}=k)$ is the component probability/proportion, $\mu_{k}$ is the mean vector, $d_{k}$ is the intrinsic dimension of the $k-$th mixture component, i.e., the number of dominant modes of shape variation, $\Lambda_{k}$ is a diagonal matrix containing the largest $d_{k}$ eigenvalues of the component’s full covariance matrix $\Sigma_{k}$, $\Psi_{k}$ is the orthonormal matrix containing the corresponding $d_{k}$ eigenvectors, and $\sigma_{k}$ is the standard deviation of the off-subspace noise. Hence the component-specific subspace $\mathbb{S}_{k}$ can be defined as follows:

[TABLE]

where $\vec{\beta}_{k}\in\mathbb{R}^{d_{k}}$ encodes the shape parameters w.r.t. the $k-$th component subspace $\mathbb{S}_{k}$.The shape probability distribution $p(\mathcal{S})$ is defined as a weighted sum of component-conditional probabilities $\{p(\mathcal{S}|\Theta_{k})\}_{k=1}^{K}$ using the mixture weights $\{\pi_{k}\}$. Notice, that the component-wise model beyond the first $d_{k}$ modes consists of an isotropic variance in the directions orthogonal to the subspace, which has the effect of penalizing point wise differences from the learned, $d_{k}$-dimensional subspace. Hence, the surface probability conditional on the $k-$th mixture component becomes a product of two marginal and indepepdent Gaussians; within-subpace and off-subspace [18].

We illustrate this training step in Fig. 1. These data are from the samples selected for training and illustrate a low-dimensional subspace embedding of every sample colored according to their class membership in a mixture model of 3 components along with example LA segmentations from each. One can observe key differences between each component such as the number and shape of pulmonary veins.

2.3 Autoencoding local intensity priors

Similar to the active shape model (ASM) approach of [16], we model local intensity profiles at the object surface consisting of the normalized first derivative of image intensities oriented along the approximate normal to the surface and centered at each landmark location, with radius of $\ell$ voxels with total length of $L=2\ell+1$ (Fig. 2 ). Corresponding profiles (i.e., those that share the same landmarks ) are combined across all training images but are treated independently of one another. The traditional approach to segmentation is to use the Mahalanobis distance of candidate profiles to the training data as a cost function to be minimized, however, it assumes the feature space to be first order linear and normally distributed. In order to overcome this limitation, we employ deep autoencoders to adaptively learn the higher order non-linear features. In this case, we use a 2 layer $\times$ 10 hidden units sparse autoencoder to generate pseudolinear representations of the feature space for each landmark. Due to the computationally expensive nature of this step we determined these parameters through a qualitative empirical analysis but there is potential for improvement via more rigorous cross-validation. To improve the statistics and capture lateral intensity changes, we append additional profiles which are parallel to the main profile and, offset by a subvoxel distance, with interpolated intensity values, to form “thick” profiles. We repeat this for multiple image resolutions and store the autoencoder as well as the encoded output.

2.4 Segmentation via conditional maximization

The marginal posterior optimization problem in (4) involves two types of latent variables, the labeling function to be estimated $\mathcal{S}$ and the component-specific shape parameters $\vec{\beta}_{k}$. To seek an optimal solution for $\mathcal{S}$, we propose a conditional maximization-based scheme [19, p. 312] that alternates between two optimization phases. In the first phase, we optimize for $\{\vec{\beta}_{k}\}_{k=1}^{K}$ given $\mathcal{S}$, and in the second phase we optimize for $\mathcal{S}$ given individual shape parameters $\{\vec{\beta}_{k}\}_{k=1}^{K}$. Starting with the coarsest image resolution, we take the mean landmark positions from the training data as the initial model and iterate through these steps.

**Localized feature-based optimization: ** For each landmark position, we generate a series of overlapping profiles of $L-$dimension that extend above and below the original profile by some specified search length $s$. For each of these candidate profiles, we generate features using pre-trained autoencoders. We then take the Mahalanobis distance of each candidate profile to the encoded training features and update the landmark position to the center of the optimal profile, oen with the smallest distance. This step ensures that landmarks are locally optimal, but they may no longer by globally optimal as there is no shape information encoded into these features.

Component-specific shape parameters optimization: For a given surface $\mathcal{S}$ (from the intensity-based step), component-specific shape parameters $\vec{\beta}_{k}$ can be obtained by projecting the surface onto the component-specific subpace, leading a closed-form solution $\vec{\beta}_{k}=\Psi_{k}^{T}(\mathcal{S}-\mu_{k})$. By limiting the model parameters to what is considered a normal contour with respect to shape, the landmark positions become globally optimal. We repeat this a predetermined number of times at each resolution, coarse to fine.

3 Results

We tested these methods on 100 3D MRI images, separated into 80 for training and 20 for testing. Fig. 3 shows a sample MRI image with the provided segmentation. The segmented volume is also separately displayed in blue along with modeled landmark positions. We first develop a shape model of 2048 homologous landmark positions modeled with a 3 component Gaussian mixture model. We then model local intensity gradients for each position across all samples using an autoencoder on profiles 11 voxels in length (optimized through cross-validation). In order to improve model statistics, we encode additional parallel information to create ”thick” profiles. This model was then used to segment the 20 test images according to the methods described above. To compare the effectiveness of our method we evaluate results both with and without autoencoded local intensity profiles as well as with and without the mixture model. Results are in Table 1 . While there appears to be some improvement in median Dice coefficient, this may not be the best metric for evaluating the efficiency of the proposed method [21]. We also evaluated the euclidean distance between final landmark positions from our segmentations to those of the ground truth. This also allows for a more detailed point by point analysis of error as shown in Fig. 4. Here we see two examples from the test data set showing clear improvement in segmentation results when using autoencoded profiles. We truncate the color scale at 10 (1 cm) and flag any higher error as an extreme outlier. Error is generally limited to image edges around the pulmonary veins when using autoencoded intensity profiles as opposed to nearly universally high error without. These results could be further improved by initializing the model with the component mean shapes as opposed to the global mean which is how we are presently doing it.

4 Conclusion

This paper proposed a shape/feature-based generative model for left atrium segmentation by modeling non-Gaussian global shape priors as mixture of Gaussians on landmark-based representations of training data and learning feature-based representations for local intensity priors. The mixture parameters were learned in the high dimensional shape space by taking into account component-specific subspaces to avoid over fitting. The method used a variant of ASM-based segmentation framework that relies on a maximum a-posteriori estimation with a marginalization over class membership within a generalized EM framework. While these results are preliminary, they suggest that ASM’s for LA segmentation can be improved through optimized prior modeling. Local optimization through autoencoding of non-linear intensity features and global optimization through mixture modeling of shape parameters increases the accuracy of the original method

Acknowledgment: This work was supported by the National Institutes of Health [grant numbers R01-HL135568-01 and P41-GM103545-19].