http://orcid.org/0000-0003-2585-6748
Figures
Abstract
Dynamic magnetic resonance imaging (DMRI) is used to noninvasively trace the movements of organs and the process of drug delivery. The results can provide quantitative or semiquantitative pathology-related parameters, thus giving DMRI great potential for clinical applications. However, conventional DMRI techniques suffer from low temporal resolution and long scan time owing to the limitations of the k-space sampling scheme and image reconstruction algorithm. In this paper, we propose a novel DMRI sampling scheme based on a golden-ratio Cartesian trajectory in combination with a compressed sensing reconstruction algorithm. The results of two simulation experiments, designed according to the two major DMRI techniques, showed that the proposed method can improve the temporal resolution and shorten the scan time and provide high-quality reconstructed images.
Citation: Li S, Zhu Y, Xie Y, Gao S (2018) Dynamic magnetic resonance imaging method based on golden-ratio cartesian sampling and compressed sensing. PLoS ONE 13(1): e0191569. https://doi.org/10.1371/journal.pone.0191569
Editor: Dalin Tang, Worcester Polytechnic Institute, UNITED STATES
Received: June 11, 2017; Accepted: January 8, 2018; Published: January 30, 2018
Copyright: © 2018 Li et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: This work is supported by the National Natural Science Foundation of China (Grant 61671026) to Song Gao, the Natural Science Foundation of Beijing, China (Grant No. 7162112) to Song Gao.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Dynamic magnetic resonance imaging (DMRI) has been widely used to noninvasively trace the movement of human tissues and organs (e.g., cardiac cine MRI [1] and dynamic contrast-enhanced MRI [2]), and the process of drug delivery. The reconstructed dynamic images obtained via DMRI can be used to extract many quantitative or semiquantitative pathology-related parameters. These parameters contain considerable biological and pathophysiological information about tissues and organs during the development and occurrence of diseases. This information can be used for clinical disease research and diagnosis; therefore, DMRI has significant potential in clinical research and applications.
The DMRI technique uses three k-space sampling schemes: sampling based on a Cartesian trajectory [3], sampling based on a golden-angle radial trajectory [4], and hybrid sampling based on both Cartesian and radial trajectories [5]. All three schemes must repeatedly acquire all or parts of the k-space data. Sampling based on a Cartesian trajectory is the traditional sampling scheme and has the advantages of a simple pulse sequence design and image reconstruction method; however, this scheme suffers from low temporal resolution. The second sampling scheme stems from the golden-ratio angle sampling scheme proposed by Winkelmann in 2007 [6]. The golden-angle radial trajectory sampling scheme improves the temporal resolution of two-dimensional (2D) DMRI to a certain extent; however, three-dimensional (3D) DMRI still suffers from low temporal resolution. Furthermore, the radial trajectory introduces aliasing artifacts into the reconstructed images [7]. The hybrid sampling scheme is a combination of linear and radial trajectories; for example, in a 3D cylindrical k-space, Cartesian slice encoding can be used in the kz direction and radial sampling can be used in the kx-ky plane. The hybrid sampling scheme improves the temporal resolution of 3D DMRI to some extent; however, the radial trajectory part still introduces artifacts. Therefore, DMRI based on a linear sampling scheme with high temporal resolution shows significant promise for future clinical research and applications.
To realize DMRI based on a linear sampling scheme with high temporal resolution, we present two basic methods in this paper. First, because the use of the golden-ratio angle proposed by Winklemann improves the temporal resolution in radial trajectory, we introduced the golden ratio into a Cartesian trajectory. Second, compressed sensing (CS), first proposed in 2006 and widely used in many scientific fields, including fast MRI [8, 9], can be used to reconstruct relatively high-quality images from undersampled k-space data.
In this paper, we first discuss the basic theories of the golden-ratio-based linear sampling scheme and CS. Second, the proposed DMRI technique based on the golden-ratio Cartesian sampling scheme and CS is described. Third, computer simulation experiments conducted to compare the proposed and traditional methods are discussed. Finally, the future development of the proposed DMRI technique is outlined.
Theory and methods
Cartesian k-space sampling scheme based on the golden ratio
The golden ratio is an irrational number defined as . It has been widely used in mathematics, physics, architecture, and other fields. Because of its special mathematical properties, Winkelmann proposed a k-space sampling scheme based on the golden-ratio angle [6]. This scheme provides approximately uniformly distributed k-space data within an arbitrary reconstruction window. Therefore, it offers a high degree of flexibility in choosing an appropriate window length and temporal resolution, in positioning the reconstruction window, and in combining adjacent time frames. This allows the application of sliding window reconstruction to variable window lengths and arbitrary window positioning. However, radial trajectories introduce aliasing artifacts into the reconstructed images, whereas Cartesian trajectories do not. We propose a k-space sampling scheme based on the golden ratio and a Cartesian trajectory. Diagrams of 2D and 3D sampling schemes are shown in Fig 1(A) and 1(B). The blue, red, and green lines in Fig 1(A) denote the first, second, and third sampling trajectories, respectively. ky1, ky2, and ky3 denote the corresponding coordinates in the phase-encoding direction, where kyn is the coordinate of the nth trajectory in the phase-encoding direction and is calculated as
(1)
where α is the proportionality coefficient, mod(a,b) denotes the remainder of a/b, and γ (≈ 1.618) is the golden ratio. In Fig 1(B), kz is the slice-encoding direction. All the trajectories with the same kz coordinate are defined as a profile, e.g., the blue, red, and green profiles. Data in the blue profile are acquired first, followed by those in the red and green profiles. The data acquisition time of each profile is mzTR, where TR is the repetition time and mz is the slice-encoding number.
Diagrams of (a) 2D and (b) 3D Cartesian k-space sampling schemes based on the golden ratio.
kyn values, calculated by using Eq (1), are distributed uniformly along the y-axis. However, the low-frequency signals are concentrated near the center of the k-space. Given that the low-frequency signals contain the main information of the image and that information can be used to trace tissue movement, the calculation of kyn can be improved as follows:
(2)
where sign(α) denotes the sign of α. Fig 2(A) and 2(B) show 64 successive 2D k-space sampling trajectories with a spatial resolution of 256 × 256 derived using Eqs (1) and (2), respectively, where the black solid lines represent the k-space sampling trajectories. The improved trajectories in Fig 2(B) are denser near the k-space center, which could be used to reconstruct the main information of the image.
Diagrams of 64 successive 2D k-space sampling trajectories based on (a) Eq (1) and (b) Eq (2).
Fig 3 shows the point spread function of the golden-angle radial trajectory and the proposed trajectory based on Eq (2).
Point spread function comparison between (a) the golden-angle radial trajectory and (b) the proposed trajectory.
DMRI reconstruction method based on CS
The CS technique can be used to precisely reconstruct the images in undersampling conditions [10]. Applying CS to MRI can significantly shorten the scanning time and provide a basis for achieving a high temporal resolution in DMRI. The fundamental concept behind DMRI image reconstruction based on CS is to solve the following optimization problem:
(3)
where x is the image, y is the measured k-space data, F is the Fourier transformation matrix, and Φ is a sparse matrix. Optimization of Eq (3) occurs when the vector z is the sparsest. Of the many algorithms that can solve Eq (3), we used a total-variation-based self-adaption algorithm [11], in which the following optimization problem must be solved:
(4)
where ∇ is the difference operator matrix, ||A||1 is the L1-norm, ||B||2 is the L2-norm, arg min{C} is the minimization of the argument C, and λ is the regularization factor.
Simulation experiments
MATLAB 7.0 (MathWorks, Natick, MA, USA) was used for 2D MRI simulation experiments in which a 2D “Modified Shepp-Logan” phantom image was used. Fig 4 shows the white ball, the movement of which simulates that of human organs or tissues. Two simulation experiments were designed to reflect the two main applications of DMRI. The periodic movement of human organs (such as cardiac motion) was imaged using three different trajectories for k-space data acquisition: a traditional Cartesian trajectory, a golden-angle-based radial trajectory, and the proposed sampling trajectory. The dynamic images were reconstructed using traditional Fourier transformation, regridding-based reconstruction, and CS-based reconstruction. The images reconstructed using these three methods were then compared.
Experiment I
In experiment I, the size of the phantom image was 256 × 256 pixels, the ball was located at the center of the phantom image, and its radius varied with time as a cosine function, with a period of 640TR. The radius r(t) is expressed as
(5)
where round(a) indicates that the expression a should be rounded. r(t) is in units of pixels and has an initial value of 12 and time t is in units of TR.
The three image reconstruction methods used in the simulation experiments are based on retrospective reconstruction [12]. Data acquired within N periods are assigned to the same phase interval according to the state of motion of the phantom and then used to reconstruct the corresponding phase image. Therefore, the undersampling factor can be written as
(6)
where tres is the temporal resolution and TR is the unit.
To compare the results of the three reconstruction methods, we reconstructed eight phase images that traced the states of motion of the ball, with tres = 10 and N = 1–25. In addition, the mean aliasing artifact power (AP) between the eight reconstructed images and the phantom images was calculated as follows:
(7)
where i and j are the pixel indexes of the 2D images and I and I' are phantom images and reconstructed images, respectively.
Experiment II
In the second experiment, the size of the phantom image was 256 × 256 pixels, the radius of the ball was 15 pixels, and the center of the ball moved linearly along the y-axis, where
(8)
where y(t) is in pixels and t is in units of TR.
To compare the results of the three reconstruction methods, we reconstructed images at eight different times that traced the states of motion of the ball using 20 different tres values in the range of 5–100. The mean AP values of the corresponding states of motion were also calculated.
Parameters of the reconstruction method
The parameters of the regridding-based reconstruction method included the Kaiser-Bessel convolution kernel [13] with a window width of 2, β = 18.5547, and an oversampling ratio of 2. For the self-adapting CS reconstruction method, we set λ = 0.05, the adaption coefficient to 0.6, the initial search step size to 1, and the maximum number of iterations to 100, and used a Wolfe line search [14].
Simulation results and discussion
Tissue and organ motion simulation results
Fig 5 shows the images reconstructed using the four methods, with tres = 10 and N = 7 (η ≈ 0.27). Fig 5(A)–5(E) show the phantom images, the images reconstructed from data acquired using the conventional Cartesian sampling method, the images reconstructed from data acquired using the golden-angle radial sampling method using the regridding reconstruction method, the images reconstructed from the data acquired using the golden-angle radial sampling method using the CS reconstruction method, and the images reconstructed using the proposed method, respectively. The images in Fig 5(B) show severe motion artifacts and the effects of undersampling. The images in Fig 5(C) have fewer motion artifacts than those in Fig 5(B); however, there are severe aliasing artifacts due to the effects of undersampling. The images for the golden-ratio angle trajectory sampling combined with the CS reconstruction method [Fig 5(D)] and the proposed method [Fig 5(E)] show better results because of the advantages of golden-ratio linear data acquisition and CS reconstruction. The results of the simulation show that both methods suppress motion artifacts and aliasing artifacts very well and yield better-quality images. However, in a clinical scan, radial trajectory sampling would suffer from the gradient delay effect caused by the MRI system, which would require additional phase correction steps [15]. The proposed method uses Cartesian k-space trajectory sampling, so no gradient delay-related phase correction is needed.
(a) phantom images, images reconstructed from (b) data from the conventional Cartesian sampling method, (c) data from golden-angle radial sampling using regridding reconstruction method, (d) data from golden-angle radial sampling using CS reconstruction method, and (e) the proposed method.
Fig 6 shows the mean AP values of the images reconstructed from the four methods, with tres = 10 and N = 1–25 (η ≈ 0.04–0.98): (a) the conventional Cartesian sampling method, (b) the golden-angle radial sampling method combined with the regridding image reconstruction method, (c) the golden-angle radial sampling method combined with the CS image reconstruction method, and (d) the proposed method. The AP values of the conventional Cartesian method are greater than 0.5 for N = 1–25, which indicates that the image quality of this method is unacceptable. The AP values of the golden-angle radial sampling method display a downward trend with increasing N. However, when η < 1, the AP value is still not ideal because the regridding image reconstruction method generally requires the acquired data to satisfy the Nyquist sampling theorem [16]. For this simulation experiment, at least 402 echoes must be acquired to meet this requirement. When tres = 10, N should be >40, which is oversampling. Fig 6(C) shows that the AP value drops to nearly zero when only five periods of data (N = 5) are used. Under the same temporal resolution condition, the proposed method can greatly reduce the scan time. Moreover, the proposed method provides the best image quality.
Drug delivery simulation results
Fig 7 shows the mean AP values of images reconstructed from the four methods but with different temporal resolutions. The AP value is higher for a lower temporal resolution for both traditional Cartesian and golden-angle radial acquisition methods. With the traditional method, there is a tradeoff between higher image quality and lower temporal resolution. However, lower temporal resolution causes more severe motion artifacts.
Fig 8 shows that when the temporal resolution is high [Fig 8(A)], the image quality is not ideal because of undersampling. The temporal resolution of the image in Fig 8(C) is the lowest of the three images, but the state of motion of the ball is the blurriest.
Reconstructed images from the golden-angle radial sampling method with temporal resolutions of (a) 10TR, (b) 30TR, and (c) 50TR.
The golden-angle radial sampling method combined with the CS image reconstruction method and the method proposed in this paper can provide high image quality under undersampling conditions. In addition, the image quality improves when the undersampling factor increases because the CS reconstruction algorithm allows high-quality reconstructed images from very little data. The proposed method and the golden-angle radial sampling method can provide approximately uniform k-space data within an arbitrary sampling window, which a conventional sampling trajectory cannot provide. Therefore, the proposed method can improve the temporal resolution.
Our results also showed that when tres was >30 [Fig 7(C)] and 40 [Fig 7(D)], the AP values increase because when selecting more data to reconstruct a frame image, movement of the organ or tissue introduces motion artifacts in the reconstructed image, thus reducing the image quality. Therefore, it is important to select the appropriate undersampling factor and set the temporal resolution appropriately.
Conclusion
Conventional DMRI techniques suffer from low temporal resolution and long scan time. Golden-angle radial trajectory sampling combined with a CS image reconstruction method provides better results but has a gradient delay effect, which imposes the need for additional steps to correct the phase. The proposed method improved the temporal resolution, shortened the scan time, and produced high-quality reconstructed images in the simulation experiments. In this study, we proposed and investigated DMRI that uses golden-ratio Cartesian trajectory sampling and a CS reconstruction algorithm. However, the proposed method has some problems. First, the state of motion of human organs or drug delivery in clinical scans is complex; however, the simulation experiments used in our study were based on only simple cosine and linear functions. Second, the design of the acquisition trajectories must be further optimized because linear acquisition is more sensitive to motion in the phase-encoding direction. Third, the CS reconstruction algorithm used in this study is a generic algorithm; the reconstruction algorithm needs improvement. Finally, the proposed method must be studied in more detail and validated after evaluation in a clinical situation.
Supporting information
S1 Matlab code. This is the major part of the simulation codes.
https://doi.org/10.1371/journal.pone.0191569.s001
(RAR)
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant No. 61671026) and the Natural Science Foundation of Beijing, China (grant No. 7162112).
References
- 1. Zhu Y, Liu J, Weinsaft J, Spincemaille P, Nguyen TD, Prince MR, et al. Free-breathing 3D imaging of right ventricular structure and function using respiratory and cardiac self-gated cine MRI. Biomed Res Int. 2015;2015:819102. pmid:26185764; PubMed Central PMCID: PMC4491385.
- 2. Liu DD, Liang D, Zhang N, Liu X, Zhang YT. Undersampling trajectory design for compressed sensing based dynamic contrast-enhanced magnetic resonance imaging. J Electron Imaging. 2015;24(1). Artn 013017 WOS:000350466100018.
- 3. Piludu F, Marzi S, Pace A, Villani V, Fabi A, Carapella CM, et al. Early biomarkers from dynamic contrast-enhanced magnetic resonance imaging to predict the response to antiangiogenic therapy in high-grade gliomas. Neuroradiology. 2015;57(12):1269–80. WOS:000365222500010. pmid:26364181
- 4. Piccini D, Feng L, Bonanno G, Coppo S, Yerly J, Lim RP, et al. Four-dimensional respiratory motion-resolved whole heart coronary MR angiography. Magn Reson Med. 2016. pmid:27052418.
- 5. Feng L, Axel L, Chandarana H, Block KT, Sodickson DK, Otazo R. XD-GRASP: Golden-angle radial MRI with reconstruction of extra motion-state dimensions using compressed sensing. Magnet Reson Med. 2016;75(2):775–88. WOS:000370597000030. pmid:25809847
- 6. Winkelmann S, Schaeffter T, Koehler T, Eggers H, Doessel O. An optimal radial profile order based on the golden ratio for time-resolved MRI. IEEE T Med Imaging. 2007;26(1):68–76. WOS:000243286800006. pmid:17243585
- 7. Crowe ME, Larson AC, Zhang Q, Carr J, White RD, Li DB, et al. Automated rectilinear self-gated cardiac cine imaging. Magnet Reson Med. 2004;52(4):782–8. WOS:000224238600012. pmid:15389958
- 8. Tolouee A, Alirezaie J, Babyn P. Compressed sensing reconstruction of cardiac cine MRI using golden angle spiral trajectories. J Magn Reson. 2015;260:10–9. WOS:000364982200002. pmid:26397216
- 9. Yang Y, Liu F, Jin ZY, Crozier S. Aliasing artefact suppression in compressed sensing MRI for random phase-encode undersampling. IEEE T Bio-Med Eng. 2015;62(9):2215–23. WOS:000360394900014. pmid:25850083
- 10. Donoho DL. Compressed sensing. IEEE T Inform Theory. 2006;52(4):1289–306. WOS:000236714000001.
- 11. Block KT, Uecker M, Frahm J. Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint. Magnet Reson Med. 2007;57(6):1086–98. WOS:000246979900013. pmid:17534903
- 12. Liu J, Spincemaille P, Codella NC, Nguyen TD, Prince MR, Wang Y. Respiratory and cardiac self-gated free-breathing cardiac CINE imaging with multiecho 3D hybrid radial SSFP acquisition. Magn Reson Med. 2010;63(5):1230–7. pmid:20432294; PubMed Central PMCID: PMC3071925.
- 13. Jackson JI, Meyer CH, Nishimura DG, Macovski A. Selection of a convolution function for fourier inversion using gridding. IEEE T Med Imaging. 1991;10(3):473–8. WOS:A1991GH74400029. pmid:18222850
- 14. Alhawarat A, Mamat M, Rivaie M, Salleh Z. An efficient hybrid conjugate gradient method with the strong Wolfe-Powell line ssearch. Math Probl Eng. 2015. Artn 103517 WOS:000358882500001.
- 15. Peters DC, Derbyshire JA, McVeigh ER. Centering the projection reconstruction trajectory: Reducing gradient delay errors. Magnet Reson Med. 2003;50(1):1–6. WOS:000183961800001 pmid:12815671
- 16. Zhu YC, Du J, Yang WC, Duan CJ, Wang HY, Gao S, et al. Reduced aliasing artifacts using shaking projection k-space sampling trajectory. Chinese Phys B. 2014;23(3). Artn 038702 WOS:000335646100102.