FIX: Simplify STC logic for too short BOLD series

.circleci/config.yml

@@ -839,7 +839,7 @@ jobs:
                 --config $PWD/nipype.cfg -w /tmp/${DATASET}/work \
                 /tmp/data/${DATASET} /tmp/${DATASET}/derivatives participant \
                 ${FASTRACK_ARG} \
-                --fs-no-reconall --use-syn-sdc \
+                --fs-no-reconall --use-syn-sdc --ignore slicetiming \
                 --dummy-scans 1 --sloppy --write-graph \
                 --output-spaces MNI152NLin2009cAsym \
                 --mem_mb 4096 --nthreads 2 -vv

docs/faq.rst

@@ -1,8 +1,14 @@
 .. include:: links.rst
 
-========================
-FAQ, Tips, and Tricks
-========================
+================================
+FAQ - Frequently Asked Questions
+================================
+
+.. contents::
+    :local:
+    :depth: 1
+
+------------
 
 Should I run quality control of my images before running *fMRIPrep*?
 --------------------------------------------------------------------
@@ -118,7 +124,7 @@ In general, please be cautious of deleting files and mindful why a file may exis
 
 Running subjects in parallel
 ----------------------------
-When running several subjects in parallel, and depending on your settings, fMRIPrep may
+When running several subjects in parallel, and depending on your settings, *fMRIPrep* may
 fall into race conditions.
 A symptomatic output looks like: ::
 
@@ -127,9 +133,9 @@ A symptomatic output looks like: ::
 If you would like to run *fMRIPrep* in parallel on multiple subjects please use
 `this method <https://neurostars.org/t/updated-fmriprep-workaround-for-running-subjects-in-parallel/6677>`__.
 
-How much CPU time and RAM should I allocate for a typical fMRIPrep run?
------------------------------------------------------------------------
-The recommended way to run fMRIPrep is to process one subject per container instance. A typical preprocessing run
+How much CPU time and RAM should I allocate for a typical *fMRIPrep* run?
+-------------------------------------------------------------------------
+The recommended way to run *fMRIPrep* is to process one subject per container instance. A typical preprocessing run
 without surface processing with FreeSurfer can be completed in about 2 hours with 4 CPUs or in about 1 hour with 16 CPUs.
 More than 16 CPUs do not translate into faster processing times for a single subject. About 8GB of memory should be
 available for a single subject preprocessing run.
@@ -140,14 +146,14 @@ CPUs and 64 GB of physical memory:
 .. figure:: _static/fmriprep_benchmark.svg
 
 **Compute Time**: time in hours to complete the preprocessing for all subjects. **Physical Memory**: the maximum of RAM usage
-used across all fMRIPrep processes as reported by the HCP job manager. **Virtual Memory**: the maximum of virtual memory used
-across all fMRIPrep processes as reported by the HCP job manager. **Threads**: the maximum number of threads per process as
+used across all *fMRIPrep* processes as reported by the HCP job manager. **Virtual Memory**: the maximum of virtual memory used
+across all *fMRIPrep* processes as reported by the HCP job manager. **Threads**: the maximum number of threads per process as
 specified with ``–omp-nthreads`` in the fMRIPrep command.
 
-The above figure illustrates that processing 2 subjects in 2 fMRIPrep instances with 8 CPUs each is approximately as fast as
-processing 2 subjects in one fMRIPrep instance with 16 CPUs. However, on a distributed compute cluster, the two 8 CPU
+The above figure illustrates that processing 2 subjects in 2 *fMRIPrep* instances with 8 CPUs each is approximately as fast as
+processing 2 subjects in one *fMRIPrep* instance with 16 CPUs. However, on a distributed compute cluster, the two 8 CPU
 instances may be allocated faster than the single 16 CPU instance, thus completing faster in practice. If more than one
-subject is processed in a single fMRIPrep instance, then limiting the number of threads per process to roughly the
+subject is processed in a single *fMRIPrep* instance, then limiting the number of threads per process to roughly the
 number of CPUs divided by the number of subjects is most efficient.
 
 .. _upgrading:
@@ -183,8 +189,8 @@ For further details, please check `its documentation section <spaces.html#templa
 
 .. _tf_no_internet:
 
-How do you use TemplateFlow in the absence of access to the Internet?
----------------------------------------------------------------------
+How do you use *TemplateFlow* in the absence of access to the Internet?
+-----------------------------------------------------------------------
 This is a fairly common situation in :abbr:`HPCs (high-performance computing)`
 systems, where the so-called login nodes have access to the Internet but
 compute nodes are isolated, or in PC/laptop environments if you are traveling.
@@ -222,9 +228,9 @@ runtime environment is able to access the filesystem, at the location of your
 *TemplateFlow home* directory.
 If you are a Singularity user, please check out :ref:`singularity_tf`.
 
-How do I select only certain files to be input to fMRIPrep?
------------------------------------------------------------
-Using the ``--bids-filter-file`` flag, you can pass fMRIPrep a JSON file that
+How do I select only certain files to be input to *fMRIPrep*?
+-------------------------------------------------------------
+Using the ``--bids-filter-file`` flag, you can pass *fMRIPrep* a JSON file that
 describes a custom BIDS filter for selecting files with PyBIDS, with the syntax
 ``{<query>: {<entity>: <filter>, ...},...}``. For example::
 
@@ -315,3 +321,25 @@ Note that any discrepancies between the pre-indexed database and
 the BIDS dataset complicate the provenance of fMRIPrep derivatives.
 If `--bids-database-dir` is used, the referenced directory should be
 preserved for the sake of reporting and reproducibility.
+
+Error in slice timing correction: *insufficient length of BOLD data after discarding nonsteady-states*
+------------------------------------------------------------------------------------------------------
+Typically, the scanner will be in a *nonsteady state* during a few initial time points of the acquisition,
+until it stabilizes.
+These *nonsteady states* (also called *dummy* scans) typically show greater T1 contrast and higher average
+intensity, and therefore potentially are detrimental if used in the interpolation of slice timing corrections.
+Hence, *nonsteady states* are discarded by the slice timing correction tool (in this case, AFNI's ``3dTShift``).
+However, ``3dTShift`` requires that at least five (5) time points are present in the target series, after
+dismissing the initial *nonsteady states*.
+
+*fMRIPrep* estimates the number of *nonsteady states* within the pipeline, unless the parameter is provided
+by the user with the argument ``--dummy-scans <num>``.
+Either way, if the number of *nonsteady states* is, say 4, then the length of the BOLD series must be greater
+than 8.
+If you encounter this error, first check that the number of *nonsteady states* is not suspiciously large
+(it typically ranges from zero to five).
+Next, if the number of *nonsteady states* is reasonable, consider why your BOLD time series are so short
+and whether slice timing correction is appropriate under these conditions.
+Finally, you can either skip the slice-timing correction with the argument ``--ignore slicetiming`` or
+enforce a number of *nonsteady states* lower than the maximum for your data with ``--dummy-scans <num>``.
+Please note that both strategies will apply to all tasks and runs that are to be processed.

fmriprep/cli/parser.py

@@ -335,7 +335,7 @@ def _bids_filter(value):
         action="store",
         default=None,
         type=int,
-        help="Number of non steady state volumes.",
+        help="Number of nonsteady-state volumes.",
     )
     g_conf.add_argument(
         "--random-seed",

fmriprep/workflows/bold/base.py

@@ -246,11 +246,10 @@ def init_func_preproc_wf(bold_file):
         # If fieldmaps are not enabled, activate SyN-SDC in unforced (False) mode
         fmaps = {'syn': False}
 
-    # Short circuits: (True and True and (False or 'TooShort')) == 'TooShort'
+    # Check whether STC must/can be run
     run_stc = (
         bool(metadata.get("SliceTiming"))
         and 'slicetiming' not in config.workflow.ignore
-        and (_get_series_len(ref_file, config.workflow.dummy_scans) > 4 or "TooShort")
     )
 
     # Build workflow
@@ -411,7 +410,7 @@ def init_func_preproc_wf(bold_file):
     final_boldref_wf.__desc__ = None  # Unset description to avoid second appearance
 
     # SLICE-TIME CORRECTION (or bypass) #############################################
-    if run_stc is True:  # bool('TooShort') == True, so check True explicitly
+    if run_stc:
         bold_stc_wf = init_bold_stc_wf(name='bold_stc_wf', metadata=metadata)
         workflow.connect([
             (initial_boldref_wf, bold_stc_wf, [
@@ -905,18 +904,6 @@ def init_func_preproc_wf(bold_file):
     return workflow
 
 
-def _get_series_len(bold_fname, skip_vols):
-    from niworkflows.interfaces.registration import _get_vols_to_discard
-    img = nb.load(bold_fname)
-    if len(img.shape) < 4:
-        return 1
-
-    if skip_vols is None:
-        skip_vols = _get_vols_to_discard(img)
-
-    return img.shape[3] - skip_vols
-
-
 def _create_mem_gb(bold_fname):
     bold_size_gb = os.path.getsize(bold_fname) / (1024**3)
     bold_tlen = nb.load(bold_fname).shape[-1]

fmriprep/workflows/bold/stc.py

@@ -27,15 +27,31 @@
 .. autofunction:: init_bold_stc_wf
 
 """
+import nibabel as nb
 from nipype.pipeline import engine as pe
 from nipype.interfaces import utility as niu, afni
+from nipype.interfaces.base import isdefined
 
 from ... import config
 
 
 LOGGER = config.loggers.workflow
 
 
+class TShift(afni.TShift):
+    """Patched version of TShift implementing the "TooShort" behavior."""
+
+    def _pre_run_hook(self, runtime):
+        ignore = self.inputs.ignore if isdefined(self.inputs.ignore) else 0
+        ntsteps = nb.load(self.inputs.in_file).shape[3]
+        if ntsteps - ignore < 5:
+            raise RuntimeError(
+                f"Insufficient length of BOLD data ({ntsteps} time points) after "
+                f"discarding {ignore} nonsteady-state (or 'dummy') time points."
+            )
+        return runtime
+
+
 def init_bold_stc_wf(metadata, name='bold_stc_wf'):
     """
     Create a workflow for :abbr:`STC (slice-timing correction)`.
@@ -89,10 +105,10 @@ def init_bold_stc_wf(metadata, name='bold_stc_wf'):
 
     # It would be good to fingerprint memory use of afni.TShift
     slice_timing_correction = pe.Node(
-        afni.TShift(outputtype='NIFTI_GZ',
-                    tr='{}s'.format(metadata["RepetitionTime"]),
-                    slice_timing=metadata['SliceTiming'],
-                    slice_encoding_direction=metadata.get('SliceEncodingDirection', 'k')),
+        TShift(outputtype='NIFTI_GZ',
+               tr=f"{metadata['RepetitionTime']}s",
+               slice_timing=metadata['SliceTiming'],
+               slice_encoding_direction=metadata.get('SliceEncodingDirection', 'k')),
         name='slice_timing_correction')
 
     copy_xform = pe.Node(CopyXForm(), name='copy_xform', mem_gb=0.1)