Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2010 July 1.
Published in final edited form as:
Med Image Comput Comput Assist Interv. 2009; 12(Pt 2): 951–959.
PMCID: PMC2895313

Atlas-based Improved Prediction of Magnetic Field Inhomogeneity for Distortion Correction of EPI data


We describe a method for atlas-based segmentation of structural MRI for calculation of magnetic fieldmaps. CT data sets are used to construct a probabilistic atlas of the head and corresponding MR is used to train a classifier that segments soft tissue, air, and bone. Subject-specific fieldmaps are computed from the segmentations using a perturbation field model. Previous work has shown that distortion in echo-planar images can be corrected using predicted fieldmaps. We obtain results that agree well with acquired fieldmaps: 90% of voxel shifts from predicted fieldmaps show subvoxel disagreement with those computed from acquired fieldmaps. In addition, our fieldmap predictions show statistically significant improvement following inclusion of the atlas.

1 Introduction

Echo-planar imaging (EPI) is the standard pulse sequence employed in functional magnetic resonance imaging (FMRI) and diffusion tensor imaging (DTI) studies due to its high temporal resolution. In order to extract meaningful information from EPI data, scientists and clinicians typically register low resolution EPI data to high resolution anatomical images. A major problem in achieving accurate registration is the inherent distortion of EPI data due to B0 field inhomogeneity.

Field inhomogeneity results in distortion because the MR image reconstruction assumes a one to one relation between spatial location and the frequency produced by known linear gradient fields. In practice, hardware constraints result in inhomogenities in the main field as high as several hundred parts per million (ppm). This can be corrected to within several ppm using shims [1]. Intrinsic magnetic susceptibility of biological materials causes additional perturbations in the field. In neuroimaging studies, soft tissue and bone have similar susceptibilities of (χt ≈ −9.1 × 10−6) [2] and (χb ≈ −11.4 × 10−6) [3] respectively, but differ significantly from the susceptibility of air (χa ≈ 0.4 × 10−6) [2]. This results in large perturbations around the air-filled sinuses and subsequent distortion of EPI data in the frontal and temporal lobes of the brain.

To correct for this, MR systems contain a set of room-temperature shims that are wound to produce fields based on the first and second order spherical harmonics [1]. Imperfect shimming and higher-order perturbations, however, result in residual subject-specific inhomogeneity, and subsequent distortion of the acquired images. Several techniques have been developed to correct for geometric distortion, the majority of which rely on the acquisition of magnetic fieldmaps [46]. A fieldmap provides a direct measure of the B0 inhomogeneity at each point in the image. In EPI, the geometric distortion is a voxel-wise translation restricted predominantly to the phase-encode direction, which can be computed directly from the fieldmap. Fieldmap techniques, however, have several limitations. First, they require additional scan time, which may be difficult to accomodate in clinical settings as well as in DTI studies where many directions need to be acquired in the same session. Second, fieldmaps suffer from low signal-to-noise ratios (SNR) at tissue/air boundaries, which can reduce their reliability in areas where the distortion is most severe. Acquisition of a single fieldmap may also be invalid for unwarping EPI data if there are significant effects due to motion or respiration in the timeseries. Finally, fieldmaps are not available in many retrospective studies.

Given these limitations, there have been several efforts to predict fieldmaps from tissue/air susceptibility distributions of the anatomy using magnetic field models [7, 8, 2]. Jenkinson et al. showed that the perturbing field can be predicted with high accuracy when CT data of the subject is used to obtain a tissue/air segmentation. Since whole-head CT is rarely available, Koch et al. proposed using registered CT from a reference subject to obtain the tissue/air susceptibility model. Both methods were limited by requiring CT to distinguish air from bone, which have similar intensities in structural MR but significant differences in magnetic susceptibility. Second, these models cannot account for the shim fields present in the scanner that partially compensate for B0 inhomogeneity during EPI acquisition. Without this additional information, applying the fieldmaps for distortion correction is not possible. In [7] it was shown that tissue/air susceptibility models could be derived from structural MRI by using an intensity-based classifier trained with CT. It was also shown that registration of the EPI and structural MR could be used to search over the unknown shim parameters allowing distortion correction of the EPI that agrees well with results obtained using acquired fieldmaps.

Variability in structural MR acquisitions, however, may limit the efficacy of an intensity-based classifier in cases where the MR intensity properties differ significantly from those of the training data. In [7], CT data sets with MR acquired on the same scanner as the subjects of interest could be used to train the classifier, but this may not be possible in many cases. Limited anatomical information below the brain may also prevent accurate estimation of the perturbing field. Therefore, obtaining more reliable susceptibility models from structural MR is critical for retrospective unwarping of EPI data sets that lack acquired fieldmaps. While previous results predicting fieldmaps from structural MR have shown good agreement with acquired fieldmaps, we hypothesized that improved segmentation methods would result in even greater accuracy.

In this paper we describe a method for using anatomical information from a set of 22 whole-head CT data sets to achieve improved, subject-specific segmentation of structural MR. In this approach a tissue/air atlas is constructed from the CT data to obtain priors on the probability of tissue or air at each location in the anatomy. The corresponding structural MR is used to train a classifier that segments the MR of the subject of interest and this is used as input to a first order perturbation field model to compute a subject-specific fieldmap. The method is evaluated by comparison of predicted fieldmaps and acquired fieldmaps. In addition, the MR classifier can be used to obtain probabilistic bone segmentations from structural MR that show promising agreement with segmented CT.

2 Methods

2.1 Atlas Construction

Automatic segmentation of neuroanatomical structures often relies on the use of probabilistic atlases [913]. These atlases are usually constructed by co-registering collections of manual segmentations or other training data. The atlas functions as a spatial prior to represent anatomical variability within a population and compensate for missing information in structural MR images [13]. Although atlas-based methods have typically been applied to the segmentation of brain structures, in this work, we construct a probabilistic tissue/air atlas from 22 CT data sets. By incorporating spatial information into the MR segmentation, we expect improved tissue/air classification in regions where bone is often mislabeled as air.

We obtained 22 datasets consisting of CT and MRI from 3 sources: the publicly available Retrospective Image Registration Evaluation (RIRE) database (17 neurosurgery patients), the Radiology department at Brigham and Women's Hospital (BWH) (4 neurosurgery patients) and the Zubal head phantom (1 subject) [14]. In the RIRE data, each CT image has 27 to 34 slices, 4 mm thick, matrix=512×512, voxel size=0.65×0.65mm. The T1-weighted MRI was acquired on a Siemens SP 1.5 Tesla scanner. MRI for 8 of the 17 subjects has 20 to 26 axial slices, 4 mm thick, no gap, matrix=256×256, voxel size= 1.25×1.25mm, TE/TR=15/650ms. T1-weighted MP-RAGE for the other 9 subjects had TE/TR=4/10ms, matrix=128×256×256 and FOV=160×250×250mm. In the BWH dataset, the CT spanned 36 slices with 512×512 in-plane voxels of size 0.46×0.46×4.8mm. 3D-SPGR MRI of these patients was obtained: slice thickness=1mm, TE/TR=3/8ms, matrix=512×512, voxel size=0.5×0.5×1mm. The Zubal data consists of CT of the head and neck: 1.2mm isotropic voxels spanning 230 slices with 256×256 in-plane voxels, and T1-weighted MRI with 90 slices of thickness 0.2cm, matrix=256×256, and 25.6×25.6cm in-plane resolution.

For each subject, the CT data was registered to its corresponding MR using 6 degrees of freedom (DOF) and mutual information as the cost function. The MR was registered to standard space using the MNI152T1 atlas as the reference, 12 DOF and normalized correlation ratio. These transformations were then applied to the co-registered CT. All registrations were carried out using FLIRT [15, 16]. Tissue/air labels were obtained by thresholding the CT data in standard space.

The wide variation in the field of view of the CT data results in highly varying amounts of data at each voxel; in particular only the zubal phantom includes observations in the neck. Probabilistic atlases are frequently constructed by counting the number of occurrences of each tissue class at each voxel and normalizing the result to obtain a probability distribution. In the two class situation, this corresponds to ML estimation of the parameter of a binomial distribution, i.e., if k ~ Binomial(n, p), then, the ML estimate of p given observed k is [p with hat] = k/n.

In the case of only one trial (n = 1), or one observation at a given voxel, this will lead to estimates of p that are zero or one, which may be unreasonably certain (i.e. when used as prior probability on tissue class in segmentation, these values would dominate any amount of data in these voxels). One way to avoid this effect, is to put a prior on p; a natural choice is the beta distribution, which is conjugate to the binomial: p ~ Beta(α, β). (The special case of α = β = 1 corresponds to a flat prior and Laplace's rule of succession). We have chosen to use α = β = ε = 0.05, and in this case, the posterior expected value of the parameter is p¯=k+εn+2ε, which avoids, by ε, the probability zero and one cases mentioned above. In addition to the tissue/air atlas, a probabilistic atlas of bone was obtained by segmenting bone from CT data and applying the same binomial model and conjugate prior. This was applied to segment bone from MR, which may be useful in other applications such as calculation of attenuation maps for absorption correction in PET or dose calculation in radiotherapy planning.

2.2 Atlas-based Segmentation

Structural MR was segmented using an MR classifier that incorporates spatially dependent prior information from the probabilistic atlas and MR intensity information (from the subject of interest) to obtain a subject-specific susceptibility model. The classifier was trained using the CT/MR training data described in section 2.1, but applied to segment MR data acquired at a separate site. The accuracy of the segmentations was evaluated by comparing fieldmaps predicted from the atlas-based segmenter to acquired fieldmaps. The fieldmaps were also compared to those predicted using intensity information alone (ie. a spatially constant prior).

We obtained T1-weighted MRI and fieldmaps of 5 subjects previously acquired at Massachusetts General Hospital on a 3T Siemens TimTrio scanner as part of the FBIRN multi-center FMRI study. The MP-RAGE spanned 160 slices, with thickness=1.2mm, matrix=256×256, voxel size=0.86×0.86mm, and TE/TR=2.94/2300ms. The gradient echo fieldmaps had 30 slices, thickness=5mm, matrix=64×64, voxel size=3.44×3.44mm, TE1/TE2/TR=3.03/5.49/500ms.

Rigid (6 DOF) registration of the T1 data to the fieldmap magnitude image was carried out using FLIRT [15] so the predicted fieldmaps would be in the same space as the acquired fieldmaps for validation. The MNI152T1 reference image was registered to the fieldmap magnitude image using 12 DOF and the resulting transformation was applied to the atlas-based probability maps.

The prior probabilities of tissue, P(T[mid ]Xn), and air, P(TC[mid ]Xn), at each voxel were obtained from the co-registered atlas. For each of the 22 T1 images in the training data, the MR was scaled according to the parameter that minimized the Kullback-Leibler Distance to the MR of the subject of interest. The data that showed the minimal distance to the subject of interest was used to compute the likelihood terms, P(Ii[mid ]T) and P(Ii[mid ]TC), using tissue/air labels from the corresponding CT. The posterior probability of tissue was computed and applied to segment the MR of the subject of interest. Intensity-based segmentation using spatially stationary priors computed from normalized intensity histograms of the training data was also carried out for comparison to the atlas-based approach.

Bone segmentations were obtained in similar fashion and evaluated using a ‘leave-one-out’ framework in which one CT was withheld as ground truth and the remaining CT data sets were used to construct the probabilistic bone atlas. A non-linear direct search was performed to solve for the thresholds that maximized the similarity of the estimated bone segmentations to the ground truth segmentation. The similarity was quantified using the dice score, where a value of 0 indicates the volumes have no overlapping voxels and a value of 1 indicates they are exactly the same.

2.3 Fieldmap Estimation

Fieldmaps are predicted from the atlas and intensity-based segmentations using the perturbation field model described in [2]. In this model, a first order perturbation solution of Maxwell's equations is calculated from a tissue/air susceptibility model, where each pixel takes continuous values between 0 (air) and 1 (tissue). For a single voxel, with the B0 field along z, the predicted field is given by: F(x)=(2Gz2)(χ1Bz(0)), where G(x) = (4πr)−1, r= x =x2+y2+z2, Bz(0) is the field strength, and χ1 = 1 is the susceptibility value of the voxel. Due to the linearity of the perturbing field solution, the total field is given by: Bz(1)(x)=xχ1(x)F(xx) where x′ are the source points (locations of the voxel centers) and x is the field point where the field is evaluated.

Current field modeling techniques, including the one described in [2], do not account for the shim fields that reduce the B0 inhomogeneity prior to fieldmap acquisition. Therefore, in order to compare an estimated fieldmap to an acquired one, the shim fields must added to the predicted fieldmaps. This is done by modeling the shim fields using the set of first and second order spherical harmonic basis functions. In addition, a global scaling of the predicted fieldmap must be estimated since the model assumes the magnetic susceptibility throughout the brain, (χt ≈ −9.1 × 10−6) [2], is constant, but this may not be accurate near bone interfaces where both partial volume effects and mis-estimation of segmentation values are most likely to occur. Furthermore, the perturbing fieldmaps are calculated assuming a perfectly homogeneous B0 field, which cannot be achieved in practice due to constraints on the hardware. The fieldmap scaling and shim parameters, θ, can be obtained by least squares fitting to the acquired fieldmap: θ^=arg minθ[BAθ]2 where A = [[B with circumflex], S1, S2, …, S8]. The column vectors B, [B with circumflex], and Si, represent the acquired fieldmap, predicted fieldmap, and shim basis functions, respectively. Once these coefficients are known, the predicted fieldmap with shim, [B with circumflex]s = A[theta w/ hat], can be compared to the acquired fieldmap.

3 Experimental Results

Results of atlas-based segmentation of structural MR is shown in Fig. 1. Fig. 1a shows T1 of the sinus region. Fig. 1b shows the limitations of using the intensity classifier to segment the MR. While it produces reasonable results for many of the voxels in the sinuses, voxels outside this region which are clearly soft tissue or bone are mislabeled with values close to zero. In contrast, using the atlas-based segmenter (as shown in Fig. 1c) achieves similar results for the highly variable subject-specific anatomy within the sinus region, while producing fewer errors in the surrounding area. The intensity and atlas-based segmentations were used as input to the perturbation field model to obtain predicted fieldmaps. The scaling and shim parameters were fit from the acquired fieldmaps as described in section 2.3. The shim fields could then be added to the predicted fieldmaps for comparison to the acquired fieldmaps as shown in Fig. 2. The first column of Fig. 2 shows fieldmaps computed from the intensity-based segmentations, which show significant differences relative to the acquired fieldmaps shown for each subject in column 3. These are especially noticable in areas that have lower signal in MR, such as in the ventricles and major sulci. Fieldmap results from the atlas-based segmentations are shown in the second column of Fig. 2 and show improved agreement with acquired fieldmaps. Quantitative analysis of the absolute error in the B0 field between these images is given in the table in Fig. 2. Since the bandwidth/pixel for the EPI data acquired in this study is 22.3 Hz, 90% of the voxels in the atlas-based fieldmaps show subvoxel error. The mean of these statistics across all five subjects is also shown for both the intensity and atlas-based classifiers. The intensity classifier shows a slight improvement over the results reported by Koch et al [8] for a single subject. The atlas-based classifier out performs both the intensity and Koch methods. Paired t-tests comparing the means of the intensity and atlas-based results shows this improvement is statistically significant (all p-values < 0.05). Results of the segmentation of bone from structural MR for a representative subject is shown in Fig. 3. The CT shown in Fig. 3a can be easily thresholded to segment bone from air and soft tissue. Fig. 3b and Fig. 3c show the results of using the intensity and atlas-based classifiers, respectively. While the intensity classifier has some success in segmenting MR into tissue/air classes, it is much less effective in segmenting bone (Fig. 3b). Inspection of the the atlas-based segmentation (Fig. 3c), however, shows good general agreement with the CT, with a dice score of 0.780 for this subject.

Fig. 1
Results of the Segmentation. The T1 MR for a representative subject is shown in (a). The tissue probability map computed using the intensity classifier (b) shows mis-classification of voxels outside the sinus region where intensities are low in MR. Using ...
Fig. 2
Results of the Fieldmap Estimation. Predicted and acquired fieldmaps for subjects 1-5 are shown in rows 1-5 respectively. Fieldmaps predicted using the intensity classifier (column 1) show significant differences relative to the acquired fieldmaps (column ...
Fig. 3
Results of the Bone Segmentation. Segmentation of bone using the intensity classifier (b) results in significant errors when compared with CT (a), while the atlas-based classifier (c) shows good overall agreement


This research is supported in part by NIH grants: P41RR13218, U54EB005149, U41RR019703, and U24RR021992. Mark Jenkinson was supported by a UK BBSRC David Phillips Fellowship.


1. Clare S, Evans J, Jezzard P. Requirements for room temperature shimming of the human brain. Magn Reson Med. 2006;55:210–214. [PubMed]
2. Jenkinson M, Wilson JL, Jezzard P. Perturbation method for magnetic field calculations of nonconductive objects. Magn Reson Med. 2004;52:471–477. [PubMed]
3. Hopkins JA, Wehrli FW. Magnetic susceptibility measurement of insoluable solids by nmr: magnetic susceptibility of bone. Magn Reson Med. 1997;37(4):494–500. [PubMed]
4. Weisskoff RM, Davis TL. Correcting gross distortion on echo planar images. Soc Magn Res Abstr. 1992;11:4515.
5. Jezzard P, Balaban RS. Correction for geometric distortion in echo planar images from b0 field variations. Magn Reson Med. 1995;34:65–73. [PubMed]
6. Hutton C, Bork A, Josephs O, Deichmann R, Ashburner J, Turner R. Image distortion correction in fmri: A quantitative evaluation. NI. 2002;16:217–240. [PubMed]
7. Poynton C, Jenkinson M, Whalen S, Golby AJ, Wells W. MICCAI LNCS. Vol. 5242. Springer; 2008. Fieldmap-free retrospective registration and distortion correction for epi-based functional imaging; pp. 271–279. [PMC free article] [PubMed]
8. Koch KM, Papademetris X, Rothman D, de Graaf RA. Rapid calculations of susceptibility-induced magnetostatic field perturbations for in vivo magnetic resonance. Phys Med Biol. 2006;51:6381–6402. [PubMed]
9. Fischl B, van der Kouwe A, Destrieux C, Halgren E, Segonne F, Salat D, Busa E, Seidman L, Goldstein J, Kennedy D, Caviness V, Makris N, Rosen B, Dale A. Automatically parcellating the human cerebral cortex. Cerebral Cortex. 2004;14:11–22. [PubMed]
10. Tosun D, Rettmann M, Han X, Tao X, Xu C, Resnick S, Pham D, Prince J. Cortical surface segmentation and mapping. NI. 2004;23:108–118. [PubMed]
11. Pohl K, Bouix S, Kikinis R, Grimson W. Anatomical guided segmentation with nonstationary tissue class distributions in an expectation-maximization framework. ISBI. 2004:81–84.
12. Heckemann RA, Hajnal J, Aljabar P, Rueckert D, Hammers A. Automatic anatomical brain mri segmentation combining label propagation and decision fusion. NI. 2006;33:115–126. [PubMed]
13. Zollei L, Shenton M, Wells W, Pohl K. The impact of atlas formation methods on atlas-guided brain segmentation. MICCAI. 2007;10(WS):39–46.
14. Zubal IG, Harrell CR, Smith EO, Rattner Z, Gindi GR, Hoffer PB. Computerized three-dimensional segmented human anatomy. Med Phys. 1994;21:299–302. [PubMed]
15. Jenkinson M, Smith SM. A global optimisation method for robust affine registration of brain images. Medical Image Analysis. 2001;5(2):143–156. [PubMed]
16. Jenkinson M, Bannister P, Brady M, Smith S. Improved optimisation for the robust and accurate linear registration and motion correction of brain images. NI. 2002;17(2):825–841. [PubMed]