[FIX] Replace FSL mcflirt with AFNI 3dVolReg

.circleci/config.yml

@@ -537,7 +537,7 @@ jobs:
                 -e FMRIPREP_DEV 1 \
                 --config $PWD/nipype.cfg -w /tmp/ds054/work \
                 /tmp/data/ds054 /tmp/ds054/derivatives participant \
-                --fs-no-reconall --sloppy \
+                --fs-no-reconall --sloppy --hmc-use-mcflirt \
                 --output-space T1w template \
                 --template-resampling-grid 2mm \
                 --mem_mb 4096 --nthreads 2 -vv
@@ -601,8 +601,8 @@ jobs:
           at: /tmp
       - restore_cache:
           keys:
-            - ds210-anat-v4-{{ epoch }}
-            - ds210-anat-v4-
+            - ds210-anat-v5-{{ epoch }}
+            - ds210-anat-v5-
       - run:
           name: Setting up test
           command: |
@@ -636,7 +636,7 @@ jobs:
                 --fs-no-reconall --sloppy --write-graph \
                 --mem_mb 4096 --nthreads 2 --anat-only -vv
       - save_cache:
-         key: ds210-anat-v4-{{ epoch }}
+         key: ds210-anat-v5-{{ epoch }}
          paths:
             - /tmp/ds210/work
             - /tmp/ds210/derivatives

docs/citing.rst

@@ -100,8 +100,8 @@ we recommend to include in your paper.
 
    <p style="font-style: italic;">
      Functional data was <span class="slicetime_text_true">slice time corrected using
-     <code>3dTshift</code> from AFNI v16.2.07 [11, RRID:SCR_005927]
-     and </span>motion corrected using <code>mcflirt</code> (FSL v5.0.9 [9]).
+     <code>3dTshift</code> (AFNI v16.2.07 [11, RRID:SCR_005927]) and </span>
+     motion-corrected using <code>3dVolreg</code> (AFNI v16.2.07 [11, RRID:SCR_005927]).
      <span class="SDC_text_TOPUP" style="display: none">Distortion correction was performed
      using an implementation of the TOPUP technique [10] using <code>3dQwarp</code> (AFNI v16.2.07 [11]).</span>
      <span class="SDC_text_FUGUE" style="display: none">Distortion correction was performed using fieldmaps

docs/outputs.rst

@@ -77,8 +77,8 @@ Derivatives related to EPI files are in the ``func`` subfolder.
 Volumetric output spaces include ``T1w`` and ``MNI152NLin2009cAsym`` (default).
 
 - ``*bold_space-<space>_brainmask.nii.gz`` Brain mask for EPI files, calculated by nilearn on the average EPI volume, post-motion correction
-- ``*bold_space-<space>_preproc.nii.gz`` Motion-corrected (using MCFLIRT for estimation and ANTs for interpolation) EPI file
-- (optional) ``*bold_space-<space>_variant-smoothAROMAnonaggr_preproc.nii.gz`` Motion-corrected (using MCFLIRT for estimation and ANTs for interpolation),
+- ``*bold_space-<space>_preproc.nii.gz`` Head-motion corrected EPI file
+- (optional) ``*bold_space-<space>_variant-smoothAROMAnonaggr_preproc.nii.gz`` Head-motion corrected,
   smoothed (6mm), and non-aggressively denoised (using AROMA) EPI file - currently produced only for the ``MNI152NLin2009cAsym`` space
 
 Surface output spaces include ``fsnative`` (full density subject-specific mesh),

docs/workflows.rst

@@ -51,6 +51,7 @@ is presented below:
         use_aroma=False,
         aroma_melodic_dim=None,
         ignore_aroma_err=False,
+        use_mcflirt=False,
     )
 
 
@@ -273,6 +274,7 @@ BOLD preprocessing
         use_aroma=False,
         aroma_melodic_dim=None,
         ignore_aroma_err=False,
+        use_mcflirt=False,
     )
 
 Preprocessing of :abbr:`BOLD (blood-oxygen level-dependent)` files is
@@ -334,11 +336,12 @@ Head-motion estimation
 
     from fmriprep.workflows.bold import init_bold_hmc_wf
     wf = init_bold_hmc_wf(
+        use_mcflirt=False,
         mem_gb=1,
         omp_nthreads=1)
 
 Using the previously :ref:`estimated reference scan <bold_ref>`,
-FSL ``mcflirt`` is used to estimate head-motion.
+AFNI's ``3dVolreg`` is used to estimate head-motion.
 As a result, one rigid-body transform with respect to
 the reference image is written for each :abbr:`BOLD (blood-oxygen level-dependent)`
 time-step.
@@ -548,7 +551,7 @@ Confounds estimation
         metadata={"RepetitionTime": 2.0,
                   "SliceTiming": [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]})
 
-Given a motion-corrected fMRI, a brain mask, ``mcflirt`` movement parameters and a
+Given a motion-corrected fMRI, a brain mask, head-motion parameters and a
 segmentation, the `discover_wf` sub-workflow calculates potential
 confounds per volume.
 

fmriprep/__about__.py

@@ -105,6 +105,7 @@
     'scikit-image',
     'versioneer',
     'pyyaml',
+    'transforms3d',
 ]
 
 LINKS_REQUIRES = [

fmriprep/cli/run.py

@@ -160,6 +160,10 @@ def get_parser():
         '--medial-surface-nan', required=False, action='store_true', default=False,
         help='Replace medial wall values with NaNs on functional GIFTI files. Only '
         'performed for GIFTI files mapped to a freesurfer subject (fsaverage or fsnative).')
+    g_conf.add_argument(
+        '--hmc-use-mcflirt', required=False, action='store_true', default=False,
+        help='Head-Motion Correction (HMC) - use FSL\'s ``mcflirt`` instead '
+        'of AFNI\'s ``3dVolreg``.')
 
     # ICA_AROMA options
     g_aroma = parser.add_argument_group('Specific options for running ICA_AROMA')
@@ -554,6 +558,7 @@ def build_workflow(opts, retval):
         use_aroma=opts.use_aroma,
         aroma_melodic_dim=opts.aroma_melodic_dimensionality,
         ignore_aroma_err=opts.ignore_aroma_denoising_errors,
+        use_mcflirt=opts.hmc_use_mcflirt,
     )
     retval['return_code'] = 0
 

fmriprep/interfaces/__init__.py

@@ -16,5 +16,5 @@
 from .utils import TPM2ROI, AddTPMs, AddTSVHeader, ConcatAffines, JoinTSVColumns
 from .fmap import FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap
 from .confounds import GatherConfounds, ICAConfounds, FMRISummary
-from .itk import MCFLIRT2ITK, MultiApplyTransforms
+from .itk import MCFLIRT2ITK, Volreg2ITK, MultiApplyTransforms
 from .multiecho import FirstEcho

fmriprep/interfaces/images.py

@@ -21,6 +21,7 @@
     traits, TraitedSpec, BaseInterfaceInputSpec, SimpleInterface,
     File, InputMultiPath, OutputMultiPath)
 from nipype.interfaces import fsl
+from fmriprep.utils.misc import remove_rotation_and_shear
 
 LOGGER = logging.getLogger('nipype.interface')
 
@@ -293,6 +294,7 @@ def _run_interface(self, runtime):
 
 class ValidateImageInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='input image')
+    remove_rotation_and_shear = traits.Bool(False, usedefault=True)
 
 
 class ValidateImageOutputSpec(TraitedSpec):
@@ -342,6 +344,7 @@ class ValidateImage(SimpleInterface):
     output_spec = ValidateImageOutputSpec
 
     def _run_interface(self, runtime):
+
         img = nb.load(self.inputs.in_file)
         out_report = os.path.join(runtime.cwd, 'report.html')
 
@@ -370,7 +373,14 @@ def _run_interface(self, runtime):
 
         # Both match, qform valid (implicit with match), codes okay -> do nothing, empty report
         if matching_affines and qform_code > 0 and sform_code > 0:
-            self._results['out_file'] = self.inputs.in_file
+            if self.inputs.remove_rotation_and_shear:
+                out_fname = fname_presuffix(self.inputs.in_file, suffix='_valid',
+                                            newpath=runtime.cwd)
+                img = remove_rotation_and_shear(img)
+                img.to_filename(out_fname)
+                self._results['out_file'] = out_fname
+            else:
+                self._results['out_file'] = self.inputs.in_file
             open(out_report, 'w').close()
             self._results['out_report'] = out_report
             return runtime
@@ -418,6 +428,9 @@ def _run_interface(self, runtime):
 </p>
 """
         snippet = '<h3 class="elem-title">%s</h3>\n%s\n' % (warning_txt, description)
+
+        if self.inputs.remove_rotation_and_shear:
+            img = remove_rotation_and_shear(img)
         # Store new file and report
         img.to_filename(out_fname)
         with open(out_report, 'w') as fobj:

fmriprep/interfaces/itk.py

@@ -9,6 +9,7 @@
 
 """
 import os
+from pathlib import Path
 from mimetypes import guess_type
 from tempfile import TemporaryDirectory
 import numpy as np
@@ -78,6 +79,49 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class Volreg2ITKInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True,
+                   desc='mat file generated by AFNI\'s 3dVolreg')
+
+
+class Volreg2ITKOutputSpec(TraitedSpec):
+    out_file = File(desc='the output ITKTransform file')
+
+
+class Volreg2ITK(SimpleInterface):
+
+    """
+    Convert an AFNI's mat file into an ITK Transform file.
+    """
+    input_spec = Volreg2ITKInputSpec
+    output_spec = Volreg2ITKOutputSpec
+
+    def _run_interface(self, runtime):
+        # Load AFNI mat entries and reshape appropriately
+        orig_afni_mat = np.loadtxt(self.inputs.in_file)
+        afni_affines = [mat.reshape(3, 4, order='C') for mat in orig_afni_mat]
+
+        out_file = Path(fname_presuffix(self.inputs.in_file, use_ext=False,
+                        suffix='_mc4d_itk.txt', newpath=runtime.cwd))
+
+        fixed_params = 'FixedParameters: 0 0 0'  # Center of rotation does not change
+        lines = ["#Insight Transform File V1.0"]
+        for i, affine in enumerate(afni_affines):
+            lines.append("#Transform %d" % i)
+            lines.append("Transform: AffineTransform_double_3_3")
+
+            ants_affine_2d = np.hstack((affine[:3, :3].reshape(1, -1),
+                                        affine[:3, 3].reshape(1, -1)))
+            params = ants_affine_2d.reshape(-1).astype('float64')
+            params_list = ["%g" % i for i in params.tolist()]
+            lines.append("Parameters: %s" % ' '.join(params_list))
+            lines.append(fixed_params)
+
+        out_file.write_text('\n'.join(lines))
+        self._results['out_file'] = str(out_file)
+        return runtime
+
+
 class MultiApplyTransformsInputSpec(ApplyTransformsInputSpec):
     input_image = InputMultiPath(File(exists=True), mandatory=True,
                                  desc='input time-series as a list of volumes after splitting'

fmriprep/utils/misc.py

@@ -7,6 +7,31 @@
 """
 
 
+def remove_rotation_and_shear(img):
+    from transforms3d.affines import decompose, compose
+    import numpy as np
+    import nibabel as nb
+
+    img = nb.as_closest_canonical(img)
+    T, _, Z, _ = decompose(img.affine)
+    affine = compose(T=T, R=np.diag([1, 1, 1]), Z=Z)
+    return img.__class__(np.asanyarray(img.dataobj), affine, img.header)
+
+
+def split_and_rm_rotshear_func(in_file):
+    import os
+    import nibabel as nb
+    from fmriprep.utils.misc import remove_rotation_and_shear
+    out_files = []
+    imgs = nb.four_to_three(nb.load(in_file))
+    for i, img in enumerate(imgs):
+        out_file = os.path.abspath('vol%04d.nii.gz' % i)
+        img = remove_rotation_and_shear(img)
+        img.to_filename(out_file)
+        out_files.append(out_file)
+    return out_files
+
+
 def fix_multi_T1w_source_name(in_files):
     """
     Make up a generic source name when there are multiple T1s

fmriprep/workflows/base.py

@@ -38,7 +38,8 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid, work_dir, output_dir, bids
                      omp_nthreads, skull_strip_template, skull_strip_fixed_seed,
                      freesurfer, output_spaces, template, medial_surface_nan, cifti_output, hires,
                      use_bbr, bold2t1w_dof, fmap_bspline, fmap_demean, use_syn, force_syn,
-                     use_aroma, ignore_aroma_err, aroma_melodic_dim, template_out_grid):
+                     use_aroma, ignore_aroma_err, aroma_melodic_dim, template_out_grid,
+                     use_mcflirt):
     """
     This workflow organizes the execution of FMRIPREP, with a sub-workflow for
     each subject.
@@ -84,7 +85,8 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid, work_dir, output_dir, bids
                               use_aroma=False,
                               ignore_aroma_err=False,
                               aroma_melodic_dim=None,
-                              template_out_grid='native')
+                              template_out_grid='native',
+                              use_mcflirt=False)
 
 
     Parameters
@@ -162,6 +164,8 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid, work_dir, output_dir, bids
         template_out_grid : str
             Keyword ('native', '1mm' or '2mm') or path of custom reference
             image for normalization
+        use_mcflirt : bool
+            Use FSL ``mcflirt`` for head motion correction, instead of AFNI ``3dVolreg``.
 
     """
     fmriprep_wf = Workflow(name='fmriprep_wf')
@@ -209,6 +213,7 @@ def init_fmriprep_wf(subject_list, task_id, run_uuid, work_dir, output_dir, bids
             use_aroma=use_aroma,
             aroma_melodic_dim=aroma_melodic_dim,
             ignore_aroma_err=ignore_aroma_err,
+            use_mcflirt=use_mcflirt,
         )
 
         single_subject_wf.config['execution']['crashdump_dir'] = (
@@ -231,7 +236,7 @@ def init_single_subject_wf(subject_id, task_id, name, reportlets_dir, output_dir
                            freesurfer, output_spaces, template, medial_surface_nan,
                            cifti_output, hires, use_bbr, bold2t1w_dof, fmap_bspline, fmap_demean,
                            use_syn, force_syn, template_out_grid,
-                           use_aroma, aroma_melodic_dim, ignore_aroma_err):
+                           use_aroma, aroma_melodic_dim, ignore_aroma_err, use_mcflirt):
     """
     This workflow organizes the preprocessing pipeline for a single subject.
     It collects and reports information about the subject, and prepares
@@ -278,7 +283,8 @@ def init_single_subject_wf(subject_id, task_id, name, reportlets_dir, output_dir
                                     template_out_grid='native',
                                     use_aroma=False,
                                     aroma_melodic_dim=None,
-                                    ignore_aroma_err=False)
+                                    ignore_aroma_err=False,
+                                    use_mcflirt=False)
 
     Parameters
 
@@ -355,6 +361,8 @@ def init_single_subject_wf(subject_id, task_id, name, reportlets_dir, output_dir
             Perform ICA-AROMA on MNI-resampled functional series
         ignore_aroma_err : bool
             Do not fail on ICA-AROMA errors
+        use_mcflirt : bool
+            Use FSL ``mcflirt`` for head motion correction, instead of AFNI ``3dVolreg``.
 
     Inputs
 
@@ -493,7 +501,8 @@ def init_single_subject_wf(subject_id, task_id, name, reportlets_dir, output_dir
                                                use_aroma=use_aroma,
                                                aroma_melodic_dim=aroma_melodic_dim,
                                                ignore_aroma_err=ignore_aroma_err,
-                                               num_bold=len(subject_data['bold']))
+                                               num_bold=len(subject_data['bold']),
+                                               use_mcflirt=use_mcflirt)
 
         workflow.connect([
             (anat_preproc_wf, func_preproc_wf,

fmriprep/workflows/bold/base.py

@@ -15,10 +15,10 @@
 import nibabel as nb
 from nipype import logging
 
-from nipype.interfaces.fsl import Split as FSLSplit
 from nipype.pipeline import engine as pe
 from nipype.interfaces import utility as niu
 
+from fmriprep.utils.misc import split_and_rm_rotshear_func
 from ...interfaces import (
     DerivativesDataSink,
     GiftiNameSource,
@@ -53,7 +53,7 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
                          use_aroma, ignore_aroma_err, aroma_melodic_dim,
                          medial_surface_nan, cifti_output,
                          debug, low_mem, template_out_grid,
-                         layout=None, num_bold=1):
+                         use_mcflirt, layout=None, num_bold=1):
     """
     This workflow controls the functional preprocessing stages of FMRIPREP.
 
@@ -86,6 +86,7 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
                                   use_aroma=False,
                                   ignore_aroma_err=False,
                                   aroma_melodic_dim=None,
+                                  use_mcflirt=False,
                                   num_bold=1)
 
     **Parameters**
@@ -146,6 +147,8 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
         template_out_grid : str
             Keyword ('native', '1mm' or '2mm') or path of custom reference
             image for normalization
+        use_mcflirt : bool
+            Use FSL ``mcflirt`` for head motion correction, instead of AFNI ``3dVolreg``.
         layout : BIDSLayout
             BIDSLayout structure to enable metadata retrieval
         num_bold : int
@@ -398,11 +401,13 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
     bold_reference_wf = init_bold_reference_wf(omp_nthreads=omp_nthreads)
 
     # Top-level BOLD splitter
-    bold_split = pe.Node(FSLSplit(dimension='t'), name='bold_split',
+    bold_split = pe.Node(niu.Function(function=split_and_rm_rotshear_func, input_names=['in_file'],
+                                      output_names=['out_files']), name='bold_split',
                          mem_gb=mem_gb['filesize'] * 3)
 
     # HMC on the BOLD
-    bold_hmc_wf = init_bold_hmc_wf(name='bold_hmc_wf',
+    bold_hmc_wf = init_bold_hmc_wf(use_mcflirt=use_mcflirt,
+                                   name='bold_hmc_wf',
                                    mem_gb=mem_gb['filesize'],
                                    omp_nthreads=omp_nthreads)