Abstract

Transformed domain sparsity of Magnetic Resonance Imaging (MRI) has recently been used to reduce the acquisition time in conjunction with compressed sensing (CS) theory. Respiratory motion during MR scan results in strong blurring and ghosting artifacts in recovered MR images. To improve the quality of the recovered images, motion needs to be estimated and corrected. In this article, a two-step approach is proposed for the recovery of cardiac MR images in the presence of free breathing motion. In the first step, compressively sampled MR images are recovered by solving an optimization problem using gradient descent algorithm. The -norm based regularizer, used in optimization problem, is approximated by a hyperbolic tangent function. In the second step, a block matching algorithm, known as Adaptive Rood Pattern Search (ARPS), is exploited to estimate and correct respiratory motion among the recovered images. The framework is tested for free breathing simulated and in vivo 2D cardiac cine MRI data. Simulation results show improved structural similarity index (SSIM), peak signal-to-noise ratio (PSNR), and mean square error (MSE) with different acceleration factors for the proposed method. Experimental results also provide a comparison between k-t FOCUSS with MEMC and the proposed method.

1. Introduction

Compressed sensing (CS) has been successfully implemented to reduce the scan time of MR images [13]. Application of CS to transformed domain MR images includes brain [1], cardiac [4], and pediatric MR imaging [5]. According to CS approach, a sparse or compressible image can be recovered, by solving norm mixed optimization problem, from randomly undersampled data using a nonlinear recovery technique [2, 6]. Since the -norm is not differentiable everywhere, an approximation to the -norm is used in the reconstruction algorithm. Different approaches [1, 7] have been used to approximate -norm using some differentiable functions. Reference [7] uses hyperbolic tangent based function as an approximation of -norm to solve the CS recovery problem for static MR images.

Different types of motion during the data acquisition process cause artifacts like ghosting and blurring in the recovered cardiac MR images. In the presence of respiratory motion, high quality MR images can be produced by combining k-space profiles of the same cardiac phases in ECG-gated MR acquisition [8]. Addition of the same cardiac phases at different respiratory motions or states produces inconsistencies in k-space which results in motion artifacts in the combined reconstructed images. Sparsity, a necessary condition for CS, is also violated by this k-space profile combination [9]. Hence, to avoid periodic breath holds during the acquisition process and to take advantage of sparse signal recovery from undersampled data using CS methods, motion artifacts may need to be corrected. The combined approach of CS and motion correction has been implemented in [811]. Otazo et al. proposed 1D translational respiratory motion correction. Usman et al. introduced a reconstruction scheme for dynamic cardiac MRI by incorporating general motion framework directly into CS reconstruction. Their method uses data binning and intensity based nonrigid registration algorithm for estimating respiratory motions. Reference [11] proposed a CS based motion correction in the free breathing environment with multiple constraints. This method uses Demon based registration to estimate the motion between reference and other respiratory states.

Interframe motion estimation and compensation for time varying features of images has been used in video compression standards [12, 13]. These standards are based on different block matching algorithms [14] for motion estimation and compensation. Similar to video compression, dynamic MR images can be predicted by exploiting temporal redundancies between the images. Asif et al. proposed an algorithm, MASTeR [15], for the breath held condition, based on interframe motion to recover different cardiac MR images. MASTeR uses motion adaptive transform that models temporal sparsity using interframe motion estimation. k-t FOCUSS [16] also uses interframe motion estimation and compensation with a fixed reference frame during the image recovery process for breath held cardiac cine MRI.

In this article, a novel framework is presented for the recovery of highly undersampled free breathing cardiac MR images. Similar cardiac phases at different respiratory states are grouped like the frames of a video sequence. The Adaptive Rood Pattern Search (ARPS) technique, based on the interframe motion, is used to estimate and correct respiratory states among the grouped images. A two-step approach is adopted for the reconstruction of dynamic MR images. In the first step, free breathing cardiac phases without motion estimation are recovered from undersampled k-space data. Next, interframe motion between the reconstructed cardiac phases is calculated using ARPS to improve the image estimates iteratively. An approximation of the -norm penalty is used in the gradient descent algorithm to recover dynamic MR images. The adjustable parameters of the -norm approximation provide an extra benefit, as it can be adjusted to reflect the changing statistics of dynamic MR images. After the application of the proposed method, the combinations of similar cardiac phases at different respiratory states are clear and accurate as compared to the combined cardiac phase without motion estimation and correction.

The rest of the paper is organized as follows. Section 2 discusses the preliminaries for interframe motion estimation in dynamic MRI and CS. Section 3 describes the methodology of the proposed algorithm. Section 4 gives the details of algorithm, Section 5 presents the simulation parameters and results followed by Section 6 that discusses the merits and demerits of the scheme, and Section 7 concludes the work.

2. Free Breathing Imaging Model and CS

Free breathing downsampled k-spaced data corrupted by motion states for cardiac phase is mathematically given as where is a two-dimensional complex MR image vector of length representing a cardiac phase at respiratory state , is a Fourier operator that transforms an image to k-space, is a random variable-density undersampling mask, different for all respiratory states and is a sensing matrix, and is a combined k-space measurement vector of length for nth cardiac phase acquired for all respiratory positions. A single cardiac phase n at respiratory state d in a specific heart cycle can be given asThe reduction or acceleration factor for MR images is given by . By increasing , the system in (1) becomes highly underdetermined. Compressed sensing solves such underdetermined system of equations effectively to recover MR images. The recovery of a sparse signal using CS can be achieved by solving the following convex optimization problem:where is the Lagrangian that provides a balance between sparsity and data consistency. The -norm keeps the solution consistent with the data and the -norm given by encourages sparsity in solution [2].

3. Methods

3.1. Smooth -Norm Approximation

CS algorithms recover sparse signals or images by solving the -norm regularized optimization problem such as that given in (3). In this article, we use gradient descent algorithm for solving the optimization problem with wavelet based penalty term. Nondifferentiability of the -norm at origin excludes the usage of mostly optimization approaches for the solution. Reference [1] uses an approximation of with the complex conjugate and a positive smoothing factor. Reference [7] proposes a new smooth function for approximating the -norm, used in this paper, which is given below

This function better approximates the absolute value and provides extra flexibility of adjusting the slope at the origin with the proper selection of and makes it more suitable for dynamic images. Mostly MR images are sparse in transformed domain so the modified version of cost function given in (3) for transformed MR images iswhere is a wavelet operator that transforms the image to a sparse domain.

In this article, we propose an iterative algorithm that uses the following approximation for the -norm penalty:where . The update equation for the algorithm, derived using the steepest descent method for a sparse vector , iswhere is positive valued step size and is the gradient operator that differentiates the cost function at th iteration. During each iteration, shrinkage given in (8) is applied in the wavelet domain after (7) to reconstruct the MR images.where is a thresholding parameter. By incorporating the approximation , the cost function can be written asThe gradient of the cost function is easy to compute:with

3.2. Respiratory Motion Based Dynamical System

Two main problems with free breathing cardiac MRI are as follows.

(1) Blurring artifacts are generated by the combination of k-space samples for the same cardiac phases at different respiratory states.

(2) The combination of k-space data in free breathing decreases the sparsity level.

In this article, we use interframe motion to estimate the respiratory states between the same cardiac phases. Video standards MPEG and H.264 [12, 13] have successfully exploited interframe motion for compression. In the dynamic MRI images, pixels are not significantly displaced in the neighboring frames. Pixel locations can be predicted using interframe motion estimation. Temporal redundancy among the frames is advantageous for the prediction of pixel locations. Let and be images having th cardiac phases at respiratory states d and , respectively. The pixel values of at location are closest to the pixel values at in . The displacement of all pixels in from to in is represented by motion vectors According to [2], cardiac phase at dth respiratory state can be generated from the cardiac phase at th respiratory state by the following equation:where is a backward transformation that uses information about the physical changes between two datasets of the same cardiac phases. The motion compensated residual is computed by taking the difference between predicted and compensated image. Using the transformation , a motion dependent linear system can be written by combining (1) and (12) as follows:To recover the cardiac phases , we solve (13a) and (13b) by exploiting sparse structure in and . The process of complete high resolution image generation is shown in Figure 1. The data is acquired in segmented fashion because MRI is a slow imaging modality. During the data scanning process, a limited k-space sample is recorded at each heart phase in all cardiac cycles. To simulate this condition, each cardiac phase at different respiratory states is multiplied with different sampling matrix .

3.3. Initial CS Reconstruction

A two-step approach is adopted to recover motion free cardiac phases. In the 1st step, images with motion effects are reconstructed from undersampled k-space data independently. For the recovery of dynamic cardiac images, the iterative algorithm, mentioned as Algorithm 1, optimizes the cost function given in (9) with the approximation given in (4). Wavelet based soft thresholding is used for the recovery of the Nth cardiac phase for each respiratory state. Daubechies-4 (db4) wavelet is used to exploit the transformed domain sparsity.

INPUTS
: k-space data
: Fourier Operator
: Sparsifying transform operator
λ = 0.005, η = 0.9, Maxiter = 50
OUTPUT
: Motion corrected final image
INITIALIZATION
,
INITIAL CS BASED RECOVERY
Step  1. For
Step  2. Update
Step  3. Shrinkage: (Using (8))
Step  4. End (i) Return
Repeat Steps1 to 4 for recovery of
MOTION ESTIMATION AND COMPENSATION
Step  5. For : Maxiter do
If
Step  6. Find motion compensated image for and
Else
Step  7. Refine motion between and
End (Else)
Step  8. Update
Step  9. Shrinkage: (Using (8))
End (j) Return refined image
Step  10.
3.4. Interframe Motion Estimation and Correction (MEMC)

In the 2nd step, interframe motion estimation and correction is performed from a pair of CS recovered images and divided into two substeps.

(a) Motion Estimation. Exploit initially CS reconstructed images to estimate or refine interframe motion and the motion transformation as follows. Cardiac phase in the 1st R-R ECG interval is taken as a P frame (frame to be predicted) and cardiac phases in subsequent R-R intervals are taken as an I frame (reference frame). The P and I frame are borrowed terminology from video compression. To estimate the motion between the 2nd cardiac phases of the 1st and 4th R-R interval, for example, we take the 2nd cardiac phase of the 1st interval as a P frame and the 2nd cardiac phase at different respiratory state in the 4th R-R interval as an I frame.

(b) Motion Correction. After finding motion vectors using ARPS block matching algorithm, we generate the corrected image from I frame and with the help of motion vectors. For the refinement of motion corrected image, solve the following optimization problem written for (13a) and (13b):where is given in (13b) and its -norm approximation isAt the final step, we generate the image of 2nd cardiac phase by combining P frame and motion corrected image to get an image with high temporal and spatial resolution.

4. Proposed Algorithm

Steps involved in the reconstruction of MR motion corrected images are given in Algorithm 1.

5. Experimental Setup and Results

The proposed method was tested on simulated data generated by the MRXCAT framework [17] and on fully sampled, free breathing, cine MRI data. The recovered images for CS-free breathing motion corrected were compared with CS-free breathing images and with CS-breath held images. All CS images were recovered in MATLAB (R2012a, MathWorks Inc., Natick, MA) using the proposed hyperbolic tangent based surrogate function to solve the nondifferentiability problem of the -norm penalty. In the gradient descent algorithm, step size in an update equation was chosen empirically. Parameter values used in the algorithm are β = 0.005, , and . The same values for , , and were used for initial CS reconstruction and for the final interframe motion estimation and compensation.

The structural similarity index (SSIM) [18], peak signal-to-noise ratio (PSNR), and mean square error (MSE) were used for quantitative comparison between CS-free breathing reconstruction with motion correction and that without motion correction. The PSNR and MSE of complete image and ROI for CS-free breathing motion corrected and CS-free breathing were calculated as where is a maximum pixel value of the current image having dimensions of and and are pixels being compared with current and reference cardiac phases, respectively. The ARPS block matching algorithm was used for interframe respiratory motion estimation between the reference image and the current image. Diastolic, middle of systolic and diastolic, and systolic heart phases at different respiratory states were used for both simulated images and clinical data.

The MRXCAT, a Matlab software for numerical simulation of cardiac MRI, is used for generating free breathing and breath held images. For the MRXCAT, the following parameters were used: reconstruction matrix size of 256 × 256, 24 cardiac phases in the presence of respiratory motion, with an image resolution of  mm3, TE = 1.5 ms, TR = 3 ms, and flip angle = 60°. In real free breathing cardiac cine MRI, fully sampled ECG-gated data was acquired on a Philips 1.5 T scanner (b-SSFP). Reconstruction matrix size of 256 × 256, 6 cardiac cycles with 24 cardiac phases in each cycle, and an image resolution of  mm3 were used.

For comparison of the proposed method with k-t FOCUSS, we used the following data: a short-axis MRI scan (images shown in Figure 5) was acquired using a GE 1.5 T Twin Speed scanner (R12M4) with a 5-element cardiac coil and a FIESTA/FastCARD cine SSFP sequence. Scan parameters were selected as follows: TE: 2.0 ms, TR: 4.1 ms, flip angle: 45°, FOV: 350 × 350 mm, slice thickness: 12 mm, 8 views per segment, 224 phase-encoding lines, 256 read-out samples, and 16 temporal frames. To emulate the estimation of sensitivity maps from a prescan, we acquired a separate scan (which we assume to be a prescan) with identical scan parameters and estimated sensitivity maps as follows: Half of the (high frequency) k-space samples from each coil were removed from the prescan via a smoothing filter followed by an inverse Fourier transform to obtain smoothed images for each coil. To estimate the sensitivity maps, we divided each smoothed coil image by the root sum of squares of all coil images.

The acquired data were retrospectively undersampled for acceleration rates (50% of samples), 3 (33% of samples), and 8 (12.5% of samples) using variable-density random undersampling method. Sampling masks for different acceleration rates are shown in Figure 2. The sampling mask randomly selects more samples from the low frequencies of the k-space data and fewer samples from the high frequencies of the k-space data.

Figure 3 provides a comparison of the CS-free breathing motion corrected (CS + MEMC), CS-free breathing (CS + no MEMC), and breath held for the short-axis MRI images generated from the MRXCAT simulation software at the reduction factors 2 and 8. Figure 3(a) illustrates frames 1, 5, and 12 out of 24 frames in a sequence, produced from fully sampled breath held k-spaced data. Most of the changes occur in the heart region, enclosed in the white box in (a), and are taken as a region of interest (ROI). Figure 3(b) shows the ROI, specifying left and right ventricles with endocardium and epicardium. Figures 3(c) and 3(e) show the proposed method recovery (CS + MEMC) at the reduction factors 2 and 8, respectively. Figures 3(d) and 3(f) represent the difference between breath held and estimated images for the proposed method. Figures 3(g) and 3(i) show CS + no MEMC at the reduction factors 2 and 8, respectively. Figures 3(h) and 3(j) represent the difference between breath held and estimated images with CS + no MEMC.

Images recovered by the proposed method show significant improvement as compared to the image recovery without MEMC at both reduction factors. Motion artifacts like ghosting and blurring can be seen in Figures 3(g) and 3(i) pointed by the black arrows. The proposed method eliminated ghosting and blurring effects and achieved high spatial and temporal accuracy as shown in Figures 3(c) and 3(e). The elimination of motion artifacts provides sharp endocardium and epicardium borders. This sharpness is very important in the clinical interpretation of ventricular dynamics. The improved recovery of the proposed method is also evident from difference images.

Figure 4 presents the comparison of the proposed method and CS + no MEMC for the short-axis MRI images at the reduction factors 3 and 8.

Figure 4(a) illustrates complete dataset of a diastolic, middle of diastolic and systolic, and systolic frames in a sequence, produced from clinically observed fully sampled k-spaced data. Figure 4(b) shows the ROI, enclosed in rectangular box in (a), representing left and right ventricle and epicardium and endocardium. Figures 4(c) and 4(e) show the proposed method recovery at reduction factors 3 and 8, respectively. The results of CS-free breathing without MEMC at reduction factors 3 and 8 are illustrated in Figures 4(g) and 4(i), respectively. For clinical data, images generated by the proposed method show significant improvement as compared to the image recovery without MEMC at both reduction factors. Motion artifacts like ghosting and blurring can be seen in Figures 4(g) and 4(i) pointed by black arrows. The systolic phase recovered by proposed method in (c) and (e) is very close and clear to fully sampled ROI in (b) as compared to (g) and (i) where not only are the images ghosted and blurred, but also the heart walls are displaced from their true location. Sharpness of epicardium and endocardium is also prominently visible in images recovered through the proposed method.

Figure 5 illustrates the comparison of the proposed method (CS + MEMC) and k-t FOCUSS with MEMC for the short-axis MRI dataset at reduction factor 4. Figure 5(a) shows frames 1, 13, and 10 (from left to right) out of the 16 frames in the sequence, calculated from the fully sampled breath held k-space data. Using k-space tutorial [19], motion-corrupted images are generated from the short-axis MRI dataset. Figure 5(b) presents the proposed technique reconstructions at reduction factor 4. The results for k-t FOCUSS with MEMC at reduction factor 4 are presented in Figure 5(c). The proposed method reconstruction shows significant improvement and less random noise than k-t FOCUSS with MEMC reconstruction. Furthermore, k-t FOCUSS with MEMC reconstructions contains motion artifacts (visible with bright regions), while the proposed method reconstructions are much cleaner. Tables 1, 2, and 3 provide a comparison of the proposed method and k-t FOCUSS for performance metrics such as SSIM, PSNR, and MSE. The numerical values of metrics for selected frames 1, 10, and 13 show that the proposed method outperforms k-t FOCUSS with MEMC.

A plot for PSNR at different reduction factors for CS-free breathing and CS-free breathing motion corrected is shown in Figure 6. It is drawn for recovered images shown in Figure 4. Dotted lines denote PSNR over the ROI and the solid line shows it over the entire image. The curves show that CS-free breathing motion corrected (the proposed method) is far better than the CS-free breathing at all reduction factors for both the full reconstruction and the reconstructions of ROI. Even at higher reduction factor like 12, PSNR for ROI is 4 db better in CS-free breathing with MEMC as compared to CS-free breathing without MEMC.

To show how the recovered images, with and without MEMC, are similar to the fully sampled images, we used SSIM. A plot for SSIM at different reduction factors for CS-free breathing and CS-free breathing motion corrected is shown in Figure 7. The plot is drawn for clinical data of Figure 4. Solid lines denote SSIM over the ROI and dotted line shows it over the entire image (full image). The curves illustrate that the images recovered with the proposed method are more similar to a gold standard as compared to CS-free breathing without MEMC at all reduction factors.

A plot for reconstruction mean square error (MSE) at different reduction factors for CS-free breathing and CS-free breathing motion corrected is illustrated in Figure 8. The plot is drawn for clinical data of Figure 4. Solid lines denote MSE over the ROI and dotted line shows it over the entire image. The curves illustrate that images recovered with proposed method have smaller MSE in comparison with CS-free breathing without MEMC at all reduction factors for both entire region and ROI.

6. Discussion

Interframe motion estimation and compensation has been used to improve cardiac images in breath held condition and at the same respiratory states. The proposed CS-free breathing motion corrected framework exploits temporal redundancy among the same cardiac phases at different respiratory states (free breathing). The source of improvement is the use of a novel smooth -norm approximation and inclusion of sparse residual in the basic cost function of (3).

This article suggests a CS-free breathing motion corrected recovery from cardiac cine MR data acquired below the Nyquist sampling rate. The proposed framework is a combination of interframe motion estimation technique and compressed sensing. The ARPS block matching algorithm is used for respiratory motion estimation and compensation. A simple gradient descent algorithm, with the -norm approximated by hyperbolic tangent function, is used to reconstruct the motion corrected images. The proposed method is implemented for 2D ECG-gated cardiac cine MRI for both simulation and clinical data. Implementation parameters are discussed in an experimental setup section.

The proposed method requires higher computations due to an iterative nature of the algorithm. To estimate motion corrected image, the algorithm needs to compute motion operators and alternatively multiple times. The ARPS introduces interpolation error during the prediction process of motion corrected image. There is a need to investigate motion estimation schemes with reduced interpolation error during the process of motion estimation and correction. This error is inherent in all estimation schemes and is not a topic in this article.

In the presented scheme, motion corrected images are produced from a fixed reference frame. In the future, motion estimation can be done from adjacent frames in both forward and backward direction and other CS recovery approaches can be used. Further research might be required to find its usage in 3D dynamic cardiac MRI. The proposed method will require modification for arrhythmic patients because heart rate variability is not considered in this work. Similar cardiac phases at different respiratory states are chosen visually from a sequence of cardiac MRI frames. Data binning might be used for selection of cardiac phases at different respiratory states.

7. Conclusion

In this article, we proposed a method for respiratory motion correction in ECG-gated free breathing cardiac MRI. Interframe motion estimation was used to estimate the respiratory motion between the same cardiac phases, but at different respiratory states. The block matching algorithm was used for MEMC. A gradient descent algorithm based on flexible -norm approximation was used for the recovery of MR images free from motion artifacts and close to the true MR images. Standard metrics like SSIM, PSNR, and MSE at different reduction factors were observed to be superior for the proposed method as compared to the results obtained without MEMC and k-t FOCUSS.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.