Probabilistic Motion Modeling from Medical Image Sequences: Application to Cardiac Cine-MRI

TL;DR

This paper introduces a probabilistic motion model for medical image sequences, enabling realistic motion prediction, interpolation, and transfer, with improved accuracy in cardiac cine-MRI compared to existing methods.

Contribution

It presents a novel probabilistic latent space and temporal dropout training scheme for motion modeling, allowing motion simulation and transfer without inter-subject registration.

Findings

Improved registration accuracy over state-of-the-art algorithms.

Effective motion interpolation from incomplete sequences.

Successful motion transfer demonstrating pathology simulation.

Abstract

We propose to learn a probabilistic motion model from a sequence of images. Besides spatio-temporal registration, our method offers to predict motion from a limited number of frames, useful for temporal super-resolution. The model is based on a probabilistic latent space and a novel temporal dropout training scheme. This enables simulation and interpolation of realistic motion patterns given only one or any subset of frames of a sequence. The encoded motion also allows to be transported from one subject to another without the need of inter-subject registration. An unsupervised generative deformation model is applied within a temporal convolutional network which leads to a diffeomorphic motion model, encoded as a low-dimensional motion matrix. Applied to cardiac cine-MRI sequences, we show improved registration accuracy and spatio-temporally smoother deformations compared to three state-of-the-art registration algorithms. Besides, we demonstrate the model's applicability to motion transport by simulating a pathology in a healthy case. Furthermore, we show an improved motion reconstruction from incomplete sequences compared to linear and cubic interpolation.

Figures

Equations

$$\sum_{t=1}^{T}\mathbb{E}_{z_{t}\sim p_{\gamma}(\cdot|\tilde{z}_{1:T},\bar{t}_{1:T})}\Big{[}\text{log }p_{\theta}(I_{t}|z_{t},I_{0})\Big{]}-\text{KL}\left[q_{\omega}(\tilde{z}_{t}|I_{0},I_{t})\|p(\tilde{z})\right].$$

$${\tilde{z}}_t=r_t\ast {\tilde{z}}_t^{\text{prior}}+(1-r_t)\ast {\tilde{z}}_t^{\text{post}}.$$

Full text

1 Introduction

In medical imaging, an important task is to analyze temporal image sequences to understand physiological processes of the human body. Dynamic organs, such as the heart or lungs, are of particular interest to study as detected motion patterns are helpful for the diagnosis and treatment of diseases. Moreover, recovering the motion pattern allows to track anatomical structures, to compensate for motion, to do temporal super-resolution and motion simulation.

Motion is typically studied by computing pairwise deformations – the registration of each of the images in a sequence with a target image. The resulting dense deformation fields track moving structures from the beginning to the end of the sequence. Providing an invertible and smooth transformation, diffeomorphic registration algorithms such as the SyN algorithm [2], the LCC-demons [11] or recent learning-based algorithms [6, 10] are especially suited for the registration of sequential images. One difficulty is to acquire temporally smooth deformations that are fundamental for tracking. That is why registration algorithms with a temporal regularizer have been proposed [7, 12, 13, 15]. In the computer vision community, temporal video super-resolution and motion compensation are of related interest [5].

However, while these methods produce accurate dense deformations, they do not aim to extract intrinsic motion parameters crucial for building a comprehensive motion model useful for analysis tasks such as motion classification or simulation. Rohé et al. [14] proposed a parameterization, the Barycentric Subspaces, as a regularizer for cardiac motion tracking. Yang et al. [16] generated a motion prior using manifold learning from low-dimensional shapes.

We propose to learn a probabilistic motion model from image sequences directly. Instead of defining a parameterization explicitly or learning from pre-processed shapes, our model captures relevant motion features in a low-dimensional motion matrix in a generic but data-driven way. This learned latent space can be used to fill gaps of missing frames (motion reconstruction), to predict the next frames in the sequence or to generate an infinite number of new motion patterns given only one image (motion simulation). Motion can be also transported by applying the motion matrix on an image of another subject.

The probabilistic motion encoding is learned by generalizing a pair-wise registration method [10] based on Bayesian inference [9] using a temporal regularizer with explicit time dependence. Furthermore, to enforce temporal consistency, we introduce a novel self-supervised training scheme called temporal dropout sampling. The framework is learned in an unsupervised fashion from image sequences of varying lengths. Smooth, diffeomorphic and symmetric deformations are ensured by applying an exponentiation layer, spatio-temporal regularization and a symmetric local cross-correlation metric. Besides motion simulation, the model demonstrates state-of-the-art registration results for diffeomorphic tracking of cardiac cine-MRI. The main contributions are as follows:

• $\bullet$

An unsupervised probabilistic motion model learned from image sequences

• $\bullet$

A generative model using explicit time-dependent temporal convolutional networks trained with self-supervised temporal dropout sampling

• $\bullet$

Demonstration of cardiac motion tracking, simulation, transport and temporal super-resolution

2 Methods

The motion observed in an image sequence with $T+1$ frames is typically described by deformation fields $\phi_{t}$ between a moving image $I_{0}$ and the fixed images $I_{t}$ with $t\in[1,T]$. Inspired by the probabilistic deformation model of [10] based on conditional variational autoencoder (CVAE) [9], we define a motion model for temporal sequences. The model is conditioned on the moving image and parameterizes the set of diffeomorphisms $\phi_{t}$ in a low-dimensional probabilistic space, the motion matrix $z\in\mathbb{R}^{d\times T}$, where $d$ is the size of the deformation encoding per image pair. Each column’s $z_{t}$-code corresponds to the deformation $\phi_{t}$. To take temporal dependencies into account, $z_{t}$ is conditioned on all past and future time steps. To learn this temporal regularization directly from data, we apply Temporal Convolutional Networks [3] with explicit time dependence and temporal dropout sampling enforcing the network to fill time steps by looking at given past and future deformations. An illustration of the model is shown in Fig. 1a.

2.0.1 Probabilistic Motion Model

Our motion model consists of three distributions. First, the encoder $q_{\omega}(\tilde{z}_{t}|I_{0},I_{t})$ maps each of the image pairs $(I_{0},I_{t})$ independently to a latent space denoted by $\tilde{z}_{t}\in\mathbb{R}^{d}$. Second, as the key component of temporal modeling, these latent vectors $\tilde{z}_{t}$ are jointly mapped to the motion matrix $z$ by conditioning them in all past and future time steps and on the normalized time $\bar{t}$: $p_{\gamma}(z|\tilde{z}_{1:T},\bar{t}_{1:T})$. Finally, the decoder $p_{\theta}(I_{t}|z_{t},I_{0})$ aims to reconstruct the fixed image $I_{t}$ by warping the moving image $I_{0}$ with the deformation $\phi_{t}$. This deformation $\phi_{t}$ is extracted from the temporally regularized $z_{t}$-codes. The decoder is conditioned on the moving image by concatenating the features at each scale with down-sampled versions of $I_{0}$.

The distributions $q_{\omega}$, $p_{\gamma}$, $p_{\theta}$ are approximated by three neural networks with trainable parameters $\omega$, $\gamma$, $\theta$. During training, a lower bound on the data likelihood is maximized with respect to a prior distribution $p(\tilde{z}_{t})$ of the latent space $\tilde{z}_{t}$ (cf. CVAE [9]). The prior $p(\tilde{z}_{t})$ is assumed to follow a multivariate unit Gaussian distribution with spherical covariance $I$: $p(\tilde{z}_{t})\sim\mathcal{N}(0,I)$. The objective function results in optimizing the expected log-likelihood $p_{\theta}$ and the Kullback-Leibler (KL) divergence enforcing the posterior distribution $q_{\omega}$ to be close to the prior $p(\tilde{z}_{t})$ for all time steps:

[TABLE]

Unlike the traditional CVAE model, the temporal regularized $z_{t}$-code is used in the log-likelihood term $p_{\theta}$ instead of the $\tilde{z}_{t}$. We model $p_{\theta}$ as a symmetric local cross-correlation Boltzmann distribution with the weighting factor $\lambda$. Encoder and decoder weights are independent of the time $t$. Their network architecture consists of convolutional and deconvolutional layers with fully-connected layers for mean and variance predictions in the encoder part [9]. We use an exponentiation layer for the stationary velocity field parameterization of diffeomorphisms [10], a linear warping layer and diffusion-like regularization with smoothing parameters $\sigma_{G}$ in spatial and $\sigma_{T}$ in temporal dimension.

2.0.2 Temporal Convolutional Networks with Explicit Time Dependence

Since the parameters of encoder $q_{\omega}$ and decoder $p_{\theta}$ are independent of time, the temporal conditioning $p_{\gamma}$ plays an important role in merging information across different time steps. In our work, this regularization is learned by Temporal Convolutional Networks (TCN). Consisting of multiple 1-D convolutional layers with increasing dilation, TCN can handle input sequences of different lengths. TCN have several advantages compared to recurrent neural networks such as a flexible receptive field and more stable gradient computations [3].

The input of the TCN is the sequence of $\tilde{z}$ concatenated with the normalized time $\bar{t}=t/T$. Providing the normalized time explicitly, provides the network with information on where each $\tilde{z}$ is located in the sequence. This supports the learning of a motion model from data representing the same type of motion with varying sequence lengths. The output of the TCN is the regularized motion matrix $z$. We use non-causal instead of causal convolutional layers to also take future time steps into account. We follow the standard implementation using zero-padding and skip connections. Each layer contains $d$ filters. A schematic representation of our TCN is shown in Fig. 1b. For cyclic sequences, one could use a cyclic padding instead of zero-padding, for example by linking $\tilde{z}_{T}$ to $\tilde{z}_{0}$. However, in case of cardiac cine-MRI, one can not assume the end of a sequence coincides with the beginning as 5-10% of the cardiac cycle are often omitted [4].

2.0.3 Training with Temporal Dropout Sampling

Using Eq. 1 for training could lead to learning the identity transform $z\approx\tilde{z}$ in the TCN $p_{\gamma}$ such that deformations of the current time step are independent of past and future time steps. To avoid this and enforce the model to search for temporal dependencies during the training, we introduce the concept of temporal dropout sampling (TDS). In TDS, some of the $\tilde{z}_{t}$ are sampled from the prior distribution $p(\tilde{z})$ instead of only sampling from the posterior distribution $q_{\omega}(\tilde{z}_{t}|I_{0},I_{t})$ as typical for CVAE. At the time steps the prior has been used for sampling, the model has no knowledge of the target image $I_{t}$ and is forced to use the temporal connections within the TCN in order to minimize the objective.

More precisely, at each time step $t$, a sample from the prior distribution $\tilde{z}_{t}^{\text{prior}}\sim p(\tilde{z}_{t})$ is selected instead of a posterior sample $\tilde{z}_{t}^{\text{post}}\sim q_{\omega}(\tilde{z}_{t}|I_{0},I_{t})$ using a binary Bernoulli random variable $r_{t}$. All independent Bernoulli random variables $r\in\mathbb{R}^{T}$ have the success probability $\delta$. The latent vector $\tilde{z}_{t}$ can be defined as:

[TABLE]

Fig. 1c illustrates the TDS procedure. At test time, for each time step independently, one can either draw $\tilde{z}_{t}$ from the prior or take the encoder’s prediction.

3 Experiments

We evaluate our motion model on 2-D cardiac MRI-cine data. First, we demonstrate accurate temporal registration by evaluating motion tracking and compensation of the cardiac sequence, taking the end-diastolic (ED) frame as the moving image $I_{0}$. Stabilization, is accomplished by warping all frames $I_{t}$ to the ED frame. Pair-wise registration results are presented for ED-ES (end-systolic) frame pairs. Second, we present motion transport, motion sampling and reconstruction with a limited number of frames.

3.0.1 Data

We used 334 short-axis sequences acquired from different hospitals including 150 sequences from the Automatic Cardiac Diagnosis Challenge 2017 (ACDC [4]). The remaining cases were obtained from the EU FP7-funded project MD-Paedigree (Grant Agreement 600932), mixing congenital heart diseases with healthy and pathological images from adults. The cine images were acquired in breath hold using 1R-R or 2R-R intervals with a retrospective or prospective gating. The sequence length $T$ varied from 13 to 35 frames. We used the 100 cases from ACDC that contain ED-ES segmentation information for testing while the remaining sequences were used for training. All slices were resampled with a spacing of 1.5 $\times$ 1.5 mm and cropped to a size of 128 $\times$ 128 pixels.

3.0.2 Implementation Details

The encoder $q_{\omega}$ consisted of 4 convolutional layers with strides (2, 2, 2, 1) and dense layers of size $d$ for mean and variance estimation of the VAE. The TCN consisted of four 1-D convolutional layers with dilations (1, 2, 4, 8), same padding, a kernel size of 3 and skip connections (cf. Fig. 1b). The decoder $p_{\theta}$ had 3 deconvolutional and 1 convolutional layer before the exponentiation and warping layers (Fig. 1a). The regularization parameters $\sigma_{G}$ and $\sigma_{T}$ were set to 3 mm respectively 1.5. The loss weighting factor $\lambda$ was chosen empirically as $6\cdot 10^{4}$. The deformation encoding size $d$ was set to 32. The dropout sampling probability $\delta$ was 0.5. We applied a first-order gradient-based method for stochastic optimization (Adam [8]) with a learning rate of 0.00015 and a batch size of one. We performed data augmentation on-the-fly by randomly shifting, rotating, scaling and mirroring images. We implemented the model in Tensorflow [1] with Keras111https://github.com/fchollet/keras. The training time was 15h on a NVIDIA GTX TITAN X GPU.

3.1 Registration: Tracking and Motion Compensation

We compare our model with and without the temporal dropout sampling (Our w/o TDS) with three state-of-the-art methods: SyN [2], the learning-based probabilistic registration (LPR, 2-D single-scale version [10]) and the b-spline-based 4D algorithm in elastix [12]. In contrast to the results in [10], in this work, LPR was trained in 2-D taking all images from a sequence into account, not only ED/ES pairs. The following results are reported for full sequences, except the metrics based on segmentations which are only reported on frame pairs with provided ground truth information (ED/ES pairs).

In Fig. 2, tracking results are visualized for the test data taking sequences of all slices where segmentation in ED and ES frames are available, resulting in 677 sequences. We report the root mean square error (RMSE), the spatial (Spatial Grad.) and temporal gradient (Temp. Grad.) of the displacement fields for evaluating the smoothness of the resulting deformations. Our model shows spatially and temporally smoother deformations. We also report DICE scores and 95%-tile Hausdorff distances (HD in mm) on five anatomical structures: myocardium (LV-Myo) and epicardium (LV) of the left ventricle, left bloodpool (LV-BP), right ventricle (RV) and LV+RV. Note, that DICE scores and HD were evaluated on ED-ES frame pairs only. The proposed method showed improved mean DICE scores and smaller mean HD of 84.6%, 6.2mm (w/o TDS: 84.7%, 6.1mm) compared to SyN, LPR, 4D-Elastix with (82.7%, 7.0mm), (82.1%, 6.6mm) respectively (83.7%, 6.3mm%). Compared to training without temporal dropout sampling, HD and DICE scores show minimal differences, indicating that using TDS does not degrade registration accuracy while improving deformation regularity. Furthermore, only TDS offers consistent motion simulation and temporal interpolation.

In the bottom right of Fig. 2, LV volume curves computed by warping the ED mask are plotted. One can see that the SyN algorithm underestimates big deformations. Both, the volume and the gradient metrics show smoother deformations for both versions of our motion model compared to the SyN, LPR and 4D-Elastix algorithms. Visual results for one case including tracking, determinant of Jacobians and compensated motion are shown in Fig. 3. Motion compensation was done by warping the $I_{t}$’s frame with inverted diffeomorphisms.

3.2 Sampling, Sequence Reconstruction and Motion Transport

We then evaluate motion sampling, reconstruction and transport. We extract simulated reconstructed motion patterns if we provide our model with different subsets of images from the original sequence. In the time steps without frames, the motion matrix $z$ is created by randomly sampling $\tilde{z}_{t}$ from the prior distribution as depicted in Fig. 1c (take same $\tilde{z}_{t}$ for all slices of one volume). Here, we choose sampling over interpolation in $\tilde{z}$-space to remain an uncertain, a probabilistic, estimation of the interpolated deformations. The left side of Fig. 4 shows LV volume curves and reconstruction errors if every second, every fifth, only the 10th, no frame (sampling) or the first 5 frames are provided besides the moving frame (ED). The LV volume errors are computed on all test sequences by taking the mean absolute differences between sampled and tracking volumes. We compare with linear and cubic interpolation of velocities (extracted from tracking) between given times. One can see that our model performs better in recovering the LV volume, compared to linear and cubic interpolation when fewer frames are provided. Given the first five frames or only one additional frame (10th frame), the model estimates the motion consistently with plausible cardiac motion patterns for the missing time steps. In the cases of providing every second and every fifth frame, our method performs equally good or marginally better. Note, the motion simulation given only the ED frame (sampling) does not overlap with the original motion, which is not intended. Nevertheless, one can see cardiac specific motion patterns such as the plateau phase before the atrial systole.

Motion can be transported by taking the motion matrix $z$ from one sequence and apply it on the ED frame $I_{0}$ of another sequence. The right side of Fig. 4 shows the transport of a pathological motion to a healthy subject and vice versa. The resulting simulated motion shows similar heart contractions and motion characteristics as the originating motion while the transported deformations are adapted to the target image without requiring explicit inter-subject registration.

4 Conclusions

In this paper, we presented an unsupervised approach for learning a motion model from image sequences. Underlying motion factors are encoded in a low-dimensional probabilistic space, the motion matrix, in which each column represents the deformation between two frames of the sequence. Our model demonstrated accurate motion tracking and motion reconstruction from missing frames, which can be useful for shorter acquisition times and temporal super-resolution. We also showed motion transport and simulation by using only one frame. Limitations of the presented approach include the support for 3-D image sequences and the generalization to other use cases such as respiratory motion.

For future work, we aim to explore these points and especially the spatial coherence between slices for 3-D applications and the influence of using different training datasets (pathological and non-pathological) on the learned motion matrix.

**Disclaimer: ** The concepts and information presented in this paper are based on research results that are not commercially available.