Some of the standard BOLD fMRI pre-processing steps include:
Each is described briefly below. Also described are two common steps that are not part of the standard sequence, but can be added optionally: spatial smoothing and normalization to standard space.
The following paragraph describes the first step in the pre-processing pathway for the installation at the University of Pennsylvania. These steps make use of proprietary General Electric code, and thus the source cannot be distributed. Instead, the VoxBo package contains the skeleton of an IDL procedure that can be replaced with idiosyncratic routines appropriate for your particular scanners and pulse sequences. The remainder of the VoxBo package will work flawlessly if, however, you choose to omit this step – the only requirement is that your data be converted at some point into the TES file format (four dimensional data files with Time, X, Y, Z array order).
At the University of Pennsylvania, our data arrive for processing in k-space (the spatial frequency domain). Thus, one of the first steps of processing is to convert the data from k-space to image space. A procedure called RawToTes, within the VoxBo_RawConvert.pro file performs this step. Several other treatments are applied to the data at this time. These include corrections for phase evolution and removal of fat signals (ask a physicist for details). Importantly, Distortion correction is also performed at this time. Static field inhomogeneities, such as those produced by air-tissue boundaries and fillings, can warp the shape of EPI images. This distortion is corrected during this step using the Chemical Shift Image (CSI) file collected immediately prior to the onset of functional imaging. The end result of the raw conversion step is a TES file that contains the distortion corrected imaging data in image space, in Time, X, Y, Z array order.
During multi-slice echoplanar scanning, each slice is acquired sequentially during the TR. As a result, data from one slice are acquired at a different point in time from data on other slices. The correction applied at this point compensates for the staggered order of slice acquisition. Without correction, the data on one slice will represent a point in time as far removed as 1/2 the TR from an adjacent slice (in the case of an interleaved sequence) (Zarahn et al., 1997b; Aguirre et al., 1998b). A procedure called AcqCorrect, within the VoxBo_RawConvert.pro file performs this step.
The correction works by lagging (shifting to the right) the time-series data on each slice using sinc-interpolation. This results in each time series having the values that would have been obtained had the slice been acquired at the beginning of each TR. To make this clear, consider a neural event (and ensuing hemodynamic response) that occurs simultaneously on two adjacent slices. Values from slice “A” are acquired starting at time zero, simultaneous to the neural event, while values from slice “B” are acquired one second later. Without correction, the “B” values will describe a hemodynamic response that will appear to have begun one second EARLIER on the “B” slice than on slice “A”. To correct for this, the “B” values need to be shifted towards the Right, i.e., towards the last value.
This correction assumes that changes of interest in the data are band-limited (i.e. there is no meaningful information present in the data at a frequency higher than that of the Nyquist). This assumption is support by the study of Josephs et al (1997, NeuroImage) that obtained event-related data at an effective TR of 166 msecs. No physiological signal change was present at frequencies higher than our typical Nyquist (0.25 HZ).
The routine accepts a keyword (GATE) that is used to indicate when the fMRI-gate pulse sequence was used to acquire data. During a gated sequence, all “dead” time within the TR is pushed to the beginning of the TR to create a brief “pause” in the scanner pulses (and thus gradient noise). The acquisition correction is slightly modified in this case to correct for the duration of the pause. The magnitude of this correction is dependent upon the TR, number of slices, magnetic field strength used. A table of corrections is built into the code for the parameters encountered in our lab and an error will be returned if a correction is needed that is not defined. NB: This correction should be the first performed (i.e., before orienting, motion correction, padding, smoothing, etc.). Additionally, it should only be performed once!
VoxBo assumes that functional imaging data are in a standard spatial arrangement. This step arranges the data in the appropriate fashion. The goal is to have the data arranged from inferior to superior, with the following array orders:
|Right - Left||X (right side at zero— the radiologic convention)|
|Anterior - Posterior||Y (posterior at zero)|
|Superior - Inferior||Z (inferior at zero)|
Modifications are performed depending upon the orientation in which the data were acquired (i.e., axial, coronal or sagittal). This information is obtained from the header of the TES file. Additionally, two keywords control the behavior of this routine. The first (INTERLEAVED) is set to correct for the acquisition of slices in an interleaved fashion. Interleaving is performed to reduce cross-talk (RF inhomogeneity) between adjacent slices. The second (INVERTED) indicates that the data were collected in a non-standard order (e.g., from superior to inferior for an axial acquisition).
A procedure called OrientTES, within the VoxBo_RawConvert.pro file performs this step.
Despite their best efforts, subjects move their heads during scanning. This step attempts to realign a series of images of the brain back to a template image. While this effort corrects for misregistration (i.e., first order) artifacts introduced by subject movement, it does not correct for so-called “spin history” effects. The algorithm attempts to minimize the squared difference between a template image and the image to be realigned through the application of six motion parameters (X, Y and Z translations and Pitch, Roll and Yaw rotations). A c program called realign accomplishes this step. REALIGN is a straight port of the motion correction algorithm included in SPM96b, written primarily by John Ashburner.
The template typically used for motion correction is the EPI scout volume acquired immediately after the T1 anatomical images. This is done to bring the functional and anatomical data in as close a register as possible.
The routine iteratively applies the correction until the residual motion estimates are below a specified threshold, or until four iterations have been attempted. By default, this threshold is one-half the voxel width along its smallest dimension (i.e., < 2 mm). The routine iterates over residual translational or rotational motion.
In addition to saving the motion corrected TES files, the routine saves a _MoveParams.ref adjacent to every corrected TES file. This file contains the parameter estimates for the first correction iteration, as well as the number of iterations required for each image. The set of MoveParams files for a given subject can easily be displayed by using the VECVIEW -m procedure at the IDL command line. Conversion to standard space
Using parameters calculated upon the anatomical localizer images, this step will “reslice” the images in a TES file into standard (i.e., Talairach) space. The bounding box and voxel sizes default to appropriate values for our data. The external program norm (written in c/c++) is used here and this code is adapted from the spatial normalization code provided in the SPM package.
This routine smooths the data in space by convolution with a 3D Gaussian kernel. The Full-Width at Half-Maximum (FWHM) of the kernel in all three dimensions is specified in units of voxels. Remember to take into account any anisotropy in voxel size when choosing a smoothing kernel to produce uniform smoothness (in mm) in all dimensions. For example, if you have an in-plane (x and y) voxel size of 3.75 mm and a between plane (z) voxel size of 5 mm, then a reasonable smoothing kernel might have a FWHM (in voxels) of 4x4x3. This would yield equal smoothness in all three dimensions (15 mm, FWHM).
This step thresholds the data to establish which voxels are located within the parenchyma of the brain. This is done so that hypotheses are only tested in voxels in which sufficient signal is present and to reduce the amount of data that must be stored in each TES file. Absolute and standard-deviation methods are available.
A final step calculates and saves a global signal and grand-average power spectrum for each TES file. These data are used to create covariates of no-interest and establish intrinsic noise models during subsequent analyses.