Merge pull request #1283 from chrisfilo/enh/3dvolreg

[ENH] Replace FSL ``mcflirt`` with AFNI ``3dVolReg``

Author: oesteban

.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
@@ -778,7 +778,7 @@ workflows:
             - get_regression_data
           filters:
             branches:
-              ignore: 
+              ignore:
                 - /docs?\/.*/
                 - /ds005\/.*/
                 - /ds054\/.*/

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
@@ -257,7 +257,7 @@ Other relevant references
   .. [Power2017] Power JD, Plitt M, Kundu P, Bandettini PA, Martin A (2017) Temporal interpolation alters
       motion in fMRI scans: Magnitudes and consequences for artifact detection. PLOS ONE 12(9): e0182939.
       doi:`10.1371/journal.pone.0182939 <https://doi.org/10.1371/journal.pone.0182939>`_.
-      
-  .. [Brett2001] Brett M, Leff AP, Rorden C, Ashburner J (2001) Spatial Normalization of Brain Images with 
-      Focal Lesions Using Cost Function Masking. NeuroImage 14(2) 
+
+  .. [Brett2001] Brett M, Leff AP, Rorden C, Ashburner J (2001) Spatial Normalization of Brain Images with
+      Focal Lesions Using Cost Function Masking. NeuroImage 14(2)
       doi:`10.006/nimg.2001.0845 <https://doi.org/10.1006/nimg.2001.0845>`_.

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),
@@ -160,7 +160,7 @@ derivative of RMS variance over voxels (or :abbr:`DVARS (D referring to differen
 standardized (``stdDVARS``), non-standardized (``non-stdDVARS``), and voxel-wise standardized (``vx-wisestdDVARS``);
 the ``FrameDisplacement`` is a quantification of the estimated bulk-head motion; ``X``, ``Y``, ``Z``, ``RotX``,
 ``RotY``, ``RotZ`` are the actual 6 rigid-body transform parameters estimated by FMRIPREP;
-the ``NonSteadyStateOutlierXX`` columns indicate non-steady state volumes with a single ``1`` value and ``0`` elsewhere (there 
+the ``NonSteadyStateOutlierXX`` columns indicate non-steady state volumes with a single ``1`` value and ``0`` elsewhere (there
 is one ``NonSteadyStateOutlierXX`` column per outlier/volume); and finally six noise components ``aCompCorXX`` calculated using
 :abbr:`CompCor (Component Based Noise Correction Method)`
 and five noise components ``AROMAaggrCompXX`` if

docs/workflows.rst

@@ -19,38 +19,40 @@ is presented below:
     :simple_form: yes
 
     from fmriprep.workflows.base import init_single_subject_wf
-    wf = init_single_subject_wf(subject_id='test',
-                                name='single_subject_wf',
-                                task_id='',
-                                longitudinal=False,
-                                t2s_coreg=False,
-                                omp_nthreads=1,
-                                freesurfer=True,
-                                reportlets_dir='.',
-                                output_dir='.',
-                                bids_dir='.',
-                                skull_strip_template='OASIS',
-                                skull_strip_fixed_seed=False,
-                                template='MNI152NLin2009cAsym',
-                                output_spaces=['T1w', 'fsnative',
-                                              'template', 'fsaverage5'],
-                                medial_surface_nan=False,
-                                cifti_output=False,
-                                ignore=[],
-                                debug=False,
-                                low_mem=False,
-                                anat_only=False,
-                                hires=True,
-                                use_bbr=True,
-                                bold2t1w_dof=9,
-                                fmap_bspline=False,
-                                fmap_demean=True,
-                                use_syn=True,
-                                force_syn=True,
-                                template_out_grid='native',
-                                use_aroma=False,
-                                aroma_melodic_dim=None,
-                                ignore_aroma_err=False)
+    wf = init_single_subject_wf(
+        subject_id='test',
+        name='single_subject_wf',
+        task_id='',
+        longitudinal=False,
+        t2s_coreg=False,
+        omp_nthreads=1,
+        freesurfer=True,
+        reportlets_dir='.',
+        output_dir='.',
+        bids_dir='.',
+        skull_strip_template='OASIS',
+        skull_strip_fixed_seed=False,
+        template='MNI152NLin2009cAsym',
+        output_spaces=['T1w', 'fsnative', 'template', 'fsaverage5'],
+        medial_surface_nan=False,
+        cifti_output=False,
+        ignore=[],
+        debug=False,
+        low_mem=False,
+        anat_only=False,
+        hires=True,
+        use_bbr=True,
+        bold2t1w_dof=9,
+        fmap_bspline=False,
+        fmap_demean=True,
+        use_syn=True,
+        force_syn=True,
+        template_out_grid='native',
+        use_aroma=False,
+        aroma_melodic_dim=None,
+        ignore_aroma_err=False,
+        use_mcflirt=False,
+    )
 
 
 T1w/T2w preprocessing
@@ -248,30 +250,32 @@ BOLD preprocessing
     :simple_form: yes
 
     from fmriprep.workflows.bold import init_func_preproc_wf
-    wf = init_func_preproc_wf('/completely/made/up/path/sub-01_task-nback_bold.nii.gz',
-                              omp_nthreads=1,
-                              ignore=[],
-                              freesurfer=True,
-                              reportlets_dir='.',
-                              output_dir='.',
-                              template='MNI152NLin2009cAsym',
-                              output_spaces=['T1w', 'fsnative',
-                                             'template', 'fsaverage5'],
-                              medial_surface_nan=False,
-                              cifti_output=False,
-                              debug=False,
-                              low_mem=False,
-                              use_bbr=True,
-                              t2s_coreg=False,
-                              bold2t1w_dof=9,
-                              fmap_bspline=True,
-                              fmap_demean=True,
-                              use_syn=True,
-                              force_syn=True,
-                              template_out_grid='native',
-                              use_aroma=False,
-                              aroma_melodic_dim=None,
-                              ignore_aroma_err=False)
+    wf = init_func_preproc_wf(
+        '/completely/made/up/path/sub-01_task-nback_bold.nii.gz',
+        omp_nthreads=1,
+        ignore=[],
+        freesurfer=True,
+        reportlets_dir='.',
+        output_dir='.',
+        template='MNI152NLin2009cAsym',
+        output_spaces=['T1w', 'fsnative', 'template', 'fsaverage5'],
+        medial_surface_nan=False,
+        cifti_output=False,
+        debug=False,
+        low_mem=False,
+        use_bbr=True,
+        t2s_coreg=False,
+        bold2t1w_dof=9,
+        fmap_bspline=True,
+        fmap_demean=True,
+        use_syn=True,
+        force_syn=True,
+        template_out_grid='native',
+        use_aroma=False,
+        aroma_melodic_dim=None,
+        ignore_aroma_err=False,
+        use_mcflirt=False,
+    )
 
 Preprocessing of :abbr:`BOLD (blood-oxygen level-dependent)` files is
 split into multiple sub-workflows described below.
@@ -332,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.
@@ -546,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,30 @@
 """
 
 
+def remove_rotation_and_shear(img):
+    from transforms3d.affines import decompose, compose
+    import numpy as np
+
+    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(
+            nb.as_closest_canonical(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