Research Articles

Download PDF (3 MB)

TOMOGRAPHY, June 2017, Volume 3, Issue 2: 79-88
DOI: 10.18383/j.tom.2017.00003

A Robust Method for Estimating B0 Inhomogeneity Field in the Liver by Mitigating Fat Signals and Phase-Wrapping

Antonis Matakos1, James M. Balter1, Yue Cao1

1Departments of Radiation Oncology;2Biomedical Engineering; and3Radiology, University of Michigan, Ann Arbor, Michigan


We developed an optimized and robust method to estimate liver B0 field inhomogeneity for monitoring and correcting susceptibility-induced geometric distortion in magnetic resonance images for precision therapy. A triple-gradient-echo acquisition was optimized for the whole liver B0 field estimation within a single-exhale breath-hold scan on a 3 T scanner. To eliminate chemical-shift artifacts, fat signals were chosen in-phase between 2 echoes with an echo time difference (∆TE) of 2.3 milliseconds. To avoid phase-wrapping, other 2 echoes provided a large field dynamic range (1/∆TE) to cover the B0 field inhomogeneity. In addition, using high parallel imaging factor of 4 and a readout-bandwidth of 1955 Hz/pixel, an ∼18-second acquisition time for breath-held scans was achieved. A 2-step, 1-dimensional regularized method for the ∆B0 field map estimation was developed, tested and validated in phantom and patient studies. Our method was validated on a water phantom with fat components and air pockets; it yielded ∆B0-field maps that had no chemicalshift and phase-wrapping artifacts, and it had a


Magnetic resonance imaging (MRI) is commonly used for planning radiation treatment (RT) for patients with liver cancers because of its superior soft tissue contrast compared with computed tomography (CT). Recent investigations (13) show the potential to generate radiation attenuation maps based on magnetic resonance (MR) images, thus eliminating the need for CT scans for RT planning, as well as the associated costs and burdens for the patients. However, ensuring geometric accuracy of MR images is crucial for its use as the sole modality for RT planning. Geometric distortion of clinical MR images needs to be assessed and corrected if necessary. System-level distortions are sufficiently corrected in state-of-the-art scanners, but the subject-induced distortions due to susceptibility and chemical-shift effects still need to be accounted for and corrected (49).

A common approach for quantifying and correcting subject-induced geometric distortions is to use the acquired B0-inhomogeneity field map from the phase difference of 2 gradient echo (GRE) images (5). Estimation of the ΔB0-field map from GRE images in the liver faces several challenges, some of which are common to all body sites such as poor signal-to-noise ratio (SNR) and phase-wrapping at locations with large susceptibility effects, and others that are more relevant for the liver such as respiratory motion, chemical-shift, acquisition time, and even geometric accuracy of the image data itself. The phase-wrapping artifact in the estimated field map is a serious challenge. Methods have been proposed to unwrap the field map (1013). Phase-unwrapping, when applied to the liver, significantly increases the computational load for image processing, and it may not guarantee that it will work in all cases. An alternative approach is to increase the available dynamic range of the field by shortening the time difference between the 2 echoes (ΔTE), with a trade-off of reduced SNR. Clinical MRI scanner hardware limits the minimum separation of echo times (TEs) within a single repetition time (TR), usually to ∼1 milliseconds or slightly less, corresponding to a field dynamic range of ±500 Hz. Furthermore, if the TEs are not tuned for fat and water to be in-phase, the field map can be corrupted by chemical-shift artifacts, particularly for a body site like the liver with substantial fat signals. However, choosing both TEs to have fat and water in-phase reduces the field dynamic range and thus increases the risk for phase-wrapping artifacts. More advanced methods for field map estimation make use of ≥3 TEs to avoid phase-wrapping (14), or estimate the field map and fat and water components (1518), thereby reducing or eliminating the chemical-shift artifact. However, modeling fat and water components requires high accuracy in the estimation to avoid errors propagating into the field map (16, 18). A joint estimation of fat and water components along with a field map (15, 17) also substantially increases the computation complexity. The multiple echo methods (19, 14) that are based on multiple separated acquisitions do not conform well to the time constraints on echo spacing posed by a single acquisition. Finally, respiratory motion is an important source of phase errors for abdominal scanning, but it can be reduced through use of breath-held acquisition (with associated limitations on acquisition time). The geometric distortion induced by the susceptibility effect in the field map also needs to be considered, which can be mitigated to a submillimeter level by using a high readout bandwidth.

In this work, we propose an optimized and efficient protocol using 3 GREs to estimate the B0-inhomogeneity field map of the whole liver. This method was optimized to overcome the challenges discussed. First, a breath-held acquisition is used to reduce respiratory motion-induced phase errors in complex images, constraining the acquisition time to ∼18 s to be suitable for the majority of patients. Next, we optimize TE spacing in a single breath-held acquisition to achieve both elimination of the chemical-shift effect of fat and extension of the dynamic range of the ΔB0 field. In addition, we minimize geometric distortion in acquired images and in the subsequently estimated B0-inhomogeneity field maps by using a high readout bandwidth. With this optimized imaging protocol, we propose a robust, 2-step, 1-dimensional (1D) regularized method to estimate the abdominal ΔB0-field map free of motion artifacts, phase-wrapping and chemical-shift effects. We tested our method in phantom and human liver studies.

Materials and Methods


The chemical-shift of fat alters the phase of MR signals of voxels that contain both fat and water. The complex MR signal ρj at voxel j from a GRE is commonly expressed as follows:

where ρfj and ρwj are, respectively, fat and water magnitude components, ωj is the B0 inhomogeneity field, Δf is the frequency shift of fat signals due to the chemical-shift effect, and TE is the TE. To mitigate the chemical-shift effect of fat, the conventional 2-GRE method requires water and fat signals to be in phase, where ΔfTE is an integer, for both echoes, which results in limited TE choices. An alternative is to estimate the field map and fat and water components simultaneously using 3-point Dixon-based approaches (15, 17). The current state-of-the-art 3-point Dixon-based methods require spacing 3 echoes such that phase shifts between fat and water are evenly distributed within 2π (15, 17, 20, 21). This requirement poses more constraints on echo spacing and is not desirable for scanning the whole liver in a single breath-hold. If the chemical-shift effect of fat is ignored, the estimated field map can have discontinuities at interfaces between fat and water.

We propose a new acquisition method to eliminate the chemical-shift effect of fat. Let us consider a simple situation wherein a voxel contains either water or fat (ie, no mixing). The phase difference of 2 GREs with a TE difference of ΔTE at voxel j is expressed as follows:

It can be seen that the phases between fat and water signals are discontinuous, dependent upon the TE difference of ΔTE. If a voxel contains a mixture of fat and water, the extra phase due to the fat signals can vary, depending upon the portion of fat in the voxel and TE used. Using such a phase map for estimating the related geometric distortion will cause errors. To remove the phase error caused by the frequency shift of fat signals, we can choose fat signals in-phase between the 2 echoes such as ΔfΔTE = k, an integer, and the phase difference of the fat component in Eq. [2] becomes the following equation:
where the phase discontinuity between fat and water is eliminated. Note that Eq. [3] holds true for a voxel containing a mixture of fat and water signals, as expressed in Eq. [1]. When the chemical-shift-induced phase is removed by having the fat signal in-phase between the 2 echoes, there is no need to estimate water and fat components in Eq. [1], which dramatically reduce the computation complexity. Most importantly, maintaining the fat signal in-phase does not require fat and water signals to be in-phase in either echoes, which provide flexibility to select TEs, for example, choosing the first echo to have a minimum TE. In addition, a minimum TE for the first echo is the first step toward shortening the total acquisition time. Under this condition, the B0 inhomogeneity field at voxel j is simply as follows:
where πΔφjππ/ΔTEωjπ/ΔTE, regardless of whether the voxel contains any fat.

Next, let us look at the issue of the dynamic range of the field map provided by the echoes, which is the root cause for the phase-wrapping. The dynamic range of the field frequency is determined by ΔTE of the 2 echoes (1/ΔTE). A short ΔTE will result in a large dynamic frequency range. Requiring fat signals to be in-phase between the 2 echoes to remove the chemical-shift effect in the phase-difference map limits the TE difference of ΔTE to be 2.3 milliseconds at 3 T and 4.5 milliseconds at 1.5 T, resulting in a dynamic range of ∼435 Hz at 3 T and 222 Hz at 1.5 T. If the B0 field inhomogeneity is greater than ±217.5 Hz at 3 T, using a ΔTE of 2.3 milliseconds will cause phase-wrapping; otherwise, it will not. To avoid the phase-wrapping artifact in the estimated field map (or avoid the estimate trapped in a local minimum), it is desirable to decrease ΔTE. An optimal solution for overcoming both the chemical-shift effect and the phase-wrapping effect is to have 2 echoes with fat signal in-phase and 2 echoes with a minimum ΔTE. In addition, a short ΔTE will result in a reduced phase difference (low signal or low SNR) between the 2 echoes, the effect of which on field mapping needs to be tested.

Our protocol for field mapping in the liver optimizes several factors, including

  1. yielding a sufficiently short total acquisition time to support imaging during a single breath-hold (short first TE and ΔTE);

  2. providing sufficient dynamic range for mapping the large B0 inhomogeneity field at the liver dome without a phase wrap (short ΔTE);

  3. avoiding phase discontinuity from the chemical-shift effect (fat signals in-phase between the 2 echoes); and

  4. yielding adequate SNR.

We implemented a triple-echo acquisition paradigm. Specifically, we selected the first and third echoes to have ΔTE of 2.3 milliseconds (fat signals in-phase). We further inserted a second echo between the first and third echoes, by which the ΔTE between the first and second echoes is short enough to provide a large dynamic range. Furthermore, we selected a minimal TE for the first echo to ensure a short total acquisition time suitable for a single breath-hold. Finally, we used a high readout bandwidth to further shorten the total acquisition time and reduce geometric distortion in the field map caused by B0 inhomogeneity to a negligible level (<0.5 mm; see Results).

To estimate the B0-inhomogenity field map ω^ based upon this protocol, we developed a 2-step regularized estimation method to eliminate both the chemical-shift effect of fat and the phase-wrapping effect. We applied a 1D regularized estimation method (19) to the complex image volumes of the first and third echoes (ρ1 and ρ3) that were free from the chemical-shift effect of fat, by minimizing the following:

where R(ω) is a quadratic regularizer that enforces field map smoothness, and λ is a parameter tuned to suppress noise without excessively compromising the spatial resolution of the field map. The quadratic regularizer is a natural choice because of the relative smoothness of the field map and computation efficiency. The intent of the use of this regularizer is to stabilize the solution. The λ choice is based on a previous investigation (8), in which the full-width–half-maximum of the point-spread function was found to be 1.02 pixels. As discussed before, the problem in Eq. [5] is nonconvex. To avoid trapping the field map solution in the wrong local minima, the initialization of the B0 inhomogeneity frequency has to have no phase-wrapping, thus ensuring that the optimization in Eq. [5] is near the global minimum. This requires that the dynamic range of the 2 echoes used for initialization has to be large enough to cover the field range. Therefore, we initialized Eq. [5] using the field map ω^1 estimated from the first 2 echoes with short ΔTE by Eq. [4]. The required ΔTE for the liver was estimated in the experiment.


First, we investigated the extent of the susceptibility-induced B0 inhomogeneity field in the liver to determine the necessary dynamic range of the field for optimal echo spacing in the triple-GRE protocol. We next investigated the effectiveness of the proposed field map estimation strategy, including the triple-GRE protocol with optimal echo spacing and the 2-step, 1D regularized method for ΔB0-field estimation. We validated the proposed method first in a phantom study with known water and fat contents, as well as geometric specifications, and subsequently, evaluated in a liver patient study for applicability.

Assessment of the Extent of Susceptibility-Induced Field Inhomogeneity in the Liver

Liver MRI scans of 13 patients who were enrolled in an IRB-approved study were acquired on a 3 T scanner (Siemens Healthineers, Erlangen, Germany). After the localizer, patient-specific automated shimming of the liver was performed based on an ultrafast, noniterative, fully automated and stable shimming procedure, provided by the vendor. The shimming setting was then locked for the imaging session. The B0 inhomogeneity field maps were acquired from the conventional phase differences of dual GREs with TE1/TE2/TR = 4.92/7.38/106 milliseconds and a readout bandwidth of 290 Hz/pixel within a single breath-hold. Both echoes had the same gradient polarity. Because of the time constraint of the breath-held acquisition, 10 coronal sections in the center of the liver dome were acquired with an in-plane resolution of 3.5–4.4 mm and section thickness of 3.75 mm. The complex image volumes were saved for the calculation of the phase differences. The phase-difference maps from the 2 echoes were unwrapped using a quantitative MRI analysis software package (FSL, FMRIB, Oxford, UK) (12). For comparison, the same images were acquired with the subject breathing freely. The resulting calculated B0 inhomogeneity maps (ΔB0 in Hz) were used to estimate the extent and location of the B0 inhomogeneity field to optimize echo spacing in the triple-GRE acquisition.

Evaluation of the Triple-Gradient-Echo Method for Liver Field Mapping

Phantom Study.

A phantom consisting of a large cylindrical compartment filled with doped water and 5 cylindrical air pockets, with known geometrical specifications, was used. To assess the chemical-shift effects on the estimated field map by our proposed method, we filled 3 of the pockets with inserts containing 100% fat (oil), a fluid mixture with 20% fat and 80% water (intralipid with 20% lipid concentration), and a fluid mixture with 10% fat and 90% water (50% intralipid diluted with 50% water). The other 2 pockets were left empty, which created large inhomogeneities in the B0 field at the air–water interfaces because of susceptibility effects, and allowed us to assess the geometric distortion of the cylinder shapes in the field map and magnitude image.

The phantom was scanned on a 3 T scanner (Siemens Healthineers, Erlangen, Germany) using a 3-dimensional (3D) triple-GRE sequence from our optimized protocol for the liver with the following parameters: field-of-view (FOV) = 280 × 280 × 141 mm3, matrix = 128 × 128 × 64, flip angle of 9°, TR of 7 milliseconds, TE1/TE2/TE3 = 1.89/2.79/4.19 milliseconds, readout bandwidth = 1955 Hz/pixel, and parallel imaging (GRAPPA) with a 4× acceleration. The first 2 echoes had a short ΔTE (0.9 milliseconds), resulting in a dynamic range of 1111 Hz, and the first and third echoes had ΔTE of 2.3 milliseconds, which kept fat signals in-phase between the 2 echoes. All complex image volumes were corrected for 3D gradient nonlinearity.

Patient Study.

The triple-GRE acquisition was optimized to have a sufficient dynamic range to cover the B0 field inhomogeneity of the liver and to avoid phase-wrapping. To estimate field maps of the whole liver, 3D triple-GRE image volumes were acquired in axial orientation on a 3 T scanner (Siemens Healthineers, Erlangen, Germany) with the following parameters: FOV = 430 × 416 × 224 mm3, matrix = 128 × 124 × 64, flip angle = 9°, TR = 6.7 milliseconds, TE1/TE2/TE3 = 1.89/3.04/4.19 milliseconds, readout bandwidth = 1955 Hz/pixel, and parallel imaging (GRAPPA) with a 4× acceleration. To insert a second echo between the first and third echoes (2.3 milliseconds), the gradient polarity of the second echo was set opposite to the first and third ones. This acquisition took ∼18 s and was done during a single exhale breath-hold. Again, the shimming was done after the localizer and locked for the whole imaging session. We also acquired 3D anatomic T1-weighted image volumes using a clinical breath-holding protocol (using a volume interpolated breath-hold examination sequence [VIBE]) in coronal orientation with FOV = 432.25 × 460 × 260 mm3, matrix = 120 × 128 × 52, flip angle = 5°, TR = 10 milliseconds, TE = 1.57 milliseconds, and readout bandwidth = 440 Hz/pixel. All complex image volumes were corrected for 3D gradient nonlinearity. Scanning was performed in 5 volunteer patients using this technique on an IRB-approved protocol.

Field Map Estimation and Image Analysis.

The field maps of the phantom and the abdominal regions of the patients, including the liver, were estimated from the triple-GRE image volumes using our proposed method, and compared with standard field maps estimated using the first 2 echoes and the first and third echoes. Because the gradient polarity of the second echo was reversed compared to the first and third echoes in the patient study, the gradient polarity mismatch resulted in a linear ramping in the phase image of the second echo (2224). This phase ramping was approximated linearly and corrected before the phase-difference maps were calculated from the first 2 echoes for initialization of Eq. [5].

The phantom was also used to assess geometric distortion in the estimated field map. The geometric distortion caused by the B0-inhomogeneity field in a 3D acquisition occurs only along the readout direction and given by Chang and Fitzpatrick (4, 5):


where ω is the B0-inhomogeneity field at a voxel, BWf is the readout-bandwidth per pixel, and Δx is the pixel size in the readout direction. The geometric distortion of the phantom at the empty air pockets estimated by this equation was compared with the geometric specification from the phantom design. In the patients, after confirming the susceptibility-induced distortion in the estimated field map to be <0.5 mm, which was considered to be negligible, we used the field maps to assess the distortion of the clinical 3D T1-weighted images that had a readout bandwidth of 440 Hz/pixel (optimized for image quality).


Extent of Susceptibility-Induced B0 Inhomogeneity Field in the Liver

The greatest B0 inhomogeneity field was observed at the dome of the liver near the lung–liver interface (Figure 1). For 12 of the 13 patients, we observed no absolute ΔB0 >434.8 Hz, for the frequency range that ΔTE of 1.15 milliseconds can provide (Table 1). For the 1 outlying patient, 0.28 cc of the imaged liver volume had ΔB0 >434.8 Hz but <514 Hz. For 11 patients, the absolute ΔB0 was <217.4 Hz for ≥98% of the voxels within the imaged liver volumes. These data indicate that 1.15 milliseconds of ΔTE provides a sufficient dynamic range for ΔB0 in the liver.

Figure 1.

A typical example of ΔB0 maps (color-coded and overlapped on gradient echo [GRE] magnitude maps). The ΔB0 maps are displayed between 250 and −250 Hz. Note that very few voxels have ΔB0 >200 Hz. Yellow contours depict outlines of the liver.

Table 1.

The Extent of ΔB0 Field Inhomogeniety in the Liver

Percentage Volume
Volume (cc)
(0, 217.4) Hz (217.4, 434.8) Hz >434.8 Hz Max ΔB0 (Hz) (0, 217.4) Hz (217.4, 434.8) Hz >434.8 Hz
Pt 1 99.9% 0.1% 0.0% 309 558.1 0.4 0.0
Pt 2 99.7% 0.3% 0.0% 273 684.7 1.8 0.0
Pt 3 99.0% 1.0% 0.0% 374 444.9 4.6 0.0
Pt 4 98.0% 1.9% 0.1% 514 577.9 11.3 0.3
Pt 5 95.4% 4.6% 0.0% 429 431.7 20.6 0.0
Pt 6 100% 0.0% 0.0% 231 573.5 0.1 0.0
Pt 7 99.8% 0.2% 0.0% 315 445.6 1.0 0.0
Pt 8 99.6% 0.4% 0.0% 346 526.4 2.2 0.0
Pt 9 99.7% 0.3% 0.0% 350 716.9 2.4 0.0
Pt 10 99.9% 0.1% 0.0% 263 661.6 0.5 0.0
Pt 11 99.2% 0.8% 0.0% 380 278.1 2.2 0.0
Pt 12 92.3% 7.7% 0.0% 368 459.7 38.6 0.0
Pt 13 98.7% 1.3% 0.0% 312 404.7 5.4 0.0

i] Note: The absolute frequency is used to define the frequency ranges. In addition, 217.4 and 434.8 Hz correspond to ΔTEs of 2.3 and 1.15 milliseconds, respectively.

Evaluation of the Proposed Field Mapping Method for the Whole Liver

Phantom Study.

Magnitude and phase images of the phantom acquired by the triple-GRE sequence are shown in Figure 2. In the phase images of all 3 echoes, phase discontinuities are visible between the main water container and the oil insert. Large B0-inhomogeneity field is also visible near the boundary of the 2 empty (air) packets. The field maps obtained using Eq. [4] from the first 2 echoes and the first and third echoes are shown in Figure 3, A and B, respectively. Note that the field map from the first 2 echoes shows no phase-wrapping, except noise, near the interface of air and water, indicating the sufficient dynamic range (∼1100 Hz) provided by ΔTE of 0.9 milliseconds to fully cover the strong B0-inhomogeneity effect induced by air; but, phase discontinuities occur at all 3 inserts with fat, even at 10% concentration (Figure 3A). In contrast, the field map from the first and third echoes shows no discontinuities in all 3 inserts with fat because of fat in-phase in the first and third echoes, but phase-wrapping artifacts near the empty air pocket because of 50% reduction in the dynamic range provided by the first and third echoes (Figure 3B). However, the field maps estimated by our proposed 2-step method show no phase-wrapping near the air–water interface and no discontinuity caused by the chemical-shift of fat (Figure 3C).

Figure 2.

Magnitude (top row) and phase (bottom row) images of 3 GREs of a water phantom with echo times (TEs) of 1.89, 2.79, and 4.19 milliseconds (from left to right). Two cylindrical inserts contained air, while the others contained 100% oil (top right pocket), 20% fat and 80% water (bottom right), and 10% fat and 90% water (bottom left). Note the phase discontinuity at the 3 insets, and large susceptibility-induced phase variation, as well as phase-wrapping near the 2 air pockets (center and top left pockets).

Figure 3.

Field maps of the phantom. The standard phase-difference map of the first 2 echoes with ΔTE = 0.9 milliseconds shows no phase-wrapping, but phase discontinuity at 3 inserts with fat (the top right insert and 2 bottom inserts) (A). The standard phase-difference map of the first and third echoes with ΔTE = 2.3 milliseconds (fat signals in-phase) shows no chemical-shift artifact at the inserts with fat, but phase-wrapping at the 2 empty pockets (center and top left) (B). The field map from our method shows no chemical-shift and phase-wrapping artifacts (C). All 3 field maps are displayed from −300 to 300 Hz. The field map in (C) converted to represent the distortion in millimeters for a readout bandwidth of 1955 Hz/pixel is color-coded and overlaid on the magnitude image from the first echo (D). The maximum resulting distortion is <0.5 mm.


The large field inhomogeneity near the air–water interface caused marginal geometric distortion in both field and magnitude maps, the estimated maximum value of <0.5 mm (Figure 3D), because of using 1995 Hz/pixel readout bandwidth in the acquisition. This level of distortion is an acceptable level for precision RT (25), and therefore, there is no need for additional geometric correction of the field map using a joint reconstruction method or other methods.

Patient Study.

Typical magnitude and phase images of 3 GREs from the abdominal region of a patient are shown in Figure 4. Phase artifacts are seen because of both the large susceptibility effect near the liver dome and the fat chemical-shift effect.

Figure 4.

Magnitude (first and third rows) and phase (second and fourth rows) images of 3 GREs of 2 sections in the abdominal region of a patient. TEs were 1.89, 3.04, and 4.19 milliseconds from left to right columns. The phase maps of the first section (second row) shows phase discontinuities (red arrow) because of the fat chemical-shift effect. The phase maps of the second section near the liver dome (fourth row) show phase-wrapping artifacts (arrows) due to the large susceptibility effect.


Similarly to the phantom study, the field maps estimated by our proposed 2-step method were compared with those produced by the standard phase difference using the first 2 echoes and first and third echoes. Figure 5 shows 2 sections from 1 patient, 1 in the lower abdomen region where substantial fat components are present and 1 on the top of the liver near the diaphragm at the lung–liver interface, where there is a strong B0-inhomogeneity field. The field maps from the first 2 echoes (ΔTE = 1.15 milliseconds) show no phase-wrapping, indicating the dynamic range of ∼870 Hz is sufficient for field inhomogeneity near the lung–liver interface, but also show clearly visible chemical-shift artifacts (discontinuities) at all the fat–water interfaces. The field maps from the first and third echoes show no chemical-shift artifacts at the fat–water interfaces because of fat signals being in-phase between the 2 echoes (ΔTE = 2.3 milliseconds), but phase-wrapping artifacts near the lung–liver interface due to the 50% reduced dynamic field range to 453 Hz by the first and third echoes. The field map estimated by our proposed 2-step method shows no phase-wrapping artifact near the lung–liver interface and no visible discontinuity caused by the fat chemical-shift (Figure 5, C and F).

Figure 5.

Patient's estimated field maps of 2 sections. Top row: an abdomen region with substantial fat components. Bottom row: the liver dome with large field inhomogeneity. Left column: the field maps from the standard phase-difference estimation of the first 2 echoes with ΔTE = 1.15 milliseconds show no phase-wrapping in (D) but fat chemical-shift artifacts in (A) (arrow). Middle column: the field maps from the standard phase-difference estimation of the first and third echoes with ΔTE = 2.3 milliseconds (fat signals in-phase) show no chemical-shift artifact in (B) but phase-wrapping artifacts at the lung–liver interface in (E). Right column: The field maps estimated by our 2-step method show neither chemical-shift nor phase-wrapping artifacts (C and E). Field maps are displayed from −350 to 350 Hz.


The field maps for all 5 patients are shown in Figure 6. The maximum field inhomogeneity in the liver was located at the dome near the interface with lung. The maximum inhomogeneity fields of the 5 patients were ∼309, 274, 222, 226, and 220 Hz. The high readout bandwidth of 1995 Hz/pixel used in the triple-GRE acquisition mitigated the geometric distortion to 0.53 ± 0.025 (STD) mm as an average for the field maps themselves. This level of distortion is considered negligible for RT applications.

Figure 6.

Color-coded ΔB0-inhomogeneity field maps of the 5 patients estimated by our proposed method. The maximum intrahepatic ΔB0-inhomogeneity field was observed at the top of the liver near the interface with the lung. The field maps are displayed between −250 and 250 Hz.


The estimated field maps were used to evaluate susceptibility-induced geometric distortion in clinical 3D T1-weighted images acquired by a VIBE sequence with 440 Hz/pixel readout bandwidth. The averaged maximum distortion at the lung–liver interface of 5 patients was 2.05 ± 0.29 (STD) mm. Most interestingly, ΔB0 decayed slowly from the interface of lung and liver (Figure 7). At a distance of 40 mm away from the dome of the liver, the average ΔB0 was 114.5 ± 16.0 Hz (STD), which could still cause 0.93 mm of distortion in the readout direction for the clinical images acquired with the VIBE sequence and parameters.

Figure 7.

ΔB0 decay curves from the interface of the liver and lung. The curves were obtained along a 75-mm-long path perpendicular to and starting from the dome of the liver. The inset shows an example of the path.



This study presented an optimized but simple and efficient protocol involving a single breath-hold scan using triple-GREs and a 2-step 1D regularized method for estimating the B0-inhomogeneity field map in the liver. Our triple-GRE acquisition has no requirement for the phase differences between fat and water for all echoes, thus reducing echo spacing constraints. We only require that the fat signals be in-phase for the first and third echoes. Our optimized acquisition and processing steps are able to mitigate effects of respiratory motion, chemical-shift of fat, phase-wrapping, and geometric distortion caused by susceptibility. Particularly, we designed 2 of the 3 echoes to have fat signals in-phase, which eliminate the chemical-shift of fat in the phase difference and thus the need for the estimation of fat and water components. As a result, the field map estimation becomes a 1D problem. We implemented a 1D regularized method to reduce noise influence. To avoid phase-wrapping, we initialized the optimization process using the phase-difference map obtained from a short TE difference of the first 2 echoes. We used a high readout bandwidth to reduce the susceptibility-induced geometric distortion to <0.5 mm in the location with the maximum susceptibility change. Overall, our approach is efficient and practical for integration into a clinical protocol for assessing the degree of susceptibility-induced distortion in clinical images and for correcting them if deemed necessary.

Breath-holding is a common approach in abdominal MRI to mitigate respiratory motion effects on images. Imaging during free breathing can cause several issues for abdominal field mapping. First, respiratory motion can lead to phase errors, making the estimated field map inaccurate (Figure 8). Second, a moving organ can cause additional magnetic-susceptibility effects. Third, free-breathing and breath-held acquisitions can create mismatched field distributions due to differences in organ locations. Using breath-holding for field mapping also matches the condition for anatomic image acquisition, and even for breath-controlled RT. In addition, it is important to acquire all echoes within a single breath-hold to avoid the liver position discrepancy across multiple breath-holds (26, 27).

Figure 8.

Magnitude (top row) and phase images (bottom row) with (left column) and without (right column) breath-holding. Note phase discrepancies caused by motion in (D) compared with that in (C) and the image quality degradation in the magnitude image in (B) compared with that in (A). The magnitude and phase images are windowed identically for both scans. The contours depict outlines of the liver.


To comfortably accommodate the breath-hold lengths of the majority of patients, acquisition of images with the triple-GREs should be accomplished within 18–19 seconds or less. We used several strategies to accelerate the acquisition, including parallel imaging (acceleration factor of 4), high readout bandwidth (1955 Hz/pixel), and modestly low spatial resolution. One of the consequences of using parallel imaging and high readout bandwidth is a lower SNR in the acquired images. However, our phantom and patient experiments indicate that the SNR realized in our investigation did not have a substantial impact on the quality of the reconstructed field map because of the noise suppression provided by the regularized method. In addition, the slow spatial variation of the field map is suitable for the regularized method. Another benefit of using the high readout bandwidth is to reduce the susceptibility-induced geometric distortion to a negligible level (∼0.5 mm) in both magnitude and phase images, eliminating the need to remove the geometric distortion using more computationally intensive joint reconstruction methods (8). Finally, the modestly low spatial resolution (voxel size = 3.5 mm3) had little impact on the field map estimation because of the slow variation in the field map.

To achieve a short acquisition time, we designed the first and third echoes to have fat signals in-phase, which requires an ΔTE of 2.3 milliseconds at 3 T, and inserted a second echo between them. However, this is only possible if the gradient polarity of the second echo is reversed compared to the first and third echoes. The gradient polarity mismatch results in a linear ramping in the phase image of the second echo (2224). If not corrected, this phase error can cause serious artifacts in the initial field map from the first 2 echoes, which are propagated into the final field map. We found that the phase error could be corrected by a linear approximation similar to the technique for Nyquist ghost artifact correction in echo planar imaging acquisition (28, 29). In addition, the phase error is largely not subject-dependent, and therefore, it can be estimated once and reused for a particular acquisition with only small residual errors. The small residual phase error in the initial field map estimated from the first 2 echoes (having the reversed gradient polarities) causes no problem in the final field map. The advantage of inserting the second echo between the 2.3 milliseconds interval is to have a short ΔTE between the first 2 echoes, which provides a large dynamic range, compared to that of the first and third echoes, to avoid phase-wrapping in the initial field map estimate (which would result in the estimate being trapped in the wrong local minima). Although the initial field map is not free from the chemical-shift effect, the final field map estimation, based upon the first and third echoes, overcomes this challenge. In addition, noise in the initial field map estimate does not affect the final field map, and therefore, there is no need to filter or regularize the initial field map.

Our initial study that assessed the extent of ΔB0 in the liver showed that the ΔTE of 1.15 milliseconds between the first 2 echoes provided a sufficient dynamic range for the B0-inhomogeneity field at the lung–liver interface. The maximum B0-inhomogeneity field at the liver dome rarely exceeded ±434 Hz with good shimming. Therefore, we optimized our triple-GRE protocol of the liver accordingly. For the subsequent study of the liver, we found the maximum ΔB0 rarely exceeded 300 Hz. However, for the body sites that have larger susceptibility effects, for example, head and neck, it may be necessary to further reduce the time difference of the 2 echoes by increasing the gradients to avoid potential phase-wrapping errors. This would result in costs to SNR. Phase-wrapping of the initial field map could also be avoided with the use of dynamic multicoil shimming that compensates for the bulk of the patient-induced B0 inhomogeneity, hence requiring only estimation of the residual inhomogeneity that should have a much lower frequency range (30).

Our method for field map estimation has no need to fit fat and water components, thus improving robustness over other techniques. In contrast, the Dixon-based methods that simultaneously estimate the field map and fat/water components from 3 echoes (14, 19) can have errors >10% in predicting the lipid concentration (16, 18). This error can propagate into the estimated field map, introducing additional artifacts. In addition, no need for a complex fat and water component model leads to high computationally efficiency for the field map estimation and avoids the additional overhead and potential inaccuracies otherwise present in the joint estimation techniques.

In summary, our proposed method estimates undistorted field maps without phase-wrapping and chemical-shift artifacts. In addition, our method depends on neither prior knowledge nor modeling of fat and water components. Our results from phantom and patient studies show that our method is robust to noise and can handle alternating gradient polarities. Our method is currently being validated in an ongoing patient study. As part of our future work, we will investigate its potential use in other body sites susceptible to motion artifacts and/or with substantial fat components, such as the pelvic region.


[2] Abbreviations:


Magnetic resonance imaging


radiation treatment


computed tomography


magnetic resonance




poor signal-to-noise ratio


repetition time


echo time






volume interpolated breath-hold examination sequence


This work was supported, in part, by NIH RO1 EB016079 and RO1 CA132834.

Disclosures: No disclosures to report.

Conflict of Interest: None reported.


    Johansson A, Karlsson M, Nyholm T. CT substitute derived from MRI sequences with ultrashort echo time. Med Phys. 2011;38(5):2708–2714.
    Hsu SH, Cao Y, Huang K, Feng M, Balter JM. Investigation of a method for generating synthetic CT models from MRI scans of the head and neck for radiation therapy. Phys Med Biol. 2013;58(23):8419–8435.
    Dowling JA, Sun J, Pichler P, Rivest-Hénault D, Ghose S, Richardson H, Wratten C, Martin J, Arm J, Best L, Chandra SS, Fripp J, Menk FW, Greer PB. Automatic substitute computed tomography generation and contouring for magnetic resonance imaging (MRI)-alone external beam radiation therapy from standard MRI sequences. Int J Radiat Oncol Biol Phys. 2015;93(5):1144–1153.
    Chang H, Fitzpatrick JM. A technique for accurate magnetic resonance imaging in the presence of field inhomogeneities. IEEE Trans Med Imaging. 1992;11(3):319–329.
    Jezzard P, Balaban RS. Correction for geometric distortion in echo planar images from B0 field variations. Magn Reson Med. 1995;34(1):65–73.
    Baldwin LN, Wachowicz K, Fallone BG. A two-step scheme for distortion rectification of magnetic resonance images. Med Phys. 2009;36(9):3917–3926.
    Wang H, Balter J, Cao Y. Patient-induced susceptibility effect on geometric distortion of clinical brain MRI for radiation treatment planning on a 3T scanner. Phys Med Biol. 2013;58(3):465–477.
    Matakos A, Balter J, Cao Y. Estimation of geometrically undistorted B(0) inhomogeneity maps. Phys Med Biol. 2014;59(17):4945–4959.
    Huang KC, Cao Y, Baharom U, Balter JM. Phantom-based characterization of distortion on a magnetic resonance imaging simulator for radiation oncology. Phys Med Biol. 2016;61(2):774–790.
    Moon-Ho Song S, Napel S, Pelc NJ, Glover GH. Phase unwrapping of MR phase images using Poisson equation. IEEE Trans Image Process. 1995;4(5):667–676.
    Cusack R, Papadakis N. New robust 3-D phase unwrapping algorithms: application to magnetic field mapping and undistorting echoplanar images. Neuroimage. 2002;16(3 Pt 1):754–764.
    Jenkinson M. Fast, automated, N-dimensional phase-unwrapping algorithm. Magn Reson Med. 2003;49(1):193–197.
    Bioucas-Dias JM, and Valadao G.. Phase unwrapping via graph cuts. IEEE Trans Image Process. 2007;16(3):698–709.
    Dagher J, Reese T, Bilgin A. High-resolution, large dynamic range field map estimation. Magn Reson Med. 2014;71(1):105–117.
    Yu H, Reeder SB, Shimakawa A, Brittain JH, Pelc NJ. Field map estimation with a region growing scheme for iterative 3-point water-fat decomposition. Magn Reson Med. 2005;54(4):1032–1039.
    Yu H, McKenzie CA, Shimakawa A, Vu AT, Brau AC, Beatty PJ, Pineda AR, Brittain JH, Reeder SB. Multiecho reconstruction for simultaneous water-fat decomposition and T2* estimation. J Magn Reson Imaging. 2007;26(4):1153–1161.
    Hernando D, Haldar JP, Sutton BP, Ma J, Kellman P, Liang ZP. Joint estimation of water/fat images and field inhomogeneity map. Magn Reson Med. 2008;59(3):571–580.
    Yu H, Shimakawa A, McKenzie CA, Brodsky E, Brittain JH, Reeder SB. Multiecho water-fat separation and simultaneous R2* estimation with multifrequency fat spectrum modeling. Magn Reson Med. 2008;60(5):1122–1134.
    Funai AK, Fessler JA, Yeo DT, Olafsson VT, Noll DC. Regularized field map estimation in MRI. IEEE Trans Med Imaging. 2008;27(10):1484–1494.
    Pineda AR, Reeder SB, Wen Z, Pelc NJ. Cramér-Rao bounds for three-point decomposition of water and fat. Magn Reson Med. 2005;54(3):625–635.
    Reeder SB, Pineda AR, Wen Z, Shimakawa A, Yu H, Brittain JH, Gold GE, Beaulieu CH, Pelc NJ. Iterative decomposition of water and fat with echo asymmetry and least-squares estimation (IDEAL): application with fast spin-echo imaging. Magn Reson Med. 2005;54(3):636–644.
    Ahn CB, Cho ZH. A new phase correction method in NMR imaging based on autocorrelation and histogram analysis. IEEE Trans Med Imaging. 1987;6(1):32–36.
    Bruder H, Fischer H, Reinfelder HE, Schmitt F. Image reconstruction for echo planar imaging with nonequidistant k-space sampling. Magn Reson Med. 1992;23(2):311–323.
    Hu X, Le TH. Artifact reduction in EPI with phase-encoded reference scan. Magn Reson Med. 1996;36(1):166–171.
    Kirkpatrick JP, Wang Z, Sampson JH, McSherry F, Herndon JE 2nd, Allen KJ, Duffy E, Hoang JK, Chang Z, Yoo DS, Kelsey CR, Yin FF. Defining the optimal planning target volume in image-guided stereotactic radiosurgery of brain metastases: results of a randomized trial. Int J Radiat Oncol Biol Phys. 2015;91(1):100–108.
    Balter JM, Dawson LA, Kazanjian S, McGinn C, Brock KK, Lawrence T, Ten Haken R. Determination of ventilatory liver movement via radiographic evaluation of diaphragm position. Int J Radiat Oncol Biol Phys. 2001;51(1):267–270.
    Dawson LA, Brock KK, Kazanjian S, Fitch D, McGinn CJ, Lawrence TS, Ten Haken RK, Balter J. The reproducibility of organ position using active breathing control (ABC) during liver radiotherapy. Int J Radiat Oncol Biol Phys. 2001;51(5):1410–1421.
    Buonocore MH, Gao L. Ghost artifact reduction for echo planar imaging using image phase correction. Magn Reson Med. 1997;38(1):89–100.
    Buonocore MH, Zhu DC. Image-based ghost correction for interleaved EPI. Magn Reson Med. 2001;45(1):96–108.
    Juchem C, Nixon TW, McIntyre S, Boer VO, Rothman DL, de Graaf RA. Dynamic multi-coil shimming of the human brain at 7 T. J Magn Reson. 2011;212(2):280–288.


Download the article PDF (3 MB)

Download the full issue PDF (137.22 MB)

Mobile-ready Flipbook

View the full issue as a flipbook (Desktop and Mobile-ready)