Quadratic species separation using balanced SSFP MRI

Abstract

MRI signals from two species such as water and fat are separated by using two SSFP acquisitions (S1, S2) with the two species signals having two different phase differences (e.g., 180°, 90°) and binomial quadratic expansion in which the species are determined based on magnitude and sign of S1 and the magnitude of S2.

Description

BACKGROUND OF THE INVENTION

This invention relates generally to magnetic resonance imaging (MRI), and more particularly the invention relates to species (e.g. water/fat) signal separation in MRI.

Magnetic resonance (MR) imaging is based on nuclear spins, which can be viewed as vectors in a three-dimensional space. During a MR experiment, each nuclear spin responds to four different effects—precession about the main magnetic field, nutation about an axis perpendicular to the main field, and both transverse and longitudinal relaxation. In steady-state MR experiments, a combination of these effects occurs periodically.

Refocused steady-state free precession (SSFP) sequences have recently gained popularity in magnetic resonance imaging, due to improved gradient hardware. SSFP imaging provides high signal and good contrast in short scan times.

As illustrated in FIGS. 1A, 1B, a refocused SSFP sequence includes a single RF excitation which is repeated periodically. All gradients used for slice selection or imaging are fully rewound over each repetitive time, TR. In the steady-state, the magnetization at points a and d is the same.

Magnetization is tipped about a traverse axis through an angle α. Between excitations, the magnetization undergoes a precession by an angle θ=2πΔfTR about the z-axis (direction of B0), where Δf is the tissue off-resonance, and also experiences both T1 and T2 relaxation.

During the sequence each spin is affected by RF pulses, relaxation and free precession. The steady-state magnetization for SSFP is a function of the sequence parameters flip angle (α), repetition time (TR) and echo time (TE) as well as the tissue parameters T1, T2, and resonant frequency shift Δf.

Balanced steady-state free precession (SSFP) imaging sequences offer relatively high signal and rapid 2D or 3D imaging, potentially overcoming two of the prime weaknesses of MR imaging.

However, clinical applications of balanced SSFP are still limited by two main factors. First, the signal is very sensitive to resonant frequency variations that result from static field inhomogeneity, susceptibility variations or chemical shift. Second, balanced SSFP produces a very bright signal from fat which obscures visualization of other structures.

Numerous methods have been presented to suppress the fat signal including methods that generate a fat suppressed steady state, methods that synthesize a fat-suppressed steady state, and methods with transient fat suppression. Alternatively, several methods separate fat and water using multiple acquisitions, or in a single acquisition.

These methods all have advantages and disadvantages when used with balanced SSFP imaging with regard to efficiency, scan time, robustness and accuracy. Many require multiple acquisitions or increased scan time. As a consequence, most methods compromise the SNR efficiency of balanced SSFP. The only exceptions all can fail when fat and water occupy a single voxel. Thus there is still a trade off that must be made when selecting an imaging sequence for a particular application. Furthermore, most of these techniques still need to be compared quantitatively and in a clinical setting.

Dixon SSFP imaging offers high image quality but requires three or more signal acquisitions with varying echo times and a challenging reconstruction.

There are numerous variations of Dixon imaging. Two applications of Dixon fat/water separation with balanced SSFP are IDEAL (iterative decomposition of water and fat using echo asymmetry and least squares estimation) which uses an iterative reconstruction that maintains efficiency by permitting the use of small echo spacings and the in-phase/out-of-phase 2-point approach that is essentially a variation of the linear combination technique. The former of these techniques provides accurate quantitation of fat and water within a voxel. The latter technique requires only 2 acquisitions, but compromises accuracy. A final “Dixon” technique that is of interest here is partially-opposed-phase (POP), which acquires 2 echoes with angular separation 0° and approximately 135°. The present invention here is similar, but with improved efficiency and robustness when applied to balanced SSFP imaging, where a refocused opposed-phase acquisition is possible.

SUMMARY OF THE INVENTION

The present invention provides an efficient method for species signal separation that is more accurate than phase-sensitive SSFP, yet significantly faster and simpler than 3-point IDEAL imaging. The technique can provide water/fat separation in balanced SSFP in two acquisitions.

In its simplest embodiment, this technique separates water and fat using two balanced SSFP acquisitions S₁ and S₂, that have water and fat 180° and 90° out of phase respectively. After phase correction of S₁ the water (W) and fat (F) signals can be determined by simply solving a quadratic equation that depends on the magnitude and sign of S₁ and the magnitude of S₂.

A second embodiment combines this technique with multiple acquisition SSFP to separate water and fat in the presence of off resonance using a total of 4 acquisitions.

The invention and object and features thereof will be more readily apparent from the following detailed description and appended claims when taken with the drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIGS. 1A, 1B illustrate a SSFP phase sequence in accordance with prior art.

FIG. 2 illustrates magnitude (top) and phase (bottom) of a balanced SSFP signal with TE=TR/2=2.3 ms for different tissues.

FIG. 3 is a geometric interpretation of signals S₁ and S₂ in a fat (F) and water (W) coordinate plane.

FIG. 4 illustrates mean and standard deviation of a water estimation error (Ŵ−S_W) as a function of water fraction for different SNR values.

FIGS. 5( a)–5(h) illustrate water/fat separation without the use of a magnetic field map.

FIG. 6 illustrates fat/water separation through a tilted projection.

FIGS. 7( a)–7(c) illustrate sagittal images of a knee showing fat (a), water (b), and recombined images (c).

FIGS. 8( a)–8(c) illustrate coronal images of a knee showing fat (a), Water (b) and recombined images (c).

FIGS. 9( a)–9(d) illustrate sagittal knee images for TR=2TE=13.8 ms.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

In balanced SSFP, the signal from spins in each band of frequencies of width 1/TR is refocused at TE=TR/2, to form a spin echo. The direction of the spin echo alternates by 180° with each frequency increase of 1/TR. This behavior is well documented and is shown in FIG. 2. If TE is not TR/2, then the signal phase is simply the phase at TR/2 plus precession between TE and TR/2.

The selection of

$\mathrm{TR}=\frac{2k+1}{\Delta f_{\mathrm{cs}}},$ where k is an integer and Δf_cs is the chemical shift frequency difference between fat and water, will guarantee that fat and water are always refocused 180° out of phase. (Selection of a TR close to such a value is often sufficient).

Two balanced SSFP images are acquired with TR chosen above, and all parameters kept constant with the exception of TE. For the first image, TE=TE₁=TR/2. For the second image,

$\mathrm{TE}={\mathrm{TE}}_2=\mathrm{TR}/2+\frac{v}{2\pi }/\Delta f_{\mathrm{cs}}.$ Denote the complex signals for the two images as S₁ and S₂ whereS₁=(W−F)e^iφ  (1)S₂=(W−e^iνF)e^iφe^{2πi(TE2−TE1)f}  (2)Here W and F are the signal magnitudes from water and fat respectively. φ is a constant signal phase term resulting from coil phase, slice profile or timing mismatches, and varies slowly in space. f is the frequency corresponding to field inhomogeneity and susceptibility variations, and can vary arbitrarily in space.

Complex images of S₁ and S₂ are reconstructed using a standard Fourier image reconstruction. First, an estimate, {circumflex over (φ)}, of the value of φ at each voxel is determined. This can be done by assuming that φ varies slowly in space, and that the sign of W−F can change arbitrarily. Beginning at a voxel known to be water or fat, the phase {circumflex over (φ)} of each voxel or block of voxels moving outwards is found, such that {circumflex over (φ)} differs by less than 90° from that of neighboring blocks or voxels. This scheme has been used for phase-sensitive SSFP and similar schemes have been presented for other Dixon reconstruction approaches.

Given the phase estimate, {circumflex over (φ)}, define the water-fat difference estimate as{circumflex over (d)}=S₁e^{−i{circumflex over (φ)}}  (3)Next, using {circumflex over (d)} and the fact that |S₂|²=W²+F²+2WF cos ν, use Eq. 3 to write|S₂|²≈({circumflex over (d)}+F)²+F²+2({circumflex over (d)}+F)F cos ν  (4)This is a simple quadratic equation in F. The desired solution of Eq. 4 gives the estimates of W and F:

$\begin{array}{cc}\hat{F}=-\frac{1}{2}\hat{d}\frac{1}{2}\sqrt{\widehat{d^2}+2\frac{{S_2}^2-{\hat{d}}^2}{1+\cos v}}& \left(5\right)\end{array}$ andŴ={circumflex over (d)}+{circumflex over (F)}  (6)The physical meaning of the lesser root of Eq. 4 is to give {circumflex over (F)} as an estimate of −W and Ŵ as an estimate of −F.

It is useful to note that in the case where ν=π/2, Eq. 5 is simply

$\begin{array}{cc}\hat{F}=-\frac{1}{2}\hat{d}\frac{1}{2}\sqrt{2{S_2}^2-{\hat{d}}^2}& \left(7\right)\end{array}$ A final, optional step is to calculate an estimate of the field inhomogeneity f (actually in Hz) as:

$\begin{array}{cc}f=\frac{1}{2\pi \left({\mathrm{TE}}_2-\mathrm{TR}/2\right)}\angle \left(S_2{e}^{-i\hat{\varphi }}\left(\hat{W}-{\mathrm{ie}}^{\mathrm{iv}}\hat{F}\right)\right)& \left(8\right)\end{array}$ However, note that this calculation is not necessary for determination of Ŵ and {circumflex over (F)}.

FIG. 3 shows a geometric interpretation of the reconstruction of F and W from signals S₁ and S₂ in a FW coordinate plane, for ν=π/2. The point {circumflex over (d)} is located on the W axis after removing phase {circumflex over (φ)} from S₁, and a line of slope 1 is drawn through this point (thin solid line). A circle centered at the origin, with radius |S₂| (dotted line) is drawn. The intersection of the line with the circle in quadrant 1, shown by the solid dot gives the water (Ŵ) and fat ({circumflex over (F)}) values, shown by hollow dots. Notice that the other intersection of the line and circle gives negative, interchanged estimates of fat and water.

It is also useful to consider that the water and fat estimates are most accurate when the line and circle intersect at 90°, when F=W. Noise in |S₂|² will enlarge/reduce the size of the circle, while noise in {circumflex over (d)} will shift the line horizontally.

FIG. 3 shows the solution to Eq. 5. The fact that the dependence on S₂ gives a circular shape is not because the phase of S₂ is unimportant in the fit. When general values of ν (other than π/2) are considered the circle becomes an ellipse that intersects the W and F axes at |S₂|. Extreme cases are ν=0 and ν=π. As νapproaches 0, the ellipse approaches a line with slope −1, representing the very simple case where {circumflex over (F)}={circumflex over (d)}/2+|S₂|/2. As ν approaches π, S₁ and S₂ become the same, and the ellipse is asymptotic to two lines of slope 1 intersecting the W axis and F axis respectively at |S₂|.

Since this method is for use with balanced SSFP imaging, it is useful to consider the effect of crossing a signal null, such as at 110 Hz in FIG. 2. For the purposes of this discussion, assume that TR=4.6 ms.

First assume the voxel is entirely water. The estimate, {circumflex over (φ)} is intentionally insensitive to phase changes of 180°, so is unaffected by the null-crossing. It follows that the sign of {circumflex over (d)} is now incorrect, so fat and water will be interchanged. Thus {circumflex over (F)}=W and Ŵ=F. Finally, consider the “water frequency” f as determined by Eq. 8. The angle of S₂ changes by 180°. Furthermore, because water and fat are interchanged, the angle of Ŵ−ie^iν{circumflex over (F)} decreases by 90°. The net result is a 90° increase in the phase, or a 220 Hz increase in calculated f.

Second, assume the voxel is entirely water. Again, {circumflex over (φ)} is unaffected by the null-crossing and fat and water will be interchanged. The angle of S₂ changes by 180°, and the angle of Ŵ−ie^iν{circumflex over (F)} increases by 90°. The net result is a 90° decrease in the phase, or a 220 Hz decrease in calculated f.

Third, assume the voxel is split between water and fat. In this case the angle of S₂ changes by 180°, while the angle of Ŵ−ie^iν{circumflex over (F)} is unchanged. Thus f changes by 440 Hz (±440 Hz are not distinguishable).

It follows that the overall frequency change is greater than 220 Hz. By simply adding 180° to {circumflex over (φ)}, fat and water will be correctly identified, and the angle change in S₂ will be cancelled.

The fact that water and fat are opposed-phase in S₁ is essential to this technique, and in absence of noise allows perfect reconstruction for any value of ν except ν=π. If, however, S₁ was acquired with water and fat in phase, the quadratic separation is not possible. This can be seen by replacing the line in FIG. 3 with a line of slope −1, and noting that this line now intersects the circle at two points in the first quadrant. Derivation of Eq. 5 for this case would show that there is an ambiguity in taking the square root, and which root to use. For the opposed-phase case, it is clear that the greater quadratic root is the appropriate solution. As previously stated, the lesser root produces negative, interchanged values of Ŵ and {circumflex over (F)}.

Using Matlab 6.5 (Natick, Mass.), a function calculates S₁(W, F, φ, f), and S₁(W, F, φ, f), given sequence parameters α (flip angle), TR, TE₁, TE₂, and the fat/water chemical shift frequency f_cs.

N “pixels” were simulated, with the values of W and F constant such that W+F=1. φ was chosen from a uniform distribution between −π/4 and π/4, and f was normally distributed with a mean of 0 Hz and standard deviation of 100 Hz.

For each pixel, S₁ and S₂ were calculated for TR=4.6 ms, TE₁=2.3 ms, TE₂=1.15 ms, α=60° and f_cs=220 Hz. For these parameters, ν=π/2. Next, four random variables, X₁, X₂, Y₁ and Y₂ were selected from a Gaussian distribution with a mean of 0 and standard deviation σ for each voxel. To simulate noise, signals {tilde over (S)}₁ and {tilde over (S)}₂ were generated as:{tilde over (S)}₁=S₁+X₁+iY₁  (9){tilde over (S)}₂=S₂+X₂+iY₂  (10)To simulate the reconstruction, {circumflex over (φ)} was estimated, assuming a small distribution, from the angle of {tilde over (S)}₁. With this small distribution of φ, the simulation does not really address the detection of φ.

Equations 3, 5 and 6 were used to determine the values of Ŵ and {circumflex over (F)} from {tilde over (S)}₁ and {tilde over (S)}₂. The mean and standard deviation of the detection errors Ŵ−S_w and {circumflex over (F)}−S_f were calculated. Here S_w is the average of |S₁| and |S₂| when W=1 and F=0, i.e. the balanced SSFP signal for an all-water voxel. S_f is, similarly the average of |S₁| and |S₂| when W=0 and F=1.

These simulations were repeated for numerous parameter combinations with 0<W<1 and SNR of 5, 15 and 50. FIG. 4 shows the mean and standard deviation of the error in water estimate (Ŵ−S_w) for SNR values of 5, 15 and 50.

Note the trend with mean—due to the fact that the radical in Eq. 5 is forced to be positive. Also, the standard deviation of the estimate is largest when the water fraction is low. This also makes sense because |S₂| is not very sensitive to W when W/F is small (as can be seen from FIG. 3) or the fact that |S₂|²=W²+F²).

The experiments to test this method include validation in a fat/water phantom, and in vivo imaging.

We verified this technique using a bottle phantom filled with tap water and vegetable oil, laid on its side. We used a custom 3D balanced SSFP pulse sequence with TR=5.8 ms, TE₁=2.9 ms, TE₂=1.75 ms, and α=60°. Other scan parameters were a 24×24×9.6 cm³ field of view (FOV) with 256×160×48 matrix, leading to 1.0×1.5×2.0 mm³ voxel size. A small linear gradient shim, approximately 0.038 mT/m, was applied along the bottle to create a linear frequency variation of 16 Hz/cm.

All scans were performed on a 1.5 T GE Signa Excite scanner with 40 mT/m gradients, 150 mT/m/ms gradient slew rates and a quadrature transmit/receive head coil (GE Healthcare, Waukesha, Wis.).

Both acquisitions were reconstructed using a standard Fourier reconstruction to produce images S₁ and S₂. {circumflex over (φ)} was estimated using a region-growing algorithm, and Ŵ and {circumflex over (F)} were calculated using Equations 3, 5 and 6 were used to determine the water and fat content.

FIG. 5 shows the magnitude and phase images of S₁ and S₂ for the fat/water phantom, followed by magnitude and phase of {circumflex over (d)}, and finally the magnitudes of Ŵ and {circumflex over (F)}. Images (g) and (h) show that fat and water are correctly separated in the central part of the bottle. In other regions of the bottle, fat and water are incorrectly determined. This is because of the off-resonance sensitivity of balanced SSFP, and emphasizes the need for a good shim and a short TR.

Projection in vertical direction of slanted block across interface. As a preliminary test of the ability of this technique to quantify fat and water content in a single voxel, vertical projections were made through a slanted slab in the images S₁ and S₂. The slab was chosen to intersect the fat/water boundary in the bottle, so that the water/fat ratio varies with position in the slab. The projected voxels from S₁ and S₂ were used to calculate Ŵ and {circumflex over (F)} as a function of position along the slab using the same reconstruction as for the 3D images.

FIG. 6 shows a tilted slab that is used to calculate projected voxels from S₁ and S₂. The fat and water content in each projected voxel is plotted as a function of position in the lower part of the figure. This is a preliminary illustration of the ability of this algorithm to separate fat and water within a voxel regardless of the water/fat ratio in the voxel.

This technique was tested in vivo by acquisition of 3D balanced SSFP images in the knee of a normal volunteer, using identical parameters as for the phantom scan, except that the matrix size was 160×160×96 matrix, leading to a 1.5×1.5×2.0 mm³ voxel size. No additional shim was applied for this scan, as no field variation is desired.

The scans were performed on a 1.5 T GE Signa Excite scanner with 40 mT/m gradients, 150 mT/m/ms gradient slew rates and a linear transmit/receive extremity coil (GE Healthcare, Waukesha, Wis.). Reconstruction procedure was identical to that used for the phantom scans.

FIGS. 7 and 8 show separated fat and water images and the recombined ({circumflex over (F)}+Ŵ) image in the sagittal and coronal planes. Fat and water are well separated, and the presence of a small amount of water in the femur and tibia shows that the technique is able to process voxels containing both fat and water.

As mentioned above, the field map estimate {circumflex over (f)} is not needed for water/fat separation. However, using Eq. 8, {circumflex over (f)} can be determined, and used to detect pixels where crossing of an SSFP signal null has caused an incorrect sign of {circumflex over (d)}. Negating {circumflex over (d)} and adding 180° to {circumflex over (φ)} at these points, gives the correct estimates of Ŵ, {circumflex over (F)} and {circumflex over (f)}.

This is of particular significance for high-field SSFP or high-resolution SSFP. In these cases the TR is usually too long to avoid SSFP signal nulls, so multiple acquisition methods can be used. The following experiment combines the proposed quadrature separation with the “maximum pixel” multiple acquisition method.

A balanced SSFP protocol with TR=2TE₁=13.8 ms and TE₂=6.75 ms, and a 30° flip angle was repeated twice, first with RF phase alternating by 180°, then with constant RF phase for a total of 4 image acquisitions. The total scan time was 15:05 for a 256×256×64 matrix over a 16×16×12.8 cm³ FOV. Ŵ, {circumflex over (F)} and {circumflex over (f)} were calculated as above for both phase cycles. At points where {circumflex over (f)} deviated by more than an empirically determined threshold of 140 Hz, {circumflex over (d)} was negated and Ŵ and {circumflex over (F)} were recalculated. Finally, at each pixel, the maximum from the two phase cycles was taken for both Ŵ and {circumflex over (F)} images.

Again, all scans were performed on a 1.5 T GE Signa Excite scanner with 40 mT/m gradients, 150 mT/m/ms gradient slew rates and a linear transmit/receive extremity coil (GE Healthcare, Waukesha, Wis.). Reconstruction procedure was identical to that used for the phantom scans.

FIG. 9 shows the original S₂ images with RF phase increments of 0° and 180° to show the SSFP signal nulls that result from the long TR of 13.8 ms. The separated water and fat images following separation and maximum-pixel combination are shown below in the figure. The quadrature separation technique correctly determines the water and fat images, while the multiple-phase-cycle acquisitions allow removal or reduction of the signal nulls.

The above experiments are proof-of-concept that quadratic water/fat separation with balanced SSFP is feasible for in vivo MR imaging. This section discusses some additional areas that will be interesting to explore.

For full-readout Cartesian imaging, selection of TE=TE₁=TR/2 is very efficient. However, keeping TR constant with a full readout becomes less efficient as the difference between TE₂ and TE₁ increases. Although this suggests making TE₂ similar to TE₁, this will decrease robustness. As indicated by Eqs. 3, 6 and 7, the reconstruction is quite simple when TE₂ are chosen so that ν=π/2. Examining Eq. 5, when π/2<ν<π, the rightmost term under the radical can become large, corresponding to a poorly-conditioned system. Conversely, reduction of ν toward 0 will improve robustness over the ν=π/2 case. Thus the choice of echo times is ultimately a tradeoff between reconstruction robustness and scan efficiency.

As was shown by FIG. 4 the propagation of noise in this reconstruction is not simple. First, the propagation of noise is improved when phase-correcting S₁, assuming the correction is accurate, since the correction removes any noise that is in quadrature with the signal. Second, using the magnitude of S₂ changes the noise statistics from Rayleigh to Rician. Third, examining FIG. 3, the quadratic solution for {circumflex over (F)} is less sensitive to noise when the line and circle intersect at right angles. Thus, the effect of noise will depend on the water/fat ratio. Finally, the selection of TE₂, and thus ν will change the circle in FIG. 3 to an ellipse, and affect this intersection angle. While a full analysis of these noise effects is beyond the scope of the current description, it is clear from simulations that the SNR is somewhere between that expected for a single acquisition and two acquisitions.

The use of phased-array coils for SNR improvement is of significance in MR. For a multicoil non-accelerated reconstruction, this method can generate phase maps for individual coils at the phase-correction stage. Once the appropriate phase maps have been generated, other steps of the reconstruction can be performed individually for each coil, and final images combined.

Parallel imaging techniques use individual coil sensitivity information to accelerate imaging. We have previously shown that phase-sensitive SSFP can be implemented following a sensitivity-encoded (SENSE) reconstruction that preserves image phase (See Vasanawala S S, et al. “Phase-sensitive SSFP parallel imaging,” Proceedings of the 12th Annual Meeting of ISMRM, Kyoto, 2004. p. 2253.). A similar reconstruction applied to individual coil signals can produce the S₁ and S₂ signals needed for quadrature reconstruction. Thus this technique will be compatible with accelerated imaging methods.

Partial k-space acquisitions can reduce scan times by close to 50%, but introduce an unwanted quatrature signal component that can be suppressed in some cases using homodyne reconstruction or other techniques. Because water and fat are 180° out of phase in S₁, partial k-space reconstruction can, in principle be performed perfectly, as the signal introduced by partial k-space acquisition is in quadrature with both water and fat. However, the same cannot be said for S₂. In fact, even with perfect field homogeneity, a partial k-space acquisition of S₂ cannot be properly reconstructed if water and fat are assumed to vary with high spatial frequency. Some advanced techniques or assumptions of slow variation between water and fat may permit a partial acquisition of S₂. However, with existing reconstruction methods, a reduction of only 25% in overall scan time is achievable using a partial k-space acquisition.

Quadratic water/fat separation is a promising new method for use with balanced SSFP imaging that exploits the refocusing properties of balanced SSFP, which previous techniques do not. Although similar to Dixon techniques, it requires only two, rather than three, acquisitions and does not require a field map to quantify water and fat.

While the invention has been described with reference to specific embodiments, the description is illustrative of the invention and is not to be construed as limiting the invention. Various applications may occur to those skilled in the art without departing from the spirit and scope of the invention as defined by the appended claims.

Claims

1. A method of separating MRI signals (Sp1, Sp2) from two species in an object being imaged comprising the steps of: a) acquiring MRI signals (S1, S2) using two balanced SSFP acquisitions in which the two species are at respective phase differences, φ and ν, and in which echo times TE1 and TE2 differ, with S1=(Sp1−Sp2)eiφS2=(Sp1−eiνSp2)eiφe2πi(TE2−TE1)f where f is a frequency corresponding to field inhomogeniety and susceptibility variations,b) estimating a phase, {circumflex over (φ)}, for voxels in a field of view, and defining species differences as {circumflex over (d)}=Sp1−Sp2=S1e−1{circumflex over (φ)}, andc) using a quadratic expression |S2|2=Sp12+Sp22+2Sp1·Sp2 cos ν

2. The method of claim 1 wherein the two species are water and fat with signal magnitude W and F, respectively.

3. The method of claim 2 where first echo time, TE1, is one-half repetition time, TR/2, and the second echo time,

4. The method of claim 3 wherein phase difference between fat and water in S1 is 180° and in S2 is 90°.

5. The method of claim 4 wherein

6. The method of claim 5 where a field map is calculated after estimating the two species.

7. The method of claim 6 where species and field map estimates are calculated using reduced-resolution images to reduce imaging time.

8. The method of claim 1 wherein S1 and S2 are acquired multiple times with different phase increments in SSFP to reduce the effects of dark bands.

9. The method of claim 1 wherein multiple coils are used, and the coil sensitivities are used to reconstruct signals S1 and S2 for different coils or combinations of coils.

10. The method of claim 1 where signals S1 and S2 are both acquired in a “multi-echo” sequence, using a single RF excitation.

11. The method of claim 1 where S1 is acquired using a pulse sequence that results in water and fat being opposed-phase.

12. The method of claim 1 where TE1=(2k+1)/(2fcs) where k is an integer and fcs is the frequency difference between species.