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³.

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.

Here, we applied three processing pipelines^2,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 SHARP^1,2 and a 3x3x3 Laplacian kernel². For Lap-Unw, non-eroded brain mask (FSL BET¹⁰ for the volunteer) and threshold¹¹ $$$\sigma = 10^{-10}$$$ were used. For Lap-Bg, the same brain mask with 2-voxel¹² erosion and threshold² $$$\sigma = 0.05$$$ were used.

Phantom

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

- Avg-Unw: Lap-Unw¹ on the phase at each TE; field map calculation:⁵ $$\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¹ on $$$\Delta B$$$.

- Avg-Bg: Lap-Bg² on the phase at each TE; field map calculation:⁷ $$ \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² on the fitted $$$\Delta B$$$^2,9.

$$$\chi$$$ maps were calculated (TKD¹⁷ with correction for underestimation² 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₁/Δ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¹⁸. 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² 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 pipelines¹³ that fit the phase over multiple echoes, should be used before applying Laplacian-based methods.