Preprint of paper to appear in Proc. SPIE Medical Imaging 2006 San Diego, CA, (Feb 2006)

ORMAT

Simulation of Susceptibility-Induced Distortions in fMRI

Ning Xu1, Yong Li1, Cynthia B. Paschal2,3, J. Christopher Gatenby3, Victoria L.Morgan3, David R.Pickens3, Benoit M.Dawant1, J. Michael Fitzpatrick1,3

1Department of Electrical Engineering and Computer Science, Vanderbilt University School of Engineering, Nashville, TN

2Department of Biomedical Engineering, Vanderbilt University School of Engineering, Nashville, TN

3Department of Radiology and Radiological Sciences, Vanderbilt University Institute of Imaging Science, Vanderbilt University, Nashville TN

ABSTRACT

It has been recently proposed that computer-simulated phantom images can be used to evaluate methods for fMRI preprocessing. It is widely recognized that Gradient-Echo Echo Planar Imaging (EPI), the most often used technique for fMRI, is strongly affected by field inhomogeneities. Accurate and realistic phantom images for use by the fMRI community for software evaluation and training must incorporate these distortions and account for the effects of head motion and respiration on the distortions. A method to generate realistic distortions caused by field inhomogeneity for the generation of an fMRI phantom is presented in this paper. Changes in field inhomogeneity due to motion are studied by means of adding motions to the brain model and calculating the induced field map numerically rather than measuring it experimentally. A fast analytic version of an MR simulation is used to generate distorted EPI images based on the calculated field maps. The new generated fMRI phantoms can be used to evaluate processing algorithms for fMRI study more accurately. We can appreciate the importance of distortions for fMRI phantom generation by simulating a distortion-free image and adding distortions afterwards. Validations are performed by comparing the calculated field maps with measured ones. In addition, we show the similarities between a simulated fMRI phantom and real EPI image from our MR scanner.

Keywords: Magnetic Resonance, Neurological, Other (simulaton)

1. introduction

Gradient-echo echo planar imaging (GE-EPI) is widely employed in fMRI due to its high sensitivity to susceptibility variations on the microscopic scale such as those due to blood oxygen level dependent (BOLD) signal and its fast imaging speed. Unfortunately, it suffers from patient-induced field inhomogeneity caused by macroscopic susceptiblity variation at air/tissue and bone/tissue interfaces, which, because of its long TE and signal readout time, results in significant geometric distortion and signal loss. As a result, both accurate registration between anatomical images and EPI images and the detection of BOLD signals are compromised.

Computer-generated phantoms provide a way of evaluating new processing algorithms for fMRI data analysis such as motion correction1. A realistic fMRI phantom should characterize noise, motion, and functional activation. It is important to model image distortions and signal loss due to magnetic field inhomogeneities as well. It is known that the variations in field inhomogeneity depend on the orientation of the head relative to the applied static field2. Thus, it is not correct simply to rotate and translate the measured field map of one image and apply it to differently rotated images in

the fMRI time series. The purpose of this work is to simulate realistic susceptibility induced distortions for fMRI phantom generation and examine the impact of motion on the changes in image distortion.

2. METHOD

Given a known distribution of tissue susceptibility, a realistic field inhomogeneity can be calculated3-7. Susceptiblity artifacts are caused primarily by the susceptiblity difference across air and tissue interfaces. Hence, to obtain an approximate susceptibility distribution, an air-tissue model needs to be generated. The typical series of EPI images exhibits varying orientations of one image relative to another. To generate a realistic phantom series, the orientation parameters of a real series can be estimated or modeled. The same orientations can then be assigned to the phantom series, with the caveat that the induced field maps should be recaculated if new images corresponding to new positions are to be generated from an initial image.

We start with the numerically calculated field inhomogeneity distribution and an initial distortion-free image with the goal to generate an image with induced image distortion and signal loss. A time series can be created after we apply the orientation parameters to the first distortion-free image, recalculate the field maps, and repeat the same process to simulate distortion and signal loss.

We first examine the image formation process, to develop an approximate method to distort the image analytically. A single shot 2D GE-EPI signal at position z1 can be expressed as

1)

where is the gyromagnetic ratio, is the proton density at , is the slice selection profile, A is a constant. are the k-space components, m and n are the corresponding sampling indices, are the sampling time intervals between adjacent k-space points in the x and y direction respectively, TE is the echo time, and is a constant k-space shift caused by system delay. T2, the transverse relaxation time constant, is used here instead of T2* since the additional “*” parts of T2* are explicitly included in the last exponential term in Eq. (1). The sign means that the corresponding exponential term changes sign for opposite k-space trajectory for odd echos and even echos. To simplify our analysis, we assume there is negligible distortion in the slice selection direction and the slice selection profile is a prefect rectangular function. Thus, can be replaced with 1.

The ideal signal would have the form, =, and, as a result, the inverse Fourier transform would get back exactly the proton density weighted by T2 decay. However, there are other factors in Eq. that perturb the result. Fundamentally, it is the introduced extra phase modulation along with the alternating k-space sampling trajectory (expressed by the sign) that cause ghosting, distortions, and signal loss artifacts in EPI images.

The alternating sign changes in the exponential term from odd echos to even echos are reflected in , which cause amplitude and phase discontinuities in the phase encoding direction8,9. This gives rise to the well known “N-over-2” ghosting artifacts in image space. Normally these artifacts can be reduced or minimized using reference scans, and thus will not be simulated in our EPI simulation. Therefore our initial distortion-free image can be simulated using the following equation.

2)

If we have as input a distribution of proton-density and T2 for each tissue, we can generate signals using Eq. . After taking the Fourier transform, we can produce a pure T2-weighted distortion-free image. We can simulate distorted EPI images by running a complete simulation using Eq. . This approach is, however, relatively time consuming for fMRI phantom generation, where large numbers of volumes need to be generated for each time series. We thus run this simulation only once, while generating the rest of images based on this initial undistorted image using the following derived analytical approximation.

Ignoring T2 decay, we will have the signal in a form as follows,

3)

In Eq. , is a phase shift at position that is proportional to the inhomogeneity at that point, which, after the inverse Fourier transform of the signal, results in a shift in the phase encoding direction, y, in image space. The term, , which is absent in SE-EPI image, causes an echo shift for each voxel10 and is the source of intravoxel dephasing. It is not a problem when the image formation is in continual space, where each reconstructed spin will be modulated with an extra position dependent phase after the Fourier transform is taken. The amplitude of the transform image will get the true image back exactly. In the discrete case, however, the image intensity for each voxel is actually the result of a weighted average of the complex image over the voxel, or, more precisely, a convolution between the complex image and a non-delta point-spread function. Since varies with position, spins in different positions within a voxel experience different fields, and, as a result they accumulate different phases during k-space sampling. The convolution will cause an in-plane intravoxel dephasing, which is reflected in the image as signal loss (meaning signal diminution). Because we are dealing with 2D multislice GE-EPI image, the measured image intensity still remains in the k-space with respect to the third dimension. However, because the image intensity is a result of an integral over slice selection direction, also causes phase dispersion in the z direction. The signal loss results from the intravoxel dephasing in all three directions, but the z direction dephasing plays a major role because of the larger voxel dimension compared with the other two dimensions11,12.

Taking the Fourier transform of the signal, the reconstructed image intensity at can be represented as

4)

where represent the shifted voxel position in phase encoding direction, is the Jacobian that expresses the change intensity that accompanies the stretching or shrinking of the image, is the gradient of along the phase encoding direction, is the field of view, are voxel size, and C is a constant.

If we assume that is varies linearly within a voxel, then at voxel can be represented as

5)

where is the constant value at , , , and are the gradients along the three respective directions.

If a rectangular truncation function is assumed in k-space, the point spread function in image space will have the form of a sinc function. If we neglect contributions to a voxel’s intensity arising from neighboring voxels, we find that

6)

where D is a constant. If the further assumptions are made that and are each constant within a voxel, we have

7)

We define

8)

in Eq. as “dephasing factors”, because they are intensity factors that account for the effect of image intravoxel dephasing. The factors are in-plane dephasing factors, and is the through-plane dephasing factor. All three of

them are real numbers in the range from 0 to 1. Given the field inhomogeneity at the center of each voxel, the pixel shift, the Jacobian, and all three dephasing factors can be calculated. The factor can be removed by taking the

magnitude. Thus, we have for the magnitude image,

9)

We can use Eq. to generate EPI images with simulated distortions and signal loss included. Rather than simulating distorted EP images using a complete MR simulation from the former description, which is very time consuming, we use the approximate analytic expressions described in Eq. to calculate the geometric distortion ( relative to ), the intensity change due to stretching or shrinking (), and signal loss caused by intravoxel dephasing for each voxel given its field inhomogeneity ().

For a healthy volunteer subject, MR scans were acquired on a Philips Intera Achieva 3.0 T scanner. A 3D T1 weighted anatomical image with matrix size of 256256170 and voxel size 1x1x1mm3 was obtained. A second volume was acquired using slice selection. We call this volume the “2D anatomical image”. The 2D anatomical image, which is also T1 weighted, has 25625628 reconstruction resolution, voxel size 0.9375mm0.9375mm4.5mm, and a 0.4mm slice gap. A measured field map was obtained based on the well known techniques for mapping the B0-field described in13,14. A single-shot GE-EPI sequence was then performed 127 times for an fMRI time series with 28 slices and 8078 scan resolution (interpolated to 128128), 1.875mm1.875mm4.5mm voxel size with a 0.4mm gap. Presumably, the 2D anatomical image has the same orientation and location as the first EPI images. We then align the 3D anatomical image with the 2D anatomical image.

Two different brain models are then constructed based on the aligned 3D anatomical image. Firstly, it is segmented into CSF, white matter, gray matter, etc. to provide a voxel-based brain model. We then take the proton density and T2 values for these tissues as inputs and run a complete simulation to generate a T2 weighted distortion free image. Secondly, the aligned the 3D anatomical image is segmented into two components air and tissue. This second segmentation when provided with the susceptibility of the components represents an approximate susceptibility distribution of the human head. We treat the susceptibility of air as zero and the susceptibility of all tissue, including bone, as equal to that of water, which is 9.05e-6. From that distribution we calculate the field inhomogeneity map using a previously published numerical method3. The calculated field map is used in Eq. .

In order to show the necessity of field map recalculation when the head is in different orientations, we calculated two different field maps with the head of the subject in differing orientations. For comparison, we also simply rotated the field map from orientation 1 to orientation 2.

3. RESULTS

Figure 1 shows a calculated field map comparision between a rotated field map and the field map of an equivalently rotated human head. The maps are not the same, as shown by the absolute difference shown in Figure 1(c). This result demonstrates that it is not correct simply to rotate and translate the field map of one image and use it to correct other rotated and translated images. The model distribution of air and tissue used for the field map calculation is shown in Figure 2. The distribution can be constructed by segmenting the MR images using a simple threshold. To assess the capability and feasibility of our field map calculation, we show comparisons between a numerically calculated field map and an experimentally measured one for the same volunteer in Figure 3. It shows that the field inhomogeneities at the air-tissue interfaces in the calculated map field agree well with those in the measured field map. A simulated distortion-free EP image, a simulated distorted EP image, and a real EPI volume are shown in Figure 4. We can observe similar image distortion and signal loss around sinuses in the simulated images compared with the real EP images. All calculations and measurements are based on the same subject.