Compressed sensing reconstruction for whole-heart imaging with 3D radial trajectories: a graphics processing unit implementation

Abstract

A disadvantage of three-dimensional (3D) isotropic acquisition in whole-heart coronary MRI is the prolonged data acquisition time. Isotropic 3D radial trajectories allow undersampling of k-space data in all three spatial dimensions, enabling accelerated acquisition of the volumetric data. Compressed sensing (CS) reconstruction can provide further acceleration in the acquisition by removing the incoherent artifacts due to undersampling and improving the image quality. However, the heavy computational overhead of the CS reconstruction has been a limiting factor for its application. In this article, a parallelized implementation of an iterative CS reconstruction method for 3D radial acquisitions using a commercial graphics processing unit is presented. The execution time of the graphics processing unit-implemented CS reconstruction was compared with that of the C++ implementation, and the efficacy of the undersampled 3D radial acquisition with CS reconstruction was investigated in both phantom and whole-heart coronary data sets. Subsequently, the efficacy of CS in suppressing streaking artifacts in 3D whole-heart coronary MRI with 3D radial imaging and its convergence properties were studied. The CS reconstruction provides improved image quality (in terms of vessel sharpness and suppression of noise-like artifacts) compared with the conventional 3D gridding algorithm, and the graphics processing unit implementation greatly reduces the execution time of CS reconstruction yielding 34-54 times speed-up compared with C++ implementation.

Copyright © 2012 Wiley Periodicals, Inc.

Figures

Figure 1

3D radial reconstruction using compressed sensing. The iterative process consists of two steps of data consistency and thresholding. The image is updated to reduce the l2-norm error between the measured data and the k-space of the image estimate in the data consistency step, and to enforce the sparsity of the image estimate in the transform domain in the thresholding step. The final image is obtained as the result of the iterative process.

Figure 2

CUDA grid hierarchy and thread assignment: A grid, which consists of multiple threads, is generated once the device kernel is invoked. Each projection line of the 3D radial trajectory is assigned to one block of threads. Each thread in a block corresponds to a 3D radial sample point in the same projection line. The total number of projections is equal to the total number of blocks. This example shows a thread assignment of a 3D radial trajectory with (Ns, Np, Ni) = (8, 3, 2).

Figure 3

Thread assignment strategies for implementation of a gridding algorithm in CUDA programming: (a) radial point driven assignment, (b) Cartesian point driven assignment. Cumulative memory writes can be observed in the radial point driven assignment. The central grid point has a larger workload than the outer grid point in the Cartesian point driven assignment.

Figure 4

Comparison of conventional 3D gridding reconstruction vs. 3D iterative CS reconstruction with different sparsity regularization (image domain and wavelet domain) for a 3D radial acquisition using four different sampling densities (40%, 30%, 20%, 10%, and 7.5%). The number of iterations were 3000 and 500 for CS with image domain sparsity and wavelet domain sparsity, respectively. For high sampling densities all three reconstruction methods yield comparable image qualities. For lower densities, both CS reconstructions provide superior image qualities compared with the gridding algorithm, while CS with image domain sparsity shows better results at sharp edges and CS with wavelet domain sparsity is better at smooth surfaces. The normalized mean-squared errors are also included at the right bottom of the images.

Figure 5

CS reconstruction with image domain regularization for a phantom imaged with 3D radial with sampling density of 7.5% at different numbers of iterations, initiated with the conventional gridding reconstruction. The streaking artifacts are gradually removed with some blurring up to 500 iterations, however, with additional iterations the streaking artifacts are suppressed with improved sharpness.

Figure 6

Example slices of axial views from 3D whole-heart images reconstructed with the conventional 3D gridding reconstruction and iterative CS reconstruction (with 1000 iterations for image domain regularization and 500 iterations for wavelet domain regularization) for different sampling densities. For all sampling densities, CS reconstructions have less high-frequency streaking artifacts, and the improvement in the image quality is more distinct at lower sampling densities.

Figure 7

Example slices of sagittal views from 3D whole-heart images reconstructed by conventional 3D gridding reconstruction and iterative CS reconstruction (with 1000 iterations for image domain regularization and 500 iterations for wavelet domain regularization) for different sampling densities. For all the sampling densities, CS reconstructions have less high-frequency streaking artifacts, and the improvement in the image quality is more distinct at lower sampling densities.

Figure 8

An example slice from 3D data-set (sampling density = 6.8%) of the coronary arteries reconstructed using CS with image domain regularization at different iterations. The high-frequency artifacts are gradually removed throughout the iterations up to 500 iterations. Slight improvement was observed after 500 iterations, but it was less prominent that the phantom case (Figure. 5).

Figure 9

Reformatted images of the RCA with isotropic resolution of (1.0 mm)3 from wholeheart 3D radial data with three sampling densities (40%, 20% and 10%) by the iterative CS reconstruction with image domain regularization and 1000 iterations on GPU. The actual scan time with sampling density of 40% was 7 minutes 28 seconds with the navigator gating efficiency of 54%. The RCA is clearly visualized with the CS reconstruction for all sampling densities, while slight blurring of the image and residual artifacts are observed at low sampling density (10%).