dMRI-Lab 1.0
atti2hydidsi.m File Reference

Functions

function atti2hydidsi (in dwifile, in gifile, in bifile, in eapfile, in dtifile, in latticefile, in Qxfile, in Qyfile, in Qzfile, in resfile, in laplfile, in loptfile, in varargin)
 

Function Documentation

◆ atti2hydidsi()

function atti2hydidsi ( in dwifile,
in gifile,
in bifile,
in eapfile,
in dtifile,
in latticefile,
in Qxfile,
in Qyfile,
in Qzfile,
in resfile,
in laplfile,
in loptfile,
in varargin )

Computes the Ensemble Average Propagator (eap) at a regular Cartesian lattice from a set of scattered q-space measurements (atti) using Hybrid Diffusion Imaging without q-space re-gridding. To that end, the R-space is first aligned with the eigenvectors of a tensor model (the z axis aligns with the maximum diffusion direction) and an appropriate lattice is designed so that the whole extent of the EAP is covered according to the computed eigenvalues. From this estimate, the q-space is cropped to avoid aliasing, meaning that some q-space samples might be discarded for some voxels. Accordingly, the tensor model estimated is necessary to interpret the computed eap (and it is returned by this function). The eap is estimated by solving a quadratic programming problem that involves the encoding matrix that relates q-space measurments to R-space Fourier coefficients (i.e. EAP samples).

Parameters
dwifileA string specifying the file containing the diffusion-weighted images.
gifileA string specifying the file containing the gradient directions.
bifileA string specifying the file containing the b-values.
eapfileA string specifying the file to save the unrolled Cartesian samples of the eap. If the lattice has radii Nx x Ny x Nz, a grand total of Nl=(2*Nx+1)(2*Ny+1)*(2*Nz+1) samples would be needed. However, due to the antipodal symmetry of the eap we need just K=(Nl+1)/2. The samples are arranged as follows:
  • [x,y,z] = meshgrid(-Nx:Nx,-Ny:Ny,-Nz:Nz);
  • x = x(:); y = y(:); z = z(:);
  • Nl = size(x,1);
  • x = x(floor(Nl/2):end);
  • y = y(floor(Nl/2):end);
  • z = z(floor(Nl/2):end);
dtifileA string specifying the file to save the estimated tensor model at each voxel.
Note
this is necessary to interpret the data in eap, since the lattice is aligned at each voxel to the eigenvectors of dti (and scaled according to its eigenvalues).
Parameters
latticefileA string specifying the file to save the radii of the lattice at each voxel. This is useful in case this parameter is automatically determined instead of fed as an input argument.
Qxfile,Qyfile,QzfileStrings specifying the files to save estimated bandwidths of the q-space signal along each dimension. They are necessary to interpret the computed eap.
resfileA string specifying the file to save the residual of the QP problem at each voxel.
laplfileA string specifying the file to save the energy of the Laplacian at each voxel.
loptfileA string specifying the file to save the optimal value of lambda computed at each voxel (makes sense only if Generalized Cross Validation is used to optimally set the value of lambda at each voxel, see below).
vararginOptional arguments passed as -name value pairs (e.g. -wlsit 5).
Returns
- eapfile: The unrolled Cartesian samples of the eap.
- dtifile: The estimated tensor model at each voxel.
- latticefile: The radii of the lattice at each voxel.
- Qxfile, Qyfile, Qzfile: The estimated bandwidths of the q-space signal along each dimension.
- resfile: The residual of the QP problem at each voxel.
- laplfile: The energy of the Laplacian at each voxel.
- loptfile: The optimal value of lambda computed at each voxel.

Mandatory Inputs:

  • dwifile: The diffusion-weighted images.
  • gifile: The gradient directions.
  • bifile: the b-values.
  • eapfile: A string specifying the file to save the unrolled Cartesian samples of the eap.
  • dtifile: A string specifying the file to save the estimated tensor model at each voxel.
  • latticefile: A string specifying the file to save the radii of the lattice at each voxel.
  • Qxfile, Qyfile, Qzfile: Strings specifying the files to save estimated bandwidths of the q-space signal along each dimension.
  • resfile: A string specifying the file to save the residual of the QP problem at each voxel.
  • laplfile: A string specifying the file to save the energy of the Laplacian at each voxel.
  • loptfile: A string specifying the file to save the optimal value of lambda computed at each voxel.
  • varargin: Optional arguments passed as -name value pairs.

Optional Arguments General Optional Parameters:

  • lambda: this is a regularization parameter that penalizes highly irregular solutions of the EAP. In precise terms, the quadratic deviation from the linear model is added lambda times the energy of the Laplacian of the EAP to compute the final cost function.
    Note
    you can pass Inf, NaN or any negative value of lambda to the fucntion. In this case, its value will be adaptively computed at each voxel using Generalized Cross Validation (GCV), BUT: this will blow up the overall computation time, since GCV implies repeated matrix inversions (default 1.0e-3).
  • tl, tu: the lower and upper thresholds, respectively, defining the range the dwi will lay within, so that tl should be close to 0 and tu should be close to 1 (default: 1.0e-7, 1-1.0e-7).
  • bcut: b-values over this value are ignored to estimate the diffusion tensor, which is used to transform the native space to the anatomical space (default: 2000 s/mm^2).
  • ADC0: estimated diffusivity of free water at body temperature. It is used to determine the lower and upper bounds of the eigenvalues of the dti and perform sanity checks (default: 3.0e-3).
  • tau: the effective diffusion time of the acquisition, necessary to determine the actual support of the q-space (default: 35e-3).
  • Rth: the threshold used to estimate the support of the EAP. The last lattice node along each ‘x’-‘y’-‘z’ axis will be placed in the position where the DTI-EAP falls to Rth times its value at the origin, which is trivially computed from the eigenvalues l1, l2, and l3 (default: 0.01).
  • lattice: the radii of the Cartesian grid where the eap will be sampled, whose size will be typically in the range 3 x 3 x 3. You may pass an empty array [] to let the function determine the proper size. In any case, the third component should be greate or equal than the first and the second (default: []).
  • const: wether (true) or not (false) constrain the estimated EAP to be positive and have unit mass. With these constraints, the estimation has a rigorous physical meaning but a quadratic problem must be solved at each voxel. Without it, we just need to solve an unconstrained least squares problem, which will be pretty faster (default: true, i.e. use constraints). Options when const is true: These options are directly passed to dmriQuadprog (check the documentation therein):
    • miters (default: 100)
    • otol (default: 1.0e-6)
    • stol (default: 1.0e-6)
    • ctol (default: 1.0e-8)

Other General Options:

  • usemex: whether (true) or not (false) use the C/MEX implementation of the method. This implementation can be notably faster (or not) depending on the CPU used. It makes extensive use of BLAS routines, so the performance is computer dependent. In POSIX systems the implementation is multi-threaded, so that you don't need a working license for the parallel computing toolbox to run the code in parallel (default: false). NOTE: if 'usemex' is set true, the 'const' option is ignored and the constrained problem is always solved.
  • maxthreads: if 'usemex' is set true, the algorithm is run as a multi-threaded mex. This is the maximum allowed number of threads, which can indeed be reduced if it exceeds the number of logical cores (default: the number of logical cores in the machine).
  • mask: a MxNxP array of logicals. Only those voxels where mask is true are processed, the others are filled with zeros (default: all trues).
  • verbose: whether (true) or not (false) show a progress bar to check the progress of the algorithm (default: false).

Example:

atti2hydidsi dwifile.nii gifile.bvec bifile.bval eapfile.nii dtifile.nii latticefile.nii Qxfile.nii Qyfile.nii Qzfile.nii resfile.nii laplfile.nii loptfile.nii -mask mask_file.nii
function atti2hydidsi(in dwifile, in gifile, in bifile, in eapfile, in dtifile, in latticefile, in Qxfile, in Qyfile, in Qzfile, in resfile, in laplfile, in loptfile, in varargin)
Note
This is the simplest use case
See also
atti2hydidsi