Method of analysing magnetic resonance imaging images

Abstract

A method of analysing the magnitude of Magnetic Resonance Imaging (MRI) data is described. The method comprising: using the magnitude only of the multi-echo MRI data of images from the subject, where images are acquired at arbitrarily timed echoes including at least one echo time where water and fat are not substantially in-phase; fitting the magnitude of said multi-echo MRI data to a single signal model to produce a plurality of potential solutions for the relative signal contributions for each of the at least two species from the model, by using a plurality of different starting conditions to generate a particular cost function value for each of the plurality of starting conditions, where said cost function values are independent of a field map term for the MRI data; analysing said cost function values to calculate relative signal separation contribution for each species at each voxel of the images.

Claims

1. A method of analysing the magnitude of Magnetic Resonance Imaging (MRI) data from acquired MRI images to determine the relative signal contributions of at least two species to each voxel of the images, the method comprising the steps of: using the magnitude of the multi-echo MRI data of images from the subject, where the images are acquired at arbitrarily timed echoes including at least one echo time where water and fat are not substantially in-phase with each other; fitting the magnitude of said acquired multi-echo MM data to a single signal model to produce a plurality of potential solutions for the relative signal contributions for each of the at least two species from the model, by using a plurality of different starting conditions to generate a particular cost function value for each of the plurality of starting conditions, where said cost function values are independent of a field map term for the MRI data; and analysing said cost function values to calculate the relative signal separation contribution for each species at each voxel of the images.

2. An image processing system arranged to analyse the magnitude of Magnetic Resonance Imaging (MRI) data from acquired MRI images to determine the relative signal contributions of at least two species to each voxel of the images, the image processing system comprising at least one processing device arranged to: use the magnitude of the multi-echo MRI data of images from the subject, where the images are acquired at arbitrarily timed echoes including at least one echo time where water and fat are not substantially in-phase with each other; fit the magnitude of said multi-echo MRI data to a single signal model to produce a plurality of potential solutions for the relative signal contributions for each of the at least two species from the model, by using a plurality of different starting conditions to generate a particular cost function value for each of the plurality of starting conditions, where said cost function values are independent of a field map term for the MRI data; and analyse said cost function values to calculate the relative signal separation contribution for each species at each voxel of the images.

3. A method according to claim 1, wherein said analysis of said cost function values comprises the step of: comparing the generated cost function values to determine which is the correct solution for said signal separation.

4. A method according to claim 3, wherein the lowest cost function value of said species is determined to be the correct solution for said signal.

5. A method according to claim 1, wherein said magnitude of said multi-echo MRI data is fitted to said single signal model using a model fitting algorithm.

6. A method according to claim 5, wherein said model fitting algorithm is an instance of at least one of the following: least squares estimation, iteratively reweighted least squares, least trimmed squares, or other robust approaches using m-estimators or s-estimators.

7. A method according to claim 5, wherein said model fitting algorithm is combined with at least one regularisation term.

8. A method according to claim 1, in which the single signal model includes a spectral model of one of the at least two species with more than one spectral component.

9. A method according to claim 1, wherein said single signal model includes at least one of the relaxation time quantities (T.sub.1, T.sub.2, T.sub.2*) to correct for signal decay.

10. A method according to claim 8, wherein the starting condition values of the relaxation time quantities are in the physically observable range.

11. A method according to claim 9, wherein the starting condition values of the relaxation time quantity T.sub.2* are between 1 and 100 ms.

12. A method according to claim 9, wherein the starting condition values of the relaxation time quantity T.sub.2* are between 10 and 30 ms.

13. A method according to claim 1, further comprising the step of using said species signal contribution to generate separate images showing the results for each species.

14. A method according to claim 13, in which one or more of the resulting images are post-processed.

15. The method according to claim 1, wherein the separated species contributions are used to estimate a field heterogeneity (‘fieldmap’) term.

16. The method according of claim 9, wherein the estimated relaxation time quantities are used to estimate a field heterogeneity (‘fieldmap’) term.

17. The method according to claim 1, wherein the at least two species include at least two of water, fat, hyperpolarized contrast elements or metabolites of such elements, markers for the presence of cancerous cells.

18. The method according to claim 1, wherein the subject includes a phantom model or animal or human tissue comprising at least one of the following organs: liver, pancreas, kidney, spleen, heart, muscle or adipose tissue.

19. The method of claim 1, in which cost function values in a certain voxel are used to update a likelihood map of the presence of at least one species in the voxel.

20. An image processing system as in claim 2, wherein the system is a non-transitory computer program product having executable program code stored therein, the program code operable for analysing the magnitude of magnetic resonance imaging (MRI) data.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Further details, aspects and embodiments of the invention will be described, by way of example only, with reference to the drawings. In the drawings, like reference numbers are used to identify like or functionally similar elements. Elements in the figures are illustrated for simplicity and clarity and have not necessarily been drawn to scale.

(2) FIG. 1 illustrates MRI PDFF images obtained using prior art techniques;

(3) FIG. 2 shows example cost function value maps;

(4) FIG. 3 is a frequency spectrum for fat;

(5) FIG. 4 illustrates magnitude related water fat ambiguity using spectral peaks and echo shifts;

(6) FIG. 5 illustrates example cost function value maps as used in an example of the invention;

(7) FIG. 6 illustrates example MRI PDFF images obtained using the method of the invention;

(8) FIG. 7 shows MRI PDFF images obtained using the method of the invention;

(9) FIG. 8 show the relationship between measured and true PDFF for images obtained under different conditions;

(10) FIG. 9 is a flow chart showing the steps involved in the method of this invention;

(11) FIG. 10 shows the agreement (absolute and Bland Altman) between prior art measurements and measurements obtained using this invention;

(12) FIG. 11 compares MRI PDFF images obtained from prior art methods, and the method of this invention;

(13) FIG. 12 shows field maps and MRI PDFF images for the prior art method, and MRI images obtained using the method of the invention;

(14) FIG. 13 shows field maps, MRI PDFF images and distributions obtained using the prior art method and the method of this invention;

(15) FIG. 14 shows field maps, MRI PDFF images and distributions obtained using the prior art method and the method of this invention.

DETAILED DESCRIPTION

(16) The present invention will now be described with reference to the accompanying drawings in which there is illustrated an example of a method and apparatus for analysing only the magnitude data from MRI images that are acquired at arbitrarily timed echoes including at least one echo time where water and fat are not substantially in-phase with each other. The phase information is not needed in this method. However, it will be appreciated that the present invention is not limited to the specific examples herein described and as illustrated in the accompanying drawings.

(17) Furthermore, because the illustrated embodiments of the present invention may for the most part, be implemented using electronic components and circuits known to those skilled in the art, details will not be explained in any greater detail than that considered necessary as illustrated below, for the understanding and appreciation of the underlying concepts of the present invention and in order not to obfuscate or distract from the teachings of the present invention.

(18) As discussed above, most methods for species separation (where the species are preferably water and fat) in MRI images are based on complex data, that is, they rely on the availability of both magnitude and phase data in order to estimate the fieldmap, which is then used to estimate water and fat components of the image.

(19) A magnetic resonance signal s.sub.i at a single voxel containing water, fat and iron (or other species) may be sampled at multiple echo times t.sub.i during relaxation. For a general complex-valued signal, the following phase-constraint model has been proposed:
s.sub.i=(ρ.sub.W+ρ.sub.F.Math.Σ.sub.pα.sub.pe.sup.j2πf.sup.p.sup.t.sup.i).Math.e.sup.j(2πψt.sup.i.sup.+ϕ.sup.0.sup.).Math.e.sup.−R.sup.2.sup.*t.sup.i+n.sub.i  Equation 1
where ρ.sub.W and ρ.sub.F are the unknown water and fat quantities, respectively, and R.sub.2*=1/T.sub.2* (s.sup.−1) is an unknown relaxation quantity. ρ.sub.W and ρ.sub.F are real-valued variables that have associated phase terms e.sup.jϕ.sup.W and e.sup.jϕ.sup.F, where ϕ.sub.W==ϕ.sub.0 may be assumed.

(20) In the above equation s.sub.i is actually s.sub.i(x), and ρ.sub.w is ρ.sub.w(x), but explicit mention of the spatial parameter (x) clutters the equation, and so established practice is that x is supressed in the equation.

(21) The field map ψ is modelled as a phase shift. The signal is further affected by noise (n.sub.i) which is typically ignored in subsequent derivations, implicitly assuming high enough signal-to-noise ratio (SNR) acquisitions. PDFF may be calculated from the water and fat amounts using

(22) $\mathrm{PDFF}=\frac{{\rho }_F}{{\rho }_W+{\rho }_F}\times 100.$

(23) The term Σ.sub.pα.sub.pe.sup.j2πf.sup.p.sup.t.sup.i is the multi-peak fat spectral model with P peaks: α.sub.p is the amplitude of the fat peak p relative to the amplitude of the water peak, and f.sub.p is the difference in precession frequency of the fat peak p with respect to the water peak in Hz (16, 38). It is usual practice to assume that the values {a.sub.p, f.sub.p} for all p are known empirically and constant throughout the image.

(24) The alleged inability of magnitude-based water-fat separation methods to determine PDFF values above 50% is known as the ‘fat-water ambiguity’ challenge, and may be explained mathematically (Bydder et al., 2008; Yu et al., 2011). With magnitude-based methods, Equation 1 above in which |a| refers to the magnitude of a, has to be optimised for a given set of echo signals s.sub.i,
|s.sub.i|=|(ρ.sub.W+ρ.sub.F.Math.Σ.sub.pα.sub.pe.sup.i2πf.sup.p.sup.t.sup.i).Math.e.sup.−R.sup.2.sup.*t.sup.i|.  (Equation 1a)

(25) As expected, the phase term vanishes, as does the field map parameter.

(26) Note that there is a set of equations for each and every voxel in an MRI image, and the goal is to estimate the fat-water fraction

(27) $\mathrm{PDFF}=\frac{{\rho }_F}{{\rho }_W+{\rho }_F}\times 100$
fraction (given as a %) for each voxel in the image. In the magnitude only formulation of the problem, there is no fieldmap term (that is familiar from complex variations), since B.sub.0 inhomogeneity is generally modelled (at each voxel) as a constant phase shift.

(28) To illustrate the fundamental problem, we assume for the moment that fat is modelled as having a spectrum with a single-peak (p=1), so that Σ.sub.pα.sub.pe.sup.i2π.sup.p.sup.t.sup.i becomes e.sup.i2πΔft.sup.i, whose magnitude |e.sup.i2πΔft.sup.i|=1. Equation 1 now becomes
|s.sub.i|=|(ρ.sub.W+ρ.sub.F).Math.e.sup.−R.sup.2.sup.*t.sup.i|.  (Equation 2)

(29) The fat-water ambiguity that is reportedly inherent to magnitude-based methods is immediately evident in Equation 2 in that exchanging ρ.sub.W and ρ.sub.F does not change the value of the equation (Yu et al., 2011). Both solutions are equally valid even though in reality there will only be one “true” solution {p.sub.W, p.sub.F}.sub.t, the other being an aliased solution {ρ.sub.W, ρ.sub.F}.sub.a.

(30) The two solutions are simply related, since PDFF.sub.a=100%−PDFF.sub.t holds. Note that R.sub.2* will be equal for both solutions. We may simulate the optimisation space for: a given field strength (1.5 T); an arbitrary set of echoes ({1.2, 3.2, 5.2, 7.2, 9.2, 11.2} ms); and true values {PDFF,R.sub.2*}.sub.t={75%, 45 s.sup.−1}. This field strength is chosen as the strength of a generally available MRI scanner, with 3 T considered a “high field”. By generating a true noise-free signal s.sub.i|{PDFF, R.sub.2*}.sub.t and then comparing the true signal with a signal generated at each possible {PDFF,R.sub.2*} point in space, a continuous map of ‘cost function values’ is produced.

(31) FIG. 2 shows examples of cost function value maps obtained under different conditions, when the optimisation space is reduced to two dimensions (ρ.sub.W and ρ.sub.F may be combined into PDFF for dimensionality reduction to allow visualisation; and R*.sub.2) when using magnitude-based estimation. FIG. 2(a) shows the simplest case where the fat is modelled as a single-peak spectrum. In this case there is unresolvable fat-water ambiguity, reflected as symmetry in the optimisation space. FIG. 2(b) is the result when fat is modelled as a multi-peak spectrum. As seen in this image, ambiguity still exists and is reflected as a local minimum but there is no longer symmetry about PDFF=50% in the image. The ambiguity in FIG. 2(b) may be resolvable using the multi-point search optimisation of this invention, with different initial conditions, where the dot in the top right-hand quadrant of the image of FIG. 2(c) indicates the chosen solution.

(32) When the spectrum of fat is assumed to consist of a single peak, there will be two minima, {ρ.sub.W, ρ.sub.F}.sub.t and {ρ.sub.W, ρ.sub.F}.sub.a, and the unresolvable ambiguity will be reflected by a symmetrical optimisation space as shown in FIG. 2(a) The choice of which of these two local minima is found by an optimisation algorithm will be determined by and is sensitive to the initial condition. Specifically, if, as is almost always the case in previously reported methods, the algorithm is initialised near to PDFF=0% then it will always converge to values below 50%.

(33) Now suppose that fat is modelled more realistically as having a multi-peak spectrum, as shown in FIG. 3, for example (p=6; Hamilton et al., 2011), then in Equation 1, |Σ.sub.pα.sub.pe.sup.i2πf.sup.p.sup.t.sup.i| will no longer be equal to 1, and exchanging ρ.sub.W and ρ.sub.F will change the value of the equation 2 (Yu et al., 2011). In the optimisation space, the ‘swapped’ solution {ρ.sub.W, ρ.sub.F}.sub.a now appears as a local minimum, where PDFF.sub.a≠100%−PDFF.sub.t; the converged R.sub.2* will now be different in both solutions (as shown in FIG. 2(b). Note that in this multi-peak situation, an optimisation algorithm initialised near PDFF=0% will still converge to values below 50%. Conversely, the same optimisation algorithm initialised near PDFF=100% will converge to the local solution greater than 50%. This ambiguity now is resolved by comparing the cost function values at both minima—achieved with at least two different sets of initial conditions—and choosing the solution with the lowest cost function value (FIG. 2(c)).

(34) In an example of the invention we use a multi-peak spectral model, preferably the spectrum is for fat, with P=6, as this is the value most commonly used in the literature, however spectral model(s) for at least one of the two species may be used as long as the model has more than one spectral component. We may now move away from the in-phase acquisition by shifting the echo sampling by a known shift δ, so that t.sub.i=(i−δ)/Δf; δ=0.5 thus meaning an opposed-phase acquisition. This is illustrated with reference to FIG. 4. The left-hand column of this figure is results for using in-phase echoes with variable amount of P, whereas the righthand column is for results using variable arbitrarily timed echoes, and P=6 in accordance with an example of this invention. That is, the magnitude only of the multi-echo MRI data is used, with MRI images acquired at arbitrarily timed echoes including at least one echo time where water and fat are not substantially in-phase with each other. The top row shows graphs for a water dominant voxel (PDFF 10%), the central row shows the results for a fat dominant voxel (PDFF 90%). While PDFF may be resolved in the 0-100% range for substantial δ, a local minimum appears for increasing δ; in general, for an asymmetric combination of echo times, there are two local minima for the physical range of PDFF and R.sub.2* (FIG. 6 graphs in right hand column). Note that optimisation algorithms initialised near 0% PDFF will converge to the correct solution for water-dominant voxels (FIG. 5 top and bottom images on right hand side), but they will converge to the local minimum for fat-dominant voxels (FIG. 5 2.sup.nd and bottom images on right hand side). Magnitude-based estimation methods previously described in the literature have consistently exhibited this behaviour, with true PDFF>50% values aliased to solutions below 50%. it is understood that the method disclosed here works so long as P>1. If 1<P<6 then the minimum we seek is less pronounced. If P>6 the extra spectral peaks, which are very low (e.g. give height of the 6th peak relative to the first one) then there is little or no additional benefit since the MRI data is discrete, relatively coarsely sampled, and subject to noise.

(35) The magnitude-related water-fat ambiguity may now be resolved because the solutions at the local minima will now have higher cost function values than the true solutions, using the presented magnitude multi-peak model. We present a multipoint method that aims to explore both possible solutions and resolve the magnitude-related water-fat ambiguity for arbitrarily timed echoes including at least one echo time where water and fat are not substantially in-phase with each other.

(36) In summary, the classical magnitude-based related fat-water ambiguity challenge may be resolvable by using: magnitude data only, to reduce the possible minima from many (using the complex-based approach that requires both magnitude and phase; see ref 24) to only two potential minima; multi-peak fat spectrum modelling, to break the symmetry and thus making the two possible minima distinguishable; and an optimisation technique that ensures both minima are explored, to compare the cost function values at both solutions and choose the solution associated with the lowest cost function value.

(37) Our method has been applied to an example range of different clinical cases, and its accuracy, precision and robustness to artefacts (especially fat-water swaps) has been evaluated.

(38) The accuracy of the magnitude-based method of this invention (MAGO) was tested in vivo against an in-house implementation of the prior art IDEAL method, for two different sets of data, where the in-house IDEAL (LMS IDEAL) method was considered the reference standard. LMS IDEAL has been previously validated against phantom data and a set of in vivo data.

(39) In order to ensure that both possible solutions are explored, the initial values of water and fat quantities need to be combined to a low PDFF in at least one run of the algorithm and to a high PDFF in at least one other run. The initial value for the relaxation quantity R.sub.2* may be set within the physiologically expected range in all runs. A given converged solution set (ρ.sub.W, ρ.sub.F, R.sub.2*) has an associated cost function value in the form of the expression cost function value (ρ.sub.W, ρ.sub.F, R.sub.2*)=Σ.sub.i.sup.N(|ŝ.sub.i|−|s.sub.i|).sup.2, where ŝ.sub.i is the estimated signal using the converged solution set (ρ.sub.W, ρ.sub.F, R.sub.2*) in Equation 2.

(40) Note that this definition of the cost function value is independent of a field map term. We choose as solution at each voxel the one with the lowest cost function value, though the other minima may be retained as alternative solutions.

(41) An example of the magnitude-based method of this invention was implemented using the Isqcurvefit function in the Matlab program (The MathWorks, Inc.). Alternatively, the method may be performed using compiled C++ routines using ITK (www.itk.org). Nonlinear fitting of the data was performed twice, each time with a different set of initial conditions for the data. As discussed above, the data is the magnitude only of the multi-echo MRI data of images from the subject, where the images are acquired at arbitrarily timed echoes including at least one echo time where water and fat are not substantially in-phase with each other. In preferred embodiments of the invention all echoes may be used in a single step when performing the fitting of the magnitude data to the single signal model.

(42) Preferably, the magnitude of the multi-echo MRI data is fitted to the single signal model using an estimation model. In a preferred example of the invention, this uses one or more of the following: regularised least squares estimation; iteratively weighted least squared estimation, m-estimators or s-estimators. Alternatively, other estimation models may be used for fitting the MRI data to the single signal model.

(43) The complex-valued vector Σ.sub.pα.sub.pe.sup.j2πf.sup.f.sup.t.sup.i was calculated a priori and the value was used for all voxels and during optimisation. The three unknown parameters ρ.sub.W, ρ.sub.F and R.sub.2* were estimated from two different sets of initial values. Both first estimates of the signal relaxation parameter were set within the physiologically meaningful range to R.sub.2*=50 s.sup.−1. In an example of the invention, at least one of the relaxation time quantities, T.sub.n*, is used to correct for signal decay, and is preferably in the physically observable range between 1 to 100 ms, still further preferably it is between 20 and 30 ms, when the operational field strength is 1.5 tesla. Alternatively, for a field strength of 3 tesla, the time quantities may be between 10 to 15 ms. If PDFF and R2* are not estimated simultaneously (as they are in the known IDEAL methods), then the PDFF may be inaccurate as it does not correct for R.sub.2* decay, and also R.sub.2* may be inaccurate as it does not correct for the presence of fat. Furthermore, PDFF and R.sub.2* may be less robust and dependent on acquisition parameters, e.g. the chosen number of echoes. Previous studies (31) have shown this effect using simulations and in vivo data.

(44) Initial estimates of water and fat amounts were as follows: {ρ.sub.W, ρ.sub.F}.sub.1={1000, 0} in the first run and the opposite in the second run, {ρ.sub.W, ρ.sub.F}.sub.2={0, 1000}, so PDFF.sub.1=0% and PDFF.sub.2=100%. The scaling of these initial conditions was chosen empirically to account for different scanner gains across acquisition settings. For each run, two sets of solutions {ρ.sub.W, ρ.sub.F, R.sub.2*} were obtained at each voxel as shown in FIG. 6.

(45) The first and second converged solutions had two associated cost function values. Cost function values were in the form of the expression cost function value (ρ.sub.W, ρ.sub.F, R.sub.2*)=Σ.sub.i.sup.N(|ŝ.sub.i|−|s.sub.i|).sup.2, where ŝ.sub.i was the estimated signal using a solution set (ρ.sub.W, ρ.sub.F, R.sub.2*) in Equation 2. Note the definition of the cost function values is also independent of a field map term. The solution set with lowest cost function value was chosen as the “correct” solution at each voxel, whereas the other solution was kept as the alternative solution, and T.sub.2*=1/R.sub.2* maps were calculated.

(46) From the first step, two converged data sets were obtained, where it is presumed that one of the data sets was the true solution. Each of the two solutions had an associated squared 2-norm of the cost function value. The final solution at each voxel was assigned based on those cost function values, where the converged set associated with the lowest cost function value was chosen. The alternative solution was stored. No additional conditions were imposed at this point, (though those versed in the art may imagine many techniques to correlate neighbouring voxels and further improve the fitting at the decision step), similarly to the fieldmap estimation problem.

(47) FIG. 6 shows intermediate and final results of an example of the magnitude only implementation of this invention on a typical case (from Uniklinik Ulm). In each case, the method uses the magnitude of the multi-echo MRI data of images from the subject, where the images are acquired at arbitrarily timed echoes including at least one echo time where water and fat are not substantially in-phase with each other.

(48) The first PDFF map (FIG. 6 a) relates to the first converged solution at each voxel independently and is evidently associated with the lower of the two possible PDFF minima in the general case of an acquisition with an arbitrary echo combination, including at least one echo time where water and fat are not substantially in-phase with each other. This first PDFF map shows the subcutaneous fat region containing PDFF values in the 0-20% range, while PDFF values within the liver are in the plausible range. Note the PDFF map from solution 1 alone is analogous to output PDFF maps from conventional magnitude-based methods. On the other hand, the second PDFF map (FIG. 6b) is associated with the higher of the two possible PDFF minima. This map shows unfeasibly high liver PDFF values, but the subcutaneous fat is estimated correctly. The cost function values for both solutions may be compared and used to choose the most likely correct solution at each voxel independently; this step allows resolving the magnitude-related water-fat ambiguity. It may be the case that two spatially separated voxels with similar PDFF values have substantially different absolute cost function values (see subcutaneous fat in FIG. 6d); it is the relative value of the cost function from solution 1 to solution 2 at the same location that is compared. For this particular acquisition, differences between cost function values may exceed an order of magnitude, so water-fat separation is highly robust and does not require any connectivity analysis. However, in some examples of the invention, the image may be processed, and connectivity analysis, or other analysis may be performed.

(49) The method of this invention was first assessed using a publicly available dataset of twenty-eight phantom acquisitions (http://dx.doi.org/10.5281/zenodo.48266). A phantom comprising a total of eleven vials with peanut oil and water mixtures (PDFF: 0%, 2.6%, 5.3%, 7.9%, 10.5%, 15.7%, 20.9%, 31.2%, 41.3%, 51.4%, 100%) was scanned at different sites (Philips, Siemens and GE Healthcare) using two different six-echo gradient echo protocols at 1.5 T and 3 T: Protocol 1 was a near in-phase/opposed-phase acquisition (TE.sub.1≈ΔTE≈2.30 ms at 1.5 T and TE.sub.1≈ΔTE≈1.00 ms at 3 T) and Protocol 2 aimed for the shortest possible echoes (TE.sub.1=1.10-1.20 ms and ΔTE≈2.00 ms at 1.5 T, and ΔTE≈1.15 ms at 3 T). Acquisitions were designed with a small flip angle (2°-3°) to minimise T.sub.1 bias and combined monopolar and bipolar readouts.

(50) Complex-valued data were available for all acquisitions, but in order to assess the magnitude only method of this invention, the phase information was discarded. The signal model used a six-peak peanut oil fat spectrum corrected for room temperature effects (22° C., relative frequencies in ppm [0.50 −0.49 −2.04 −2.70 −3.50−3.90], relative amplitudes [0.048 0.039 0.004 0.128 0.694 0.087]). 15 mm diameter ROIs were placed manually and were used to extract a median value for each phantom vial from the central slice. Median MRI-PDFF values were plotted against designed phantom concentrations and linear regression was performed for comparison to the available IDEAL results.

(51) FIG. 7 shows the extracted median PDFF values from the ROIs placed on all 11 phantom vials, for all sites, acquisition protocols and field strengths. FIG. 7(a) is the PDFF result according to the prior art method of Hernando, image (b) is using the method of this invention, and image (c) is the difference between the images. Both Protocol 1 and Protocol 2 in the phantom acquisitions used at least one echo time where water and fat were not substantially in-phase with each other: in Protocol 1, these would be any of the opposed-phase echoes; in Protocol 2, these would be most echoes present, as the echoes were not timed according to water and fat resonance frequency shifts.

(52) Linear regression results are shown qualitatively in FIG. 8 and quantitatively in Table 1 in terms of slope, intercept, and R-squared agreement between ground truth values and designed phantom concentrations.

(53) TABLE-US-00001 TABLE 1 MAGO PDFF 1.5T Protocol 1 Hernando PDFF 1.5T Protocol 1 Site R.sup.2 Slope [95% CI] Intercept [95% CI] R.sup.2 [Slope [95% CI] Intercept [95% CI] 1 0.999 1.00 [0.97 1.03]  0.41 [−0.67 1.48] 0.999 1.00 [0.97 1.02]  0.42 [−0.66 1.50] 2 0.997 1.04 [0.99 1.08]  0.81 [−0.93 2.55] 1 1.02 [1.00 1.03]  0.73 [0.15 1.32] 3 0.999 1.01 [0.99 1.03]  0.31 [−0.52 1.14] 0.999 1.01 [0.99 1.04]  0.33 [−0.55 1.22] 4 0.997 0.99 [0.95 1.04] −0.57 [−2.27 1.13] 0.997 0.99 [0.95 1.04] −0.55 [−2.24 1.14] 5 0.999 1.00 [0.99 1.02]  0.18 [−0.50 0.87] 0.999 1.01 [0.99 1.02]  0.21 [−0.46 0.88] 6 0.998 1.00 [0.96 1.03] −0.15 [−1.37 1.06] 0.998 1.00 [0.96 1.03] −0.14 [−1.34 1.07] MAGO PDFF 1.5T Protocol 2 Hernando PDFF 1.5T Protocol 2 Site R.sup.2 Slope [95% CI] Intercept [95% CI] R.sup.2 [Slope [95% CI] Intercept [95% CI] 1 0.998 1.02 [0.99 1.06]  0.08 [−1.25 1.40] 0.998 1.02 [0.99 1.06]  0.09 [−1.21 1.39] 2 0.999 1.02 [1.00 1.04]  0.99 [0.28 1.71] 0.999 1.02 [1.00 1.04]  1.01 [0.29 1.74] 3 0.998 1.01 [0.97 1.04] −0.38 [−1.58 0.82] 0.999 1.00 [0.98 1.03]  0.39 [−1.46 0.68] 4 0.998 0.97 [0.94 1.00]  0.08 [−1.06 1.21] 0.998 0.97 [0.94 1.00]  0.09 [−1.04 1.21] 5 1 0.96 [0.95 0.98]  1.26 [0.68 1.84] 1 0.98 [0.96 0.99]  1.08 [0.67 1.50] 6 0.995 1.01 [0.96 1.06] −0.74 [−2.72 1.25] 0.995 1.01 [0.96 1.06] −0.71 [−2.69 1.26] MAGO PDFF 3T Protocol 1 Hernando PDFF 3T Protocol 1 Site R.sup.2 Slope [95% CI] Intercept [95% CI] R.sup.2 [Slope [95% CI] Intercept [95% CI] 1 0.998 1.00 [0.97 1.03] −0.03 [−1.20 1.14] 0.998 1.00 [0.97 1.03]  0.01 [−1.17 1.15] 2 0.999 1.01 [0.99 1.03]  0.83 [0.11 1.55] 0.999 1.01 [0.99 1.03]  0.85 [0.12 1.58] 3 0.999 1.01 [0.99 1.03]  0.36 [−0.40 1.11] 0.999 1.01 [0.99 1.03]  0.38 [−0.38 1.13] 4 0.997 1.00 [0.96 1.04] −0.13 [−1.71 1.44] 0.997 1.00 [0.96 1.04] −0.12 [−1.68 1.45] 5 0.999 1.00 [0.98 1.01]  0.56 [−0.14 1.26] 0.999 1.00 [0.98 1.01]  0.57 [−0.12 1.26] 6 0.998 0.99 [0.96 1.02] −0.45 [−1.74 0.84] 0.998 0.99 [0.96 1.02] −0.43 [−1.71 0.85] MAGO PDFF 3T Protocol 2 Hernando PDFF 3T Protocol 2 Site R.sup.2 Slope [95% CI] Intercept [95% CI] R.sup.2 [Slope [95% CI] Intercept [95% CI] 1 0.999 0.98 [0.96 1.00]  0.39 [−0.36 1.15] 0.999 0.98 [0.96 1.00]  0.40 [−0.34 1.15] 2 1 0.98 [0.97 0.99]  0.42 [−0.12 0.96] 1 0.98 [0.97 0.99]  0.43 [−0.11 0.97] 3 0.999 0.97 [0.95 1.00]  1.18 [0.12 2.24] 0.999 0.97 [0.95 1.00]  1.18 [0.15 2.21] 4 0.999 0.96 [0.94 0.99]  0.84 [−0.13 1.81] 0.999 0.96 [0.94 0.99]  0.84 [−0.13 1.82] 5 0.999 0.98 [0.96 1.00]  0.69 [−0.15 1.53] 0.999 0.98 [0.96 1.01]  0.61 [−0.25 1.48] 6 0.998 0.97 [0.94 1.00]  0.20 [−1.01 1.41] 0.998 0.97 [0.94 1.00]  0.20 [−1.00 1.41]

(54) Linear regression was performed distinguishing between field strengths (1.5 T and 3 T) and protocols (Protocol 1 and Protocol 2) and averaging across sites. A representative acquisition is shown in FIG. 7 a, with PDFF maps obtained from existing methods and also calculated with the method of this invention, reporting excellent voxel-wise agreement between methods (note the −1 to 1% scaling), and similar noise performance, even though MAGO used half the available information (magnitude data only). Table 1 compares results from the downloaded PDFF maps against the calculated ones and reports both excellent agreement between both methods in-vitro and high accuracy (note high R-squared coefficients, slope close to 1 and intercept close to 0). This encourages the feasibility of the MAGO method for six-echo acquisition protocols.

(55) In this example of the invention the method is able to resolve the magnitude-related water-fat ambiguity after PDFF=50% in all cases, mainly reflected in the 100% phantom vial results, and without compromising accuracy over the 0-50% PDFF range. In general, higher agreement between methods and between MAGO and ground truth values was reported on Protocol 2 data over Protocol 1 data, and also on 3 T data over 1.5 T data.

(56) Research Subjects

(57) In an example of the method of the invention two different sets of test data were used for in vivo testing, each of the test sets having a different purpose. The first test set was a validation set, with more controlled, reliable magnitude and phase data from healthy volunteers to assess the accuracy of the new method against a standard, the complex-based LMS IDEAL, which served as the reference method. The second test set explored the robustness of the two methods to a breadth of realistic acquisition conditions, which may output less reliable data and lead to the presence of artefacts.

(58) An initial test set used in an example of the invention, consisted of single-slice protocol acquisitions from UK Biobank, using a Siemens 1.5 T MAGNETOM Aera, and consisted of N=186 nominally healthy volunteers, with expected low iron and low fat content values. The mentioned protocol had 6 echoes, with echo times={1.2, 3.2, 5.2, 7.2, 9.2, 11.2} ms. 8 cases were discarded due to not having an automatic segmentation mask (which will become relevant in data analysis).

(59) FIG. 9 shows the data analysis pipeline as carried out for the test set. At 101 the magnitude-based algorithm implementation was applied to each case in each data set. At 102 liver segmentation masks were used to draw a median value of PDFF/T.sub.2* within the liver. This is similar to how data analysts draw median values from manually placed ROIs. Reported values were compared to the reference processing method to assess bias and precision of the magnitude-based method (Bland-Altman plot). Then, at 103 automatic liver segmentation masks were used to draw PDFF and T2* histogram-type distributions within the liver mask. In a similar way to how PDFF/T2* values are reported by data analysts (median value of 3 manually placed ROIs), the median values of both PDFF and T2* distributions were reported at 104. Those values were then compared at 105 with the references from LMS IDEAL, also median values drawn from the same liver segmentation.

(60) Results

(61) Validation (Biobank Test Set)

(62) Median values drawn from the magnitude-based PDFF and T.sub.2* output maps were compared to the reference values. FIG. 10 shows absolute comparisons (top images) and Bland-Altman plots (lower images) for both PDFF and T.sub.2* reported values. The 95% confidence intervals (CI) were calculated for the joined analysed cases. For PDFF, 95% CI were (0.052±0.102)%. For T.sub.2*, 95% CI were (−0.053±0.291)%. Difference between complex based LMS IDEAL and the magnitude-based method was less than 0.5% PDFF points in all cases, and less than 1 ms T.sub.2* in all cases.

(63) A typical set of images from the Biobank cohort are shown in FIG. 11, along with a Bland Altman plot for the PDFF measurements. The reference LMS IDEAL method image (left hand image a) is shown with the overlapped automatic segmentation mask. The magnitude-based method of this invention resulted in image (b). The difference image is shown as image (c). illustrating small (<1% PDFF) differences between the complex-based and the magnitude-based methods in the subcutaneous fat regions, demonstrating that the magnitude-based method of this invention was able to accurately estimate PDFF within the full dynamic range (0-100%). Note the smoothness prior on the fieldmap causes LMS IDEAL images to look smoother compared to the magnitude-based method images. The Bland-Altman plot (image 11(d)) comparing all N=178 median PDFF measures shows excellent unbiased agreement between reconstruction methods.

(64) FIG. 12 shows MRI PDFF images and fieldmaps obtained from a Philips Ingenia 3 T, multi-slice case from the Leiden site that presented substantial fat-water swap artefacts in LMS IDEAL caused by incorrect field map estimation, possibly due to unreliable phase information. This figure shows the effect of inconsistent phase information on field map estimation and subsequent PDFF maps calculation for the complex-based method LMS IDEAL.

(65) Field map piece-wise convergence to local minima (row a) causes observable fat-water swap artefacts in the LMS IDEAL PDFF maps that are propagated throughout the images (row b), notably in the liver region and subcutaneous fat, but also in the spleen, the spine and the descending aorta.

(66) PDFF maps from MAGO (the method of this invention) show no evidence of fat-water swaps and are still able to resolve the magnitude-related water-fat ambiguity over the full dynamic range 0-100% (row c). The magnitude only method of this invention (bottom row) looks robust to such errors and was able to resolve the magnitude-related water-fat ambiguity in all slices.

(67) FIG. 13 shows fieldmaps, MRI PDFF images, and distributions obtained from a Siemens TrioTim 3 T slice from the Coimbra site. This presented a fat-water swap artefact in the posterior area of subcutaneous fat in LMS IDEAL. Phase errors also seem to propagate into the liver, shown by the difference image (d) LMS IDEAL-MAGO, but may be imperceptible in many cases. We would expect a normally distributed histogram of PDFF values within the liver mask (FIG. 13 d, e), and a noisy distribution to present higher spread around the median value. While the reported median PDFF values from the automatic whole-liver segmentation agree (11.11% PDFF in LMS IDEAL, 11.07% PDFF in MAGO).

(68) FIG. 13 also shows the effect of field inhomogeneity on the estimation of PDFF maps using LMS IDEAL and the magnitude only method of this invention. The field map image from the complex-based LMS IDEAL method (a) shows aliasing in the posterior area that reflects as a substantial fat-water swap artefact in the LMS IDEAL PDFF map (b), mainly affecting subcutaneous fat but also muscles and the left arm of the individual (right in the image). The PDFF map for the magnitude only method of this invention (c) shows no signs of fat-water swap artefacts and resolves PDFF over the full dynamic range 0-100%. Phase errors at a high field strength are reflected in the PDFF difference image (d) through a speckle pattern and a wider distribution of LMS IDEAL PDFF values within the liver mask (e) compared to the magnitude only distribution (f). A median whole-liver PDFF value is still quite robust to the observed variability in this case, demonstrated by the reported PDFF measures of the two methods (LMS IDEAL PDFF=11.11%, magnitude only method of the invention PDFF=11.07%).

(69) FIG. 14 show the results from a superior slice of a RADIcAL case (Uniklinik, Siemens Skyra 3 T). The MRI PDFF images show the impact of field inhomogeneity and phase errors on the reported PDFF measures by the LMS IDEAL prior art method (images a and b) and the MAGO (magnitude only) method of this invention (image c).

(70) High field map values from the complex-based LMS IDEAL method (a) cause a noisy pattern within the liver region in the LMS IDEAL PDFF map (b), which is less evident from the MAGO PDFF map (c). The method of this invention resolves the PDFF over the full dynamic range 0-100% but the plots (FIGS. 14(e) and (f)) have been rescaled to 0-50% to appreciate the local artefactual texture in the liver, which also shows in the difference image (d) and the histogram distributions of voxel intensities (e, f). The difference image (d) shows agreement where the field map has converged to lower values. It can be seen that the two distributions have the same total count of voxels, since they were drawn from the same automatic segmentation mask. The borders of the automatic segmentation mask were manually corrected in this case. The median whole-liver PDFF values differ substantially (LMS IDEAL PDFF=9.35%, MAGO PDFF=10.97%).

(71) Note the magnitude only method of this invention does not apply smoothing at any stage of the processing, contrarily to smoothness constraints on the LMS IDEAL field map estimation step. In the prior art methods ambiguities may be approached by regularisation over a particular domain, whereas this is not needed for the method of this invention. Of course, this does not mean that spatial regularisation should not be used in association with the MAGO method of the invention; for example, it may be useful in low signal to noise cases. The results presented above for this invention have not used spatial regularisation in order to provide a “base case” of performance. Evidently, a more sophisticated spatial regularisation method could be added in the MAGO PDFF image, based for example on Markov Random Fields.

(72) Generally fat (or other species) may be homogeneously distributed in an organ such as the liver, or may be heterogeneously distributed. If the species distribution is homogeneous, then an MRI PDFF image of a single slice through the organ will be representative for that species throughout the organ. However, if the species is heterogeneously distributed in the organ, then this may have clinical implications, and it is important to obtain accurate MRI PDFF images throughout the organ if the image will be used for clinical or surgical decisions. The method of this invention can be used to provide MRI PDFF images for multiple slices of an organ where the fat (or other species) is heterogeneously distributed.

(73) In summary, the method of this invention has the following features. Firstly, the method only uses magnitude data: this reduces the number of variables to estimate, and the local minima to just two within the physiologically meaningful search space. The method also uses a multi-peak spectral model of at least one species: this enables resolution of the species ambiguity by displacing the aliased solution and reducing its associated cost function value. Also, the method uses a multi-point search step and comparing the cost function values at least two solutions: this enables exploration of both solutions using at least two sets of initial conditions.

(74) MRI CSE methods have become increasingly important clinically for (a) robust water-fat separation—the inclusion of complete multi-peak fat models, compared to conventional Fat-Sat which are only able to target the main fat peak (70% relative amplitude)—, and (b) accurate liver fat fraction quantification for many applications. The non-invasiveness of CSE methods avoids the need for painful expensive biopsies and allow for imaging heterogeneous disease. This invention is for a magnitude-only CSE method which embodies a multi-peak spectrum for fat (or other species) and which can use flexible echo combinations to estimate PDFF across the entire dynamic range (0-100%). Unlike field map estimation, as used in complex-based PDFF estimates where the search algorithm has to contend with multiple local minima, and where an incorrect choice typically gives rise to fat-water swaps, we have shown that in general that this method has to choose between two local minima placed about 50% PDFF. Using 6-point phantom and 6- to 12-point clinical data, we have shown that the “correct” solution can be determined from the cost function values to two runs of the algorithm, for example one starting at 0%, the other at 100%.

(75) As has been shown from the theory and simulated data described above, three necessary conditions are needed for the MAGO method of this invention to work. The first necessary condition implies the use of magnitude data only in order to reduce the multiple local minima that result from complex-based (phase and magnitude) methods to generally two local minima. It has been shown that the fieldmap search space from complex-based methods is not always periodic in the plausible range of fieldmap values (24), so using a multi-point method in complex-based methods may generally be less effective. Furthermore, convergence to wrong fieldmap solutions may not be readily apparent in PDFF maps, as ‘double fat-water swaps’ have been described, where reported PDFF values are incorrect but still in the feasible range. In general, using magnitude data alone ensures only two local minima have to be explored, and the appearance of possible mis-identification is more evident. Also, the use of magnitude data allows for direct estimation of PDFF, without the need of a field map estimation step nor typically used smoothness assumptions (29), which may not always hold. The second necessary condition involves the use of a multi-peak spectral model for one of the species (in this case, fat) in order to break the symmetry in the search space, so the two local minima have different cost function values and the magnitude-related water-fat ambiguity may be resolved (17). The third necessary condition relates to the use of a search space method to explore the two minima; a multipoint search technique was used hereby since previous information on the optimisation space is available: there will be a PDFF local minimum below 50% and another one above 50% in the general case. This ensures correct convergence when actual PDFF values are high and prevents convergence to local minima that has been observed in traditional magnitude-based methods. We note that there are many other potential search space techniques that may be used.

(76) The availability of public phantom data and results allowed for the comparison of the new method to the implementation of the prior art IDEAL method. These experiments also enabled the assessment of the accuracy of the new magnitude-based method directly against ground truth phantom concentrations. The results show comparable accuracy of the MAGO method with respect to the prior art IDEAL method, in terms of slope, intercept, and r-squared agreement, and also showing overall reproducibility for the full dynamic range of PDFF values. The MAGO method of this invention is able to accurately resolve the magnitude-related fat-water ambiguity over the vials at 51.4% PDFF and at 100% PDFF. These reproducibility results suggest the potential of the MAGO method for in-vivo standardization across scanner manufacturers, acquisition protocols and field strengths.

(77) We have noted throughout the experiments that the SNR benefits of high field MRI translate in the method of this invention to more robust and accurate PDFF estimates. This seems to be in contrast to complex-based PDFF estimation methods, whose apparently higher SNR is not found in clinical practice due to phase errors and high field variations. One reason why higher field strength benefits the method of this invention is that it enables more echoes to be used, resulting in higher confidence in the difference in the cost function values for the runs initialised at 0% and 100%. This in turn results in greater resolvability and accuracy of PDFF estimation.

(78) As discussed above, previous methods using phase and magnitude data necessitates fieldmap estimation, and there are many local minima when using complex data and estimating the fieldmap (24), so a multi-point search approach is less feasible. In addition, If the exhaustive search step is not used, the converged solution of water and fat estimates will be dependent on their initial values in the iterative optimisation. In the case PDFF=0 is used as initial conditions, the result will always take PDFF<50% values (15) Finally, using a single-peak model will cause unresolvable ambiguity between the fat and water components.

(79) A new magnitude-only method is presented that shows effectiveness in resolving fat-water ambiguity above 50% PDFF, and its accuracy is validated with multiple Biobank cases against a reference complex-based method, which uses both magnitude and phase information. The new magnitude-only based method is tested against a more challenging cohort, demonstrating similar accuracy and precision to the complex-based in cases where phase information is reliable. Furthermore, the new method presents increased robustness to errors (most of them in the phase images) that often cause complex-based methods in general to fail in clinical practice.

(80) The method of this invention also allows spatial regularisation: generally, the fieldmap, which may be estimated following the estimation of the species (e.g. fat, water) using the method of this application, varies smoothly. The fieldmap can be used to assess image quality, including quantification of any artefacts.

(81) The method of this invention provided MRI species separation. The formulation given in Equation (1) applies to many practical cases, the most of important of which is fat/water (proton density fat fraction estimation). The method described can also be used for artefact detection and estimation, and for estimation of iron content of the liver.

(82) The present invention has been described with reference to the accompanying drawings. However, it will be appreciated that the present invention is not limited to the specific examples herein described and as illustrated in the accompanying drawings. Furthermore, because the illustrated embodiments of the present invention may for the most part, be implemented using electronic components and circuits known to those skilled in the art, details will not be explained in any greater extent than that considered necessary as illustrated above, for the understanding and appreciation of the underlying concepts of the present invention and in order not to obfuscate or distract from the teachings of the present invention.

(83) The invention may be implemented in a computer program for running on a computer system, at least including code portions for performing steps of a method according to the invention when run on a programmable apparatus, such as a computer system or enabling a programmable apparatus to perform functions of a device or system according to the invention.

(84) A computer program is a list of instructions such as a particular application program and/or an operating system. The computer program may for instance include one or more of: a subroutine, a function, a procedure, an object method, an object implementation, an executable application, an applet, a servlet, a source code, an object code, a shared library/dynamic load library and/or other sequence of instructions designed for execution on a computer system.

(85) The computer program may be stored internally on a tangible and non-transitory computer readable storage medium or transmitted to the computer system via a computer readable transmission medium. All or some of the computer program may be provided on computer readable media permanently, removably or remotely coupled to an information processing system.

(86) A computer process typically includes an executing (running) program or portion of a program, current program values and state information, and the resources used by the operating system to manage the execution of the process. An operating system (OS) is the software that manages the sharing of the resources of a computer and provides programmers with an interface used to access those resources. An operating system processes system data and user input, and responds by allocating and managing tasks and internal system resources as a service to users and programs of the system.

(87) The computer system may for instance include at least one processing unit, associated memory and a number of input/output (I/O) devices. When executing the computer program, the computer system processes information according to the computer program and produces resultant output information via I/O devices.

(88) In the foregoing specification, the invention has been described with reference to specific examples of embodiments of the invention. It will, however, be evident that various modifications and changes may be made therein without departing from the scope of the invention as set forth in the appended claims. Those skilled in the art will recognize that the boundaries between logic blocks are merely illustrative and that alternative embodiments may merge logic blocks or circuit elements or impose an alternate decomposition of functionality upon various logic blocks or circuit elements. Thus, it is to be understood that the architectures depicted herein are merely exemplary, and that in fact many other architectures can be implemented which achieve the same functionality.

(89) Any arrangement of components to achieve the same functionality is effectively ‘associated’ such that the desired functionality is achieved. Hence, any two components herein combined to achieve a particular functionality can be seen as ‘associated with’ each other such that the desired functionality is achieved, irrespective of architectures or intermediary components. Likewise, any two components so associated can also be viewed as being ‘operably connected,’ or ‘operably coupled,’ to each other to achieve the desired functionality.

(90) Furthermore, those skilled in the art will recognize that boundaries between the above described operations merely illustrative. The multiple operations may be combined into a single operation, a single operation may be distributed in additional operations and operations may be executed at least partially overlapping in time. Moreover, alternative embodiments may include multiple instances of a particular operation, and the order of operations may be altered in various other embodiments.

(91) However, other modifications, variations and alternatives are also possible. The specifications and drawings are, accordingly, to be regarded in an illustrative rather than in a restrictive sense. Unless stated otherwise, terms such as ‘first’ and ‘second’ are used to arbitrarily distinguish between the elements such terms describe. Thus, these terms are not necessarily intended to indicate temporal or other prioritization of such elements. The mere fact that certain measures are recited in mutually different claims does not indicate that a combination of these measures cannot be used to advantage.

REFERENCES

(92) 1. Younossi Z M, Koenig A B, Abdelatif D, Fazel Y, Henry L, Wymer M. Global epidemiology of nonalcoholic fatty liver disease-Meta-analytic assessment of prevalence, incidence, and outcomes. Hepatology 2016. DOI: 10.1002/hep.28431. 2. Williams C D, Stengel J, Asike M I, Torres D M, Shaw J, Contreras M, Landt C L, Harrison S A. Prevalence of nonalcoholic fatty liver disease and nonalcoholic steatohepatitis among a largely middle-aged population utilizing ultrasound and liver biopsy: a prospective study. Gastroenterology 2011. DOI: 10.1053/j.gastro.2010.09.038. 3. Caussy C, Reeder S B, Sirlin C B, Loomba R. Non-invasive, quantitative assessment of liver fat by MRI-PDFF as an endpoint in NASH trials. Hepatology 2018. DOI: 10.1002/hep.29797. 4. Szczepaniak L S, Nurenberg P, Leonard D, Browning J D, Reingold J S, Grundy S, Hobbs H H, Dobbins R L. Magnetic resonance spectroscopy to measure hepatic triglyceride content: prevalence of hepatic steatosis in the general population. Am J Physiol Endocrinol Metab 2005. DOI: 10.1152/ajpendo.00064.2004. 5. Dixon W. Simple proton spectroscopic imaging. Radiology 1984; 153: 189-194. 6. Glover G. Multipoint Dixon technique for water and fat proton and susceptibility imaging. J Magn Reson Imaging. 1991 September-October; 1(5):521-30. [PMID: 1790376]. 7. Glover G H, Schneider E. Three-point Dixon technique for true water/fat decomposition with B0 inhomogeneity correction. Magn Reson Med. 1991 April; 18(2):371-83. [PMID: 2046518]. 8. Xiang Q S, An L. Water-fat imaging with direct phase encoding. J Magn Reson Imaging. 1997 November-December; 7(6):1002-15. [PMID: 9400843]. 9. Ma J, Singh S K, Kumar A J, Leeds N E, Broemeling L D. Method for efficient fast spin echo Dixon imaging. Magn Reson Med. 2002 December; 48(6):1021-7. DOI: 10.1002/mrm.10306. 10. Reeder S B, Wen Z, Yu H, Pineda A R, Gold G E, Markl M, Pelc N J. Multicoil Dixon chemical species separation with an iterative least-squares estimation method. Magn Reson Med. 2004 January; 51(1):35-45. DOI: 10.1002/mrm.10675. 11. Reeder S B, Pineda A R, Wen Z, Shimakawa A, Yu H, Brittain J H, Gold G E, Beaulieu C H, Pelc N J. Iterative decomposition of water and fat with echo asymmetry and least-squares estimation (IDEAL): application with fast spin-echo imaging. Magn Reson Med. 2005 September; 54(3):636-44. DOI: 10.1002/mrm.20624. 12. Hussain H K, Chenevert T L, Londy F J, Gulani V, Swanson S D, McKenna B J, Appelman H D, Adusumilli S, Greenson J K, Conjee-varam H S. Hepatic fat fraction: M R imaging for quantitative measurement and display—early experience. Radiology. 2005 December; 237(3):1048-55. DOI: 10.1148/radiol.2373041639. 13. Yu H, McKenzie C A, Shimakawa A, Vu A T, Brau A C, Beatty P J, Pineda A R, Brittain J H, Reeder S B. Multiecho reconstruction for simultaneous water-fat decomposition and T2* estimation. J Magn Reson Imaging. 2007 October; 26(4):1153-61. DOI: 10.1002/jmri.21090. 14. Liu C Y, McKenzie C A, Yu H, Brittain J H, Reeder S B. Fat quantification with IDEAL gradient echo imaging: correction of bias from T(1) and noise. Magn Reson Med. 2007. DOI: 10.1002/mrm.21301. 15. Bydder M, Yokoo T, Hamilton G, Middleton M S, Chavez A D, Schwimmer J B, Lavine J E, Sirlin C B. Relaxation effects in the quantification of fat using gradient echo imaging. Magn Reson Imaging. 2008 April; 26(3):347-59. DOI: 10.1016/j.mri.2007.08.012. 16. Yu H, Shimakawa A, McKenzie C A, Brodsky E, Brittain J H, Reeder S B. Multiecho water-fat separation and simultaneous R2* estimation with multifrequency fat spectrum modeling. Magn Reson Med. 2008. DOI: 10.1002/mrm.21737. 17. Yu H, Shimakawa A, Hines C D, McKenzie C A, Hamilton G, Sirlin C B, Brittain J H, Reeder S B. Combination of complex-based and magnitude-based multiecho water-fat separation for accurate quantification of fat-fraction. Magn Reson Med. 2011. DOI: 10.1002/mrm.22840. 18. Cassidy F H, Yokoo T, Aganovic L, Hanna R F, Bydder M, Middleton M S, Hamilton G, Chavez A D, Schwimmer J B, Sirlin C B. Fatty liver disease: M R imaging techniques for the detection and quantification of liver steatosis. Radiographics 2009. DOI: 10.1148/rg.291075123. 19. Bley T A, Wieben O, François C J, Brittain J H, Reeder S B. Fat and water magnetic resonance imaging. J Magn Reson Imaging 2010 DOI: 10.1002/jmri.21895. 20. Haufe W M, Wolfson T, Hooker C A, Hooker J C, Covarrubias Y, Schlein A N, Hamilton G, Middleton M S, Angeles J E, Hernando D, Reeder S B, Schwimmer J B, Sirlin C B. Accuracy of PDFF estimation by magnitude-based and complex-based MRI in children with M R spectroscopy as a reference. J Magn Reson Imaging 2017. DOI: 10.1002/jmri.25699. 21. Noureddin M, Lam J, Peterson M R, Middleton M, Hamilton G, Le T A, Bettencourt R, Changchien C, Brenner D A, Sirlin C, Loomba R. Utility of magnetic resonance imaging versus histology for quantifying changes in liver fat in nonalcoholic fatty liver disease trials. Hepatology 2013. DOI: 10.1002/hep.26455. 22. Tyagi A, Yeganeh O, Levin Y, Hooker J C, Hamilton G C, Wolfson T, Gamst A, Zand A K, Heba E, Loomba R, Schwimmer J, Middleton M S, Sirlin C B. Intra- and inter-examination repeatability of magnetic resonance spectroscopy, magnitude-based MRI, and complex-based MRI for estimation of hepatic proton density fat fraction in overweight and obese children and adults. Abdom Imaging. 2015. DOI: 10.1007/s00261-015-0542-5. 23. Yokoo T, Bydder M, Hamilton G, Middleton M S, Gamst A C, Wolfson T, Hassanein T, Patton H M, Lavine J E, Schwimmer J B, Sirlin C B. Nonalcoholic fatty liver disease: diagnostic and fat-grading accuracy of low-flip-angle multiecho gradient-recalled-echo M R imaging at 1.5 T. Radiology 2009. DOI: 10.1148/radiol.2511080666. 24. Yu H, Reeder S B, Shimakawa A, Brittain J H, Pelc N J. Field map estimation with a region growing scheme for iterative 3-point water-fat decomposition. Magn Reson Med. 2005 October; 54(4):1032-9. DOI: 10.1002/mrm.20654. 25. Hernando D, Haldar J P, Sutton B P, Ma J, Kellman P, Liang Z P. Joint estimation of water/fat images and field inhomogeneity map. Magn Reson Med. 2008 March; 59(3):571-80. DOI: 10.1002/mrm.21522. 26. Hernando D, Kellman P, Haldar J P, Liang Z-P. Robust water/fat separation in the presence of large field inhomogeneities using a graph cut algorithm. Magn Reson Med. 2010 January; 63(1):79-90. DOI: 10.1002/mrm.22177. 27. Lu W, Hargreaves B A. Multiresolution field map estimation using golden section search for water-fat separation. Magn Reson Med. 2008 July; 60(1):236-44. DOI: 10.1002/mrm.21544. 28. Lu W, Lu Y. JIGSAW: Joint inhomogeneity estimation via global segment assembly for water-fat separation. IEEE Trans Med Imaging. 2011 July; 30(7):1417-26. DOI: 10.1109/TMI.2011.2122342. 29. Soliman A S, Yuan J, Vigen K K, White J A, Peters T M, McKenzie C A. Max-IDEAL: a max-flow based approach for IDEAL water/fat separation. Magn Reson Med. 2014 August; 72(2):510-21. DOI: 10.1002/mrm.24923. 30. Sharma S D, Artz N S, Hernando D, Horng D E, Reeder S B. Improving chemical shift encoded water-fat separation using object-based information of the magnetic field inhomogeneity. Magn Reson Med. 2015 February; 73(2):597-604. DOI: 10.1002/mrm.25163. 31. Glover et al “three-point Dixon Technique for True Water/Fat Decomposition with B0 inhomogeneity correction” Magnetic Resonance in Medicine 1991, 18 p 371-383. 32. Yu et al., Robust Multipoint Water-Fat Separation Using Fat Likelihood Analysis 2012, Magnetic Resonance in Medicine 67, 1062-1076. 33. Hernando, Diego, et al. “R2* estimation using “in-phase” echoes in the presence of fat: the effects of complex spectrum of fat.” Journal of magnetic resonance imaging 37.3 (2013): 717-726. 34. Zhong, Xiaodong, et al. “Liver fat quantification using a multi-step adaptive fitting approach with multi-echo GRE imaging.” Magnetic resonance in medicine 72.5 (2014): 1353-1365.