Non-isotropic voxels; What to do and why



A recent question came up on the forum that was concerned with isotropic versus non-isotropic voxel sizes:

Hi, I am trying to process some DTI data using FSL for TBSS or Voxel based comparison. I am following the basic guideline from http://www.cabiatl.com/Resources/Course/tutorial/html/dti.html and this works good.

However, I found my DTI data set is not isotrophic, about 0.9 x 0.9 x 3 mm.
My question is

  1. Does the DTI voxel size matter for TBSS or VBM? Can I run those statistical using this asymmetric voxel data?
  2. If it matters, how can I overcome the problem? I have FSL on my macbook and Matlab.
  3. How about the effect of this asymmetric voxel for tractography?

My response:

That is a great question, I think it is indeed important to use isotropic voxels for DTI. You can run VBM and TBSS on non-isotropic voxels, interpolation will not likely alter your results, although isotropic acquisition might have.

It is a different story for tractography, where it will be important to interpolate to isotropic voxels. Non-isotropic voxels can negatively influence fiber tracking results, eg. you are less likely to find specific smaller/bendier tracts. My personal experience when running different tests, is that higher (isotropic) resolution (even when through interpolation) will improve tractography results. Although you'll have to keep an eye on your computing power, it is easy to "break" your processor/RAM by going overboard with higher voxel resolution which will make visualizations tedious (although that doesn't become a problem till 0.5mm^3 and higher resolutions).

I tried to find some references to support this, and although I know they must be out there, I couldn't quickly find any conclusive quotes. Somehow in most simulations that I found the authors assume isotropic voxels, while most clinical articles use non-isotropic voxels. Which is unfortunate.
Some quotes:
"FA values that are measured in regions containing crossing fibers are underestimated when using nonisotropic DTI."
http://www.ajnr.org/content/28/6/1102.short
"The use of isotropic voxels is recommended to ensure that the accuracy of the tractography scheme is independent of fiber direction."
http://onlinelibrary.wiley.com/doi/10.1002/1522-2594(200010)44:4%3C625::AID-MRM17%3E3.0.CO;2-O/full
"When selecting the acquisition parameters, one should strive for isotropic resolution since different in-plane and between-plane resolutions would lead to differential averaging of fiber orientations, leading to complicated modeling requirements if this is to be accounted for."
http://www.ncbi.nlm.nih.gov/pubmed/22846632
In addition to the above references Alexander Leemans directed my attention to this work:
http://cds.ismrm.org/protected/11MProceedings/files/1945.pdf

What to do:

These answers however do not address how you can alter you voxel size if you indeed would like to change your voxel dimensions to be more isotropic. Sam Hurley gave the great suggestions of using FSL's FLIRT to make this change. The nice thing about this tool is that you can stay within the FSL environment, but the down side it that it assumes you have a suitable reference image available with the correct dimensions. If you do not, you can instead decide to use other tools, like for example the ones developed by Gary Zhang for DTI-TK; SVResample and TVResample. The first will be able to resample scalar images (eg. the DWI images), while the latter is able to resample tensor images. Which tool you should use depends on what phase the data is in and when you want to alter the voxel dimensions.

The trick with either tool is to be considerate of the fact that resampling any dimension of your voxel space will inadvertently influence the size of the brain within the specified dimensions of the image. If you increase the voxel size, e.g. in the x direction from 1mm to 0.5mm, the space occupied by the brain will include twice as many slices. So where your image dimension used to be 248 slices it should now be 512. If you do not accommodate for this change the brain might get chopped off at the edges of the image.
Below is an worked out example of how to calculate the correct image dimension when altering specific image dimensions: 



x = (256x1)/1 = 512
y = (256x1)/0.5 = 512
z = (64x4)/0.5 = 512

The below code is just an example, there are more options available both within this tool and the TVResample tool that you can explore further:

SVResample -in 019_bravo.nii.gz -out 019_final.nii.gz -xs 512 -ys 512 -zs 512 -xv 0.5 -yv 0.5 -zv 0.5
Methods

DTI processing with Brainvoyager



BrainVoyager QX is a powerful neuroimaging software package. It started as a tool for the analysis of anatomical and functional MRI data sets but has evolved over the years into a multi-modal analysis tool for fMRI, DTI, TMS, EEG and MEG data. The software is highly optimized and user friendly running on all major computer platforms including Windows (XP/Vista/7), Linux (i.e. Ubuntu, SUSE, Fedora) and Mac OS X (10.6 or higher).
Software tools

Neuroscience Information Framework for DTI resources




The Neuroscience Information Framework (NIF) is an initiative by NIH blueprint to bring together all knowledge on neuroscience. Of course DTI is part of this. Find everything NIF put together on DTI here:
https://neuinfo.org/mynif/search.php?q=dti&t=registry
Software tools

Calculate the diffusion tensor from diffusion weighted images: A mathematical example




A guest post by Rodrigo Dennis Perea on the computational models for estimating the diffusion tensor from diffusion weighted images: 
"When I first started working with DTI, I had a hard time understanding the derivation of the diffusion tensor from what was acquired on the scanner, the diffusion weighted image (DWI). Once I understood this process a little better I decided to create a brief review. I hope this will help you! It might be lacking information or it might be too implicit so please feel free to ask me anything if you need more detailed information."
                                                                                               
I have to thank Peter B. Kingsley for his very helpful publication on this topic: http://onlinelibrary.wiley.com/doi/10.1002/cmr.a.20048/abstract 
First, it is worth noting that this example will highlight the computations in a single voxel (e.g. x=76, y=64, z=21). Essentially, a similar approach will be applied to every voxel in the entire MRI scan. The figure below illustrates our example where we have six DWI's, and the sample voxel is denoted by a small color coded square: 
Diffusion weighted images

The simplified diffusion matrix will look like this:
Diffusion matrix
Which will translate into this tensor:
Tensor model

Next the Stejskal-Tanner equation is of importance, this equation shows how there will be a reduction in the signal due to the application of a pulse gradient, this change will be proportional to the amount of diffusion: 

Image[1] 
It is important to notice that we will solve for tensor D given the other parameters from the MRI acquisition. Sk is the intensity at the "single voxel" when a specific gk gradient direction is applied (e.g. a 1x3 vector for the x,y,z direction).  b0 is the no diffusion signal and S0 is the intensity value you get from the no-diffusion signal. So lets begin...
A random MRI-DWI sequence (lets say sample_001) is used as an example. Given the following gradient vectors with their specific intensity values at an specific voxel location in the brain, we will derive equation [1]  to solve for the simplest case with only 7 volumes (1 3D volume with no diffusion and 6 3D volumes with pulse gradient directions).
A simple approach: The “H” approach: Six gradient directions (Sk) images “plus” a no-diffusion image (S0)
Given the following gradient direction pulses and intensities with a b-value of 800 s/mm2:
          
Given equation Stejskal-Tanner equation:
Image[1] 
with gk­:
Image [2]                 
And D:
Image    [3]
We can solve expand equation [1] in terms of D:
Image  [4]
Where the subscript “i” denotes each gradient direction pulse ("i" in this case will go from 1 to 6).  One approach to solve this system of equation is to use matrix algebra. Thus let express D as a six-element column vector, d:
Image      [5]
With a six-element row matrix H containing a large M x 6 matrix, where M is the number of gradient directions (in this case M=6):
H_6[6]
Lastly, let’s define a Y matrix for the left side of equation [4]:
Image[7]
Thus we can express equation [4] with the following 6 directions as:
Image    [8]
With exactly six directions, there is an exact analytic solution and can be solved using standard methods such as the Cramer’s rule:
Image   [9]
Using these equation, I was able to create functions in matlab. I also compared my accurate results with FSL dtifit, which gave me similar results. 
I'll be more than happy to share my matlab functions and explain these in more detailed if needed. 
101
Subscribe to: Posts ( Atom )