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 (1–3) 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 (4–9).
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 (10–13). 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 (15–18), 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: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:2] becomes the following equation: 3] holds true for a voxel containing a mixture of fat and water signals, as expressed in Eq. . 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. , 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:
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
yielding a sufficiently short total acquisition time to support imaging during a single breath-hold (short first TE and ΔTE);
providing sufficient dynamic range for mapping the large B0 inhomogeneity field at the liver dome without a phase wrap (short ΔTE);
avoiding phase discontinuity from the chemical-shift effect (fat signals in-phase between the 2 echoes); and
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: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.  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.  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.  using the field map estimated from the first 2 echoes with short ΔTE by Eq. . 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
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.
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 (22–24). This phase ramping was approximated linearly and corrected before the phase-difference maps were calculated from the first 2 echoes for initialization of Eq. .
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.
Evaluation of the Proposed Field Mapping Method for the Whole Liver
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.  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).
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.
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.
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).
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.
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.
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).
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 (22–24). 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.