inversion_recovery: Compute a T1 map using Inversion Recovery dataĀ¶

https://mybinder.org/badge_logo.svg

Contents

% This m-file has been automatically generated using qMRgenBatch(inversion_recovery)
    % 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 inversion_recovery on CLI.
    %
    % Demo files are downloaded into inversion_recovery_data folder.
    %
    % Written by: Agah Karakuzu, 2017
    % =========================================================================
    

I- DESCRIPTION

qMRinfo('inversion_recovery'); % Describe the model
    
 inversion_recovery: Compute a T1 map using Inversion Recovery data

    Assumptions:
    (1) Gold standard for T1 mapping
    (2) Infinite TR

    Inputs:
    IRData      Inversion Recovery data (4D)
    (Mask)      Binary mask to accelerate the fitting (OPTIONAL)

    Outputs:
    T1          transverse relaxation time [ms]
    b           arbitrary fit parameter (S=a + b*exp(-TI/T1))
    a           arbitrary fit parameter (S=a + b*exp(-TI/T1))
    idx         index of last polarity restored datapoint (only used for magnitude data)
    res         Fitting residual


    Protocol:
    IRData  [TI1 TI2...TIn] inversion times [ms]

    Options:
    Method          Method to use in order to fit the data, based on whether complex or only magnitude data acquired.
    'complex'         RD-NLS (Reduced-Dimension Non-Linear Least Squares)
    S=a + b*exp(-TI/T1)
    'magnitude'      RD-NLS-PR (Reduced-Dimension Non-Linear Least Squares with Polarity Restoration)
    S=|a + b*exp(-TI/T1)|

    Example of command line usage (see also a href="matlab: showdemo inversion_recovery_batch"showdemo inversion_recovery_batch/a):
    Model = inversion_recovery;  % Create class from model
    Model.Prot.IRData.Mat=[350.0000; 500.0000; 650.0000; 800.0000; 950.0000; 1100.0000; 1250.0000; 1400.0000; 1700.0000];
    data = struct;  % Create data structure
    data.MET2data ='IRData.mat';  % Load data
    data.Mask = 'Mask.mat';
    FitResults = FitData(data,Model); %fit data
    FitResultsSave_mat(FitResults);

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

    Author: Ilana Leppert, 2017

    References:
    Please cite the following if you use this module:
    A robust methodology for in vivo T1 mapping. Barral JK, Gudmundson E, Stikov N, Etezadi-Amoli M, Stoica P, Nishimura DG. Magn Reson Med. 2010 Oct;64(4):1057-67. doi: 10.1002/mrm.22497.
    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 inversion_recovery


    

II- MODEL PARAMETERS

a- create object

Model = inversion_recovery;
    

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.
    

III- FIT EXPERIMENTAL DATASET

a- load experimental data

         |- inversion_recovery object needs 2 data input(s) to be assigned:
    |-   IRData
    |-   Mask
data = struct();

    % IRData.mat contains [128  128    1    9] data.
    load('inversion_recovery_data/IRData.mat');
    % Mask.mat contains [128  128] data.
    load('inversion_recovery_data/Mask.mat');
    data.IRData= double(IRData);
    data.Mask= double(Mask);
    

b- fit dataset

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

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');
    qMRshowOutput(FitResults_old,data,Model);
    

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);
    Model.saveObj('inversion_recovery_Demo.qmrlab.mat');
    
Warning: Directory already exists.
    

V- SIMULATIONS

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

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.T1 = 600;
    x.rb = -1000;
    x.ra = 500;
    % Set simulation options
    Opt.SNR = 50;
    Opt.T1 = 600;
    Opt.M0 = 1000;
    Opt.TR = 3000;
    Opt.FAinv = 180;
    Opt.FAexcite = 90;
    Opt.Updateinputvariables = false;
    % 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
      %              T1            rb            ra
    OptTable.st = [6e+02         -1e+03        5e+02]; % nominal values
    OptTable.fx = [0             1             1]; %vary T1...
    OptTable.lb = [0.0001        -1e+04        0.0001]; %...from 0.0001
    OptTable.ub = [5e+03         0             1e+04]; %...to 5000
    Opt.SNR = 50;
    Opt.Nofrun = 5;
    % run simulation
    SimResults = Model.Sim_Sensitivity_Analysis(OptTable,Opt);
    figure('Name','Sensitivity Analysis');
    SimVaryPlot(SimResults, 'T1' ,'T1' );