Dixon techniques^[1] offers robust water/fat separation in the presence of static field inhomogeneity. Two-point Dixon imaging with flexible echo times is desirable because of its high scan efficiency and flexibility^[2,3]. A major challenge in two-point Dixon imaging is how to estimate the phase error resulting from field inhomogeneity. Region growing^[4] and regional iterative phasor extraction (RIPE)^[5] are popular algorithms to resolve the phase ambiguity. Graph cut^[6] and tree-reweighted message-passing^[2] formulate a global optimization to determine the phase error, and have shown to be robust in areas of large field inhomogeneity. However, solving the global optimization can be computationally intensive, especially for 3D imaging. In this work, we formulate a binary quadratic optimization and propose a fast projected power method to resolve the phase ambiguity in two-point Dixon imaging.

Following the notation in [3], two-point Dixon imaging can be modeled as:

$$S_1=(W+c_1F)e^{i\phi_1},S_2=(W+c_2F)e^{i\phi_2}$$

where $$$W$$$ and $$$F$$$ are water and fat signals in image space, $$$S_1$$$ and $$$S_2$$$ are complex composite signals acquired at two echo times, $$$c_1$$$ and $$$c_2$$$ are the dephasing factors with respect to water (usually known $$$a~priori$$$) corresponding to the echo times using a single peak or multiple peak fat model. The main problem here is to estimate the phase error (or phasor) between two echoes, $$$P=e^{i(\phi_2-\phi_1)}=e^{i\Delta \phi}$$$.

Two pairs of solutions $$$(W_1,F_1)$$$ and $$$(W_2,F_2)$$$ can be calculated analytically^[3], yielding two potential solutions for $$$P$$$:

$$P_1=\frac{S_1^*S_2}{(W_1+c_1^*F_1)(W_1+c_2F_1)},P_2=\frac{S_1^*S_2}{(W_2+c_1^*F_2)(W_2+c_2F_2)}$$

Because field inhomogeneity is usually smooth, the problem of selecting the correct phasor can be formulated as a global optimization^[2,6]:

$$\textrm{minimize}~\sum_{r}\sum_{{s}\in\mathcal{N}_{r}}v(P({r}),P({s}))$$

$$\textrm{subject}~\textrm{to}~P({r})\in\left\{P_1({r}),P_2({r})\right\},\forall{r}$$

where $$${r}$$$ is an index of image pixels, $$${s}$$$ presents pixels in the neighboring set ($$$\mathcal{N}_{r}$$$) of $$${r}$$$, and $$$v$$$ is a discontinuity cost function for phasor pair $$$(P({r}),P({s}))$$$. For simplicity, second-order neighborhood^[6] (26 neighboring pixels in 3D) was adopted in this work to promote spatial smoothness: $$$v(P({r}),P({s}))=|P({r})-P({s})|$$$.

To choose the optimal $$$P$$$, we can reformulate the previous problem into a binary quadratic optimization. As demonstrated in Fig. 1, define variable $$$X({r})=\left(\begin{array}{c}X_1(r)\\X_2(r)\end{array}\right)\in\left\{{\left(\begin{array}{c}0\\1\end{array}\right),\left(\begin{array}{c}1\\0\end{array}\right)}\right\}$$$ for each pixel, and stack them all into a vector $$$X$$$. Then the previous problem can be reformulated as^[7]:

$$\textrm{minimize}~f(X)=X^{T}VX$$

$$\textrm{subject}~\textrm{to}~X({r})\in\left\{{\left(\begin{array}{c}0\\1\end{array}\right),\left(\begin{array}{c}1\\0\end{array}\right)}\right\},\forall{r}$$

where $$$V$$$ contains the possible values of $$$v$$$ with respect to $$$X$$$. Note that $$$V$$$ is large but sparse, and $$$VX$$$ can be implicitly represented as simple pixel-wise operations. Here, an iterative projected power method is used to solve this problem. The major steps are summarized below:

(1) Set $$$P$$$ to unity and initialize $$$X$$$ based on:

$$X({r})=\begin{cases}\left(\begin{array}{c}1\\0\end{array}\right),&|P({r})-P_1({r})|<|P({r})-P_2({r})|\\\left(\begin{array}{c}0\\1\end{array}\right),&\textrm{otherwise}\end{cases},\forall{r}$$

(2) Update $$$X\leftarrow{VX}$$$, then set $$$X(r)\leftarrow{\left(\begin{array}{c}\frac{X_2(r)}{X_1(r)+X_2(r)}\\{\frac{X_1(r)}{X_1(r)+X_2(r)}}\end{array}\right)}$$$. Go to (4).

(3) Update $$$X\leftarrow{VX}$$$, then set $$$X({r})=\begin{cases}\left(\begin{array}{c}1\\0\end{array}\right),&X_1(r)<X_2(r)\\\left(\begin{array}{c}0\\1\end{array}\right),&\textrm{otherwise}\end{cases}$$$. Go to (4).

(4) If the iteration count is less than 50, repeat (2); otherwise repeat (3) until $$$X$$$ converges.

The final phasor is $$$P({r})=X(r)^T[P_1(r);P_2(r)],\forall{r}$$$. Water and fat can be separated by solving:

$$\begin{cases}W'+c_1F'=S_1\\W'+c_2F'=S_2P^*\end{cases}$$

where $$$W'=We^{i\phi_1}$$$, and $$$F'=Fe^{i\phi_1}$$$.

The proposed method as well as the method described in [3] using RIPE were performed on a 3D $$$in~vivo$$$ pelvic dataset acquired on a 3T GE MR750 scanner. The acquisition parameters were TE1/TE2/TR 2.2/3.3/5.9 ms, BW $$$\pm$$$ 100 kHz, FOV 30$$$\times$$$30$$$\times$$$23 cm$$$^3$$$, and matrix 354$$$\times$$$320$$$\times$$$230. Both single peak and multiple peak (six peaks) fat models were compared for each method. The dataset was downsampled first to a 3$$$\times$$$3$$$\times$$$3 mm$$$^3$$$ voxel size before phasors were calculated and selected using different methods. The resolved phasor was then upsampled to the original voxel size for the final water/fat separation.

Water and fat images are shown in Figs. 2 and 3. The total reconstruction time in Matlab with 2 Intel Xeon E5-2640 processors was 13.4 seconds for RIPE and 98.7 seconds for the proposed method. A complete global water/fat swap was observed for RIPE with both fat models, which indicated the wrong phasor candidate was selected. Water and fat images were manually switched and labeled as RIPE1. No global swap was observed with the proposed method. For complete comparison, RIPE was also performed with the phasor manually set to the other phasor candidate, and the results were referred to as RIPE2. As shown by the arrows, the proposed method also reduced local water/fat swaps in places with large field inhomogeneity. When the phase error was correctly resolved, better water/fat separation was obtained using the multiple peak fat model.In this work, we have formulated a binary quadratic optimization and developed a fast projected power method to resolve the phase ambiguity in two-point Dixon imaging. Compared to RIPE, improved water/fat separation can be achieved. A quadratic penalty^[2,6] of $$$v$$$ with spatial weighting can also be applied for the proposed method, but it may require more iterations to converge.