The phase of the complex-valued MRI gradient echo data contains information on local magnetic field changes. The phase is implicitly wrapped into the range $$$\small{\left[-\pi,\pi\right)}$$$. In order to recover the field information, the phase needs to be unwrapped. Laplacian-based phase unwrapping [1-3] is a commonly used technique, particularly as a pre-processing step in Quantitative Susceptibility Mapping [4]. The method unwraps phase by calculating the Laplacian of the unwrapped phase, $$$\small{\phi_u}$$$, using the wrapped phase, $$$\small{\phi_w}$$$:

$$\phi_u = \nabla^{-2} \left( \cos \phi_w \nabla^2 \left( \sin \phi_w \right) - \sin \phi_w \nabla^2 \left( \cos \phi_w \right) \right), \tag{1}$$

where $$$\small{\nabla^2}$$$ and $$$\small{\nabla^{-2}}$$$ are the Laplacian and inverse Laplacian operations, computed using the time derivative property of the Fourier transform. Whilst this is analytically correct, it was derived with the assumption of a continuous signal and employing the continuous Fourier transform. These assumptions do not translate well to discrete signals and the use of the discrete Fast Fourier Transform (FFT), which are present in the context of MRI phase data and computation. In fact, the FFT of the discrete Laplacian operator is commonly used to calculate the Laplacian operations. Yet, even when employing this modification, the resulting unwrapped phase (for 1D) becomes

$$\phi\left[x\right] = \nabla^{-2} \left[ \sin\left(\phi_{w}\left[x+1\right]-\phi_{w}\left[x\right]\right) - \sin\left(\phi_{w}\left[x\right]-\phi_{w}\left[x-1\right]\right) \right] ,\tag{2}$$

where $$$\small{x}$$$ is the discrete image space coordinate. From this, it can be seen that a change of 0 is indistinguishable from a change of $$$\small{\pi}$$$,resulting in incorrect unwrapped phase values. We present an analytically correct alternative to this method, by correctly assuming the discrete nature of the data and computations.

We derived a discrete version of Laplacian-based unwrapping, where the second order partial derivative is

$$\frac{\delta^2\phi_{u}}{\delta t^2} = \tan^{-1}\frac{\sin\left(\phi_{w}\left[t+1\right]-\phi_{w}\left[t\right]\right)}{\cos\left(\phi_{w}\left[t+1\right]-\phi_{w}\left[t\right]\right)}-\tan^{-1}\frac{\sin\left(\phi_{w}\left[t\right]-\phi_{w}\left[t-1\right]\right)}{\cos\left(\phi_{w}\left[t\right]-\phi_{w}\left[t-1\right]\right)},\tag{3}$$

the 3D Laplacian is defined as$$ \nabla^2 \phi_{u} = \frac{\delta^2\phi_{u}}{\delta x^2} + \frac{\delta^2\phi_{u}}{\delta y^2} + \frac{\delta^2\phi_{u}}{\delta z^2}\tag{4}$$and the inverse Laplacian is calculated using the Discrete Cosine Transform:

$$\phi_{u}\left[x,y,z\right] = \text{DCT}^{-1}\left(\left[2\cos\left(\frac{\pi}{N}\left(k-1\right)\right)-2\right]^{-1}\text{DCT}\left(\nabla^2 \phi_{u}\right)\right).\tag{5}$$

We shall refer to the original algorithm as Laplacian-based Uwrapping using Continuous K-space (LUCK) and the new algorithm as Fast Unwrapping using Discrete Gradient Evaluation (FUDGE). We applied both methods to two numerical simulation datasets and an experimental in vivo MRI gradient echo data. The numerical simulations comprised: a phase image with spheres of varying phase values; a phase image containing a sphere (diameter=5mm, susceptibility=0.3ppm, B0=3T, TE=20ms). In vivo MRI brain data from a healthy volunteer (35yo, female) was acquired on a 7T research MRI scanner (Siemens, Erlangan, Germany) in accordance with local ethics (fully flow-compensated GRE, TE=7.65ms,15.3ms, TR=18ms, FA=13$$$^{\circ}$$$, iPAT=3, matrix=366$$$\scriptsize{\times}$$$316$$$\scriptsize{\times}$$$224, voxel=0.6mm isotropic, TA=7:42min). iLSQR [5] was used to calculate QSM for the second numerical simulation.

LUCK produced erroneous phase in the first numerical simulation (Fig. 1), whilst FUDGE produced phase values exact to within numerical floating point error. The second numerical simulation (Fig. 2) demonstrated similar results. The susceptibility of the sphere was estimated by the iLSQR method at 0.281ppm, 0.239ppm and 0.281ppm for the true phase, LUCK unwrapped phase and FUDGE unwrapped phase, respectively. For LUCK, this equated to 85.2% of the susceptibility calculated using the true phase. The unwrapped in vivo phase images (Fig. 3) demonstrated much less error in the FUDGE image than the LUCK image, particularly in the vicinity of the periventricular vasculature. FUDGE was 4.9 times faster to compute than LUCK on the in vivo data.The simulation results reflected the errors predicted by the analysis of the mathematical formulation of LUCK. Similarly, the lack of error in the FUDGE results reflects the expected accuracy when using a mathematical formulation based on the correct assumptions of discrete data and computation. The inaccuracies in unwrapped phase produced with LUCK have expected effects on subsequent QSM maps, with the 15% error in susceptibility being substantial, particularly for a quantitative method. The in vivo results demonstrated greater accuracy with FUDGE compared to LUCK. The presence of open-ended fringe lines, a phase feature common at high field (e.g. 7T), resulted in some residual error in the phase unwrapped with FUDGE, although to a much lesser extent than in the LUCK image. FUDGE was also much faster, due to there being only 2 FFTs in the algorithm, compared to 6 FFTs in the LUCK algorithm.The analytical, numerical and experimental results demonstrate that relying on LUCK will result in erroneous phase data and QSM maps, however data processed with FUDGE produces numerically accurate results.