Mapping and Correction of Inhomogeneity in Magnetic Resonance Imaging Magnetic Field

Abstract

A system and method of mapping and correcting the inhomogeneity of a magnetic field within an object using an Magnetic Resonance Imaging (MRI) system where there is a single dominant resonance. The method includes acquiring at least three MRI images, each at different echo times (TE). At least two ΔTE images (ΔTE.sub.i=1 . . . N) are generated based on the at least three MRI images, wherein the subscripts I=1 N refer to images with sequentially increasing ΔTE times. Aliasing in the ΔTE.sub.1 image is permitted. The ΔTE times of ΔTE.sub.1 and ΔTE.sub.2 are set such that the alias points at which wrapping occurs in ΔTE.sub.1 does not overlap with the alias points of ΔTE.sub.2. Each ΔTE image is unwrapped. A final B0 map is set to the unwrapped ΔTE.sub.N image.

Claims

1. A method for mapping and correcting the inhomogeneity of a magnetic field within an object using an Magnetic Resonance Imaging (MRI) system where there is a single dominant resonance, the method comprising: acquiring at least three MRI images, each at different echo times (TE); generating at least two ΔTE images (ΔTE.sub.i=1 . . . N) based on the at least three MRI images, wherein the subscripts i=1 . . . N refer to images with sequentially increasing ΔTE times; wherein aliasing in the ΔTE.sub.1 image is permitted, and wherein the ΔTE times of ΔTE.sub.1 and ΔTE.sub.2 are set such that the alias points at which wrapping occurs in ΔTE.sub.1 does not overlap with the alias points of ΔTE.sub.2; and unwrapping each ΔTE image; and setting a final B.sub.0 map to the unwrapped ΔTE.sub.N image.

2. The method according to claim 1, further comprising: determining shim currents and/or passive shims based on the final B.sub.0 map, so as to reduce B.sub.0 inhomogeneity.

3. The method according to claim 2, further comprising applying the shim currents and/or passive shims to the sample to reduce B.sub.0 inhomogeneity in the sample.

4. The method according to claim 2, wherein the shim currents and/or passive shims generate a B.sub.0 field opposite in sign to the inhomogeneity of the magnetic field.

5. The method according to claim 1, wherein the at least two ΔTE images differ in their echo time (TE) such that the difference in resulting phase (Δϕ) of any two ΔTE images varies according to Δϕ=ΔTE.Math.Δv, where Δv is the difference in frequency at a voxel in the image from a central frequency, and ΔTE is the difference in echo time between any pair of images.

6. The method according to claim 1, wherein unwrapping each ΔTE image to form a final B.sub.0 map includes: unwrapping aliasing in the ΔTE.sub.1 image is based on information in the next longest ΔTE image (the ΔTE.sub.2 image), to form an unwrapped ΔTE.sub.1 image, and using the unwrapped ΔTE.sub.1 image to unwrap the next longest ΔTE image (the ΔTE.sub.2 image).

7. The method according to claim 6, wherein unwrapping aliasing in the ΔTE.sub.1 and ΔTE.sub.2 images includes calculating a logistic temporal unwrapping parameter LTUP, where LTUP images the fractional difference measured between two ΔTE images into a set of values describing the extent of aliasing.

8. The method according to claim 6, wherein unwrapping each ΔTE image to form a final B.sub.0 map further includes: for i=3 to N, using the previously unwrapped ΔTE.sub.i-1 image to unwrap the ΔTE.sub.i image.

9. The method according to claim 6, wherein using the unwrapped ΔTE.sub.1 image to unwrap the next longest ΔTE image (the ΔTE.sub.2 image) includes using a numerical estimate of the difference between the unwrapped value in the ΔTE.sub.1 image with values representing different degrees of aliasing in the ΔTE.sub.2.

10. The method according to claim 1, wherein each of the MRI images is acquired using gradient echo, asymmetric spin echo, and/or stimulated echo acquisitions.

11. The method according to claim 1, wherein unwrapping each ΔTE image includes a single pattern matching algorithm to simultaneously unwrap each ΔTE images.

12. The method according to claim 1, wherein acquiring the at least three MRI images includes using simultaneous multi-slice methods, K-space undersampling in the phase encoding direction, and/or use of 2 and 3 dimensional readout trajectories, to decrease acquisition times.

13. A system for mapping and correcting the inhomogeneity of the magnetic field within an object using Magnetic Resonance Imaging (MRI) where there is a single dominant resonance, the system comprising: an MRI imaging device for acquiring at least three MRI images, each at different echo times (TE), wherein the MRI imaging device includes a controller configured to: generate at least two ΔTE images (ΔTE.sub.i=1 . . . N) based on the at least three MRI images, wherein the subscripts I=1 . . . N refer to images with sequentially increasing ΔTE times; wherein aliasing in the ΔTE.sub.1 image is permitted, and wherein the ΔTE times of ΔTE.sub.1 and ΔTE.sub.2 are set such that the alias points at which wrapping occurs in ΔTE.sub.1 does not overlap with the alias points of ΔTE.sub.2, unwrap each ΔTE image; and set a final B.sub.0 map to the unwrapped ΔTE.sub.N image.

14. The system according to claim 13, wherein the controller is further configured to determine shim currents and/or passive shims based on the final B.sub.0 map, so as to reduce B.sub.0 inhomogeneity.

15. The system according to claim 14, wherein the controller is further configured to apply the shim currents and/or passive shims to the sample to reduce B.sub.0 inhomogeneity in the sample.

16. The system according to claim 14, wherein the shim currents and/or passive shims generate a B.sub.0 field opposite in sign to the inhomogeneity of the magnetic field.

17. The system according to claim 13, wherein the at least two ΔTE images differ in their echo time (TE) such that the difference in resulting phase (Δϕ) of any two ΔTE images varies according to Δϕ=ΔTE.Math.Δv, where Δv is the difference in frequency at a voxel in the image from a central frequency, and ΔTE is the difference in echo time between any pair of images.

18. The system according to claim 13, wherein to unwrap each ΔTE image to form a final B.sub.0 map, the controller is configured to: unwrap aliasing in the ΔTE.sub.1 image is based on information in the next longest ΔTE image (the ΔTE.sub.2 image), to form an unwrapped ΔTE.sub.1 image, and use the unwrapped ΔTE.sub.1 image to unwrap the next longest ΔTE image (the ΔTE.sub.2 image).

19. The system according to claim 18, wherein to unwrap aliasing in the ΔTE.sub.1 and ΔTE.sub.2 images, the controller is configured to include calculating a logistic temporal unwrapping parameter LTUP, where LTUP images the fractional difference measured between two ΔTE images into a set of values describing the extent of aliasing.

20. The system according to claim 18, wherein to unwrap each ΔTE image to form a final B.sub.0 map, the controller is configured to: for i=3 to N, use the previously unwrapped ΔTE.sub.i-1 image to unwrap the ΔTE.sub.i image.

21. The system according to claim 18, wherein to use the unwrapped ΔTE.sub.1 image to unwrap the next longest ΔTE image (the ΔTE.sub.2 image), the controller is further configured to use a numerical estimate of the difference between the unwrapped value in the ΔTE.sub.1 image with values representing different degrees of aliasing in the ΔTE.sub.2.

22. The system according to claim 13 wherein the MRI device acquired each of the MRI images using gradient echo, asymmetric spin echo, and/or stimulated echo acquisitions.

23. The system according to claim 13 wherein to unwrap each ΔTE image, the controller is configured to use a single pattern matching algorithm to simultaneously unwrap each ΔTE images.

24. The system according to claim 13, wherein the MRI device is configured to acquire the at least three MRI images using simultaneous multi-slice methods, K-space undersampling in the phase encoding direction, and/or use of 2 and 3 dimensional readout trajectories, so as to decrease acquisition times.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0017] The foregoing features of embodiments will be more readily understood by reference to the following detailed description, taken with reference to the accompanying drawings, in which:

[0018] FIG. 1 shows a schematic of a Magnetic Resonance Imaging (MRI) system (prior art);

[0019] FIG. 2 shows a conventional gradient echo pulse sequence for B.sub.0 mapping;

[0020] FIGS. 3A and 3B show Aliasing and Uncertainty in B.sub.0 Mapping;

[0021] FIG. 4 shows an example of logistical temporal unwrapping, in accordance with embodiment of the invention;

[0022] FIG. 5 shows the effects of uncertainty on LTUP in the ΔTE=2 ms image, in accordance with an embodiment of the invention;

[0023] FIG. 6 shows the effects of uncertainty in the ΔTE=2 ms image and the ΔTE=9 ms image, in accordance with an embodiment of the invention;

[0024] FIG. 7 compares images and B.sub.0 maps acquired from a healthy subject using conventional methodology; and using the Logistical Temporal Unwrapping methodology in accordance with embodiments of the invention;

[0025] FIG. 8 is an illustrative flow diagram of a Logistical Temporal Unwrapping methodology, in accordance with an embodiment of the invention; and

[0026] FIG. 9 shows three wrapped ΔTE maps and the final unwrapped B.sub.0 map obtained from a spherical phantom, in accordance with various embodiments of the invention.

DETAILED DESCRIPTION OF SPECIFIC EMBODIMENTS

[0027] Definitions. As used in this description and the accompanying claims, the following terms shall have the meanings indicated, unless the context otherwise requires:

[0028] The term “image” refers to a dataset in which a parameter(s) such as, without limitation, phase, B.sub.0 field variation etc. . . . , is provided as a two or three dimensional map reflecting the spatial variation of that parameter.

[0029] In illustrative embodiments of the invention, a system and method is provided for mapping and correcting the inhomogeneity of a magnetic field within an object using Magnetic Resonance Imaging (MRI) or Nuclear Magnetic Resonance (NMR) systems. The system and method includes logistical temporal unwrapping that enables the use of longer ΔTE images to unwrap aliasing in the shortest ΔTE image, and utilizes the patterns of aliasing in all images to unwrap the most accurate longer ΔTE image. Aliasing in the shortest ΔTE image is advantageously allowed, minimizing the need to accurately predict the maximum inhomogeneity present. Additionally, the ΔTE values used can lengthen by factors larger than two, enabling greater accuracy with fewer ΔTE images and thus shorter acquisition times. Details are described below.

[0030] FIG. 1 shows a conventional Magnetic Resonance Imaging (MRI) system. The MRI system includes a magnet, gradient coil (used for spatial localization), shim coil (used to correct B.sub.0 inhomogeneity) and RF coil (used for excitation and detection of signals). A controller, which may include a computer or other processing elements is used for control of the system. Shim and RF amplifiers may be used to drive their respective coils in the magnet. While the below described embodiments refer to MRI systems, it is to be understood that the invention is applicable to NMR or other imaging systems known in the art as well.

[0031] Performance Considerations for B.sub.0 Mapping

[0032] FIG. 2 shows a conventional gradient echo pulse sequence for B.sub.0 mapping. The sequence generally includes: a slice selective RF pulse (RF, SLICE) followed by a variable delay, A, an incremented phase encoding gradient (PE) and a readout gradient (RO). The acquisition is successively looped over N slices, M TE times, and finally P phase encode steps.

[0033] B.sub.0 mapping is commonly achieved using gradient echo images acquired with two or more different echo times. For samples or tissues with a single dominant resonance (e.g. water) the phase of the acquired images is given by the product of the difference in frequency of the dominant resonance from the carrier frequency (Δv) and the echo time (TE) and a spatially varying constant (ϕ.sub.0(r)) dependent on a number of parameters including RF coil system, reconstruction method and other instrumental factors.

ϕ=2π.Math.Δv.Math.TE+ϕ.sub.0(r)  Eq. 1

[0034] Thus the phase difference, Δϕ, between any two images is given by [2], such that the only unknown is the frequency offset Δv.

Δϕ=2π.Math.Δv.Math.ΔTE  Eq. 2.

[0035] Thus, assuming a single dominant resonance such that the distribution of frequencies in the image are solely related to the distribution of B.sub.0 field present, the spatial variation in Δϕ provides a direct measure of the inhomogeneity (ΔB.sub.0) present in the sample or tissue.

[0036] For a given level of phase noise (∈.sub.phase) and constant signal level, the measured phase (Δϕ.sub.msd) and true phase (Δϕ.sub.true) and measured (Δv.sub.msd) and true frequencies (Δv.sub.true) are given by Eqs. 3 and Eq. 4.

[00001]$\begin{array}{cc}\Delta {\varphi }_{\mathrm{msd}}=\Delta {\varphi }_{\mathrm{true}}+{ϵ}_{\mathrm{phase}}& \mathrm{Eq}.3\end{array}$$\begin{array}{cc}\Delta v_{\mathrm{msd}}=\Delta v_{\mathrm{true}}+\frac{{ϵ}_{\mathrm{phase}}}{2\pi .Math.\Delta \mathrm{TE}}& \mathrm{Eq}.4\end{array}$

[0037] Thus the accuracy of the measurement generally increases with increasing ΔTE.

[0038] However, since the phase Δϕ.sub.msd is periodic, with periodicity 2π, such that for long ΔTEs (i.e. more accurate measurements) if Δv.sub.true>1/(2.Math.ΔTE) the phase will alias such that the true phase and frequency is given by Eq. 5 and Eq. 6 with n being an integer describing the amount of wrapping.

[00002]$\begin{array}{cc}\Delta {\varphi }_{\mathrm{true}}=\Delta {\varphi }_{\mathrm{msd}}+2n\pi +{ϵ}_{\mathrm{phase}}& \mathrm{Eq}.5\end{array}$$\begin{array}{cc}\Delta v_{\mathrm{true}}=\Delta v_{\mathrm{msd}}+\frac{n}{\Delta \mathrm{TE}}+\frac{{ϵ}_{\mathrm{phase}}}{2\pi .Math.\Delta \mathrm{TE}}& \mathrm{Eq}.6\end{array}$

[0039] Thus Δv.sub.true may not be uniquely determined in the presence of wrapping and noise/uncertainty. Thus the most accurate (longest ΔTE) is prone to large systematic errors due to aliasing. For example, FIG. 3A shows the measured values (Δv.sub.msd) along the vertical axis as a function of the true value (Δv.sub.true) along the horizontal axis for a ΔTE=2 ms image over the range of +500 Hz. The wrapping points are at ±250 Hz. Displayed in FIG. 3B (lines added for visualization) is an expansion of the region from −270 to −240 Hz (box in FIG. 3a) displaying Δv.sub.msd for the ΔTE=2 ms image and a ΔTE=2 ms with an “error” of +10 Hz. For values between −270 and −260 Hz, both images alias and their difference is −10 Hz. However from −260 to −250 Hz, the ΔTE=2 ms image with a+10 Hz error term does not alias, while the ΔTE=2 ms images does resulting a difference of +490 Hz.

[0040] Methods of Unwrapping

[0041] To overcome this limitation two broad methods, spatial and temporal unwrapping have been described to unalias the data. See, for example: Jenkinson, M., Fast, automated, N-dimensional phase-unwrapping algorithm. Magn Reson Med, 2003. 49(1): p. 193-7; Hetherington, H. P., et al., Robust fully automated shimming of the human brain for high-field 1H spectroscopic imaging. Magnetic resonance in medicine, 2006. 56(1): p. 26-33; Robinson, S., H. Schodl, and S. Trattnig, A method for unwrapping highly wrapped multi-echo phase images at very high field: UMPIRE. Magn Reson Med, 2014. 72(1): p. 80-92; Geiger, Y. and A. Tal, Optimal echo times for multi-gradient echo-based B0 field-mapping. NMR Biomed, 2020. 33(7): p. e4316; and Dagher, J., T. Reese, and A. Bilgin, High-resolution, large dynamic range field map estimation. Magn Reson Med, 2014. 71(1): p. 105-17, each of which is incorporated herein by reference in its entirety.

[0042] Spatial unwrapping relies on the spatial continuity of the data to identify discontinuities in the B.sub.0 data as indicators of aliasing. See Jenksinson et al. These methods perform acceptably under conditions of moderate wrapping (i.e. moderate values of ΔTE), such that the accuracy of the measurements is limited by the LITE used. For temporal unwrapping, multiple ΔTE values are used, with the shortest ΔTE chosen to eliminate aliasing (see Hetherington et al.), i.e. ΔTE.sub.min<1/BW.sub.max=1/(2.Math.|Δv.sub.max|). Under this condition, the values in the shortest ΔTE are used to “correct” aliased values in longer ΔTE images. Unfortunately, if the maximum uncertainty in the shortest ΔTE image ∈.sup.max,

[00003]${ϵ}^{\max }=k.Math.\frac{{ϵ}_{\mathrm{phase}}}{2\pi .Math.\Delta \mathrm{TE}},$

where k is a scaling value such that all “outliers” are included, exceeds the bandwidth of the longer ΔTE image (1/ΔTE) the extent of aliasing will be ambiguous. Thus the “longer” ΔTE image cannot use an arbitrarily long ΔTE. For a normal distribution of noise/uncertainty, with

[00004]$\frac{{ϵ}_{\mathrm{phase}}}{2\pi .Math.\Delta \mathrm{TE}}$

being me standard deviation of the noise/uncertainty in the frequency domain, k would reflect number of standard deviations to account for the maximum noise value present in the sample

[0043] To overcome this limitation previously, a multi-ΔTE acquisition used a “boot strap” approach where: 1) the initial ΔTE value is chosen to be sufficiently short to eliminate aliasing; 2) each subsequent ΔTE image uses a factor of 2 in increasing duration and 3) each subsequent ΔTE image is used as the reference image to unwrap the next longest ΔTE image. Thus as the accuracy increases and bandwidth decreases in each subsequent image, the uncertainty in the “reference” ΔTE image decreases preserving the ability to unwrap the data unambiguously. However this method is limited by 1) the requirement to capture all inhomogeneity present (i.e. preclude aliasing) in the shortest ΔTE image and 2) the number of images required to get to arbitrarily long ΔTE values for a desired level of accuracy. These limitations manifest as longer acquisition times (more ΔTE images) for both increasing initial inhomogeneity and final accuracy.

[0044] Logistical Temporal Unwrapping

[0045] In illustrative embodiments of the invention, logistical temporal unwrapping provides the following unique advantages: 1) aliasing in the shortest ΔTE image is allowed, minimizing the need to accurately predict the maximum inhomogeneity present; and 2) the ΔTE values used can lengthen by factors larger than 2, enabling greater accuracy with fewer ΔTE images and thus shorter acquisition times. This methodology achieves these advantages by 1) using the longer ΔTE images to unwrap aliasing in the shortest ΔTE image and 2) utilizing the patterns of aliasing in all images to unwrap the most accurate image, i.e. the longest ΔTE image.

[0046] Unwrapping Aliasing in the Shortest TE Image

[0047] Aliasing results in discontinuities in the measured phase at well-defined frequency values in the bandwidth of the measurement. By comparing the measured values in the shortest ΔTE image with that measured in the longer ΔTE images, aliasing can be detected. Specifically a logistic temporal unwrapping parameter (LTUP) can be calculated, which enables the detection of wrapping in the shortest ΔTE image.

LTUP=(Δv.sub.short−Δv.sub.long)/(BW.sub.long/2)

[0048] FIG. 4 displays an example of the process graphically using ΔTE values of 2 and 9 ms (based on three MRI images with TE times of 4, 6 and 13 ms) with up to a factor of 2 in aliasing (±500 Hz, maximal expected bandwidth) present in the 2 ms image, in accordance with an embodiment of the invention. When the value of the LTUP is (LTUP=0,±2, or ±4, i.e., even numbers) the measured value in the shortest ΔTE image is unaliased, whereas if the value is (LTUP=±1,±3,±5, i.e., odd numbers) the value is aliased and can be corrected. For example, in the past conventional methodology, to achieve a true bandwidth of ±500 Hz, and achieve reasonable accuracy 4 ΔTEs would be acquired (i.e. 1, 2, 4 and 8 ms). Whereas, under the current method, only 2 LTEs would be acquired yielding equivalent maximal bandwidth, but with superior accuracy (9 vs 8 ms) with fewer ΔTEs (2 versus 4).

[0049] The effects of noise/uncertainty in the measurements is most pronounced at the aliasing points, where a few Hz of “noise” can dramatically change the measured values. For example, for a true value of 250 Hz, if ∈.sup.max=10 Hz in the ΔTE=2 ms image, the “true” value will map into a set of values between +240 Hz and +260 Hz. However, due to aliasing this set of true values {240,244, . . . 260} maps into the values {240,241, . . . 250,−249,−248, . . . −240}. This is graphically represented in FIG. 5, which shows the effects of uncertainty in the LTE=2 ms image on LTUP, in accordance with an embodiment of the invention, by two lines 2 ms(−10 Hz) and 2 ms(+10 Hz) representing the extremes of the numerical range. Assuming ∈.sup.max˜0 in the ΔTE=9 ms image, the true value of 250 Hz maps into a value of 27.8 Hz in the ΔTE=9 ms image. However, due to the uncertainty in the ΔTE=2 ms image, the LTUP can take on not only non-integral values, e.g. (240−27.7)/55.6=+3.82, but also very different values when the frequency is aliased, e.g. (−240/−27.7)/55.6˜−4.81.

[0050] The non-integral values of LTUP can be addressed by calculating LTUP*, where LTUP*=round(LTUP). After correction for non-integral values, however LTUP*=−5 or +4. Similarly for frequencies within ±∈.sup.max of −250 Hz, the LTUP* can take on values of −4 or +5. However, the set of LTUP* values over these two regions are distinct, such that the LTUP* can still be used to determine if the value has aliased or not. For example if the true frequency is 245 Hz, the ±10 Hz limits will place the measured value between 235 and 255 Hz, yielding value of 235 to 250 and −250 to −245 Hz in ΔTE=2 ms image with two possible LTUP* values, +4 and −5. This range of values in the ΔTE=2 ms image could also have been generated by a true frequency of −255 Hz, with LTUP* values of −4 and +5. Outside of the potential aliasing regions, the calculated LTUP* values for the extremes of uncertainty in the ΔTE=2 ms image are coincident (FIG. 5) and no ambiguity for aliasing exists. Thus the LTUP*, even in the presence of noise provides unambiguous unwrapping over these regions.

[0051] Unwrapping the ΔTE.sub.long Image

[0052] Once the ΔTE=2 ms image has been unwrapped and Δv.sub.short.sup.unwrap is outside of ±∈.sub.2 (the uncertainty/noise), of the alias point of the image±1/(2ΔTE.sub.long), Δv.sub.long.sup.msd can also be unwrapped. FIG. 6 assists in showing the unwrapping the ΔTE.sub.long Image, in accordance with an embodiment of the invention (FIG. 6 includes experimental uncertainty in both the ΔTE 2 ms and 9 ms images). Specifically, for each Δv.sub.long.sup.msd in the ΔTE=9 ms image there are nine theoretical potential values in the ΔTE=9 ms image, Δv.sub.short.sup.potential(n), which are given by:

Δv.sub.short.sup.potential(n)=Δv.sub.long.sup.msd+n(1/ΔTE.sub.long) where n=0,±1,±2,±3,±4  Eq. 8

forming the pairs

{(Δv.sub.short.sup.potential(−4),Δv.sub.long.sup.msd) . . . (Δv.sub.short.sup.potential(0),Δv.sub.long.sup.msd) . . . (Δv.sub.short.sup.potential(+4),Δv.sub.long.sup.msd)}

Thus the unwrapped value ΔTE=9 ms image, Δv.sub.long.sup.unwrap, is achieved by selecting the nth pair of values, minimizing ε.sub.short(n) where

ε.sub.short(n)=|Δv.sub.short.sup.potential(n)−Δv.sub.short.sup.unwrap|.  Eq. 9

[0053] However, if the value of Δv.sub.long.sup.msd is within ±∈.sub.long (the uncertainty/noise) of the alias point of the image±1/(2ΔTE.sub.long), “false” aliasing due to noise may have occurred. If “false” aliasing due to noise does occur, ε.sub.short(n) will be larger than ∈.sub.short, the uncertainty in the ΔTE.sub.short image. In this case the value Δv.sub.short.sup.unwrap is used for Δv.sub.long.sup.unwrap at this location to unwrap the ΔTE.sub.long image.

[0054] Alternatively, the process of unwrapping can be more generally viewed from the perspective of pattern matching, i.e. the finding the best numerical agreement between pairs of measured values (i.e. {(Δv.sub.short.sup.msd, Δv.sub.long.sup.msd)} and then sets of theoretical potential pairs i.e. {(Δv.sub.short.sup.potential(n), Δv.sub.long.sup.potential(n)), . . . } of values, where the number of pairs is dependent upon the extent of aliasing and maximum bandwidth.

[0055] Thus, in general, in the presence of noise or uncertainty in the ΔTE.sub.long image, unambiguous unwrapping can be achieved if uncertainties in the frequency region about the alias points in the image are spanned by continuous sets of differing values from the ΔTE.sub.short image. Continuity of these values in the ΔTE.sub.short image is achieved when the specific ΔTE values for both the ΔTE.sub.short and ΔTE.sub.long images are chosen such that the alias points±∈.sup.max do not overlap. Once the ΔTE.sub.short image is unwrapped, the ΔTE.sub.long image can be unwrapped.

[0056] As described, we have allowed for noise and uncertainty in the ΔTE.sub.short image but not the ΔTE.sub.long image(s). The presence of noise or uncertainty in the ΔTE.sub.long image pixels with true values near the alias points of the ΔTE.sub.long image can result in multiple LTUP* values. To unalias the ΔTE.sub.short image the regions of uncertainty about the alias points in the short ΔTE image are extended by the noise in the ΔTE.sub.long image, i.e. ∈.sub.max=∈.sub.short.sup.max+∈.sub.long.sup.max. As long as these bands, alias points±∈.sup.max, are within the continuous region of the ΔTE.sub.long image, the ΔTE.sub.short image can be unaliased as described previously (see FIG. 5) using the LTUP*. Since the ΔTE.sub.short image is continuous over these regions, the ΔTE.sub.long image can be unwrapped as described previously by creating the set of paired values {(Δv.sub.short.sup.potential(n),Δv.sub.long.sup.msd)} and selecting the nth pair such that ε.sub.short(n) is a minimum or Δv.sub.short.sup.unwrap when near the wrap points of the ΔTE.sub.long image and ε.sub.short(n) is large.

[0057] Creating Very High Accuracy Maps with Additional ΔTE Images.

[0058] As described, we have utilized 2 ΔTE images (ΔTE, =ΔTE.sub.short, ΔTE.sub.2=ΔTE.sub.long) to generate a B.sub.0 map. However, the ΔTE times chosen are limited by the maximum bandwidth of the sample, the aliasing patterns and the uncertainty in the individual ΔTE maps. Thus in general arbitrarily long ΔTE times for very high accuracy may not be possible under conditions of significant uncertainty and large sample bandwidth. To overcome this limitation additional ΔTE times (ΔTE.sub.3, ΔTE.sub.4, . . . ) can be acquired with increasing length. Each additional ΔTE image then forms a new “pair” of images, with the long ΔTE image from the previous pair (e.g. ΔTE.sub.2) now serving as the short ΔTE image and the new image (e. g. ΔTE.sub.3) serving as ΔTE.sub.long. Since the new ΔTE.sub.short image has already been unwrapped the requirement for non-overlapping alias point regions, can be relaxed, and the new ΔTE chosen must only satisfy the criteria that the maximal uncertainty in the new pair of images (e.g for (ΔTE.sub.2, ΔTE.sub.3), ∈.sup.max=∈.sub.2.sup.max+∈.sub.3.sup.max, be less than the bandwidth of the new ΔTE image, (e.g. 1/ΔTE.sub.3). Thus the process becomes recursive for each ΔTE image. Notably, as described previously, the multiplicative factor used in increasing the length of subsequent ΔTE times is not constrained by a factor of 2, but rather the actual noise/uncertainty in the images. This allows for very efficient increases in accuracy per additional ΔTE image acquired. At some point, the noise/uncertainty in the longer ΔTE image may not decrease inversely with increasing duration.

[0059] FIG. 7 compares scout images and B.sub.0 maps acquired from a healthy subject using conventional methodology using ΔTE times of 1, 2, 4 and 8 ms; and using the Logistical Temporal Unwrapping methodology in accordance with embodiments of the invention, with ΔTE times of 2 and 9 ms. The maps display every 4th slice in the image set, from a set of 33 slices. The images are plotted on a scale of ±100 Hz. The target region of interest for shimming is shown in red on the scout images (top row) and black on the B.sub.0 maps. Despite fewer ΔTE times, the Logistical Temporal Unwrapping methodology shows a small but significant reduction (10%, 20.2 versus 22.5) in standard deviation of B.sub.0 over the entire brain in comparison to the conventional method. This is consistent with the lengthening of the maximum echo time from 8 ms to 9 ms (12.5%).

[0060] FIG. 8 is an illustrative flow diagram of a Logistical Temporal Unwrapping methodology, in accordance with an embodiment of the invention. At step 801, the varying echo times (TE) are selected based on BW.sub.max and accuracy, which may be determined, for example, by ΔTE long. The B.sub.0 mapping image data is then acquired, step 803. Illustratively, at least two ΔTE images (ΔTE.sub.i=1 . . . N) may be generated based on at least three MRI images, wherein the subscripts i=1 . . . N refer to images with sequentially increasing ΔTE times. Advantageously, aliasing in the ΔTE.sub.1 image is permitted, and the ΔTE times of ΔTE.sub.1 and ΔTE.sub.2 are set such that the alias points at which wrapping occurs in ΔTE.sub.1 does not overlap with the alias points of ΔTE.sub.2. The phase differences in the ΔTE images are calculated in the ΔTE maps are calculated and converted to frequency, step 805. Each ΔTE image is then unwrapped. More specifically, as an example, LTUP* is calculated and the ΔTE.sub.1 (short image) is unwrapped, step 807. Then, the ΔTE.sub.2 (long image) is unwrapped using the unwrapped ΔTE.sub.1 (short image), step 809. If there are more ΔTE images, step 811, set the unwrapped ΔTE (short image) to the unwrapped ΔTE (long image), and the ΔTE (long image) to the next longest ΔTE image, step 813. If there are no more ΔTE images at step 811, set the final B.sub.0 map to the last unwrapped ΔTE.sub.N (long image), step 815. The shim currents and/or passive shims may then be calculated to generate the B.sub.0 correction field, step 817, which may then be applied to improve the B.sub.0 homogeneity.

[0061] Illustratively, the following is a more specific algorithm for unwrapping and setting of the final B0 map, in accordance with an embodiment of the invention. [0062] 1. Select a set of evolution times, [0063] a. A minimum of 3 evolution times to generate at least 2 ΔTE images [0064] b. Select the individual ΔTE times in the first two ΔTE image such that the frequency regions about the alias points (2/ΔTE), including the maximal uncertainty (∈.sub.short.sup.max+∈.sub.long.sup.max), do not overlap over the range of maximum expected frequencies (BW.sub.max) in the sample or tissue. [0065] c. Select additional ΔTE times such that (∈.sub.short.sup.max+∈.sub.long.sup.max)<BW.sub.long, where the subscripts short and long refer to images with sequentially increasing ΔTE times. [0066] d. Aliasing in the ΔTE.sub.short image is permitted over BW.sub.max. [0067] 2. Acquire the data. [0068] 3. Calculate the phase in the ΔTE images and convert to frequency.

[00005]$\Delta v_{\mathrm{msd}}=\frac{\Delta {\varphi }_{\mathrm{msd}}}{2\pi \Delta \mathrm{TE}}$ [0069] 4. Unwrap the ΔTE.sub.short image: [0070] a. Divide the bandwidth of the ΔTE.sub.short image into two conceptual regions: [0071] i. Frequencies within ∈.sub.short.sup.max+∈.sub.long.sup.max of the alias points of the ΔTE.sub.short.sup.image [0072] ii. Frequencies outside of this region [0073] b. For frequencies within region (ii) [0074] i. Calculate the sets of theoretical values for LTUP* over BW.sub.max as a function of extent of aliasing (n=0,±1, . . . ) [0075] ii. Calculate LTUP* from the measured data LTUP*=ROUND((Δv.sub.short.sup.msd−Δv.sub.long.sup.msd)/(BW.sub.long/2)) [0076] iii. Use the measured LTUP* to look up the value of n from the theoretical values in 4.b.i [0077] iv. Calculate Δv.sub.short.sup.unwrap=Δv.sub.short.sup.msd+n.Math.BW.sub.short, [0078] c. For frequencies within region (i) [0079] i. Perform the same calculation as in 4.b.i over the range of frequencies spanned by the alias point±(∈.sub.short.sup.max+∈.sub.long.sup.max) [0080] ii. Calculate LTUP* from the measured data LTUP*=ROUND((Δv.sub.short.sup.msd−Δv.sub.long.sup.msd)/(BW.sub.long/2)) [0081] iii. Use the measured LTUP* to look up the value of n from the set of theoretical values in 4.c.i. Each aliased region will have multiple but distinct LTUP* values for a given value of n. [0082] iv. Calculate Δv.sub.short.sup.unwrap=Δv.sub.short.sup.msd+n.Math.BW.sub.short, [0083] 5. Unwrap the long ΔTE image: [0084] a. Using the measured value in the ΔTE.sub.long image calculate the potential pairs of frequency values for the (Δv.sub.short.sup.potential(n),Δv.sub.long.sup.msd) over BW.sub.max {(Δv.sub.short.sup.potential(n), Δv.sub.long.sup.msd) . . . } where Δv.sub.short.sup.potential(n)=Δv.sub.long.sup.msd+n(1/ΔTE.sub.long) [0085] b. Calculate ε.sub.short(n)=|Δv.sub.short.sup.potential(n)−Δv.sub.short.sup.unwrap| for each pair [0086] c. Select the value of n that minimizes ε.sub.short(n). [0087] d. If ε.sub.short(n)<∈.sub.short.sup.max Assign Δv.sub.long.sup.unwrap=Δv.sub.long.sup.msd+n(1/ΔTE.sub.long) [0088] If ε.sub.short(n)>∈.sub.short.sup.max Assign Δv.sub.long.sup.unwrap=Δv.sub.short.sup.unwrapped [0089] 6. Unwrap additional longer ΔTE images (ΔTE.sub.longer) if present by . . . . [0090] a. Substituting the parameters for the ΔTE.sub.short image with those from the ΔTE.sub.long image e.g. Δv.sub.short.sup.unwrap=Δv.sub.long.sup.unwrap [0091] b. Substituting the parameters for the ΔTE.sub.long image with those from the ΔTE.sub.longer image e.g. Δv.sub.long.sup.msd=msdΔv.sub.longer.sup.msd [0092] c. Repeating the calculation described in 5 [0093] 7. After all images have been unwrapped, use the values from the longest ΔTE image to generate the final map. [0094] 8. Using the final B.sub.0 map described in 7 calculate a set of shim currents and/or determine passive shimming for generating a correction B.sub.0 field of approximately equal amplitude but opposite sign. [0095] 9. Apply the calculated shim currents and/or passive shims to reduce the B.sub.0 inhomogeneity in the sample.

[0096] FIG. 9 shows three wrapped ΔTE maps and the final unwrapped B.sub.0 map obtained from a spherical phantom, in accordance with various embodiments of the invention. More particularly, the figure shows nine slices from a set of 33 slices acquired from the spherical phantom in which the shims were purposely offset to generate a significant amount of inhomogeneity. The bottom three rows show the wrapped maps from the: ΔTE=2 ms image (bottom row); 9 ms image; (2.sup.nd row from bottom); and 21 ms image, plotted on their scales of +/−BW/2 where BW=1/ΔTE. The top row shows the final B.sub.0 map obtained from the unwrapped last ΔTE=21 ms map, using the above-described methodology and embodiments of the invention.

[0097] Embodiments of the invention may be located implemented in part in any conventional computer programming language. For example, preferred embodiments may be implemented in a procedural programming language (e.g., “C”) or an object oriented programming language (e.g., “C++”, Python). Alternative embodiments of the invention may be implemented as pre-programmed hardware elements, other related components, or as a combination of hardware and software components.

[0098] Embodiments also can be implemented in part as a computer program product for use with a computer system—for example, the controller of the MRI system described above. Such implementation may include a series of computer instructions fixed either on a tangible medium, such as a computer readable medium (e.g., a diskette, CD-ROM, ROM, or fixed disk) or transmittable to a computer system, via a modem or other interface device, such as a communications adapter connected to a network over a medium. The medium may be either a tangible medium (e.g., optical or analog communications lines) or a medium implemented with wireless techniques (e.g., microwave, infrared or other transmission techniques). The series of computer instructions embodies all or part of the functionality previously described herein with respect to the system. Those skilled in the art should appreciate that such computer instructions can be written in a number of programming languages for use with many computer architectures or operating systems. Furthermore, such instructions may be stored in any memory device, such as semiconductor, magnetic, optical or other memory devices, and may be transmitted using any communications technology, such as optical, infrared, microwave, or other transmission technologies. It is expected that such a computer program product may be distributed as a removable medium with accompanying printed or electronic documentation (e.g., shrink wrapped software), preloaded with a computer system (e.g., on system ROM or fixed disk), or distributed from a server or electronic bulletin board over the network (e.g., the Internet or World Wide Web). Of course, some embodiments of the invention may be implemented as a combination of both software (e.g., a computer program product) and hardware. Still other embodiments of the invention are implemented as entirely hardware, or entirely software (e.g., a computer program product).

[0099] Although various exemplary embodiments of the invention have been disclosed, it should be apparent to those skilled in the art that various changes and modifications can be made which will achieve some of the advantages of the invention without departing from the true scope of the invention.