There exist several 3-point Dixon reconstruction methods. However a common problem is so called water-fat swaps, especially at low SNR. This can limit correct clinical assessments and also quantitative measurements. The purpose of the method suggested herein is to give correct results even at low SNR.

The signal model used is described in eq.1.

$$S_1=W+a_1F\hspace{20mm}\\S_2=(W+a_2F)b\hspace{10mm}[1]\\S_3=(W+a_3F)b^2\hspace{13mm}$$

where S_n is the expected signal of echo n. The echoes have a constant time spacing, Δt.

W and F are the complex signals coming from water and fat respectively at the time of excitation multiplied by a phasor. A 9-peak fat model¹ is used to calculate a_n. Field inhomogeneities and signal decay are modeled by b, which can be written as:

$$b=e^{(iω-R_2^*)Δt}\hspace{20mm}[2]$$

The value of b for each voxel is unknown, but two possible values can be calculated analytically using a formula derived from eq.1, the two solutions are each given a label, 1 or 2.

$$b_{1,2}=\dfrac{(a_1-a_3)S_2}{2(a_2-a_3 )S_1 }±\sqrt{{\Bigg(\dfrac{(a_1-a_3)S_2}{2(a_2-a_3 )S_1 }\Bigg)}^2-\dfrac{(a_1-a_2)S_3}{(a_2-a_3 )S_1 }}\hspace{10mm}[3]$$

Finding the correct label is done using a global image optimization with unary and binary terms. For each label of b a corresponding R (=W/F) is calculated:

$$R_1=R(b_1,b_2)=\dfrac{a_3(a_1-a_2)b_1+a_1(a_3-a_2)b_2}{(a_2-a_1)b_1+(a_2-a_3 )b_2}\hspace{20mm}[4]$$

R₂ is calculated by swapping b₁ and b₂ in eq.4.

Since the phases of W and F are expected to be equal in each voxel the phase of R should be close to 0. Using this knowledge, a unary cost can be assigned for each voxel.

Binary costs for each of the four possible combinations of labels between neighbors are also calculated. Based on the idea that the phase map should be spatially smooth, the cost gets lower the closer the phases are to each other. The unary and binary costs are both taken into account, and a single parameter, λ, is used to set the relative importance between them, resulting in a single energy function (eq.5) of the entire image that is to be minimized.

$$E=\sum_{p∈V}-λ\cos(∠R_l^p)+\sum_{(p,q)∈\mathcal{E}}-\dfrac{m_p m_q}{d_{p,q}}\cos(∠b_l^p-∠b_l^q)\hspace{10mm}[5]$$

In eq.5 E is the energy to be minimized, p and q are voxel indices, V is the set of all voxels, $$$\mathcal{E}$$$ is the set of all 6-neighborhood voxels, m is the magnitude weight of the voxels², d is the Euclidian distance between two voxels in mm and l is the label of the voxel.

The magnitude weights are used to lessen the effect of the background noise. Eq.5 is minimized using QPBO³ which is a graph-cut method. The parameter λ can be automatically set. After finding the correct label, estimates of W and F can be calculated using a least-squares estimate. $$$R_2^*$$$ can be calculated from b.

The method presented will be referred to as analytical 3-point Dixon graph cut (AGC).

The resulting images were compared to those from a reference method, ASR², which uses region-growing to find the correct label. Noise performance was tested on the images by adding noise to set SNR to different levels (1 2 4 8 16). Images were compared both visually and automatically. The automatic method considers a voxel to be incorrectly labelled if the phase map value is significantly different from the phase map value when reconstructed at full SNR.

The properties of the images used can be found in table 1. The images were whole-body scans of 36 adult humans.

In Fig.1 reconstructed images using the two methods for different SNR can be seen. In Fig.2 the estimated percentage of voxels with the incorrect labels are presented along with p-values for the difference in the percentage of errors between the methods. For SNR 8 the p-value shows no significance, although visual inspection clearly shows more errors for ASR than for AGC. For SNR 16 no significant difference was detected automatically or visually.

The computation time for AGC was less than 1 minute on a modern desktop computer for the whole-body test images at full SNR.

The results show that the suggested method is superior to the reference method at low SNR when it comes to the amount of swaps, while maintaining a low computation time. The method only uses one parameter, which can be set automatically.

The method used for estimating if the incorrect label has been assigned in a voxel is not perfect. For example the automatic method shows no difference in performance between the methods for SNR 8, in contrast to visual inspection.

AGC reconstructs data with low SNR correctly which demonstrates the method's robustness. It can potentially give clinically and quantitatively relevant information where the reference method fail. The short reconstruction time should allow for online reconstruction.