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:
[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:
[1]
with gk:
[2]
And D:
[3]
We can solve expand equation [1] in terms of D:
[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:
[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):
[6]
Lastly, let’s define a Y matrix for the left side of equation [4]:
[7]
Thus we can express equation [4] with the following 6 directions as:
[8]
With exactly six directions, there is an exact analytic solution and can be solved using standard methods such as the Cramer’s rule:
[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.
