dti: Compute a tensor from diffusion dataĀ¶



% This m-file has been automatically generated using qMRgenBatch(dti)
    % Command Line Interface (CLI) is well-suited for automatization
    % purposes and Octave.
    % Please execute this m-file section by section to get familiar with batch
    % processing for dti on CLI.
    % Demo files are downloaded into dti_data folder.
    % Written by: Agah Karakuzu, 2017
    % =========================================================================


qMRinfo('dti'); % Describe the model
 dti: Compute a tensor from diffusion data

    Anisotropic Gaussian diffusion tensor
    Valid at relatively low b-value (i.e. ~ 2000 s/mm2)

    DiffusionData       4D DWI
    (SigmaNoise)        map of the standard deviation of the noise per voxel
    (Mask)              Binary mask to accelerate the fitting

    D                   [Dxx Dxy Dxz Dxy Dyy Dyz Dxz Dyz Dzz] Diffusion Tensor
    L1                  1rst eigenvalue of D
    L2                  2nd eigenvalue of D
    L3                  3rd eigenvalue of D
    FA                  Fractional Anisotropy: FA = sqrt(3/2)*sqrt(sum((L-L_mean).^2))/sqrt(sum(L.^2));
    S0_TEXX             Signal at b=0 at TE=XX
    (residue)           Fitting residuals

    At least 2 shells (e.g. b=1000 and b=0 s/mm2)
    diffusion gradient direction in 3D

    DiffusionData       Array [NbVol x 7]
    Gx                Diffusion Gradient x
    Gy                Diffusion Gradient y
    Gz                Diffusion Gradient z
    Gnorm (T/m)         Diffusion gradient magnitude
    Delta (s)         Diffusion separation
    delta (s)         Diffusion duration
    TE (s)            Echo time

    fitting type
    'linear'                              Solves the linear problem (ln(S/S0) = -bD)
    'non-linear (Rician Likelihood)'      Add an additional fitting step,
    using the Rician Likelihood.
    Rician noise bias                       only for non-linear fitting
    SigmaNoise map is prioritary.
    'Compute Sigma per voxel'             Sigma is estimated by computing the STD across repeated scans.
    'fix sigma'                           Use scd_noise_std_estimation to measure noise level. Use 'value' to fix Sigma.

    Example of command line usage (see a href="matlab: web(which('dti_batch.html'))" dti_batch.html/a):
    Model = dti
    %% LOAD DATA
    data.DiffusionData = load_nii_data('DiffusionData.nii.gz');
    data.SigmaNoise = load_nii_data('SigmaNoise.nii.gz');
    data.Mask = load_nii_data('Mask.nii.gz');
    % Specific voxel:         datavox = extractvoxel(data,voxel);
    % Interactive selection:  datavox = extractvoxel(data);
    voxel       = round(size(data.DiffusionData(:,:,:,1))/2); % pick FOV center
    datavox     = extractvoxel(data,voxel);
    FitResults  = Model.fit(datavox);
    Model.plotModel(FitResults, datavox); % plot fit results
    %% FIT all voxels
    FitResults = FitData(data,Model);
    % SAVE results to NIFTI
    FitResultsSave_nii(FitResults,'DiffusionData.nii.gz'); % use header from 'DiffusionData.nii.gz'

    For more examples: a href="matlab: qMRusage(dti);"qMRusage(dti)/a

    Author: Tanguy Duval, 2016

    Please cite the following if you use this module:
    Basser, P.J., Mattiello, J., LeBihan, D., 1994. MR diffusion tensor spectroscopy and imaging. Biophys. J. 66, 259?267.
    In addition to citing the package:
    Karakuzu A., Boudreau M., Duval T.,Boshkovski T., Leppert I.R., Cabana J.F.,
Gagnon I., Beliveau P., Pike G.B., Cohen-Adad J., Stikov N. (2020), qMRLab:
Quantitative MRI analysis, under one umbrella doi: 10.21105/joss.02343

    Reference page in Doc Center
    doc dti



a- create object

Model = dti;

b- modify options

         |- This section will pop-up the options GUI. Close window to continue.
    |- Octave is not GUI compatible. Modify Model.options directly.
Model = Custom_OptionsGUI(Model); % You need to close GUI to move on.


a- load experimental data

         |- dti object needs 3 data input(s) to be assigned:
    |-   DiffusionData
    |-   SigmaNoise
    |-   Mask
data = struct();
    % DiffusionData.nii.gz contains [74   87   50  109] data.
    % Mask.nii.gz contains [74  87  50] data.

b- fit dataset

           |- This section will fit data.
FitResults = FitData(data,Model,0);
=============== qMRLab::Fit ======================
    Operation has been started: dti
    Elapsed time is 1.655054 seconds.
    Operation has been completed: dti

c- show fitting results

         |- Output map will be displayed.
    |- If available, a graph will be displayed to show fitting in a voxel.
    |- To make documentation generation and our CI tests faster for this model,
    we used a subportion of the data (40X40X40) in our testing environment.
    |- Therefore, this example will use FitResults that comes with OSF data for display purposes.
    |- Users will get the whole dataset (384X336X224) and the script that uses it for demo
    via qMRgenBatch(qsm_sb) command.
FitResults_old = load('FitResults/FitResults.mat');

d- Save results

         |-  qMR maps are saved in NIFTI and in a structure FitResults.mat
    that can be loaded in qMRLab graphical user interface
    |-  Model object stores all the options and protocol.
    It can be easily shared with collaborators to fit their
    own data or can be used for simulation.
FitResultsSave_nii(FitResults, 'dti_data/DiffusionData.nii.gz');
Warning: Directory already exists.


   |- This section can be executed to run simulations for dti.

a- Single Voxel Curve

         |- Simulates Single Voxel curves:
    (1) use equation to generate synthetic MRI data
    (2) add rician noise
    (3) fit and plot curve
      x = struct;
    x.L1 = 2;
    x.L2 = 0.7;
    x.L3 = 0.7;
    Opt.SNR = 50;
    % run simulation
    figure('Name','Single Voxel Curve Simulation');
    FitResult = Model.Sim_Single_Voxel_Curve(x,Opt);

b- Sensitivity Analysis

         |-    Simulates sensitivity to fitted parameters:
    (1) vary fitting parameters from lower (lb) to upper (ub) bound.
    (2) run Sim_Single_Voxel_Curve Nofruns times
    (3) Compute mean and std across runs
      %              L1            L2            L3
    OptTable.st = [2             0.7           0.7]; % nominal values
    OptTable.fx = [0             1             1]; %vary L1...
    OptTable.lb = [0             0             0]; %...from 0
    OptTable.ub = [5             5             5]; %...to 5
    Opt.SNR = 50;
    Opt.Nofrun = 5;
    % run simulation
    SimResults = Model.Sim_Sensitivity_Analysis(OptTable,Opt);
    figure('Name','Sensitivity Analysis');
    SimVaryPlot(SimResults, 'L1' ,'L1' );