Dynamic imaging using a deep generative SToRM (Gen-SToRM) model

A. Dynamic MRI from undersampled data: problem setup

Our main focus is to recover a series of images x₁, ..x_M from their undersampled multichannel MRI measurements. The multidimensional dataset is often compactly represented by its Casoratti matrix

Each of the images is acquired by different multichannel measurement operators

where n_i is zero mean Gaussian noise matrix that corrupts the measurements.

B. Smooth manifold models for dynamic MRI

The smooth manifold methods model images x_i in the dynamic time series as points on a smooth manifold $\mathcal{M}$. These methods are motivated by continuous domain formulations that recover a function f on a manifold from its measurements as

where the regularization term involves the smoothness of the function on the manifold. This problem is adapted to the discrete setting to solve for images lying on a smooth manifold from its measurements as

where L is the graph Laplacian matrix. L is the discrete approximation of the Laplace-Beltrami operator on the manifold, which depends on the structure or geometry of the manifold. The manifold matrix L is estimated from k-space navigators. Different approaches, ranging from proximity-based methods [8] to kernel low-rank regularization [10] and sparse optimization [12], have been introduced.

The results of earlier work [10], [30] show that the above manifold regularization penalties can be viewed as an analysis prior. In particular, these schemes rely on a fixed non-linear mapping φ of the images. The theory shows that if the images x₁, ..x_M lie in a smooth manifold/surface or union of manifolds/surfaces, the mapped points live on a subspace or union of subspaces. The low-dimensional property of the mapped points φ(x₁), ..φ(x_M) is used to recover the images from undersampled data or derive the manifold using a kernel low-rank minimization scheme:

This nuclear norm regularization scheme is minimized using an iterative reweighted algorithm, whose intermediate steps match (4). The non-linear mapping φ may be viewed as an analysis operator that transforms the original images to a low-dimensional latent subspace, very similar to analysis sparsity-based approaches used in compressed sensing.

C. Unsupervised learning using Deep Image Prior

The recent work of DIP uses the structure of the network as a prior [32], enabling the recovery of images from ill-posed measurements without any training data. Specifically, DIP relies on the property that CNN architectures favor image data more than noise. The regularized reconstruction of an image from undersampled and noisy measurements is posed in DIP as

where $\mathbf{x}={\mathcal{G}}_{\mathbf{\theta }*}[\mathbf{z}]$ is the recovered image, generated by the CNN generator ${\mathcal{G}}_{\mathbf{\theta }*}$ whose parameters are denoted by θ. Here, z is the random latent variable, which is chosen as random noise and kept fixed.

The above optimization problem is often solved using stochastic gradient descent (SGD). Since CNNs are efficient in learning natural images, the solution often converges quickly to a good image. However, when iterated further, the algorithm also learns to represent the noise in the measurements if the generator has sufficient capacity, resulting in poor image quality. The general practice is to rely on early termination to obtain good results. This approach was recently extended to the dynamic setting by Jin et al. [33], where a sequence of random vectors was used as the input.

A. Generative model

We model the images in the time series as

where ${\mathcal{G}}_{\theta }$ is a non-linear mapping, which is termed as the generator. Inspired by the extensive work on generative image models [32], [36], [37], we represent ${\mathcal{G}}_{\theta }$ by a deep CNN, whose weights are denoted by θ. The parameters z_i are the latent vectors, which live in a low-dimensional subspace. The non-linear mapping ${\mathcal{G}}_{\theta }$ may be viewed as the inverse of the image-to-latent space mapping φ, considered in the SToRM approach.

We propose to estimate the parameters of the network θ as well as the latent vectors z_i by fitting the model to the undersampled measurements. The main distinction of our framework with DIP, which is designed for a single image, is that we use the same generator for all the images in the dynamic dataset. The latent vector z_i for each image is different and is also estimated from the measurements. This strategy allows us to exploit non-local information in the dataset. For example, in free-breathing cardiac MRI, the latent vectors of images with the same cardiac and respiratory phase are expected to be similar. When the gradient of the network is bounded, the output images at these time points are expected to be the same. The proposed framework is hence expected to learn a common representation from these time-points, which are often sampled using different sampling trajectories. Unlike conventional manifold methods [8], [10], [12], the use of the CNN generator also offers spatial regularization.

It is often impossible to acquire fully-sampled training data in many free-breathing dynamic imaging applications, and a key benefit of this framework over conventional neural network schemes is that no training data is required. As discussed previously, the number of parameters of the model in (7) is orders of magnitude smaller than the number of pixels in the dataset. The dramatic compression offered by the representation, together with the mini-batch training provides a highly memory-efficient alternative to current manifold based and low-rank/tensor approaches. Although our focus is on establishing the utility of the scheme in 2-D settings in this paper, the approach can be readily translated to higher dimensional applications. Another benefit is the implicit spatial regularization brought in by the convolutional network as discussed for DIP. We now introduce novel regularization priors on the network and the latent vectors to further constrain the recovery to reduce the need for manual early stopping.

B. Distance/Network regularization

As in the case of analysis SToRM regularization [8], [10], our interest is in generating a manifold model that preserves distances. Specifically, we would like the nearby points in the latent space to map to similar images on the manifold. With this interest, we now study the relation between the Euclidean distances between their latent vectors and the shortest distance between the points on the manifold (geodesic distance).

We consider two points z₁ and z₂ in the latent space, which are fed to the generator to obtain $\mathcal{G}\left({\mathbf{z}}_1\right)$ and $\mathcal{G}\left({\mathbf{z}}_2\right)$, respectively. We have the following result, which relates the the Euclidean distance ${\Vert {\mathbf{z}}_1-{\mathbf{z}}_2\Vert }^2$ to the geodesic distance ${\mathrm{dist}}_{\mathcal{M}}\left(\mathcal{G}\left({\mathbf{z}}_1\right),\mathcal{G}\left({\mathbf{z}}_2\right)\right)$, which is the shortest distance on the manifold. The setting is illustrated in Fig. 2, where the geodesic distance is indicated by the red curve.

C. Latent vector regularization penalty

The time frames in a dynamic time series have extensive redundancy between adjacent frames, which is usually capitalized by temporal gradient regularization. Directly penalizing the temporal gradient norm of the images requires the computation of the entire image time series, which is difficult when the entire image time series is not optimized in every batch.

We consider the norm of the finite differences between images specified by ${\Vert {\nabla }_p\mathbf{G}\left[{\mathbf{z}}_p\right]\Vert }^2$. Using Taylor series expansion, we obtain ${\nabla }_p\mathbf{G}\left[{\mathbf{z}}_p\right]=J_{\mathbf{z}}(\mathcal{G}[\mathbf{z}]){\nabla }_p\mathbf{z}+\mathcal{O}(p)$. We thus have

Since $J_{\mathbf{z}}(\mathcal{G}[\mathbf{z}])$ is small because of the distance regularization, we propose to add a temporal regularizer on the latent vectors. For example, when applied to free-breathing cardiac MRI, we expect the latent vectors to capture the two main contributors of motion: cardiac motion and respiratory motion. The temporal regularization encourages the cardiac and respiratory phases change slowly in time.

D. Proposed optimization criterion

Based on the above analysis, we derive the parameters of the network θ and the low-dimensional latent vectors z_i; i = 1, .., M from the measured data by minimizing:

with respect to z and θ. We use the ADAM optimization to determine the optimal parameters, and random initialization is used for the network parameters and latent variables.

A potential challenge with directly solving (13) is its high computational complexity. Unlike supervised neural network approaches that offer fast inference, the proposed approach optimizes the network parameters based on the measured data. This strategy will amount to a long reconstruction time when there are several image frames in the time series.

E. Strategies to reduce computational complexity

To minimize the computational complexity, we now introduce some approximation strategies.

1). Approximate data term for accelerated convergence:

When the data is measured using non-Cartesian sampling schemes, M non-uniform fast Fourier transform (NUFFT) evaluations are needed for the evaluation of the data term, where M is the number of frames in the dataset. Similarly, M inverse non-uniform fast Fourier transform (INUFFT) evaluations are needed for each back-propagation step. These NUFFT evaluations are computationally expensive, resulting in slow algorithms. In addition, most non-Cartesian imaging schemes over-sample the center of k-space. Since the least-square loss function in (5) weights errors in the center of k-space higher than in outer k-space regions, it is associated with slow convergence.

To speed up the intermediate computations, we propose to use gridding with density compensation, together with a projection step for the initial iterations. Specifically, we will use the approximate data term

instead of ${\sum }_i{\Vert {\mathcal{A}}_i\left(\mathcal{G}\left[{\mathbf{z}}_i\right]\right)-{\mathbf{b}}_i\Vert }^2$ in early iterations to speed up the computations. Here, g_i are the gridding reconstructions

where, $\mathcal{W}$ are diagonal matrices corresponding to multiplication by density compensation factors. The operators ${\mathcal{P}}_i$ in (14) are projection operators:

We note that the term $\left({\mathcal{A}}_i^H\mathcal{W}{\mathcal{A}}_i\right)\text{x}$ can be efficiently computed using Toeplitz embedding, which eliminates the need for expensive NUFFT and INUFFT steps. In addition, the use of the density compensation serves as a preconditioner, resulting in faster convergence. Once the algorithm has approximately converged, we switch the loss term back to (5) since it is optimal in a maximum likelihood perspective.

2). Progressive training-in-time:

To further speed up the algorithm, we introduce a progressive training strategy, which is similar to multi-resolution strategies used in image processing. In particular, we start with a single frame obtained by pooling the measured data from all the time frames. Since this average frame is well-sampled, the algorithm promptly converges to the optimal solution. The corresponding network serves as a good initialization for the next step. Following convergence, we increase the number of frames. The optimal θ parameters from the previous step are used to initialize the generator, while the latent vector is initialized by the interpolated version of the latent vector at the previous step. This process is repeated until the desired number of frames is reached.

This progressive training-in-time approach significantly reduces the computational complexity of the proposed algorithm. In this work, we used a three-step algorithm. However, the number of steps (levels) of training can be chosen based on the dataset. This progressive training-in-time approach is illustrated in Fig. 3.

A. Structure of the generator

The structure of the generator used in this work is given in Table. I. The output images have two channels, which correspond to the real and imaginary parts of the MR images. Note that we have a parameter d in the network. This user-defined parameter controls the size of the generator or, in other words, the number of trainable parameters in the generator. We also have a number ℓ(z) as a user-defined parameter. This parameter represents the number of elements in each latent vector. In this work, it is chosen as ℓ(z) = 2 as we have two motion patterns in cardiac images. We use leaky ReLU for all the non-linear activations, except at the output layer, where it is tanh activation.

B. Datasets

This research study was conducted using data acquired from human subjects. The Institutional Review Board at the local institution (The University of Iowa) approved the acquisition of the data, and written consents were obtained from all subjects. The experiments reported in this paper are based on datasets collected in the free-breathing mode using the golden angle spiral trajectory. We acquired eight datasets on a GE 3T scanner. One dataset was used to identify the optimal hyperparameters of all the algorithms in the proposed scheme. We then used the hyperparameters to generate the experimental results for all the remaining datasets reported in this paper. The sequence parameters for the datasets are: TR = 8.4 ms, FOV = 320 mm× 320 mm, flip angle = 18°, slice thickness = 8 mm. The datasets were acquired using a cardiac multichannel array with 34 channels. We used an automatic algorithm to pre-select the eight best coils, that provide the best signal to noise ratio in the region of interest. The removal of the coils with low sensitivities provided improved reconstructions [39]. We used a PCA-based coil combination using SVD such that the approximation error < 5%. We then estimated the coil sensitivity maps based on these virtual channels using the method of Walsh et al. [40] and assumed they were constant over time.

For each dataset in this research, we binned the data from six spiral interleaves corresponding to 50 ms temporal resolution. If a Cartesian acquisition scheme with TR = 3.5ms were used, this would correspond to ≈14 lines/frame; with a 340×340 matrix, this corresponds roughly to an acceleration factor of 24. Moreover, each dataset has more than 500 frames. During reconstruction, we omit the first 20 frames in each dataset and use the next 500 frames for SToRM reconstructions; this is then used as the simulated ground truth for comparisons. The experiments were run on a machine with an Intel Xeon CPU at 2.40 GHz and a Tesla P100-PCIE 16GB GPU. The source code for the proposed GenSToRM scheme can be downloaded from this link: https://github.com/qing-zou/Gen-SToRM.

C. Quality evaluation metric

In this work, the quantitative comparisons are made using the Signal-to-Error Ratio (SER) metric (in addition to the standard Peak Signal-to-Noise Ratio (PSNR) and the Structural Similarity Index Measure (SSIM)) defined as:

Here x_orig and x_recon represent the ground truth and the reconstructed image. The unit for SER is decibel (dB).

The SER metric requires a reference image, which is chosen as the SToRM reconstruction with 500 frames. However, we note that this reference may be imperfect and may suffer from blurring and related artifacts. Hence, we consider the Blind/referenceless Image Spatial Quality Evaluator (BRISQUE) [41] to evaluate the score of the image quality. The BRISQUE score is a perceptual score based on the support vector regression model trained on an image database with corresponding differential mean opinion score values. The training image dataset contains images with different distortions. A smaller score indicates better perceptual quality.

D. State-of-the-art methods for comparison

We compare the proposed scheme with the recent state-of-the-art methods for free-breathing and ungated cardiac MRI. We note that while there are many deep learning algorithms for static MRI, those methods are not readily applicable to our setting.

• Analysis SToRM [9], [10], published in 2020: The manifold Laplacian matrix is estimated from k-space navigators using kernel low-rank regularization, followed by solving for the images using (4).

• Time-DIP [33] implementation based on the arXiv form at the submission of this article: This is an unsupervised learning scheme, that fixes the latent variables as noise and solves for the generator parameters. For real-time applications, Time-DIP chooses a preset period, and the noise vectors of the frames corresponding to the multiples of the period were chosen as independent Gaussian variables [33]. The latent variables of the intermediate frames were obtained using linear interpolation. We chose a period of 20 frames, which roughly corresponds to the period of the heart beats.

• Low-rank [2]: The image frames in the time series are recovered using the nuclear norm minimization.

E. Hyperparameter tuning

We used one of the acquired datasets to identify the hyperparameters of the proposed scheme. Since we do not have access to the fully-sampled dataset, we used the SToRM reconstructions from 500 images (acquisition time of 25 seconds) as a reference. The smoothness parameter λ of this method was manually selected as λ = 0.01 to obtain the best recovery, as in the literature [9]. All of the comparisons relied on image recovery from 150 frames (acquisition time of 7.5 seconds). The hyperparameter tuning approach yielded the parameters d = 40, λ₁ = 0.0005, and λ₂ = 2 for the proposed approach. We demonstrate the impact of tuning d in Fig. 6, while the impact of choosing λ₁ and λ₂ is shown in Fig. 4. The hyperparameter optimization of SToRM from 150 frames resulted in the optimal smoothness parameter λ = 0.0075. For Time-DIP, we follow the design of the network shown by Jin et al. [33], where the generator consists of multiple layers of convolution and upsampling operations. To ensure fair comparison, we used a similar architecture, where the base size of the network was tuned to obtain the best results.

We use a three-step progressive training strategy. In the first step, the learning rate for the network is 1 × 10⁻³ and 1000 epoches are used. For the second step of training, the learning rate for the network is 5 × 10⁻⁴ and the learning rate for the latent variable is 5×10⁻³. In this stage, 600 epoches are used. In the final step of training, the learning rate for the network is 5×10⁻⁴, the learning rate for the latent variable is 1×10⁻³, and 700 epoches are used.

A. Impact of different regularization terms

We first study the impact of the two regularization terms in (13). The parameter d corresponding to the size of the network (see Table I) was chosen as d = 24 in this case. In Fig. 4 (a), we plot the reconstruction performance with respect to the number of epoches for three scenarios: (1) using both regularization terms; (2) using only latent regularization; and (3) using only distance/network regularization. In the experiment, we use 500 frames of SToRM (~ 25 seconds of acquisition) reconstructions, which is called “SToRM500”, as the reference for SER computations. We tested the reconstruction performance for the three scenarios using 150 frames, which corresponds to around 7.5 seconds of acquisition. From the plot, we observe that without using the network regularization, the SER degrades with increasing epoches, which is similar to that of DIP. In this case, an early stopping strategy is needed to obtain good recovery. The latent vectors corresponding to this setting are shown in (c), which shows mixing between cardiac and respiratory waveforms. When latent regularization is not used, we observe that the SER plot is roughly flat, but the latent variables show quite significant mixing, which translates to blurred reconstructions. By contrast, when both network and latent regularizations are used, the algorithm converges to a better solution. We also note that the latent variables are well decoupled; the blue curve captures the respiratory motion, while the orange one captures the cardiac motion. We also observe that the reconstructions agree well with the SToRM reconstructions. The network now learns meaningful mappings, which translate to improved reconstructions when compared to the reconstructions obtained without using the regularizers.

B. Benefit of progressive training-in-time approach

In Fig. 5, we demonstrate the significant reduction in runtime offered by the progressive training strategy described in Section III-E2. Here, we consider the recovery from 150 frames with and without the progressive strategy. Both regularization priors were used in this strategy, and d was chosen as 24. We plot the reconstruction performance, measured by the SER with respect to the running time. The SER plots show that the proposed scheme converges in around ≈ 200 seconds, while the direct approach takes more than 2000 seconds. We also note from the SER plots that the solution obtained using progressive training is superior to the one without progressive training.

C. Impact of size of the network

The architecture of the generator ${\mathcal{G}}_{\theta }$ is given in Table I. Note that the size of the network is controlled by the user-defined parameter d, which dictates the number of convolution filters and hence the number of trainable parameters in the network. In this section, we investigate the impact of the user-defined parameter d on the reconstruction performance. We tested the reconstruction performance using d = 8, 16, 24, 32, 40, and 48, and the obtained results are shown in Fig. 6. From the figure, we see that when d = 8 or d = 16, the generator network is too small to capture the dynamic variations. When d = 8, the generator is unable to capture both cardiac motion and respiratory motion. When d = 16, part of the respiratory motion is recovered, while the cardiac motion is still lost. The best SER scores with respect to SToRM with 500 frames is obtained for d = 24, while the lowest Brisque scores are obtained for d = 40. We also observe that the features including papillary muscles and myocardium in the d = 40 results appear sharper than those of SToRM with 500 frames, even though the proposed reconstructions were only performed from 150 frames. We use d = 40 for the subsequent comparisons in the paper.

D. Comparison with the state-of-the-art methods

In this section, we compare our proposed scheme with several state-of-the-art methods for the reconstruction of dynamic images.

In Fig. 7, we compare the region of interest for SToRM500, SToRM with 150 frames (SToRM150), the proposed method with two different d values, the unsupervised Time-DIP approach, and the low-rank algorithm. From Fig. 7, we observe that the proposed scheme can significantly reduce errors in comparison to SToRM150. Additionally, the proposed scheme is able to capture the motion patterns better than Time-DIP, while the low-rank method is unable to capture the motion patterns. From the time profile in Fig. 7, we notice that the proposed scheme is capable of recovering the abrupt change in blood-pool contrast between diastole and systole. This is due to inflow effects associated with gradient echo (GRE) acquisitions. In particular, the blood from regions outside the slice enters the heart, which did not experience any of the former slice-selective excitation pulses; the differences in magnetization of the blood with no magnetization history, and that was within the slice, results in the abrupt change in intensity. We note that some of the competing methods such as Time-DIP and low-rank, blur these details.

We also perform the comparisons on a different dataset in Fig. 8. We compare the proposed scheme with SToRM500, SToRM150, Time-DIP, and the low-rank approach. The results are shown in Fig. 8. From the figure, we see that the proposed reconstructions appear less blurred than those of the conventional schemes.

We also compared the proposed scheme with SToRM500, SToRM150, and the unsupervised Time-DIP approach quantitatively. We omit the low-rank method here because low-rank approach often failed in some datasets. The quantitative comparisons are shown in Table II. We used SToRM500 as the reference for SER, PSNR, and SSIM calculations. The quantitative results are based on the average performance from six datasets.

Finally, we illustrate the proposed approaches in Fig. 9 and Fig. 10, respectively. The proposed approach decoupled the latent vectors corresponding to the cardiac and respiratory phases well, as shown in the representative examples in Fig. 9 (a) and Fig. 10 (a).