Distortion correction and robust tensor estimation for MR diffusion imaging

Abstract

This paper presents a new procedure to estimate the diffusion tensor from a sequence of diffusion-weighted images. The first step of this procedure consists of the correction of the distortions usually induced by eddy-current related to the large diffusion-sensitizing gradients. This correction algorithm relies on the maximization of mutual information to estimate the three parameters of a geometric distortion model inferred from the acquisition principle. The second step of the procedure amounts to replacing the standard least squares-based approach by the Geman–McLure M-estimator, in order to reduce outlier-related artefacts. Several experiments prove that the whole procedure highly improves the quality of the final diffusion maps.

Introduction

There is currently considerable interest in the use of MRI for imaging the apparent diffusion of water in biological tissues (Le Bihan et al. 2001). The physical process underlying this diffusion is the random walk motion of the molecules in a fluid: due to thermal agitation, the molecules are constantly moving and colliding with neighbors. When the fluid is embedded into the complex geometry of biological tissue, however, the collisions with cell membranes and macromolecules and the restriction to various compartments highly influence this process. Hence, by probing the microscopic motion of tissular water, diffusion-weighted imaging provides a unique in vivo tool for studying the structure of biological tissue. In particular, this imaging modality gives access to various information about the brain microstructures that could be used to improve the interpretation of functional imaging studies.

Diffusion-weighted MRI relies on a phenomenon that could have been considered as a technical problem only. When protons are placed into a static magnetic field (B0), they begin to precess (i.e. their magnetic vector rotates around B0). A key point underlying MRI is the linear relationship between the local field strength and the frequency of the proton precession motion, which is the frequency of the signal produced in the receiving antenna. Adding a gradient to B0 encodes the spin localization into this frequency, which leads to images after Fourier transform. Without special preparation, the protons (spins) precessing in a static magnetic field (B0) do not produce signal in the receiving antenna (in x–y plane) because of lack of coherence between the individual precessions (they are all out of phase and hence have no net transverse component). By applying a 90 degree radio-frequency (RF) pulse, the frequency of which matches the frequency of precession of protons, the spins can be made to be in phase and have a net transverse component, producing signal in an antenna. After the 90 degree RF pulse the spins will again go out of phase, mainly because of the effect of external field inhomogeneities. For static spins, the dephasing caused by external field inhomogeneities can be eliminated with a 180 degree pulse leading to what is called a spin echo. This is not possible for spins undergoing diffusion because they are not static (their position fluctuates randomly because of the random character of the thermal spin motion). The result is diffusion-related signal attenuation.

While with standard MR sequences, the diffusion-related signal attenuation is negligible, diffusion imaging sequences increase this effect with the addition of two strong diffusion sensitizing gradient pulses (Stejskal and Tanner, 1965). These additional gradients increase the attenuation of the signal produced by the spins that move along the gradient direction. Within a simple isotropic medium like a glass of water, the attenuation is related to an exponential of the medium property called the diffusion coefficient D (the standard coefficient of Fick’s law). By using, for instance, an image without diffusion weighting and one diffusion-weighted image, we can calculate a D value for each voxel. As a consequence of their spatial structure, however, many substances and biological tissues exhibit anisotropic diffusion behavior: the computed diffusion coefficient depends on the direction of the sensitizing gradient. Therefore, when anisotropy of the 3D diffusion process is of interest, for instance for fiber bundle tracking (Poupon et al., 2001), a symmetric diffusion tensor D has to be calculated for each voxel from a series of diffusion-weighted volumes (Basser et al., 1994a, Basser et al., 1994b). Each such volume is acquired with a different applied diffusion-sensitizing gradient (Stejskal and Tanner, 1965). These gradients are applied in order to vary a symmetric matrix b (s/mm²) that depends on the gradient direction, strength and timing (Mattiello et al., 1994). The diffusion-sensitizing gradient affects the signal intensity of any given voxel in a manner that can be described by the linear equation

$$\text{ln}\text{S(b)=}\text{ln}\text{S(0)−D}{}_{\mathrm{xx}}\text{b}{}_{\mathrm{xx}}\text{−2D}{}_{\mathrm{xy}}\text{b}{}_{\mathrm{xy}}\text{−2D}{}_{\mathrm{xz}}\text{b}{}_{\mathrm{xz}}\text{−D}{}_{\mathrm{yy}}\text{b}{}_{\mathrm{yy}}\text{−2D}{}_{\mathrm{yz}}\text{b}{}_{\mathrm{yz}}\text{−D}{}_{\mathrm{zz}}\text{b}{}_{\mathrm{zz}}\text{,}$$

where S denotes the signal of the selected voxel. When a sufficient number of different b matrices is used (related to at least six different gradient directions), the diffusion tensor D can be estimated. Such calculations are simple if each voxel in the different volumes represents the same point in the anatomy of the subject, but can be impractical if different volumes of the series are distorted relative to each other. Diffusion-weighted images, however, are often acquired using echo-planar imaging (EPI), to reduce acquisition time and artefacts related to physiological motions. Unfortunately, this fast acquisition scheme is highly sensitive to eddy currents induced by the large diffusion gradients (Haselgrove and Moore, 1996). These eddy currents can cause significant distortions in the phase-encoding direction where the image bandwidth is quite low (see Fig. 1). Since the degree and nature of this artefact typically vary both with the strength and orientation of the diffusion-sensitizing gradient, distortions can dramatically change the estimated diffusion tensor.

The methods for reducing the effects of eddy currents may be divided into three categories. The first one simply consists of modifications of the gradient sequences (Alexander et al., 1997). This approach, however, seems insufficient to get completely rid of artefacts. A second family of approaches, which rely deeply on MR physics, require additional experimental data (Jezzard et al., 1998, Horsfield, 1999). Since the eddy-current distortions do not rely on the subject’s head geometry, these cumbersome additional acquisitions can be done on phantoms only during a calibration operation. Unfortunately, the obtained correction scheme has to be updated on a regular basis because of some slow variations of the magnet (Bastin and Armitage, 2000).

The last kind of approaches are purely retrospective and can be considered as registration methods. They use a distortion geometric model inferred from the acquisition principle, which leads to estimate a few parameters using a standard similarity measure like cross-correlation (Haselgrove and Moore, 1996, Calamante et al., 1999, Bastin, 1999). Such simple similarity measures, however, are not sufficient to perfectly take into account the complex dependencies embedded in Eq. (1) (Bastin and Armitage, 2000). In this paper we propose to estimate the few parameters of the distortion geometric model using the mutual information as similarity measure in order to achieve a more robust correction scheme (Maes et al., 1997, Wells III et al., 1996).

It should be understood that EPI imaging leads to another kind of distortion induced by susceptibility artefacts (Jezzard and Clare, 1999) that is not addressed in this paper because it is not dependent on the diffusion gradient. These distortions are non-linear and depend on the subject’s head geometry. They have to be corrected to relate functional MRI experiments performed with EPI with standard high resolution anatomical images. Several registration schemes have been proposed for this purpose using either a free deformation model with a high number of parameters (Hellier and Barillot, 2000), or a more constrained model taking into account a priori knowledge about the main distortion direction and a difference of squares based similarity measure (Kybic et al., 2000).

This paper proposes a second improvement of the standard calculation of the diffusion tensor D. The linearity of Eq. (1) usually leads to a least squares-based regression method (Basser et al., 1994a). This approach, however, is not robust to the various kinds of noises that can be observed in diffusion-weighted data (Bastin et al., 1998, Basser and Pajevic, 2000, Skare et al., 2000, Anderson, 2001). Non-Gaussian noise can stem for instance from physiological motions (brain beat), subject motions (Atkinson et al., 2000, Clark et al., 2000) or residual distortions. While careful acquisition schemes including cardiac gating (Dietrich et al., 2000) and navigator echo (Butts et al., 1996, Clark et al., 2000) may reduce some of these problems, some weaknesses of the tensor diffusion model lead to other regression problems: each voxel includes several water compartments endowed with different diffusion processes that are mixed up in the data (Clark and LeBihan, 2000). Furthermore, the choice of the gradient directions used by the MR sequence can lead to very different estimation situations (Papadakis et al., 2000). Sophisticated restoration schemes dedicated to diffusion-weighted data are bound to be developed in the future to improve the situation using for instance anisotropic smoothing (Parker et al., 2000). In our opinion, however, these restoration methods cannot be perfect because of the poor quality of the raw data. Hence, in order to overcome the influence of outliers on the tensor estimation, we propose the use of a standard robust M-estimator (Meer et al., 1991). A comparison of the behaviour of both regression methods in the presence of various levels of corrupted data proves the interest of the robust approach. An earlier shorter version of this paper was published in the MICCAI proceedings (Mangin et al., 2001).