Merge branch 'master' into normalizeBeforeMAPCA

.github/ISSUE_TEMPLATE/config.yml

@@ -0,0 +1,4 @@
+contact_links:
+    - name: Usage question
+      url: https://neurostars.org/tag/tedana
+      about: Please ask questions about using tedana here.

tedana/combine.py

@@ -14,7 +14,7 @@
 @due.dcite(Doi('10.1002/(SICI)1522-2594(199907)42:1<87::AID-MRM13>3.0.CO;2-O'),
            description='T2* method of combining data across echoes using '
                        'monoexponential equation.')
-def _combine_t2s(data, tes, ft2s):
+def _combine_t2s(data, tes, ft2s, report=True):
     """
     Combine data across echoes using weighted averaging according to voxel-
     (and sometimes volume-) wise estimates of T2*.
@@ -27,6 +27,8 @@ def _combine_t2s(data, tes, ft2s):
         Echo times in milliseconds.
     ft2s : (M [x T] X 1) array_like
         Either voxel-wise or voxel- and volume-wise estimates of T2*.
+    report : bool, optional
+        Whether to log a description of this step or not. Default is True.
 
     Returns
     -------
@@ -42,14 +44,15 @@ def _combine_t2s(data, tes, ft2s):
       Medicine: An Official Journal of the International Society
       for Magnetic Resonance in Medicine, 42(1), 87-97.
     """
-    RepLGR.info("Multi-echo data were then optimally combined using the "
-                "T2* combination method (Posse et al., 1999).")
-    RefLGR.info("Posse, S., Wiese, S., Gembris, D., Mathiak, K., Kessler, "
-                "C., Grosse‐Ruyken, M. L., ... & Kiselev, V. G. (1999). "
-                "Enhancement of BOLD‐contrast sensitivity by single‐shot "
-                "multi‐echo functional MR imaging. Magnetic Resonance in "
-                "Medicine: An Official Journal of the International Society "
-                "for Magnetic Resonance in Medicine, 42(1), 87-97.")
+    if report:
+        RepLGR.info("Multi-echo data were then optimally combined using the "
+                    "T2* combination method (Posse et al., 1999).")
+        RefLGR.info("Posse, S., Wiese, S., Gembris, D., Mathiak, K., Kessler, "
+                    "C., Grosse‐Ruyken, M. L., ... & Kiselev, V. G. (1999). "
+                    "Enhancement of BOLD‐contrast sensitivity by single‐shot "
+                    "multi‐echo functional MR imaging. Magnetic Resonance in "
+                    "Medicine: An Official Journal of the International Society "
+                    "for Magnetic Resonance in Medicine, 42(1), 87-97.")
     n_vols = data.shape[-1]
     alpha = tes * np.exp(-tes / ft2s)
     if alpha.ndim == 2:
@@ -71,7 +74,7 @@ def _combine_t2s(data, tes, ft2s):
 @due.dcite(Doi('10.1002/mrm.20900'),
            description='PAID method of combining data across echoes using just '
                        'SNR/signal and TE.')
-def _combine_paid(data, tes):
+def _combine_paid(data, tes, report=True):
     """
     Combine data across echoes using SNR/signal and TE via the
     parallel-acquired inhomogeneity desensitized (PAID) ME-fMRI combination
@@ -83,6 +86,8 @@ def _combine_paid(data, tes):
         Masked data.
     tes : (1 x E) array_like
         Echo times in milliseconds.
+    report : bool, optional
+        Whether to log a description of this step or not. Default is True.
 
     Returns
     -------
@@ -99,16 +104,17 @@ def _combine_paid(data, tes):
       International Society for Magnetic Resonance in Medicine,
       55(6), 1227-1235.
     """
-    RepLGR.info("Multi-echo data were then optimally combined using the "
-                "parallel-acquired inhomogeneity desensitized (PAID) "
-                "combination method.")
-    RefLGR.info("Poser, B. A., Versluis, M. J., Hoogduin, J. M., & Norris, "
-                "D. G. (2006). BOLD contrast sensitivity enhancement and "
-                "artifact reduction with multiecho EPI: parallel‐acquired "
-                "inhomogeneity‐desensitized fMRI. "
-                "Magnetic Resonance in Medicine: An Official Journal of the "
-                "International Society for Magnetic Resonance in Medicine, "
-                "55(6), 1227-1235.")
+    if report:
+        RepLGR.info("Multi-echo data were then optimally combined using the "
+                    "parallel-acquired inhomogeneity desensitized (PAID) "
+                    "combination method.")
+        RefLGR.info("Poser, B. A., Versluis, M. J., Hoogduin, J. M., & Norris, "
+                    "D. G. (2006). BOLD contrast sensitivity enhancement and "
+                    "artifact reduction with multiecho EPI: parallel‐acquired "
+                    "inhomogeneity‐desensitized fMRI. "
+                    "Magnetic Resonance in Medicine: An Official Journal of the "
+                    "International Society for Magnetic Resonance in Medicine, "
+                    "55(6), 1227-1235.")
     n_vols = data.shape[-1]
     snr = data.mean(axis=-1) / data.std(axis=-1)
     alpha = snr * tes
@@ -215,17 +221,19 @@ def make_optcom(data, tes, adaptive_mask, t2s=None, combmode='t2s', verbose=True
     data = data[mask, :, :]  # mask out unstable voxels/samples
     tes = np.array(tes)[np.newaxis, ...]  # (1 x E) array_like
     combined = np.zeros((data.shape[0], data.shape[2]))
+    report = True
     for echo in np.unique(adaptive_mask[mask]):
         echo_idx = adaptive_mask[mask] == echo
 
         if combmode == 'paid':
             combined[echo_idx, :] = _combine_paid(data[echo_idx, :echo, :],
-                                                  tes[:echo])
+                                                  tes[:echo], report=report)
         else:
             t2s_ = t2s[mask, ..., np.newaxis]  # mask out empty voxels/samples
 
             combined[echo_idx, :] = _combine_t2s(
-                data[echo_idx, :echo, :], tes[:, :echo], t2s_[echo_idx, ...])
+                data[echo_idx, :echo, :], tes[:, :echo], t2s_[echo_idx, ...], report=report)
+        report = False
 
     combined = unmask(combined, mask)
     return combined

tedana/decay.py

@@ -61,31 +61,39 @@ def monoexponential(tes, s0, t2star):
     return s0 * np.exp(-tes / t2star)
 
 
-def fit_monoexponential(data_cat, echo_times, adaptive_mask):
+def fit_monoexponential(data_cat, echo_times, adaptive_mask, report=True):
     """
     Fit monoexponential decay model with nonlinear curve-fitting.
 
     Parameters
     ----------
     data_cat : (S x E x T) :obj:`numpy.ndarray`
         Multi-echo data.
-    echo_times
+    echo_times : (E,) array_like
         Echo times in milliseconds.
-    adaptive_mask
+    adaptive_mask : (S,) :obj:`numpy.ndarray`
         Array where each value indicates the number of echoes with good signal
         for that voxel.
+    report : bool, optional
+        Whether to log a description of this step or not. Default is True.
 
     Returns
     -------
-    t2s_limited, s0_limited, t2s_full, s0_full
+    t2s_limited, s0_limited, t2s_full, s0_full : (S,) :obj:`numpy.ndarray`
+        T2* and S0 estimate maps.
+
+    Notes
+    -----
+    This method is slower, but more accurate, than the log-linear approach.
     """
-    RepLGR.info("A monoexponential model was fit to the data at each voxel "
-                "using nonlinear model fitting in order to estimate T2* and S0 "
-                "maps, using T2*/S0 estimates from a log-linear fit as "
-                "initial values. For each voxel, the value from the adaptive "
-                "mask was used to determine which echoes would be used to "
-                "estimate T2* and S0. In cases of model fit failure, T2*/S0 "
-                "estimates from the log-linear fit were retained instead.")
+    if report:
+        RepLGR.info("A monoexponential model was fit to the data at each voxel "
+                    "using nonlinear model fitting in order to estimate T2* and S0 "
+                    "maps, using T2*/S0 estimates from a log-linear fit as "
+                    "initial values. For each voxel, the value from the adaptive "
+                    "mask was used to determine which echoes would be used to "
+                    "estimate T2* and S0. In cases of model fit failure, T2*/S0 "
+                    "estimates from the log-linear fit were retained instead.")
     n_samp, n_echos, n_vols = data_cat.shape
 
     # Currently unused
@@ -156,7 +164,42 @@ def fit_monoexponential(data_cat, echo_times, adaptive_mask):
 
 
 def fit_loglinear(data_cat, echo_times, adaptive_mask, report=True):
-    """
+    """Fit monoexponential decay model with log-linear regression.
+
+    The monoexponential decay function is fitted to all values for a given
+    voxel across TRs, per TE, to estimate voxel-wise :math:`S_0` and :math:`T_2^*`.
+    At a given voxel, only those echoes with "good signal", as indicated by the
+    value of the voxel in the adaptive mask, are used.
+    Therefore, for a voxel with an adaptive mask value of five, the first five
+    echoes would be used to estimate T2* and S0.
+
+    Parameters
+    ----------
+    data_cat : (S x E x T) :obj:`numpy.ndarray`
+        Multi-echo data. S is samples, E is echoes, and T is timepoints.
+    echo_times : (E,) array_like
+        Echo times in milliseconds.
+    adaptive_mask : (S,) :obj:`numpy.ndarray`
+        Array where each value indicates the number of echoes with good signal
+        for that voxel.
+    report : :obj:`bool`, optional
+        Whether to log a description of this step or not. Default is True.
+
+    Returns
+    -------
+    t2s_limited, s0_limited, t2s_full, s0_full: (S,) :obj:`numpy.ndarray`
+        T2* and S0 estimate maps.
+
+    Notes
+    -----
+    The approach used in this function involves transforming the raw signal values
+    (:math:`log(|data| + 1)`) and then fitting a line to the transformed data using
+    ordinary least squares.
+    This results in two parameter estimates: one for the slope  and one for the intercept.
+    The slope estimate is inverted (i.e., 1 / slope) to get  :math:`T_2^*`,
+    while the intercept estimate is exponentiated (i.e., e^intercept) to get :math:`S_0`.
+
+    This method is faster, but less accurate, than the nonlinear approach.
     """
     if report:
         RepLGR.info("A monoexponential model was fit to the data at each voxel "
@@ -215,7 +258,7 @@ def fit_loglinear(data_cat, echo_times, adaptive_mask, report=True):
     return t2s_limited, s0_limited, t2s_full, s0_full
 
 
-def fit_decay(data, tes, mask, adaptive_mask, fittype):
+def fit_decay(data, tes, mask, adaptive_mask, fittype, report=True):
     """
     Fit voxel-wise monoexponential decay models to `data`
 
@@ -234,6 +277,8 @@ def fit_decay(data, tes, mask, adaptive_mask, fittype):
         given sample
     fittype : {loglin, curvefit}
         The type of model fit to use
+    report : bool, optional
+        Whether to log a description of this step or not. Default is True.
 
     Returns
     -------
@@ -254,18 +299,11 @@ def fit_decay(data, tes, mask, adaptive_mask, fittype):
 
     Notes
     -----
-    1.  Fit monoexponential decay function to all values for a given voxel
-        across TRs, per TE, to estimate voxel-wise :math:`S_0` and
-        :math:`T_2^*`:
-
-        .. math::
-            S(TE) = S_0 * exp(-R_2^* * TE)
-
-            T_2^* = 1 / R_2^*
-
-    2.  Replace infinite values in :math:`T_2^*` map with 500 and NaN values
-        in :math:`S_0` map with 0.
-    3.  Generate limited :math:`T_2^*` and :math:`S_0` maps by doing something.
+    This function replaces infinite values in the :math:`T_2^*` map with 500 and
+    :math:`T_2^*` values less than or equal to zero with 1.
+    Additionally, very small :math:`T_2^*` values above zero are replaced with a floor
+    value to prevent zero-division errors later on in the workflow.
+    It also replaces NaN values in the :math:`S_0` map with 0.
     """
     if data.shape[1] != len(tes):
         raise ValueError('Second dimension of data ({0}) does not match number '
@@ -285,10 +323,10 @@ def fit_decay(data, tes, mask, adaptive_mask, fittype):
 
     if fittype == 'loglin':
         t2s_limited, s0_limited, t2s_full, s0_full = fit_loglinear(
-            data_masked, tes, adaptive_mask_masked)
+            data_masked, tes, adaptive_mask_masked, report=report)
     elif fittype == 'curvefit':
         t2s_limited, s0_limited, t2s_full, s0_full = fit_monoexponential(
-            data_masked, tes, adaptive_mask_masked)
+            data_masked, tes, adaptive_mask_masked, report=report)
     else:
         raise ValueError('Unknown fittype option: {}'.format(fittype))
 
@@ -355,12 +393,14 @@ def fit_decay_ts(data, tes, mask, adaptive_mask, fittype):
     t2s_full_ts = np.copy(t2s_limited_ts)
     s0_full_ts = np.copy(t2s_limited_ts)
 
+    report = True
     for vol in range(n_vols):
         t2s_limited, s0_limited, t2s_full, s0_full = fit_decay(
-            data[:, :, vol][:, :, None], tes, mask, adaptive_mask, fittype)
+            data[:, :, vol][:, :, None], tes, mask, adaptive_mask, fittype, report=report)
         t2s_limited_ts[:, vol] = t2s_limited
         s0_limited_ts[:, vol] = s0_limited
         t2s_full_ts[:, vol] = t2s_full
         s0_full_ts[:, vol] = s0_full
+        report = False
 
     return t2s_limited_ts, s0_limited_ts, t2s_full_ts, s0_full_ts

tedana/workflows/tedana.py

@@ -355,18 +355,17 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
     # Removing handlers after basicConfig doesn't work, so we use filters
     # for the relevant handlers themselves.
     log_handler.addFilter(ContextFilter())
+    logging.root.addHandler(log_handler)
     sh = logging.StreamHandler()
     sh.addFilter(ContextFilter())
+    logging.root.addHandler(sh)
 
     if quiet:
-        logging.basicConfig(level=logging.WARNING,
-                            handlers=[log_handler, sh])
+        logging.root.setLevel(logging.WARNING)
     elif debug:
-        logging.basicConfig(level=logging.DEBUG,
-                            handlers=[log_handler, sh])
+        logging.root.setLevel(logging.DEBUG)
     else:
-        logging.basicConfig(level=logging.INFO,
-                            handlers=[log_handler, sh])
+        logging.root.setLevel(logging.INFO)
 
     # Loggers for report and references
     rep_handler = logging.FileHandler(repname)
@@ -375,7 +374,7 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
     ref_handler.setFormatter(text_formatter)
     RepLGR.setLevel(logging.INFO)
     RepLGR.addHandler(rep_handler)
-    RepLGR.setLevel(logging.INFO)
+    RefLGR.setLevel(logging.INFO)
     RefLGR.addHandler(ref_handler)
 
     LGR.info('Using output directory: {}'.format(out_dir))
@@ -665,10 +664,16 @@ def tedana_workflow(data, tes, out_dir='.', mask=None,
     report += '\n\nReferences:\n\n' + references
     with open(repname, 'w') as fo:
         fo.write(report)
-    os.remove(refname)
 
-    for handler in logging.root.handlers[:]:
-        logging.root.removeHandler(handler)
+    log_handler.close()
+    logging.root.removeHandler(log_handler)
+    sh.close()
+    logging.root.removeHandler(sh)
+    for local_logger in (RefLGR, RepLGR):
+        for handler in local_logger.handlers[:]:
+            handler.close()
+            local_logger.removeHandler(handler)
+    os.remove(refname)
 
 
 def _main(argv=None):