dMRI-Lab 1.0
demo_MiSFIT_developer.m File Reference

This demo illustrates most of the functions and utilities in the micro/ folder, i.e, the module related to the computational analysis of multi-shell dMRI data. More...

Detailed Description

This demo illustrates most of the functions and utilities in the micro/ folder, i.e, the module related to the computational analysis of multi-shell dMRI data.

It is based on a linear convolutional model of the so-called ODF with an elemental impulse response corresponding to a prolate tensor with eigenvalues lpar and lperp.

Estimate the diffusion tensor

This section illustrates the usage of the atti2micro, which is the core of micro-structure computation. It returns f (the volume fraction of confined, i.e. non-free water), lpar, and lperp ( the two eigenvalues of the elemental tensor in the linear convolution model). It is always assured that mlp<=lperp<=lpar<=ADC0 and that 0<=f<=1. In the simplest use case, the function reduces to the original Kaden's method, so that f=1 (except for the background) and no regularization is imposed to lpar and lperp:

atti2micro dwifile.nii gifile.bvec bifile.bval lambdaparfile.nii lambdaperpfile.nii ffile.nii -mask maskfile.nii -usef 0 -bth 100 -mlperp 0.0 -mu 0.0 -verbose 1
function atti2micro(in dwifile, in gifile, in bifile, in lambdaparfile, in lambdaperpfile, in ffile, in varargin)
Note
This is the simplest use case
atti2micro dwifile.nii gifile.bvec bifile.bval lambdaparfile.nii lambdaperpfile.nii ffile.nii -mask maskfile.nii -usef 1 -bth 100 -mlperp 0.0 -mu 0.0 -verbose 1
Note
We have changed the 'usef' opt. argument

ODFs estimation

Once the micro-structure model is available, computing the ODFs field is a spherical deconvolution problem. Since we are using Spherical Harmonics, this is a trivial task that reduces to a linear LS fitting of SH coefficients. The present section illustrates how to use the function micro2shodf. This may take a good while since we invert one matrix per voxel:

micro2shodf dwifile.nii gifile.bvec bifile.bval lambdaparfile.nii lambdaperpfile.nii ffile.nii shfile.nii -mask maskfile.nii -lambda 0.0005 -L 8 -optimal 1
function micro2shodf(in dwifile, in gifile, in bifile, in lambdaparfile, in lambdaperpfile, in ffile, in shfile, in varargin)

The color-by-orientation, FA, and ODFs field look quite reasonable even without any denoising (this is not always the case with HARDI). To avoid the large computation time associated to the inversion of the LS matrix at each voxel, we may think of a sub-optimal scheme where only one matrix is pre-computed and inverted for the whole data set:

micro2shodf dwifile.nii gifile.bvec bifile.bval lambdaparfile.nii lambdaperpfile.nii ffile.nii shfile.nii -mask maskfile.nii -lambda 0.0005 -L 8 -optimal 0
Note
'optimal' is now false and we use 'reweight'

The results look pretty similar to those obtained with the optimal estimation method, both for the color-coded/FA maps and the ODFs field.

Characterization of diffusion via moments computation. Full moments

The next three sections illustrate the computation of arbitrary order moments to characterize the diffusion process, i.e., they illustrate the usage of micro2moments. We start with the so-called 'full moments', computed as integrals in the whole 3-D space of either the attenaution signal or the EAP. This kind of moments depend just on the micro-structure (lpar and lperp) and not on the ODF. For the attenaution signal, E(q):

micro2moments shfile.nii lambdaparfile.nii lambdaperpfile.nii ffile.nii mufile.nii -mask maskfile.nii -nu -2
function micro2moments(in shfile, in lambdaparfile, in lambdaperpfile, in ffile, in mufile, in varargin)
Note
Different moments provide different contrast features

RTOP

micro2moments shfile.nii lambdaparfile.nii lambdaperpfile.nii ffile.nii mufile.nii -mask maskfile.nii -nu 0 -clean 100
Note
This is just RTOP

QMSD

micro2moments shfile.nii lambdaparfile.nii lambdaperpfile.nii ffile.nii mufile.nii -mask maskfile.nii -nu 2 -clean 100
Note
This is just QMSD

MSD

micro2moments shfile.nii lambdaparfile.nii lambdaperpfile.nii ffile.nii mufile.nii -mask maskfile.nii -nu 2
Note
This is just MSD

Axial moments

Axial moments depend on a direction provided at each voxel, typically the direction of maximum diffusion. But they could be otherwise computed over vector fields computed from tractography, or over any other meaningful vector field:

RTPP

micro2moments shfile.nii lambdaparfile.nii lambdaperpfile.nii ffile.nii mufile.nii -mask maskfile.nii -nu 0 -clean 100
Note
This is just RTPP

Planar moments

Planar moments depend on a direction provided at each voxel, typically the direction of maximum diffusion. But they could be otherwise computed over vector fields computed from tractography, or over any other meaningful vector field. Planar moments will depend on both the micro-structure model and the ODF. Begin with E(q):

RTAP

micro2moments shfile.nii lambdaparfile.nii lambdaperpfile.nii ffile.nii mufile.nii -mask maskfile.nii -nu 0 -clean 100
Note
This is just RTAP
Note
Finally note that for nu=0 there is a direct relationship between planar moments computed over E(q) and axial moments compute over P(R) (and vice-versa) by virtue of the central section theorem. More concretely:
  • The axial moment of order nu=0 computed over E(q) will equal the planar moment of order nu=0 computed over P(R) (this is just RTPP).
  • The planar moment of order nu=0 computed over E(q) will equal the axial moment of order nu=0 computed over P(R) (this is just RTAP).

Propagator Anisotropy

The PA is related to the distance from the EAP (or, equivalently, E(q)) to its isotropic equivalent defined as its spherical average. In our framework, this reduces once again to a manipulation of the SH coefficients. In this case there is not an analytical closed-form available, so the calculations rely on precomputed integrals implemented in dmri_compute_PA_weights. However, it is unlikely that you need to explicitly call that function, since it will be invoked as needed by micro2pa, which is the function we illustrate in this section. In the simplest case,

micro2pa shfile.nii lambdaparfile.nii lambdaperpfile.nii ffile.nii pafile.nii -mask maskfile.nii -micro 1
function micro2pa(in shfile, in lambdaparfile, in lambdaperpfile, in pafile, in varargin)

However, the raw values of the PA are not well distributed over the range [0,1], so that a gamma correction is usually applied to improve its visual performance. You can fix this by properly tunning the optional input parameter epsilon:

micro2pa shfile.nii lambdaparfile.nii lambdaperpfile.nii ffile.nii pafile.nii -mask maskfile.nii -epsilon 0.4 -micro 1

Non-gaussianity

The NG is related to the distance from the EAP (or, equivalently, E(q)) to its Gaussian equivalent, which may be estimated in several ways. In our framework, this reduces once again to a manipulation of the SH coefficients. In this case there is not an analytical closed-form available, so the calculations rely on precomputed integrals implemented in dmri_compute_PA_weights. However, it is unlikely that you need to explicitly call that function, since it will be invoked as needed by micro2ng, which is the function we illustrate in this section. In the simplest case,

micro2atti shfile.nii lambdaparfile.nii lambdaperpfile.nii ffile.nii gifile.bvec bifile.bval attifile.nii -mask maskfile.nii
function micro2atti(in shfile, in lparfile, in lperpfile, in ffile, in gifile, in bifile, in attifile, in varargin)