Preprocessing of diffusion data

> Preprocessing of diffusion data
This page gives an example of how to use the FSL tools “topup” and “eddy” to preprocess diffusion data.

This is all discussed in greater detail on the FSL wiki, and I recommend you read those pages as well before performing any actual analysis of your own data.

Acquisition:

First of all, we need data acquired with at least two different phase encoding directions: for example, AP and PA, or RL and LR. At a minimum, we need at least two B0 volumes with opposite phase encoding directions. We could either run the entire diffusion acquisition twice with opposite phase encoding, or just rerun the B0 volumes with opposite phase encoding, or run two completely different acquisitions, as long as we wind up with at least two B0 volumes, each with an opposite phase encoding direction.

For our example, we have two identical acquisitions consisting of a single B0 volumes followed by 64 diffusion directions arranged on a sphere. Series 1 has RL phase encoding; series 2 has LR phase encoding.

Preparing for topup:

The first step is to convert our data to nifti format. We are using MRIConvert, and we always make sure we use the latest version. We know that later we’ll need the total readout time for our acquisition, and we can get that now from the MRIConvert “series information” screen: 0.063 s. This information is also output in a text file for each converted series. (This parameter is discussed in more depth on the “field map” section of this wiki.) Next, we need to extract the B0 images from the nifti files and combine them into a single volume:

fslroi series1 b0_rl 0 -1 0 -1 0 -1 0 1
fslroi series2 b0_lr 0 -1 0 -1 0 -1 65 1
fslmerge -a b0_rl_lr b0_rl b0_lr

Now we create an acquisition parameters file. This file has one row for each b0 volume (2 in our example), and four columns. The first three columns indicate the phase encoding directions (+RL, +PA, +HF) and column four is for the total readout time. Our file, which we will name acqparams.txt, looks like this:

1 0 0 0.063
-1 0 0 0.063

The topup FAQ (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TOPUP/Faq) has a helpful illustration of different phase encode directions and proper values for this file.

Running topup

We’re going to use the default parameters and configuration file. We need the output parameters and the unwarped images. We are also going to output the field map, because we’re going to use it to correct some other EPI images later (see here for how we might use this file):

topup --imain=b0_rl_lr --datain=acqparams.txt --config=b02b0.cnf --out=topup_results --iout=b0_unwarped --fout=fieldmap_Hz

Decision time

We now have the topup results (go read the FSL topup page if you want to know what just happened) and a pair of corrected B0 images. There are two ways to proceed from here. One is to use “applytopup” to apply those results to the remaining data. The other is to use the FSL tool “eddy” together with the topup results to correct the data. Although considerably slower, “eddy” will probably do a better job, especially if there is subject motion present. If you decide to use “applytopup” you may want to run the tool “eddy_correct” first. You should NOT run both “applytopup” and “eddy”, however. Since we meet the requirements for “eddy” (sampling on the full sphere) we decide to use that.

Running eddy

We need a brain mask, and we’re going to get it from the undistorted b0 data:

fslmaths b0_unwarped -Tmean b0_unwarped_mean
bet b0_unwarped_mean b0_unwarped_brain -make

Then we need a single 4D nifti file with all the diffusion data. The first volume of this file should be identical to the first b0 volume fed to “topup” earlier:

fslmerge -a diffusion_data series1 series2

We need to combine the two bvec files from series1 and series2 into a single file, and we need to do the same for the bval files. We decide to do this with the linux “paste” command:

paste -d"" series1_bvecs series2_bvecs >> bvecs
paste -d"" series1_bvals series2_bvals >> bvals

We also an index file that tells eddy which line of the “acqparams.txt” file to use for each volume of data. In our example, that will be a file (we’ll call it “index.txt”) that has a single line with 65 1’s followed by 65 2’s, separated by spaces:

indx=""
for ((i=0; i<65; ++i)); do indx="$indx 1"; done
for ((i=0; i<65; ++i)); do indx="$indx 2"; done
echo $indx > index.txt

Finally, we can run eddy:

eddy --imain=diffusion_data --mask=b0_unwarped_brain_mask --acqp=acqparams.txt --index=index.txt --bvecs=bvecs --bvals=bvals --topup=topup_results --out=eddy_corrected_data

After some time (possibly hours), we will have our eddy_corrected_data file, ready to be used in dtifit or whatever diffusion processing pipeline we intend to use.

Category: Tags: