Application of Laplacian-based Methods to Multi-echo Phase Data for Accurate Susceptibility Mapping

Emma Biondetti^{1}, David L. Thomas^{2}, and Karin Shmueli^{1}

In Susceptibility Mapping (SM), Laplacian-based techniques, e.g. SHARP^{1,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$$$ estimates^{3}.

Several studies^{4-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
studies^{7,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, others^{2,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.

Laplacian-based phase unwrapping (*Lap-Unw*)
or background-field removal (*Lap-Bg*)
were implemented with SHARP^{1,2} and a 3x3x3 Laplacian kernel^{2}. For *Lap-Unw*, non-eroded brain mask (FSL BET^{10} for the volunteer) and threshold^{11}
$$$\sigma = 10^{-10}$$$ were used. For *Lap-Bg*, the same brain mask with 2-voxel^{12} erosion and threshold^{2}
$$$\sigma = 0.05$$$ were used*.*

Phantom

Wrapped phase ($$$\phi_0 = 0$$$) was simulated at five echoes (TE_{1}/ΔTE
= 10/10 ms) from a ground-truth susceptibility distribution^{13,14} (Zubal phantom^{15}, Figure 1), using a Fourier-based forward model^{16}.
Background-field free field maps were then calculated using three distinct pipelines:

- **Avg-Unw**: *Lap-Unw*^{1}
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-Bg*^{1} on $$$\Delta B$$$.

- **Avg-Bg**: *Lap-Bg*^{2}
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-Bg*^{2} on the fitted $$$\Delta B$$$^{2,9}.

$$$\chi$$$ maps were calculated (TKD^{17} with correction for underestimation^{2} 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.

Volunteer

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 (TE_{1}/Δ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 PRELUDE^{18}. The mean and SD
of the processed phase were calculated in three ROIs (Figure 3a) drawn on the
fifth-echo magnitude image.

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 SHARP^{2} 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 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. http://weill.cornell.edu/mri/pages/qsm.html

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 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)

1547