Nothing Special   »   [go: up one dir, main page]
Skip to main content
PMC home page
NIHPA Author Manuscripts logoLink to NIHPA Author Manuscripts
. Author manuscript; available in PMC: 2014 May 23.
Published in final edited form as: Inf Process Med Imaging. 2009;21:664–675. doi: 10.1007/978-3-642-02498-6_55

A Variational Image-Based Approach to the Correction of Susceptibility Artifacts in the Alignment of Diffusion Weighted and Structural MRI

Ran Tao 1, P Thomas Fletcher 1, Samuel Gerber 1, Ross T Whitaker 1
PMCID: PMC4031680  NIHMSID: NIHMS576341  PMID: 19694302

Abstract

This paper presents a method for correcting the geometric and greyscale distortions in diffusion-weighted MRI that result from in-homogeneities in the static magnetic field. These inhomogeneities may due to imperfections in the magnet or to spatial variations in the magnetic susceptibility of the object being imaged—so called susceptibility artifacts. Echo-planar imaging (EPI), used in virtually all diffusion weighted acquisition protocols, assumes a homogeneous static field, which generally does not hold for head MRI. The resulting distortions are significant, sometimes more than ten millimeters. These artifacts impede accurate alignment of diffusion images with structural MRI, and are generally considered an obstacle to the joint analysis of connectivity and structure in head MRI. In principle, susceptibility artifacts can be corrected by acquiring (and applying) a field map. However, as shown in the literature and demonstrated in this paper, field map corrections of susceptibility artifacts are not entirely accurate and reliable, and thus field maps do not produce reliable alignment of EPIs with corresponding structural images. This paper presents a new, image-based method for correcting susceptibility artifacts. The method relies on a variational formulation of the match between an EPI baseline image and a corresponding T2-weighted structural image but also specifically accounts for the physics of susceptibility artifacts. We derive a set of partial differential equations associated with the optimization, describe the numerical methods for solving these equations, and present results that demonstrate the effectiveness of the proposed method compared with field-map correction.

1 Introduction

Echo-planar imaging (EPI), used in virtually all diffusion-weighted MRI acquisition protocols, including diffusion tensor imaging (DTI), assumes a homogeneous static magnetic field, which generally does not hold for head MRI. These field inhomogeneities are caused in large part by spatial variation in the magnetic susceptibility of the objects being scanned. Thus, these susceptibility artifacts are image dependent, and will be different for each subject being scanned.

The effect of static field inhomogeneities is described by a physical model of the distortions in k-space. These distortions result in both a spatial warp of the image geometry as well as a change in the image intensities. The geometric distortion is a nonlinear warping along the phase-encoding direction, which is typically the same for the baseline (b = 0) and diffusion weighted images (DWIs). Conventionally, this direction is along the coronal axis of the patient, and is denoted as the y coordinate. The image intensities are preserved as densities in this transformation, and thus the intensities of EPIs are affected locally, proportional to the Jacobian of the corresponding transformation.

The effects of susceptibility artifacts on clinical images are significant. The geometric distortions can be more than a centimeter, and the effects on EPI-acquired intensities in brain tissue can be more than 100% of the original signal. The example in Figure 2 shows a pair of images, an EPI (baseline from a DTI sequence) of a brain and a coregistered, structural T2 image of the same patient. Here we see the obvious distortion in frontal cortex (and associated changes in intensity), but the entire head is also distorted (shortened) in this case. These artifacts impede accurate alignment of diffusion images with structural MRI, and are generally considered an obstacle to the joint analysis of connectivity and structure in head MRI. This joint analysis is one of the primary motivations for the work in this paper. However, even analysis of DTI alone is impaired by these artifacts, because tensor estimates depend on the distorted intensities of the each of the DWIs (with an inhomogeneous relationship), and subsequent geometric analysis, such as tractography, depends on the distorted geometric relationships between nearby tensors. Thus, susceptibility artifacts are an important problem in the analysis of diffusion weighted MRI.

Fig. 2.

Fig. 2

Open in a new tab

EPI Correction. (a) T2 Image. (b) B0. (c) field map corrected EPI image. (d) Proposed method corrected EPI image.

Susceptibility artifacts can be corrected, to a certain extent, by the acquisition of a corresponding field map, as described in [1]. However, the literature shows [2] and we will further demonstrate in this paper, that field map corrections of susceptibility artifacts are not entirely accurate and reliable, and thus field maps do not produce reliable alignment of EPI to corresponding structural MRI. The reasons for this discrepency are not entirely understood, but noise in the field map reconstructions as well as small inaccuracies in the associated physical model can lead to significant artifacts. Fundamentally, however, any attempt to correct EPIs so that they align with structural data must explicitly account for geometry of the structural data (e.g., T1, T2, PD) in order to overcome inevitable noise and inaccuracies in the physical model as well as distortions in the structural data itself. This is the motivation for the image-driven EPI correction proposed in this paper.

The literature does show several examples of image-driven EPI correction [3,4]. However, the methods in [3] employ low-dimensional, relatively smooth warps, which are not adequate for many of the examples in this paper. The methods presented in [4], while using higher-order deformations, do not explicitly account for the intensity transformations in the signal that result from local expansion or contraction of the corresponding transformation. Thus, these methods are not adquate for the very nonlinear distortions that we have observed in practice, particularly in the frontal cortex.

The proposed method uses an image-based correction for susceptibility art-facts. We formulate the correction as an optimization of a penalty that captures the intensity differences between the Jacobian corrected EPI baseline images and a corresponding T2-weighted structural image. The penalty includes a penalty on the dervative of the transformation, to ensure smoothness in the final warp. We derive the first variation of this expression, including the Jacobian term, and show that this formulation reduces to a novel PDE for intensity-conserving image registration. We present a numerical method that allows the method to solve for the optimal registration in practical times and we show results on phantoms and real head data that demonstrates the effectiveness of the method relative to corrections based solely on the field map.

2 Related Work

There is a significant body of work in the area of correcting susceptibility artifacts in echo-planar MRI. Here we review some of that work, for context, with an emphasis on approaches that are technically and strategically related to the proposed method.

The physical model of susceptibility-related distortion is described by Jezzard et al. [1], and they present a method for correcting these distortions using an associated field map. The correction of EPI distortion by field maps is, in principle, correct. However, a variety of researchers have noted that such field maps are not always available, and, when they are, often fail to completely correct geometric distortions. This was most recently observed by Wu et a. [2], and will be confirmed by results presented later in this paper.

Therefore, a variety of researchers have proposed image-based methods for correcting EPI artifacts. For instance, Kybic et al. [5] present a method that uses a deformable registration, represented by B-splines and minimize the mean-squared-difference between a B0 (baseline) EPI image and a T2 weighted, anatomically correct MRI. They optimize by a gradient descent using analytical derivatives and coarse-to-fine multiresolution framework. They do not account for the significant changes in image intensity that occur, due to the conservation of the EPI signal when deformed, which is particularly important with more recent imaging acquisitions that use higher field strengths. Studholme et al. [3] use a spline-based deformation and include the signal amplitude correction (which is proportional to the Jacobian of the transformation). They use the log of the signal to accentuate the impact of low-amplitude image regions, and optimize mutual information using a gradient descent with a finite difference scheme for the updates. Wu et al. [2] also propose a B-spline registration between EPI and structural MRI, and make quantitative comparisons with field maps.

With modern imaging systems, which employ larger magnetic fields to improve noise characteristics, EPI distortions can be quite severe. Thus, there is a need for EPI-structural registration technologies that incorporate a more versatile set of transformations, such as those described by dense displacement fields. For instance, Hellier et al. [6] describe an elastically-regularized displacement field for EPI-structural registration, based on mutual information, but without a correction for signal conservation. Likewise, Tao et al. [4] describe a dense, regularized warp between structural images and EPI, without a Jacobian correction.

The contribution of this paper is to describe a complete variational approach to the registration of EPI images to structural data which gives a gradient descent on dense displacement fields with a proper correction for the conservation of the EPI signal. The result is novel partial differential equation for signal-preserved image registration with some unique characteristics. We present the derivation, describe a numerical scheme for solving this equation, and demonstrate its effectiveness at EPI susceptibility-artifact correction relative to a state-of-the art field-map correction.

3 EPI Susceptibility Artifact Correction

3.1 Variational Formulation

The correction of susceptibility artifacts, which is the subject of this paper, is a necessarily part of a processing pipeline [2] for DWI processing, which corrects for head motion, eddy current artifacts, and coordinate system inconsistencies between scans. We begin with the methodology for susceptibility correction, and then describe the entire pipeline which we use for the subsequent experiments.

Generally, we denote the transformation between the EPI and structural coordinate systems as xE = ϕ(x), where x = (x, y, z) is a shorthand for vectors, and the “E” subscript indicates the EPI coordinate system. The coordinate transformation associated with field inhomogeneities is usually approximated as a displacement along the phase direction of the scan (e.g. as described in [3]). For most protocols this is along the posterior-anterior axis, Thus we denote, without a loss of generality, this displacement direction as y. The displacement is υ(x, y, z) and EPI image coordinates are

(xE,yE,zE)≈(x,y+υ(x,y,z),z). (1)

The relationship between the uncorrected and corrected EPI images must account for the conservation of signal [3], and thus the expression of the EPI image, IE, in structural coordinates is:

IE−S(x,y,z)=Jϕ(x,y,z)IE(x,y+υ(x,y,z),z), (2)

where Jϕ is the Jacobian of ϕ, which, in this case, has a conventient form:

Jϕ(x,y,z)=1+∂υ(x,y,z)∂dy. (3)

The strategy we use in this paper is to register the B0 image in an DWI acquisition with a T2-weighted structural image. Because these two images are correlated, we contruct a squared-error penality between the T2 and the corrected EPI, and include an elastic penalty on υ. That is,

ɛυ=12∫D[(1+∂υ∂y)IE(x,y+υ,z)−IS(x,y,z)]2dxdydz+λ2∫D|∇υ|2dxdydz,

where we have left out the formal dependency of υ on x to keep the expressions concise. We optimize with gradient descent, which we derive from the first variation. Thus we have

δɛυ=IE(x,y+υ,z)∂∂y[IE(x,y+υ,z)−IS(x,y,z)]−IE(x,y+υ,z)∂∂y[IE(x,y+υ,z)∂υ∂y]−λ∇⋅∇υ.

We introduce a dummy variable, t, and let υ descend the energy ευ as a function of t, and thus we have the PDE which describes the evolution of υ as we optimize the EPI-structural alignment. This gives

∂υ∂t=I∼E∂∂y[IS−I∼E]+I∼E∂∂y[I∼E∂υ∂y]+λ∇⋅∇υ, (4)
=E(υ)+F(υ)+G(x), (5)

where ĨE = ĨE(x,y z) = IE(x, y + υ, z) is a shorthand for the unwarped EPI image, without the Jacobian correction.

We have factored the first variation into three terms; E, F, and G; so that Equation 5 has a particular form that demonstrates its unique nature. The equation differs from many other elastic-registration PDEs in two ways. The first term (on the right hand side), which is often considered as a forcing term that drives the two images to correspond, depends on the derivative of the image residual rather than the image residual itself, as is typical with elastically regularized registration methods. Second, besides the regularization in the third term, there is an extra heterogeneous diffusion in the second term, which is a result of the Jacobian correction. This second term has important implications for the numerical implementation of this PDE.

3.2 Numerical Implementation

The numerical approach is to iterate with discrete time steps, Δt, on the transformation υ until ∂υ/∂t falls below some threshold, and we have reached some approximate steady state. The time steps should be chosen so that the updates in υ are smaller than the grid spacing and the solution does not jump over image data in our attempt to find a minimum. We use a coarse-to-fine multiresolution minimization strategy in order to overcome this computational burden of these relatively small steps and to alleviate the tendency toward local minima.

With this limitation on Δt, the first term in Equation 5, E(υ), can be treated as a reaction term and implemented with finite forward differences. The third term, G(υ), is a diffusion term, but it is stationary and is usually implemented as convolution. The second term, however, must be treated with care. The F term is a nonstationary diffusion. Because this diffusion is proportional to the square of the image intensity, it results in a very large diffusion number for this equation. Thus, the time steps required for this equation to be stable in a forward difference scheme are so small it becomes impractical.

For this reason, we solve for the second term using a semi-implicit method. Here we describe briefly the numerical scheme, in one dimension; that is, the y axis along which the diffusion takes places. Subsequently we describe how the regularization (third term in Eq. 5), which is the only part of this equation that depends on x − z, is solved for the 3D case.

Let υik be the ith grid point of the displacement for the kth time step. We denote the discretely sampled input images, the EPI and the structural images, as ei and si, respectively. Let e∼ik denote the deformed EPI image, which is computed by sampling the discrete image ei at position i=υik using a linear interpolation. We denote the discrete forms of intermediate quantities in Eq. 5 as Ei, Fi, and Gi.

Now we consider the three terms in Eq. 5 one by one. For the first term we have:

Ei≈e∼ik2((si+1−si−1)−(e∼i+1k−e∼i−1k)) (6)

The second and third terms are linear operations on vk+1, which forms the implicit part of the update:

Fi+Gi≈+(λ+e∼ike∼ik+e∼i−1k2)υi−1k+1−(2λ+(e∼ik)2+e∼ike∼i+1k+e∼i−1k2)υik+1+(λ+e∼ike∼i+1k+e∼ik2)υi+1k+1 (7)

We treat this as a linear system, letting the grid quantities without subscripts denote the vector forms of the corresponding grid data and the linear operator in Eq. 7 be A. For the update of υ, we have the following

υk+1=(I−ΔtA)−1(υk+ΔtE). (8)

The matrix (I − ΔtA) is tridiagonal, and we can efficiently solve this system using a two-pass Gaussian elimination. We choose Δt so that max[ΔtE] < 0.4.

For three dimensions, the discrete method is the same, except that we have two extra terms in Eq. 7 that are proportional to second derivatives in x and z, which result from the 3D Laplacian. There are several ways to include these. Because λ is usually small relative to max(Ĩ2), this regularization term is not time-step limiting, and we can include the second derivatives from regularization in the explicit part of the equation, which we do in the results that follow. Alternatively, one could treat all of these second order terms implicitly, and use some type of splitting scheme to solve this system [7].

3.3 EPI-Structural Processing Pipeline

The proposed EPI distortion correction is, in practive, the final step in a combined diffusion and structural image processing pipeline. This full pipeline for aligning a set o DWI images with structural data consists of the following steps:

  1. Preprocess structural images (T1 and T2) to remove skull, correct bias field, normalize intensities, and segment tissue classes (to provide a white matter mask) [8].

  2. Correct diffusion weighted images for eddy currents and head motion [9].

  3. Rigidly align the T1, T2, and segmentation label images to the B0 image based on mutual information metric.

  4. Calculate the deformation field υ using the proposed method.

  5. Concatenate the eddy current transformation and the deformation field υ and apply the resulting transformation to all DWIs. By sampling the original images using the composed transformations, we interpolate the images only once.

4 Results

In this section we present results of the proposed method for EPI distortion correction on real data. The first experiment is applying the correction to phantom data, and the second experiment is to apply the method to diffusion imagery of the human brain. In each experiment we compare the results of our method to field-map correction, which was done using the Statistical Parametric Mapping (SPM) field map toolbox [1,10,11,12] 1.

4.1 Phantom Data Results

We applied our method and field-map correction to a phantom dataset provided with the SPM field map toolbox 2. The phantom is a gel phantom with a spherical geometry. Figure 1 shows a slice from the original phantom dataset. Notice that the geometry of the phantom has been globally squashed as well as a smaller local change on the left. The results from the field-map correction (Figure 1(c)) show that the field map does a good job correcting the global squashing in the image, but it does not do a good job correcting the local “dimple” on the left side. The results from our proposed method (Figure 1(d)) show a correction of both global and local distortions. The field map correction may be unable to handle small scale geometric distortions because the field map must be blurred to remove noise.

Fig. 1.

Fig. 1

Open in a new tab

Phantom data EPI correction. (a) Structural image. (b) squashed EPI image. (c) field map corrected EPI image. (d) Proposed method corrected EPI image.

4.2 Data Acquisition

Head data were acquired on a 3 Tesla Siemens scanner. Diffusion weighted images were acquired with a single-shot spin-echo EPI sequence with high-resolution (2×2×2.5mm3), which was performed using bipolar gradients with dual-echo refocusing to reduce eddy currents. This consisted of one image with b = 0s/mm2 and 12 images with b = 1000s/mm2 with different gradient orientations. An undistorted T2 weighted scan was acquired with a turbo spin echo (TSE) sequence. For field mapping, two gradient-echo images with different echo times (TE = 4 and 6.46 ms) were collected.

4.3 Validation Methods

A common problem in validating results from non-rigid image registration algorithms, is the lack of ground truth. Therefore, we are limited to using indirect measures to establish the reliability of our variational approach. Since a primary motivation for this work is to facilitate joint analysis of diffusion and structural imagery, we will use the fit of the diffusion images to the structural images as a validation measure. In each experiment the result of the correction presented in this paper is compared to the correction from application of a field map. Since our method explicitly minimizes the mean-squared error between the EPI and structural T2 image, we use two independent measures for comparison: visual assessment of edges and mutual information between the DWIs and a T1 structural image. First, we show the results of our EPI distortion correction method alongside the field-corrected image in Figure 2.

Visual Assessment of Edges

As an initial check, we confirm through visual inspection that the DWIs and the corresponding T2 image are well aligned after registration by superimposing the contour of the EPIs onto the T2. Figure 3 shows the resulting edge map of the corrected baseline image overlaid onto the T2 image. Also shown for comparison is the edge map for the field-map corrected image. Notice that the result of the proposed method shows better alignment to the T2 than the field map result.

Fig. 3.

Fig. 3

Open in a new tab

EPI Correction. (a)field map overlay T2. (b) Proposed method corrected EPI image overlay T2.

Mutual information of DWIs and T1

We compared the corrected DWIs with the structural T1 image using mutual information. The T1 image has been coregistered with the T2 image in step 3 of the pipeline described in Section 3.3. Mutual information to T1 was chosen because it gives a good metric of alignment of the diffusion images to the structural image, but it is independent of the registration metric used in the proposed method (mean-squared error to T2). Table 1 gives the results of the comparison for both our method and for field map correction. The table shows a comparison of mutual information of the T1 to corrected B0 image as well as the average mutual information of the T1 to each of the diffusion-weighted images. Notice that the mutual information improves for our method in all cases, while it actually decreases on average for the field-map corrected DWIs.

Table 1. Mutual Information Comparison.
Measure B0 Average DWIs
Before 0.725 0.731
Proposed 0.839 0.814
Fieldmap 0.736 0.718
Open in a new tab

Comparison to Low-Order Deformation

In order to test the need for the high-order deformation (dense displacement field) that we use, we compared the results to a B-spline representation of the deformation. To do this we projected the final deformation field onto a cubic B-spline basis with control point spacings of increasing resolution. These lower-order deformations were then applied to the diffusion images in the same way. The results of a mutual information comparison to the T1 image are shown in Table 2. The high-order deformation of the proposed method achieves the best mutual information, and mutual information increases monotonically with increased resolution. This suggests that the local deformations present in the EPI are best handled with a dense displacement field.

Table 2. B-spline Mutual Information Comparison.
Measure B0 Average DWIs
B-spline (8x10x15) 0.764 0.762
B-spline (16x20x10) 0.781 0.780
B-spline (24x30x15) 0.795 0.786
Open in a new tab

5 Discussion

We have presented a variational image registration framework for correcting the geometric distortions from susceptibility artifacts in EPI, which has application to the alignment of diffusion and structural imaging. We have shown results on phantom data as well as a clinical DTI dataset that our method compares favorably to field map correction. Because of the variability of field maps and the dependence of field distortions on the geometry of the objects being imaged, this method will need to be tested on further data to thoroughly compare it to field maps. Another area of future investigation is the interaction between the various geometric corrections required in DTI, including head motion, eddy current, and susceptibility artifacts. Currently, we perform a head motion and eddy current correction jointly, followed by the EPI distortion correction presented in this paper. Since performing the head motion and eddy current corrections changes the phase-encoding direction, this should be accounted for in the EPI distortion correction. Future work will investigate a joint correction model that includes all three models of geometric distortion inherent in DTI. Finally, another useful improvement would be a formulation that uses mutual information rather than image correlation as an image match function. This would allow registration directly to a T1 image in cases where an undistorted structural T2 image is not available.

Acknowledgments

We would like to thank Molly DuBray and Dr. Janet Lainhart for providing the image data. This work was supported by an Autism Speaks Mentor-based Postdoctoral Fellowship, Grant Number RO1 MH080826 from the National Institute Of Mental Health, and the National Alliance for Medical Image Computing (NAMIC): NIH Grant U54 EB005149.

Footnotes

1

Statistical Parametric Mapping : http://www.fil.ion.ucl.ac.uk/spm/

References

  • 1.Jezzard P, Balaban R. Correction for geometric distortion in echo planar images from B0 field variations. Magnetic Resonance in Medicine. 1995;34:65–73. doi: 10.1002/mrm.1910340111. [DOI] [PubMed] [Google Scholar]
  • 2.Wu M, Chang LC, Walker L, Lemaitre H, Barnett AS, Marenco S, Pierpaloi C. Comparison of EPI distortion correction methods in diffusion tensor mri using a novel framework. In: Metaxas D, Axel L, Fichtinger G, Székely G, editors. MICCAI 2008, Part II LNCS. Vol. 5242. Springer; Heidelberg: 2008. pp. 321–329. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 3.Studholme C, Constable R, Duncan J. Accurate alignment of functional epi data to anatomical mri using a physics based distortion model. IEEE Trans Med Imaging. 2000:1115–1127. doi: 10.1109/42.896788. [DOI] [PubMed] [Google Scholar]
  • 4.Tao G, He R, Poonawalla AH, Narayana P. The correction of epi-induced geometric distortions and their evaluation. Proc ICIP. 2007 [Google Scholar]
  • 5.Kybic J, Nirkko A, Unser M. Unwarping of unidirectionally distorted epi images. IEEE Trans on Medical Imaging. 2000;19:80–93. doi: 10.1109/42.836368. [DOI] [PubMed] [Google Scholar]
  • 6.Hellier P, Barillot C. Multimodal non-rigid warping for correction of distortions in functional MRI. In: Delp SL, DiGoia AM, Jaramaz B, editors. MICCAI 2000 LNCS. Vol. 1935. Springer; Heidelberg: 2000. pp. 512–520. [Google Scholar]
  • 7.Tai XC, Lie KA, Chan T, Osher S. Image Processing Based on Partial Differential Equations. Springer; Heidelberg: 2007. [Google Scholar]
  • 8.Leemput KV, Maes F, Vandermeulen D, Suetens P. Automated model-based tissue classification of MR images of the brain. IEEE Trans on Medical Imaging. 1999;18(10):897–908. doi: 10.1109/42.811270. [DOI] [PubMed] [Google Scholar]
  • 9.Rohde G, Barnett A, Basser P, Marenco S, Pierpaoli C. Comprehensive approach for correction of motion and distortion in diffusion-weighted MRI. Magnetic Resonance in Medicine. 2004;51:103–114. doi: 10.1002/mrm.10677. [DOI] [PubMed] [Google Scholar]
  • 10.Andersson J, Hutton C, Ashburner J, Turner R, Friston K. Modelling geometric deformations in epi time series. NeuroImage. 2001;13:903–919. doi: 10.1006/nimg.2001.0746. [DOI] [PubMed] [Google Scholar]
  • 11.Hutton C, Bork A, Josephs O, Deichmann R, Ashburner J, Turner R. Image distortion correction in fmri: A quantitative evaluation. NeuroImage. 16:217–240. doi: 10.1006/nimg.2001.1054. [DOI] [PubMed] [Google Scholar]
  • 12.Jenkinson M. Fast, automated, N-dimensional phase-unwrapping algorithm. Magnetic Resonance in Medicine. 49:193–197. doi: 10.1002/mrm.10354. [DOI] [PubMed] [Google Scholar]

RESOURCES