Rapid two-step dipole inversion for susceptibility mapping with sparsity priors
Introduction
Quantitative susceptibility mapping (QSM) is a novel reconstruction method that uses gradient echo (GRE) phase data to estimate the underlying tissue magnetic susceptibility (Deistung et al., 2016, Haacke et al., 2015, Li, 2001, Salomir et al., 2003, Schweser et al., 2016, Shmueli et al., 2009, Wang and Liu, 2015, Wharton et al., 2010). The phase contrast is related to the field perturbations caused by the tissues' susceptibility distributions, and is modulated by acquisition parameters, e.g. echo time (Wharton and Bowtell, 2012), and the tissues' micro- and macroscopic architectures (He and Yablonskiy, 2009, Shmueli et al., 2011, Zhong et al., 2008). Extracting the magnetic susceptibility from the GRE phase within each voxel is further complicated by the non-local relationship between the phase and the underlying magnetic susceptibility. In order to obtain tissue susceptibility maps, the field-to-source inversion has to be solved once the additional pre-processing steps of phase unwrapping and local field estimations have been performed. This inverse problem, known as dipole inversion or dipole deconvolution, is challenged by regions in the proximity of the conical surfaces of the dipole in k-space, where the Fourier coefficients decay to zero. The ill-conditioning can be mitigated by sampling multiple orientations, as performed by COSMOS (calculation of susceptibility through multiple orientation sampling) (Liu et al., 2009), the current gold-standard of susceptibility mapping. However, multiple orientation sampling is uncomfortable for the subjects, more time-consuming than single orientation scans, and may require the use of larger coils, often resulting in inferior signal-to-noise ratios.
To overcome these challenges, a number of groups have developed single-orientation regularization algorithms addressing the ill-posed dipole deconvolution (Bilgic et al., 2014, Liu et al., 2012, Liu et al., 2011, Schweser et al., 2012, Wang and Liu, 2015, Wu et al., 2012). Some approaches make use of priors, e.g. information obtained from the GRE-magnitude used in MEDI (Morphology Enabled Dipole Inversion) (Liu et al., 2012, Liu et al., 2011) or a combination of binary masks obtained from the gradient and laplacian of the phase data and the gradient of the magnitude data as used in HEIDI (Homogeneity Enabled Incremental Dipole Inversion) (Schweser et al., 2012) while other approaches utilize reconstruction strategies used in compressed sensing (CSC) (Wu et al., 2012). Deconvolution approaches that are k-space based are computationally cheap and offer fast reconstruction times at the cost of streaking artifacts and underestimation of the susceptibility values (Schweser et al., 2013, Wang and Liu, 2015, Wharton et al., 2010). Iterative image-space optimization methods, on the other hand, either produce high quality, accurate susceptibility maps in phantoms at long reconstruction times, (15 min to 2–3 h, e.g. MEDI, HEIDI, CSC) (Liu et al., 2012, Liu et al., 2011, Schweser et al., 2012, Wu et al., 2012), or provide fast reconstruction times (<5 min, e.g. total variation Split-Bregman (TVSB)) (Bilgic et al., 2014), but cause over-smoothing that could render visual inspection of small structures infeasible. Generally, in vivo data are reconstructed with lower accuracy compared to COSMOS than techniques achieve in phantom experiments (Wang and Liu, 2015).
Herein, we propose a rapid two-step dipole inversion algorithm, which adopts the incremental inversion scheme previously presented by Wu et al. (2012), Li et al. (Li et al., 2011, Li et al., 2015), and Schweser et al. (2012). without the use of a priori information to guide the inversion.
Assuming an isotropic susceptibility model, the underlying tissue magnetic susceptibility, χ, can be calculated from the measured local field perturbations, φ, by inverting: , where F represents the Fourier Transform, and D is the dipole kernel in k-space. Because D decays to zero in the proximity of two conical surfaces (Marques and Bowtell, 2005), the inversion becomes problematic in the presence of noise and non-harmonic non-local field contributions. In order to address the ill-conditioned k-space regions during the inversion, we adopt an incremental inversion strategy dividing k-space into a well-conditioned and an ill-conditioned region. Although such sub-domain algorithms have been previously applied to the deconvolution problem (Li et al., 2011, Schweser et al., 2012, Wu et al., 2012), each method utilizes different approaches for recovering the susceptibility information in each of the domains. In the 2016 QSM reconstruction challenge the proposed algorithm achieved the lowest root-mean-square error (RMSE) with respect to the provided susceptibility tensor imaging (STI) reconstruction across 27 submissions (Langkammer et al., 2017). Our two-step method is detailed below.
The well-conditioned k-space region is reconstructed by invertingusing an LSMR solver (Fong and Saunders, 2011) and exiting the iterations early, hence utilizing the implicit regularization properties of Krylov subspace methods (Wang et al., 2013, Wang and Liu, 2015), in order to avoid streaking artifacts.
Reconstructing the ill-conditioned k-space region is achieved by solving the following unconstrained TV-model:where is an operator that applies the binary mask, which sets the Fourier coefficients in the ill-conditioned regions to zero. μ is the regularization parameter, is the isotropic TV norm and is the discrete gradient operator in the three spatial directions. This model is equivalent to a CS-minimization, where the total variation is taken as the sparsifying transform (Lustig et al., 2007). CS-minimization has been previously applied to the dipole deconvolution problem by Wu et al. (2012) using a wavelet transform as well as a TV penalty. Incorporating two sparsifying transforms takes advantage of the sparsity of the image in both transform domains, but comes at the cost of increased reconstruction times and additional parameters which need to be optimized and may affect the solution. A TV approach similar to Equation (1) is also used in HEIDI for the reconstruction of the ill-conditioned k-space region. In contrast to HEIDI's TV, our TV approach does not include a priori information as weights.
Here, the minimization was solved using the alternating direction method of multipliers (ADMM), which is known as the Split-Bregman method when applied to l1-regularization problems (Goldstein et al., 2014). This method is also used in TVSB with the difference being that we apply the regularization parameter μ to the data fidelity term instead of the regularization term.
Introducing the constraint and the dual variable p, we can write the augmented Lagrangian as:where α is the regularization parameter for the quadratic penalty term. The ADMM technique allows to find a saddle point of, which is also the solution of our original problem described in Equation (1), by iteratively performing alternating minimizations of Equation (2) over χ and y. These ADMM iterations are given by:
Equation (3) can be solved by expanding the norms, differentiating with respect to χ, and setting the result to zero to get its normal equation:
By taking advantage of the identity and noting that under periodic boundary conditions the resulting matrix is block circulant, we can diagonalize the matrix by a discrete Fourier transform (DFT), i.e. where L is a diagonal matrix, and subsequently solve for χi+1 using only two DFTs:
The solution to Equation (4) can be obtained by using a generalized shrinkage formula (Wang et al., 2007):where
We further adopt an update scheme for the regularization parameter α that has shown to improve convergence rates in ADMM methods (Chan et al., 2011). Introducing the variables 0 < λ < 1 and γ > 1, the following update scheme ensures that the residue of the quadratic penalty term decreases asymptotically:
Consequently, the complete algorithm for the reconstruction of the ill-conditioned k-space region is given by:
Section snippets
Methods
The performance of the proposed algorithm was evaluated versus COSMOS using a single orientation from a multiple orientation COSMOS dataset acquired in one healthy volunteer. In a second volunteer dataset, different parallel imaging acceleration factors and image resolutions were tested in contrast to COSMOS. Both volunteer dataset were acquired at our centre. This study was approved by the Clinical Research Ethics Board of the University of British Columbia and is in accord with the
Parameter selection
Fig. 1 demonstrates the effect that the iteration number, chosen for exiting Step 1, has on the obtained susceptibility maps and the determined error values, R2 and RMSE, respectively. The displayed values for all δ′s are averages across the range of μ′s tested with error bars representing one standard deviation. We compared the effect on the full brain and the deep grey matter ROI separately in two datasets. Dataset A is the online data provided by the lab at Cornell University, dataset B
Discussion
We presented a rapid two-step dipole inversion algorithm for the calculation of susceptibility maps. By limiting the LSMR solver in Step 1 to a few iterations, the implicit regularization property of Krylov subspace methods ensures that reconstruction errors and noise in the ill-conditioned region do not result in streaking artifacts in image space. Although other Krylov subspace methods, e.g. CG or LSQR, in lieu of LSMR, produce similar error values and susceptibility maps, it was found
Conclusion
This study presented a two-step dipole inversion algorithm with reconstruction times 50 times faster than current state-of-the-art methods. When compared to COSMOS, the proposed method has shown to produce comparable in vivo susceptibility maps with a lower RMSE and R2 than the susceptibility maps produced by existing methods.
Acknowledgements
The authors would like to thank Yi Wang and Tian Liu for generously sharing their dataset. This study was supported by NSERC (RGPIN-2016-05371). VW is supported by a graduate studentship award from the MS Society of Canada (EGID 2002). AR is supported by Canada Research Chairs (230363). The UBC MRI Research Centre thanks Philips Medical for continuing support.
References (39)
- et al.
Quantitative susceptibility mapping: current status and future directions
Magn. Reson. Imaging
(2015) - et al.
Improved optimization for the robust and accurate linear registration and motion correction of brain images
NeuroImage
(2002) - et al.
A global optimisation method for robust affine registration of brain images
Med. Image Anal.
(2001) - et al.
A method for estimating and removing streaking artifacts in quantitative susceptibility mapping
NeuroImage
(2015) - et al.
Magnetic susceptibility anisotropy of human brain in vivo and its molecular underpinnings
NeuroImage
(2012) - et al.
Quantitative susceptibility mapping of human brain reflects spatial variation in tissue composition
NeuroImage
(2011) - et al.
Morphology enabled dipole inversion for quantitative susceptibility mapping using structural consistency between the magnitude image and the susceptibility map
NeuroImage
(2012) - et al.
Foundations of mri phase imaging and processing for quantitative susceptibility mapping (QSM)
Z. Für Med. Phys.
(2016) - et al.
Quantitative susceptibility mapping for investigating subtle susceptibility variations in the human brain
NeuroImage
(2012) - et al.
The molecular basis for gray and white matter contrast in phase imaging
NeuroImage
(2008)