Time-Dependent Deep Image Prior for Dynamic MRI

I Introduction

The aim of dynamic magnetic resonance imaging (MRI) is to capture the dynamics associated with moving organs, which requires a fast imaging process. A typical approach is to accelerate data acquisition by a partial sampling of the k-space. The resulting partial loss of data must then be compensated to maintain the image quality. Several methods have addressed this by exploiting spatial or temporal redundancy, including parallel MRI , k-t acceleration methods , compressed sensing (CS) MRI , low‐rank methods , and many others. In the specific case of cardiac applications, the current state-of-the-art methods further improve the reconstruction by exploiting the fact that the heart motion is approximately cyclic. They typically use electrocardiograms or self-gating techniques . However, all of these methods are limited by constraints over the signal-to-noise ratios (SNR), restrictions in the coil design, hand-picked priors, multiple processing steps, or inefficient algorithms in their deployment of the standard convex-optimization techniques.

More recently, inspired by the development of deep-learning techniques in various imaging modalities , supervised-learning approaches have been applied to the fast and accurate reconstruction of partially sampled MRI . These methods, however, heavily depend on a training dataset, especially on ground-truth data (i.e., fully sampled measurements), which are typically unavailable for dynamic MRI. Unlike the direct deep-learning approaches, the model-based deep-learning framework of formulates the image recovery as an optimization scheme. By unrolling an iterative algorithm, it minimizes a cost function that combines data consistency and a deep-learned prior. Because the learned prior incorporates patient-specific noise patterns into the algorithm, this approach successfully recovers images with fast reconstruction and acceptable quality. However, it still requires ground-truth data to train the denoising network.

In this paper, we propose an unsupervised learning framework in which a generative network is optimized to reconstruct a sequence of golden-angle radial lines in k-space, also called spokes. Inspired by deep image priors (DIP) , we use a convolutional neural networks (CNN) architecture as an implicit structural prior that constrains the search space of the optimization problem. In addition, to learn the temporal dependencies of the dynamic measurements, we impose a one-dimensional manifold parameterized by time. Aided by this explicit cue, the network then learns to encode the temporal variations of the sequential images into the spatial closeness of the samples on the imposed manifold. This simple temporal coupling already enables our model to outperform the other CS algorithms without bells and whistles—note that our approach is purely unsupervised and optimized in an end-to-end manner. We further improve the reconstruction by introducing a mapping network (MapNet) that brings more flexibility to our latent space . MapNet consists of a few fully connected layers with nonlinear activations; it learns to map the fixed manifold into a more expressive latent space. This allows the subsequent generative network to adapt its input to a given dataset, thereby improving image quality (Figure 1).

In short, our generative model takes the latent variables from MapNet and reconstructs dynamic images by exploiting its powerful structural prior. With the extensive analyses in Section IV and experimental results in Section V, we show that both the manifold design and MapNet are essential to achieve good reconstructions. To the best of our knowledge, this is the first unsupervised deep-learning-based method that can reconstruct the full temporal frames of dynamic MRI sequences with high spatial resolution.

I-B Related Work

Unsupervised Learning. Starting from the seminal work of DIP , there have been several studies that applied unsupervised learning to medical imaging, such as MRI and positron emission tomography , albeit both cases address the reconstruction of static images.

The closest work to ours is the one that used DIP for video compression , which also considers a sequence of latent inputs. However, unlike our goal (the reconstruction of an image sequence), theirs is to find compact codes for the representation of video frames. To find such codes, they optimize both the network weights and latent variables. Without any constraint on the latent space, however, the latent codes may diverge to an arbitrary space. To prevent this, they imposed either low-rank or similarity constraints on the latent sequence. However, this optimization not only requires additional effort to tune hyper-parameters but also entails a singular-value decomposition at each iteration, which severely increases the computational burden. By contrast, we solve this by simply inputting an explicit manifold and letting a mapping network adapt the manifold to the given data. This makes the training much easier and yields the one-dimensional manifold in latent space that is adapted to the given data. In addition, their forward model is an identity operator, while ours is an MR measurement operator with severe under-sampling.

II Methods

We first briefly recapitulate the content of deep image prior (Section II-A) as well as the physics of dynamic MRI (Section II-B). Then, we describe our method based on DIP with a mapping network and on the learning of the underlying latent manifolds (Section II-C).

II-B Dynamic MRI

The measurement process that relates the image at time $t$ and the $k$-space measurements with angle $\vartheta=\vartheta(t)$ is linear and formally described by the relation

where ${\bf{H}}(\vartheta)$ is the $M_{0}\times N$ system matrix that represents the combined effect of taking the 2D Fourier transform of ${\bf{x}}$ and resampling along a radial line with direction $\vartheta$. The type of measurement provided by (2) is referred to as an angular spoke. In practice, we acquire a series of $K$ spokes taken at regularly spaced time point $t_{k}=t_{0}+k\Delta t$, $k=0,\dots,(K-1)$ with step size $\Delta t$. The spoke orientations follow the golden-angle strategy

II-B2 Spoke-Sharing

To further describe this pooling process, we introduce the augmented measurement vector {\mathbf{y}}_{{k}}=\big{(}\mathbf{y}(t_{m})\big{)}_{m=k-(n_{\rm s}-1)/2}^{k+(n_{\rm s}-1)/2} of size $M=n_{\rm s}\times M_{0}$. Correspondingly, we define the column-wise concatenated system matrix ${\mathbf{H}}_{{k}}=({\mathbf{H}}(\vartheta_{k}))_{m=k-(n_{\rm s}-1)/2}^{k+(n_{\rm s}-1)/2}$, whose time dependence is indicated by the index $k$. This results in the forward imaging model

II-C Proposed Framework

To address the dynamic MRI reconstruction problem, we first modify the original DIP so that it takes a sequence of input and output pairs (Figure 1). More specifically, we optimize an untrained neural network $f_{\boldsymbol{\theta}}$ to map a sequence of inputs $\{{\bf{z}}_{k}\}_{k=0}^{K-1}$ to the spoke-shared measurements $\{{\bf{y}}_{k}\}_{k=0}^{K-1}$, thereby reconstructing the sequence of images $\{{\bf{x}}(t_{k})\}_{k=0}^{K-1}$ by searching for

leading to ${\bf{x}}^{*}(t_{k})=f_{\boldsymbol{\theta}^{*}}({\bf{z}}_{k})$. Note that the optimization is done in the measurement domain. This enforces the image sequence to be consistent with the measurements, while the modified DIP scheme regularizes the reconstructed images.

Manifold Design. To fully exploit the characteristics of dynamic MRI, the underlying model must be able to effectively encode the temporal variations of the measurements while preserving the structure of the individual frames. To this end, we propose to design a manifold $\mathcal{Z}$, thereby effectively injecting a specific prior into the network. For example, an ordered sequence $\{{\bf{z}}_{k}\}$ from a straight-line manifold will guide the network to associate spatial closeness of input variables with temporal closeness of images. This encourages the network to reconstruct an image sequence with temporally similar attributes. For a quasi-periodic signal such as the cardiac motion, we can encode the expected behavior by letting the manifold take the structure of a three-dimensional helix.

Mapping Network (MapNet). Although a careful choice of temporally meaningful manifolds typically results in an excellent performance, the fact that the design is hand-crafted may also sometimes limit the performance of the network . To add flexibility to our model and to exploit the rich representation power of the network, we introduce a mapping network (MapNet). In our design, MapNet $g_{\phi}$ involves a few fully connected layers with nonlinearities. It learns to map a fixed manifold into the more expressive latent space $\mathcal{W}=g_{\phi}(\mathcal{Z})$. More specifically, our model $f_{\boldsymbol{\theta}}$ now has a hierarchical architecture that consists of the MapNet $g_{\phi}$ followed by CNN $h_{\psi}$ so that $f_{\boldsymbol{\theta}}=h_{\psi}\circ g_{\phi}$ and $\boldsymbol{\theta}=\{\phi,\psi\}$ (Figure 1 (B)). This leads us to replace (5) by

The role of $g_{\phi}$ is to appropriately warp the input manifold to facilitate $h_{\psi}$ in its reconstruction of the true dynamics. Overall, the insertion of $g_{\phi}$ provides better flexibility to our model and lets us efficiently exploit the representation power of neural networks, resulting in a good reconstruction.

Final Algorithm. Our optimization scheme is given in Algorithm 1. We minimize the loss function (6) using standard gradient-descent methods for $n_{\rm iter}$ iterations. At each iteration, instead of (6), a batch loss $L_{B}(\boldsymbol{\theta})$ is updated where a batch $\{k_{0},\ldots,k_{B-1}\}$ of size $B$ is randomly sampled from the index set $\{0,\ldots,K-1\}$. The corresponding input variables $\{{\bf{z}}_{k_{b}}\}_{b=0}^{B-1}$ are fed to the network and its parameters are updated using the gradient with respect to $\boldsymbol{\theta}$.

III Experiments

In this section, we describe the datasets, baseline methods, cardiac-cycle estimation, evaluation setups, and implementation details.

All experimental datasets are breath-hold MR images. We assume a twofold upsampling of measurements for every dataset. Therefore, the size of the reconstructed fields of view is half that of the first dimension of the measurements.

A cardiac cine dataset was acquired using a 3T whole-body MRI scanner (Siemens; Tim Trio) equipped with a 32-element cardiac coil array. The acquisition sequence was bSSFP and prospective cardiac gating was used. The imaging parameters were as follows: FOV=$(300\times 300)\,\mbox{mm}^{2}$, acquisition matrix size=$(128\times 128)$, TE/TR=$1.37/2.7\,\mbox{ms}$, receiver bandwidth=$1184\,\mbox{Hz/pixel}$, and flip angle=$40^{\circ}$. The number of frames was $23$ and the temporal resolution was $43.2\,\mbox{ms}$. The resulting fully sampled Cartesian trajectories are used as ground-truth. To retrospectively simulate the radial sampling, we implemented the forward model using the golden-angle strategy with NuFFThttps://github.com/marchdf/python-nufft. Sinograms are obtained as shown in Figure 1. The number of spokes per frame is $n_{s}=13$. For a single-cycle simulation, the dimension of sinograms is $(K\times n_{s}\times M_{\mathbf{\omega}}\times C)=(23\times 13\times 256\times 32)$. For a multicycle simulation, we acquire $p=13$ cycles, which results in $K=13\cdot 23=299$ frames.

III-A2 Fetal Cardiac Dataset

Fetal cardiac MRI data were acquired on a 1.5 T clinical MR scanner (MAGNETOM Aera, Siemens AG, Healthcare Sector, Erlangen, Germany) with an 18-channel body array coil and a 32-channel spine coil for signal reception. We used an untriggered continuous 2D bSSFP sequence that was modified to acquire radial readouts with a golden-angle trajectory . The acquisition parameters were: FOV = (260 $\times$ 260) $\text{mm}^{2}$, acquisition matrix size = (256 $\times$ 256) pixels, slice thickness = 4.0 mm, TE/TR = 1.99/4.1 ms, RF excitation angle = 70∘, radial readouts = 1400, acquisition time = 6.7 s, and bandwidth = 1028 Hz/pixel.

III-B Baseline Methods

Back Projection (BP) is a zero-filled discrete Fourier transform, which is the most basic baseline one can think of.

GRASP is a golden-angle radial sparse parallel MRI algorithm, which extends the idea of k-t SPARSE-SENSE to volumetric golden-angle radial acquisitions. Here, the spoke-sharing strategy is not applied.

Reordering Method (RD) is a three-step algorithm. RD first reconstructs real-time images of limited image quality and uses these images to reorder or self-gate the measurements, which in turn are used for the final reconstruction with k-t SPARSE-SENSE . In the retrospective experiment, where we know the phase indices, we use the exact order of frames for self-gating.

III-C Estimation of cardiac Cycles

For the processing of the fetal cardiac dataset, RD and our algorithm both require a rough estimate of the number of cardiac cycles seen over the whole duration of a sequence of data acquisition. It can be typically obtained from k-space. Simple techniques to estimate the cardiac cycles from radial data have been previously reported by . Radial acquisition schemes sample the center of k-space at every readout, which supports the extraction of physiological motion signals. The central k-space coefficient of a radial readout (i.e., the echo peak) corresponds to the complex sum of the transverse magnetization across the entire image volume. In the presence of moving structures such as a beating heart, changes in the overall transverse magnetization due to motion will induce a modulation of the consecutive echo peaks. (Trajectory imperfections and eddy currents can also modulate echo peaks, but their frequency responses differ from the physiological motion frequencies and, thus, can be filtered out.) The resulting signal can then be used to estimate the number of cardiac cycles and to inform the manifold network. For our fetal cardiac dataset, we find that the time-course has approximately 13 periods so that we finally set $p=13$.

III-D Evaluation Metric

We use the regressed SNR as a quantitative metric. With the oracle $\mathbf{x}$ and the reconstructed image ${\mathbf{x}}^{*}$, RSNR is given by

where a higher RSNR corresponds to a better reconstruction.

III-E Implementation Details

We use an Intel i7-7820X (3.60GHz) CPU and an NVIDIA Titan X (Pascal) GPU. Pytorch 1.0.0 on Python 3.6 is used to implement our generative modelWe shall provide a link to the repository upon paper acceptance.. The network is optimized until $n_{\text{iter}}=10,000$ with $B=1$ using Adam optimizer of default setting and the learning rate of $10^{-3}$.

III-F Architectures

The mapping network $g_{\phi}$ is two consecutive fully connected layers of 512 hidden dimension with ReLU in between. It outputs $L=64$-dimensional latent vector, which is reshaped to $\left(8\times 8\right)$ for the following generative network $h_{\psi}$ (Table I). The generative network consists of convolutional layers, batch normalization layers, ReLU, and nearest-neighbor interpolations. We apply zero-padding before convolution to let the size of the output mirror that of the input. At the last layer, ReLU is not used. The output has two channels because MRI images take complex values.

IV Design of the latent space

In this section, we analyze the individual components of our model and compare the performance with baselines. We first demonstrate the simplest setup that reconstructs a single heart cycle. We then move on to a more complicated dataset that has multiple heart cycles.

Although simple, this configuration already outperforms the other baseline methods and successfully reconstructs the dynamics for a single cycle dataset (Table II).

IV-B Manifolds for Multiple Heart Cycles

In practice, the measurements generally span several heart cycles. To better exploit the fact that the cardiac movement has a quasi-periodic behavior, it is of interest to explore more sophisticated manifolds.

Effect of the Manifolds. In Figure 2, we show the reconstructed (y-t) images of the cross section that is denoted by a white line in GT (y-x) imageFor display purposes, we show only one cycle of our cross section.. When we use a straight-line manifold, the network fails to capture the heart movement and outputs the same static image over all frames. This is natural since most of the pixels are static and the dynamic parts are localized in a small area. Thus, the network easily finds a local minimum that corresponds to an image that remains constant over all frames. However, as soon as we switch to “periodic-like” manifold designs, the network starts to reconstruct the movement (Table III). For example, when we use a line with 13 segments as an input, the performance is better than the RD that uses the same information. Using circles with 13 repetitions as input manifold, we improve even further. However, the helix input manifold gives the best performance among the others without MapNet because the heartbeat is a quasi-periodic signal.

Effect of the Mapping Network. In addition to the choice of its manifold, our method has another design choice: its mapping network. By introducing MapNet, the network can adapt its input manifold to a given dataset, which allows us to further improve the reconstruction (Table III). This can be clearly seen in the t-SNE visualization of the mapped latent space (Figure 6), which we discuss in Section VI.

In summary, our analysis shows that a careful design of the manifold and the use of a mapping network are both necessary to achieve the best performance. Based on these, from now on, we use ‘Helix+MapNet’ as our default setup.

V Results

We first show results on the retrospective dataset, where the desired behaviors of the reconstruction methods are well-defined. We then illustrate on the fetal cardiac dataset that the observations extend well to a real scenario.

The benefits of our method are evident in both the (y-t) view (Figure 2) and (y-x) view (Figure 3) of each frame. In Figure 2, both GRASP and RD reconstruct the movement of the heart. RD shows better performance than GRASP, which was expected because it takes advantage of the period information that is estimated while reordering the frames. However, as can be seen in the residuals, GRASP and RD show significant errors in the reconstruction of the dynamics. In the (y-x) view of Figure 3, GRASP leads to blurring artifacts, while the residual image reveals errors around the wall of the heart in both GRASP and RD reconstructions. By contrast, our method gives better results with fewer artifacts.

V-B Fetal Cardiac Dataset

Having demonstrated the superior behavior of our method on the retrospective dataset, we now assess our model on real data. In the absence of ground-truth, we shall take the static image that is generated from all spokes as pseudo-gold standard—note that it is of high quality only in the regions that are not moving.

Like in the retrospective experiments, both GRASP and RD are able to reconstruct multiple cardiac phases. RD gives better reconstructions, especially in the dynamic region (Figure 4 (A)). However, RD shows a spurious artifact at the edge area (Figure 4 (B)) and fails to find the detailed structures of the static background (Figure 4 (C)). By contrast, our method produces better-resolved features in both dynamic and static areas (particularly for the hyperintense dot-like structures in Figure 4 (A)), while it does not suffer from artifacts at the edges and recovers the low-intensity background areas as well (Figure 4 (C)).

In Figure 5, it is apparent that BP completely fails in capturing the fetal cardiac beats. The GRASP reconstruction is less noisy but still far from satisfactory. RD fares better; unfortunately, its reordering process can lead it to superpose in the same frame spokes that belong to different phases of the cardiac cycle. By contrast, our method reconstructs each frame with data from just a few neighboring spokes, thus avoiding the mingling of different cycles. The reconstructed systolic phase captures the true motion of the heart better. The cross section from our method is similar to that of RD but the motion is smoother in our case, which is the expected behavior of a beating heart.

We provide in Figure 5 (Bottom row) our whole reconstructed sequence of cardiac cycles. The quasi-periodicity of the cardiac motion is clearly visible along the temporal axis, while motion variations can still be discerned from cycle to cycle. Note that this is a unique benefit of our method that the other algorithms cannot provide.

VI Discussion

To assess the extent of structural change as a function of time, we used t-stochastic neighborhood embeddings (t-SNE) which capture the underlying manifold by projecting the high-dimensional entities onto a three-dimensional space. In Figure 6 (A), we show the t-SNE result of the original manifold when the variables are generated according to (10) with $L=64$ and $p=13$. Unsurprisingly, this recovers a helix with 13 cycles. In Figure 6 (B), we show the embedding of the 64-dimensional mapped variables $\{g_{\phi^{*}}({\bf{z}}_{k})\}_{k=0}^{K-1}$, where the ${\bf{z}}_{k}$ are generated according to (10) with $L=3$ and $p=13$. Again, we recover a helical geometry with 13 periods, although the height of the helix is now shortened—the first and the last cycles become closer. This shows that MapNet successfully recovers some similarity between the different cycles which, in turn, translates into better reconstructions. It warps the given manifold in adaptive fashion, while retaining the prior information that we inject via the manifold geometry. As shown in Section IV, this design (Helix+MapNet) outperforms the ‘Helix’ with $L=64$ and $p=13$ which is fed directly to the vanilla CNN. In Figure 6 (C), we display the projected manifold of the $(256\times 256)$-dimensional reconstructed images. It shows a helical structure with 13 local folds, each of which corresponds to a single cycle of the cardiac motion. This also shows that the quasi-periodic characteristic of the data is well represented by the network.

VI-B Benefits of Our Approach

Continuous Dynamic Reconstruction. One major benefit of our approach is that it lets us reconstruct temporally continuous dynamic images. We showed that the network $f_{\boldsymbol{\theta}^{*}}$ successfully captures the underlying nonlinear dynamics of the image manifold, and the input variable ${\bf{z}}_{k}$ lets us reconstruct the image at the corresponding time stamp (Figure 6). Because our method represents images as a learned parametric function $f_{\boldsymbol{\theta}^{*}}$, we can recover nontrivial intra-frame images by navigating between two consecutive input variables, which would not be possible with other standard interpolation methods such as temporal bilinear interpolation.

Memory Savings. In the methods based on compressed sensing (CS), the gradient updates of the iterative optimization process necessitate memory that is large enough to hold the target reconstruction volume. For example, the reconstruction of $5,\!000$ frames with spatial size $\left(256\times 256\right)$ would need one to handle data of size $\left(256\times 256\times 5,000\right)$, which demands for over a gigabyte of memory. Our approach, by contrast, requires much less memory. It optimizes the neural network using batches, which requires the simultaneous handling of only those frames that correspond to the batch size. In short, the fact that our proposed approach handles few 2D images whereas CS handles a 2D+t extended sequence leads to substantial savings, particularly for golden-angle dynamic MRI with many frames. In our approach, we only store a 2D generative model; for example, its memory demands for the spatial size $\left(256\times 256\right)$ are about half-a-dozen megabytes. This cost is negligible compared to that of the CS approach.

Efficient Reconstruction. Our model visits each frame about seven times during training (10,000 iterations / 1400 frames $\approx 7$ ), while GRASP sees all frames during the entire iterations (24 outer iterations). Regarding the execution time, the major bottleneck of our method is the slow forward model. It depends on the NuFFT package which, in its current implementation, does not benefit from a GPU and is a major cause for slowdown. Indeed, NuFFT takes 47 % of the entire running time of our algorithm per each iteration; the average processing time for 100 repetitions is 6.55 s for back and forth NuFFTs, and 3.08 s for the remaining parts. With a more efficient implementation, our algorithm could be substantially accelerated.

Because our model is fully automated, it leads to a simpler optimization task with fewer hyperparameters than the conventional methods. For instance, k-t SENSE requires three interdependent hyperparameters whose optimal values are found only after some substantial grid-search effort, while the two hyperparameters of our approach are easier to interpret since they trivially consist of just an initial learning rate, along with a number of iterations.

VII Conclusion

In this paper, we proposed an unsupervised deep-learning-based algorithm for dynamic MRI reconstruction that provides high spatial resolution with access to the sub-frame—or even continuous—temporal control of dynamic images. By designing a one-dimensional manifold, combined with the mapping network, our generative network model fully exploits the representation power of the network as well as its structural priors. Our study showed that the proposed method successfully reconstructs dynamic MRI in an end-to-end manner and outperforms the state-of-the-art CS approaches by 3.8 dB. To the best of our knowledge, this is the first unsupervised-learning approach in accelerated dynamic MRI.

Acknowledgements

The authors thank Prof. Jong Chul Ye at KAIST for providing the bSSFP cardiac MRI k-space dataset (retrospective dataset).

References