rf: Calculate RMSD from motion transforms (#3427)

Reconstructed from:

*
https://git.fmrib.ox.ac.uk/fsl/mcflirt/-/blob/master/mcflirt.cc#L450-452
*
https://git.fmrib.ox.ac.uk/fsl/miscmaths/-/blob/master/miscmaths.cc#L1065-1092

Author: effigies

fmriprep/interfaces/confounds.py

@@ -90,6 +90,58 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _FSLRMSDeviationInputSpec(BaseInterfaceInputSpec):
+    xfm_file = File(exists=True, mandatory=True, desc='Head motion transform file')
+    boldref_file = File(exists=True, mandatory=True, desc='BOLD reference file')
+
+
+class _FSLRMSDeviationOutputSpec(TraitedSpec):
+    out_file = File(desc='Output motion parameters file')
+
+
+class FSLRMSDeviation(SimpleInterface):
+    """Reconstruct FSL root mean square deviation from affine transforms."""
+
+    input_spec = _FSLRMSDeviationInputSpec
+    output_spec = _FSLRMSDeviationOutputSpec
+
+    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)
+
+        center = 0.5 * (np.array(boldref.shape[:3]) - 1) * boldref.header.get_zooms()[:3]
+
+        # 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])
+
+        diff = fsl_matrix[1:] @ np.linalg.inv(fsl_matrix[:-1]) - np.eye(4)
+        M = diff[:, :3, :3]
+        t = diff[:, :3, 3] + M @ center
+        Rmax = 80.0
+
+        rmsd = np.concatenate(
+            [
+                [np.nan],
+                np.sqrt(
+                    np.diag(t @ t.T)
+                    + np.trace(M.transpose(0, 2, 1) @ M, axis1=1, axis2=2) * Rmax**2 / 5
+                ),
+            ]
+        )
+
+        params = pd.DataFrame(data=rmsd, columns=['rmsd'])
+        params.to_csv(self._results['out_file'], sep='\t', index=False, na_rep='n/a')
+
+        return runtime
+
+
 class _FSLMotionParamsInputSpec(BaseInterfaceInputSpec):
     xfm_file = File(exists=True, desc='Head motion transform file')
     boldref_file = File(exists=True, desc='BOLD reference file')
@@ -139,7 +191,7 @@ def _run_interface(self, runtime):
             columns=['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z'],
         )
 
-        params.to_csv(self._results['out_file'], sep='\t', index=False)
+        params.to_csv(self._results['out_file'], sep='\t', index=False, na_rep='n/a')
 
         return runtime
 
@@ -172,7 +224,7 @@ def _run_interface(self, runtime):
 
         fd = pd.DataFrame(diff.abs().sum(axis=1, skipna=False), columns=['FramewiseDisplacement'])
 
-        fd.to_csv(self._results['out_file'], sep='\t', index=False)
+        fd.to_csv(self._results['out_file'], sep='\t', index=False, na_rep='n/a')
 
         return runtime
 
@@ -200,7 +252,9 @@ def _run_interface(self, runtime):
         )
 
         metadata = pd.read_csv(self.inputs.in_file, sep='\t')
-        metadata[metadata.retained].to_csv(self._results['out_file'], sep='\t', index=False)
+        metadata[metadata.retained].to_csv(
+            self._results['out_file'], sep='\t', index=False, na_rep='n/a'
+        )
 
         return runtime
 
@@ -263,13 +317,15 @@ def _run_interface(self, runtime):
         final_components = components.rename(columns=dict(zip(c_orig, c_new, strict=False)))
         final_components.rename(columns=dict(zip(w_orig, w_new, strict=False)), inplace=True)
         final_components.rename(columns=dict(zip(a_orig, a_new, strict=False)), inplace=True)
-        final_components.to_csv(self._results['components_file'], sep='\t', index=False)
+        final_components.to_csv(
+            self._results['components_file'], sep='\t', index=False, na_rep='n/a'
+        )
 
         metadata.loc[c_comp_cor.index, 'component'] = c_new
         metadata.loc[w_comp_cor.index, 'component'] = w_new
         metadata.loc[a_comp_cor.index, 'component'] = a_new
 
-        metadata.to_csv(self._results['metadata_file'], sep='\t', index=False)
+        metadata.to_csv(self._results['metadata_file'], sep='\t', index=False, na_rep='n/a')
 
         return runtime
 

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

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

fmriprep/interfaces/tests/test_confounds.py

@@ -33,6 +33,26 @@ def test_FilterDropped(tmp_path, data_dir):
     assert filtered_meta == target_meta
 
 
+def test_FSLRMSDeviation(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'
+    timeseries = data_dir / f'{base}_desc-motion_timeseries.tsv'
+
+    rmsd = pe.Node(
+        confounds.FSLRMSDeviation(xfm_file=str(xfms), boldref_file=str(boldref)),
+        name='rmsd',
+        base_dir=str(tmp_path),
+    )
+    res = rmsd.run()
+
+    orig = pd.read_csv(timeseries, sep='\t')['rmsd']
+    derived = pd.read_csv(res.outputs.out_file, sep='\t')['rmsd']
+
+    # RMSD is nominally in mm, so 0.1um is an acceptable deviation
+    assert np.allclose(orig.values, derived.values, equal_nan=True, atol=1e-4)
+
+
 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'

fmriprep/workflows/bold/base.py

@@ -675,7 +675,6 @@ def init_bold_wf(
             ('outputnode.bold_mask', 'inputnode.bold_mask'),
             ('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

@@ -40,6 +40,7 @@
     FMRISummary,
     FramewiseDisplacement,
     FSLMotionParams,
+    FSLRMSDeviation,
     GatherConfounds,
     RenameACompCor,
 )
@@ -124,8 +125,6 @@ def init_bold_confs_wf(
         BOLD series mask
     motion_xfm
         ITK-formatted head motion transforms
-    rmsd_file
-        Root mean squared deviation as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
     skip_vols
         number of non steady state volumes
     t1w_mask
@@ -225,7 +224,6 @@ def init_bold_confs_wf(
                 'bold_mask',
                 'hmc_boldref',
                 'motion_xfm',
-                'rmsd_file',
                 'skip_vols',
                 't1w_mask',
                 't1w_tpms',
@@ -269,7 +267,8 @@ def init_bold_confs_wf(
     motion_params = pe.Node(FSLMotionParams(), name='motion_params')
 
     # Frame displacement
-    fdisp = pe.Node(FramewiseDisplacement(), name='fdisp', mem_gb=mem_gb)
+    fdisp = pe.Node(FramewiseDisplacement(), name='fdisp')
+    rmsd = pe.Node(FSLRMSDeviation(), name='rmsd')
 
     # Generate aCompCor probseg maps
     acc_masks = pe.Node(aCompCorMasks(is_aseg=freesurfer), name='acc_masks')
@@ -373,12 +372,6 @@ def init_bold_confs_wf(
         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
@@ -524,6 +517,8 @@ def _select_cols(table):
                             ('bold_mask', 'in_mask')]),
         (inputnode, motion_params, [('motion_xfm', 'xfm_file'),
                                     ('hmc_boldref', 'boldref_file')]),
+        (inputnode, rmsd, [('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'),
@@ -567,7 +562,6 @@ def _select_cols(table):
         (merge_rois, signals, [('out', 'label_files')]),
 
         # Collate computed confounds together
-        (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')]),
@@ -577,7 +571,7 @@ def _select_cols(table):
         (rename_acompcor, concat, [('components_file', 'acompcor')]),
         (crowncompcor, concat, [('components_file', 'crowncompcor')]),
         (motion_params, concat, [('out_file', 'motion')]),
-        (add_rmsd_header, concat, [('out_file', 'rmsd')]),
+        (rmsd, concat, [('out_file', 'rmsd')]),
         (add_dvars_header, concat, [('out_file', 'dvars')]),
         (add_std_dvars_header, concat, [('out_file', 'std_dvars')]),
 

fmriprep/workflows/bold/fit.py

@@ -178,8 +178,6 @@ def init_bold_fit_wf(
     boldref2fmap_xfm
         Affine transform mapping from BOLD reference space to the fieldmap
         space, if applicable.
-    rmsd_file
-        Root mean squared deviation as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
     dummy_scans
         The number of dummy scans declared or detected at the beginning of the series.
 
@@ -285,7 +283,6 @@ def init_bold_fit_wf(
                 'motion_xfm',
                 'boldref2anat_xfm',
                 'boldref2fmap_xfm',
-                'rmsd_file',
             ],
         ),
         name='outputnode',
@@ -299,9 +296,7 @@ def init_bold_fit_wf(
         name='hmcref_buffer',
     )
     fmapref_buffer = pe.Node(niu.Function(function=_select_ref), name='fmapref_buffer')
-    hmc_buffer = pe.Node(
-        niu.IdentityInterface(fields=['hmc_xforms', 'rmsd_file']), name='hmc_buffer'
-    )
+    hmc_buffer = pe.Node(niu.IdentityInterface(fields=['hmc_xforms']), name='hmc_buffer')
     fmapreg_buffer = pe.Node(
         niu.IdentityInterface(fields=['boldref2fmap_xfm']), name='fmapreg_buffer'
     )
@@ -354,7 +349,6 @@ def init_bold_fit_wf(
         (fmapreg_buffer, outputnode, [('boldref2fmap_xfm', 'boldref2fmap_xfm')]),
         (hmc_buffer, outputnode, [
             ('hmc_xforms', 'motion_xfm'),
-            ('rmsd_file', 'rmsd_file'),
         ]),
         (inputnode, func_fit_reports_wf, [
             ('bold_file', 'inputnode.source_file'),
@@ -444,9 +438,6 @@ def init_bold_fit_wf(
                 ('bold_file', 'inputnode.bold_file'),
             ]),
             (bold_hmc_wf, ds_hmc_wf, [('outputnode.xforms', 'inputnode.xforms')]),
-            (bold_hmc_wf, hmc_buffer, [
-                ('outputnode.rmsd_file', 'rmsd_file'),
-            ]),
             (ds_hmc_wf, hmc_buffer, [('outputnode.xforms', 'hmc_xforms')]),
         ])  # fmt:skip
     else:

fmriprep/workflows/bold/hmc.py

@@ -71,8 +71,6 @@ def init_bold_hmc_wf(mem_gb: float, omp_nthreads: int, name: str = 'bold_hmc_wf'
     -------
     xforms
         ITKTransform file aligning each volume to ``ref_image``
-    rmsd_file
-        Root mean squared deviation as measured by ``fsl_motion_outliers`` [Jenkinson2002]_.
 
     """
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
@@ -89,27 +87,19 @@ def init_bold_hmc_wf(mem_gb: float, omp_nthreads: int, name: str = 'bold_hmc_wf'
     inputnode = pe.Node(
         niu.IdentityInterface(fields=['bold_file', 'raw_ref_image']), name='inputnode'
     )
-    outputnode = pe.Node(niu.IdentityInterface(fields=['xforms', 'rmsd_file']), name='outputnode')
+    outputnode = pe.Node(niu.IdentityInterface(fields=['xforms']), name='outputnode')
 
     # Head motion correction (hmc)
-    mcflirt = pe.Node(
-        fsl.MCFLIRT(save_mats=True, save_plots=True, save_rms=True),
-        name='mcflirt',
-        mem_gb=mem_gb * 3,
-    )
+    mcflirt = pe.Node(fsl.MCFLIRT(save_mats=True), name='mcflirt', mem_gb=mem_gb * 3)
 
     fsl2itk = pe.Node(MCFLIRT2ITK(), name='fsl2itk', mem_gb=0.05, n_procs=omp_nthreads)
 
-    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, outputnode, [(('rms_files', _pick_rel), 'rmsd_file')]),
         (fsl2itk, outputnode, [('out_file', 'xforms')]),
     ])  # fmt:skip