Elsevier

NeuroImage

Volume 167, 15 February 2018, Pages 276-283
NeuroImage

Rapid two-step dipole inversion for susceptibility mapping with sparsity priors

https://doi.org/10.1016/j.neuroimage.2017.11.018Get rights and content

Highlights

  • A rapid two-step dipole inversion algorithm is proposed.

  • 50 times faster when compared to MEDI and HEIDI in one healthy volunteer data set.

  • High-quality susceptibility maps from low resolution clinical data set.

Abstract

Quantitative susceptibility mapping (QSM) is a post-processing technique of gradient echo phase data that attempts to map the spatial distribution of local tissue magnetic susceptibilities. To obtain these maps, an ill-posed field-to-source inverse problem must be solved to remove non-local magnetic field perturbations. Current state-of-the-art algorithms which aim to solve the dipole inversion problem are plagued by the trade-off between reconstruction speed and accuracy. A two-step dipole inversion algorithm is proposed to bridge this gap. Our approach first addresses the well-conditioned k-space region, which is reconstructed using a Krylov subspace solver. Then the ill-conditioned k-space region is reconstructed by solving a constrained l1-minimization problem. The proposed pipeline does not incorporate a priori information, but utilizes sparsity constraints in the second step. We compared our method to well-established QSM algorithms with respect to COSMOS in in vivo volunteer datasets. Compared to MEDI and HEIDI the proposed algorithm produces susceptibility maps with a lower root-mean-square error and a higher coefficient of determination, with respect to COSMOS, while being 50 times faster. Our two-step dipole inversion algorithm without a priori information yields improved QSM reconstruction quality at reduced computation times compared to current state-of-the-art methods.

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: ϕ=F1DFχ, 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 invertingϕ=F1DFχwellusing 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:χ=argminχχTV+μ2Mχχwell22where M=F1RFis an operator that applies the binary maskR=|D|>δ, which sets the Fourier coefficients in the ill-conditioned regions to zero. μ is the regularization parameter, χTV=χ2is 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 y=χ and the dual variable p, we can write the augmented Lagrangian as:L(χ,y,p)=y1+μ2Mχχwell22pT(yχ)+α2yχ22where α is the regularization parameter for the quadratic penalty term. The ADMM technique allows to find a saddle point ofL(χ,y,p), 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:χi+1=argminχμ2Mχχwell22piT(yiχ)+α2yiχ22yi+1=argminyy1piT(yχi+1)+α2yχi+122pi+1=piα(yi+1χi+1)

Equation (3) can be solved by expanding the norms, differentiating with respect to χ, and setting the result to zero to get its normal equation:(μMTM+αT)χ=μMTχwell+T(αyp)

By taking advantage of the identity T=Δ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. Δ=F1LFwhere L is a diagonal matrix, and subsequently solve for χi+1 using only two DFTs:χ=F1(μMTMαL)1F(μMTχwell+T(αyp))

The solution to Equation (4) can be obtained by using a generalized shrinkage formula (Wang et al., 2007):y=shrink(χ+1αp,1α)=max(s1α,0)χ+1αpswheres=|xχ+1αpx|2|yχ+1αpy|2|zχ+1αpz|2

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:{γα,ifyi+1-χi+1λyi-χiα,otherwise

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)

  • W. Bian et al.

    A serial in vivo 7T magnetic resonance phase imaging study of white matter lesions in multiple sclerosis

    Mult. Scler. J.

    (2013)
  • B. Bilgic et al.

    Fast quantitative susceptibility mapping with L1-regularization and automatic parameter selection: fast QSM with L1-Regularization

    Magn. Reson. Med.

    (2014)
  • S.H. Chan et al.

    An augmented lagrangian method for total variation video restoration

    IEEE Trans. Image Process

    (2011)
  • A. Deistung et al.

    Overview of quantitative susceptibility mapping

    NMR Biomed.

    (2016)
  • C. Denk et al.

    Susceptibility weighted imaging with multiple echoes

    J. Magn. Reson. Imaging

    (2010)
  • A. Eissa et al.

    Detecting lesions in multiple sclerosis at 4.7 tesla using phase susceptibility-weighting and T2-weighting

    J. Magn. Reson. Imaging

    (2009)
  • D. Fong et al.

    LSMR: an iterative algorithm for sparse least-squares problems. SIAM j

    Sci. Comput.

    (2011)
  • T. Goldstein et al.

    Fast alternating direction optimization methods

    SIAM J. Imaging Sci.

    (2014)
  • X. He et al.

    Biophysical mechanisms of phase contrast in gradient echo MRI

    Proc. Natl. Acad. Sci.

    (2009)
  • Cited by (0)

    View full text