Magnetic Resonance Imaging xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Magnetic Resonance Imaging journal homepage: www.mrijournal.com

Original contribution

Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing Rajikha Raja ⁎, Neelam Sinha International Institute of Information Technology, Bangalore, India

a r t i c l e

i n f o

Article history: Received 27 March 2013 Revised 14 October 2013 Accepted 1 December 2013 Available online xxxx Keywords: DCE-MRI CS-based image reconstruction Gradient priors Adaptive k-space sampling

a b s t r a c t The critical challenge in dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is the tradeoff between spatial and temporal resolution due to the limited availability of acquisition time. To address this, it is imperative to under-sample k-space and to develop specific reconstruction techniques. Our proposed method reconstructs high-quality images from under-sampled dynamic k-space data by proposing two main improvements; i) design of an adaptive k-space sampling lattice and ii) edgeenhanced reconstruction technique. A high-resolution data set obtained before the start of the dynamic phase is utilized. The sampling pattern is designed to adapt to the nature of k-space energy distribution obtained from the static high-resolution data. For image reconstruction, the well-known compressed sensing-based total variation (TV) minimization constrained reconstruction scheme is utilized by incorporating the gradient information obtained from the static high-resolution data. The proposed method is tested on seven real dynamic time series consisting of 2 breast data sets and 5 abdomen data sets spanning 1196 images in all. For data availability of only 10%, performance improvement is seen across various quality metrics. Average improvements in Universal Image Quality Index and Structural Similarity Index Metric of up to 28% and 24% on breast data and about 17% and 9% on abdomen data, respectively, are obtained for the proposed method as against the baseline TV reconstruction with variable density random sampling pattern. © 2014 Elsevier Inc. All rights reserved.

1. Introduction Dynamic magnetic resonance imaging (MRI) is used in various clinical applications to study the dynamics involved in various physiological processes, such as cardiac imaging, functional MRI, dynamic contrast-enhanced magnetic resonance imaging (DCEMRI). This could involve finding changes in tissue manifestations, such as signal intensity, over time. In dynamic imaging, a time series of MRI volumes are acquired, hence capturing the required dynamics for study. Such studies require images with high spatio-temporal resolution, which is however, not feasible due to the inherent limitations of MRI. In order to address the problem outlined, only a fraction of k-space data is acquired during the dynamic phase of the acquisition. It is designed such that the time saved by reducing acquired data would lead to adequate temporal resolution for the study. Various methods have been proposed to reduce the number of sampling points by designing optimal sampling patterns that lead to increased quality of image reconstruction [1–14]. Most of these sampling patterns, utilize prior knowledge of the signal of interest. ⁎ Corresponding author. Tel.:+91 9611703878. E-mail addresses: [email protected] (R. Raja), [email protected] (N. Sinha).

However the direct image reconstruction with reduced k-space data leads to undesirable poor spatial resolution. To overcome this, several reconstruction techniques [15,16] have been reported to improve the spatial resolution of the reconstructed images, using reduced k-space data. The main issues in dynamic MRI are as follows: (1) Design of optimal k-space sampling pattern. (2) Optimal image reconstruction technique. (3) Incorporation of prior information from the data available before the dynamic phase. The success of the techniques used to address the above mentioned issues is measured by comparing images that are obtained using all k-space data as against images reconstructed using reduced k-space data. Clearly, high-quality image reconstruction is desired with as reduced data as possible. Quantitative measures that check the quality of the image for a given time instant, such as Mean Square Error (MSE), Structural Similarity Index Metric (SSIM) [17], and Universal Image Quality Index (UIQI) [18] are typically used. MSE is a most widely used objective metric which calculates the error metric pixel-wise but fails to model the human visual system. SSIM and UIQI can provide a better approximation to

http://dx.doi.org/10.1016/j.mri.2013.12.022 0730-725X/© 2014 Elsevier Inc. All rights reserved.

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

2

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

the perceived image distortion as they measure the change in structural information and contrast. While MSE measures any departure from the original indiscriminately, SSIM and UIQI measure the visually perceived distortions. Additionally, in DCE-MRI, the intensity uptake curves through a region of interest (ROI), typically a tumor, are compared. In the present study we examine aspects that affect the image quality in scenarios of dynamic contrast-enhanced MRI. In the proposed work, we attempt to design a customized sampling pattern and to develop an optimal image reconstruction technique. The proposed sampling pattern follows the k-space energy distribution of the static high-resolution data to select random samples from k-space data. Next, the proposed reconstruction technique attempts to incorporate the gradient information from the high-resolution data thus modifying the regions of sharp transitions in the compressed sensing (CS)-based total variation (TV) minimization image reconstruction [19]. The gradient information and the time elapsed in dynamic acquisition are used to derive a weighting factor, “ζ(x, y, z, t)” which decides on the pixels to be used from the reference image to sharp the edges. An overview of the proposed method is shown in Fig. 1. 2. Related work The importance of dynamic MR imaging in clinical practice is evident in the volume of work carried out over the last two decades. Various methods have been proposed with the common objective of addressing the main issue in dynamic MR imaging of reducing the data acquisition time: through the design of optimal sampling lattice, efficient reconstruction technique and incorporating prior information from the available data. The initial methods were focused on the optimal selection of samples in k-space; this was followed up with conventional linear reconstruction techniques. However, over the last decade several non-linear reconstruction techniques were proposed to estimate the unacquired data more reliably. Alongside, methods have been proposed onhow to incorporate prior information from the data available in pre-dynamic phase to obtain highquality image reconstruction. The main intent of designing an optimal sampling pattern involves: simultaneous maximum energy compaction and minimum data acquisition. Reduced data encoding methods, extensively used [1,4], work by capturing a significant fraction of total k-space energy with less number of samples at an adequate temporal resolution. This is done to achieve satisfactory reconstruction of the physiological dynamics. As k-space energy is concentrated at lower frequencies, most of the sampling patterns prioritize low-frequency samples as

against that of high frequency. In case of dynamic imaging, the main assumption is that no significant changes happen in the high frequencies during the dynamic phase and that most of the information is mainly contained in the data at low frequencies. Therefore, if prior information related to high frequencies is collected during predynamic phase, only less number of samples corresponding to the low frequencies is required to be collected during dynamic phase. Methods such as keyhole imaging [1] and reduced encoding imaging by generalized series reconstruction (RIGR) [2] follow reduced data encoding approach for reduced data acquisition. In both the reconstruction methods, high-resolution dynamic images are obtained with the use of high-resolution reference image which is acquired before the dynamic portion of the study begins. During the dynamic phase, only a subset of samples from a rectangular window symmetric about kx axis is collected, to achieve higher temporal resolution. In keyhole imaging, images are reconstructed from the under-sampled dynamic frames by simply substituting the samples from reference image for the unacquired portion of dynamic data and then performing a linear reconstruction. On the other hand, in RIGR, the missing samples in under-sampled dynamic frames are recovered using a generalized series model, where the basis functions are based on reference data and finally thedynamic images are reconstructed using linear reconstruction technique. While the keyhole method suffers from artifacts due to discontinuities in k-space incurred by direct substitution, RIGR method also exhibits truncation artifacts due to partial encoding of the data. Many variations and improvements have been made to both the keyhole and RIGR techniques, such as keyhole with elliptical and rhomboid windows, radial keyhole [3] and Two-reference RIGR (TRIGR) [4]. In several other methods, [5,9], categorized as datasharing methods, the missing data can be obtained from adjacent frames by assuming that the dynamics in an image sequence change only by a small amount from frame to frame. Thus the acquired data are interleaved to obtain high-quality image reconstruction with adequate temporal resolution. BRISK [7] and TRICKS [9] indeed follow the same principle of acquiring low-frequency k-space data more frequently than high-frequency k-space data. BRISK was adapted for cardiac acquisitions while TRICKS was applied for 3D DCE-MRI acquisitions. To explore higher achievable under-sampling factors, without compromising image quality, new algorithms that use non-linear reconstruction techniques have been proposed in recent times. Compressed Sensing theory [20,21] provides a novel way to accomplish this goal. Various applications of CS theory to dynamic MRI such as k-t SPARSE [15] and k-t FOCUSS [22], have been successfully demonstrated.

Fig. 1. Flowchart of the proposed method.

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

3

Fig. 2. k-Space energy adaptive sampling pattern algorithm (the size of the volume for every time point is referred as fr_en x ph_en x no_slices, where fr_en denotes the number of frequency encodes, ph_en denotes the number of phase encodes and no_slices denotes the number of slices per volume). Spatial coordinates are referred as (x, y, z) and k-space coordinates are referred as (kx, ky, kz).

A popular MR image reconstruction approach that utilizes undersampled k-space is the CS-based TV minimization reconstruction, explored in Ref. [19]. In the method proposed in Ref. [19], random sampling is carried out and image reconstruction is obtained by solving the constrained optimization problem, Minimizekψmk1 such that k F u m−yk2 bϵ

ð1Þ

where m is the reconstructed image, y is the measurements in k-space, Fu is the under-sampled Fourier transform, ψ is the sparsifying transform used. Variations of CS-based TV minimization [19] have been explored in Refs. [10,14]. The variations involve (i) modifying sampling patterns based on the nature of k-space energy distribution, (ii) incorporating prior information, and (iii) using spatio-temporal correlations of the dynamic data. One of the main factors influencing the quality of CS-based image reconstruction is the choice of sampling pattern. Sampling patterns could be broadly categorized as (i) uniform sampling and (ii) random sampling schemes. Cartesian sampling is most commonly used, but has less incoherence and hence not ideal for CS MRI reconstruction. On the other hand non-Cartesian

trajectories like radial, spiral yield good results but have implementational limitations. In order to have higher incoherence, various sampling strategies to obtain random sampling patterns have been proposed for CS-based reconstruction. The most popular is the variable density sampling scheme where a greater number of high-energy low-frequency components are selected compared to low-energy higher-frequency components, to reduce aliasing. Lustig et al. [19] proposed variable density random sampling which determines sampling mask based on a probability density function (PDF). However the model used for the PDF is nonadaptive and lacks the ability to follow the corresponding k-space energy distribution. The k-space corresponding to MR images of different anatomies show different characteristics. These specific characteristics can be exploited in optimal k-space sampling leading to customized k-space trajectories for the respective imaging scenario. Thesampling trajectory proposed in Ref. [14], exploits the characteristic energy distribution of k-space for a particular anatomy. This is accomplished by considering the mean of the k-space energy in a sub-region to calculate thedensity of samples to be selected in the sub-region. However, it would be ideal if the sampling pattern could be customized for the dynamic images at hand, on a case by case basis, with no generalization. On the other hand, Ravishankar and Bresler [11] proposed an adaptive under-sampling design to

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

4

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

Fig. 3. Blocks of k-space energy distribution matrix.

redistribute the phase encodings based on k-space errors between the reconstructed and full-sampled data. However, this method may not be effective when the initial trajectory included enough central phase encodings. This method was improved in Ref. [12] by avoiding the dependence on initial trajectory and the k-space errors used in Ref. [11] were substituted with Bayesian Posterior entropy [13]. Incorporation of prior information has been handled differently across research groups. While methods like keyhole imaging directly substituted values from pre-dynamic data, RIGR

utilized it in the modeling of the generalized series. CS-based reconstruction of DCE-MRI was also demonstrated in Ref. [10], where the sparse nature of DCE data set in the finite difference transform along the time dimension is used. However, in a recently proposed method, EESTCR [16], edge information is used to enhance the TV minimization constraint and to improve the sharpness of the edges by adding an edge matching function based on a reference image. Improvement in results was shown with EESTCR for an under-sampling factor of 4.

Fig. 4. TV-gradient-kspace reconstruction algorithm (the size of the volume for every time point is referred as fr_en x ph_en x no_slices, where fr_en denotes the number of frequency encodes, ph_en denotes the number of phase encodes and no_slices denotes the number of slices per volume). Spatial coordinates are referred as (x, y, z) and k-space coordinates are referred as (kx, ky, kz).

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

5

Fig. 5. Representative images from the data sets used. (a) Slice 21 of first volume of breast data. (b) Slice 11 of first volume of abdomen data (arrow points at tumors in panels (a) and (b)).

3. Proposed method In the proposed method, we work on the following two aspects: (1) customizing the k-space sampling pattern to the imaging scenario and (2) incorporation of gradient information obtained from the high-resolution static data. The method mainly hinges on the assumption that the nature of the k-space does not vary much during the course of the dynamic phase [1,4]. 3.1. Designing of k-space sampling pattern In the proposed method, an imaging scenario-customized sampling pattern is designed which adapts to the energy distribution in the k-space and hence is named as “k-space Energy Adaptive Sampling Pattern.” The k-space of the static high-resolution data is considered for designing the sampling pattern. The algorithm is detailed in Fig. 2. The size of the volume for every time point is referred as fr_en x ph_en x no_slices, where fr_en denotes the number of frequency encodes, ph_en denotes the number of phase encodes and no_slices denotes the number of slices per volume. The coordinates (x, y, z) referred in this paper denotes the spatial coordinates and (kx, ky, kz) denotes the k-space coordinates. The division of k-space energy distribution matrix into “N” annular blocks is shown in Fig. 3 with each block coloured differently.

directions. In the proposed reconstruction technique, initially the gradient values are calculated slice-wise for the high-resolution static volume to get a gradient image for each of the slice. In order to get the full range of direction, gradient images in both horizontal and vertical directions are computed. Thus, gradient magnitudes measured in horizontal (GRx) and vertical (GRy) directions are summed up to get the gradient GR of an image. The gradient images are computed so that large gradient values become possible edge pixels and this in turn gives more weighting for calculating ζ which is explained below in Eq. (2). The gradient of the static image, IRef (x, y, z) is computed as, GRðx; y; zÞ ¼ absðGRx ðI Ref ðx; y; zÞÞÞ þ absðGRy ðIRef ðx; y; zÞÞÞ

The information in GR(x, y, z) is used to compute a weighting factor “ζ” that is designed to be high at those image regions where the gradient GR(x, y, z) is high. Besides, ζ also includes a time factor “tζ” that varies inversely with time elapsed. A variable “time_elapsed” is defined, whose value ranges from 1,2,3…T, where T is the number of time points for which the data are acquired. Reciprocal of time elapsed is added with the image gradient in order to decrease the weighting factor “ζ” accorded to the static image, as the time since the start of the dynamic phase increases.

t ζ ðtime elapsedÞ ¼ 1=time elapsed

3.2. Gradient-based reconstruction technique The term “gradient” used in this context refers to the image gradient which is used to measure the change in intensity with direction by computing the derivatives in horizontal and vertical

Table 1 Data description. Data set

Anatomy

No. of time points

Size

1 2 3 4 5 6 7

Breast Breast Abdomen Abdomen Abdomen Abdomen Abdomen

4 4 5 4 9 5 5

256 256 256 512 256 256 256

× × × × × × ×

256 256 256 512 256 256 256

× × × × × × ×

36 20 36 20 36 20 20

ð2Þ

ð3Þ

where time_elapsed = 1,2,3…T. Thus, the weighting factor ζ(x, y, z, t) is given as: ζ ðx; y; z; t Þ ¼ GRðx; y; zÞ þ t ζ ðt Þ

ð4Þ

The k-space in dynamic phase is acquired using the variable density random sampling pattern as outlined in Ref. [19]. The CS-based TV minimization image reconstruction (ITV) is utilized to obtain an estimate of the images in the dynamic phase. However, these pixel values are altered to incorporate the gradient prior as well as time elapsed, in the proposed image reconstruction (ITV _ Gradient) as, ITV

Gradient ðx; y; z; t Þ

¼ ζ ðx; y; z; t Þ  IRef ðx; y; zÞ þ ½1−ζ ðx; y; z; t Þ  ITV ðx; y; z; t Þ

ð5Þ

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

6

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

Fig. 6. Representative image from the synthetic data. (a) Single slice of synthetic data at first time point (tumors are numbered). (b) Uptake curves of 6 tumors over 8 time points.

Fig. 7. Comparison of reconstructed slice from the third volume of synthetic data using various methods for under-sampling factor of 5. (a) Keyhole method. (b) DynCS. (c) TVVDRS. (d) TV-gradient. (e) Comparison of uptake curves for tumor 2 for different reconstruction methods. (f) Comparison of mean tumor intensities at second time point to show the robustness to noise of various methods. Comparison of visual quality of tumor 3 in the reconstructed slice (shown in yellow boxes in panel (a)–(d)) from the third volume of synthetic data for under-sampling factor of 5 (zoomed inside the yellow box). (g) Keyhole method. (h) DynCS. (i) TV-VDRS. (j) TV-gradient.

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

7

Fig. 8. Comparison of ζ maps for breast data set (slice 21) at different time points with the corresponding sampling masks for different under-sampling factors. (a) Static highresolution image. (b) Sampling mask for under-sampling factor of 20. (c) Sampling mask for under-sampling factor of 15. (d) Sampling mask for under-sampling factor of 10. (e)– (h) Fully sampled images at successive time points from 1 to 4. (i)–(l) ζ maps at successive time points from 1 to 4 for an under-sampling factor of 10.

The pseudo-code for the gradient-based reconstruction algorithm along with the proposed sampling pattern (hereafter referred to as TV-gradient-kspace reconstruction) is shown in Fig. 4.

points. As mentioned earlier, in our experiments, the first data volume is taken as the high-resolution acquisition. To simulate the dynamic under-sampled k-space data, the respective k-space is under-sampled by the corresponding under-sampling factor.

4. Data used 5. Results All simulations are carried out in MATLAB (MathWorks, Natick, MA). The proposed method is tested on seven real volumes of data of which two are breast and five are abdomen data sets. The data sets are acquired as seven scans on patients after informed consent on GE 1.5T Excite Scanner with the following scan parameters: CA = 10–15 cc gadolinium contrast; TE/TR = 1.4 ms/3.9 ms; acquisition matrix = 256 × 192 × 20 (36); FOV = 20/36 cm (breast/abdomen); and slice thickness = 3/4.2 mm. 1) Breast data: One of the time series is available for 4 time points, where the 21st slice has a tumor, as shown in Fig. 5(a). The other data set is available at 4 time points. 2) Abdomen data: These data sets are time series of more than 4 time points, of different sizes. Table 1 details the data set sizes. A representative slice of this volume is shown in Fig. 5(b). As proof of concept, the experiment is carried out on synthetic data of a single slice of size 256 × 256, containing six tumors of various sizes and shapes that enhance at different rates over 8 time points. Fig. 6(a) shows the original synthetic image at first time point, with the tumors numbered from 1 to 6 and Fig. 6(b) shows the uptake curves for the 6 tumors from the original image over 8 time

The experiments are carried out at under-sampling factors of 5, 8 and 10, corresponding to data availability of 20%, 12.5% and 10%, respectively. In the trials reported below, data availability of 10% and 12.5% is considered. In each of the cases, we present the comparison between the following methods, (i) Baseline TV reconstruction with variable density random sampling [19] (hereafter referred to as TV-VDRS) (ii) TV-gradient-kspace (iii) Keyhole method [1] (iv) DynCS method [10] Reconstruction of the real and synthetic data sets using TV-VDRS is implemented using Lustig's SparseMRI V0.2 set of functions available in Ref. [23] and the reconstructed volume from TV-VDRS is further enhanced using proposed methods. The reconstruction is quantitatively and qualitatively compared against TV-VDRS to show improvements. Additionally, relevant methods such as keyhole and DynCS which are proposed for reconstruction of DCE-MRI data set using highresolution reference images are also implemented and the outputs of which are compared with that of the proposed method. A comparison

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

8

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

Fig. 9. Representative results on breast data volume (size: 256 × 256 × 36, for 4 time points) Comparison of reconstructed images (at time point 3) obtained using 10% of k-space data. (1) Reconstructed images. (2) SSIM maps. (3) UIQI maps. (4) Tumor region (shown inside yellow box) zoomed in.

Fig. 10. (a) Comparison of scan line profiles for third volume (x-axis: scan line index; y-axis: signal intensity) for original, keyhole, DynCS, TV-VDRS and TV-gradient-kspace. (b) Comparison of uptake curves through the tumor (x-axis: dynamic imaging time; y-axis: mean signal intensity) for original, keyhole, DynCS, TV-VDRS and TV-gradient-kspace (both panels (a) and (b) using 10% of k-space data for slice 21 of breast data).

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

9

Fig. 11. Representative results on abdomen data volume (size: 256 × 256 × 36, for 5 time points). Comparison of reconstructed images (at time point 2) obtained using 12.5% of k-space data. (1) Reconstructed images. (2) SSIM maps. (3) UIQI maps. (4) Region of interest (shown inside the yellow box) zoomed in.

of the scan lines and the mean signal intensity curves through the ROI is also presented. Validation of the proposed method is carried out using the following quantitative measures: (a) Mean Square Error, MSE (for 8-bit grayscale image); (2) Structural Similarity Index Metric, SSIM as described in Ref. [17] (range is 0–1); and (3) Universal Image Quality Index as described in Ref. [18] (range is 0–1) 5.1. Synthetic data analysis Different levels of Gaussian noise are added to the synthetic data and the performance of TV-gradient reconstruction is evaluated by comparing the output image with that of other reconstruction methods such as keyhole method, DynCS and TVVDRS. Noise levels of 2%, 3% and 5% are tested for evaluating the robustness of the proposed algorithm for noise. Five percent noise added to the synthetic data which is reconstructed with various methods is compared in Fig. 7(a)–(d). Signal uptake curves for tumor 2 over the entire dynamic phases for various reconstruction methods are compared with that of the original as shown in Fig. 7(e). In the uptake curves, TV-VDRS and TV-gradient reconstruction seems to follow the original curve closely as compared to that of the other methods. The robustness of TV-Gradient method to different levels of noise is evaluated which is illustrated in Fig. 7(f). For different noise levels

of 2%,3% and 5%, the mean signal intensity of tumor 2 in the TVgradient reconstruction matches closely with the original image whereas keyhole and DynCS deviate more from the original and this proves the consistent performance irrespective of increased noise levels. One of the tumors is considered (tumor 3) for analyzing the visual quality of reconstruction as shown in Fig. 7(g)–(j). Tumor 3 which is the most complicated shaped tumor is zoomed in this case to show the better quality in TV-gradient reconstruction. While keyhole and TV-VDRS reconstruction lost the complete structure of the tumor and affected by noise, DynCS regained the structure of the tumor but still affected by noise. In this case, TV-gradient performs better in reproducing the structure of the tumor and less affected by noise. Moreover edges are well preserved for all the six tumors with TV-gradient whereas edges of tumor 1, tumor 2 and tumor 6 are blurred in keyhole, DynCS and TV-VDRS respectively. In keyhole reconstruction, the entire image shows aliasing artifacts while DynCS and TV-VDRS are affected by noise. 5.2. Sampling masks and ζ maps The sampling masks generated using the proposed sampling pattern are illustrated in Fig. 8(a)–(d). For a given k-Space of a highresolution static image, the corresponding sampling mask is

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

10

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

Fig. 12. (a) Comparison of scan line profiles for fourth volume (x-axis: scan line index; y-axis: signal intensity) for original, keyhole, DynCS, TV-VDRS and TV-gradient-kspace. (b) Comparison of uptake curves through the ROI (x-axis: dynamic imaging time; y-axis: mean signal intensity) for original, keyhole, DynCS, TV-VDRS and TV-gradient-kspace (both panels (a) and (b) using 12.5% of k-space data for slice 11 of abdomen data).

generated for the required under-sampling factors. Fig. 8(a) shows the slice 21 of high-resolution static volume of breast data and Fig. 8(b)–(d) show the corresponding sampling masks for under-sampling factors of 20, 15 and 10. For different under-sampling factors, the density of samples varies based on the energy of the k-space. The samples selected are then reconstructed using the CS-based TV minimization reconstruction. The reconstructed image is then enhanced using the calculated ζ values. The ζ maps calculated for every time points are shown in Fig. 8(j)–(l) corresponding to the slice 21 of breast volume at 4 time points (shown in Fig. 8(e)–(h)). The ζ values are high for boundary regions and edges which is well intended. 5.3. Real-data analysis For the breast data set, the images reconstructed using keyhole, DynCS and baseline TV-VDRS are compared with TV-gradient-kspace in Fig. 9(a)–(d) with the corresponding SSIM maps and UIQI maps in Table 2 Image quality metrics for breast and abdomen data reconstructed with 10% k-space data. Average of 2 breast data sets TV_VDRS TV_kSpace TV_Gradient_kSpace Average of 5 abdomen data sets TV_VDRS TV_kSpace TV_Gradient_kSpace Average of 7 data sets TV_VDRS TV_kSpace TV_Gradient_kSpace

MSE

SSIM

UIQI

247.051 172.31 97.7785

0.590572 0.7283 0.86965

0.37291 0.511517 0.691721

MSE

SSIM

UIQI

41.1002 22.5963 27.4371

0.845899 0.927003 0.960808

0.42619 0.525038 0.617757

MSE

SSIM

UIQI

99.9433 65.3715 47.5347

0.772949 0.870231 0.934763

0.410967 0.521175 0.638889

Values for breast data is averaged 2 time series (180 + 100 reconstructed images). Values for abdomen data is averaged over 5 time series (216 + 100 + 360 + 120 + 120 reconstructed images). Last row shows averaged metrics over all the seven real data sets.

Fig. 9(e)–(h) and (i)–(l) respectively. Fig. 10 illustrates the comparison of scan lines and the intensity uptake curves. Similarly, for the abdomen data set, the images reconstructed using baseline TV-VDRS are compared with TV-gradient-kspace in Fig. 11(a)–(d) with the corresponding SSIM maps and UIQI maps in Fig. 11(e)–(h) and (i)–(l) respectively. Fig. 12 illustrates the comparison of scan lines and the intensity uptake curves. Table 2 records the quantitative evaluation of the reconstructed images, over the well-known quality indices: MSE, SSIM and UIQI. 5.4. Computational complexity All implementations were done on a Windows machine with 3.2 GHz Intel Xeon Processor and 8-GB RAM. The computational complexity of the proposed algorithm is the summation of the computational requirements of its three sub-parts: 1. k-space adaptive sampling pattern algorithm, which is of complexity O(v) for “v” voxels in the sampling plane. 2. CS-based TV minimization reconstruction by Lustig [23] 3. Computation of gradient-based reconstruction, which is of complexity O(v) for a single slice. The time requirements of the reconstruction methods are compared in Table 3. The proposed method is built upon CS-based TV minimization reconstruction in which the optimization is done via iterative non-linear conjugate gradient algorithm which is computationally intensive.

Table 3 Comparison of complexities of reconstruction methods. Method

TV-gradient- TV-VDRS kspace

Time taken for reconstructing 1.1573 sec a single slice (machine configuration as given in paper)

DynCS

Keyhole

1.0128 sec 1.0779 sec 0.1385 sec

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

11

Fig. 13. Comparison of the improvement in image quality indices of TV-gradient-kspace as against TV-VDRS across breast and abdomen data sets (x-axis: anatomies; y-axis: improvement in SSIM/UIQI) using 10% of k-space data.

The sparsifying transform used for the CS-based TV minimization algorithm is finite differences in the spatial domain. As Lustig et al. [19] pointed out in their paper, the number of samples required for a good reconstruction should be roughly two to five times the number of sparse coefficients. The reconstructed images shown in this article are representative results of the data sets. Detailed results for the entire data volumes are made available online at www. resultstvgradientkspace.wordpress.com. 6. Discussion The proposed method aims to improve the well-known TV minimization-based reconstruction for DCE-MRI applications. The proposed improvement is carried out by incorporating information from the available high-resolution static image. In this work, we utilize the a priori information in k-space energy distribution of the highresolution static image to design the adaptive sampling pattern. We also introduce a weighting factor that decides, for a dynamic time point, how the current reconstruction can be combined with the available prior. 6.1. Quantitative analysis measures The reconstruction is analyzed through quantitative measures such as MSE,SSIM and UIQI. The performance improvement is shown across quality metrics to prove the consistency of the reconstruction performance. The obtained quality indices such as MSE, SSIM and UIQI, are shown in Table 2. MSE objectively quantifies the error signal, giving the average performance. However, outliers could adversely affect the quality index. SSIM and UIQI penalize the distortions in brightness, contrast and structure of the reconstructed images. Both SSIM and UIQI measure image distortion as a combination of three factors: loss of correlation, luminance distortion, and contrast distortion. It must be seen that the values show a significant improvement across the various indices, consistently over 1196 (180 + 100 breast images and 216 + 100 + 360 + 120 + 120 abdomen images) images. These values are obtained by taking the average over all slices across all the dynamic sampling time instants for each of the categories. The performance is consistent in spite of diversity in details of these imaged anatomies. 6.2. Qualitative analysis of real-data reconstruction The proposed reconstruction is qualitatively analyzed for different factors such as aliasing, blurring, robustness to noise, preserva-

tion of edges, homogeneous regions, texture and contrast. With increased under-sampling factors keyhole methods showed pronounced blurring in both the breast and abdomen data sets as this method just substitutes the missing k-space values by using those from the reference data. The SSIM and UIQI maps indicate the ability of the reconstruction method to conserve structure and contrast. The values in SSIM and UIQI maps range from 0 to 1 where 0 denotes the least similarity and 1 denotes the best similarity. The SSIM and UIQI maps in Figs. 9(2) and 11(2) show high similarity for TV-gradientkspace as compared to other methods. Specifically, the SSIM map in Fig. 11(e) shows high dissimilarity for the background region in abdomen data with keyhole reconstruction. Regions zoomed in Figs. 9(4) and 11(4) to illustrate the visibly better spatial resolution obtained in reconstruction using the proposed method, as compared to the other methods. The tumor region zoomed in Fig. 9(q) for TVgradient-kspace seems to preserve the structure and edges better than other methods. In DynCS, the difference between the undersampled k-space data and the static image is reconstructed using CSbased TV minimization and finally combined with the static data to get the full image. Moreover, under-sampling is done via variable density random sampling. The tumor region showed in Fig. 9(m)–(q) indicates the structure at the center of the tumor is well preserved in the reconstruction method using the proposed methodas compared to that obtained using DynCS. In addition, the proposed method is seen to be robust to the noise levels in the data. It is illustrated through the performance of the proposed method when higher levels of noise are added to the synthetic data. The ability of the proposed method to follow intensity transitions is illustrated in a representative scan line comparison shown in Figs. 10(a) and 12(a), for the breast and abdomen data, respectively. This improvement in structural fidelity is manifested in the improvement of about 24% in the SSIM index and about 28% in the values of the UIQI, on the breast data sets. In the abdomen data set, similar improvement in the structural indices can be seen, at about 9%for SSIM index and about 17% for the UIQI. Figs. 10(b) and 12(b) compare a representative plot of mean signal intensity values through the specified regions of interest. In the breast data, the tumor is the region through which the uptake curve is plotted. It is clear that the fidelity in reconstruction is reflected in the plots. While both the reconstruction techniques follow the trends manifested by the original mean intensity curve, it is the curve generated through the proposed reconstruction that deviates far less than the curve generated using keyhole, DynCS and TV-VDRS reconstruction. It is also observed that as the size of the

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

12

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

Fig. 14. Comparison of sampling mask generated for under-sampling factor of 10 of k-space data for slice 21 of first volume of breast data along with the respective reconstructed images and SSIM maps. (a) k-Space energy adaptive sampling mask (regions enclosed inside curly brackets shows the difference in the density of the samples in the proposed method). (b) Image reconstructed using TV-gradient-kspace. (c) SSIM map corresponding to TV-gradient-kspace. (d) Customized k-space trajectory mask. (e) Image reconstructed using customized TV reconstruction. (f) SSIM map corresponding to customized TV reconstruction. (g) VDRS mask (h) image reconstructed using TV-VDRS reconstruction. (i) SSIM Map corresponding to TV-VDRS reconstruction.

specified region increases, the differences between the curves generated using the baseline TV-VDRS and the one generated using the proposed method are similar. This is due to the smoothening nature of the TV-VDRS reconstruction.

the time required for acquiring one data set with the conventional reconstruction, could be used to acquire 10 data sets with the proposed method. 6.4. Comparison of proposed sampling design with customized sampling design

6.3. Performance of proposed method w.r.t. anatomies The improvement in the reconstruction quality is analyzed with respect to different anatomies as shown in Fig. 13. It can be seen that the improvement in results is high for breast data as compared to the abdomen data. Breast data contains more edgeinformation as compared to abdomen data which looks homogeneous. By using our proposed reconstruction method, we have demonstrated a potential reduction in acquisition time of up to 90% which translates to a saving of up to 90% of acquisition time. This in turn means that

In the customized sampling method proposed in Ref. [14], the authors attempt to determine a customized sampling pattern for each structure, separately. They work under the assumption that, images of the same structure exist with not much variation, and hence they can all follow a similar k-space distribution. They further illustrate that the k-space distribution is not typically Gaussian, as assumed by Lustig et al. In our sampling design, the pattern is customized for the dynamic images to be reconstructed with no generalization by adapting to the k-space energy distribution of the

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

13

Fig. 15. Impact of sampling schemes: comparison of reconstructed images using the proposed sampling scheme and customized TV sampling pattern under similar conditions of a fixed reconstruction procedure and a fixed under-sampling factor (representative results on breast data for the first volume obtained using 10% of k-space data). (a) TV reconstruction with k-space energy adaptive sampling (features shown in arrow marks which are perceptibly preserved better by TV reconstruction with k-space energy adaptive sampling as compared to utilizing customized TV). (b) TV with customized k-space trajectory.

corresponding high-resolution data. As a means of comparison we implemented the method proposed in Ref. [14], which we will refer to as customized TV reconstruction. In Fig. 14, the sampling mask generated by customized k-space trajectory is compared against the sampling mask proposed in this work. As can be seen, the density of selected samples is more in the regions enclosed within brackets in the proposed mask which is not in the case of customized mask. Also, the images reconstructed with the respective masks are shown in Fig. 15 for a single slice of breast data set. Improvement in spatial resolution is seen with the proposed sampling pattern as compared with the customized trajectory. This is also reflected in the SSIM maps for the same shown in Fig. 14(c) and (f). The quality indices for the reconstructed images of Fig. 14(b) and (c) are shown in Table 4 for an under-sampling factor of 10%. As seen in Table 3, the quality indices are consistently higher for the proposed method.

7. Conclusion In this paper, we have proposed enhancements of the wellknown CS-based TV minimization constrained image reconstruction and the variable density random sampling for dynamic MR imaging. The high-resolution static data acquisition is utilized to improve the regions of sharp transitions in the dynamic images. The performance of the algorithm is illustrated on real data sets of abdomen and breast image volumes. The reconstruction quality, at 10% data availability, when assessed by quality metrics shows an average improvement of 28% in breast data and 17% in abdomen data, in UIQI, as against the baseline TV reconstruction, spanning over 1196 images, in all.

Table 4 Image quality metrics for the reconstructed images shown in Fig. 14. Method

MSE

SSIM

UIQI

TV_kspace TV_Customized

0.00266 0.00316

0.60581 0.55281

0.52037 0.46394

References [1] Van Vaals JJ, Brummer ME, Thomas Dixon W, et al. keyhole method for accelerating imaging of contrast agent uptake. J Magn Reson Imaging 1993;3 (4):671–5. [2] Liang Z-P, Lauterbur PC. An efficient method for dynamic magnetic resonance imaging. IEEE Trans Med Imaging 1994;13(4):677–86. [3] Mistry NN, Pollaro J, Song J, De Lin M, Johnson GA. Pulmonary perfusion imaging in the rodent lung using dynamic contrast-enhanced MRI. Magn Reson Med 2008;59(2):289–97. [4] Hanson JM, Liang Z-P, Wiener EC, Lauterbur PC. Fast dynamic imaging using two reference images. Magn Reson Med 1996;36(1):172–5. [5] Parrish T, Hu X. Continuous update with random encoding (cure): a new strategy for dynamic imaging. Magn Reson Med 2005;33(3):326–36. [6] Riederer SJ, Tasciyan T, Farzaneh F, Lee JN, Wright RC, Herfkens RJ. MR fluoroscopy: technical feasibility. Magn Reson Med 2005;8(1):1–15. [7] Doyle M, Walsh EG, Blackwell GG, Pohost GM. Block regional interpolation scheme for k-space (brisk): a rapid cardiac imaging technique. Magn Reson Med 2005;33(2):163–70. [8] Zaitsev M, Zilles K, Shah N. Shared k-space echo planar imaging with keyhole. Magn Reson Med 2001;45(1):109–17. [9] Korosec FR, Frayne R, Grist TM, Mistretta CA. Time-resolved contrast-enhanced 3D MR angiography. Magn Reson Med 2005;36(3):345–51. [10] Ji J, Lang T. Dynamic MRI with compressed sensing imaging using temporal correlations. 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro, 2008. ISBI 2008. IEEE; 2008. p. 1613–6. [11] Ravishankar S, Bresler Y. Adaptive sampling design for compressed sensing MRI. 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBC. IEEE; 2011. p. 3751–5. [12] Duan-duan Liu XL, Liang Dong, Ting Zhang Y. Under-sampling trajectory design for compressed sensing MRI. 2012 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBC. IEEE; 2012. p. 73–6. [13] Seeger M, Nickisch H, Pohmann R, Schölkopf B. Optimization of k-space trajectories for compressed sensing by Bayesian experimental design. Magn Reson Med 2010;63(1):116–26. [14] Ilievska ES, Ivanovski ZA. Customized k-space trajectory for compressed sensing MRI. 2011 19th Telecommunications Forum (TELFOR). IEEE; 2011. p. 631–4. [15] Lustig M, Santos JM, Donoho DL, Pauly JM. k-t sparse: high frame rate dynamic MRI exploiting spatio-temporal sparsity. Proceedings of the 13th Annual Meeting of ISMRM, Seattle; 2006. p. 2420. [16] Iyer SK, DiBella EV, Tasdizen T. Edge enhanced spatio-temporal constrained reconstruction of undersampled dynamic contrast enhanced radial MRI. 2010 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. IEEE; 2010. p. 704–7. [17] Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process 2004;13 (4):600–12. [18] Wang Z, Bovik AC. A universal image quality index. IEEE Signal Process Lett 2002;9(3):81–4.

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

14

R. Raja, N. Sinha / Magnetic Resonance Imaging xxx (2014) xxx–xxx

[19] Lustig M, Donoho D, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. Magn Reson Med 2007;58(6):1182–95. [20] Candes E, Romberg J. Signal recovery from random projections. Proc. SPIE, 5674; 2005. p. 76–86. [21] Donoho DL. Compressed sensing. IEEE Trans Inf Theory 2006;52(4):1289–306.

[22] Jung H, Sung K, Nayak KS, Kim EY, Ye JC. k-t focuss: a general compressed sensing framework for high resolution dynamic MRI. Magn Reson Med 2008;61 (1):103–16. [23] Lustig M. http://www.eecs.berkeley.edu/mlustig/software/sparseMRI-v0.2.tar. gz; 2008 . [TV-VDRS].

Please cite this article as: Raja R, Sinha N, Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing, Magn Reson Imaging (2014), http://dx.doi.org/10.1016/j.mri.2013.12.022

Adaptive k-space sampling design for edge-enhanced DCE-MRI using compressed sensing.

The critical challenge in dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is the trade-off between spatial and temporal resolution due ...
5MB Sizes 1 Downloads 3 Views