Monday, August 4, 2014

DTI Analysis, Steps 1 & 2: Distortion Correction and Eddy Correction

DTI data comes into this world just as we do - raw, unprocessed, unthresholded, and disgusting. Today's parents are much too easy on their children; back in my day, boys and girls were expected to correct their own erratic motion artifacts, to strip their own skulls when too much cranium started sprouting around their cortex, and reach significance through their own diligence and hard work. "No dinner for you until you've gotten rid of those filthy distortion currents on your hands, young man!" was a common saying around my household.

Unfortunately, these days we live in a society increasingly avoiding taking responsibility for their own actions, filled with people more and more unwilling to strip their own skulls. (Even though that would save neuroimaging researchers countless hours of preprocessing.) This trend is the result of several complex cultural and sociological factors, by which I mean: Lawyers. Things have gotten so bad that you could, to take a random example, potentially sue someone posting advice on an unlicensed website if following those instructions somehow messed up your data.

This actually happened to a friend of mine who publishes the nutrition blog tailored for Caucasian males, White Man's Secret. After one of his readers formed kidney stones as a result of following advice to drink a two-liter of Mello Yello every day, the author was sued for - and I quote from the court petition - "Twelve Bazillion Dollars."

In any case, DTI data takes some time to clean up before you can do any statistics on it. However, the preprocessing steps are relatively few, and are analogous to what is done in a typical FMRI preprocessing pipeline. The steps are to identify and correct any distortions during the acquisition, and then to correct for eddy currents.

The first step, distortion correction, can be done with FSL's topup command. Note that this will require two diffusion weighted images, one with the encoding direction in anterior-to-posterior direction (AP) and one encoded in the posterior-to-anterior direction (PA). You will also need a textfile containing information about the encoding directions and the total readout time (roughly the amount of time between the first echo and the last echo); for example, let's say we have a text file, acqparams.txt, containing the following rows:

0 -1 0 0.0665
0 1 0 0.0665

Here a -1 means AP, and 1 means PA. The last column is the readout time. If you have any questions about this, ask your scanner technician. "What if he's kind of creepy?" I don't care, girl; you gotta find out somehow.

While you're at it, also ask which image from the scan was the AP scan, and which was the PA scan. Once you have identified each, you should combine the two using fslmerge:

fslmerge -t AP_PA_scan AP_scan PA_scan 

Once you have that, you have all you need for topup:

topup --imain=AP_PA_scan --datain=acqparams.txt --out=topup_AP_PA_scan

This will produce an image, topup_AP_PA_scan, of where the distortions are located in your DTI scan. This will then be called by the eddy tool to undo any distortions at the same time that it does eddy correction.

The only other preparation you need for the eddy correction is an index list of where the acquisition parameters hold (which should usually just be a 1 for each row for each scan), and a mask excluding any non-brain content, which you can do with bet:

bet dwi_image betted_dwi_image -m -f 0.2

Then feed this into the eddy command (note that a bvecs and bvals image should be automatically generated if you imported your data using a tool such as dcm2niigui):

eddy --imain=dwi_image --mask=betted_dwi_image --index=index.txt --acqp=acqparams.txt --bvecs=bvecs --bvals=bvals --fwhm=0 --topup=topup_AP_PA_scan --flm=quadratic --out=eddy_unwarped_images

Note that the --fwhm and --flm arguments I chose here are the defaults that FSL recommends.

The differences before and after can be seen in fslview:




Of course, the distortion correction step isn't strictly necessary; and you can still do a much simpler version of eddy correction by simply selecting FDT Diffusion from the FSL menu and then selecting the Eddy Correction tab, which only requires input and output arguments. All I'm saying is that, if you're feeling particularly litigious, you should probably be safe and do both steps.



25 comments:

  1. Hi Andy, Do you have any experience or suggestions for the distortion correction when the acq_params file is not available as in the case with the public dataset like ADNI or PPMI. I read that you can use some registration based methods but I didn't find a good tool yet. Do you know about some registration based tools for distortion corrections ?

    ReplyDelete
    Replies
    1. Hey Vikash,

      Unfortunately, no, I don't know of any registration based tools for distortion correction; my understanding is that in order to do it correctly, you need some information from the scanner technicians that isn't otherwise available from the header information in the data.

      If I find out about anything though, I'll let you know. There are other DTI analysis packages, such as TORTOISE, that might be worth looking into (although I haven't used it yet).

      Best,

      -Andy

      Delete
    2. exploredti

      Delete
    3. There is also Brainsuite which is freely downloadable and it offers registration based distortion correction:

      http://brainsuite.org/processing/diffusion/distortion-correction-and-co-registration/

      Delete
  2. u can use exploredti to do registeration based look at manual

    ReplyDelete
  3. This comment has been removed by the author.

    ReplyDelete
  4. Hello Andy, Can you tell me how long it takes for the eddy command to produce the unwarped images output? My terminal is taking a very long time to process it, and I'm not sure whether it's because I've written something uncorrectly or whether it is actually a long process. Thank you.

    ReplyDelete
    Replies
    1. Hey there,

      It's difficult to say, given that there are many differences between computers doing the processing. If you've written something incorrectly, it should either terminate with an error, or give you crazy output. Either way, I'm assuming that you've processed at least one person by now, so I'm curious how long it took?

      -Andy

      Delete
    2. Thank you for getting back to me! Initially, I terminated the process, because I thought that I did something wrong. However, I did rerun one person, and it ended up taking 1 hour and 20 minutes to run that person. I do have one more question. Our lab initially assumed that we should extract all 7 of our b0 volumes when creating the nodif image from our dwidata file, because we thought that would give us a more precise and clean image. However, entering only the first two b0 values, which are also the only b0 values not interrupted by any other value, is the only measure entered into the terminal which still actually provides the image of a brain in fslview. Adding any other b0 volume to this code: fslroi dwidata nodif 0 2 gives us an error. Therefore, is there no use for all of the b0 values together? Do they not produce better results when extracted all together? Thanks!

      Delete
    3. I've been wondering the same thing about multiple b0 images. I have a dataset where there are 4 b0 AP scans and 1 b0 PA scan. What I've done is just pair each AP b0 scan with the 1 PA b0 scan (so my AP_PA_b0 file has 8 scans in it) and run top-up on that. Not sure if this is the right thing to do, though...

      Delete
  5. Hi Andy,

    I've completed topop successfully, but now I need to run a process through freesurfer called tracula, on this new combined AP/PA file. In order to run tracula I need the correct number of bvals and bvecs. I have these files, but they are for the AP and PA images separately. When I try to run the program, it tells me the number of images doesn't match the bvals and bvecs. I've tried just combining the files, copying AP into the same file as PA and running it that way but I still run into the same error. Do you know how I could get the proper bvals and bvecs for the new combined file?

    Thanks

    ReplyDelete
    Replies
    1. Hi Marissa,

      I'm not sure why that's happening; these days I use TORTOISE to preprocess my DTI data before running the analyses in Tracula, and I haven't had issues with that. Give that program a try and see if it fixes your problem.


      Best,

      -Andy

      Delete
    2. Which export option did you use from TORTOISE? I tend to use the FSL (sorted) one for TRACULA. Haven't had any problems... BUT, you may have to transpose the b-vectors and b-values to get them to play nice.

      Delete
  6. Hi Andy,
    What do you advise me to do to correct motion distortion? My patients have a lot of motion artifact ( because they are elderly). I am going to preprocess DTI data and don´t implement TOP UP so.
    I would like to know if I can use any tool to better correct motion.

    ReplyDelete
    Replies
    1. Hi Danielle,

      If you don't want to use topup, you can use another DTI preprocessing package like TORTOISE. The results are pretty good compared to FSL's preprocessing.

      -Andy

      Delete
  7. Hi Andy,
    Can you tell me how to extract Mean Diffusivity (MD) values of each subject with the help of FSL tool ?

    ReplyDelete
    Replies
    1. Hey there,

      I would check out the featquery tool, either from the practical here:

      http://fsl.fmrib.ox.ac.uk/fslcourse/lectures/practicals/feat1/

      Or the video here:

      https://www.youtube.com/watch?v=SNVt7smHLm8

      It should be similar to extracting parameter estimates like you would with any FMRI experiment; just give it a mask and the maps to extract from, and you should get the average MD across all the voxels in the mask.

      Best,

      -Andy

      Delete
  8. Hi Andy,
    Will you be making a video or how-to tutorial on TORTOISE?

    ReplyDelete
  9. This comment has been removed by the author.

    ReplyDelete
  10. Hi Andy,
    If I am using a multishell DTI scheme with 6 shells each with one b0 and 10 directions all acquired A-P is it still possible or beneficial to use topup? Do I need to reorganize my b-val file and volumes so that topup will extract all 6 b0 values first (such that the bval lists 000000,500,1000, etc)? Or does it require just one b0 followed by the number of volumes (in my case 6)?

    ReplyDelete
    Replies
    1. Hi there,

      I haven't used multishell DTI, so I can't say. However, my understanding is that topup requires images in opposite phase encoding directions - for example, one set of images collected A-P, and another set P-A. I don't think it would be appropriate in your case.

      -Andy

      Delete
    2. Thanks Andy!
      So in referencing my PAR files, is "phase encoding direction" synonymous with "preparation direction" designated by the PAR file?

      Or, must you really look at the FSL time series as you did in your video to visualize the distortions themselves in order to determine this?

      Delete
    3. Rather, if I have one scan out of 6 in the shell with a different calculated readout time in for the acqp file for topup (6 rows reading 0 1 0 x) for anterior to posterior aquisition but one row has a different x value then I should still be able to use TOPUP as one of the aquisitions was different?

      Delete
    4. Not necessarily; my understanding is that the directions need to be opposite to each other, not for there to be different readtimes.

      I'm not an expert in this, so I would double-check with the FSL listserv to make sure.

      -Andy

      Delete
    5. I see. Well thank you for the insight!!

      Delete