1. Introduction
Automatic shimming for optimizing magnetic field uniformities is highly desirable in MR spectroscopy. Objects are often heterogeneous and contain intrinsically unshimmable field variations due to rapid susceptibility changes, which can lead to distortions of the lineshape obtained from the volume [1] [2] [3]. Several shimming techniques using volumes of interest (VOIs) have been proposed in order to improve the B0 field homogeneity [3] - [8]. For example, Holtz et al. (1988) used a surface coil [3] and the signal integral of the free induction decay (FID) over a VOI, iteratively, for field optimization [7] - [12]; however, this technique is time-consuming and impractical for many in vivo applications. Moreover, the FID (or the spectral peak amplitudes) is sensitive to changes in shim settings [7] [9] [10] [11] [12] [13].
The use of linear shim coils is highly advantageous in MR imaging [10] [11]. The use of second, or higher-order, shim coils can introduce nonlinear interactions in the B0 field; specifically, whenever the origin of the VOI is offset from the isocenter [9] [10] [11] [12]. The fast automatic shimming technique, by mapping along projections (FASTMAP) [13] [14] [15], has been very effective in improving B0 field inhomogeneity [5] [16] [17] [18] [19]. This method computes the corrective first and second-order shim currents by mapping the B0 field along six projection directions (or columns). FASTMAP, however, is restricted to selected cubic VOIs [5] [18] [20] [21] [22].
FASTMAP works well over reasonably homogeneous volumes with moderate field inhomogeneity [20] [23]. This technique performs well in applications probing smaller volumes (e.g., single voxel spectroscopy) [24], but not larger ones. For example: during human brain imaging studies at high-fields where VOIs are extended into the frontal and inferior brain regions; where off-resonance may be present, or whenever fields rapidly change.
The FASTMAP technique incorrectly assumes that shim coil fields can be fully characterized by a minimal set of spherical harmonics [25] [26]. Therefore, a shimming technique less susceptible to signal voids than projection based methods, and capable of handling arbitrarily shaped VOIs is highly desirable.
In this paper, we follow the same general principals outlined in FASTMAP but propose several improvements. In brief, we propose combining spherical harmonic functions and linear least squares fitting for estimating field inhomogeneity. This method entails the computation of 3D phase images and the determination of first and second-order spherical harmonic coefficients for specific shim currents, by changing the Digital to Analog Converter (DAC) settings, which control voltages across different shim coils. The spherical harmonic calibration constants are then determined by computing the gradients between spherical harmonic coefficients and the DAC values of each coil—followed by a first order correction [see Equation (5)]. Our analysis is numerically robust and completely flexible when selecting VOIs for shimming. A performance analysis comparing our technique with FASTMAP, on a phantom and a human brain, demonstrates how our proposed method outperforms the FASTMAP technique in terms of B0 homogeneity.
2. Methods
2.1. Equipment
All experiments were performed on a 4T whole-body Varian INOVA (Palo Alto, CA) MRI scanner located in Cincinnati, Ohio. The system was equipped with the following resistive shim coils: X, Z, Y (n = 1, m = 1, 0, −1) and second-order [X2 - Y2, ZX, Z2, ZY, XY (n = 2, m = −2, −1, 0, 1, 2)]. A TEM volume head coil was used for RF transmission and reception.
2.2. Imaging
The imaging protocol employed a modified 3D gradient-echo pulse sequence (see FigureS1 in supplementary material), which was used to obtain B0 field maps. Frequency distortion correction (along the read-out direction) was performed on B0 field maps. All acquisitions used a 256 × 256 ×256 mm field of view; a 128 × 64 × 64 acquisition matrix; a 10˚ pulse flip angle; a repetition time (TR) of 16 ms, and echo times of 5.25 ms and 7 ms. The data was acquired in the axial orientation, with a slab-selective pulse used for excitation.
After acquisition, inverse Fourier transformation was performed on the acquired 3D k-space data. Subsequently, 3D phase unwrapping was performed on the resultant phase images as necessary. Frequency maps were then computed from the difference of the two phase images (acquired at different echo times) with the following equation:
(1)
After calculation of the 3D frequency maps, voxels corresponding to the selected VOI were extracted. All image reconstruction steps were performed in Matlab (Mathworks, Natick, MA).
Images were obtained in both a phantom and in-vivo. The phantom was a water sphere with a diameter of 178 mm. In-vivo images of a human head were obtained from a single subject. Consent was obtained with an IRB protocol approved by the University of Cincinatti School of Medicine. The VOI for shimming was defined as the entire spherical phantom and the brain only, respectively (see supplementary material for details). B0 field maps were acquired both prior to, and after, the shimming procedure outlined below.
2.3. Constructing Calibration Tables for Active Shimming
A one-time procedure was performed to construct shim calibration tables for active shimming. B0 field maps were acquired upon each of the system’s 8 shim coils at different shim current levels. Specifically, the shim current was varied from −15,000 to 15,000 by increments of 5000 per acquisition. Thus, 7 field maps were acquired per shim coil. A spherical phantom (d = 178 mm) was used as the reference object for this calibration procedure. After reconstruction of the 3D phase images for each shim coil, and shim current setting, frequency distribution maps were computed. The matrix representation of
is given by:
(2)
where ηnm are the coefficients of spherical harmonics, and Fn,m is the Cartesian spherical harmonic spatial dependence function (see FigureS4). Using the linear least-squares method, the optimized spherical harmonic coefficients of the first- and second-order shim coils over the selected VOI can be estimated. The frequency distributions of all shim coils (at each DAC step) can be projected onto the spherical harmonics by using Equation (2). We assume that the
of each shim coil is linearly varying with the DAC values.
(3)
Here,
is the calibration constant for each spherical harmonic. These
values can be estimated using the following expression:
(4)
The
values for all 8 shim coils are obtained from Figures S4-S6. Finally, the spherical harmonic calibration constants are computed by the gradient between spherical harmonic coefficients and the DAC values of each coil; this can be used to update the DAC settings. The R2 of this linear fit was also computed (see supplementary material for details).
2.4. Correction Procedure for 1st Order Shims
Generally, first order coils should produce orthogonal fields that correspond to first order spherical harmonics. The second order coils could potentially produce fields that correspond to first and second-order spherical harmonics. Therefore, we propose the following correction when computing optimal DAC settings of first-order shims, in order to counter the contributions of second-order shims:
(5)
Here,
is the shim setting of the 1st order m1th degree coil (X, Y, or Z), and
is the setting for the 2nd order m2th degree shim coil for correct shimming of an object.
and
are the 1st order m1th degree, and the 2nd order m2th degree calibration coefficients of coils, respectively. The term
is the contribution of the second-order coil to the first-order spherical harmonics. Multiplying the term
by a proportionality
constant
, then using Equation (4), we can compute an updated
setting. This new setting has effectively subtracted the contributions of the second order coil from the first order coil (or first order spherical harmonics).
3. Results
Spherical harmonic calibration constants and corresponding R2 values of linear fits (for all shims) are tabulated in Table 1 and Table 2, which are used to compute optimal DAC settings for shimming an object. Second-order shims seem to exhibit higher R2 values in spherical harmonic calibration constants for first-order shims (Table 2). R2 values that are ≥0.9 are highlighted in light blue in Table 2. For example, changes in DAC values of the xy coil influence coefficients of some
Table 1. Spherical harmonic calibration constants. Spherical harmonic coefficients corresponding to frequency distribution maps of objects are multiplied by calibration constants to obtain DAC settings for optimal shimming.
Table 2. The R2 values of linear fits for respective shims. Numbers highlighted in green (diagonal) indicate almost perfect correlation with DAC settings. Numbers highlighted in light blue (off diagonal) indicate strong cross influences between coils with changing DAC settings. In other words, a given second order shim coil can produce undesired field components that project onto the entire spherical harmonic coefficient set.
second-order harmonics (i.e., x2 - y2, xz, and zy) in addition to first-order harmonics. On the other hand, both Table 1 & Table 2 suggest that DAC changes in first-order shims are relatively independent and only influence the first three spherical harmonic calibration constants (e.g., A11, A10, A1-1).
3.1. Active Shimming of a Phantom
Figure 1(A) and Figure 1(B) show B0 field distribution in the phantom before and after active shimming. The histograms of magnetic field distributions (over the entire phantom), before and after active shimming, are shown in Figure 2(A) and Figure 2(B). The full width at the half maximum (FWHM) value of the field distribution after active shimming is reduced by approximately 94.8% (Figure 2(B)). Note: Figure 1(B) and Figure 2(B) are similar to what can be achieved with the proposed shimming method, i.e., using Table 1. These results show that our method improves B0 homogeneity significantly within the phantom.
3.2. Comparison of FASTMAP and Corrected B0 Field Maps Using the proposed Method in a Human Brain
B0 maps following FASTMAP and active shimming methods are shown in Figure 3. Significant field inhomogeneity can be observed after FASTMAP shimming (Figure 3(A)). The introduction of first- and second-order field corrections improved B0 homogeneity (Figure 3(B) and Figure 3(C)), although, the prefrontal cortex and regions near the nasal sinus still remained inhomogeneous. The large susceptibility variations made shimming these regions difficult, whenever the VOI includes the whole brain.
Figure 4 shows the FWHMs after respective shim procedures. With FASTMAP, the FWHM is about 127.1 Hz. This value was reduced to 91.9 Hz after optimal first- and second-order shimming which is a 28% improvement in the field homogeneity (Figure 4(B)). This value was further improved (by approximately 38%) after incorporating the corrections shown in Equation (5), (Figure 4(C)).
Figure 1. Comparison of B0 magnetic field homogeneity within the phantom before (A) and after (B) active shimming. The calibration table is used to achieve the optimization of the field homogeneity.
Figure 2. The histograms of the magnetic field variation of the phantom before (A) and after (B) active shimming. A narrower histogram indicates a highly uniform magnetic field over the entire phantom.
Figure 3. B0 field maps of a subject’s brain at the center slice (A) after FASTMAP shimming and (B) after adjusting optimal first- and second-order corrections and (C) after incorporating adjustments given by Equation (5) to optimal first and second-order corrections.
Figure 4. Histograms of the magnetic field variations in the brain corresponding to field maps shown in Figure 3. A narrower histogram indicates a more uniform magnetic field over the entire brain.
4. Discussion and Conclusion
Performance analyses of phantom results and in vivo results of a human brain showed that our proposed method can significantly outperform FASTMAP. When field maps are derived using all data points within a VOI, B0 homogeneity can be improved by countering the contributions, or effects, of higher-order shims on first-order shims. First order shims play significant roles in B0 homogeneity within small VOIs. Accordingly, taking into account the contributions of higher order shims within small VOIs can be important for many MR spectroscopy applications. Specifically, our method highlights the advantage of using spherical harmonic expansion corrections for shimming spherical volumes.
Our method, however, could not improve the magnetic field homogeneity near regions of the nasal sinus to a satisfying degree: these regions are known for significant susceptibility variations. Future research, focusing on combining active and passive shimming, must be pursued in order to further improve field homogeneity in the frontal brain [27]. Combining these two shimming techniques could be very important for high field MR setups which inherently require higher second-order shim fields [8] [28] [29].
Magnetic field gradient pulses can produce eddy-currents in conductive brain regions [17] [30] [31], affecting the accuracy of field map calculations. These effects can be mitigated by fixing the relative timing of gradient pulses immediately preceding excitation pulses or acquisition windows during δ1 and δ2 (see FigureS1 for details).
There may be instances where simultaneous shimming of arbitrary volumes (with differing levels of field uniformity) becomes necessary. For example: to establish a shim over a particular organ, with a tight B0 range, while maintaining a coarser uniformity over the entire abdominal slice to prevent frequency-based fat-suppression techniques from failing. Thus, our method provides greater flexibility and can be advantageous for shimming arbitrary volumes over FASTMAP.
Here, we followed the method of projecting shim maps onto spherical harmonics: an a priori basis set to represent field maps. Due to some arguments suggesting that the use of spherical harmonics may be sub-optimal [22], Webb et al. (1991) used shim maps themselves as basis sets to produce highly uniform B0 fields over large volumes [31] [32]. A performance analysis comparing our techniques with theirs should be the focus of future research.
Supplementary Materials
Inhomogeneous magnetic fields in the MRI scanner can be corrected by adjusting shim coils to produce additional magnetic fields. These shim coils generate unique magnetic field distributions which are modelled using orthogonal spherical harmonic functions [5] [14] [15] [16] [17] [19].
Below, we present the theory and methods to: 1) numerically estimate inhomogeneous magnetic fields by varying shim settings; 2) derive calibration tables, and (3) determine appropriate shim currents for the first and second-order shim coils.
S.1. Modeling the B0 Static Magnetic Field
Assuming a current density of zero (
), the static inhomogeneous magnetic field
in a region of interest is given by Laplace’s equation (S1).
(S1)
The solution to this equation
can be expressed as a sum of spherical harmonics [14] [20] [29].
(S2)
Here,
and
are the spherical coordinates. n and m are integers satisfying the conditions
; n is the order and m is the degree of a given spherical harmonic.
are the coefficients of spherical harmonic functions. The
is Ferrer’s associated Legendre polynomial [26], and
can be expressed in terms of Cartesian coordinates using TableS1.
(S3)
Table S1. Converting first, and second-order, spatially-dependent, spherical harmonic functions to Cartesian coordinates.
S.2. Experimentally Determining
We performed a phantom study using the pulse sequence shown in FigureS1 to compute the B0 field maps.
was computed by comparing two phase images with different echo times.
At each voxel, the relationship between
, phase evolution
, and echo time (DTE) is given by Equation (S4):
. (S4)
Here
is the gyromagnetic ratio in radian/s/T for proton
. Since the phase can only have magnitudes between
, phase unwrapping must be performed on an as needed basis. At each voxel, the distribution of the precessional frequency
is related to
by Equation (S5):
. (S5)
These
maps were computed for DAC values: (A) −15,000, (B) −10,000, (C) −5000, (D) 0, (E) 5000, (F) 10,000, and (G) 15,000.
S.3. Phantom Study
This procedure was repeated on a water phantom to compute frequency distribution maps. FigureS2 and FigureS3 show
maps for the water phantom, for all 8 shim coils.
S.3.1. Computing Calibration Tables
By combining Equation (S3) and Equation (S5) we obtain the following expression
(S6)
Here,
and
represent the 0th and higher-order coefficients of spherical harmonics. The matrix representation of
is given by:
Figure S1. The modified 3D-gradient-echo pulse sequence.
Figure S2.
maps for the first order shim coils: (I) X, (II) Y and (III) Z. These maps correspond to DAC settings of −15,000 to 15,000, moving by increments of +5000 for each step.
Figure S3.
field maps for the second-order shim coils: (I) X2 - Y2, (II) XZ, (III) Z2C, (IV) YZ and (V) XY. These maps correspond to DAC settings of −15,000 to 15,000, moving by increments of +5000 for each step.
. (S7)
where ηnm are the coefficients of spherical harmonics. Using the linear least-squares method, the optimized spherical harmonic coefficientsof the first and second-ordershim coils can be estimated. The frequency distributions of all shim coils (at each DAC step) can be projected onto spherical harmonics by using Equation (S7). We assume that the
of each shim coil is linearly varying with the DAC values, i.e.,
. (S9)
Here
is the calibration constant for each spherical harmonic. These
values can be estimated using the following expression:
. (S10)
The values for all 8 shim coils are computed using tables and Figures S4-S6 shown below. Finally, the spherical harmonic calibration constants are computed by the gradients between spherical harmonic coefficients and the DAC values of each coil. This can be used to update the DAC settings.
Figure S4. (A) The first- and second-order spherical harmonic coefficients for the X coil (Hz∙cm−1 and Hz∙cm−2) at DAC values of −15,000, −10,000, −5000, 0, 5000, 10,000, and 15,000. For each DAC value,
maps were computed and spherical harmonic coefficients were estimated using the least-squares technique. (B) Using linear regression, the spherical harmonic calibration constant (Cnm,g) for the X coil (highlighted in orange) is estimated as the gradient (slope) between spherical harmonic coefficients (yellow) and the DAC values. The linear fit of this regression (R2) is also computed. The procedure is explained in the graph above. The same procedure was repeated for other shims.
The calibration constants and their corresponding R2 are obtained from FigureS4 & FigureS5 and used in Table 1 and Table 2.
Figure S5. The gradient [slope: spherical harmonic calibration constant (Cnm,g)] estimates for Y and Z coils. The linear fits of regression (R2) are computed using figures similar to Figure S4. The same procedure described in Figure S4 was followed.
Figure S6. The gradient [slope: spherical harmonic calibration constant (Cnm,g)] are estimates for X2 - Y2, XZ, Z2C, ZY and XY coils. The linear fits of regression (R2) are computed using figures similar to Figure S4. The same procedure described in Figure S4 was followed.
Abbreviation
B0: Static main magnetic field
FASTMAP: Fast automatic shimming technique, by mapping along projections
DAC: Digital to Analog Converter
n: Order of a spherical harmonic
m: Degree of a spherical harmonic
ηnm: Coefficients of spherical harmonics
: Cartesian spherical harmonic spatial dependence function
: Calibration constant.
R2: Linear fit
: Current density
: Static inhomogeneous magnetic field
: Ferrer’s associated Legendre polynomial
: Phase evolution
ΔTE: Echo time
: Gyromagnetic ratio in radian/s/T.
: Precessional frequency distribution at each voxel
: 0th coefficients of spherical harmonic
ηnm: Coefficients of spherical harmonics.
: Calibration constant for each spherical harmonic.
VOI: Volume of Interest