It is of clinical interest for many applications to get accurate fat-fraction quantification or separated fat/water-images [1,2]. To determine an accurate fat/water fraction by chemical shift based methods (Dixon [3]) using multi-echo GRE acquisition it is necessary to estimate the underlying B₀-field as well. Because of strong B₀-gradients or even discontinuities at some tissue interfaces, this can be very challenging. Many approaches were already presented to overcome this problem like IDEAL [4] or variants of it [5-7]. Most of them do not allow discontinuities in the B₀-field or use rough models that still could lead to artefacts near tissue boundaries. In this work we propose a novel approach which uses a global piecewise smoothness constraint for the underlying B₀-field.

The signal $$$s_q\left(TE_n\right)$$$ at echo time $$$TE_n$$$ in voxel $$$q$$$ can be described by

$$s_q\left(TE_n\right)=\left(\rho_{W,q}+\rho_{F,q}\sum_{m=1}^M\alpha_me^{j2\pi\Delta\,\!f_mTE_n}\right)e^{j2\pi\gamma\Delta\,\!B_{0,q}TE_n}e^{-R_{2,q}^*TE_n}$$

where $$$\rho_{W,q}$$$ is the water-signal fraction, $$$\rho_{F,q}$$$ is the fat-signal fraction, $$$\Delta\,\!B_{0,q}$$$ is the static field inhomogeneity, $$$R_{2,q}^*$$$ is the reciprocal of the transverse-relaxation-time-constant $$$T_{2,q}^*$$$ in voxel $$$q$$$ and $$$\gamma$$$ is the gyromagnetic ratio. As a fat-signal model we used the multi-peak model of the ISMRM fat/water-challenge 2012, where $$$\Delta\,\!f_m$$$ is the chemical shift of peak $$$m$$$ and $$$\alpha_m$$$ is its normalized signal contribution $$$\left(\sum_{m=1}^M\alpha_m=1\right)$$$, which are assumed to be known.

To solve this problem we propose a 3-step procedure: First, we estimate $$$R_{2,q}^*$$$ roughly by a linear fit on $$$\log\vert\,\!s_q\left(TE_n\right)\vert$$$ in each voxel. As described in [6] a rough estimate for $$$T_2^*$$$ is usually sufficient for this purpose. In the second step a spatially piecewise smooth solution for $$$\Delta\,\!B_{0,q}$$$ is calculated. If we require a constant echo spacing $$$\Delta\,\!TE$$$

$$TE_n=TE_1+\left(n-1\right)\Delta\,\!TE$$

and take $$$R_{2,q}^*$$$ out of step1, we can rewrite the least squares residual

$$r_q^2\left(\Delta\,\!B_{0,q}\right)=\min_{\rho_{W,q},\rho_{F,q}}\sum_{n=1}^N\left\|s_q\left(TE_n\right)-\left(\rho_{W,q}+\rho_{F,q}\sum_{m=1}^M\alpha_me^{j2\pi\Delta\,\!f_mTE_n}\right)e^{j2\pi\gamma\Delta\,\!B_{0,q}TE_n}e^{-R_{2,q}^*TE_n}\right\|_2^2$$

as a trigonometric polynomial in $$$z$$$ and $$$z^{-1}$$$ by the following substitution:

$$z_q=e^{j2\pi\gamma\Delta\,\!B_{0,q}\Delta\,\!TE}\qquad\,r_q^2\left(\Delta\,\!B_{0,q}\right)\rightarrow\,r_q^2\left(z_q\right).$$

This formulation was first described in [7] in the case of three echoes, while our proposed solution is able to deal with an arbitrary number of echoes. All local minima of this polynomial can be determined easily by root-finding in each voxel and because of their finite number, also the global one. Due to B₀-inhomogeneity and noise the global minimizer is not always correct, which can lead to fat/water-swaps. We now reformulate the given problem into a variational problem extended by a regularization term. Piecewise smoothness is enforced by the total-generalized-variation (TGV) functional introduced in [8] which has proven to be beneficial for MRI-applications in [9]. The optimization problem can be written as

$$\min_z\int_\Omega\,r_q\left(z_q\right)\,dq+TGV_\alpha^2\left(z\right)\qquad\,s.t.\,z_q\,\,\text{local}\,\,\text{minimum}\,\,\text{of}\,\,r_q$$

meaning that the integral over the residuals restricted to the local minima of the polynomial on the unit-circle is minimized s.t. TGV constraints. The overall goal is to select the local minimum which gives the globally best estimate (Fig.1A). Because of the combinatorial complexity, the convex relaxation of the problem is solved, where we extend the solution space and each functional $$$r_q$$$ to its convex-hull (i.e., all solutions in-between) (Fig.1B)[10].

$$\min_z\int_\Omega\,r_{q,relax}\left(z_q\right)\,dq+TGV_\alpha^2\left(z\right).$$

This problem can be solved easily using the primal-dual algorithm described in [11]. The solution is then back-projected to the nearest local minimum on the unit-circle. Finally, in the third and last step we solve the signal-equation using a least-squares-fit to determine $$$\rho_{W,q}$$$ and $$$\rho_{F,q}$$$ in each voxel.

This method was applied to datasets of the ISMRM fat/water-challenge 2012 and datasets acquired on a clinical 3T system (Skyra, Siemens, Erlangen Germany) with 4 equally-spaced echoes. The acquisition was performed with the following sequence-parameters: FOV=300mm, TR=100ms, $$$\alpha_{flip}=25°$$$, matrix-size=256x256, slice-thickness=3mm, $$$TE_1=1.86ms$$$ and $$$\Delta\,\!TE=3.32ms$$$. The results are compared to an IDEAL-implementation according to [12].

Fig.2 shows fat-fraction, field-map, fat- and water-image of the acquired dataset using our method compared to that achieved with the IDEAL-implementation. Fig.3 and Fig.4 show the same quantities for dataset 7 and 8 of the ISMRM fat/water-challenge 2012, respectively.

The separation accuracy is considerably improved compared to IDEAL, where the results are corrupted by wide fat/water-swaps, rendering the result unusable. The same can be seen on the challenge results, especially dataset 8 (Fig.4), where IDEAL produces fat/water-swaps nearly periodically over the whole FOV, which is completely removed with our method. Further, the B₀-field is estimated uncorrectly in wide areas in all shown datasets using IDEAL, which cause the described fat/water-swaps. We showed the capability of our method to deal with huge B₀-inhomogeneities and its ability to avoid fat/water-swaps. Further, we are able to consider small susceptibility-changes and the corresponding jumps in the B₀-field on tissue-boundaries parallel to the main-field orientation.