Multiresolution reconstruction of real-time MRI with motion compensated compressed sensing: Application to 2D free-breathing cardiac MRI

Real-time MRI is a novel noninvasive imaging technique that allows the visualization of physiological processes with both good spatial and temporal resolutions. However, the reconstruction of images from highly undersampled data, needed to perform real-time imaging, remains challenging. Recently, the combination of Compressed Sensing theory with motion compensation techniques has shown to achieve better results than previous methods. In this work we describe a real-time MRI algorithm based on the acquisition of the k-space data following a Golden Radial trajectory, Compressed Sensing reconstruction and a groupwise temporal registration algorithm for the estimation and compensation of the motion in the image, all this embedded within a temporal multiresolution scheme. We have applied the proposed method to the reconstruction of free-breathing acquisition of short axis views of the heart, achieving a temporal resolution of 25ms, corresponding to an acceleration factor of 28 with respect to fully sampled Cartesian acquisitions.


INTRODUCTION
The recent development of signal analysis based acceleration techniques in dynamic MRI has significantly pushed the achievable acceleration factors to unforeseen levels. Recently, the combination of Compressed Sensing (CS) theory [1] with motion estimation (ME) and motion compensation (MC) techniques has been successfully applied to CS reconstruction of dynamic MRI based on the assumption that sparser representations of the signal can be obtained if some knowledge of the motion in the image is incorporated during the reconstruction [2,3].
This work was partially supported by the Spanish Ministerio de Ciencia e Innovación and the European Regional Development Fund (ERDF-FEDER) under Research Grant TEC2014-57428-R, the Spanish Ministerio de Ciencia e Innovación under Research Grant TEC2013-44194-P, the Spanish Junta de Castilla y León under Grant VA136U13 and by the Universidad de Valladolid-Banco de Santander grant program.
In real-time (RT) MRI, the high acceleration factors needed to achieve the desired spatio-temporal resolution make the reconstruction problem highly ill-conditioned, becoming very sensitive to noise in the data and model inconsistencies. Moreover, ME/MC techniques lack the true dynamic image to estimate the motion from. A common approach is to perform the ME step on an initial reconstruction without MC which, however, may be corrupted by the undersampling artifacts that cannot be corrected in the first step. These artifacts, in turn, would hinder the estimation of the true motion in the dynamic image. In a previous work [3], a groupwise (GW) temporal registration procedure, which showed to be more robust to artifacts than its pairwise counterparts, was proposed for ME/MC reconstruction of breath-hold cine acquisitions. However, the acceleration factors achieved could not enable its direct application to free-breathing RT imaging, due to the limitations of the Cartesian acquisition employed.
To overcome the previous limitations, in this work we propose a reconstruction algorithm with three key components. First, we adopt the golden-angle radial (GR) trajectory [4] for continuous data acquisition in which the angular step between consecutive spokes is given by the Golden Ratio. For a given undersampling factor, GR provides a more efficient coverage of the k-space than the Cartesian scheme used in [3]. Second, the GR allows to retrospectively set up an arbitrary window length (i.e. temporal resolution), which enables us to introduce a multiresolution approach in which finer temporal resolution reconstructions are obtained iteratively from coarser ones. Third, for ME/MC we resort to the GW registration method employed in [3] to iteratively obtain MEs from coarser to finer temporal resolutions, avoiding the need of estimating the motion directly from highly undersampled data..
The proposed method is applied to the reconstruction of free-breathing, 2D short axis views of the heart of three healthy volunteers. Results are compared with those obtained with the non linear inversion (NLINV) method [5], in which both the images and coil sensitivity maps are jointly estimated. In NLINV, temporal regularization is introduced by enforcing the similarity between consecutive frames in a sequential way.

THEORY
The reconstruction of accelerated MRI data can be generally formulated as solving a linear inversion problem y = Em, where y stands for the acquired multi-coil k-t space data, m represents the sequence of images to be reconstructed and the encoding operator E comprises the multiplication by the coils sensitivity maps and the frame-by-frame undersampled spatial Fourier transform. Under the MC-CS framework [6] both the image and the motion information are jointly reconstructed. In our proposal, we split the problem in two subproblems (reconstruction and ME/MC) and alternate between both following a multiresolution scheme.

Motion Compensated Reconstruction
The reconstruction problem, assuming motion is known beforehand, is formulated as: with Ψ a spatio-temporal transform and λ a trade off parameter between data fidelity and sparsity promotion. T Θ is a warping operator governed by the set of parameters Θ designed to compensate the motion in the image.

Groupwise temporal registration for ME/MC
During the ME/MC step a GW dissimilarity metric that provides a unique global measure of the differences between the images comprised within the sequence is minimized. In the current implementation, a simple sum of squared differences metric (SSD) [7] is used: where T Θ m stands for the temporal average of the registered sequence. C(T Θ ) is a regularization term related to the elastic energy of a bending plate [8] which penalizes non-smooth deformations. As it can be derived from this expression, there is no need for a fixed reference frame. The SSD metric relies on the assumption that the intensity of each pixel is preserved along time. Other metrics can be incorporated for different imaging modalities where such a condition is not met [9]. The deformation model used in the GW registration algorithm is based on a non-rigid 2D+t free form deformations (FFD) with cubic B-splines, which has been widely used in practice and has shown to be flexible enough to describe the motion and deformation of anatomical structures [8]. Eq. (2) is solved with a nonlinear conjugate gradient method.

Multiresolution reconstruction of Real-Time MRI
As discussed in the introduction, a major difficulty is to obtain both good initial reconstructions and accurate ME from highly undersampled data. To this end, we propose a multiresolution scheme (see Fig. 1) in which coarser reconstructions serve as initialization of the finer ones. An initial, low temporal resolution sequence is obtained by a regular CS reconstruction equivalent to solving Eq. 1 with T Θ set to the identity. Since no motion information is available at this point, Ψ is set to a spatial wavelet transform. For the next levels temporal total variation is used. At each level, the method alternates until convergence between reconstruction (Eq. 1) and ME (Eq. 2). In order to provide a good initialization to the next resolution level, frames could be interpolated as the weighted average of the surrounding frames. However, high inter-frame motion is present such straightforward interpolation can lead to corrupted images. To improve the quality of the interpola- tion, we make use of the motion information in the estimated deformation. Following the approach in [7], for each new frame to be interpolated at instant t j we obtain the transfor- which aligns a (moving) frame at instant t i to the frame at instant t j (fixed image). The originally described B-splines registration method can be extended to guarantee the existence of T −1 Θ (x, t j ) [10]. In our current implementation, the spatio-temporal regularization term and proper spacing on the grid of control points favor the invertibility and the inverse transform is therefore found numerically.

APPLICATION TO 2D FREE-BREATHING CARDIAC MRI
Three healthy volunteers were scanned with a 32-element cardiac coil and a golden radial trajectory on a 1.5T Philips scanner. Other relevant scan parameters include: b-SSFP sequence, TR/TE/α = 2.9ms/1.44ms/60 • , FOV = 320x320mm 2 , spatial resolution = 2x2mm 2 , slice thickness = 8mm. A single slice was continuously imaged for about 12 seconds. Images were reconstructed following the scheme in Fig.  1 with three resolution levels and 15 and 9 radial spokes per frame used in the last level (135, 45 and 15 spokes per frame in each level for the first case and 81, 27 and 9 in the second). Corresponding final temporal resolutions are 44 and 26 ms. Parameter λ was set to 0.1 attending to visual inspection of the results from one dataset and maintained for the other two. NLINV reconstructions were performed on the same data using the code available at http://www.user.gwdg.de/∼muecker/realtime.html after griding the data to a Cartesian lattice with a Kaisser-Bessel window of size L = 6 and shape parameter β = 13.8551. The number of Newton iterations was set to 3 and the scale factor for regularization to 0.9 and, as proposed, a median temporal filter of length 5 is applied after reconstruction to remove remaining artifacts (see Ref. [5] for details).
Reconstructions were run offline on a server with two Intel Xeon E5-2695 v3 CPUs @ 2.30 GHz and 64 GB of RAM using MATLAB (R2015a, The MathWorks, Natick, MA). Fig. 2 shows reconstructed images from one dataset at systole and diastole and the temporal evolution of a horizontal line across the left and right ventricles. The results with the proposed method (GW-MR) show sharper edges, and finer details can be observed in the interior of the left ventricle both at 15 and 9 spokes per frame. Differences in the reconstructions become more evident in their temporal evolution. NLINV reconstructions show poorer contrast between blood and myocardium in the left ventricle as well as severe remaining striking artifacts, specially at 9 spokes per frame, in spite of the temporal median filter required by NLINV method. In Fig. 3, NLINV results before temporal filtering are included for the other two volunteers at 9 spokes per frame. Remaining artifacts severely degrade the temporal evolution of the dynamic image and are not completely removed by the temporal filter. On the contrary, artifacts are better removed by GW-MR without additional post-processing steps while preserving spatial detail. Results for the three volunteers are shown jointly in the video provided as supplementary material.

RESULTS AND DISCUSSION
Average reconstruction times are plotted in Fig. 4. GW-MR approximately doubles the computational demand of NLINV with both algorithms running in MATLAB without further optimizations. 14-17 seconds per frame are required by GW-MR, what hinders its practical application. However, the algorithm is highly parallelizable and GPU implementa-tions will be investigated as future work.
One limitation of our study is the lack of a ground truth in RT imaging to obtain a quantitative measurement of the quality of the reconstructions. In future work, the usage of a numerical phantom will provide such a quality metric. Functional parameters such as ventricular volumes and ejection fraction computed both from RT and standard segmented acquisitions will be compared as well.

CONCLUSIONS
We have presented an algorithm based on continuous GR data acquisition and a robust GW registration method integrated in a temporal multiresolution scheme for ME/MC CS retrospective reconstruction of 2D RT MRI data. We have shown the application of this method to free-breathing, RT imaging of the heart employing only 9 spokes per frame, obtaining a final temporal resolution of 26 ms without view sharing or post-processing of the obtained images.