Alan Seth Barnett^{1}, Mustafa Okan Irfanoglu^{1}, Baxter P Rogers^{2}, Bennett A Landman^{2}, and Carlo Pierpaoli^{1}

Gradient nonlinearity causes systematic errors in the measurement of DTI parameters. These errors can be greatly reduced if the actual fields generated by the gradient coils is known. Although the gradient field maps are known to the manufacturers, many users do not have access to them. We describe a method for measuring the gradient field maps using a set of diffusion weighted images of an isotropic diffusion PVP phantom. We use the field maps to analyze a DTI study of the phantom and compare the results to analysis performed without the field maps. The results show that the method works well.

The signal S^{q} in voxel q of a diffusion weighted image of an isotropic diffusion phantom with is

(1) $$$S^q=S_0^q exp(-D\ tr(\bf b))$$$

where $$$S_0^q$$$ is the signal in voxel q without diffusion weighting, D is the diffusivity of the phantom, and tr(**b**) is the trace of the b-matrix.

For an ideal gradient set, the b-matrix is typically expressed as a dyad

(2) $$$b_{ij}=\kappa G_iG_j$$$

where κ is a constant that depends on the timing of the pulse sequence and G_{i} is the prescribed gradient sensitization in the i(=x,y,z)-direction. The effect of non-linear gradients can be taken into account by writing^{1}

(3) $$$G_i={{\partial B^i}\over {r_j}}G^p_j$$$

where $$$G^p_j$$$ is the prescribed gradient, B^{i} is the field produced by gradient coil i (=x,y,z), and r_{i} = x, y, or z for i=1, 2, or 3, respectively.

The field produced by the gradient coils is

(4) $$$B^i(r,\theta,\phi)=\sum_{l,m} \alpha_{lm}(c_{lm}^ir^lP_{lm}cos(\theta))cos(\phi)+s_{lm}^ir^lP_{lm}(cos(\theta))sin(\phi))$$$

where (r, θ, φ) are the spherical coordinates of an image point, P_{lm}(x) is the associated Legendre polynomial, α_{lm} is a normalization constant, and s_{lm} and c_{lm} are expansion coefficients. Manufacturers also use expansion (4) to describe the fields produced by the gradient coils, but a source of confusion and possible error is different values for the normalization constants.

To estimate the expansion coefficients, we acquire a set of diffusion weighted images and find the values of c^{i}_{lm} and s^{i}_{lm} that minimizes χ^{2}. In addition to the c_{lm} and s_{lm}, best fit values of diffusivity D and unattenuated signal S_{0} also have to be computed. Due to B_{0} inhomogeneity, S_{0}
is a function of position, so it cannot be described by a single
parameter. To reduce the number of parameters and thereby increase the stability of the fit, we divide the diffusion
weighted images by the unattenuated reference image and fit the
attenuation function

(5) $$$A={S\over S_0}$$$.

(6) $$$\chi^2=\sum_{q}(A^{qp}_{measured}-A^{qp}(c_{lm}^i,s_{lm}^i))^2 $$$

where the superscript q labels the voxel and the superscript p labels the image, and A is defined by Eqs 1-5, and $$$A^{qp}_{measured}$$$ is the normalized diffusion image p.

The problem of minimizing χ^{2} is not well defined because Eq 1 depends only on the product of the expansion coefficients and the diffusivity D. We therefore further assume that the gradients are properly calibrated and set the terms that describe the "ideal" gradient to 1.

The minimum number of images required to determine the expansion coefficients c_{lm} and s_{lm} for a single gradient is 2: an unweighted image and an image with diffusion weighting in the desired direction (x, y, or z only). In practice, the contribution of the imaging gradients to the diffusion weighting causes a systematic error in the derived fit parameters. To reduce the effect of the imaging gradients, we acquire a second image with the polarity of the diffusion pulses reversed. To compute expansion coefficients for all three gradients we therefore need 7 images.

1. Bammer R1, Markl M, Barnett A, Acar B, Alley MT, Pelc NJ, Glover GH, Moseley ME. *Analysis and generalized correction of the effect of spatial gradient field distortions in diffusion-weighted imaging.* Magn Reson Med. 2003 Sep;50(3):560-9.

2. Rogers BP, Blaber J, Newton AT, Hansen CB, Welch EB, Anderson AW, Luci JJ,
Pierpaoli C, Landman BA. *Phantom-based field maps for gradient nonlinearity
correction in diffusion imaging.* Proc SPIE Int Soc Opt Eng. 2018 Mar;10573. pii:
105733N. doi: 10.1117/12.2293786. PMID: 29887658; PMCID:
PMC5990280.

3. Carlo Pierpaoli1, Joelle Sarlls1, Uri Nevo1,2, Peter J. Basser1, Ferenc Horkay1 Polyvinylpyrrolidone (PVP) Water Solutions as Isotropic Phantoms for Diffusion MRI Studies, ISMRM 2009, page: 1414

4. United States Patent 9,603,546 Horkay et al. Phantom for diffusion MRI imaging March 28, 2017

Figure 1. tr(**D**) maps of PVP phantom. The first row are axial slices, the second row are sagittal slices, and the third row are coronal slices. The first column are uncorrected, the second column are corrected using the manufacturers b-field model, the third column are corrected using parameters from a B_{0} map, and the fourth column are corrected using the diffusion data. All three correction methods produce similar results. The diffusivity values in the axial slice are very well corrected; the diffusivity variations in the z-direction (horizontal in the second row, vertical in the third row) are over corrected.

Figure 2. Fractional Anisotropy maps. The rows and columns are the same as for Figure 1. FA for an isotropic phantom is 0. In the absence of gradient nonlinearity, noise causes computed FA to increase. The largest FA values in these maps are artifacts due to the meniscus at the top of the phantom. The bright regions near the other edges are artifacts caused by eddy current induced image distortion. As is the case for tr(**D**), the three correction methods produce very similar FA maps, with residual error in off center slices.

Figure 3. tr(**D**) histogram. tr(**D**) is uniform across the phantom. Gradient inhomogeniety causes the peak in the histogram to broaden and have a peculiar shape. The effect of the correction is obvious; the corrected histograms are narrower and more symmetrical. The histograms of the data corrected using different methods to determine the field map lie on top of each other.

Figure 4. FA histogram. The actual FA of the phantom is everywhere 0. Both gradient inhomogeneity and noise cause the histogram to broaden and the peak to move to the right. The effect of the correction is obvious; the corrected histograms are narrower and peak at a smaller value than the uncorrected histogram. The histograms of the data corrected using different methods to determine the field map lie on top of each other.

Figure 5. DEC maps. The rows and columns are the same as Figure 1. For a data set with no systematic errors, the directions of the eigenvectors would be random. Gradient inhomogeneity introduces a preferred direction, creating regions of uniform color. The corrections greatly reduce the extent and brightness of the regions of uniform color. The off center slices (bottom of bottom row, right side of middle row) contain regions where the correction is not perfect.