rf: Reconstruct motion confounds from minimal derivatives (#3424)

After an exploration in nipreps/fmripost-aroma#94, we now know how to reconstruct the original motion parameters (to at least 0.001 mm/rad precision) from the ITK transforms we generate.

To ensure consistency between single-shot and fit+apply runs, I propose to use this function in place of the original outputs.

fmriprep/interfaces/confounds.py

@@ -34,6 +34,7 @@
 import re
 
 import nibabel as nb
+import nitransforms as nt
 import numpy as np
 import pandas as pd
 from nipype import logging
@@ -50,6 +51,8 @@
 from nipype.utils.filemanip import fname_presuffix
 from nireports.reportlets.modality.func import fMRIPlot
 from niworkflows.utils.timeseries import _cifti_timeseries, _nifti_timeseries
+from scipy import ndimage as ndi
+from scipy.spatial import transform as sst
 
 LOGGER = logging.getLogger('nipype.interface')
 
@@ -87,6 +90,93 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _FSLMotionParamsInputSpec(BaseInterfaceInputSpec):
+    xfm_file = File(exists=True, desc='Head motion transform file')
+    boldref_file = File(exists=True, desc='BOLD reference file')
+
+
+class _FSLMotionParamsOutputSpec(TraitedSpec):
+    out_file = File(desc='Output motion parameters file')
+
+
+class FSLMotionParams(SimpleInterface):
+    """Reconstruct FSL motion parameters from affine transforms."""
+
+    input_spec = _FSLMotionParamsInputSpec
+    output_spec = _FSLMotionParamsOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['out_file'] = fname_presuffix(
+            self.inputs.boldref_file, suffix='_motion.tsv', newpath=runtime.cwd
+        )
+
+        boldref = nb.load(self.inputs.boldref_file)
+        hmc = nt.linear.load(self.inputs.xfm_file)
+
+        # FSL's "center of gravity" is the center of mass scaled by zooms
+        # No rotation is applied.
+        center_of_gravity = np.matmul(
+            np.diag(boldref.header.get_zooms()),
+            ndi.center_of_mass(np.asanyarray(boldref.dataobj)),
+        )
+
+        # Revert to vox2vox transforms
+        fsl_hmc = nt.io.fsl.FSLLinearTransformArray.from_ras(
+            hmc.matrix, reference=boldref, moving=boldref
+        )
+        fsl_matrix = np.stack([xfm['parameters'] for xfm in fsl_hmc.xforms])
+
+        # FSL uses left-handed rotation conventions, so transpose
+        mats = fsl_matrix[:, :3, :3].transpose(0, 2, 1)
+
+        # Rotations are recovered directly
+        rot_xyz = sst.Rotation.from_matrix(mats).as_euler('XYZ')
+        # Translations are recovered by applying the rotation to the center of gravity
+        trans_xyz = fsl_matrix[:, :3, 3] - mats @ center_of_gravity + center_of_gravity
+
+        params = pd.DataFrame(
+            data=np.hstack((trans_xyz, rot_xyz)),
+            columns=['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z'],
+        )
+
+        params.to_csv(self._results['out_file'], sep='\t', index=False)
+
+        return runtime
+
+
+class _FramewiseDisplacementInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, desc='Head motion transform file')
+    radius = traits.Float(50, usedefault=True, desc='Radius of the head in mm')
+
+
+class _FramewiseDisplacementOutputSpec(TraitedSpec):
+    out_file = File(desc='Output framewise displacement file')
+
+
+class FramewiseDisplacement(SimpleInterface):
+    """Calculate framewise displacement estimates."""
+
+    input_spec = _FramewiseDisplacementInputSpec
+    output_spec = _FramewiseDisplacementOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['out_file'] = fname_presuffix(
+            self.inputs.in_file, suffix='_fd.tsv', newpath=runtime.cwd
+        )
+
+        motion = pd.read_csv(self.inputs.in_file, delimiter='\t')
+
+        # Filter and ensure we have all parameters
+        diff = motion[['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z']].diff()
+        diff[['rot_x', 'rot_y', 'rot_z']] *= self.inputs.radius
+
+        fd = pd.DataFrame(diff.abs().sum(axis=1, skipna=False), columns=['FramewiseDisplacement'])
+
+        fd.to_csv(self._results['out_file'], sep='\t', index=False)
+
+        return runtime
+
+
 class _FilterDroppedInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, desc='input CompCor metadata')
 

fmriprep/interfaces/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motion_timeseries.tsv

@@ -0,0 +1,17 @@
+trans_x	trans_y	trans_z	rot_x	rot_y	rot_z	framewise_displacement
+1.19425e-05	0.0443863	0.00472985	0.000356176	-0.000617445	0.0	n/a
+-2.57666e-05	0.0463662	0.0623273	-0.000208795	-0.00012937	0.0	0.1122673591
+-2.64055e-05	-0.00438628	0.067513	-2.59508e-05	-0.00012937	0.000173904	0.0737762289
+0.0161645	-0.0226134	0.0630764	-2.59508e-05	-0.000199844	0.000279081	0.0476371754999999
+0.0161497	-0.0263834	0.0464668	0.000161259	-0.00012937	0.000573335	0.04799129
+0.0161482	-0.0226144	0.0345415	6.52323e-05	-7.276439999999998e-05	0.000573335	0.023327415
+0.0121946	-0.00426109	0.0671039	-2.59508e-05	-0.00012937	0.000312581	0.075296445
+0.0121556	-0.0175135	0.04042	-2.59508e-05	-0.00012937	0.000166363	0.04728621
+0.0126526	-0.000813328	0.0778061	-2.59508e-05	-0.00012937	0.000166363	0.054583272
+0.012614	0.0250656	0.106248	-0.000320333	0.000149271	0.000166363	0.0830105879999999
+0.0126599	-0.0252459	0.0731423999999999	2.99512e-05	-0.000204037	0.000166363	0.11864261
+-0.00608005	0.0349207	0.110289	-0.000485119	-0.00012937	8.57718e-05	0.1495695699999999
+-0.00607796	-0.00933714	0.0796319999999999	-0.000126125	-0.00012937	0.000166363	0.09689619
+0.0124531	0.00996903	0.0986678	-0.000266813	-0.000207949	0.000166363	0.06783638
+0.010915	-0.00933714	0.0986667	-0.000126125	-0.000248281	0.000166363	0.02989637
+0.0119349	0.00531637	0.0808524	-0.000209587	-0.0	0.000226933	0.0531033599999999

fmriprep/interfaces/tests/data/sub-01_task-mixedgamblestask_run-01_from-orig_to-boldref_mode-image_desc-hmc_xfm.txt

@@ -0,0 +1,80 @@
+#Insight Transform File V1.0
+#Transform 0
+Transform: AffineTransform_float_3_3
+Parameters: 1 2.19652e-07 -0.000617 2.19652e-07 1 0.000356 0.000617 -0.000356 0.999999 -0.0198273 0.0558417 -0.00925011
+FixedParameters: 0 0 0
+
+#Transform 1
+Transform: AffineTransform_float_3_3
+Parameters: 1 -2.6961e-08 -0.000129 -2.6961e-08 1 -0.000209 0.000129 0.000209 1 -0.00408853 0.0396441 -0.0579032
+FixedParameters: 0 0 0
+
+#Transform 2
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000173997 -0.000129005 -0.000174003 1 -2.59776e-05 0.000128995 2.60224e-05 1 -0.000785439 -0.00589913 -0.0665673
+FixedParameters: 0 0 0
+
+#Transform 3
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000278995 -0.000200007 -0.000279005 1 -2.59442e-05 0.000199993 2.60558e-05 1 -0.0172875 -0.0245076 -0.0618151
+FixedParameters: 0 0 0
+
+#Transform 4
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000573021 -0.000128908 -0.000572979 1 0.000161074 0.000129092 -0.000160926 1 -0.00936206 -0.0233705 -0.0490914
+FixedParameters: 0 0 0
+
+#Transform 5
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000573005 -7.29627e-05 -0.000572995 1 6.50418e-05 7.30372e-05 -6.49581e-05 1 -0.0076044 -0.0226739 -0.0354961
+FixedParameters: 0 0 0
+
+#Transform 6
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000312997 -0.000129008 -0.000313003 1 -2.59596e-05 0.000128992 2.60404e-05 1 -0.0103974 -0.00632903 -0.0661615
+FixedParameters: 0 0 0
+
+#Transform 7
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000165997 -0.000129004 -0.000166003 1 -2.59786e-05 0.000128996 2.60214e-05 1 -0.0130841 -0.0189498 -0.0394762
+FixedParameters: 0 0 0
+
+#Transform 8
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000165997 -0.000129004 -0.000166003 1 -2.59786e-05 0.000128996 2.60214e-05 1 -0.0135735 -0.00224875 -0.0768618
+FixedParameters: 0 0 0
+
+#Transform 9
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000166048 0.000148947 -0.000165952 1 -0.000320025 -0.000149053 0.000319975 1 -0.00466364 0.014244 -0.100674
+FixedParameters: 0 0 0
+
+#Transform 10
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000166006 -0.000203995 -0.000165994 1 3.00339e-05 0.000204005 -2.99661e-05 1 -0.0160123 -0.0248835 -0.0729365
+FixedParameters: 0 0 0
+
+#Transform 11
+Transform: AffineTransform_float_3_3
+Parameters: 1 8.59374e-05 -0.000129042 -8.60625e-05 1 -0.000484989 0.000128958 0.000485011 1 0.00358861 0.0190405 -0.100594
+FixedParameters: 0 0 0
+
+#Transform 12
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000165984 -0.000129021 -0.000166016 1 -0.000125979 0.000128979 0.000126021 1 0.00515558 -0.0139722 -0.0767701
+FixedParameters: 0 0 0
+
+#Transform 13
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000165944 -0.000208044 -0.000166056 1 -0.000266965 0.000207956 0.000267034 1 -0.0159206 0.000794425 -0.0928171
+FixedParameters: 0 0 0
+
+#Transform 14
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000165969 -0.000248021 -0.000166031 1 -0.000125959 0.000247979 0.000126041 1 -0.0156607 -0.0139734 -0.0953506
+FixedParameters: 0 0 0
+
+#Transform 15
+Transform: AffineTransform_float_3_3
+Parameters: 1 0.000227 -4.767e-08 -0.000227 1 -0.00021 -4.767e-08 0.00021 1 -0.00762153 -0.00231912 -0.0768985
+FixedParameters: 0 0 0

fmriprep/interfaces/tests/test_confounds.py

@@ -1,5 +1,7 @@
 from pathlib import Path
 
+import numpy as np
+import pandas as pd
 from nipype.pipeline import engine as pe
 
 from fmriprep.interfaces import confounds
@@ -29,3 +31,53 @@ def test_FilterDropped(tmp_path, data_dir):
     target_meta = Path.read_text(data_dir / 'component_metadata_filtered.tsv')
     filtered_meta = Path(res.outputs.out_file).read_text()
     assert filtered_meta == target_meta
+
+
+def test_FSLMotionParams(tmp_path, data_dir):
+    base = 'sub-01_task-mixedgamblestask_run-01'
+    xfms = data_dir / f'{base}_from-orig_to-boldref_mode-image_desc-hmc_xfm.txt'
+    boldref = data_dir / f'{base}_desc-hmc_boldref.nii.gz'
+    orig_timeseries = data_dir / f'{base}_desc-motion_timeseries.tsv'
+
+    motion = pe.Node(
+        confounds.FSLMotionParams(xfm_file=str(xfms), boldref_file=str(boldref)),
+        name='fsl_motion',
+        base_dir=str(tmp_path),
+    )
+    res = motion.run()
+
+    derived_params = pd.read_csv(res.outputs.out_file, sep='\t')
+    # orig_timeseries includes framewise_displacement
+    orig_params = pd.read_csv(orig_timeseries, sep='\t')[derived_params.columns]
+
+    # Motion parameters are in mm and rad
+    # These are empirically determined bounds, but they seem reasonable
+    # for the units
+    limits = pd.DataFrame(
+        {
+            'trans_x': [1e-4],
+            'trans_y': [1e-4],
+            'trans_z': [1e-4],
+            'rot_x': [1e-6],
+            'rot_y': [1e-6],
+            'rot_z': [1e-6],
+        }
+    )
+    max_diff = (orig_params - derived_params).abs().max()
+    assert np.all(max_diff < limits)
+
+
+def test_FramewiseDisplacement(tmp_path, data_dir):
+    timeseries = data_dir / 'sub-01_task-mixedgamblestask_run-01_desc-motion_timeseries.tsv'
+
+    framewise_displacement = pe.Node(
+        confounds.FramewiseDisplacement(in_file=str(timeseries)),
+        name='framewise_displacement',
+        base_dir=str(tmp_path),
+    )
+    res = framewise_displacement.run()
+
+    orig = pd.read_csv(timeseries, sep='\t')['framewise_displacement']
+    derived = pd.read_csv(res.outputs.out_file, sep='\t')['FramewiseDisplacement']
+
+    assert np.allclose(orig.values, derived.values, equal_nan=True)

fmriprep/workflows/bold/base.py

@@ -673,7 +673,8 @@ def init_bold_wf(
         ]),
         (bold_fit_wf, bold_confounds_wf, [
             ('outputnode.bold_mask', 'inputnode.bold_mask'),
-            ('outputnode.movpar_file', 'inputnode.movpar_file'),
+            ('outputnode.hmc_boldref', 'inputnode.hmc_boldref'),
+            ('outputnode.motion_xfm', 'inputnode.motion_xfm'),
             ('outputnode.rmsd_file', 'inputnode.rmsd_file'),
             ('outputnode.boldref2anat_xfm', 'inputnode.boldref2anat_xfm'),
             ('outputnode.dummy_scans', 'inputnode.skip_vols'),

fmriprep/workflows/bold/confounds.py

@@ -38,6 +38,8 @@
 from ...interfaces.confounds import (
     FilterDropped,
     FMRISummary,
+    FramewiseDisplacement,
+    FSLMotionParams,
     GatherConfounds,
     RenameACompCor,
 )
@@ -120,8 +122,8 @@ def init_bold_confs_wf(
         when available.
     bold_mask
         BOLD series mask
-    movpar_file
-        SPM-formatted motion parameters file
+    motion_xfm
+        ITK-formatted head motion transforms
     rmsd_file
         Root mean squared deviation as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
     skip_vols
@@ -221,7 +223,8 @@ def init_bold_confs_wf(
             fields=[
                 'bold',
                 'bold_mask',
-                'movpar_file',
+                'hmc_boldref',
+                'motion_xfm',
                 'rmsd_file',
                 'skip_vols',
                 't1w_mask',
@@ -262,8 +265,11 @@ def init_bold_confs_wf(
         mem_gb=mem_gb,
     )
 
+    # Motion parameters
+    motion_params = pe.Node(FSLMotionParams(), name='motion_params')
+
     # Frame displacement
-    fdisp = pe.Node(nac.FramewiseDisplacement(parameter_source='SPM'), name='fdisp', mem_gb=mem_gb)
+    fdisp = pe.Node(FramewiseDisplacement(), name='fdisp', mem_gb=mem_gb)
 
     # Generate aCompCor probseg maps
     acc_masks = pe.Node(aCompCorMasks(is_aseg=freesurfer), name='acc_masks')
@@ -367,12 +373,6 @@ def init_bold_confs_wf(
         mem_gb=0.01,
         run_without_submitting=True,
     )
-    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',
@@ -518,12 +518,13 @@ def _select_cols(table):
             if not col.startswith(('a_comp_cor_', 't_comp_cor_', 'std_dvars'))
         ]
 
-    # fmt:off
     workflow.connect([
         # connect inputnode to each non-anatomical confound node
         (inputnode, dvars, [('bold', 'in_file'),
                             ('bold_mask', 'in_mask')]),
-        (inputnode, fdisp, [('movpar_file', 'in_file')]),
+        (inputnode, motion_params, [('motion_xfm', 'xfm_file'),
+                                    ('hmc_boldref', 'boldref_file')]),
+        (motion_params, fdisp, [('out_file', 'in_file')]),
         # Brain mask
         (inputnode, t1w_mask_tfm, [('t1w_mask', 'input_image'),
                                    ('bold_mask', 'reference_image'),
@@ -566,7 +567,6 @@ def _select_cols(table):
         (merge_rois, signals, [('out', 'label_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')]),
@@ -576,7 +576,7 @@ def _select_cols(table):
                             ('pre_filter_file', 'cos_basis')]),
         (rename_acompcor, concat, [('components_file', 'acompcor')]),
         (crowncompcor, concat, [('components_file', 'crowncompcor')]),
-        (add_motion_headers, concat, [('out_file', 'motion')]),
+        (motion_params, 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')]),
@@ -617,8 +617,7 @@ def _select_cols(table):
         (concat, conf_corr_plot, [('confounds_file', 'confounds_file'),
                                   (('confounds_file', _select_cols), 'columns')]),
         (conf_corr_plot, ds_report_conf_corr, [('out_file', 'in_file')]),
-    ])
-    # fmt: on
+    ])  # fmt: skip
 
     return workflow