[RTM] ENH: Merge T1w images to unbiased template if N < 3 or --longitudinal is passed

docs/environment.yml

@@ -28,4 +28,4 @@ dependencies:
     - svgutils
     - nitime
     - nilearn
-    - niworkflows>=0.1.4
+    - git+https://github.com/poldracklab/niworkflows.git@0a71b0f16acec10520b6eb824649d19519e65b59#egg=niworkflows-0.1.5-dev

docs/links.rst

@@ -6,7 +6,7 @@
 .. _workflows: workflows.html
 .. _FreeSurfer: https://surfer.nmr.mgh.harvard.edu/
 .. _`submillimeter reconstruction`: https://surfer.nmr.mgh.harvard.edu/fswiki/SubmillimeterRecon
-.. _`mri_robust_register`: https://surfer.nmr.mgh.harvard.edu/fswiki/mri_robust_template
+.. _`mri_robust_template`: https://surfer.nmr.mgh.harvard.edu/fswiki/mri_robust_template
 .. _GIFTI: https://www.nitrc.org/projects/gifti/
 .. _`Connectome Workbench`: https://www.humanconnectome.org/software/connectome-workbench.html
 .. _`fmriprep-docker`: https://pypi.python.org/pypi/fmriprep-docker

docs/workflows.rst

@@ -20,6 +20,7 @@ slice-timing information and no fieldmap acquisitions):
     wf = init_single_subject_wf(subject_id='test',
                                 name='single_subject_wf',
                                 task_id='',
+                                longitudinal=False,
                                 omp_nthreads=1,
                                 freesurfer=True,
                                 reportlets_dir='.',
@@ -60,13 +61,14 @@ T1w/T2w preprocessing
                                              'template', 'fsaverage5'],
                               skull_strip_ants=True,
                               freesurfer=True,
+                              longitudinal=False,
                               debug=False,
                               hires=True)
 
 The anatomical sub-workflow begins by constructing a template image by
 :ref:`conforming <conformation>` any T1-weighted images to RAS orientation and
 a common voxel size, and, in the case of multiple images, merges them into a
-single template using `mri_robust_register`_.
+single template (see `Longitudinal processing`_).
 This template is then skull-stripped, and the white matter/gray
 matter/cerebrospinal fluid segments are found.
 Finally, a non-linear registration to the MNI template space is estimated.
@@ -86,6 +88,19 @@ Finally, a non-linear registration to the MNI template space is estimated.
 
     Animation showing T1w to MNI normalization (ANTs)
 
+Longitudinal processing
+~~~~~~~~~~~~~~~~~~~~~~~
+In the case of multiple sessions, T1w images are merged into a single template
+image using FreeSurfer's `mri_robust_template`_.
+This template may be *unbiased*, or equidistant from all source images, or
+aligned to the first image (determined lexicographically by session label).
+For two images, the additional cost of estimating an unbiased template is
+trivial and is the default behavior, but, for greater than two images, the cost
+can be a slowdown of an order of magnitude.
+Therefore, in the case of three or more images, ``fmriprep`` constructs
+templates aligned to the first image, unless passed the ``--longitudinal``
+flag, which forces the estimation of an unbiased template.
+
 Surface preprocessing
 ~~~~~~~~~~~~~~~~~~~~~
 :mod:`fmriprep.workflows.anatomical.init_surface_recon_wf`

fmriprep/cli/run.py

@@ -81,6 +81,9 @@ def get_parser():
         choices=['fieldmaps', 'slicetiming'],
         help='ignore selected aspects of the input dataset to disable corresponding '
              'parts of the workflow')
+    g_conf.add_argument(
+        '--longitudinal', action='store_true',
+        help='treat dataset as longitudinal - may increase runtime')
     g_conf.add_argument('--bold2t1w-dof', action='store', default=9, choices=[6, 9, 12], type=int,
                         help='Degrees of freedom when registering BOLD to T1w images. '
                              '9 (rotation, translation, and scaling) is used by '
@@ -253,6 +256,7 @@ def create_workflow(opts):
                                        ignore=opts.ignore,
                                        debug=opts.debug,
                                        anat_only=opts.anat_only,
+                                       longitudinal=opts.longitudinal,
                                        omp_nthreads=omp_nthreads,
                                        skull_strip_ants=opts.skull_strip_ants,
                                        reportlets_dir=reportlets_dir,

fmriprep/workflows/anatomical.py

@@ -32,7 +32,7 @@
 
 #  pylint: disable=R0914
 def init_anat_preproc_wf(skull_strip_ants, output_spaces, template, debug, freesurfer,
-                         omp_nthreads, hires, reportlets_dir, output_dir,
+                         longitudinal, omp_nthreads, hires, reportlets_dir, output_dir,
                          name='anat_preproc_wf'):
     """T1w images preprocessing pipeline"""
 
@@ -73,12 +73,11 @@ def bidsinfo(in_file):
     t1_merge = pe.Node(
         # StructuralReference is fs.RobustTemplate if > 1 volume, copying otherwise
         StructuralReference(auto_detect_sensitivity=True,
-                            initial_timepoint=1,
-                            fixed_timepoint=True,     # Align to first image
+                            initial_timepoint=1,      # For deterministic behavior
                             intensity_scaling=True,   # 7-DOF (rigid + intensity)
-                            no_iteration=True,
                             subsample_threshold=200,
-                            ), name='t1_merge')
+                            ),
+        name='t1_merge')
 
     # 2. T1 Bias Field Correction
     # Bias field correction is handled in skull strip workflows.
@@ -129,11 +128,23 @@ def bidsinfo(in_file):
         name='mni_tpms'
     )
 
+    def set_threads(in_list, maximum):
+        return min(len(in_list), maximum)
+
+    def len_above_thresh(in_list, threshold, longitudinal):
+        if longitudinal:
+            return False
+        return len(in_list) > threshold
+
     workflow.connect([
         (inputnode, bids_info, [(('t1w', fix_multi_T1w_source_name), 'in_file')]),
         (inputnode, t1_conform, [('t1w', 't1w_list')]),
-        (t1_conform, t1_merge, [('t1w_list', 'in_files'),
-                                (('t1w_list', add_suffix, '_template'), 'out_file')]),
+        (t1_conform, t1_merge, [
+            ('t1w_list', 'in_files'),
+            (('t1w_list', set_threads, omp_nthreads), 'num_threads'),
+            (('t1w_list', len_above_thresh, 2, longitudinal), 'fixed_timepoint'),
+            (('t1w_list', len_above_thresh, 2, longitudinal), 'no_iteration'),
+            (('t1w_list', add_suffix, '_template'), 'out_file')]),
         (t1_merge, skullstrip_wf, [('out_file', 'inputnode.in_file')]),
         (skullstrip_wf, t1_seg, [('outputnode.out_file', 'in_files')]),
         (skullstrip_wf, outputnode, [('outputnode.bias_corrected', 't1_preproc'),

fmriprep/workflows/base.py

@@ -26,7 +26,7 @@
 
 
 def init_fmriprep_wf(subject_list, task_id, run_uuid,
-                     ignore, debug, anat_only, omp_nthreads,
+                     ignore, debug, anat_only, longitudinal, omp_nthreads,
                      skull_strip_ants, reportlets_dir, output_dir, bids_dir,
                      freesurfer, output_spaces, template, hires,
                      bold2t1w_dof, fmap_bspline, fmap_demean, use_syn, force_syn,
@@ -48,6 +48,7 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid,
                                                    ignore=ignore,
                                                    debug=debug,
                                                    anat_only=anat_only,
+                                                   longitudinal=longitudinal,
                                                    omp_nthreads=omp_nthreads,
                                                    skull_strip_ants=skull_strip_ants,
                                                    reportlets_dir=reportlets_dir,
@@ -81,7 +82,7 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid,
 
 
 def init_single_subject_wf(subject_id, task_id, name,
-                           ignore, debug, anat_only, omp_nthreads,
+                           ignore, debug, anat_only, longitudinal, omp_nthreads,
                            skull_strip_ants, reportlets_dir, output_dir, bids_dir,
                            freesurfer, output_spaces, template, hires,
                            bold2t1w_dof, fmap_bspline, fmap_demean, use_syn, force_syn,
@@ -122,6 +123,7 @@ def init_single_subject_wf(subject_id, task_id, name,
                                            output_spaces=output_spaces,
                                            template=template,
                                            debug=debug,
+                                           longitudinal=longitudinal,
                                            omp_nthreads=omp_nthreads,
                                            freesurfer=freesurfer,
                                            hires=hires,

fmriprep/workflows/confounds.py

@@ -29,7 +29,7 @@ def init_discover_wf(bold_file_size_gb, use_aroma, ignore_aroma_err, metadata,
     Calculates segment regressors and aCompCor
         from the fMRI and a white matter/gray matter/CSF segmentation ('inputnode.t1_seg'), after
         applying the transform to the images. Transforms should be fsl-formatted.
-    Calculates noise components identified from ICA_AROMA (if ``use_aroma=True``)
+    Calculates noise components identified from ICA-AROMA (if ``use_aroma=True``)
     Saves the confounds in a file ('outputnode.confounds_file')'''
 
     inputnode = pe.Node(utility.IdentityInterface(
@@ -310,7 +310,7 @@ def aroma_add_header_func(np_arr, col_base, comp_nums):
 
         return os.path.abspath(str(col_base) + "AROMAConfounds.tsv")
 
-    # load the txt files from ICA_AROMA
+    # load the txt files from ICA-AROMA
     melodic_mix = os.path.join(ica_out_dir, 'melodic.ica/melodic_mix')
     motion_ics = os.path.join(ica_out_dir, 'classified_motion_ICs.txt')
 
@@ -378,8 +378,8 @@ def init_ica_aroma_wf(name='ica_aroma_wf', ignore_aroma_err=False):
 
     Steps:
     1) smooth data using SUSAN
-    2) run melodic outside of ICA_AROMA to generate the report
-    3) run ICA_AROMA
+    2) run melodic outside of ICA-AROMA to generate the report
+    3) run ICA-AROMA
     4) print identified motion components (aggressive) to tsv
     5) pass classified_motion_ICs and melodic_mix for user to complete nonaggr denoising
     '''
@@ -462,7 +462,7 @@ def getusans_func(image, thresh):
         (smooth, ica_aroma, [('smoothed_file', 'in_file')]),
         (inputnode, ica_aroma, [('movpar_file', 'motion_parameters')]),
         (melodic, ica_aroma, [('out_dir', 'melodic_dir')]),
-        # geneerate tsvs from ICA_AROMA
+        # generate tsvs from ICA-AROMA
         (ica_aroma, ica_aroma_confound_extraction, [('out_dir', 'ica_out_dir')]),
         # output for processing and reporting
         (ica_aroma_confound_extraction, outputnode, [('aroma_confounds', 'aroma_confounds'),

requirements.txt

@@ -1 +1 @@
-git+https://github.com/poldracklab/niworkflows.git@1d67e9463f8b7871d75e6f41e4b2c9ccdca8bf93#egg=niworkflows-0.1.5-dev
+git+https://github.com/poldracklab/niworkflows.git@0a71b0f16acec10520b6eb824649d19519e65b59#egg=niworkflows-0.1.5-dev

test/workflows/test_base.py

@@ -17,6 +17,7 @@ def test_single_subject_wf(self, _):
                                          ignore=[],
                                          debug=False,
                                          anat_only=False,
+                                         longitudinal=False,
                                          omp_nthreads=1,
                                          skull_strip_ants=False,
                                          reportlets_dir='.',