MR images hold complex-valued information whose phase component $\phi$ needs to be examined in certain scenarios. Computing $\phi$ by the standard function atan2(Im,Re) only provides phase principal values, $-\pi\lt\phi\leq+\pi$, while the dynamic range of the true phase values may extend well beyond this range. Phase unwrapping is necessary to remove artefactual phase jumps associated with principal value estimation, while being robust to jumps of physical origin. Number of techniques have been proposed to handle this problem (e.g. [1-2]). They generally rely on complex methods and often require a priori knowledge on the underlying data. To accelerate the process of multidimensional phase unwrapping, we devised a new approach based on a linear algebraic formulation and well adapted to modern parallel processors.

$\mathbf{\text{Principle}}$ Phase derivatives robust to phase jumps can be computed from the real and imaginary components of a signal (see below). We use this property, followed by integration to recover phase estimates robust to phase jumps. The integral is however only determined up to an unknown integration constant, which needs to be estimated. For multidimensional signals, this approach can be applied independently along each dimension. Ensuring consistency between multidimensional phase maps estimated by differentiating-then-integrating along orthogonal axes provides a means to determine the unknown integration constants introduced in the process.

$\mathbf{\text{Robust phase gradient}}$ Normalized complex numbers can be represented as $e^{\imath\,\phi}$. For such numbers, we have: $$\nabla_\alpha\left(e^{i\,\phi}\right)\,=\,i\,\nabla_\alpha\,{\phi}\,\cdot\,e^{i\,\phi}\quad\Longleftrightarrow\quad\nabla_\alpha\,{\phi}\,=\,\frac{-i\,\nabla_\alpha\left(e^{i\,\phi}\right)}{e^{i\,\phi}}\qquad(1)$$ where $\nabla_\alpha$ represents the gradient operator along axis $\alpha$. Normalizing each element of a complex MR image $I(\overrightarrow{r})$ does not affect the phase and it is possible to estimate $\nabla_\alpha\,{\phi}$ from the normalized complex image $\frac{I(\overrightarrow{r})}{|I(\overrightarrow{r})|}\,\equiv\,e^{i\,\phi(\overrightarrow{r})}$ and its corresponding gradient $\nabla_\alpha{(I/\left|I\right|)}$.

$\phi$ is real, and so should be $\nabla_\alpha\,{\phi}$. Eq.1 for $\nabla_\alpha\,{\phi}$ involves complex quantities but the imaginary part of the result should be negligible (Figs.3-1,2). We use this principle as a test to detect remaining phase discontinuities of physical origin, that should be allowed to subsist in a phase unwrapped map.

$\mathbf{\text{Linear system of equations}}$ We detail the process for a 2D, $N_x\,\times\,N_y$ image. A version of the original image normalized per element,$\frac{I(\overrightarrow{r})}{|I(\overrightarrow{r})|}$, is produced. For each of its $N_x$ rows, the unidirectional $\nabla_y\,{\phi}$ gradient along $N_y$ columns is computed and integrated, $\hat{\phi}_y\,=\,\int\,\nabla_y\,{\phi}\,dy$, see Fig. 4-1, with one unknown integration constant $x_{\text{uic}}$ per row. For each element $\overrightarrow{r}$ where the imaginary part of $\nabla_y\,{\phi}$ is above some threshold while the intensity $|I(\overrightarrow{r})|$ of the original image is above some other SNR-related threshold, we consider an additional unknown integration constant $x_{\text{auic}}$ must be introduced. The process is similarly repeated for unidirectional $\nabla_x\,{\phi}$ along $N_x$ rows for each of the $N_y$ columns, see Fig.4-2.

$\hat{\phi}_x\,=\,\int\,\nabla_x\,{\phi}\,dx$ and $\hat{\phi}_y\,=\,\int\,\nabla_y\,{\phi}\,dy$ should provide equal estimates of the phase. Eventually, the difference $\delta=\hat{\phi}_x-\hat{\phi}_y$ (Fig.4-5) can be expressed as a function of the $N_{\textrm{ic}}\,=\,N_x+N_y+N_j$ integration constants, with $N_j$ the total number of additional constants $x_{\text{auic}}$. This problem can be formulated as a linear system $\mathbf{A}\,x\,=\,\delta'$, where: (1) $\delta'=vec(\delta)$ is the $(N_xN_y)\,\times\,1$ vectorized difference map; (2) $x$ is a $N_{\textrm{ic}}\times\,1$ vector holding the integration constants; (3) $\mathbf{A}$ is a $(N_xN_y)\times\,N_{\textrm{ic}}$ sparse matrix with +1 or -1 non-zero entries specifying how the integration constants contribute to the observed difference $\delta'$.

$\mathbf{\text{Weighted least-square estimation}}$ Solving for $x$ in the least squares sense gives the integration constants, used to correct the integrated phase gradient maps (Figs.4-3,4) and make them match (Fig.4-5). To bring robustness against experimental random noise in the data, the least squares problem is weighted using magnitude or preferably power information.

Fig.1 represents the average magnitude of two GE sagittal images acquired at different echo times, at 3T. Fig.2 is the phase difference map between the two echoes, obtained using atan2() and affected by phase wraps. Figs.3-1,2 show the across-columns and across-rows phase gradients imaginary part, the norm of which is an indicator of possible phase jump necessitating to introduce an additional integration constant. Figs.4-1,2 show the integrated across-columns and across-rows phase gradient maps. Figs.4-3,4 are their modified versions after solving and correcting for the unknown integration constants, while Figs.4-5,6 display the difference between them, and their average, respectively. Fig.4-7 is an adjusted version of Fig.4-6 where the unwrapped values have been rounded to ensure modulo 2π equivalence with Fig.2.Execution time is dominated by the resolution of the large linear problem. Exploiting the sparsity of the matrix $\mathbf{A}$ results in dramatic speed-up. Using the SuiteSparse library [3] with a NVIDIA Tesla K20xm allows to solve the 2D problem illustrated in a few ms.We have developped, implemented and demonstrated a new algorithm for multi-dimensional phase unwrapping.