Merge pull request #2239 from emdupre/fix/meepi-t1-trans

REF/FIX: Simplify transform aggregation in resampling, pass identity transforms for multi-echo cases

Author: effigies

fmriprep/workflows/bold/base.py

@@ -340,8 +340,6 @@ def init_func_preproc_wf(bold_file):
     # apply BOLD registration to T1w
     bold_t1_trans_wf = init_bold_t1_trans_wf(name='bold_t1_trans_wf',
                                              freesurfer=freesurfer,
-                                             use_fieldwarp=bool(fmaps),
-                                             multiecho=multiecho,
                                              mem_gb=mem_gb['resampled'],
                                              omp_nthreads=omp_nthreads,
                                              use_compression=False)
@@ -398,10 +396,12 @@ def init_func_preproc_wf(bold_file):
                                        debug=config.execution.debug)
 
     # MULTI-ECHO EPI DATA #############################################
-    if multiecho:
+    if multiecho:  # instantiate relevant interfaces, imports
         from niworkflows.func.util import init_skullstrip_bold_wf
         skullstrip_bold_wf = init_skullstrip_bold_wf(name='skullstrip_bold_wf')
 
+        split_opt_comb = bold_split.clone(name='split_opt_comb')
+
         inputnode.inputs.bold_file = ref_file  # Replace reference w first echo
 
         join_echos = pe.JoinNode(niu.IdentityInterface(fields=['bold_files']),
@@ -415,13 +415,6 @@ def init_func_preproc_wf(bold_file):
                                        omp_nthreads=omp_nthreads,
                                        name='bold_t2smap_wf')
 
-        workflow.connect([
-            (skullstrip_bold_wf, join_echos, [
-                ('outputnode.skull_stripped_file', 'bold_files')]),
-            (join_echos, bold_t2s_wf, [
-                ('bold_files', 'inputnode.bold_file')]),
-        ])
-
     # MAIN WORKFLOW STRUCTURE #######################################################
     workflow.connect([
         (inputnode, t1w_brain, [('t1w_preproc', 'in_file'),
@@ -450,8 +443,6 @@ def init_func_preproc_wf(bold_file):
             ('t1w_aparc', 'inputnode.t1w_aparc')]),
         (t1w_brain, bold_t1_trans_wf, [
             ('out_file', 'inputnode.t1w_brain')]),
-        # unused if multiecho, but this is safe
-        (bold_hmc_wf, bold_t1_trans_wf, [('outputnode.xforms', 'inputnode.hmc_xforms')]),
         (bold_reg_wf, outputnode, [
             ('outputnode.itk_bold_to_t1', 'bold2anat_xfm'),
             ('outputnode.itk_t1_to_bold', 'anat2bold_xfm')]),
@@ -470,7 +461,6 @@ def init_func_preproc_wf(bold_file):
             ('outputnode.ref_image_brain', 'inputnode.epi_brain'),
             ('outputnode.bold_mask', 'inputnode.epi_mask')]),
         (bold_sdc_wf, bold_t1_trans_wf, [
-            ('outputnode.out_warp', 'inputnode.fieldwarp'),
             ('outputnode.epi_mask', 'inputnode.ref_bold_mask'),
             ('outputnode.epi_brain', 'inputnode.ref_bold_brain')]),
         (bold_sdc_wf, bold_bold_trans_wf, [
@@ -516,20 +506,34 @@ def init_func_preproc_wf(bold_file):
                 ('outputnode.bold', 'inputnode.bold')]),
             (bold_split, bold_t1_trans_wf, [
                 ('out_files', 'inputnode.bold_split')]),
+            (bold_hmc_wf, bold_t1_trans_wf, [
+                ('outputnode.xforms', 'inputnode.hmc_xforms')]),
+            (bold_sdc_wf, bold_t1_trans_wf, [
+                ('outputnode.out_warp', 'inputnode.fieldwarp')])
         ])
-    else:  # for meepi, create and use optimal combination
+    else:  # for meepi, use optimal combination
         workflow.connect([
             # update name source for optimal combination
             (inputnode, func_derivatives_wf, [
                 (('bold_file', combine_meepi_source), 'inputnode.source_file')]),
             (bold_bold_trans_wf, skullstrip_bold_wf, [
                 ('outputnode.bold', 'inputnode.in_file')]),
+            (skullstrip_bold_wf, join_echos, [
+                ('outputnode.skull_stripped_file', 'bold_files')]),
+            (join_echos, bold_t2s_wf, [
+                ('bold_files', 'inputnode.bold_file')]),
             (bold_t2s_wf, bold_confounds_wf, [
                 ('outputnode.bold', 'inputnode.bold')]),
-            (bold_t2s_wf, bold_t1_trans_wf, [
-                ('outputnode.bold', 'inputnode.bold_split')]),
+            (bold_t2s_wf, split_opt_comb, [
+                ('outputnode.bold', 'in_file')]),
+            (split_opt_comb, bold_t1_trans_wf, [
+                ('out_files', 'inputnode.bold_split')]),
         ])
 
+        # Already applied in bold_bold_trans_wf, which inputs to bold_t2s_wf
+        bold_t1_trans_wf.inputs.inputnode.fieldwarp = 'identity'
+        bold_t1_trans_wf.inputs.inputnode.hmc_xforms = 'identity'
+
     if fmaps:
         from sdcflows.workflows.outputs import init_sdc_unwarp_report_wf
         # Report on BOLD correction
@@ -628,7 +632,6 @@ def init_func_preproc_wf(bold_file):
             spaces=spaces,
             name='bold_std_trans_wf',
             use_compression=not config.execution.low_mem,
-            use_fieldwarp=bool(fmaps),
         )
         workflow.connect([
             (inputnode, bold_std_trans_wf, [
@@ -637,14 +640,10 @@ def init_func_preproc_wf(bold_file):
                 ('bold_file', 'inputnode.name_source'),
                 ('t1w_aseg', 'inputnode.bold_aseg'),
                 ('t1w_aparc', 'inputnode.bold_aparc')]),
-            (bold_hmc_wf, bold_std_trans_wf, [
-                ('outputnode.xforms', 'inputnode.hmc_xforms')]),
             (bold_reg_wf, bold_std_trans_wf, [
                 ('outputnode.itk_bold_to_t1', 'inputnode.itk_bold_to_t1')]),
             (bold_bold_trans_wf, bold_std_trans_wf, [
                 ('outputnode.bold_mask', 'inputnode.bold_mask')]),
-            (bold_sdc_wf, bold_std_trans_wf, [
-                ('outputnode.out_warp', 'inputnode.fieldwarp')]),
             (bold_std_trans_wf, outputnode, [('outputnode.bold_std', 'bold_std'),
                                              ('outputnode.bold_std_ref', 'bold_std_ref'),
                                              ('outputnode.bold_mask_std', 'bold_mask_std')]),
@@ -664,18 +663,22 @@ def init_func_preproc_wf(bold_file):
         if not multiecho:
             workflow.connect([
                 (bold_split, bold_std_trans_wf, [
-                    ('out_files', 'inputnode.bold_split')])
+                    ('out_files', 'inputnode.bold_split')]),
+                (bold_sdc_wf, bold_std_trans_wf, [
+                    ('outputnode.out_warp', 'inputnode.fieldwarp')]),
+                (bold_hmc_wf, bold_std_trans_wf, [
+                    ('outputnode.xforms', 'inputnode.hmc_xforms')]),
             ])
         else:
-            split_opt_comb = bold_split.clone(name='split_opt_comb')
             workflow.connect([
-                (bold_t2s_wf, split_opt_comb, [
-                    ('outputnode.bold', 'in_file')]),
                 (split_opt_comb, bold_std_trans_wf, [
-                    ('out_files', 'inputnode.bold_split')
-                ])
+                    ('out_files', 'inputnode.bold_split')])
             ])
 
+            # Already applied in bold_bold_trans_wf, which inputs to bold_t2s_wf
+            bold_std_trans_wf.inputs.inputnode.fieldwarp = 'identity'
+            bold_std_trans_wf.inputs.inputnode.hmc_xforms = 'identity'
+
         # func_derivatives_wf internally parametrizes over snapshotted spaces.
         workflow.connect([
             (bold_std_trans_wf, func_derivatives_wf, [
@@ -693,7 +696,6 @@ def init_func_preproc_wf(bold_file):
                 mem_gb=mem_gb['resampled'],
                 metadata=metadata,
                 omp_nthreads=omp_nthreads,
-                use_fieldwarp=bool(fmaps),
                 err_on_aroma_warn=config.workflow.aroma_err_on_warn,
                 aroma_melodic_dim=config.workflow.aroma_melodic_dim,
                 name='ica_aroma_wf')

fmriprep/workflows/bold/confounds.py

@@ -548,7 +548,6 @@ def init_ica_aroma_wf(
     err_on_aroma_warn=False,
     name='ica_aroma_wf',
     susan_fwhm=6.0,
-    use_fieldwarp=True,
 ):
     """
     Build a workflow that runs `ICA-AROMA`_.
@@ -601,8 +600,6 @@ def init_ica_aroma_wf(
     susan_fwhm : :obj:`float`
         Kernel width (FWHM in mm) for the smoothing step with
         FSL ``susan`` (default: 6.0mm)
-    use_fieldwarp : :obj:`bool`
-        Include SDC warp in single-shot transform from BOLD to MNI
     err_on_aroma_warn : :obj:`bool`
         Do not fail on ICA-AROMA errors
     aroma_melodic_dim : :obj:`int`
@@ -628,8 +625,6 @@ def init_ica_aroma_wf(
         BOLD series mask in template space
     hmc_xforms
         List of affine transforms aligning each volume to ``ref_image`` in ITK format
-    fieldwarp
-        a :abbr:`DFM (displacements field map)` in ITK format
     movpar_file
         SPM-formatted motion parameters file
 

fmriprep/workflows/bold/registration.py

@@ -175,8 +175,8 @@ def _bold_reg_suffix(fallback, freesurfer):
     return workflow
 
 
-def init_bold_t1_trans_wf(freesurfer, mem_gb, omp_nthreads, multiecho=False, use_fieldwarp=False,
-                          use_compression=True, name='bold_t1_trans_wf'):
+def init_bold_t1_trans_wf(freesurfer, mem_gb, omp_nthreads, use_compression=True,
+                          name='bold_t1_trans_wf'):
     """
     Co-register the reference BOLD image to T1w-space.
 
@@ -196,10 +196,6 @@ def init_bold_t1_trans_wf(freesurfer, mem_gb, omp_nthreads, multiecho=False, use
     ----------
     freesurfer : :obj:`bool`
         Enable FreeSurfer functional registration (bbregister)
-    use_fieldwarp : :obj:`bool`
-        Include SDC warp in single-shot transform from BOLD to T1
-    multiecho : :obj:`bool`
-        If multiecho data was supplied, HMC already performed
     mem_gb : :obj:`float`
         Size of BOLD file in GB
     omp_nthreads : :obj:`int`
@@ -327,38 +323,18 @@ def init_bold_t1_trans_wf(freesurfer, mem_gb, omp_nthreads, multiecho=False, use
     # Generate a reference on the target T1w space
     gen_final_ref = init_bold_reference_wf(omp_nthreads, pre_mask=True)
 
-    if not multiecho:
-        # Merge transforms placing the head motion correction last
-        nforms = 2 + int(use_fieldwarp)
-        merge_xforms = pe.Node(niu.Merge(nforms), name='merge_xforms',
-                               run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB)
-        if use_fieldwarp:
-            workflow.connect([
-                (inputnode, merge_xforms, [('fieldwarp', 'in2')])
-            ])
-
-        workflow.connect([
-            # merge transforms
-            (inputnode, merge_xforms, [
-                ('hmc_xforms', 'in%d' % nforms),
-                ('itk_bold_to_t1', 'in1')]),
-            (merge_xforms, bold_to_t1w_transform, [('out', 'transforms')]),
-            (inputnode, bold_to_t1w_transform, [('bold_split', 'input_image')]),
-        ])
-
-    else:
-        from nipype.interfaces.fsl import Split as FSLSplit
-        bold_split = pe.Node(FSLSplit(dimension='t'), name='bold_split',
-                             mem_gb=DEFAULT_MEMORY_MIN_GB)
-
-        workflow.connect([
-            (inputnode, bold_split, [('bold_split', 'in_file')]),
-            (bold_split, bold_to_t1w_transform, [('out_files', 'input_image')]),
-            (inputnode, bold_to_t1w_transform, [('itk_bold_to_t1', 'transforms')]),
-        ])
+    # Merge transforms placing the head motion correction last
+    merge_xforms = pe.Node(niu.Merge(3), name='merge_xforms',
+                           run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB)
 
     workflow.connect([
         (inputnode, merge, [('name_source', 'header_source')]),
+        (inputnode, merge_xforms, [
+            ('hmc_xforms', 'in3'),  # May be 'identity' if HMC already applied
+            ('fieldwarp', 'in2'),   # May be 'identity' if SDC already applied
+            ('itk_bold_to_t1', 'in1')]),
+        (inputnode, bold_to_t1w_transform, [('bold_split', 'input_image')]),
+        (merge_xforms, bold_to_t1w_transform, [('out', 'transforms')]),
         (gen_ref, bold_to_t1w_transform, [('out_file', 'reference_image')]),
         (bold_to_t1w_transform, merge, [('out_files', 'in_files')]),
         (merge, gen_final_ref, [('out_file', 'inputnode.bold_file')]),

fmriprep/workflows/bold/resampling.py

@@ -13,7 +13,6 @@
 
 from nipype.pipeline import engine as pe
 from nipype.interfaces import utility as niu, freesurfer as fs
-from nipype.interfaces.fsl import Split as FSLSplit
 import nipype.interfaces.workbench as wb
 
 
@@ -164,7 +163,6 @@ def init_bold_std_trans_wf(
     spaces,
     name='bold_std_trans_wf',
     use_compression=True,
-    use_fieldwarp=False,
 ):
     """
     Sample fMRI into standard space with a single-step resampling of the original BOLD series.
@@ -215,8 +213,6 @@ def init_bold_std_trans_wf(
         Name of workflow (default: ``bold_std_trans_wf``)
     use_compression : :obj:`bool`
         Save registered BOLD series as ``.nii.gz``
-    use_fieldwarp : :obj:`bool`
-        Include SDC warp in single-shot transform from BOLD to MNI
 
     Inputs
     ------
@@ -335,13 +331,8 @@ def init_bold_std_trans_wf(
     mask_merge_tfms = pe.Node(niu.Merge(2), name='mask_merge_tfms', run_without_submitting=True,
                               mem_gb=DEFAULT_MEMORY_MIN_GB)
 
-    nxforms = 3 + use_fieldwarp
-    merge_xforms = pe.Node(niu.Merge(nxforms), name='merge_xforms',
+    merge_xforms = pe.Node(niu.Merge(4), name='merge_xforms',
                            run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB)
-    workflow.connect([(inputnode, merge_xforms, [('hmc_xforms', 'in%d' % nxforms)])])
-
-    if use_fieldwarp:
-        workflow.connect([(inputnode, merge_xforms, [('fieldwarp', 'in3')])])
 
     bold_to_std_transform = pe.Node(
         MultiApplyTransforms(interpolation="LanczosWindowedSinc", float=True, copy_dtype=True),
@@ -361,8 +352,9 @@ def init_bold_std_trans_wf(
                                  ('templates', 'keys')]),
         (inputnode, mask_std_tfm, [('bold_mask', 'input_image')]),
         (inputnode, gen_ref, [(('bold_split', _first), 'moving_image')]),
-        (inputnode, merge_xforms, [
-            (('itk_bold_to_t1', _aslist), 'in2')]),
+        (inputnode, merge_xforms, [('hmc_xforms', 'in4'),
+                                   ('fieldwarp', 'in3'),
+                                   (('itk_bold_to_t1', _aslist), 'in2')]),
         (inputnode, merge, [('name_source', 'header_source')]),
         (inputnode, mask_merge_tfms, [(('itk_bold_to_t1', _aslist), 'in2')]),
         (inputnode, bold_to_std_transform, [('bold_split', 'input_image')]),
@@ -431,7 +423,6 @@ def init_bold_preproc_trans_wf(mem_gb, omp_nthreads,
                                name='bold_preproc_trans_wf',
                                use_compression=True,
                                use_fieldwarp=False,
-                               split_file=False,
                                interpolation='LanczosWindowedSinc'):
     """
     Resample in native (original) space.
@@ -459,9 +450,6 @@ def init_bold_preproc_trans_wf(mem_gb, omp_nthreads,
         Save registered BOLD series as ``.nii.gz``
     use_fieldwarp : :obj:`bool`
         Include SDC warp in single-shot transform from BOLD to MNI
-    split_file : :obj:`bool`
-        Whether the input file should be splitted (it is a 4D file)
-        or it is a list of 3D files (default ``False``, do not split)
     interpolation : :obj:`str`
         Interpolation type to be used by ANTs' ``applyTransforms``
         (default ``'LanczosWindowedSinc'``)
@@ -518,19 +506,26 @@ def init_bold_preproc_trans_wf(mem_gb, omp_nthreads,
         niu.IdentityInterface(fields=['bold', 'bold_mask', 'bold_ref', 'bold_ref_brain']),
         name='outputnode')
 
+    merge_xforms = pe.Node(niu.Merge(2), name='merge_xforms',
+                           run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB)
+
     bold_transform = pe.Node(
         MultiApplyTransforms(interpolation=interpolation, float=True, copy_dtype=True),
         name='bold_transform', mem_gb=mem_gb * 3 * omp_nthreads, n_procs=omp_nthreads)
 
-    merge = pe.Node(Merge(compress=use_compression), name='merge',
-                    mem_gb=mem_gb * 3)
+    merge = pe.Node(Merge(compress=use_compression), name='merge', mem_gb=mem_gb * 3)
 
     # Generate a new BOLD reference
     bold_reference_wf = init_bold_reference_wf(omp_nthreads=omp_nthreads)
     bold_reference_wf.__desc__ = None  # Unset description to avoid second appearance
 
     workflow.connect([
+        (inputnode, merge_xforms, [('fieldwarp', 'in1'),
+                                   ('hmc_xforms', 'in2')]),
+        (inputnode, bold_transform, [('bold_file', 'input_image'),
+                                     (('bold_file', _first), 'reference_image')]),
         (inputnode, merge, [('name_source', 'header_source')]),
+        (merge_xforms, bold_transform, [('out', 'transforms')]),
         (bold_transform, merge, [('out_files', 'in_files')]),
         (merge, bold_reference_wf, [('out_file', 'inputnode.bold_file')]),
         (merge, outputnode, [('out_file', 'bold')]),
@@ -540,37 +535,6 @@ def init_bold_preproc_trans_wf(mem_gb, omp_nthreads,
             ('outputnode.bold_mask', 'bold_mask')]),
     ])
 
-    # Input file is not splitted
-    if split_file:
-        bold_split = pe.Node(FSLSplit(dimension='t'), name='bold_split',
-                             mem_gb=mem_gb * 3)
-        workflow.connect([
-            (inputnode, bold_split, [('bold_file', 'in_file')]),
-            (bold_split, bold_transform, [
-                ('out_files', 'input_image'),
-                (('out_files', _first), 'reference_image'),
-            ])
-        ])
-    else:
-        workflow.connect([
-            (inputnode, bold_transform, [('bold_file', 'input_image'),
-                                         (('bold_file', _first), 'reference_image')]),
-        ])
-
-    if use_fieldwarp:
-        merge_xforms = pe.Node(niu.Merge(2), name='merge_xforms',
-                               run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB)
-        workflow.connect([
-            (inputnode, merge_xforms, [('fieldwarp', 'in1'),
-                                       ('hmc_xforms', 'in2')]),
-            (merge_xforms, bold_transform, [('out', 'transforms')]),
-        ])
-    else:
-        def _aslist(val):
-            return [val]
-        workflow.connect([
-            (inputnode, bold_transform, [(('hmc_xforms', _aslist), 'transforms')]),
-        ])
     return workflow
 
 

setup.cfg

@@ -34,7 +34,7 @@ install_requires =
     pybids >= 0.10.2
     pyyaml
     requests
-    sdcflows ~= 1.3.1
+    sdcflows @ git+https://github.com/nipreps/sdcflows.git@master
     smriprep @ git+https://github.com/nipreps/smriprep.git@master
     tedana >= 0.0.9a1, < 0.0.10
     templateflow ~= 0.6