Merge pull request #2088 from rciric/fsl_motion_outliers

ENH: Add RMS framewise displacement following Jenkinson et al. (2002)

Author: oesteban

docs/outputs.rst

@@ -275,6 +275,8 @@ frames with sudden and large motion or intensity spikes.
 
 - ``framewise_displacement`` - is a quantification of the estimated bulk-head motion calculated using
   formula proposed by [Power2012]_;
+- ``rmsd`` - is a quantification of the estimated relative (frame-to-frame) bulk head motion
+  calculated using the :abbr:`RMS (root mean square)` approach of [Jenkinson2002]_;
 - ``dvars`` - the derivative of RMS variance over voxels (or :abbr:`DVARS (derivative of
   RMS variance over voxels)`) [Power2012]_;
 - ``std_dvars`` - standardized :abbr:`DVARS (derivative of RMS variance over voxels)`;
@@ -529,6 +531,10 @@ See implementation on :mod:`~fmriprep.workflows.bold.confounds.init_bold_confs_w
      NeuroImage. 2013. doi:`10.1016/j.neuroimage.2013.05.116
      <https://doi.org/10.1016/j.neuroimage.2013.05.116>`_
 
+  .. [Jenkinson2002] Jenkinson M, Bannister P, Brady M, Smith S. Improved optimization for the
+     robust and accurate linear registration and motion correction of brain images. Neuroimage.
+     2002. doi:`10.1016/s1053-8119(02)91132-8 <https://doi.org/10.1016/s1053-8119(02)91132-8>`__.
+
   .. [Lindquist2019] Lindquist, MA, Geuter, S, and Wager, TD, Caffo, BS,
      Modular preprocessing pipelines can reintroduce artifacts into fMRI data.
      Human Brain Mapping. 2019. doi: `10.1002/hbm.24528

fmriprep/interfaces/confounds.py

@@ -29,6 +29,7 @@ class GatherConfoundsInputSpec(BaseInterfaceInputSpec):
     dvars = File(exists=True, desc='file containing DVARS')
     std_dvars = File(exists=True, desc='file containing standardized DVARS')
     fd = File(exists=True, desc='input framewise displacement')
+    rmsd = File(exists=True, desc='input RMS framewise displacement')
     tcompcor = File(exists=True, desc='input tCompCorr')
     acompcor = File(exists=True, desc='input aCompCorr')
     cos_basis = File(exists=True, desc='input cosine basis')
@@ -82,6 +83,7 @@ def _run_interface(self, runtime):
             dvars=self.inputs.dvars,
             std_dvars=self.inputs.std_dvars,
             fdisp=self.inputs.fd,
+            rmsd=self.inputs.rmsd,
             tcompcor=self.inputs.tcompcor,
             acompcor=self.inputs.acompcor,
             cos_basis=self.inputs.cos_basis,
@@ -135,7 +137,7 @@ def _run_interface(self, runtime):
 
 
 def _gather_confounds(signals=None, dvars=None, std_dvars=None, fdisp=None,
-                      tcompcor=None, acompcor=None, cos_basis=None,
+                      rmsd=None, tcompcor=None, acompcor=None, cos_basis=None,
                       motion=None, aroma=None, newpath=None):
     r"""
     Load confounds from the filenames, concatenate together horizontally
@@ -186,6 +188,7 @@ def _adjust_indices(left_df, right_df):
                            (std_dvars, 'Standardized DVARS'),
                            (dvars, 'DVARS'),
                            (fdisp, 'Framewise displacement'),
+                           (rmsd, 'Framewise displacement (RMS)'),
                            (tcompcor, 'tCompCor'),
                            (acompcor, 'aCompCor'),
                            (cos_basis, 'Cosine basis'),

fmriprep/workflows/bold/base.py

@@ -473,7 +473,8 @@ def init_func_preproc_wf(bold_file):
         (inputnode, bold_confounds_wf, [('t1w_tpms', 'inputnode.t1w_tpms'),
                                         ('t1w_mask', 'inputnode.t1w_mask')]),
         (bold_hmc_wf, bold_confounds_wf, [
-            ('outputnode.movpar_file', 'inputnode.movpar_file')]),
+            ('outputnode.movpar_file', 'inputnode.movpar_file'),
+            ('outputnode.rmsd_file', 'inputnode.rmsd_file')]),
         (bold_reg_wf, bold_confounds_wf, [
             ('outputnode.itk_t1_to_bold', 'inputnode.t1_bold_xform')]),
         (bold_reference_wf, bold_confounds_wf, [

fmriprep/workflows/bold/confounds.py

@@ -98,6 +98,8 @@ def init_bold_confs_wf(
         BOLD series mask
     movpar_file
         SPM-formatted motion parameters file
+    rmsd_file
+        Framewise displacement as measured by ``fsl_motion_outliers``.
     skip_vols
         number of non steady state volumes
     t1w_mask
@@ -140,6 +142,9 @@ def init_bold_confs_wf(
 Several confounding time-series were calculated based on the
 *preprocessed BOLD*: framewise displacement (FD), DVARS and
 three region-wise global signals.
+FD was computed using two formulations following Power (absolute sum of
+relative motions, @power_fd_dvars) and Jenkinson (relative root mean square
+displacement between affines, @mcflirt).
 FD and DVARS are calculated for each functional run, both using their
 implementations in *Nipype* [following the definitions by @power_fd_dvars].
 The three global signals are extracted within the CSF, the WM, and
@@ -173,8 +178,8 @@ def init_bold_confs_wf(
 were annotated as motion outliers.
 """.format(fd=regressors_fd_th, dv=regressors_dvars_th)
     inputnode = pe.Node(niu.IdentityInterface(
-        fields=['bold', 'bold_mask', 'movpar_file', 'skip_vols',
-                't1w_mask', 't1w_tpms', 't1_bold_xform']),
+        fields=['bold', 'bold_mask', 'movpar_file', 'rmsd_file',
+                'skip_vols', 't1w_mask', 't1w_tpms', 't1_bold_xform']),
         name='inputnode')
     outputnode = pe.Node(niu.IdentityInterface(
         fields=['confounds_file', 'confounds_metadata']),
@@ -258,6 +263,9 @@ def init_bold_confs_wf(
     add_motion_headers = pe.Node(
         AddTSVHeader(columns=["trans_x", "trans_y", "trans_z", "rot_x", "rot_y", "rot_z"]),
         name="add_motion_headers", mem_gb=0.01, run_without_submitting=True)
+    add_rmsd_header = pe.Node(
+        AddTSVHeader(columns=["rmsd"]),
+        name="add_rmsd_header", mem_gb=0.01, run_without_submitting=True)
     concat = pe.Node(GatherConfounds(), name="concat", mem_gb=0.01, run_without_submitting=True)
 
     # CompCor metadata
@@ -382,6 +390,7 @@ def _pick_wm(files):
 
         # Collate computed confounds together
         (inputnode, add_motion_headers, [('movpar_file', 'in_file')]),
+        (inputnode, add_rmsd_header, [('rmsd_file', 'in_file')]),
         (dvars, add_dvars_header, [('out_nstd', 'in_file')]),
         (dvars, add_std_dvars_header, [('out_std', 'in_file')]),
         (signals, concat, [('out_file', 'signals')]),
@@ -390,6 +399,7 @@ def _pick_wm(files):
                             ('pre_filter_file', 'cos_basis')]),
         (acompcor, concat, [('components_file', 'acompcor')]),
         (add_motion_headers, concat, [('out_file', 'motion')]),
+        (add_rmsd_header, concat, [('out_file', 'rmsd')]),
         (add_dvars_header, concat, [('out_file', 'dvars')]),
         (add_std_dvars_header, concat, [('out_file', 'std_dvars')]),
 

fmriprep/workflows/bold/hmc.py

@@ -54,6 +54,8 @@ def init_bold_hmc_wf(mem_gb, omp_nthreads, name='bold_hmc_wf'):
         ITKTransform file aligning each volume to ``ref_image``
     movpar_file
         MCFLIRT motion parameters, normalized to SPM format (X, Y, Z, Rx, Ry, Rz)
+    rms_file
+        Framewise displacement as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
 
     """
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
@@ -68,15 +70,18 @@ def init_bold_hmc_wf(mem_gb, omp_nthreads, name='bold_hmc_wf'):
 `mcflirt` [FSL {fsl_ver}, @mcflirt].
 """.format(fsl_ver=fsl.Info().version() or '<ver>')
 
-    inputnode = pe.Node(niu.IdentityInterface(fields=['bold_file', 'raw_ref_image']),
-                        name='inputnode')
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=['bold_file', 'raw_ref_image']),
+        name='inputnode')
     outputnode = pe.Node(
-        niu.IdentityInterface(fields=['xforms', 'movpar_file']),
+        niu.IdentityInterface(
+            fields=['xforms', 'movpar_file', 'rmsd_file']),
         name='outputnode')
 
     # Head motion correction (hmc)
-    mcflirt = pe.Node(fsl.MCFLIRT(save_mats=True, save_plots=True),
-                      name='mcflirt', mem_gb=mem_gb * 3)
+    mcflirt = pe.Node(
+        fsl.MCFLIRT(save_mats=True, save_plots=True, save_rms=True),
+        name='mcflirt', mem_gb=mem_gb * 3)
 
     fsl2itk = pe.Node(MCFLIRT2ITK(), name='fsl2itk',
                       mem_gb=0.05, n_procs=omp_nthreads)
@@ -85,13 +90,17 @@ def init_bold_hmc_wf(mem_gb, omp_nthreads, name='bold_hmc_wf'):
                                name="normalize_motion",
                                mem_gb=DEFAULT_MEMORY_MIN_GB)
 
+    def _pick_rel(rms_files):
+        return rms_files[-1]
+
     workflow.connect([
         (inputnode, mcflirt, [('raw_ref_image', 'ref_file'),
                               ('bold_file', 'in_file')]),
         (inputnode, fsl2itk, [('raw_ref_image', 'in_source'),
                               ('raw_ref_image', 'in_reference')]),
         (mcflirt, fsl2itk, [('mat_file', 'in_files')]),
         (mcflirt, normalize_motion, [('par_file', 'in_file')]),
+        (mcflirt, outputnode, [(('rms_files', _pick_rel), 'rmsd_file')]),
         (fsl2itk, outputnode, [('out_file', 'xforms')]),
         (normalize_motion, outputnode, [('out_file', 'movpar_file')]),
     ])