Implements MTsat function

scripts/sct_compute_mtsat.py

@@ -0,0 +1,98 @@
+#!/usr/bin/env python
+# -*- coding: utf-8
+#########################################################################################
+#
+# Compute MT saturation map and T1 map from a PD-weigthed, a T1-weighted and MT-weighted FLASH images
+#
+# Reference paper:
+#    Helms G, Dathe H, Kallenberg K, Dechent P. High-resolution maps of magnetization transfer with inherent correction
+#    for RF inhomogeneity and T1 relaxation obtained from 3D FLASH MRI. Magn Reson Med 2008;60(6):1396-1407.
+
+# ---------------------------------------------------------------------------------------
+# Copyright (c) 2018 Polytechnique Montreal <www.neuro.polymtl.ca>
+# Author: Julien Cohen-Adad
+#
+# About the license: see the file LICENSE.TXT
+#########################################################################################
+
+
+import argparse
+
+
+def get_parser():
+    parser = argparse.ArgumentParser(
+        description='Compute MTsat and T1map.'
+                    'Reference: Helms G, Dathe H, Kallenberg K, Dechent P. High-resolution maps of magnetization '
+                    'transfer with inherent correction for RF inhomogeneity and T1 relaxation obtained from 3D FLASH '
+                    'MRI. Magn Reson Med 2008;60(6):1396-1407.')
+    parser.add_argument("-mt",
+                        help="Image with MT_ON",
+                        required=True)
+    parser.add_argument("-pd",
+                        help="Image PD weighted (typically, the MT_OFF)",
+                        required=True)
+    parser.add_argument("-t1",
+                        help="Image T1-weighted",
+                        required=True)
+    parser.add_argument("-trmt",
+                        help="TR [in ms] for mt image.",
+                        type=float,
+                        required=True)
+    parser.add_argument("-trpd",
+                        help="TR [in ms] for pd image.",
+                        type=float,
+                        required=True)
+    parser.add_argument("-trt1",
+                        help="TR [in ms] for t1 image.",
+                        type=float,
+                        required=True)
+    parser.add_argument("-famt",
+                        help="Flip angle [in deg] for mt image.",
+                        type=float,
+                        required=True)
+    parser.add_argument("-fapd",
+                        help="Flip angle [in deg] for pd image.",
+                        type=float,
+                        required=True)
+    parser.add_argument("-fat1",
+                        help="Flip angle [in deg] for t1 image.",
+                        type=float,
+                        required=True)
+    parser.add_argument("-b1map",
+                        help="B1 map",
+                        default=None)
+    parser.add_argument("-omtsat",
+                        help="Output file for MTsat. Default is mtsat.nii.gz",
+                        default=None)
+    parser.add_argument("-ot1map",
+                        help="Output file for T1map. Default is t1map.nii.gz",
+                        default=None)
+    parser.add_argument("-v",
+                        help="Verbose: 0 = no verbosity, 1 = verbose (default).",
+                        choices=('0', '1'),
+                        type=int,
+                        default=1)
+    return parser
+
+
+def run_main(args):
+    import sct_utils as sct
+    from spinalcordtoolbox.mtsat import mtsat
+
+    sct.start_stream_logger()
+
+    fname_mtsat, fname_t1map = mtsat.compute_mtsat_from_file(
+        args.mt, args.pd, args.t1, args.trmt, args.trpd, args.trt1, args.famt, args.fapd, args.fat1,
+        fname_b1map=args.b1map, fname_mtsat=args.omtsat, fname_t1map=args.ot1map, verbose=1)
+
+    sct.display_viewer_syntax([fname_mtsat, fname_t1map],
+                              colormaps=['gray', 'gray'],
+                              minmax=['-10,10', '0, 3'],
+                              opacities=['1', '1'],
+                              verbose=args.v)
+
+
+if __name__ == '__main__':
+    parser = get_parser()
+    arguments = parser.parse_args()
+    run_main(arguments)

scripts/sct_utils.py

@@ -319,8 +319,8 @@ def display_viewer_syntax(files, colormaps=[], minmax=[], opacities=[], mode='',
     sct.display_viewer_syntax([file1, file2], colormaps=['gray', 'red'], minmax=['', '0,1'], opacities=['', '0.7'])
     """
     list_viewer = ['fsleyes', 'fslview_deprecated', 'fslview']  # list of known viewers. Can add more.
-    dict_fslview = {'gray': 'Greyscale', 'red-yellow': 'Red-Yellow', 'blue-lightblue': 'Blue-Lightblue', 'red': 'Red', 'random': 'Random-Rainbow'}
-    dict_fsleyes = {'gray': 'greyscale', 'red-yellow': 'red-yellow', 'blue-lightblue': 'blue-lightblue', 'red': 'red', 'random': 'random'}
+    dict_fslview = {'gray': 'Greyscale', 'red-yellow': 'Red-Yellow', 'blue-lightblue': 'Blue-Lightblue', 'red': 'Red', 'random': 'Random-Rainbow', 'hsv': 'hsv'}
+    dict_fsleyes = {'gray': 'greyscale', 'red-yellow': 'red-yellow', 'blue-lightblue': 'blue-lightblue', 'red': 'red', 'random': 'random', 'hsv': 'hsv'}
     selected_viewer = None
 
     # find viewer

spinalcordtoolbox/mtsat/mtsat.py

@@ -0,0 +1,139 @@
+# coding: utf-8
+# This is the interface API to compute MTsat
+# Code is based on QMRLab: https://github.com/neuropoly/qMRLab
+# Author: Julien Cohen-Adad
+# Copyright (c) 2018 Polytechnique Montreal <www.neuro.polymtl.ca>
+# About the license: see the file LICENSE.TXT
+
+import os
+import sct_utils as sct
+from msct_image import Image
+import numpy as np
+
+
+def compute_mtsat(nii_mt, nii_pd, nii_t1,
+                  tr_mt, tr_pd, tr_t1,
+                  fa_mt, fa_pd, fa_t1,
+                  nii_b1map=None, verbose=1):
+    """
+    Compute MTsat and T1 map based on FLASH scans
+    :param nii_mt:
+    :param nii_pd:
+    :param nii_t1:
+    :param tr_mt:
+    :param tr_pd:
+    :param tr_t1:
+    :param fa_mt:
+    :param fa_pd:
+    :param fa_t1:
+    :param nii_b1map:
+    :param verbose:
+    :return:
+    """
+    # params
+    nii_t1map = None  # it would be possible in the future to input T1 map from elsewhere (e.g. MP2RAGE). Note: this
+                      # T1map needs to be in s unit.
+    b1correctionfactor = 0.4  # empirically defined in https://www.frontiersin.org/articles/10.3389/fnins.2013.00095/full#h3
+
+    # convert all TRs in s
+    tr_mt *= 0.001
+    tr_pd *= 0.001
+    tr_t1 *= 0.001
+
+    # Convert flip angles into radians
+    fa_mt_rad = np.radians(fa_mt)
+    fa_pd_rad = np.radians(fa_pd)
+    fa_t1_rad = np.radians(fa_t1)
+
+    # ignore warnings from division by zeros (will deal with that later)
+    seterr_old = np.seterr(over='ignore', divide='ignore', invalid='ignore')
+
+    # check if a T1 map was given in input; if not, compute R1
+    if nii_t1map is None:
+        # compute R1
+        sct.printv('Compute T1 map...', verbose)
+        r1map = 0.5 * np.true_divide((fa_t1_rad / tr_t1) * nii_t1.data - (fa_pd_rad / tr_pd) * nii_pd.data,
+                                     nii_pd.data / fa_pd_rad - nii_t1.data / fa_t1_rad)
+        # remove nans and clip unrelistic values
+        r1map = np.nan_to_num(r1map)
+        ind_unrealistic = np.where(r1map < 0.01)  # R1=0.01 s^-1 corresponds to T1=100s which is reasonable to clip
+        r1map[ind_unrealistic] = np.inf  # set to infinity so that these values will be 0 on the T1map
+        # compute T1
+        nii_t1map = nii_mt.copy()
+        nii_t1map.data = 1. / r1map
+    else:
+        sct.printv('Use input T1 map.', verbose)
+        r1map = 1. / nii_t1map.data
+
+    # Compute A
+    sct.printv('Compute A...', verbose)
+    a = (tr_pd * fa_t1_rad / fa_pd_rad - tr_t1 * fa_pd_rad / fa_t1_rad) * np.true_divide(np.multiply(nii_pd.data, nii_t1.data, dtype=float), tr_pd * fa_t1_rad * nii_t1.data - tr_t1 * fa_pd_rad * nii_pd.data)
+
+    # Compute MTsat
+    sct.printv('Compute MTsat...', verbose)
+    nii_mtsat = nii_mt.copy()
+    nii_mtsat.data = tr_mt * np.multiply((fa_mt_rad * np.true_divide(a, nii_mt.data) - 1), r1map, dtype=float) - (fa_mt_rad ** 2) / 2.
+    # sct.printv('nii_mtsat.data[95,89,14]' + str(nii_mtsat.data[95,89,14]), type='info')
+    # remove nans and clip unrelistic values
+    nii_mtsat.data = np.nan_to_num(nii_mtsat.data)
+    ind_unrealistic = np.where(np.abs(nii_mtsat.data) > 1)  # we expect MTsat to be on the order of 0.01
+    nii_mtsat.data[ind_unrealistic] = 0
+    # convert into percent unit (p.u.)
+    nii_mtsat.data *= 100
+
+    # Apply B1 correction to result
+    # Weiskopf, N., Suckling, J., Williams, G., Correia, M.M., Inkster, B., Tait, R., Ooi, C., Bullmore, E.T., Lutti, A., 2013. Quantitative multi-parameter mapping of R1, PD(*), MT, and R2(*) at 3T: a multi-center validation. Front. Neurosci. 7, 95.
+    if not nii_b1map is None:
+        nii_mtsat.data = np.true_divide(nii_mtsat.data * (1 - b1correctionfactor), (1 - b1correctionfactor * nii_b1map.data))
+
+    # set back old seterr settings
+    np.seterr(**seterr_old)
+
+    return nii_mtsat, nii_t1map
+
+
+def compute_mtsat_from_file(fname_mt, fname_pd, fname_t1, tr_mt, tr_pd, tr_t1, fa_mt, fa_pd, fa_t1, fname_b1map=None,
+                            fname_mtsat=None, fname_t1map=None, verbose=1):
+    """
+    Compute MTsat and T1map.
+    :param fname_mt:
+    :param fname_pd:
+    :param fname_t1:
+    :param tr_mt:
+    :param tr_pd:
+    :param tr_t1:
+    :param fa_mt:
+    :param fa_pd:
+    :param fa_t1:
+    :param fname_b1map:
+    :param fname_mtsat:
+    :param fname_t1map:
+    :param verbose:
+    :return: fname_mtsat: file name for MTsat map
+    :return: fname_t1map: file name for T1 map
+    """
+    # load data
+    sct.printv('Load data...', verbose)
+    nii_mt = Image(fname_mt)
+    nii_pd = Image(fname_pd)
+    nii_t1 = Image(fname_t1)
+    if fname_b1map is not None:
+        nii_b1map = Image(fname_b1map)
+
+    # compute MTsat
+    nii_mtsat, nii_t1map = compute_mtsat(nii_mt, nii_pd, nii_t1, tr_mt, tr_pd, tr_t1, fa_mt, fa_pd, fa_t1,
+                                         nii_b1map=nii_b1map, verbose=verbose)
+
+    # Output MTsat and T1 maps
+    # by default, output in the same directory as the input images
+    sct.printv('Generate output files...', verbose)
+    if fname_mtsat is None:
+        fname_mtsat = os.path.join(nii_mt.path, "mtsat.nii.gz")
+    nii_mtsat.setFileName(fname_mtsat)
+    nii_mtsat.save()
+    if fname_t1map is None:
+        fname_t1map = os.path.join(nii_mt.path, "t1map.nii.gz")
+    nii_t1map.setFileName(fname_t1map)
+    nii_t1map.save()
+
+    return fname_mtsat, fname_t1map