Application of Laplacian-based Methods to Multi-echo Phase Data for Accurate Susceptibility Mapping
Emma Biondetti1, David L. Thomas2, and Karin Shmueli1

1Department of Medical Physics and Biomedical Engineering, University College London, London, United Kingdom, 2Leonard Wolfson Experimental Neurology Centre, UCL Institute of Neurology, University College London, London, United Kingdom


In Susceptibility Mapping (SM) using multi-echo gradient-echo phase data, unwrapping and/or background-field removal is often performed using Laplacian-based methods. However, SM pipelines in the literature have applied these methods at different stages. Here, using simulated and acquired images, we compared the performance of three pipelines that apply Laplacian-based methods at different stages. We showed that Laplacian-based methods alter the linearity of the phase over time. We demonstrated that only a processing pipeline that takes this into account, i.e. by fitting the multi-echo data over time to correctly estimate a field map before applying Laplacian-based methods, gives accurate susceptibility values.


In Susceptibility Mapping (SM), Laplacian-based techniques, e.g. SHARP1,2, have been used to perform unwrapping and/or background-field removal of multi-echo gradient-echo phase data. The unwrapped gradient-echo phase at time $$$t$$$ and location $$$\bf r$$$ is $$\phi(t, {\bf r}) = \gamma \cdot \Delta B({\bf r}) \cdot t + \phi_0(0,{\bf r}) \quad (1),$$ where $$$\gamma$$$ is the proton gyromagnetic ratio, $$$\phi_0$$$ the $$$t = 0$$$ phase offset and $$$\Delta B$$$ the total field variation from local ($$$\Delta B_{loc}$$$) and background ($$$\Delta B_{bg}$$$) sources. In SM, multi-echo acquisitions are preferable because they allow fitting to Equation (1) to give $$$\phi_0$$$ and increase the accuracy of $$$\Delta B$$$ estimates3.

Several studies4-6 have used Laplacian unwrapping at each echo time (TE) as an initial step. Assuming $$$\phi_0 = 0$$$, a field map has then been calculated by scaling the Laplacian-unwrapped phase according to (1) at each TE and averaging the results. Similarly, other studies7,8 have used simultaneous Laplacian unwrapping and background-field removal at each TE followed by scaling according to (1), again assuming $$$\phi_0 = 0$$$. However, averaging the processed phase images assumes that Laplacian-based methods preserve the linear phase-time dependence (1). In contrast, others2,9 have fitted the unwrapped phase over TEs to calculate $$$\Delta B$$$ (and $$$\phi_0 = 0$$$) and have then removed $$$\Delta B_{bg}$$$ from the fitted $$$\Delta B$$$ using a Laplacian-based technique.


Here, we applied three processing pipelines2,4-9 to phase images of a numerical phantom and a healthy volunteer. We investigated the effect of using Laplacian-based methods at different stages of the SM pipeline on the phase-time linearity (1) and the accuracy of $$$\chi$$$ estimation.


Laplacian-based phase unwrapping (Lap-Unw) or background-field removal (Lap-Bg) were implemented with SHARP1,2 and a 3x3x3 Laplacian kernel2. For Lap-Unw, non-eroded brain mask (FSL BET10 for the volunteer) and threshold11 $$$\sigma = 10^{-10}$$$ were used. For Lap-Bg, the same brain mask with 2-voxel12 erosion and threshold2 $$$\sigma = 0.05$$$ were used.


Wrapped phase ($$$\phi_0 = 0$$$) was simulated at five echoes (TE1/ΔTE = 10/10 ms) from a ground-truth susceptibility distribution13,14 (Zubal phantom15, Figure 1), using a Fourier-based forward model16. Background-field free field maps were then calculated using three distinct pipelines:

- Avg-Unw: Lap-Unw1 on the phase at each TE; field map calculation:5 $$\Delta B = \frac{\sum_{echo = 1}^5 \phi_{echo} (TE_{echo}, {\bf r})}{\gamma \sum_{echo = 1}^5 TE_{echo}} \qquad (2)$$

then Lap-Bg1 on $$$\Delta B$$$.

- Avg-Bg: Lap-Bg2 on the phase at each TE; field map calculation:7 $$ \Delta B = \frac{1}{5} \sum_{echo = 1}^5 \frac{\phi(TE_{echo}, {\bf r})}{\gamma \cdot TE_{echo}} \qquad (3).$$

- Fit: linear fit of the simulated (non-wrapped) phase over TEs; Lap-Bg2 on the fitted $$$\Delta B$$$2,9.

$$$\chi$$$ maps were calculated (TKD17 with correction for underestimation2 and $$$\delta = 2/3$$$) to assess the effect of each pipeline on the $$$\chi$$$ values. $$$\chi$$$ mean, standard deviation (SD) and Root Mean Square Error (RMSE) were calculated in all the regions of interest (ROIs) shown in Figure 1.


3D gradient-echo brain images of a healthy volunteer were acquired on a Philips Achieva 3T scanner with a 32-channel head coil, 1-mm isotropic resolution, 7 echoes (TE1/ΔTE = 3.7/6.9 ms), TR = 50 ms, SENSE factor = 2 and flip angle = 10°.

The effect of Lap-Unw and Lap-Bg was tested. For comparison, the phase at each TE was also unwrapped with PRELUDE18. The mean and SD of the processed phase were calculated in three ROIs (Figure 3a) drawn on the fifth-echo magnitude image.

Results and Discussion

Laplacian-based processing altered the linearity of the phase-time relationship (1) in both the numerical phantom (Figure 2) and the volunteer (Figure 3). Such alterations of (1) were expected, because SHARP2 involves non-linear operations. These findings suggest that Laplacian unwrapping does not only unwrap the phase but also removes some $$$\Delta B_{bg}$$$ from $$$\Delta B$$$, even with a very small $$$\sigma$$$.

In the phantom, scaling and averaging the SHARP-processed phase caused inaccuracies in the estimated field maps (Figure 4) and, therefore, errors in $$$\chi_{Avg-Unw}$$$ and $$$\chi_{Avg-Bg}$$$ versus $$$\chi_{Fit}$$$ (Figure 5). Unlike $$$\chi_{Fit}$$$, $$$\chi_{Avg-Unw}$$$ and $$$\chi_{Avg-Bg}$$$ underestimated $$$\chi_{True}$$$ in the caudate nucleus and globus pallidus, whereas all mean $$$\chi$$$ estimates were similar in the thalamus and white matter. $$$\chi_{Fit}$$$ had the lowest SDs or the smallest RMSE values in all ROIs except the globus pallidus, in which, however, $$$\chi_{Fit}$$$ had the lowest RMSE as a percentage of $$$\chi_{Fit}$$$.


We demonstrated that Laplacian-based techniques alter the phase-time linearity (1). We also showed that Fit, therefore, gave the most accurate $$$\chi$$$ results, suggesting that Fit, or analogous pipelines13 that fit the phase over multiple echoes, should be used before applying Laplacian-based methods.


We would like to acknowledge our healthy volunteers and Dr Marios Yiannakas for his assistance with MRI scanning.

Emma Biondetti is supported by the Engineering and Physical Sciences Research Council (EP/M506448/1).

Dr David L. Thomas is supported by the UCL Leonard Wolfson Experimental Neurology Centre (PR/ylr/18575).


1. F. Schweser, A. Deistung, B. W. Lehr et al. Quantitative imaging of intrinsic magnetic tissue properties using MRI signal phase: An approach to in vivo brain iron metabolism? NeuroImage 2011; 54(4): 2789–2807.

2. F. Schweser, A. Deistung, K. Sommer et al. Toward online reconstruction of quantitative susceptibility maps: Superfast dipole inversion. Magn. Reson. Med. 2013; 69(6): 1582–1594.

3. G. Gilbert, G. Savard, C. Bard et al. Quantitative comparison between a multiecho sequence and a single-echo sequence for susceptibility-weighted phase imaging. Magn. Reson. Imaging 2012; 30(5): 722–730.

4. B. Bilgic, A. P. Fan, J. R. Polimeni et al. Fast quantitative susceptibility mapping with L1-regularization and automatic parameter selection. Magn. Reson. Med. 2014; 72(5): 1444–1459.

5. W. Li, C. Langkammer, Y.-H. Chou et al. Association between Increased Magnetic Susceptibility of Deep Gray Matter Nuclei and Decreased Motor Function in Healthy Adults. NeuroImage 2015; 105: 45–52.

6. H. Wei, R. Dibb, Y. Zhou et al. Streaking artifact reduction for quantitative susceptibility mapping of sources with large dynamic range. NMR Biomed. 2015; 28(10): 1294–1303.

7. D. Qiu, R. C. Brown, B. Sun et al. Abnormal Iron Levels in the Brain of Pediatric Sickle Cell Disease Patients: a Study using Quantitative Susceptibility Mapping (QSM). Proc. Intl. Soc. Mag. Reson. Med. 2014; 22: 0897.

8. R. C. Brown, D. Qiu, L. Hayes et al. Quantitation of Iron in Brain Structures of Children with Sickle Cell Disease and Transfusion Hemosiderosis. Blood 2014; 124(21): 1393.

9. H. Sun, A. J. Walsh, R. M. Lebel, et al. Validation of quantitative susceptibility mapping with Perls’ iron staining for subcortical gray matter. NeuroImage 2015; 105: 486–92.

10. S. M. Smith. Fast robust automated brain extraction. Hum. Brain Mapp. 2002; 17(3): 143–55.

11. MEDI toolbox, MATLAB function: unwrapLaplacian.m.

12. F. Schweser, K. Sommer, M. Atterbury et al. On the impact of regularization and kernel type on SHARP-corrected GRE phase images. Proc. Intl. Soc. Mag. Reson. Med. 2011; 19: 2667.

13. T. Liu, C. Wisnieff, M. Lou et al. Nonlinear formulation of the magnetic field to source relationship for robust quantitative susceptibility mapping. Magn. Reson. Med. 2013; 69(2): 467–476.

14. H. Sun and A. H. Wilman. Background field removal using spherical mean value filtering and Tikhonov regularization. Magn. Reson. Med. 2014; 71(3): 1151–1157.

15. I. G. Zubal, C. R. Harrell, E. O. Smith et al. Two dedicated software, voxel-based, anthropomorphic (torso and head) phantoms. Proceeding Int. Conf. Natl. Radiol. Prot. Board 19.

16. J. P. Marques and R. Bowtell. Application of a Fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility. Concepts Magn. Reson. Part B Magn. Reson. Eng. 2005; 25B(1): 65–78.

17. K. Shmueli, J. A. de Zwart, P. van Gelderen et al. Magnetic susceptibility mapping of brain tissue in vivo using MRI phase data. Magn. Reson. Med. 2009; 62(6): 1510–1522.

18. M. Jenkinson. Fast, automated, N-dimensional phase-unwrapping algorithm. Magn. Reson. Med. 2003; 49(1): 193–197.


Figure 1. Susceptibility phantom and $$$\chi$$$ values. $$$\chi$$$ values were taken from Ref. 13 for white matter, caudate nucleus, putamen, globus pallidus, thalamus, superior sagittal sinus and other brain regions and from Ref. 14 for air/non-brain regions.

Figure 2. Phantom: effect of Avg-Unw (b) and Avg-Bg (c) on the phase simulated from the ground-truth $$$\chi$$$ distribution (a). Mean and SD of the phase in four ROIs (caudate nucleus, thalamus, globus pallidus and white matter) are plotted against TEs.

Figure 3. Healthy volunteer: effect of PRELUDE unwrapping (b), Lap-Unw (c) and Lap-Bg (d) on the measured phase. Mean and SD of the phase in three ROIs drawn on the fifth-echo magnitude image (a) are plotted against TEs.

Figure 4. Phantom: transverse and coronal slices of field maps calculated with Avg-Unw (a), Avg-Bg (b) and Fit (c).

Figure 5. Phantom: transverse and coronal slices of $$$\chi$$$ maps calculated with Avg-Unw (b), Avg-Bg (c) and Fit (d) versus the true $$$\chi$$$ (a).

The table (e) shows the true and calculated $$$\chi$$$ [ppm] (mean ± SD) and the RMSEs [ppm] for each method in all ROIs.

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)