{'fraction_non_alphanumeric': 0.041401955885502176, 'fraction_numerical': 0.015918338211149977, 'mean_word_length': 4.378119001919385, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Latent linear dynamics in spatiotemporal medical data

Spatiotemporal imaging is common in medical imaging, with applications in e.g. cardiac diagnostics, surgical guidance and radiotherapy monitoring. In this paper, we present an unsupervised model that identifies the underlying dynamics of the system, only based on the sequential images. The model maps the input to a low-dimensional latent space wherein a linear relationship holds between a hidden state process and the observed latent process. Knowledge of the system dynamics enables denoising, imputation of missing values and extrapolation of future image frames. We use a Variational Auto-Encoder (VAE) for the dimensionality reduction and a Linear Gaussian State Space Model (LGSSM) for the latent dynamics. The model, known as a Kalman Variational Auto-Encoder, is end-to-end trainable and the weights, both in the VAE and LGSSM, are simultaneously updated by maximizing the evidence lower bound of the marginal log likelihood. Our experiment, on cardiac ultrasound time series, shows that the dynamical model provide better reconstructions than a similar model without dynamics. And also possibility to impute and extrapolate for missing samples.


Introduction
Today, most medical imaging modalities support some form of time-resolved imaging. In some modalities, like ultrasound, it is the default mode of operation, while in others it is known under different names depending on the application, e.g. fluoroscopy [28], 4DCT [22] and 4D flow MRI [18]. It is used for analysis, in e.g. cardiac diagnostics [1] and for guidance, monitoring and control in e.g. interventional and intraoperative percutaneous procedures and surgeries [2,31] and image-guided radiation therapy [20,26].
In this paper we want to uncover the dynamics in a medical image time series based on nothing but the images themselves. Our aim is to learn a parameterized mathematical function, f , such that the image y t at time t is described by the dynamical system y t = f (y 1:t−1 ) + ǫ t , where y 1:t−1 = {y 1 , . . . y t−1 } are the previously observed images and ǫ t is noise. This representation makes it possible to predictŷ t = f (y 1:t−1 ) based on observations up to time T , i.e. to impute (t < T ) or extrapolate (t > T ) images. In the medical domain, there are many examples of motion modeling based on side information [19,27] such as from X-ray tracking of implanted markers or surface tracking. But, to the best of our knowledge, there are no examples where the motion model is based directly and exclusively on the temporal medical images under study. Temporal medical images contains high dimensional data where the dynamics, due to cyclic and deformable motion, is nonlinear. While linear assumption is preferable for several reasons, like tractable filtering and smoothing posteriors, it is inappropriate for the raw image series. The idea of this paper is to reduce the high dimensional data to a low-dimensional latent space wherein the linearity assumptions are valid. For this we combine techniques from generative and dynamic modelling and train the model end-to-end. In literature, this is referred to as Deep State-Space Models (DSSM) [5,12,15].

Method
The most common and well-established approach for dynamical systems is to use a Linear Time Invariant (LTI) State Space model [4]. An LTI state space model is characterized by a linear relationship between a hidden state process z t and an observed state process x t . The hidden state process is assumed to be a first-order Markov process, i.e. p(z t+1 | z t , . . . , z 1 ) = p(z t+1 | z t ), and the observed state process is independent given the state process. In discrete time where A is called the transition matrix and C the emission matrix, w t and v t are process and measurement noise in the model respectively.
One of the main advantages of an LTI state space model is the highly efficient methods for filtering and smoothing it enables. However, there are two major issues that prevent a direct application of an LTI state space model on temporal medical images: the high dimensionality of the observations and the linear dynamics. Even for small 2D images, like 64 × 64 pixels, the observed vector space lies in R 4096 . Since smoothing algorithms often involve matrix inverses of the emission matrix, this quickly becomes prohibitive. Moreover, temporal medical images often undergo non-linear deformations, which violate the assumption of linear dynamics. In this paper, we explain the high dimensional image space, y ∈ Y, in a lower dimensional latent space wherein the latent observations, x ∈ X , and dynamics, z ∈ Z, can be represented linearly.
We use a Kalman Variational Auto-Encoder [5] where the mappings to and from the latent representation are given by a Variational Auto-Encoder (VAE) [13] and the latent dynamics are described by a Linear Gaussian State Space Model (LGSSM) [11].

Variational Auto-Encoder
Similar to traditional auto-encoders [14], VAEs embeds the input y to a lower latent representation x using an encoder and reconstruct the original input with a decoder. They differ in that VAEs are generative and reconstruct the data distribution p θ (y) instead of the input sample y. In this case the true posterior p θ (x | y) is intractable. By approximating the variational posterior as a multivariate Gaussian, q φ (x | y) = N (x | µ enc , Σ enc ), where µ enc and Σ enc are learned by the encoder, it is possible to sample from the variational approximation. From the KL divergence between the approximate and the real posterior distribution one can obtain a lower bound on the true likelihood (see (7) and (8) in Supplementary A.1), where the prior distribution over the latent space is usually chosen to be a multivariate Gaussian p θ (x) = N (x | 0, I). This lower bound is called the evidence lower bound (ELBO) and can be estimated by sampling [13].

Linear Gaussian State Space Model
In a VAE, each sample, x t ∈ R p , from the approximate posterior distribution is normally distributed with mean, µ enc t , and covariance Σ enc t . If we assume Gaussian noise, it follows that the state space vector, z t , in our linear state space model also follows a normal distribution. More precisely, we have a Linear Gaussian State Space Model (LGSSM), where R and Q are covariance matrices for the process and measurement noise, respectively. Given an initial guess z 1 ∼ N (z 1 | µ 1|0 , Σ 1|0 ) the joint probability distribution of an LGSSM can be expressed using the LGSSM model (from (4)), Given observations x = {x t } T t=1 the mean and covariance of the state space variables can be estimated efficiently using a Kalman filter [11], µ t|t , Σ t|t , and a Kalman smoother [23], µ t|T , Σ t|T , see Algorithm 1 and Algorithm 2 in Supplementary B.

Kalman Variational Auto-Encoder
For the model we use a VAE to reduce the dimension of the image time series distribution wherein we assume that the dynamics can be represented linearly using a LGSSM in the latent space. An illustration of the model is given by Figure 1. The goal is to describe the dynamics of the system in the latent space with an LGSSM and use the decoder to reconstruct the image time series.
Illustration of a Kalman VAE. Similar to a regular VAE a lower-dimensional latent distribution is learnt from the data.
Similar to a regular VAE, the ELBO can be derived from the KL divergence between the approximate and real posterior distribution. For a single time series y = {y t } T t=1 the ELBO is given by (see Supplementary A.2 for derivation), where φ and θ are the encoder and decoder parameters, respectively, and γ = { A, C, R, Q, µ 1|0 , Σ 1|0 } are the LGSSM parameters. It is possible to sample (x,z) by first samplingx ∼ q θ (x | y) and then conditional samplingz ∼ p γ (z |x). Notice that p γ (z |x) is tractable using the Kalman smoother algorithm and the joint distribution p γ (x, z) is given by (5). With re-parameterization tricks [13] the model can be trained end-to-end to minimize the negative ELBO using e.g. stochastic gradient descent.

Dataset
In our experiments we use the EchoNet-Dynamic dataset [21]. This dataset includes 10, 030 2D echocardiogram time series where 7, 465 are allocated for training, 1, 277 for testing and 1, 288 for validation. We use a fixed time horizon of 50 time steps (T = 50) where each time step represents every second step in the original temporal resolution. This makes each sequence cover approximately two cardiac cycles. Sequences shorter than the fixed horizon are removed (245 in training, 40 in testing and 30 in validation). We use one sequence per time series and its start position is selected randomly. We downsample each image to an image size of 64 × 64 pixels and normalize the intensity to the range [−1, 1].

Implementation Details
The encoder consists of four 2-stride convolutional layers with [64, 128, 156, 512] filters respectively, where each layer is followed by batch normalization [9] and a leaky ReLU [17] activation. The output of the last convolutional layer is flattened and two dense layers are used-one for the mean, {µ enc t } T t=1 and one for the logarithm of the variances, {log σ enc t 2 } T t=1 . In the LGSSM we use a dimension of 16 for our observed state process, x ∈ R 16 , and 32 for the state space process, z ∈ R 32 . The decoder layers mirrors the encoder and we use a superresolution procedure [24], similar to the original Kalman VAE [5]. We assume the output distribution to be a multivariate Gaussian so that the log-likelihood is proportional to the mean-squared error, log p θ (y | x) ∝ ||y −ŷ|| 2 , whereŷ is the output from the decoder.
The model is trained end-to-end using importance sampling to maximize the ELBO in (6). We focus the training in the latent space by weighting the reconstruction term log p θ (y | x) with a factor 0.3 and we use KL cost annealing [3] for the KL divergence, log p θ (x)/q φ (x | y). Initially, in the first 5 epochs, we update only the VAE parameters, θ, φ. Then we include the LGSSM parameters, γ = [A, C, R, Q, µ 1|0 , Σ 1|0 ], in the training procedure and jointly update all parameters. For the LGSSM parameters we estimate the transition and emission matrices, the mean of the initial state µ 1|0 and the variances of the initial state Σ 1|0 , process-Q and measurement noise R, i.e. we assume the initial state and the noise to be uncorrelated. As optimizer we use Adam with exponential decay with factor 0.85 every 20 epoch and an initial learning rate of 10 −4 . We use a batch size of 16 image time series and train the model for 100 epochs on a single Nvidia GeForce GTX 1080 Ti graphic card. Our implementation uses Tensorflow and the code is publicly available 3 .

Result
The trained KVAE model was compared with a regular VAE of the same architecture as the KVAE. For comparison we calculated the structural similarity index (SSIM) [30] between the true validation image time series and the estimated/imputed ones. For the Kalman variational autoencoder we removed every 2th, 3th, 5th, 7th and 10th time sample and used Kalman smoothing and filtering to impute the missing values. Table 1 shows the overall result and in Figure 2 how the performance differs over time. In Figure 3 we show the smooth, p γ (x t | z 1:T ) = N x t | µ t|T , Σ t|T and filtered prediction,  p γ (x t | z 1:t−1 ) = N x t | µ t|t−1 , Σ t|t−1 distribution of the observed latent process for three scenarios: all observed values are known, imputation with every 3th sample observed and extrapolation for the last 5. An example of the true and imputed images based on every 5th sample is shown in Figure 4. The model is generative and as last example we generate a time series in the latent space by sample from the initial distribution N (µ 1|0 , Σ 1|0 ), extrapolate forward and decode the samples with the decoder. The generated image time series is shown in Figure 5.

Discussion
We have described an unsupervised method for extracting latent linear dynamics directly from a medical image time series, and shown how this linearity can be used to perform efficient filtering and smoothing. We found that imputation using smoothing led to marginally higher structural similarity indices, although it deteriorated for longer time horizons. We also demonstrated the ability to generate new time series from the model by randomly sampling states at t = 1 and extrapolating forward. However, in Figure 3 one can notice that the observed latent process is only locally linear and extrapolation over long horizons comes with a large uncertainty.
In related work [5,12] non-linear dynamics were handled by a neural network, e.g. several estimated state matrices where weighted based on the output of a recurrent dynamics parameter network [5]. In other approaches like extended   Kalman filter [10] and unscented Kalman filter [29] the linearity assumption can be removed, at the cost of approximate filtering and smoothing posteriors.
Comparisons with such models would be an interesting topic for future work.
Since we used a VAE with an explicit distribution, chosen to be multivariate Gaussian, the decoded images are quite blurry. Crisper reconstructions have been shown [16,8], in particular when VAEs have been combined with generative adversarial networks (GANs) [6].
By analyzing the transition matrix A, we have observed that the eigenvalues of A are within the unit circle-real, in fact-hence the system is stable.
While the focus of this project was to infer the dynamics based on the image intensity, another useful and interested topic would be the dynamics of the displacement field between two images. With known dynamics of the displacement fields, one could easily propagate auxiliary information, like structural regions, over time. This, we believe would have an impact in e.g. tracking. The displacement fields are, however, unknown in most cases and would have to be estimated using e.g. deformable image registration [7,25].