Patent ID: 12190483

DETAILED DESCRIPTION

FIG.1is a block diagram of one example of an imaging system. The imaging system may include a scanner12. The scanner may further include a sequence controller14, a gantry16, and a bed18. An operation unit20is shown which includes a computing device22. The operation unit20is configured to receive measurements24from the scanner12. The computing device22may be configured to use the measurements24in determining parameters36and latent variables28to use in a generative manifold model30. The generative manifold model30may be represented by a series of instructions which may be executed by the computing device22. One or more input devices40may be operatively connected to the storage unit20. One or more displays42may be operatively connected to the operations unit. The display42may be used to display images based on the measurements24and application of the generative manifold model30. Thus, dynamic images of three or more dimensions may be displayed on the display42. Storage44is operatively connected to the operation unit which may be used to store representations of the measurements24or dynamic imagery.

FIG.2is one example of a method. In step50a set of measurements are acquired such as may be acquired using a magnetic resonance imaging (MRI) scanner. In step52, parameters of a generator are estimated. In step54, latent variables are estimated. In step56, modeling is performing using a generator to map to a manifold. In step58results may be displayed and/or stored.

FIG.3further illustrates on example of a computing device22. The computing device22may have a memory60and a processor62. The processor62may be configured to execute machine readable instructions64which include instructions for implementing a generative manifold model30.

The present disclosure sets forth several different sets of examples. The first is an embodiment showing dynamic imaging using a deep generative StoRM (Gen-SToRM) model. The second is an embodiment showing dynamic imaging using Motion-Compensated SmooThness Regularization on Manifolds (MoCo-SToRM). The third is an embodiment showing variational manifold learning from incomplete data: application to multislice dynamic MRI.

Part I: Dynamic Imaging Using a Deep Generative SToRM (Gen-SToRM) Model

I. Introduction

The imaging of time-varying objects at high spatial and temporal resolution is key to several modalities, including MRI and microscopy. A central challenge is the need for high resolution in both space and time [1], [2]. Several computational imaging strategies have been introduced in MRI to improve the resolution, especially in the context of free-breathing and ungated cardiac MRI. A popular approach pursued by several groups is self-gating, where cardiac and respiratory information is obtained from central k-space regions (navigators) using bandpass filtering or clustering [3]-[7]. The data is then binned to the respective phases and recovered using total variation or other priors. Recently, approaches using smooth manifold regularization have been introduced. These approaches model the images in the time series as points on a high-dimensional manifold [8]-[12]. Manifold regularization algorithms, including the smoothness regularization on manifolds (SToRM) framework [8]-[10], have shown good performance in several dynamic imaging applications. Since the data is not explicitly binned into specific phases as in the self-gating methods, manifold algorithms are less vulnerable to clustering errors than navigator-based corrections. Despite the benefits, a key challenge with the current manifold methods is the high memory demand. Unlike self-gating methods that only recover specific phases, manifold methods recover the entire time series. The limited memory on current GPUs restricts the number of frames that can be recovered simultaneously, which makes it challenging to extend the model to higher dimensionalities. The high memory demand also makes it difficult to use spatial regularization priors on the images using deep learned models.

One main focus is to capitalize on the power of deep convolutional neural networks (CNN) to introduce a memory efficient generative or synthesis formulation of SToRM. CNN based approaches are now revolutionizing image reconstruction, offering significantly improved image quality and fast image recovery [13]-[19]. In the context of MRI, several novel approaches have been introduced [20], [21], including transfer learning [22], domain adaptation [23], learning-based dynamic MRI [24]-[26], and generative-adversarial models [27]-[29]. Unlike many CNN-based approaches, the proposed scheme does not require pre-training using large amounts of training data. This makes the approach desirable in free-breathing applications, where the acquisition of fully sampled training data is infeasible. We note that the classical SToRM approach can be viewed as an analysis regularization scheme (seeFIG.4A). Specifically, a non-linear injective mapping is applied on the images such that the mapped points of the alias-free images lie on a low-dimensional subspace [10], [30], [31]. When recovering images from undersampled data, the nuclear norm prior is applied in the transform domain to encourage their non-linear mappings to lie in a subspace. Unfortunately, this analysis approach requires the storage of all the image frames in the time series, which translates to high memory demand. The proposed generative SToRM formulation offers quite significant compression of the data, which can overcome the above challenge. Both the relation between the analysis and synthesis formulations and the relation of the synthesis formulation to neural networks were established in earlier work [31]. We assume that the image volumes in the dataset are smooth non-linear functions of a few latent variables, i.e., xt=θ(zt), where ztare the latent vectors in a low-dimensional space. xtis the t-th generated image frame in the time series. This explicit formulation implies that the image volumes lie on a smooth non-linear manifold in a high-dimensional ambient space (seeFIG.4B). The latent variables capture the differences between the images (e.g., cardiac phase, respiratory phase, contrast dynamics, subject motion). We model the G using a CNN, which offers a significantly compressed representation. Specifically, the number of parameters required by the model (CNN weights and latent vectors) are several orders of magnitude smaller than required for the direct representation of the images. The compact model proportionately reduces the number of measurements needed to recover the images. In addition, the compression also enables algorithms with much smaller memory footprint and computational complexity. We propose to jointly optimize for the network parameters θ and the latent vector zt, based on the given measurements. The smoothness of the manifold generated byθ(z) depends on the gradient ofθwith respect to its input. To enforce the learning of a smooth image manifold, we regularize the norm of the Jacobian of the mapping ∥Jzθ∥2. We experimentally observe that by penalizing the gradient of the mapping, the network is encouraged to learn meaningful mappings. Similarly, the images in the time series are expected to vary smoothly in time. Hence, we also use a Tikhonov smoothness penalty on the latent vectors ztto further constrain the solutions. We use the ADAM optimizer with stochastic gradients, where random batches of ztand biare chosen at iteration to determine the parameters. Unlike traditional CNN methods that are fast during testing/inference, the direct application of this scheme to the dynamic MRI setting is computationally expensive. We use approximations, including progressive-intime optimization and an approximated data term that avoids non-uniform fast Fourier transforms, to significantly reduce the computational complexity of the algorithm.

The proposed approach is inspired by deep image prior (DIP), which was introduced for static imaging problems [32], as well as its extension to dynamic imaging [33]. The key difference of the proposed formulation is the joint optimization of the latent variables z and G. The work of Jin et al. [33] was originally developed for CINE MRI, where the latent variables were obtained by linearly interpolating noise variables at the first and last frames. Their extension to real-time applications involved setting noise latent vectors at multiples of a preselected period, followed by linearly interpolating the noise variables. This approach is not ideally suited for applications with free breathing, when the motion is not periodic. Another key distinction is the use of regularization priors on the network parameters and latent vectors, which encourages the mapping to be an isometry between latent and image spaces. Unlike DIP methods, the performance of the network does not significantly degrade with iterations. While we call our algorithm “generative SToRM”, we note that our goal is not to generate random images from stochastic inputs as in generative-adversarial networks (GAN). In particular, we do not use adversarial loss functions where a discriminator is jointly learned as in the literature [34], [35].

II. Background

A. Dynamic MRI From Undersampled Data: Problem Setup

Our main focus is to recover a series of images x1, . . . xMfrom their undersampled multichannel MRI measurements. The multidimensional dataset is often compactly represented by its Casoratti matrix
X=[x1, . . . xM].  (1)

Each of the images is acquired by different multichannel measurement operators
bi=Ai(xi)+ni,  (2)

where niis zero mean Gaussian noise matrix that corrupts the measurements.

B. Smooth Manifold Models for Dynamic MRI

The smooth manifold methods model images xiin the dynamic time series as points on a smooth manifold M. These methods are motivated by continuous domain formulations that recover a function ƒ on a manifold from its measurements as

f=arg⁢minf∑if⁡(xi)-bi2+λ⁢∫M∇Mf2⁢d⁢x(3)

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

X=arg⁢minX∑i=1MA⁡(xi)-bi2+λ⁢trace(XL⁢XH),(4)

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 x1, . . . xMlie 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 φ(x1), . . . φ(xM) is used to recover the images from undersampled data or derive the manifold using a kernel low-rank minimization scheme:

X*=arg⁢minX∑i=1MA⁡(xi)-bi2+λ⁢[φ⁡(x1),…⁢φ⁡(xN)]*.(5)

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 illposed 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

{θ*}=arg⁢minθA⁡(x)-b2⁢such⁢that⁢x=𝒢θ[z](6)

where x=θ*[z] is the recovered image, generated by the CNN generatorθ* 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.

III. Deep Generative Storm Model

We now introduce a synthesis SToRM formulation for the recovery of images in a time series from undersampled data (seeFIG.4B). Rather than relying on a non-linear mapping of images to a low-dimensional subspace [10] (seeFIG.4A), we model the images in the time series as non-linear functions of latent vectors living in a low-dimensional subspace.

FIG.4AandFIG.4Billustrate Analysis SToRM and Generative SToRM, respectively. Analysis SToRM considers a non-linear (e.g. exponential) lifting of the data. If the original points lie on a smooth manifold, the lifted points lie on a low dimensional subspace. The analysis SToRM cost function in (5) is the sum of the fit of the recovered images to the undersampled measurements and the nuclear norm of the lifted points. A challenge with analysis SToRM is its high memory demand and the difficulty in adding spatial regularization. The proposed method models the images as the non-linear mappingθof some latent vectors zi, which lie in a very low-dimensional space. Note that the same generator is used to model all the images in the dataset. The number of parameters of the generator and the latent variables is around the size of a single image, which implies a highly compressed representation. In addition, the structure of the CNN offers spatial regularization as shown in DIP. The proposed algorithm in (13) estimates the parameters of the generator and the latent variables from the measured data. A distance regularization prior is added to the generator to ensure that nearby points in the latent subspace are mapped to nearby points on the manifold. Similarly, a temporal regularization prior is added to the latent variables. The optimization is performed using ADAM with batches of few images.

A. Generative Model

We model the images in the time series as
xi=θ(zi),i=1, . . . ,M,(7)

whereθ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θby a deep CNN, whose weights are denoted by zi. The parameters ziare the latent vectors, which live in a low-dimensional subspace. The non-linear mappingθ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 ziby 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 zifor 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 z1and z2in the latent space, which are fed to the generator to obtain(z1) and(z2), respectively. We have the following result, which relates the Euclidean distance ∥z1−z2∥2to the geodesic distance distM((z1),(z2)) which is the shortest distance on the manifold. The setting is illustrated inFIG.5, where the geodesic distance is indicated by the red curve.

Proposition 1. Let z1, z2∈nbe two nearby points in the latent space, with mappings denoted by(z1),(z2)∈M. Here, M={G(z)|z∈n}. Then, the geodesic distance on the manifold satisfies:
distM((z1),(z2))≤∥z1−z2∥F∥Jz((z1))∥F.  (8)

Proof. The straight-line between the latent vectors is denoted by c(s), s∈[0,1] with c(0)=z1and c(1)=z2. We also assume that the line is described in its curvilinear abscissa, which implies ∥c′(s)∥=1; ∀s∈[0,1]. We note thatmay map to the black curve, which may be longer than the geodesic distance. We now compute the length of the black curve[c(s)] as
d=∫01∥∇s[c(s)]∥ds.(9)

Using the chain rule and denoting the Jacobian matrix ofby Jz, we can simplify the above distance as

d=∫01Jz(𝒢)⁢c′(s)F⁢ds≤∫01Jz(𝒢)F⁢c′(s)F︸1⁢ds=Jz(𝒢[z1])F⁢∫01ds︸z1-z2.(10)

We used the Cauchy-Schwartz inequality in the second step and in the last step, we use the fact that Jz()(c(t))=Jz(z1)+(t) when the points z1and z2are close. Since the geodesic distance is the shortest distance on the manifold, we have distM((z1),(z2))≤d and hence we obtain (8).

The result in (8) shows that the Frobenius norm of the Jacobian matrix ∥Jz∥ controls how far apartmaps two vectors that are close in the latent space. We would like points that are close in the latent space map to nearby points on the manifold. We hence use the gradient of the map:
Rdistance=∥Jz((z))∥F2(11)

as a regularization penalty. We note that the above penalty will also encourage the learning of a mappingsuch that the length of curve(c(t)) is the geodesic distance. We note that the above penalty can also be thought of as a network regularization. Similar gradient penalties are used in machine learning to improve generalization ability and to improve the robustness to adversarial attacks [38]. The use of gradient penalty is observed to be qualitatively equivalent to penalizing the norm of the weights of the network.

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 ∥∇pG[zp]∥2. Using Taylor series expansion, we obtain ∇pG[zp]=Jz([z])∇pz+(p). We thus have
∥∇p[zp]∥≈∥Jz([z])∇pz∥≤∥Jz([z])∥∥∇pz∥.(12)

Since Jz([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 zi; i=1, . . . , M from the measured data by minimizing:

C⁡(z,θ)=∑i=1NAi(Gθ[zi])-b2︸data⁢term+λ1⁢Jz⁢𝒢θ[zi]2︸distance⁢regularization+λ2⁢∇tzt2︸latent⁢regularization(13)

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
D(z,θ)=Σi=1M|i(θ[zi]−gi)∥2(14)

instead of Σi=1∥Ai([zi]−bi)∥2in early iterations to speed up the computations. Here, giare the gridding reconstructions
gi=(AiHAi)†AiHbi≈AiHb,(15)

where,are diagonal matrices corresponding to multiplication by density compensation factors. The operators Piin (14) are projection operators:
Pix=(AiHAi)†(AiHAi)x≈(AiHAi)x(16)

We note that the term (AiHAi) 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 inFIG.6.

IV. Implementation Details and Datasets

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 tan h activation.

TABLE 1ARCHITECTURE OF THE GENERATORθ·(z) MEANS THE NUMBER OFELEMENTS IN EACH LATENT VECTOR.Input sizefilter sz# filtersPaddingStrideOutput size1 × 1 ×(z)1 × 1100011 × 1 × 1001 × 1 × 1003 × 38d013 × 3 × 8d3 × 3 × 8d3 × 38d015 × 5 × 8d5 × 5 × 8d4 × 44d1210 × 10 × 4d10 × 10 × 4d4 × 44d1220 × 20 × 4d20 × 20 × 4d3 × 34d0241 × 41 × 4d41 × 41 × 4d5 × 52d1285 × 85 × 2d85 × 85 × 2d4 × 4d12170 × 170 × d170 × 170 × d4 × 4d12340 × 340 × d340 × 340 × d3 × 3212340 × 340 × 2
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 3 T 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:5 ms 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 16 GB GPU. The source code for the proposed Gen-SToRM 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:

SER=20·log1⁢0⁢xo⁢r⁢i⁢gxorig-xrecon.

Here xorigand xreconrepresent 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, λ1=0.0005, and λ2=2 for the proposed approach. We demonstrate the impact of tuning dinFIG.9, while the impact of choosing λ1and λ2is shown inFIG.7. 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−3and 1000 epoches are used. For the second step of training, the learning rate for the network is 5×10−4and the learning rate for the latent variable is 5×10−3. In this stage, 600 epoches are used. In the final step of training, the learning rate for the network is 5×10−4, the learning rate for the latent variable is 1×10−3, and 700 epoches are used.

V. Experiments and Results

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.7A, 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 inFIG.7C, 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

InFIG.8, 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θ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 inFIG.9. 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.

InFIG.10, 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. FromFIG.10, 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.

TABLE IIQUANTITATIVE COMPARISONS BASEDON SIX DATASETS: WE USED SIXDATASETS TO OBTAIN THE AVERAGESER, PSNR, SSIM, BRISQUE SCORE, ANDTIME USED FOR RECONSTRUCTION.MethodsSToRM500SToRM150PropsedTime-DIPSER (dB)NA17.318.216.7PSNR (dB)NA32.733.532.0SSIMNA0.860.890.87Brisque35.240.237.142.9Time (min)47131757

From the time profile inFIG.10, 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 inFIG.11. We compare the proposed scheme with SToRM500, SToRM150, Time-DIP, and the low-rank approach. The results are shown inFIG.11. 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 inFIGS.12A-12EandFIGS.13A-13F, respectively. The proposed approach decoupled the latent vectors corresponding to the cardiac and respiratory phases well, as shown in the representative examples inFIG.12AandFIG.13A.

VI. Conclusion

In this work, we introduced an unsupervised generative SToRM framework for the recovery of free-breathing cardiac images from spiral acquisitions. This work assumes that the images are generated by a non-linear CNN-based generatorθ, which maps the low-dimensional latent variables to high resolution images. Unlike traditional supervised CNN methods, the proposed approach does not require any training data. The parameters of the generator and the latent variables are directly estimated from the undersampled data. The key benefit for this generative model is its ability to compress the data, which results in a memory-effective algorithm. To improve the performance, we introduced a network/distance regularization and a latent variable regularization. The combination of the priors ensures the learning of representations that preserve distances and ensure the temporal smoothness of the recovered images; the regularized approach provides improved reconstructions while minimizing the need for early stopping. To reduce the computational complexity, we introduced a fast approximation of the data loss term as well as a progressive training-in-time strategy. These approximations result in an algorithm with computational complexity comparable to our prior SToRM algorithm.

The main benefits of this scheme are the improved performance and considerably reduced memory demand. While our main focus in this work was to establish the benefits of this work in 2D, we plan to extend this work to 3D applications in the future.

Part II: Dynamic Imaging Using Motion-Compensated SmooThness Regularization on Manifolds (MoCo-SToRM)

1. Introduction

Magnetic resonance imaging (MRI) is an attractive imaging modality for patient groups that require serial follow up because it does not use ionizing radiation. Ultra-short echo-time MRI [42] methods are capable of significantly reducing the T2* losses, mitigating some of the main challenges associated with lung MRI. However, MRI is a slow imaging modality, which makes it challenging to image moving organs such as the lung. For lung MRI, the respiratory motion can be frozen by breath-holds. However, subjects usually are unable to hold their breath for a long time, which will significantly limit the achievable spatial resolution and coverage. Besides, there are several patient groups (e.g., patients with chronic obstructive pulmonary disease (COPD), pediatric patients, and neonates) who cannot hold their breath even for a short duration [43,44,45]. For these reasons, some of these patients need to be sedated for MRI exams or are often not eligible for MRI scans.

Several gating approaches were introduced to eliminate the need for breath-holding in pulmonary MRI [46,47,48,49]. For instance, classical methods (e.g., [50]) rely on respiratory bellows or self-gating signals to bin the data to specific phases. Prospective methods only acquire the data during a specific respiratory phase, while retrospective methods continuously acquire the data but only use the data from a specific phase. Self-gating approaches such as XD-GRASP [4] use the information from the central k-space samples to estimate the motion signal, which is used to bin the acquired data into several respiratory phases. After the binning, a compressed-sensing approach is used to jointly reconstruct the phase images (SeeFIG.14Afor illustration), which is more data efficient than traditional binned acquisitions. These approaches are often called motion-resolved methods. A challenge with these methods is the potential sensitivity to bulk motion during the scan. In particular, the subjects may move abruptly during the scan. Because XD-GRASP and similar gating methods rely on low-pass filtering to estimate the pseudo-periodic motion signal, the bulk motion effects are often filtered out. In addition, these approaches are only able to recover the respiratory phase images, which correspond to the averaged data over several minutes, and not the true dynamics. Another challenge associated with the motion-resolved scheme is the trade-off between residual aliasing and blurring resulting from intra-bin motion. For instance, increasing the number of bins can reduce intra-bin motion artifacts. However, this will come at the expense of k-t space data available for each bin, which will translate to residual alias artifacts. While manifold approaches [9,11,51,52], which perform soft-gating as opposed to explicit binning, offer improved trade-offs but are also vulnerable to these challenges. These schemes use machine learning algorithms to perform soft-binning of the data using the manifold structure of images in the dataset. These unsupervised machine learning methods have been shown to offer improved performance and robustness to different motion patterns over explicit binning strategies.

Motion compensation (MoCo) is often used to further improve the data efficiency and to reduce residual aliasing and noise in the reconstructed images. Many of the approaches require a high-resolution reference image. The recovered images are then registered to the reference image to obtain the motion fields [53]. Another approach is to estimate the motion-maps between the phase images reconstructed by XD-GRASP; the different motion phases are registered together and averaged to obtain a MoCo volume [48]. Recently, Zhu et al. used a non-linear compressed sensing algorithm to directly recover the MoCo volume from the k-t space data, using the deformation maps estimated from XD-GRASP [54]. This approach, named iMoCo, is shown to significantly improve the image quality. The main challenge with this multi-step strategy (binning based on motion estimation, followed by XD-GRASP reconstruction, followed by the final MoCo reconstruction) is the dependence of the image quality on the intermediate steps. In particular, this approach inherits the sensitivity of XD-GRASP to bulk motion because it is dependent on motion estimates from XD-GRASP. For example, when the data is corrupted with a single bulk motion effect, the data with bulk motion is removed in [54] to obtain good reconstructions; this approach is not readily applicable to settings where there are multiple bulk motion events during the acquisition. This multi-step strategy was recently extended to dynamic PET [55]. The main distinction of this scheme from [54] is the use of a deep learning algorithm to estimate the deformation maps between the different motion phases. The different motion phases are registered to a fixed state using a deep network, followed by a reconstruction scheme similar to iMoCo [54].

The main focus of this work is to introduce a novel unsupervised deep-learning MoCo reconstruction scheme, which can be readily applied for free-breathing pulmonary MRI. This method is the generalization of the previous motion-resolved generative manifold methods [51,52] to the MoCo setting; we hence call the proposed approach motion-compensated smoothness regularization on manifolds (MoCo-SToRM). Unlike [51,52], which assume the images to be on a smooth manifold, we assume that the motion deformation maps at different time instants are living on a manifold, parameterized by low-dimensional latent vectors. We assume the deformation maps to be the output of a convolutional neural network (CNN) based generator, whose inputs are time-dependent low-dimensional latent vectors that capture the motion information (SeeFIG.14Bfor illustration). The generated deformation maps are used to deform a learned template image, which corresponds to the image volume frame in the time series. A multi-channel non-uniform Fourier transform (NUFFT) is used to generate the k-space measurements of the images. Unlike prior MoCo approaches that use a series of algorithms for binning, reconstruction, motion estimation, and reconstruction, we formulate the joint recovery of the latent vectors, deformation maps, and the template image directly from the measured k-t space data as a single non-linear optimization scheme. The cost function is the squared error between the multi-channel Fourier measurements of the image volumes and the actual measurements acquired from the specific subject. We note that the deformation maps are smooth and are less complex than the images themselves; we expect the proposed scheme to be less data-demanding than motion-resolved approaches [4,51,52].

The proposed framework, built using modular blocks, is highly explainable. In particular, the learned latent vectors capture the intrinsic temporal variability in the time series, including respiratory and bulk motion as seen from our experiments. Moreover, the smooth deformation maps capture the spatial deformation of the image template and can be visualized. More importantly, the reconstructions can be viewed as a movie, allowing one to visualize the images at respective time-frames, unlike binning-based approaches that only recover the phase images. Unlike current deep-learning strategies that pre-learn the CNN parameters from example data, the proposed scheme learns all the parameters from the data of the specific subject. We also note that the MoCo approach enables us to minimize the trade-off between intra-phase motion and the data available for reconstruction. We do not make any assumptions on the latent vectors, which allows it to learn all the motion events during the acquisition, including bulk motion, which is challenging for traditional methods.

2. Methods

2.1. Brief Background on Motion-Resolved and Motion-Compensated Reconstruction

Several self-navigated motion-resolved free-breathing MRI schemes [4,56], which use 3D radial ultra-short-echo (UTE) acquisition, were recently introduced for lung imaging. These schemes rely on a combination of low-pass filtering and/or clustering to derive the self-gating signals, which are used to bin the k-space data into different motion phases. Once the data is binned, these schemes (e.g., XD-GRASP [4]) perform the motion-resolved joint reconstruction of the phases by solving the following:

F*=arg⁢minF𝒜⁡(F)-B22+λs⁢Ψ⁢F1+λt⁢T⁢V⁡(F).

The first term is the data consistency term that compares the multi-channel measurements of the phase images F={f1, . . . , fN} with the binned data B={b1, . . . , bN}. Here N is the number of bins and hence the number phases in the XD-GRASP reconstruction. The second term is a spatial sparsity e1-wavelet penalty term, in which ψ is the wavelet transform. The third term is the total variation penalty along the motion phases. The above scheme is usually called XD-GRASP-type motion-resolved reconstruction, which is illustrated inFIG.14A.

Deep manifold-based approaches [51,57] offer an alternate route for motion-resolved recovery. In particular, the images in the time series are modeled as ft=θ(zt), whereθis a deep CNN generator that is shared across different image frames. The parameters of the generator denoted by θ and the low-dimensional latent vectors ztare learned such that the cost function Σt∥tθ(zt)−bt∥2is minimized. Here,tdenotes the forward model corresponding to the tthimage frame, and btare the corresponding k-space measurements. The proposed approach is a generalization of these deep manifold models [51, 57] to the MoCo setting.

Many of the early MoCo methods rely on a high-resolution static reference image [53]. Recent approaches [48,54] rely on motion-resolved XD-GRASP reconstructions. Once the motion-resolved reconstructions are obtained, one of the phases (usually the exhalation phase) will be chosen as the reference. Then the other motion phases are registered to the reference phase to obtain the deformations Φl, . . . Φp. The iMoCo approach solves the following optimization scheme to obtain the MoCo reconstruction fimoco:

fi⁢m⁢o⁢c⁢o*=argf⁢mini⁢m⁢o⁢c⁢o∑p=1N𝒜t(ϕp(fi⁢m⁢o⁢c⁢o))-bp22+λ⁢T⁢G⁢V⁡(fi⁢m⁢o⁢c⁢o).

The first term is the data consistency term, and the second is a spatial total generalized variation sparsity regularizer [58].

2.2. Proposed Approach

In this work, we extend the motion-resolved deep manifold methods [51,57] to the MoCo reconstruction setting. The proposed framework is a end-to-end self-supervised deep-learning algorithm involving explainable learning modules. We model the image volumes ft in the time series as the deformed versions of a single image template f:
ft(x,y,z)=f(x−ϕx(t),y−ϕy(t),z−ϕz(t)):=(f,ϕ(t)).  (17)

Here ϕ(t)={ϕx(t), ϕy(t), ϕz(t)} is the motion/deformation map at the time instant t. We implement D as a differentiable interpolation layer.

We propose to jointly estimate the deformation maps and the single image template directly from the k-t space data. We note that it is impossible for us to acquire the k-t space data at the Nyquist-Shannon sampling rate. Therefore, the joint estimation problem is highly ill-posed. In order to regularize the deformation maps, we use the manifold assumption. In other words, we assume that the deformation map for each image frame ft is living on a smooth manifold, parameterized by low-dimensional latent vectors zt, t=1, . . . M that capture the dynamics (e.g., respiratory motion, bulk motion in lung imaging). We model the non-linear mapping between the low-dimensional latent vectors and the high-dimensional deformation maps by a CNN generator:
ϕ(t)=θ(zt),  (18)

whose input is the low-dimensional latent vector zt∈Rd. As the dominant motion in free-breathing lung imaging is the respiratory motion, we set d=1 in this work; we will consider higher dimensional latent space in our future work. Combining (17) and (18), each image frame ft in the time series is modeled as:

ft(r)=𝒟⁢(f,𝒢θ(zt)︸)ϕ⁡(t).(19)

Here, r=(x, y, z) is the spatial coordinate. SeeFIG.14Bfor an illustration.

We note that the generated image at each time instant t is dependent on the image template f, the parameters of the deep CNN generator θ, and the low-dimensional latent vectors z=[zl, . . . , zM]. Here, M is the number of image frames in the time series. When golden-angle or bit-reversed radial acquisitions are used, M can be a user-defined parameter. We propose to jointly solve for the above unknowns directly from the k-t space data of the specific subjects as the optimization problem:
C(z,θ,f)=Σt=1M∥t(ft−bt∥2+λ1∥∇tz∥l1+λ2∥∇rf∥l1,  (20)

where ftis related to the static image f by (19). Here, Atare the forward operators that are performed on each of the time points. We implement Atas multi-channel NUFFT [59] operators using the k-space trajectory at the tthtime instant using the SigPy package [60]. btare the multi-channel k-space measurements acquired from the subject. The second term in (20) is a smoothness penalty on the latent vector that captures the dynamics (e.g., respiratory and/or bulk motion). If this term is not added, the learned latent vectors will learn high-frequency oscillations. To minimize this risk, we added a total-variation penalty on the latent vector z along the time direction to encourage the latent vectors to learn piecewise smooth motion. The last term in (20) is the spatial total variation penalty on the static image, which enables us to further reduce alias artifacts in the learned static image.

The proposed self-supervised scheme offers several benefits. First of all, the reconstruction relies only on the undersampled data acquired from the specific subject. Unlike most deep-learning strategies, the proposed framework does not require fully sampled training datasets, which are not available in our setting, to train the networks. Secondly, the proposed scheme does not require physiological monitors such as respiratory belts or dedicated k-space navigators. It also eliminates the need for band-pass filtering or clustering to estimate the phase information, which will filter out bulk motion. Finally, unlike binned approaches that recover average images over the acquisition duration within respective respiratory phases, the proposed scheme enables the recovery of the natural dynamics of the lung

2.3 Approaches to Minimize Computational Complexity

We use ADAM optimization [61] with a batch size of one time-frame to find the optimal z, θ and f. The small memory footprint enables us to use this scheme for high-resolution 3D+time problems. The network and optimization scheme was implemented in PyTorch. The motion generator is implemented using an eight-layer network. The first seven layers are 3D convolutional layers with 200 features per layer. The last layer is an up-sampling layer, which uses tri-linear interpolation to interpolate the deformation maps from lower resolution to high resolution. The final interpolation step allows us to account for the prior knowledge that the deformation maps are smooth functions. ReLU activation function [62] is used for all the convolutional layers.

Directly solving the above optimization problem (20) is computationally expensive, especially when the image resolution is high. To minimize the computational complexity, we use a progressive strategy. In particular, we first solve for (20) for very low-resolution images using the corresponding region in the central k-t space. These latent vectors, motion fields, and the images are used to initialize the network at a higher spatial resolution. We use two progressive steps to refine the resolution, until the final resolution is reached. We observe that this progressive strategy significantly improves the convergence rate. In particular, few iterations are needed at the highest resolution, compared to the setting where the parameters of the motion network and the image are initialized randomly.

The above joint optimization strategy offers good estimates of the latent vectors with few iterations, even at the lowest resolution. By contrast, the stochastic gradients with a batch size of one can result in low convergence rates for the static image f and the CNN parameters. To further accelerate the convergence rate, we additionally use a binning strategy similar to motion-resolved schemes shown inFIG.14A, assuming z to be fixed. We bin the latent vector to P phases, based on the latent vectors we estimated. We use 25 phases in the adult subjects with less extensive motion and 150 in the neonatal intensive care unit (NICU) patient with extensive motion. This approach allows us to bin the data to different phases. We update θ and f using the optimization strategy:

C⁡(θ,f)=∑p=1P𝒜p⁢𝒟⁡(f,𝒢θ(zp))︸fp-bp2+λ2⁢∇rfℓ1,(21)

Here, zpis the latent vector at the pthbin, and ϕ(t)=θ(zp) is the motion vector for the pthbin. Here fpis the image in the pthbin, obtained by deforming f with ϕp.

We only use the above binning-based optimization in (21) at the highest resolution level. Note that (21) only solves for θ and f, and not z. We hence alternate between (21) and (20) at the highest resolution.

3. Datasets and Evaluation

3.1. Experimental Datasets

The datasets used in the experiments in this work were acquired using an optimized 3D UTE sequence with variable-density readouts to oversample the k-space center [63]. We used four datasets acquired from two adult subjects. Two of them are from a healthy subject (pre-contrast and post-contrast) and another two are from a fibrotic subject (pre-contrast and post-contrast). We also used one dataset acquired from a female subject with severe bronchopulmonary dysplasia (BPD), who was admitted to the NICU. The gestational age of the patient at birth is 24 weeks, and the MRI is prescribed at the chronological age of 15 weeks and 3 days. The weight of the patient at MRI is 3.18 Kg. The datasets from the healthy subject were acquired on a 1.5 T GE scanner using 8 coils. The variable-density readouts help retain signal-to-noise ratio (SNR), and oversampling reduces aliasing artifacts. A bit-reversed ordering was used during the data acquisition. The prescribed field of view (FOV)=32×32×32 cm3. The matrix size is 256×256×256. The data were acquired with 90K radial spokes with TR≈3.2 ms and 655 samples/readout, corresponding to an approximately five-minute acquisition.

The datasets from the diseased subject were acquired on a 3T GE scanner using 32 coils. The prescribed FOV=32×32×32 cm3. The matrix size is 256×256×256. The data were acquired with 91K radial spokes with TR≈2.8 ms and 655 samples/readout, corresponding to an approximately four-minute acquisition. When we were processing the datasets with 32 coils, we used a PCA-based coil combination [64] using SVD to keep only eight virtual coils.

The dataset from the neonatal subject was acquired on a 1.5 T small footprint Mill scanner [44] located in the NICU. The data was acquired using the built-in body coil, which translates to poor SNR compared to the adult scans. In addition, the scan is made challenging because of the extensive bulk motion by the subject during the scan. For this dataset, the prescribed FOV=18×18×18 cm3. The matrix size is 256×256×256. The data was acquired with 200K radial spokes with TR≈5 ms and 1013 samples/readout, corresponding to an approximately 16-minute acquisition.

This research study was conducted using human subject data. Approval was granted by the Ethics Committees of the institutions where the data was acquired.

3.2. Numerical Phantom to Validate MoCo-SToRM

We note that we do not have accurate ground truth datasets to evaluate the quantitative accuracy of the proposed scheme and its potential impact to bulk motion. We hence constructed a high-resolution numerical phantom using the XD-GRASP and iMoCo reconstructions of the pre-contrast dataset from a healthy subject. Specifically, we registered the XD-GRASP exhalation phase to the inhalation phase to obtain the deformation maps. Then we modulated the deformation maps by a periodic triangular function with a specific frequency and a DC off-set to simulate the motion from the exhalation phase to the inhalation phase. To study the impact of bulk motion, we also consider additional random translational motion (move about 2-5 pixels), at random time points as shown inFIG.16. We deform the static iMoCo reconstruction with the above deformation maps to generate the high-resolution images at different time points. The multi-channel NUFFT of these images were used as the measurements, with an additive white Gaussian noise of 0.5%. The creation of the simulation data is illustrated inFIG.15.

3.3. Figures of Merit for Quantitative Evaluation

For image quality comparisons, we compare the proposed MoCo-SToRM reconstruction with XD-GRASP and iMoCo. To quantitatively compare the image quality, we use three images metrics in this work.Diaphragm maximum derivative (DMD) [54,65]: the DMD will be used to measure the sharpness of the lung-liver diaphragm. It is defined as:

D⁢M⁢D=Max(∂I)Mean⁢(Il⁢i⁢v⁢e⁢r),where Max(∂I) is the maximum intensity change between the lung-liver interface, which is computed by choosing the maximum value of the image gradient. Mean(Iliver) is the mean intensity in the chosen liver region. A higher DMD implies sharper edges.Signal-to-noise ratio (SNR) [66,67,68]: The SNR is computed as

SNR=20⁢log⁡(μsσn),where μsis the mean of the intensity of the chosen region of interest and σnis the standard deviation of the intensity of a chosen noise region. A higher SNR usually means better image quality. In our study, we manually choose the regions of interest.Contrast-to-noise ratio (CNR) [54,69]: The CNR is computed as

CNR=20⁢log⁡(μA-μBσn),where μAand μBare the mean of the intensity of two regions within the region of interest and σnis the standard deviation of the intensity of a chosen noise region. The higher CNR usually means better image quality. In our study, we manually choose the regions of interest.

4. Results

4.1. Numerical Simulation Experiments

The results of the simulation study are shown inFIG.16. In this simulation study, we investigate the impact of bulk motion events. Specifically, we create the simulation data with no, two, four, and ten bulk motion events. We quantitatively compare the reconstructions using the metrics of PSNR, SSIM, relative error of the reconstruction, and relative error of the deformation maps. The quantitative results of the simulation study are summarized inFIG.16A. InFIG.16BandFIG.16C, we show some results of the simulation study. We show the learned latent vectors, the time profiles of the reconstructed image volumes, and the deformation maps. The comparisons of the reconstructions are also shown in the figure. From the simulation study, we see that the proposed MoCo-SToRM approach works reasonably when there are four bulk motion events. By contrast, when there are ten bulk motion events, the performance of the proposed scheme degrades in a graceful fashion. In particular, the limited number of radial k-space spokes translates to imperfect estimation of motion, which in turn translates to blurry reconstructions.

4.2. Experimental Datasets

InFIG.17, we show the learned latent vectors, time profiles, and example estimated deformation maps and corresponding time profiles from the pre-contrast dataset from the healthy subject. We show the estimated latent vectors from the first 200 frames inFIG.17A. We also show the time profile of the reconstructed images inFIG.17Band the time profile of the deformation maps inFIG.17C, corresponding to the blue lines in the images. From the two profiles, we see that the motion patterns coincide with the learned latent vectors. InFIG.17B, we show the estimated deformation maps from two time points, indicated by red and green dots inFIG.17A, corresponding to the inhalation phase and the exhalation phase. The results show that the latent vectors closely capture the dynamics of the motion.

4.3. Comparison with State-of-the-Art Methods

In this section, we compare the results of the proposed scheme with XD-GRASP and iMoCo. InFIG.18, we show the visual comparisons of the methods on post-contrast datasets. FromFIG.18, we observe that the MoCo-SToRM reconstructions can reduce the noise and capture more details when compared to the motion-resolved XD-GRASP reconstructions. Furthermore, the MoCo-SToRM reconstructions are less blurred than those of the motion-compensated iMoCo reconstructions. We note that the post-contrast dataset from the diseased subject had a bulk motion event (seeFIG.18), which translates to blurred iMoCo reconstructions.

The quantitative comparison of the proposed scheme with the competing methods on four datasets (two from healthy adult subjects and two from diseased adult subjects) are shown inFIG.19. We first measure the DMD on 15 sagittal slices in each dataset, and the quantitative results are shown inFIG.19A. From the DMD results, we see that the proposed MoCo-SToRM scheme is able provide comparable results. In particular, the motion-compensated methods (i-MoCO and MoCO-SToRM) are observed to yield marginally higher DMD than XD-GRASP, implying reduced blurring, as shown inFIG.19A.

In addition to DMD, we also report the SNR and CNR of the aortic arch and lung parenchyma region, as shown inFIG.19B. We see that the proposed MoCo-SToRM results are comparable with the results obtained from the motion-compensated iMoCo.

4.4. Impact of Bulk Motion

The acquisition time in pulmonary MRI is around four minutes for adult subjects and around 16 minutes for NICU subjects. The relatively long scan time makes current approaches vulnerable to bulk motion artifacts, especially during the imaging of diseased and pediatric patients. If they are not compensated, these bulk motion errors translate to residual blurring.

In cases with significant bulk motion, existing methods use additional image-based approaches to detect and reject sections of data with bulk motion [54]. These approaches are readily applicable to cases with very few bulk motion events; for instance, a case with a single bulk motion event was considered in [54]. If multiple events are in the dataset, this approach may severely restrict the available k-t space data and hence translate to significantly degraded image quality.

An advantage of the proposed MoCo-SToRM scheme is its ability to directly account for bulk motion during the scan. In the proposed MoCo-SToRM scheme, the latent vectors and the CNN-based generator have the ability to capture the bulk motion in the data and account for it during the reconstruction. In particular, we observe sudden jumps in the learned latent vectors. The non-linear nature of the learned generator allows the generation of deformation maps for each motion state, depending on the value of the latent vectors.

4.4.1. Adult study: In the post-contrast dataset acquired from a diseased subject, we detected one bulk motion, which is shown inFIG.20AandFIG.20B. InFIG.20A, we show the latent vectors that are learned using MoCo-SToRM. As we mentioned before, the sudden jump in the latent vectors indicates the bulk motion. We highlight two time regions, one without bulk motion (red box) and one with bulk motion (green box). We also show the time profiles of the yellow line indicated inFIG.20A. From the plots of the time profiles, we see that when the latent vectors have no sudden jump, then no bulk motion can be seen. However, when sudden jump happens in the latent vectors, we clearly see a bulk motion event from the time profile. In (I) ofFIG.20A, we zoomed the reconstructed image corresponding to the yellow cross shown in the latent vectors, and in (II) ofFIG.20A, we zoomed the reconstructed image corresponding to the purple cross shown in the latent vectors. From the red line marker, we can see that the subject moved the shoulder during the scan.

InFIG.20B, we also show the deformation maps estimated by the proposed scheme at four different time points indicated by the brown, yellow, blue, and purple crosses inFIG.20A. From the red ellipses in the images, we see that the subject moved the shoulder during the scan. The deformation map will be very different from the deformation maps when there is no bulk motion. InFIG.18B, we compared the reconstructions from the proposed scheme and the iMoCo and XD-GRASP reconstructions. From the figure, we see that the proposed scheme is able to deal with bulk motion and offer improved reconstructions.

4.4.2. Feasibility study in neonatal imaging: Neonatal subjects often suffer from several developmental lung disorders; the non-ionizing nature of MRI radiation makes it the ideal modality to image the lung of these subjects [43,44,45]. A low-field, low-footprint MRI system was considered in [43,45]. By imaging the neonatal subjects within the NICU, this approach minimizes the risk of infection. The subject was imaged using a body coil, which offers limited SNR compared to the multi-coil array used in the adult setting. One of the main challenges with neonatal MRI is bulk motion, especially when the subjects are awake. In this work, we study the feasibility of the proposed scheme to offer motion compensation in a challenging subject with extensive bulk motion, which was challenging for the conventional methods. We show the results inFIG.21AtoFIG.21C. InFIG.21A, we show the motion signal (latent vectors) estimated from the proposed MoCo-SToRM scheme. We note that there are several discontinuities in the latent signal, which correspond to bulk motion events. Two of the bulk motion events are highlighted in the red and purple boxes, respectively. Besides the two examples of bulk motion, we also zoom into a section with no bulk motion. We also show the reconstructed images corresponding to the marked positions in the latent vectors in each of the sub-series in the respective boxes. The bulk motions in the two examples can be clearly seen from the boundary of the body, denoted by the red dotted curve. Furthermore, we see that when the patient was in different positions, the shape of the lung was different, as indicated by the yellow dotted curves. InFIG.21B, we show the comparisons with iMoCo and XD-GRASP based on two slices from the axial view. From the figure, we see that MoCo-SToRM is able to reduce the motion artifacts, as indicated by the red arrows in the figures. Furthermore, we can see that the boundaries in iMoCo and XD-GRASP are blurred out due to bulk motions, as shown by the yellow arrows. However, MoCo-SToRM is able to reconstruct the boundaries. Also, we can see that some details in the lung region are captured in the MoCo-SToRM reconstructions, as indicated by the green arrows.FIG.21Cshows the sagittal view of the comparisons.

4.5. Maximum Intensity Projections of the Reconstructions

In this section, we show some results that we obtained from the proposed MoCo-SToRM reconstructions. InFIG.22, we show the reconstructions obtained from two post-contrast datasets, one from a healthy subject and anther from a diseased subject. Maximum intensity projection (MIP) [70] is used to generate the results. MIP is known to have the benefit that the vascular structures can be clearly seen as tubular and branching structures in MIP images [71]. By showing the three views for the reconstruction using MIP, the lung structure and vascular structures for each subject can be seen in a direct way, which can be readily used by doctors in clinics. For each sub-figure inFIG.22, 20 slices are used for MIP images.

More results, including some movies on the reconstructions and showcases of the bulk motions in both the adult subject and the neonatal subject, can be found on our website: https://sites.google.com/view/qing-zou/blogs/moco-storm.

5. Discussion & Conclusion

In this work, we proposed an unsupervised motion-compensated scheme using smoothness regularization on manifolds for the reconstruction of high-resolution free-breathing lung MRI. The proposed algorithm jointly estimates the latent vectors that capture the motion dynamics, the corresponding deformation maps, and the reconstructed motion-compensated images from the raw k-t space data of each subject. Unlike current motion-resolved strategies, the proposed scheme is more robust to bulk motion events during the scan, which translates to less blurred reconstructions in datasets with extensive motion. The proposed approach may be applicable to pediatric and neonatal subjects that are often challenging to image using traditional approaches. In this study, we restricted our attention to 1-D latent vectors. In our future work, we will consider its extension using higher-dimensional latent space, which will allow improved robustness to different motion components, including cardiac motion and bulk motion. The challenge with the direct extension of the proposed scheme to this setting is the increased computational complexity.

A difference between the proposed scheme and the motion-resolved reconstruction is that rather than just resolve some phases, the proposed scheme is able to get the temporal-resolved reconstruction with 0.1 s temporal resolution. This means that we are able to have more intermediate motion states (˜25 states) between the exhalation state and the inhalation state. The proposed MoCo-SToRM scheme is also able to deal with bulk motions as discussed in the previous section. This offers the possibility of using the proposed scheme for patient groups such as pediatric patients.

Part III: Variational Manifold Learning From Incomplete Data: Application to Multislice Dynamic MRI

I. INTRODUCTION

Deep generative models [72] that rely on convolutional neural networks (CNNs) are now widely used to rep-resent data living on nonlinear manifolds. For instance, the variational autoencoder (VAE) [73] represents the data points as CNN mappings of the latent vectors, whose parameters are learned using the maximum likelihood formulation. Since the exact log-likelihood of the data points is intractable, VAE relies on the maximization of a lower bound of the likelihood, involving an approximation for the conditional density of the latent variable represented by an encoder neural network. The

VAE framework offers several benefits over the vanilla autoencoder [72], including improved generalization [74] and ability to disentangle the important latent factors [76], [76]. Unfortunately, most of the current generative models are learned from fully sampled datasets. Once learned, the probability density of the data can be used as a prior for various applications, including data imputation [77, 78]. Unfortunately, fully-sampled datasets are often not available in many high-resolution structural and dynamic imaging applications to train autoencoder networks.

The main focus of this paper is to introduce a variational framework to learn a deep generative manifold directly from undersampled/incomplete measurements. The main application motivating this work is the multislice free-breathing and ungated cardiac MRI. Breath-held CINE imaging, which pro-vides valuable indicators of abnormal structure and function, is an integral part of cardiac MRI exams. Compressed sensing [79]-[82] and deep learning methods have emerged as powerful options to reduce the breath-hold duration, with excellent performance [83]-[87]. Despite these advances, breath-held CINE imaging is challenging for several subject groups, including pediatric and chronic obstructive pulmonary disease (COPD) subjects. Several authors have introduced self-gating [3, 4, 5, 39, 88, 89] and manifold approaches [11]-[12], [90]-[92] to enable free-breathing and ungated single-slice cardiac MRI. For instance, the smoothness regularization on manifolds (SToRM) ap-proach [8]-[10] models the images as points on a low-dimensional manifold whose structure is exploited using a kernel low-rank formulation [8], [10] to recover the im-ages from highly undersampled measurements. Recently, deep learning-based manifold models were introduced [51, 52, 57] to further improve the performance; these schemes learn a deep generative network and its latent variables directly from the measured k-space data using a non-probabilistic formulation.

All of the previously described free-breathing cardiac MRI reconstruction approaches (e.g., compressed sensing-based approaches, manifold approaches, and deep learning-based approaches) independently recover the data from each slice. Cardiac Mill often relies on slice-by-slice acquisition to preserve myocardium to blood pool contrast, resulting from the in-flow of blood from unexcited regions to the slice of interest; the improved contrast facilitates segmentation. The above-mentioned 2D self-gating and manifold methods are thus unable to exploit the extensive redundancies between adjacent slices, which could offer improved performance. Note that the respiratory and cardiac motion during the acquisition of the different slices is often very different; this makes the direct 3D extension of the 2D self-gating and manifold methods impossible. Another challenge with the approaches mentioned above is the need for post-processing methods to determine matching slices at specific cardiac/respiratory phases for estimation of cardiac parameters (e.g., ejection fraction, strain). Several post-processing methods have been introduced to align the data post reconstruction [91], [93]-[96]. Because these methods require fully sampled data, they will not facilitate the exploitation of the inter-slice redundancies during image recovery.

We introduce a novel variational framework for the joint recovery and alignment of multislice data from the entire heart. This approach combines the undersampled k-t space data from different slices, possibly acquired with multiple cardiac and respiratory motion patterns, to recover the 3D dynamic MRI dataset. We use a 3D CNN generative model, which takes in a latent vector and outputs a 3D image volume. The time-varying latent vectors capture the intrinsic variability in the dataset, including cardiac and respiratory motion. The latent variables and the parameters of the 3D CNN are jointly learned from the multislice k-t space data using a maximum likelihood formulation. Since the likelihood is not tractable, we maximize its variational lower bound involving a model for the conditional distribution of the latent variables, which is conceptually similar to the VAE approach [73]. The VAE scheme uses an encoder network to derive the conditional probabilities of the latent vectors from fully sampled data [73]. This approach is not directly applicable in our setting because each data sample is measured using a different measurement operator. We hence model the conditional densities as a Gaussian distribution whose parameters are learned from the undersampled data directly using back-propagation. We use a Gaussian prior on the latent variables while deriving the evidence-based lower bound (ELBO); the Gaussian prior ensures that the latent variables from different slices have similar distributions, facilitating the alignment of the slices. We note that the direct extension of our previous generative manifold model [51], [57] to the 3D setting does not have any constraint on the latent variables; this extension results in poor alignment of the slices and degradation in image quality in the 3D setting. We also use smoothness priors on the latent variables to further improve the performance. Once learned, the representation can be used to generate matching 3D volumes with any desired cardiac/respiratory phase by exciting the generator with appropriate latent vectors. This approach of learning a generative model of the entire heart may thus be viewed as a paradigm shift from conventional slice-by-slice image-recovery algorithms [3]-[5], [88]-[92], [8]-[12], [39].

II. BACKGROUND ON DYNAMIC MRI

A. Multislice Free-Breathing MRI: Problem Statement

The main application considered in this paper is the recovery of 3D cardiac volumes of the heart from undersampled 2D multislice k-t space data acquired in the free-breathing and ungated setting. In particular, we consider the recovery of the time series x(r, tz), where r=(x, y, z) represents the spatial coordinates and tzdenotes the time frame during the acquisition of the zthslice. We model the acquisition of the
b(tz)=Atz(x(r,tz))+ntz,  (22)

where b(tz) is the k-t space data of the zthslice at the tthtime frame. Here, Atzare the time-dependent measurement operators, which evaluate the multi-channel single-slice Fourier measurements of the 3D volume x(r, tz) on the trajectory Ktzcorresponding to the time point t. Specifically, tzextracts the zthslice from the volume x(r, tz) and evaluates its single-slice measurements. ntzrepresents the noise in the measurements.

B. CNN-Based Generative Manifold Models in Dynamic MRI

CNN-based generative models were recently introduced for single-slice dynamic MRI [51]. This scheme models the 2-D images in the time series as the output of a CNN generator Dθ:
xi=Dθ(ci),i=1, . . . ,M.

The input ciis the latent vector, which lives in a low-dimensional subspace. The recovery of the images in the time series involves the minimization of the criterion

C⁡(c,θ)=∑i=1NAi(Dθ(ci))-bi2︸data⁢intern+λ1⁢Jc⁢Dθ(c)2︸net⁢reg.+λ2⁢∇ici2︸latent⁢reg..(23)

The first term in the cost function is a measure of data consistency, while the second term is a network regularization term that controls the smoothness of the generated manifold [51]. The last term is the temporal smoothness of the latent variables, which is used to further improve the performance.

III. VARIATIONAL MANIFOLD LEARNING

We now introduce a novel variational formulation to learn a manifold from undersampled measurements, which is the generalization of the seminal VAE approach [73] to the under-sampled setting. We will first present the proposed approach in a simple and general setting for simplicity and ease of understanding. The use of this variational manifold model for the joint alignment and recovery of 3D images from 2-D multislice MRI data will be described in Section IV.

A. General Problem Statement and Intuition

We assume that the images in the time series, indexed by i, live on a smooth manifold M and hence can be modeled as the output of a CNN-based generator:
xi=Dθ(ci),  (24)

where ciis the low-dimensional latent variable corresponding to xi. Here, θ denotes the weights of the generator, which is shared for all the images.

Most generative models consider the learning of the above model from fully sampled data. By contrast, we consider the recovery from incomplete measurements
bi=Ai(xi)+ni(25)

Here is an undersampled measurement operator corresponding to the ithimage frame. Here, ni∈(0, σ2I) are noise vectors. Note that the measurement operators for each xiare different. If the same sampling operators are used for all the data points, it is impossible to recover the images without additional prior information. We assume that the sampling operators satisfy the following properties:

We assume Aito be a rectangular sub-matrix, obtained by picking specific rows of an orthonormal measurement operator (e.g., Fourier transform).

We assume that the measurement operators A˜S are drawn from a distribution and satisfy
A˜S[ATA]=I,(26)

which is the identity operator. The above condition guarantees diversity in the measurement operators.

We now provide some intuition about why the learning of the model with the above settings will succeed under the restrictive assumptions on the measurement operators described above. In the noiseless setting, we consider the learning of the latent variables ciand the weights θ by minimizing the empirical error:

{θ*,ci*}=arg⁢minθ,ci∑iAi(xi-𝒟θ(ci))2.︸ℒ(27)

Here, xiare the fully sampled data points. When A˜S, this empirical sum approximates

ℒ=𝔼x~M⁢𝔼A~S⁢𝒜⁡(x-𝒟θ(c))2=𝔼x~M⁢𝔼A~S⁢〈x-𝒟θ(c),AH⁢A⁡(x-𝒟θ(c))〉=𝔼x~M⁢〈x-𝒟θ(c),𝔼A~S[AH⁢A]︸I⁢(x-𝒟θ(c))〉=argminθ,c𝔼x~M⁢x-𝒟θ(c)2.

The above result follows from (26) and the orthonormality of the full measurement operator. This result shows that the recovery of the true manifold is feasible from undersampled data when the sampling operators satisfy the properties listed above.

B. Proposed Algorithm

We consider the recovery of the images xi from their measurements (25) by maximizing their likelihood, specified by

p⁡(bi)=p⁡(bi,ci)p⁡(ci|bi)(28)

We note that the posterior p(ci|bi) is not tractable. Following the VAE approach in [73], we use a surrogate distribution to approximate p(ci|bi). The VAE formulation uses an encoder network to model p(ci|xi) from the fully sampled data (bi=xi). Unfortunately, this approach is not directly applicable in our setting since bi is the undersampled data, measured using Aithat vary with i.

We propose to use a Gaussian model qi(ci)≈p(ci|bi), parameterized by its mean μiand diagonal covariance matrix Σi, and to estimate these parameters using back-propagation. Following a similar argument as in [73], we show in the Appendix that the likelihood term in (28) can be lower-bounded as

log⁢p⁡(bi)≥-12⁢σ2⁢𝔼ci~qi(ci)[Ai⁢𝒟θ(ci)-bi2]︸data⁢term⁢KL[qi(ci)⁢❘"\[LeftBracketingBar]"❘"\[RightBracketingBar]"⁢p⁡(ci)]︸L⁡(qi):latent⁢regularization.(29)

Here, p(ci) is a prior on the latent variables. In this work, we assume p(ci)=(0, I), where I is the identity matrix. In this case, the KL divergence can be explicitly evaluated as

L⁡(ci)=-log[det(∑)]-n+trace⁢(∑)+μT⁢μ2,

where we assume a latent space of dimension n.

We hence solve for the unknown weights of the generator θ as well as the parameters of qidenoted by μiand Σiby minimizing the negative of the lower bound in (29).

Following [73], we use a Monte-Carlo approach to approximate the expectation in the data term. In particular, at each epoch of the training loop, we derive the samples cias
ci=μi+Σiϵ,  (30)

where ϵ is a zero-mean unit variance Gaussian random variable. At each iteration, the estimation process thus involves the minimization of the criterion

C(θ,{μi,∑i}︸qi)=∑i=1N⁢d⁢a⁢t⁢a(𝒜i⁢𝒟θ(ci)-bi2+σ2⁢L⁡(qi)),(31)

with respect to the unknowns θ, μiand Σi.

C. Illustration Using MNIST Data

We provide a simple example for the illustration of the above variational model from undersampled data of the digit 1 in the MNIST dataset [97]. The images used are scaled to the range [−1, 1].

The generator we used here is a simple CNN with three layers. ReLU activation function is used for the first two layers and tan h is used for the last layer. The dimension of the latent space is chosen as 2. In this example, all the trainable parameters are initialized as small random numbers, and the hyper-parameter for the latent regularization L(ci) is chosen as 1. We used 1,000 epoches to train the CNN generator.

We first trained the model from the fully sampled data (i=I), whose results are shown in the first row ofFIG.23. Then we trained the model from undersampled noisy data. In the example, 70% of the pixel values in each image are missing, while Gaussian white noise with standard deviation 0.05 is added to the known 30% pixel values. The recovered images are shown in the second row ofFIG.23. We report the peak signal-to-noise ratio (PSNR) and the structural similarity index measure (SSIM) for the results.

IV. APPLICATION TO DYNAMIC MRI

We first describe the application of the algorithm in the single-slice free-breathing and ungated data, which is the setting considered in [51]. We then generalize the approach to the novel setting of the joint alignment and recovery of 3D MRI from multislice free-breathing data in Section IV-C.

A. Acquisition Scheme and Pre-Processing of Data

The datasets used in this work are acquired using a 2D (GRE) sequence with golden angle spiral readouts in the free-breathing and ungated setting on a MR750W scanner (GE Healthcare, Waukesha, WI, USA). The sequence parameters for the datasets are: FOV=320 mm×320 mm, flip angle=18°, slice thickness=8 mm. The datasets were acquired using a cardiac multi-channel array with 34 channels. The Institutional Review Board at the University of Iowa approved the acquisition of the data, and written consents were obtained from the subjects. The number of slices acquired for different subjects varies.

We used an algorithm developed in house to pre-select the coils that provide the best signal-to-noise ratio in the region of interest. A PCA-based coil combination scheme was then used such that the approximation error was less than 5%. We then estimated the coil sensitivity maps based on these virtual channels using ESPIRiT [98] and assumed them to be constant over time.

A total of 3,192 spirals were acquired for each slice in the subjects with TR=8.4 ms, which corresponds to an acquisition time of 27 seconds. Among the 3,192 spirals, every sixth spiral was acquired with the same angle; these spirals were used for self-navigation in the reconstruction methods that require self-navigation. We binned the data from six spiral interleaves corresponding to 50 ms temporal resolution for each frame.

B. Single-Slice Variational SToRM Algorithm

Based on the analysis in the previous sections, we use the following scheme for the recovery of single-slice dynamic MRI. We use a re-parameterization layer to obtain the latent variables c(t) from the time-varying probability distributions q(c(t)) with parameters μtand Σt. These latent variables are fed to the CNN generatorθ, which generates the reconstructed volumes x(t)=θ(c(t). The multi-channel, non-uniform, Fourier transform-based forward operators are applied on the reconstructed images, which are then compared to the actual noisy measurements bi. The illustration of this scheme is shown inFIG.24A. The parameters in the generator and the μiand Σiare updated based on the loss function
(θ,{μt,Σt})=C(θ,{μt,Σt})+λ1∥θ∥12+λ2∥∇μt∥2.   (32)

Here, C(θ, {μt, Σt}) is defined in (31), which is the lower bound for maximum likelihood estimation. The second term in (32) is a regularization penalty on the generator weights. It has been shown in [51] that adding this term makes the training of the decoder more stable. The third term involves the temporal gradients of the latent vectors, which enforces the latent vectors to capture the smooth nature of motion patterns in the dynamic images. We use the ADAM optimization to determine the optimal parameters. We also adopt the progressive-in-time training strategy introduced in [51] to realize a computationally efficient reconstruction. We term this dynamic MRI reconstruction scheme as single-slice variational SToRM.

C. Multislice Variational SToRM Algorithm

We now generalize the single-slice variational SToRM scheme for the joint alignment and recovery of multislice dynamic MRI. We assume that the image volume at the time point t during the acquisition of the zthslice, denoted by x(r, tz), as the output of the generator:
x(r,tz)=θ(c(tz)).

Here, c(tz) are the low-dimensional latent vectors corresponding to slice z at the time point t, which is formed by the reparameterization layer. We note that the generator θ is shared across all slices and time points; this approach facilitates the exploitation of the spatial redundancies between the slices and time points. We propose to jointly align and reconstruct the multislice MRI by jointly estimating the parameters θ, μ(tz) and Σ(tz) from the measured multislice data by minimizing the following cost function:
M S(θ,μ(tz),Σ(tz))=CM S(θ,μ(tz),Σ(tz))+λ1∥θ∥12+λ2∥∇tzμ(tz)∥2.  (33)

where

CM S=Σz=1NsliceΣt=1Ndata|tz[θ(c(tz))]−btz∥2+σ2L(q(tz)) is the lower bound for maximum likelihood as the first term in (32). The illustration of this scheme is given inFIG.24B. The parameters of the shared 3D generatorθare jointly learned in an unsupervised fashion from the measured k-t space data using the ADAM optimization algorithm. After the training process is complete, we will generate the image time series by feeding the generator with the latent variables of any specific slice. Following successful learning, we expect the volumes of the multislice reconstructions to have the same motion patterns characterized by the latent variables of that particular slice. We refer to this dynamic MRI reconstruction scheme as multislice variational SToRM, or V-SToRM.

D. Comparison with State-of-the-Art (SOTA) Methods

We compare the proposed V-SToRM approach with the following existing methods.

Analysis SToRM [9]: The analysis SToRM model uses a kernel low-rank formulation, which involves the esti-mation of the manifold Laplacian matrix from the k-space navigators using kernel low-rank regularization. This Laplacian is then used to solve for the images. We note that the analysis SToRM approach has been demonstrated to yield improved performance over state-of-the-art self-gated methods, as shown in our prior work [9,10]. We refer to this approach as A-SToRM.

Single-slice generative SToRM [51]: The single-slice generative SToRM approach uses a CNN generator to generate the single-slice image series from the highly undersampled k-t space data. This scheme does not rely on a variational formulation. It performs the independent recovery of each slice and hence fails to exploit the inter-slice redundancies. We refer to this approach as G-SToRM:SS.

Multislice generative SToRM: We extended the single-slice generative SToRM approach without the variational framework to the multislice setting. In particular, we use the CNN generator to produce the image volume; the generator parameters and the latent vectors for each slice are jointly learned. Finally, we feed the latent variables of a particular slice into the generator to obtain the aligned multislice reconstruction. We refer to this approach as G-SToRM:MS.

For the quantitative comparisons, in addition to the SSIM metric, we also use the signal-to-error ratio (SER) defined as

S⁢E⁢R=20·log10⁢xrefxref-xrecon.

Here, xrefand xreconrepresent the reference and the recon-structed images, respectively. The unit for SER is decibel (dB). In our free-breathing and ungated cardiac MRI setting, we usually do not have access to the ground truth. Therefore, in our work, we employ the analysis SToRM method using 25 seconds of data for the reconstruction as the simulated ground truth.

V. EXPERIMENTS AND RESULTS

A. Implementation Details

In this work, we use deep CNN to build the generator. The number of generator output channels is dependent on the specific datasets. For the experiments using the MNIST dataset, the channel is chosen as 1. By contrast, a two-channel output corresponding to the real and imaginary parts of the MR images is used for the rest of the experiments. In the MRI setting, we use a generator of 10 layers. The total number of trainable parameters is about 6 times the size of the image volume. For the convolutional layers in the generator, the activation function is chosen as leaky ReLU [99] except for the final layer, where tan h is used as the activation function. Random initialization is used to initialize the generator network.

The algorithm has three free parameters, σ2, λ1, and λ2. For each method, we optimize these parameters as well as the architecture of the generator on a single dataset such that the reconstructions closely match the 25-second A-SToRM reconstructions. Once the optimal parameters are determined, they are kept fixed for the remaining datasets. Our experiments showed that two latent vectors were sufficient for the good recovery of the single-slice datasets, which correspond to the cardiac and respiratory phases. In the multislice case, we required three to obtain good reconstructions. In this case, two of the three latent vectors captured cardiac and respiratory motion, respectively. The third latent vector seemed to capture a harmonic of the respiratory motion.

B. Single-Slice V-SToRM and Comparisons

In this section, we focus on single-slice V-SToRM; the reconstructions of a dataset and its latent vectors are shown inFIG.25. We trained the variational model using the data of one slice. The latent vectors we obtained are shown at the bottom ofFIG.25. Four different phases in the time series are shown in the figure, and their corresponding latent vectors are indicated in the plot of the latent vectors.

The comparisons between the single-slice V-SToRM and the state-of-the-art methods on a different dataset are shown inFIG.26. In these experiments, we compare the region of interest for A-SToRM, G-SToRM, and V-SToRM reconstructions using the 7.5 seconds of data. We use A-SToRM reconstructions from 25 seconds of data as the reference. FromFIG.26, we see that G-SToRM (7.5 s) and V-SToRM (7.5) are able to reduce errors and noise in the images when compared to A-SToRM (7.5 s). The proposed V-SToRM (7.5 s) is able to provide sharper edges than G-SToRM (7.5 s). These observations are further confirmed by the quantitative results shown at the bottom of the figure.

C. Joint Alignment and Recovery of Multislice Data

In this section, we show the results of the joint alignment and recovery of multislice data using the proposed multislice V-SToRM scheme. We also compare the alignment results obtained from the straightforward multislice extension of the G-SToRM scheme. The results are shown inFIG.5. More results are shown in the supplementary material.

The dataset used inFIG.27was acquired with eight slices that covered the whole heart. We trained the variational model based on the undersampled k-t space data and fed the latent vectors corresponding to the second slice to the generator, which produces the aligned multislice reconstructions. Shown in the figures are four time points based on the different phases identified by the latent variables. The rows inFIG.27Acorrespond to diastole in End-Inspiration, diastole in End-Expiration, systole in End-Inspiration, and systole in End-Expiration for each slice obtained using the proposed multislice V-SToRM scheme. FromFIG.27A, we see that the proposed multislice V-SToRM scheme is able to jointly reconstruct and align the multislice free-breathing and ungated cardiac MRI. We note that all the slices in each row have the same cardiac phase and respiratory phase.

InFIG.27B, we show the corresponding results for the direct extension of the multislice G-SToRM approach. In particular, we trained the model using the undersampled k-t space data and fed the latent vectors corresponding to the second slice into the generator to produce the aligned multislice reconstructions. FromFIG.27B, we see that the multislice G-SToRM approach has some ability to align the multislice reconstructions. However, we find that the image quality for some of the frames (e.g., slices 5-8) is poor. For example, the diastole phases for the G-SToRM:MS reconstructions are blurred and the cardiac boundaries are missing.

The reason for the poor reconstructions offered by multislice G-SToRM and the improved performance of V-SToRM can be easily appreciated from the distribution of the latent vectors shown inFIG.27CandFIG.27D, respectively. The use of the variational formulation in V-SToRM encouraged the latent variables of the slices to approximate a Gaussian distribution. We also reported the KL divergence value compared to(0, I) for each set of the latent vector in the figure. We note that the V-SToRM scheme offers low KL divergence values, indicating that the latent distribution of all the slices are roughly similar to a unit Gaussian. By contrast, the G-SToRM scheme cannot guarantee that the latent variables follow any distribution. We note from the top rows ofFIG.27Dthat the distribution of the latent variables of the second slice is very different from that of the other slices. When we feed the latent vectors of the second slice into the generator, the generator is only able to generate reasonable results for the second slice.

D. Comparison of Image Quality With State-of-the-Art Methods

We compare the image quality of the multislice V-SToRM reconstructions with the image quality of the reconstructions from the state-of-the-art methods, including single-slice methods, inFIGS.28A and28B. Note that the motion patterns of the slices recovered by the single-slice methods may be very different. For comparison, we manually matched the images of the slices of the single-slice and multislice methods by their cardiac and respiratory phases. The quantitative comparisons of the slices are shown at the bottom of each sub-figure. We also show more results using another dataset in the supplementary material.

The single-slice A-SToRM and G-SToRM:SS comparisons roughly match the observations inFIG.28AandFIG.28Band the results in [31]. The results show that the multislice V-SToRM approach is able to offer reconstructions that are less blurred and have fewer alias artifacts when compared to single-slice analysis methods (A-SToRM and G-SToRM:SS). The improved performance is also evidenced by the higher SER and SSIM values. We attribute the improved performance to the exploitation of the redundancies across slices, enabled by V-SToRM. We also note that the G-SToRM:MS method offers poor performance, evidenced by image blurring and missing details on the myocardium. The poor performance of G-SToRM:MS can be understood in terms of the differences in distribution of the latent vectors, shown inFIG.27A through27D.

VI. DISCUSSION AND CONCLUSION

In this work, we introduced an approach for the variational learning of a CNN manifold model from undersampled measurements. This work generalized the traditional VAE scheme to the undersampled setting. Unlike the traditional VAE scheme that uses an encoder to learn the conditional distribution from the images, we propose to learn the parameters of the distribution from the measurements using back-propagation. The application of the framework to multislice cardiac MRI data enabled the joint alignment and recovery from highly undersampled measurements. Unlike current single-slice methods that perform independent recovery of the slices, the proposed approach aligns the acquisitions and jointly recovers the images from the undersampled k-t space data. In addition to facilitating the exploitation of inter-slice redundancies, this approach also eliminates the need for post-processing schemes to match the phases of the slices.

Our results show that the joint alignment and recovery of the slices offer reduced blurring and reduction of artifacts compared to the direct generalization of G-SToRM to the multislice setting. In particular, the variational framework encourages the latent variables of different slices to have the same distribution. By contrast, the G-SToRM framework cannot guarantee the similarity of the probability distributions; the improper alignment translates to image blurring and other artifacts. Similarly, the use of the CNN generator offers implicit spatial regularization, resulting in improved recovery over A-SToRM.

A benefit with the proposed scheme is that it does not re-quire fully sampled data to train the CNN. The subject-specific CNN parameters and the latent vectors are learned directly from the undersampled data. We note that the acquisition of fully sampled data to train neural networks is not always possible, especially in the high-resolution and dynamic MRI settings considered in this work. In this context, direct learning from undersampled data is desirable. However, a challenge of the proposed scheme when compared to pretrained deep learning methods that offer super-fast inference is the higher computational complexity. We will explore training strategies, including transfer learning and meta-learning, to reduce the run time in the future.

VII. APPENDIX

In this appendix, we show that the likelihood term in (7) can be lower-bounded by (8). According to (7) and using the result of joint probability, we obtain

p⁡(bi)=p⁡(bi,ci)qi(ci)⁢qi(ci)p⁡(ci|bi)=p⁡(bi,ci)pi(ci)︸p⁡(bi|ci)⁢p⁡(ci)qi(ci)⁢q⁡(ci)pi(ci|bi).(1⁢3)

Taking the logarithm on both sides of (13), we have

log⁢p⁡(bi)=log⁢p⁡(bi|ci)-log⁢qi(ci)pi(ci)+log⁢qi(ci)pi(ci|bi).(14)

Next, we take the expectation with respect to ci˜qi(ci) of both sides of (14), and realizing thatci˜qi(ci) log p(bi)=log p(bi), we obtain

log⁢p⁡(bi)=𝔼ci~qi(ci)⁢log⁢qi(ci)p⁡(ci|bi)︸data⁢term-𝔼ci~qi(ci)⁢log⁢qi(ci)p⁡(ci)︸K⁢L[qi(ci)⁢p⁡(c1)]+𝔼ci~qi(ci)⁢log⁢qi(ci)p⁡(ci|bi)︸K⁢L[qi(ci)⁢p⁡(ci|bi)]>0.(15)

The last term is always greater than zero. The first term is the conditional density of the measurements bigiven the images xi=θ(ci). With the measurement model specified by (4), we obtain

𝔼ci~qi(ci)⁢log⁢p⁡(bi❘ci)=-12⁢σ2⁢𝔼ci~qi(ci)⁢𝒜i⁢𝒟⁡(ci)-bi2+c,

Where c is a constant independent of the parameters of interest. Ignoring the constant c and pluggingci˜qi(ci)log p(bi|ci) back into (15), we obtain the desired lower bound (8).

Although specific ensembles are shown and described, the present invention contemplates various options and alternatives. The methods described herein or aspects thereof may be incorporated into software in the form of instructions stored on a non-transitory computer or machine readable medium. Thus, it is contemplated that existing imaging systems may be reprogrammed to perform methods described herein. It is further to be understood that the methods described here may allow for portable systems to provide improved performance and thus imaging systems such as for performing MRI may be constructed specifically to utilize the methods shown and described herein.

Throughout this specification, plural instances may implement components, operations, or structures described as a single instance. Although individual operations of one or more methods are illustrated and described as separate operations, one or more of the individual operations may be performed concurrently, and nothing requires that the operations be performed in the order illustrated. Structures and functionality presented as separate components in example configurations may be implemented as a combined structure or component. Similarly, structures and functionality presented as a single component may be implemented as separate components. These and other variations, modifications, additions, and improvements fall within the scope of the subject matter herein.

Certain embodiments may be described herein as implementing mathematical methodologies including logic or a number of components, modules, or mechanisms. Modules may constitute either software modules (e.g., code embodied on a machine-readable medium or in a transmission signal) or hardware modules. A hardware module is tangible unit capable of performing certain operations and may be configured or arranged in a certain manner. In example embodiments, one or more computer systems (e.g., a standalone, client or server computer system) or one or more hardware modules of a computer system (e.g., a processor or a group of processors) may be configured by software (e.g., an application or application portion) as a hardware module that operates to perform certain operations as described herein.

In various embodiments, a hardware module may be implemented mechanically or electronically. For example, a hardware module may comprise dedicated circuitry or logic that is permanently configured (e.g., as a special-purpose processor, such as a field programmable gate array (FPGA) or an application-specific integrated circuit (ASIC)) to perform certain operations. A hardware module may also comprise programmable logic or circuitry (e.g., as encompassed within a general-purpose processor or other programmable processor) that is temporarily configured by software to perform certain operations. It will be appreciated that the decision to implement a hardware module mechanically, in dedicated and permanently configured circuitry, or in temporarily configured circuitry (e.g., configured by software) may be driven by cost and time considerations.

Accordingly, the term “hardware module” should be understood to encompass a tangible entity, be that an entity that is physically constructed, permanently configured (e.g., hardwired), or temporarily configured (e.g., programmed) to operate in a certain manner or to perform certain operations described herein. As used herein, “hardware-implemented module” refers to a hardware module. Considering embodiments in which hardware modules are temporarily configured (e.g., programmed), each of the hardware modules need not be configured or instantiated at any one instance in time. For example, where the hardware modules comprise a general-purpose processor configured using software, the general-purpose processor may be configured as respective different hardware modules at different times. Software may accordingly configure a processor, for example, to constitute a particular hardware module at one instance of time and to constitute a different hardware module at a different instance of time.

Hardware modules can provide information to, and receive information from, other hardware modules. Accordingly, the described hardware modules may be regarded as being communicatively coupled. Where multiple of such hardware modules exist contemporaneously, communications may be achieved through signal transmission (e.g., over appropriate circuits and buses) that connect the hardware modules. In embodiments in which multiple hardware modules are configured or instantiated at different times, communications between such hardware modules may be achieved, for example, through the storage and retrieval of information in memory structures to which the multiple hardware modules have access. For example, one hardware module may perform an operation and store the output of that operation in a memory device to which it is communicatively coupled. A further hardware module may then, at a later time, access the memory device to retrieve and process the stored output. Hardware modules may also initiate communications with input or output devices, and can operate on a resource (e.g., a collection of information).

The various operations of example methods described herein may be performed, at least partially, by one or more processors that are temporarily configured (e.g., by software) or permanently configured to perform the relevant operations. Whether temporarily or permanently configured, such processors may constitute processor-implemented modules that operate to perform one or more operations or functions. The modules referred to herein may, in some example embodiments, comprise processor-implemented modules. Where the term “processor” is used, it is to be understood that it encompasses one or more processors whether located together or remote from one other.

Similarly, the methods described herein may be at least partially processor-implemented. For example, at least some of the operations of a method may be performed by one or processors or processor-implemented hardware modules. The performance of certain of the operations may be distributed among the one or more processors, not only residing within a single machine, but deployed across a number of machines. In some example embodiments, the processor or processors may be located in a single location (e.g., within a hospital, clinic, or medical office environment), while in other embodiments the processors may be distributed across a number of locations.

The one or more processors may also operate to support performance of the relevant operations in a “cloud computing” environment or as a “software as a service” (SaaS). For example, at least some of the operations may be performed by a group of computers (as examples of machines including processors), these operations being accessible via a network (e.g., the Internet) and via one or more appropriate interfaces (e.g., application program interfaces (APIs).)

The performance of certain of the operations may be distributed among the one or more processors, not only residing within a single machine, but deployed across a number of machines. In some example embodiments, the one or more processors or processor-implemented modules may be located in a single geographic location (e.g., within a hospital, an office environment, or a server farm). In other example embodiments, the one or more processors or processor-implemented modules may be distributed across a number of geographic locations.

Some portions of this specification are presented in terms of algorithms or symbolic representations of operations on data stored as bits or binary digital signals within a machine memory (e.g., a computer memory). These algorithms or symbolic representations are examples of techniques used by those of ordinary skill in the data processing arts to convey the substance of their work to others skilled in the art. As used herein, an “algorithm” is a self-consistent sequence of operations or similar processing leading to a desired result. In this context, algorithms and operations involve physical manipulation of physical quantities. Typically, but not necessarily, such quantities may take the form of electrical, magnetic, or optical signals capable of being stored, accessed, transferred, combined, compared, or otherwise manipulated by a machine. It is convenient at times, principally for reasons of common usage, to refer to such signals using words such as “data,” “content,” “bits,” “values,” “elements,” “symbols,” “characters,” “terms,” “numbers,” “numerals,” or the like. These words, however, are merely convenient labels and are to be associated with appropriate physical quantities.

Unless specifically stated otherwise, discussions herein using words such as “processing,” “computing,” “calculating,” “determining,” “presenting,” “displaying,” or the like may refer to actions or processes of a machine (e.g., a computer) that manipulates or transforms data represented as physical (e.g., electronic, magnetic, or optical) quantities within one or more memories (e.g., volatile memory, non-volatile memory, or a combination thereof), registers, or other machine components that receive, store, transmit, or display information.

As used herein any reference to “one embodiment” or “an embodiment” means that a particular element, feature, structure, or characteristic described in connection with the embodiment is included in at least one embodiment. The appearances of the phrase “in one embodiment” in various places in the specification are not necessarily all referring to the same embodiment.

As used herein, the terms “comprises,” “comprising,” “includes,” “including,” “has,” “having” or any other variation thereof, are intended to cover a non-exclusive inclusion. For example, a process, method, article, or apparatus that comprises a list of elements is not necessarily limited to only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Further, unless expressly stated to the contrary, “or” refers to an inclusive or and not to an exclusive or. For example, a condition A or B is satisfied by any one of the following: A is true (or present) and B is false (or not present), A is false (or not present) and B is true (or present), and both A and B are true (or present).

In addition, use of the “a” or “an” are employed to describe elements and components of the embodiments herein. This is done merely for convenience and to give a general sense of the disclosure. This description should be read to include one or at least one and the singular also includes the plural unless it is obvious that it is meant otherwise.

The terms “first,” “second,” “third,” “fourth,” and the like in the description and in the claims, if any, are used for distinguishing between similar elements and not necessarily for describing a particular sequential or chronological order. It is to be understood that the terms so used are interchangeable under appropriate circumstances such that the embodiments described herein are, for example, capable of operation in sequences other than those illustrated or otherwise described herein. Similarly, if a method is described herein as comprising a series of steps, the order of such steps as presented herein is not necessarily the only order in which such steps may be performed, and certain of the stated steps may possibly be omitted and/or certain other steps not described herein may possibly be added to the method.

As used herein, a plurality of items, structural elements, compositional elements, and/or materials may be presented in a common list for convenience. However, these lists should be construed as though each member of the list is individually identified as a separate and unique member. Thus, no individual member of such list should be construed as a de facto equivalent of any other member of the same list solely based on their presentation in a common group without indications to the contrary.

Reference throughout this specification to “an example” means that a particular feature, structure, or characteristic described in connection with the example is included in at least one embodiment. Thus, appearances of the phrases “in an example” in various places throughout this specification are not necessarily all referring to the same embodiment or example.

The techniques presented and claimed herein are referenced and applied to material objects and concrete examples of a practical nature that demonstrably improve the present technical field and, as such, are not abstract, intangible or purely theoretical. Further, if any claims appended to the end of this specification contain one or more elements designated as “means for [perform]ing [a function or “step for [perform]ing [a function] . . . ”, it is intended that such elements are to be interpreted under 35 U.S.C. 112 § (f). However, for any claims containing elements designated in any other manner, it is intended that such elements are not to be interpreted under 35 U.S.C. § 112(f).

The invention is not to be limited to the particular embodiments described herein. In particular, the invention contemplates numerous variations in the specific methodology used with respect to the deep manifold learning algorithms. The foregoing description has been presented for purposes of illustration and description. It is not intended to be an exhaustive list or limit any of the invention to the precise forms disclosed. It is contemplated that other alternatives or exemplary aspects are considered included in the invention. The description is merely examples of embodiments, processes, or methods of the invention. It is understood that any other modifications, substitutions, and/or additions can be made, which are within the intended spirit and scope of the invention.

REFERENCES

The following references are hereby incorporated by reference in their entireties.[1] J. Tsao, P. Boesiger, and K. P. Pruessmann, “k-t blast and k-t sense: dynamic mri with high frame rate exploiting spatiotemporal correlations,” Magnetic Resonance in Medicine, vol. 50, no. 5, pp. 1031-1042, 2003.[2] S. G. Lingala, Y. Hu, E. DiBella, and M. Jacob, “Accelerated dynamic mri exploiting sparsity and low-rank structure: kt slr,” IEEE Transactions on Medical Imaging, vol. 30, no. 5, pp. 1042-1054, 2011.[3] L. Feng et al., “Golden-angle radial sparse parallel mri: combination of compressed sensing, parallel imaging, and golden-angle radial sampling for fast and flexible dynamic volumetric mri,” Magnetic Resonance in Medicine, vol. 72, no. 3, pp. 707-717, 2014.[4] L. Feng, L. Axel, H. Chandarana, K. T. Block, D. K. Sodickson, and R. Otazo, “Xd-grasp: golden-angle radial mri with reconstruction of extra motion-state dimensions using compressed sensing,” Magnetic Resonance in Medicine, vol. 75, no. 2, pp. 775-788, 2016.[5] A. G. Christodoulou et al., “Magnetic resonance multitasking for motion-resolved quantitative cardiovascular imaging,” Nature Biomedical Engineering, vol. 2, no. 4, pp. 215-226, 2018.[6] C. Prieto et al., “Highly efficient respiratory motion compensated free-breathing coronary mra using golden-step cartesian acquisition,” Journal of Magnetic Resonance Imaging, vol. 41, no. 3, pp. 738-746, 2015.[7] A. Bustin, N. Fuin, R. M. Botnar, and C. Prieto, “From compressed-sensing to artificial intelligence-based cardiac mri reconstruction,” Frontiers in Cardiovascular Medicine, vol. 7, pp. 17, 2020.[8] S. Poddar and M. Jacob, “Dynamic mri using smoothness regularization on manifolds (storm),” IEEE Transactions on Medical Imaging, vol. 35, no. 4, pp. 1106-1115, 2015.[9] A. H. Ahmed, R. Zhou, Y. Yang, P. Nagpal, M. Salerno, and M. Jacob, “Free-breathing and ungated dynamic mri using navigator-less spiral storm,” IEEE Transactions on Medical Imaging, 2020.[10] S. Poddar, Y. Q. Mohsin, D. Ansah, B. Thattaliyath, R. Ashwath, and M. Jacob, “Manifold recovery using kernel low-rank regularization: Application to dynamic imaging,” IEEE Transactions on Computational Imaging, vol. 5, no. 3, pp. 478-491, 2019.[11] U. Nakarmi, Y. Wang, J. Lyu, D. Liang, and L. Ying, “A kernel-based low-rank (klr) model for low-dimensional manifold recovery in highly accelerated dynamic mri,” IEEE Transactions on Medical Imaging, vol. 36, no. 11, pp. 2297-2307, 2017.[12] U. Nakarmi, K. Slavakis, and L. Ying, “Mls: Joint manifold-learning and sparsity-aware framework for highly accelerated dynamic magnetic resonance imaging,” in 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018). IEEE, 2018, pp. 1213-1216.[13] G. Wang, “A perspective on deep imaging,” IEEE Access, vol. 4, pp. 8914-8924, 2016.[14] G. Wang, M. Kalra, and C. G. Orton, “Machine learning will transform radiology significantly within the next 5 years,” Medical Physics, vol. 44, no. 6, pp. 2041-2044, 2017.[15] G. Dardikman-Yoffe and Y. C. Eldar, “Learned sparcom: Unfolded deep super-resolution microscopy,” arXiv preprint arXiv:2004.09270, 2020.[16] J. C. Ye, Y. Han, and E. Cha, “Deep convolutional framelets: A general deep learning framework for inverse problems,” SIAM Journal on Imaging Sciences, vol. 11, no. 2, pp. 991-1048, 2018.[17] K. H. Jin, M. T. McCann, E. Froustey, and M. Unser, “Deep convolutional neural network for inverse problems in imaging,” IEEE Transactions on Image Processing, vol. 26, no. 9, pp. 4509-4522, 2017.[18] V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” arXiv preprint arXiv:1912.10557, 2019.[19] A. Pramanik, H. K. Aggarwal, and M. Jacob, “Deep generalization of structured low-rank algorithms (deep-slr),” IEEE Transactions on Medical Imaging, 2020.[20] S. Wang et al., “Accelerating magnetic resonance imaging via deep learning,” in 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI). IEEE, 2016, pp. 514-517.[21] S. Wang et al., “Deepcomplexmri: Exploiting deep residual network for fast parallel mr imaging with complex convolution,” Magnetic Resonance Imaging, vol. 68, pp. 136-147, 2020.[22] S. UH Dar, M. O ̈zbey, A. B. C, atli, and T. C, ukur, “A transfer-learning approach for accelerated mri using deep neural networks,” Magnetic resonance in medicine, vol. 84, no. 2, pp. 663-685,2020.[23] Y. Han et al., “Deep learning with domain adaptation for accelerated projection-reconstruction mr,” Magnetic resonance in medicine, vol. 80, no. 3, pp. 1189-1205, 2018.[24] T. Sanchez et al., “Scalable learning-based sampling optimization for compressive dynamic mri,” in ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 8584-8588.[25] S. Wang et al., “Dimension: dynamic mr imaging with both k-space and spatial prior knowledge obtained via multi-supervised network training,” NMR in Biomedicine, p. e4131, 2019.[26] S. Wang et al., “Lantern: Learn analysis transform network for dynamic magnetic resonance imaging,” Inverse Problems & Imaging, p. 1, 2020.[27] S. UH Dar et al., “Prior-guided image reconstruction for accelerated multi-contrast mri via generative adversarial networks,” IEEE Journal of Selected Topics in Signal Processing, vol. 14, no. 6, pp. 1072-1087, 2020.[28] S. UH Dar et al., “Image synthesis in multi-contrast mri with conditional generative adversarial networks,” IEEE transactions on medical imaging, vol. 38, no. 10, pp. 2375-2388, 2019.[29] M. Yurt, S. U H Dar, A. Erdem, E. Erdem, and T. C, ukur, “mustgan: Multi-stream generative adversarial networks for mr image synthesis,” arXiv preprint arXiv:1909.11504, 2019.[30] Q. Zou, S. Poddar, and M. Jacob, “Sampling of planar curves: Theory and fast algorithms,” IEEE Transactions on Signal Processing, vol. 67, no. 24, pp. 6455-6467, 2019.[31] Q. Zou and M. Jacob, “Recovery of surfaces and functions in high dimensions: sampling theory and links to neural networks,” SIAM Journal on Imaging Sciences, in press.[32] D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Deep image prior,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2018, pp. 9446-9454.[33] K. H. Jin, H. Gupta, J. Yerly, M. Stuber, and M. Unser, “Time-dependent deep image prior for dynamic mri,” arXiv preprint arXiv:1910.01684, 2019.[34] A. Bora, E. Price, and A. G. Dimakis, “Ambientgan: Generative models from lossy measurements,” in International Conference on Learning Representations, 2018.[35] K. Kazuhiro et al., “Generative adversarial networks for the creation of realistic artificial brain magnetic resonance images,” Tomography, vol. 4, no. 4, pp. 159, 2018.[36] I. Goodfellow et al., “Generative adversarial nets,” in Advances in neural information processing systems, 2014, pp. 2672-2680.[37] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein gan,” arXiv preprint arXiv:1701.07875, 2017.[38] D. Varga, A. Csisz{acute over ( )}arik, and Z. Zombori, “Gradient regularization improves accuracy of discriminative models,” arXiv preprint arXiv:1712.09936, 2017.[39] R. Zhou et al., “Free-breathing cine imaging with motion-corrected reconstruction at 3 t using spiral acquisition with respiratory correction and cardiac self-gating (spares),” Magnetic resonance in medicine, vol. 82, no. 2, pp. 706-720, 2019.[40] D. O. Walsh, A. F. Gmitro, and M. W. Marcellin, “Adaptive reconstruction of phased array mr imagery,” Magnetic Resonance in Medicine, vol. 43, no. 5, pp. 682-690, 2000.[41] A. Mittal, A. K. Moorthy, and A. C. Bovik, “No-reference image quality assessment in the spatial domain,” IEEE Transactions on Image Processing, vol. 21, no. 12, pp. 4695-4708, 2012.[42] A.-K. Bracher, C. Hofmann, A. Bornstedt, S. Boujraf, E. Hell, J. Ulrici, A. Spahr, B. Haller, and V. Rasche, “Feasibility of ultra-short echo time (ute) magnetic resonance imaging for identification of carious lesions,” Magnetic resonance in medicine, vol. 66, no. 2, pp. 538-545, 2011.[43] N. S. Higano, A. D. Hahn, J. A. Tkach, X. Cao, L. L. Walkup, R. P. Thomen, S. L. Merhar, P. S. Kingma, S. B. Fain, and J. C. Woods, “Retrospective respiratory self-gating and removal of bulk motion in pulmonary ute mri of neonates and adults,” Magnetic resonance in medicine, vol. 77, no. 3, pp. 1284-1295, 2017.[44] A. D. Hahn, N. S. Higano, L. L. Walkup, R. P. Thomen, X. Cao, S. L. Merhar, J. A. Tkach, J. C. Woods, and S. B. Fain, “Pulmonary mri of neonates in the intensive care unit using 3d ultrashort echo time and a small footprint mri system,” Journal of Magnetic Resonance Imaging, vol. 45, no. 2, pp. 463-471, 2017.[45] N. S. Higano, D. R. Spielberg, R. J. Fleck, A. H. Schapiro, L. L. Walkup, A. D. Hahn, J. A. Tkach, P. S. Kingma, S. L. Merhar, S. B. Fain, et al., “Neonatal pulmonary magnetic resonance imaging of bronchopulmonary dysplasia predicts short-term clinical outcomes,” American journal of respiratory and critical care medicine, vol. 198, no. 10, pp. 1302-1311, 2018.[46] F. Gibiino, L. Sacolick, A. Menini, L. Landini, and F. Wiesinger, “Free-breathing, zero-te mr lung imaging,” Magnetic Resonance Materials in Physics, Biology and Medicine, vol. 28, no. 3, pp. 207-215, 2015.[47] S. Kumar, R. Rai, A. Stemmer, S. Josan, L. Holloway, S. Vinod, D. Moses, and G. Liney, “Feasibility of free breathing lung mri for radiotherapy using non-cartesian k-space acquisition schemes,” The British journal of radiology, vol. 90, no. 1080, p. 20170037, 2017.[48] W. Jiang, F. Ong, K. M. Johnson, S. K. Nagle, T. A. Hope, M. Lustig, and P. E. Larson, “Motion robust high resolution 3d free-breathing pulmonary mri using dynamic 3d image self-navigator,” Magnetic resonance in medicine, vol. 79, no. 6, pp. 2954-2967, 2018.[49] D. P. Capaldi, R. L. Eddy, S. Svenningsen, F. Guo, J. S. Baxter, A. J. McLeod, P. Nair, D. G. McCormack, G. Parraga, and C. R. R. Network, “Free-breathing pulmonary mr imaging to quantify regional ventilation,” Radiology, vol. 287, no. 2, pp. 693-704, 2018.[50] Z. Berecova, V. Neuschl, P. Boruta, J. Masura, and E. Ghersin, “A complex pulmonary vein varix-diagnosis with ecg gated mdct, mri and invasive pulmonary angiography,” Journal of radiology case reports, vol. 6, no. 12, p. 9, 2012.[51] Q. Zou, A. H. Ahmed, P. Nagpal, S. Kruger, and M. Jacob, “Dynamic imaging using a deep generative storm (gen-storm) model,” IEEE Transactions on Medical Imaging, 2021.[52] J. Yoo, K. H. Jin, H. Gupta, J. Yerly, M. Stuber, and M. Unser, “Time-dependent deep image prior for dynamic mri,” IEEE Transactions on Medical Imaging, 2021.[53] N. R. Huttinga, T. Bruijnen, C. A. van den Berg, and A. Sbrizzi, “Nonrigid 3d motion estimation at high temporal resolution from prospectively undersampled k-space data using low-rank mr-motus,” Magnetic resonance in medicine, vol. 85, no. 4, pp. 2309-2326, 2021.[54] X. Zhu, M. Chan, M. Lustig, K. M. Johnson, and P. E. Larson, “Iterative motion-compensation reconstruction ultra-short to (imoco ute) for high-resolution free-breathing pulmonary mri,” Magnetic resonance in medicine, vol. 83, no. 4, pp. 1208-1221, 2020.[55] T. Li, M. Zhang, W. Qi, E. Asma, and J. Qi, “Motion correction of respiratory-gated pet images using deep learning based image registration framework,” Physics in Medicine & Biology, vol. 65, no. 15, p. 155003, 2020.[56] L. Feng, J. Delacoste, D. Smith, J. Weissbrot, E. Flagg, W. H. Moore, F. Girvin, R. Raad, P. Bhattacharji, D. Stoffel, et al., “Simultaneous evaluation of lung anatomy and ventilation using 4d respiratory-motion-resolved ultrashort echo time sparse mri,” Journal of Magnetic Resonance Imaging, vol. 49, no. 2, pp. 411-422, 2019.[57] Q. Zou, A. H. Ahmed, P. Nagpal, S. Kruger, and M. Jacob, “Deep generative storm model for dynamic imaging,” in 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI), pp. 114-117, IEEE, 2021.[58] F. Knoll, K. Bredies, T. Pock, and R. Stollberger, “Second order total generalized variation (tgv) for mri,” Magnetic resonance in medicine, vol. 65, no. 2, pp. 480-491, 2011.[59] J. A. Fessler and B. P. Sutton, “Nonuniform fast fourier transforms using min-max interpolation,” IEEE transactions on signal processing, vol. 51, no. 2, pp. 560-574, 2003.[60] F. Ong and M. Lustig, “Sigpy: a python package for high performance iterative reconstruction,” Proceedings of the International Society of Magnetic Resonance in Medicine, Montr{acute over ( )}eal, QC, vol. 4819, 2019.[61] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in 3rd International Conference on Learning Representations (ICLR), 2015.[62] A. F. Agarap, “Deep learning using rectified linear units (relu),” arXiv preprint arXiv:1803.08375, 2018.[63] K. M. Johnson, S. B. Fain, M. L. Schiebler, and S. Nagle, “Optimized 3d ultrashort echo time pulmonary mri,” Magnetic resonance in medicine, vol. 70, no. 5, pp. 1241-1250, 2013.[64] H. Pedersen, S. Kozerke, S. Ringgaard, K. Nehrke, and W. Y. Kim, “k-t pca: temporally constrained k-t blast reconstruction using principal component analysis,” Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, vol. 62, no. 3, pp. 706-716, 2009.[65] K. Breuer, C. B. Meyer, F. A. Breuer, A. Richter, F. Exner, A. M. Weng, S. Str ̈ohle, B. Polat, P. M. Jakob, O. A. Sauer, et al., “Stable and efficient retrospective 4d-mri using non-uniformly distributed quasi-random numbers,” Physics in Medicine & Biology, vol. 63, no. 7, p. 075002, 2018.[66] O. Dietrich, J. G. Raya, S. B. Reeder, M. F. Reiser, and S. O. Schoenberg, “Measurement of signal-to-noise ratios in mr images: influence of multichannel coils, parallel imaging, and reconstruction filters,” Journal of Magnetic Resonance Imaging: An Official Journal of the International Society for Magnetic Resonance in Medicine, vol. 26, no. 2, pp. 375-385, 2007.[67] M. Gutberlet, T. F. Kaireit, A. Voskrebenzev, F. Lasch, J. Freise, T. Welte, F. Wacker, J. M. Hohlfeld, and J. Vogel-Claussen, “Free-breathing dynamic 19f gas mr imaging for mapping of regional lung ventilation in patients with copd,” Radiology, vol. 286, no. 3, pp. 1040-1051, 2018.[68] N. Kobayashi, “Magnetic resonance imaging with gradient sound respiration guide,” Plos one, vol. 16, no. 7, p. e0254758, 2021.[69] J. L. Prince and J. M. Links, Medical imaging signals and systems. Pearson Prentice Hall Upper Saddle River, 2006.[70] S. Napel, M. P. Marks, G. D. Rubin, M. D. Dake, C. H. McDonnell, S. M. Song, D. R. Enzmann, and R. Jeffrey Jr, “Ct angiography with spiral ct and maximum intensity projection,” Radiology, vol. 185, no. 2, pp. 607-610, 1992.[71] J. F. Gruden, S. Ouanounou, S. Tigges, S. D. Norris, and T. S. Klausner, “Incremental benefit of maximum-intensity-projection images on observer detection of small pulmonary nodules revealed by multidetector ct,” American Journal of Roentgenology, vol. 179, no. 1, pp. 149-157, 2002.[72] G. E. Hinton and R. R. Salakhutdinov, “Reducing the dimensionality of data with neural networks,” Science, vol. 313, no. 5786, pp. 504-507, 2006.[73] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” arXiv preprint arXiv:1312.6114, 2013.[74] I. Higgins et al., “Early visual concept learning with unsupervised deep learning,” arXiv preprint arXiv:1606.05579, 2016.[75] I. Higgins et al., “beta-vae: Learning basic visual concepts with a constrained variational framework,” 2016.[76] J. Chung et al., “A recurrent latent variable model for sequential data,” Advances in Neural Information Processing Systems, vol. 28, pp. 2980-2988, 2015.[77] Y. Li, Z. Wang, R. Sun, and F. Lam, “Separation of metabolites and macromolecules for short-te 1 h-mrsi using learned component-specific representations,” IEEE Transactions on Medical Imaging, vol. 40, no. 4, pp. 1157-1167, 2021.[78] M. Mani, V. A. Magnotta, and M. Jacob, “qmodel: A plug-and-play model-based reconstruction for highly accelerated multi-shot diffusion mri using learned priors,” Magnetic Resonance in Medicine, vol. 86, no. 2, pp. 835-851, 2021.[79] T. Kido et al., “Compressed sensing real-time cine cardiovascular magnetic resonance: accurate assessment of left ventricular function in a single-breath-hold,” Journal of Cardiovascular Magnetic Resonance, vol. 18, no. 1, pp. 1-11, 2016.[80] L. Axel and R. Otazo, “Accelerated mri for the assessment of cardiac function,” The British Journal of Radiology, vol. 89, no. 1063, pp. 20150655, 2016.[81] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse mri: The application of compressed sensing for rapid mr imaging,” Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182-1195, 2007.[82] T. Kustner et al., “Isotropic 3D Cartesian single breath-hold CINE MRI with multi-bin patch-based low-rank reconstruction,” Magnetic Resonance in Medicine, vol. 84, no. 4, 2020.[83] C. Qin et al., “Convolutional Recurrent Neural Networks for Dynamic MR Image Reconstruction,” IEEE Transactions on Medical Imaging, vol. 38, no. 1, pp. 280-290, 2019.[84] A. Bustin, N. Fuin, R. M. Botnar, and C. Prieto, “From compressed-sensing to artificial intelligence-based cardiac mri reconstruction,” Fron tiers in Cardiovascular Medicine, vol. 7, no. 17, 2020.[85] T. Kustner et al., “Cinenet: deep learning-based 3d cardiac cine mri reconstruction with multi-coil complex-valued 4d spatio-temporal convolutions,” Scientific Reports, vol. 10, no. 1, pp. 1-13, 2020.[86] C. M. Sandino, P. Lai, S. S. Vasanawala, and J. Y. Cheng, “Accelerating cardiac cine mri using a deep learning-based espirit reconstruction,” Magnetic Resonance in Medicine, vol. 85, no. 1, pp. 152-167, 2021.[87] T. Wang et al., “Ica-unet: Ica inspired statistical unet for real-time 3d cardiac cine mri segmentation,” in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2020, pp. 447-457.[88] Z. Deng et al., “Four-dimensional mri using three-dimensional radial sampling with respiratory self-gating to characterize temporal phase-resolved respiratory motion in the abdomen,” Magnetic Resonance in Medicine, vol. 75, no. 4, pp. 1574-1585, 2016.[89] S. Rosenzweig, N. Scholand, H. C. M. Holme, and M. Uecker, “Cardiac and respiratory self-gating in radial mri using an adapted singular spectrum analysis (ssa-fary),” IEEE Transactions on Medical Imaging, vol. 39, no. 10, pp. 3029-3041, 2020.[90] M. Usman, D. Atkinson, C. Kolbitsch, T. Schaeffter, and C. Prieto, “Manifold learning based ECG-free free-breathing cardiac CINE MRI,” Journal of Magnetic Resonance Imaging, vol. 41, no. 6, 2015.[91] X. Chen et al., “High-Resolution Self-Gated Dynamic Abdominal MRI Using Manifold Alignment,” IEEE Transactions on Medical Imaging, vol. 36, no. 4, pp. 960-971, 2017.[92] G. N. Shetty et al., “Bi-linear modeling of data manifolds for dynamic-mri recovery,” IEEE Transactions on Medical Imaging, vol. 39, no. 3, pp. 688-702, 2019.[93] M. Gori, M. Maggini, and L. Sarti, “Exact and approximate graph matching using random walks,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 27, no. 7, pp. 1100-1111, 2005.[94] C. F. Baumgartner et al., “High-resolution dynamic mr imaging of the thorax for respiratory motion correction of pet using groupwise manifold alignment,” Medical Image Analysis, vol. 18, no. 7, pp. 939-952, 2014.[95] C. F. Baumgartner et al., “Self-aligning manifolds for matching disparate medical image datasets,” in International Conference on Information Processing in Medical Imaging. Springer, 2015, pp. 363-374.[96] X. Chen et al., “Dynamic volume reconstruction from multi-slice abdominal mri using manifold alignment,” in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2016, pp. 493-501.[97] Y. LeCun, “The mnist database of handwritten digits,” http://yann.lecun.com/exdb/mnist/, 1998.[98] M. Uecker et al., “Espirit-an eigenvalue approach to autocalibrating parallel mri: where sense meets grappa,” Magnetic Resonance in Medicine, vol. 71, no. 3, pp. 990-1001, 2014.[99] B. Xu, N. Wang, T. Chen, and M. Li, “Empirical evaluation of rectified activations in convolutional network,” arXiv preprint arXiv:1505.00853, 2015.