To propose a general phase unwrapping method that exploits smoothness of the field map in 3-D and provides the globally optimal solution without over-smoothing effect. The proposed scheme aims to produce a robust solution for field map estimation that will further improve the quality of QSM.

Phase unwrapping is an important topic in many fields of study. It can be mathematically expressed as recovering the unwrapped phase $$$ω_{0}({\mathbf r})$$$ from the wrapped phase $$$ω({\mathbf r})$$$ at each location $$${\mathbf r}$$$ of an object (e.g. an image):

$$ω_{0}({\mathbf r}) = ω({\mathbf r}) + 2\pi k({\mathbf r}) \space \space \space [1]$$

$$$k$$$ is a positive integer and $$$ω({\mathbf r}) \in [0, 2\pi]$$$. For the B0 field map estimation in MRI with the presence of fat and/or water, state-of-the-art methods (e.g. region-growing, region-merging or iterative graph cut) recover the unwrapped field map value by incorporating the smoothness of the field map into either an energy minimization problem formulation or a morphological method. These methods have proven effective in many cases, but fail to provide the global minimum in challenging cases commonly due to the presence of a large field inhomogeneity, anatomical discontinuities or low SNR in certain regions of the map (e.g. boundaries). The GlObally Optimal Surface Estimation (GOOSE)¹ and 3D GOOSE² introduced previously by the same authors are able to provide the globally optimal solutions, but could also possibly cause over-smoothing due to the strict constraint imposed on the graph. Computational cost is also a concern for 3D-GOOSE when the dataset contains a large number of slices in the 3rd dimension such as in QSM. Now we propose to exploit the smoothness of the field map in 3-D by minimizing the difference between all neighboring voxels while allowing full interactions among voxels in neighboring columns in the graph. The graph is constructed with only keeping in the bottom layer and $$$k$$$ layers above. The estimate of the field map is formulated as the following optimization scheme:

$${\mathbf k} = \arg\min_{\mathbf k} \sum_{\mathbf r} \nabla | ω_0({\mathbf r}) | = \arg\min_{\mathbf k} \sum_{\mathbf r} \nabla | ω({\mathbf r}) + 2\pi k({\mathbf r}) | \space \space \space [2]$$

$$${\mathbf k}$$$ is the vector matrix for $$$k({\mathbf r})$$$ at all locations. $$$ω({\mathbf r})$$$ is the B0 field value that minimizes a predefined likelihood measure for the field map at location $$${\mathbf r}$$$ in one period: $$$[0, 2\pi]$$$. For example it is obtained using the VARPRO formulation in GOOSE. The difference between field map values at adjacent voxels is modeled using a linear-convex constraint. An illustration of the interaction among the nodes in the graph can be seen in Fig. 1. Eq. [2] is solved by using the globally optimal surface search algorithm as introduced in Wu's work³. in which the theoretical proof of the global minimum of Eq. [2] is achieved. Meanwhile, the over-smoothing effect is well-prevented by the new formulation since the smoothness assumption is built on the independently-precalculated field values. The new formulation greatly reduces the number of layers (usually to single digit) in the graph therefore enables a great saving of computation during graph search.

The proposed method is run with the Cornell QSM 3T datasets and compared with the region-growing (marked in Fig. 2) method⁴ in the Cornell QSM Toolbox (http://weill.cornell.edu/mri/pages/qsm.html). As we can see in Fig. 2, the field maps obtained from the region-growing method and the proposed method are almost identical for the 28th slice. However, for the 21st slice, the region-growing method produces artifacts at multiple regions with high field value variations (near boundary), whereas the proposed method delivers a homogeneous field map. One possible explanation is that the region-growing scheme fails to plant seeds in the correct locations in regions with low SNRs but a high field variation. The difference can also be visualized in the QSM maps. The artifacts near boundary might undermine the capability of QSM providing useful information to physicians when it comes to measuring the cavity after tumor resection, as well as observing tumors at the edge of the brain. The proposed algorithm also provides the computational gain in the experiment. We observe it on average takes $$$k = $$$ 2 or 3 to obtain a fully recovered 3D field map which costs ~5mins for a dataset of size 256*256*64. A fast, linear-convex constrained 3-D graph search algorithm for phase unwrapping is proposed. The scheme shows promises with current MRI datasets. The method needs to be further trained on other challenging datasets such as torso, abdomen or in high magnetic field.