A Fast and Globally Optimal 3-D Graph Search Algorithm for Phase Unwrapping in MRI with Applications in Quantitative Susceptibility Mapping (QSM)

Chen Cui^{1}, Abhay Shah^{1}, Cameron M. Cushing^{2}, Vince Magnotta^{3}, Xiaodong Wu^{1,2}, and Mathews Jacob^{1}

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)^{1} and 3D GOOSE^{2} 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^{3}. 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.

[1] Cui, C., Wu, X., Newell, J. D., & Jacob, M. (2015). Fat water decomposition using globally optimal surface estimation (GOOSE) algorithm. Magnetic Resonance in Medicine, 73(3), 1289-1299.

[2] Cui C., Wu, X., Newell, J. D., & Jacob, 3D Globally optimal surface estimation algorithm for fat and water. ISMRM Proceedings, Milan, Italy, 2014.

[3] Wu, X., & Chen, D. Z. (2002). Optimal net surface problems with applications. In Automata, languages and programming (pp. 1029-1042). Springer Berlin Heidelberg.

[4] Wang, Y., & Liu, T. (2015). Quantitative susceptibility mapping (QSM): decoding MRI data for a tissue magnetic biomarker. Magnetic Resonance in Medicine,73(1), 82-101.

Figure 1: 3D optimal
graph surface search algorithm illustration. At each slice, a 3-D graph (x, y, k) is built in which node a
is allowed to interact with 1) all the nodes (marked in red) in the neighboring
columns in both x and y directions, 2) with all nodes (marked in red) in the
corresponding columns of neighboring slices (z-1, z+1) using the linear convex
constraint. S1 is the optimal surface for the slice z. The 3D graph search will eventually generate a 3-D volume of field maps.

Figure 2: Field maps
and QSM results using both the proposed 3-D method and the region-growing
method in the Cornell Toolbox on a brain dataset, with ∆TE = 5.9ms.

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)

3272