# Presentation: 095266

## Diffusion Tensor Estimation Method

- Robert Cox and Daniel Glen
- Scientific and Statistical Computing Core
- NIMH

## DTI Overview

- MRI imaging method that shows relative rates and directions of diffusion corresponding to flow principally along fiber tracts
- Applications for diseases involving white matter including perinatal brain injury, stroke, multiple sclerosis, tumors, ....
- Acquisition of MRI volume followed by at least six or more volumes (typically many more) with an additional gradient applied from non-collinear  directions
- I(q) =  J e –b(q) **.** D
- Where
- 	 I(q) = the image voxel intensity for each gradient q
-      J = the ideal image intensity without applied gradient
-          usually taken as J = I(0)
-      b(q) = “b-matrix” = 2GiGj2(-/3) for the qth encoding gradient
-           = 267.5 rad/ms.mT,
- = time lag between starts of gradient,
- = duration of gradient
-      D = the Diffusion tensor, Diffusion in 6 principal directions,
-           Dxx, Dxy, Dxz, Dyy, Dyz, Dzz
- ln(I(0) / I(q)) = b(q) **.** D
- D = ln(I(0) / I(q) ) x b-1

## Eigenvalue calculations	

- D V =  V
- Where
-    D = diffusion tensor in a symmetric, square matrix form (3x3)
-    V = the eigenvector, a vector corresponding to an orientation (3x1)
-     = the eigenvalue, a scalar constant
- For a 3x3 matrix, there are 3 sets of orthogonal eigenvector and eigenvalue solutions
- det(D – 
- (D – V 
- Solved with f2c converted eispack routine in AFNI using a tridiagonal reduction followed by a QL2 solution for the eigenvalues and vectors

## Fractional Anisotropy (FA)
  FA =

Lattice Index (LI) =  anLIn /  an

Mean Diffusivity 

-   FA =
- Lattice Index (LI) =  anLIn /  an
- Mean Diffusivity 

- General Measures

- ~

-  

- LIn =

## DWI, DT, Eigenvalue Samples

- DWI

- Dxx

- 1

- DWI q1

## DTI images

- Color coded EV1

- A-P green
- L-R red
- I-S blue

## Negative Eigenvalues

- Skare, et. al, Mag. Res. Imaging, 18, 2000

## Negative Eigenvalue Solutions

- Negative values more likely with increasing FA and with increasing noise (Skare, et al, Basser)
- Calculate FA with negative eigenvalues reset to 0 (GE, Hopkins DTI Studio)
- Use another index, e.g. Lattice Anisotropy Index (Basser)
- Spatial smoothing of DWI images (Hahn, et al)
- Temporal smoothing with repeat experiments (Skare, et al)
- Calculate D, limiting D to positive definite matrix (Tschumperlé and Deriche, Mangin et al)

## Direct Least Squares for D

- I(q) =  J e –b(q) **.** D + noise
- Noise – Johnson noise, Eddy currents, motion including periodic beat of CSF with blood flow, partial volume effect
- Find the symmetric non-negative definite matrix D that minimizes the Error functional (cost function)
-    E(D, J) = 1⁄2  w(q) (J e –b(q) **.** D - I(q))2
-                           q
- Generally, estimate D, adjust through a gradient step to find new estimate for D until D converges. Construct gradient step to guarantee D is always positive definite (no negative eigenvalues).

## Initial Estimate of D,  D0

- Estimate for D the traditional way and find eigenvalues, vectors too
- Limit eigenvalues to positive ones by setting 2, eigenvalues to be at least 0.2*1
- Recompute D0 = U  UT
- Where U = [u1, u2, u3] matrix of eigenvectors
-  = diag()

## Compute J (Ideal image voxel value)

- E(D, J) = 1⁄2  w(q) (J e –b(q) **.** D - I(q))2
- E/J =  w(q) (J e –b(q) **.** D - I(q)) e –b(q) **.** D
-    J = ( w(q) I(q) e –b(q) **.** D) /
-           ( w(q) e –2b(q) **.** D)

## Modified Gradient Descent

- Gradient of E (error) with respect to D = F
- F =  w(q) (J e –b(q) **.** D-I(q)) b(q)
- We want to change D to minimize E the fastest. From definition of gradient,
- D/ = -F   (we don’t use this)
- Where  is pseudo-time in the descent.
- But this doesn’t prevent D from becoming non-positive definite, so instead ...

## Modified Gradient Descent

- D/  = -(FD2 + D2F)
- This can be shown to be the fastest descent while remaining positive definite
- We will do the descent with finite steps
- If D/  = -(ND + DN) with N as a constant matrix, it can also be shown
- D() = e(-N) D() e(-N)
- Let N=FD and approximately constant over  step
- D() = e(-FD) D() e(-FD)

## Modified Gradient Descent

- We can replace e-x with the Padé
-  approximant (similar to a Taylor series expansion)
- e-x ~ (1 - x/2) / (1 + x/2)
- Similarly, for a matrix exponential,
- e-M ~ (I – 1⁄2 M) (I + 1⁄2 M)-1
- If we let
- H± = I ± 1⁄2  FD
- Then
- D() = H-() H+()-1 D() H+()-1 H-()
- D() = A() D() A()T
- where A = H- H+ , which will always be symmetric and positive definite

## , pseudo-time step size calculation

- Initial  conservatively estimated
-  = 0.01  |D0|  /   |G|
- where G = FD2 + D2F
- Start with initial calculation of J, E(D,J) (Cost function)
- Take trial step of 
- D() = A() D() A()T
- Recalculate J and E(D, J)
- If the new E(D, J) is less than the previous E(D,J), use this time step

## Modified Gradient Descent Algorithm

- Compute D traditional linear way
- Compute eigenvalues and adjust
- Compute D based on new eigenvalues
- Calculate Ed = E(D,J) error
- Compute Initial 
- Take trial steps until convergence
- Find acceptable trial step  that gives lower Ed by halving the initial  up to 10 times
- Try step size of 2, 1⁄2. Compute corresponding Ed, D for each and pick the , D that gives the lowest Ed. Use  as the initial time step in the next convergence loop
- Test convergence (starting with second step)
- |Dnew – Dold|  /  |Dnew| < 10-4
- Optionally recompute convergence loop with new weight factors

## Weight Factor Computation	

- The initial weight factors were all set to 1
- Recompute weight factors to downweight data points (gradients) that don’t fit well (outliers)
- Compute residual at each gradient level from (0..n) as
- rq = J e –b(q) **.** D-I(q)
- Estimate Std. Dev. as
-  = [1/Nq  rq2] 1⁄2
- wq = [1 / sqrt( 1+ (rq /)2)] * Nq /
-           [1 / sqrt( 1+ (rq /)2)]

## FA - non-linear

- FA - linear

## Ratio of final E(D,J) to E(D,J) at =0 

**Ratio of final E(D,J) to E(D,J) at ********=0 **

## Number steps to convergence without reweighting

## Usage: 3dDWItoDT [options] gradient-file dataset
Computes 6 principle direction tensors from multiple gradient vectors
 and corresponding DTI image volumes.
 The program takes two parameters as input :  
    a 1D file of the gradient vectors with lines of ASCII floats Gxi,Gyi,Gzi.
    Only the non-zero gradient vectors are included in this file (no G0 line).
    a 3D bucket dataset with Np+1 sub-briks where the first sub-brik is the
    volume acquired with no diffusion weighting.
 Options:
   -automask =  mask dataset so that the tensors are computed only for
    high-intensity (presumably brain) voxels.  The intensity level is
    determined the same way that 3dClipLevel works.

   -nonlinear = compute iterative solution to avoid negative eigenvalues.
    This is the default method.

   -linear = compute simple linear solution

   -reweight = recompute weight factors at end of iterations and restart

   -max_iter n = maximum number of iterations for convergence (Default=10)
    values can range from -1 to any positive integer less than 101.
    A value of -1 is equivalent to the linear solution.
    A value of 0 results in only the initial estimate of the diffusion tensor
    solution adjusted to avoid negative eigenvalues.

- Computes 6 principle direction tensors from multiple gradient vectors
-  and corresponding DTI image volumes.
-  The program takes two parameters as input :
-     a 1D file of the gradient vectors with lines of ASCII floats Gxi,Gyi,Gzi.
-     Only the non-zero gradient vectors are included in this file (no G0 line).
-     a 3D bucket dataset with Np+1 sub-briks where the first sub-brik is the
-     volume acquired with no diffusion weighting.
-  Options:
-    -automask =  mask dataset so that the tensors are computed only for
-     high-intensity (presumably brain) voxels.  The intensity level is
-     determined the same way that 3dClipLevel works.
-    -nonlinear = compute iterative solution to avoid negative eigenvalues.
-     This is the default method.
-    -linear = compute simple linear solution
-    -reweight = recompute weight factors at end of iterations and restart
-    -max_iter n = maximum number of iterations for convergence (Default=10)
-     values can range from -1 to any positive integer less than 101.
-     A value of -1 is equivalent to the linear solution.
-     A value of 0 results in only the initial estimate of the diffusion tensor
-     solution adjusted to avoid negative eigenvalues.

- -max_iter_rw n = max number of iterations after reweighting (Default=5)
-     values can range from 1 to any positive integer less than 101.
-    -eigs = compute eigenvalues, eigenvectors and fractional anisotropy in
-     sub-briks 6-18. Computed as in 3dDTeig
-    -debug_briks = add sub-briks with Ed (error functional), Ed0 (original error)
-     and number of steps to convergence
-    -cumulative_wts = show overall weight factors for each gradient level
-     May be useful as a quality control
-    -verbose nnnnn = print convergence steps every nnnnn voxels that survive to
-     convergence loops (can be quite lengthy)
-    -drive_afni = show convergence graphs every nnnnn voxels that survive to convergence loops. AFNI must have NIML communications on (afni -niml).
-  Example:
-   3dDWItoDTnoise -prefix rw01 -automask -reweight -max_iter 10 \
-             -max_iter_rw 10 tensor25.1D grad02+orig.
-  The output is a 6 sub-brick bucket dataset containing Dxx,Dxy,Dxz,Dyy,Dyz,Dzz. Additional sub-briks may be
- appended with the -eigs and -debug_briks options.
-  These results are appropriate as the input to the 3dDTeig program.

## Slide 22

## Other Methods Comparison

- Tschumperlé and Deriche, Variational Framework
- Simultaneous spatial smoothing
- Complicated cost function versus x2
- (ln (I(0)/I(q)) - gkT D gk) +   (|D|)
- where (s) = log(1+s2) and (s) = sqrt(1+s2)
- We use (s) = s2, = 0 (no spatial contribution)
- Slightly more complicated gradient function
- G = (F+FT)D2 + D2 (F+FT)
- Linearized with ln(I(0) / I(q)) versus non-linear relationship,
- 	I(q) =  J e –b(q) **.** D + noise
- No reweighting
- Mangin, et al, Robust Tensor Estimation
- Cost function, Geman-McLure M-estimator, made to remove outliers, i2 / (i2+C2) where C=1.48 mediani |i| (we use reweighting)
- Similar to traditional method, does not enforce positive definiteness on D

## Future Directions

- Create and show fiber tracts in SUMA and AFNI
- Test model with computed DWI and artificial noise
- Add other indices (Lattice Index, Mean diffusivity,...)
- Refine model and algorithm
- Respond to AFNI user requests

## Acknowledgements

- Rick Reynolds
- Rich Hammett
- Ziad Saad
- Peter Basser
- Wolfgang Gaggl, Fernanda Tovar-Moll