A General Algorithm for Compensation of Trajectory Errors: Application to Radial Imaging
Correspondence to :
L420 Pappajohn Biomedical Discovery Building
169 Newton Road
Iowa City, Iowa, 52242
phone number: (319) 335-5183.
Word count : about 2800
Figures+ tables count : 5
Running title: Annihilating filter based k-space shift correction
Purpose: To reconstruct artifact-free images from measured k-space data, when the actual k-space trajectory deviates from the nominal trajectory due to gradient imperfections.
Methods: Trajectory errors arising from eddy currents and gradient delays introduce phase inconsistencies in several fast scanning MR pulse sequences, resulting in image artifacts. The proposed algorithm provides a novel framework to compensate for this phase distortion. The algorithm relies on the construction of a multi-block Hankel matrix, where each block is constructed from k-space segments with the same phase distortion. In the presence of spatially smooth phase distortions between the segments, the complete block-Hankel matrix is known to be highly low-rank. Since each k-space segment is only acquiring part of the k-space data, the reconstruction of the phase compensated image from their partially parallel measurements is posed as a structured low-rank matrix optimization problem, assuming the coil sensitivities to be known.
Results: The proposed formulation is tested on radial acquisitions in several settings including partial Fourier and golden-angle acquisitions. The experiments demonstrate the ability of the algorithm to successfully remove the artifacts arising from the trajectory errors, without the need for trajectory or phase calibration. The quality of the reconstruction was comparable to corrections achieved using the Trajectory Auto-Corrected Image Reconstruction (TrACR) for radial acquisitions.
Conclusion: The proposed method provides a general framework for the recovery of artifact-free images from radial trajectories without the need for trajectory calibration.
Keywords: trajectory correction, annihilating filter, calibration-free, radial, EPI, MUSSELS, structured low rank
Several different k-space sampling strategies, each with unique benefits, are used in MRI. Echo-planar imaging (EPI) trajectories that scans the k-space in a single sweep, is routinely used in high temporal sampling applications such as functional MRI and diffusion MRI(1, 2, 3). Likewise, radial schemes, that sample the k-space center during each projection is ideally suited for accelerated imaging (4). Center-out radial methods with their short read-outs has enabled zero echo-time applications such as lung and sodium imaging (5, 6). Despite their benefits, a key problem that affects the performance of the above fast-scan methods is gradient fluctuations. Several hardware imperfections exist in the MR system that limit the exact realization of a prescribed gradient waveform on the MR hardware. Among these, eddy currents and gradient timing delays are understood as the leading source of gradient fluctuations (7, 8). The gradient fluctuations result in the actual k-space encoding to be different from the expected trajectory leading to reconstruction errors. The manifestation of the artifacts in the images will vary based upon the trajectory (9).
The classical approach to compensate for the trajectory errors is to measure the deviated k-space trajectory either using specialized hardwares (10) or employing calibration scans (11). The calibration-based methods have been proposed for radial (12) and EPI trajectory correction (13). However, these methods are not effective in dealing with temporal variations of the trajectory arising from gradient heating. Subject motion between the calibration scans and the imaging experiment can also lead to poor corrections. Self-calibrating trajectory correction methods are immune to these problems.
Self-calibrating methods based on iterative phase cycling were originally introduced for EPI (14, 15, 16) and later extended to radial acquisitions (17, 18). The EPI phase cycling methods reconstruct a series of images using different choices of phase errors between the odd and even lines and the phase which minimized some criteria was used in the final reconstruction. Likewise, Deshmane et al. (17) iteratively shift the radial data until it maximizes the root sum-of-squares (SOS) DC signal of each individual projection, while Wech et al. (18) iteratively shift the radial trajectory in a direction that provides the best data fit. Recently, these methods were generalized by an optimization scheme for multiple trajectories; this approach relies on a basis expansion for eddy current and gradient non-linearities, along with a gradient descent optimization of the trajectory parameters (10). A challenge with the above approach is the non-quadratic nature of the optimization scheme, associated with the need for careful initialization of the parameters to ensure convergence.
Recently, we proposed a self-calibrating Nyquist ghost correction strategy for EPI reconstruction (19, 20). Here, we recovered the missing k-space samples in the odd and the even blocks after the odd and the even lines of the EPI k-space were separated from each other. The images corresponding to the odd/even datasets differ only by a phase error, induced by the k-space shifts. We observed that these phase relations manifest as annihilation conditions which enabled the formulation of a block-Hankel matrix, created from these datasets, that is heavily low-rank. We rely on Hankel structured low-rank matrix completion (SLRMC) to fill in the missing entries in the datasets. This phase-compensating reconstruction is termed as MUSSELS(21, 19). In the current work, we generalize the above idea to trajectories beyond EPI. We show the applicability of the method for correction of trajectory errors in multiple radial-based acquisitions where we do not know the trajectory shifts a priori. A preliminary account of this work was presented at the 2017 ISMRM (22).
Modeling of trajectory errors as phase modulations
Gradient fluctuations arising from hardware imperfections result in trajectory shifts in several acquisitions. The proposed method is based on the idea that, in many cases, the effect of k-space shifts can be modeled as phase modulations in the data. Let us denote the ideal image that is not affected by trajectory errors and its k-space data by and respectively for an arbitrary trajectory; and being the spatial and k-space position index vectors. In the presence of trajectory errors, the k-space lines of will be shifted by various amounts. We propose that the shifted k-space data can be segmented into segments, where the data in each segment has the same shift.
Here, ; denote the discrete shifts of each segment. An illustration with (N=) 4 different shifts affecting the k-space lines of is provided in Figure 1a and our proposed splitting is illustrated in Figure 1b. Note that each of the segments are not fully sampled in k-space. If they were fully sampled, we can recover the images by simply taking the inverse Fourier transform of the full k-space data as:
Note that the images are phase-modulated versions of the underlying image . Since each of the datasets are not fully sampled, we denote the measurement process as
Here, includes the coil sensitivity weighting, Fourier transform, as well as the sampling operator that selects the acquired samples for the segment. Similarly, denotes the acquired measurements of the dataset corresponding to the segment. We propose to recover the images using the phase relations between them.
Trajectory Corrected Reconstruction based on Structured Low-Rank Matrix Completion
In (21), we showed that there exists annihilation relations between the Fourier samples of the phase-modulated images . Here, * denotes convolution and and are 2-D finite impulse response filters. These annihilation conditions can be compactly expressed as null-space conditions on the below block-Hankel matrix:
Here, each block is a Hankel matrix, formed out of the fully sampled k-space data of the phase-modulated images (23). is the operator that converts the k-space data of the phase-modulated images into a multi-block Hankel matrix (please see supporting information for details). The null-space conditions induced by the annihilation relations imply that is a low-rank matrix.
Since each of the phase-modulated images are not fully sampled, several entries of the matrix are unknown. These entries can be recovered using the SLRMC approach as in (23, 24, 25, 26). However, in contrast to the parallel imaging setting in (23, 24, 25, 26), there are no shared samples of k-space between (Figure 1b). Hence, the direct extension of the approaches in (23, 24, 25, 26) to our setting will be ill-posed. Specifically, one could fill the missing entries in each matrix with the samples from other ones. This will also result in a heavily low-rank matrix since the k-space data for all blocks are the same (27). To avoid this pitfall, the MUSSELS algorithm is formulated as:
The reconstruction given in Eq. 5 with the forward model operator appropriately defined for the specific trajectory gives the general trajectory correction method based on the idea of phase compensation.
The key step in achieving effective trajectory correction based on the above method is grouping the k-space lines that experience similar errors. While this step is straightforward for some trajectories (e.g. EPI), a well-defined separation may not be easily available for several trajectories, as in the case of conventional or golden angle radial acquisitions. The challenge is in identifying this grouping. For both conventional and golden angle radial trajectories, it has been shown that the trajectory shifts vary as a function of the angle of the spoke (28, 10). Therefore, we propose a grouping that is based on the angle of the radial spokes for these trajectories. Specifically, the radial spokes with similar angles can be assumed to have similar trajectory errors. This enables us to generalize the phase-compensated reconstruction based on SLRMC for correction of trajectory errors in radial-based acquisitions also.
Denoting the measured radial k-space data on the non-Cartesain grid by , we split the radial spokes of a given frame into segments based on their angles. To keep consistent notations for Cartesian and non-Cartesian trajectories, we denote the images and the k-space data on the Cartesian lattice corresponding to the N segments of the radial trajectories as and , respectively. The operator, , for the case of radial trajectories, is given by where represent coil sensitivities and represent the non-uniform fast Fourier transform (NUFFT) to the grid corresponding to . is computed using the Fourier samples of on the Cartesian grid.
The proposed reconstruction in Eq. 5 was implemented in MATLAB 2016a (The Mathworks, Natick, MA) on a desktop PC with an Intel i7-4770, 3.4 GHz CPU with 8 GB RAM. The different steps involved in the proposed reconstruction are illustrated in Figure 2, which is solved using an augmented Lagrangian scheme (derived in the supporting information). The NUFFT is implemented as in (29). The conjugate gradient updates use 20 iterations or a tolerance of 1e-15, with no pre-conditioning.
Simulation study: We first demonstrate the trajectory correction employing MUSSELS (referred to as T-MUSSELS) using simulations. We used a 32-channel brain data (available as part of the ESPIRiT toolbox) for this experiment. We simulated trajectory errors into a nominal radial trajectory by explicitly shifting the spokes of the trajectory. The shifts were generated as a function of the spoke angle. The k-space data to be used for testing was synthesized by performing a NUFFT of the above brain data using the corrupted trajectory. Using this simulated experiment, we also demonstrate that it is possible to recover the corrected image for a given trajectory error using more than one scheme of segmentation for the case of radial trajectories.
Experimental datasets: Trajectory correction using T-MUSSELS was also validated using phantoms and in-vivo data from cardiac and neuroimaging studies. To demonstrate the generalizability of the above method, the experimental data was collected using a variety of different radial acquisitions. The cardiac data was acquired using a short-axis radial scan on a Siemens 3T scanner with 3-channel acquisition. The golden angle (GA) radial brain data is distributed as part of the TrACR package (10). This 14-channel data consist of 201 radial spokes with 256 points per spoke. The phantom data was collected with a partial Fourier radial acquisition on a GE 3T scanner. This 32-channel data consist of 256 radial spokes with 144 points per spoke.
Coil sensitivity maps are required for all the reconstructions proposed here. The coil sensitivity maps were estimated from the data itself in all experiments, by taking the data from a small region around the center of k-space using the ESPIRiT method (30).
Simulations: The results of the simulated experiment are provided in Figure 3. The shifted radial trajectory is shown in Figure 3a. A CG-NUFFT reconstruction (29) of the shifted k-space data with the nominal radial trajectory show severe artifacts (Figure 3c). We used T-MUSSELS to correct for the trajectory errors using the two schemes shown in Figure 3b: by splitting the radial data into (i) 8 segments, and (ii) 6-segments. The resulting reconstructions are shown in Figure 3c. As noted from the results, it is possible to achieve good reconstruction with more than one splitting scheme for the case of the radial trajectory using the proposed method.
In-vivo results: Figure 4 shows the reconstruction from two frames of a cardiac data that was acquired using a radial trajectory. Figure 4a show respectively the nominal radial trajectory and the trajectory that was split into 6 segments to be used with the proposed reconstruction. The left column in Figure 4b-c shows the CG-NUFFT reconstruction using the nominal trajectory from two cardiac frames. A 4x4 filter was used for the 6-block Hankel matrix construction. The trajectory corrected reconstructions using T-MUSSELS from this three-channel acquisition show reduction in streaking and improved reconstruction in many regions. For comparison, the correction achieved using the TrACR method is also provided. Notably, the T-MUSSELS reconstruction provide improved contrast in the heart region (boxed region).
Figure 5a shows the reconstruction of a GA radial data. Since the consecutive spokes of a GA radial are not close in angle, we reordered the radial data based on the angle of the spokes so that the data obtained from spokes with similar angles are grouped together. The new grouping, which segmented the spokes into 8 segments are shown in figure 5a(ii). Figure 5a(iii) shows the CG-NUFFT reconstruction that assumes the nominal trajectory shown in Figure 5a(i). Figure 5a(iv)-(v) show the reconstructions using the TrACR and T-MUSSELS respectively, which show improved reconstructions.
Figure 5b shows T-MUSSELS reconstruction validated on a phantom data that was acquired using a partial Fourier radial trajectory. Here, the spokes were split into 4 segments based on their angles as shown in Figure 5b(ii). Figure 5b(iii) shows the CG-NUFFT reconstruction that assumes the nominal trajectory given in Figure 5b(i). The effect of trajectory errors are visible in this reconstruction. Figure 5b(iv) shows the trajectory corrected reconstructions using the TrACR method. Even though there is significant reduction in artifacts, residual errors can still be seen in the resulting reconstruction. Figure 5b(v) shows the proposed reconstruction using the nominal trajectory which shows improved results compared to the TrACR method. To quantify the difference, we computed the background-to-object signal intensity ratio similar to the ghost-to-signal ratio (GSR) that is typically computed for EPI based reconstructions (26). As observed, the GSR is much lower in T-MUSSELS compared to TrACR.
The proposed method provides a general self-calibrating two-dimensional phase compensation framework for the correction of trajectory errors. This work adds to the recent literature on structured low-rank based calibration-free strategies, which are emerging as powerful tools in a variety of applications including parallel MRI (31, 23, 32, 24, 25, 33, 26) and correction of artifacts in MRI (25, 21, 26, 19). As far as we know, this is the first work that use the idea of SLRMC for radial trajectory correction. The robustness of the reconstruction was tested using simulations and in-vivo data and compared to state-of-the-art methods. The proposed method can be extended to 3D and 4D radial trajectory corrections also. However, the high memory demands of the present method would require the adaption of memory-efficient structured low rank methods such as (34, 35) for efficient implementation.
Due to their data-driven nature, the self-calibrating methods are more robust to temporal variations in trajectory errors compared to calibration-based correction methods. There are two potential ways in which these methods can be used for dynamic scan corrections. In the first approach, the trajectory correction can be performed on the first time point and the null space learned from this can be applied to the other time points of the dynamic scan. This approach may be prone to errors if the null space learned from the first time point may not be valid for a later time point; however it can be re-learned from a closer scan, similar to the idea of interleaved calibration scans for dynamics scanning (12). A second approach is to perform independent trajectory correction for each time point at the expense of increased computation time.
Another advantage of the proposed scheme as a general phase compensation strategy is that this formulation allows for compensation of phase errors regardless of their source. Hence parametrization of the phase errors is not required in our approach. Regardless, the proposed method is shown to work well for in-vivo data, where we do not know the characteristics of the phase errors present in the data beforehand. Our only assumption is that radial spokes with similar angles exhibit similar errors. While the partial Fourier data provided good correction with a 4-segment split, the conventional radial and golden angle full Fourier data provided good correction with a 6-segment and an 8-segment split respectively. Based on simulation studies, we observe that it is possible to get good reconstructions with more than one scheme of segmentation. It is intuitive to see that if the number of segments chosen are a multiple of the number of phases simulated, the proposed phase compensation constraint should still hold. Thus, as long as the trajectory shifts are within the band-limits of the implicit filter (parametrized by the size of the window chosen to form the Hankel matrix), the phase compensation can be effectively achieved. For small trajectory shifts this should be valid and we might be able to find band-limited annihilating filters for several choices of segments.
Other structured low-rank based self-calibrating ghost correction methods have been proposed for EPI trajectories (26, 27). In addition to the low-rank constraint, (26) makes use of a sparsity constraint to recover ghost-corrected EPI data on a channel-by-channel basis. In contrast, our formulation imposes a parallel-imaging based data consistency to recover images. While this makes the proposed method applicable to only multi-coil acquisitions, the cardiac data in Figure 4 use only 3 channels in the reconstruction. Hence the number of coils is not critical for good performance of the algorithm. At the same time, this makes it possible to work with accelerated acquisitions also efficiently. (27) has shown the utility of non-convex penalties for improved performance for EPI correction.
In conclusion, we proposed a simple and generalizable method that formulates trajectory corrected reconstruction as a phase compensation problem. The method does not require explicit calibration of trajectory shifts or phase calibration and works for a variety of trajectories and accelerated acquisitions.
- Poustchi-Amin et al. (2013) Poustchi-Amin M, Mirowitz SA, Brown JJ, McKinstry RC, and Li T. Principles and applications of echo-planar imaging: a review for the general radiologist. Radiographics, 21(3):767–779, 2013.
- Poser and Norris (2007) Poser BA and Norris DG. Fast spin echo sequences for BOLD functional MRI. Magma, 20(1):11–7, 2007.
- Jaermann et al. (2004) Jaermann T, Crelier G, Pruessmann KP, Golay X, Netsch T, van Muiswinkel AMC, Mori S, van Zijl PCM, Valavanis A, Kollias S, and Boesiger P. SENSE-DTI at 3 T. Magnetic Resonance in Medicine, 51(2):230–6, 2004.
- Peters et al. (2006) Peters DC, Rohatgi P, Botnar RM, Yeon SB, Kissinger KV, and Manning WJ. Characterizing radial undersampling artifacts for cardiac applications. Magnetic Resonance in Medicine, 55(2):396–403, 2006.
- Nagel et al. (2009) Nagel AM, Laun FB, Weber MA, Matthies C, Semmler W, and Schad LR. Sodium MRI using a density-adapted 3D radial acquisition technique. Magnetic Resonance in Medicine, 62(6):1565–1573, 2009.
- Johnson et al. (2013) Johnson KM, Fain SB, Schiebler ML, and Nagle S. Optimized 3D ultrashort echo time pulmonary MRI. Magnetic Resonance in Medicine, 70(5):1241–1250, 2013.
- Boesch et al. (1991) Boesch C, Gruetter R, and Martin E. Temporal and spatial analysis of fields generated by eddy currents in superconducting magnets: Optimization of corrections and quantitative characterization of magnet/gradient systems. Magnetic Resonance in Medicine, 20:268–284, 1991.
- Aldefeld and Bornert (1998) Aldefeld B and Bornert P. Effects of gradient anisotropy in MRI. Magnetic Resonance in Medicine, 39:606–614, 1998.
- Smith and Nayak (2010) Smith TB and Nayak KS. MRI artifacts and correction strategies. Imaging in Medicine, 2(4):445–457, 2010.
- Ianni and Grissom (2016) Ianni JD and Grissom WA. Trajectory Auto-Corrected image reconstruction. Magnetic Resonance in Medicine, 76(3):757–768, 2016.
- Xu et al. (2010) Xu D, King KF, Zur Y, and Hinks RS. Robust 2D phase correction for echo planar imaging under a tight field-of-view. Magnetic Resonance in Medicine, 64(6):1800–13, 2010.
- Block and Uecker (2011) Block KT and Uecker M. Simple Method for Adaptive Gradient-Delay Compensation in Radial MRI. In Proceedings of the 19th Annual Meeting of ISMRM, 2011.
- Bruder et al. (1992) Bruder H, Fischer H, Reinfelder HE, and Schmitt F. Image reconstruction for echo planar imaging with nonequidistant k-space sampling. Magnetic Resonance in Medicine, 23:311–323, 1992.
- Clare (2003) Clare S. Iterative Nyquist Ghost Correction for Single and Multi-shot EPI using an Entropy Measure. In Proceedings of the 11th Annual Meeting of ISMRM, page 1411, 2003.
- Foxall et al. (1999) Foxall DL, Harvey PR, and Huang J. Rapid iterative reconstruction for echo planar imaging. Magnetic Resonance in Medicine, 42(3):541–7, 1999.
- Chen et al. (2011) Chen NK, Avram AV, and Song AW. Two-dimensional phase cycled reconstruction for inherent correction of echo-planar imaging nyquist artifacts. Magnetic Resonance in Medicine, 66(4):1057–66, 2011.
- Deshmane et al. (2016) Deshmane A, Blaimer M, Breuer F, Jakob P, Duerk J, Seiberlich N, and Griswold M. Self-calibrated trajectory estimation and signal correction method for robust radial imaging using GRAPPA operator gridding. Magnetic Resonance in Medicine, 75(2):883–896, 2016.
- Wech et al. (2015) Wech T, Tran-Gia J, Bley TA, and Köstler H. Using self-consistency for an iterative trajectory adjustment (SCITA). Magnetic Resonance in Medicine, 73(3):1151–7, 2015.
- Mani et al. (2016a) Mani M, Magnotta V, Kelley D, and Jacob M. Comprehensive reconstruction of multi-shot multi-channel diffusion data using MUSSELS. In Proc of Engineering in Medicine and Biology Society, Orlando, 2016a.
- Mani et al. (2017a) Mani M, Jacob M, Yang B, and Magnotta V. Comprehensive correction of Motion and Nyquist Ghost Artifacts for Multi-shot Diffusion Imaging. In Proceedings of the 25th Annual Meeting of ISMRM, 2017a.
- Mani et al. (2016b) Mani M, Jacob M, Kelley D, and Magnotta V. Multi-shot sensitivity-encoded diffusion data recovery using structured low-rank matrix completion (MUSSELS). Magnetic Resonance in Medicine, 78:494–507, 2016b.
- Mani et al. (2017b) Mani M, Poddar S, Magnotta V, and Jacob M. Trajectory Correction of Radial Data Using MUSSELS. In Proceedings of the 25th Annual Meeting of ISMRM, 2017b.
- Shin et al. (2014) Shin PJ, Larson PEZ, Ohliger MA, Elad M, Pauly JM, Vigneron DB, and Lustig M. Calibrationless parallel imaging reconstruction based on structured low-rank matrix completion. Magnetic Resonance in Medicine, 72(4):959–70, 2014.
- Haldar and Zhuo (2015) Haldar JP and Zhuo J. P-LORAKS: Low-rank modeling of local k-space neighborhoods with parallel imaging data. Magnetic Resonance in Medicine, 2015.
- Jin et al. (2017) Jin KH, Lee J, Lee D, and Ye JC. MRI Artifact Correction Using Sparse+Low-Rank Decomposition of Annihilating Filter-Based Hankel Matrix. Magnetic Resonance in Medicine, 78:327–340, 2017.
- Lee et al. (2016) Lee J, Jin KH, and Ye JC. Reference-free single-pass EPI Nyquist ghost correction using annihilating filter-based low rank Hankel matrix (ALOHA). Magnetic resonance in medicine, 2016.
- Lobos et al. (2017) Lobos RA, Kim TH, Hoge WS, and Haldar JP. Navigator-free EPI ghost correction using low-rank matrix modeling : Theoretical insights and practical improvements. In Proceedings of the 25th Annual Meeting of ISMRM, 2017.
- Peters et al. (2003) Peters DC, Derbyshire JA, and McVeigh ER. Centering the projection reconstruction trajectory: Reducing gradient delay errors. Magnetic Resonance in Medicine, 50:1–6, 2003.
- Yang and Jacob (2009) Yang Z and Jacob M. Efficient NUFFT algorithm for non-Cartesian MRI reconstruction. In 2009 IEEE International Symposium on Biomedical Imaging, pages 117–120, 2009. ISBN 978-1-4244-3931-7.
- Uecker et al. (2014) Uecker M, Lai P, Murphy MJ, Virtue P, Elad M, Pauly JM, Vasanawala SS, and Lustig M. ESPIRiT - An eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA. Magnetic Resonance in Medicine, 71(3):990–1001, 2014.
- Morrison et al. (2007) Morrison R, Jacob M, and Do M. Multichannel Estimation Of Coil Sensitivities In Parallel MRI. In 2007 IEEE International Symposium on Biomedical Imaging, pages 117–120. IEEE, 2007. ISBN 1-4244-0671-4.
- Haldar (2014) Haldar JP. Low-rank modeling of local k-space neighborhoods (LORAKS) for constrained MRI. IEEE Transactions on Medical Imaging, 33(3):668–81, 2014.
- Ye et al. (2015) Ye JC, Kim JM, Jin KH, and Lee K. Compressive Sampling using Annihilating Filter-based Low-Rank Interpolation. IEEE Transactions on Information Theory, 63(2):777 – 801, 2015.
- Ongie and Jacob (2017) Ongie G and Jacob M. A Fast Algorithm for Convolutional Structured Low-rank Matrix Recovery. IEEE Transactions on Computational Imaging, (99), 2017.
- Ongie and Jacob (2016) Ongie G and Jacob M. A fast algorithm for structured low-rank matrix recovery with applications to undersampled MRI reconstruction. In 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI), pages 522–525. IEEE, 2016. ISBN 978-1-4799-2349-6.