From d3c1a41bdb9ea9f3586d0229550f062a564ae6ee Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 10 Feb 2025 14:25:02 -0500 Subject: [PATCH 1/6] feat: Add implementation of FSL RMS Deviation (rmsd) --- fmriprep/interfaces/confounds.py | 52 ++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/fmriprep/interfaces/confounds.py b/fmriprep/interfaces/confounds.py index ec4cf8933..dae07939e 100644 --- a/fmriprep/interfaces/confounds.py +++ b/fmriprep/interfaces/confounds.py @@ -90,6 +90,58 @@ def _run_interface(self, runtime): return runtime +class _FSLRMSDeviationInputSpec(BaseInterfaceInputSpec): + xfm_file = File(exists=True, desc='Head motion transform file') + boldref_file = File(exists=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') From 50e4b9226a5350400bde9ab6cf73d0c87dfa6e1a Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 10 Feb 2025 14:36:28 -0500 Subject: [PATCH 2/6] test: Update test timeseries.tsv --- ...blestask_run-01_desc-motion_timeseries.tsv | 34 +++++++++---------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/fmriprep/interfaces/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motion_timeseries.tsv b/fmriprep/interfaces/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motion_timeseries.tsv index f8371b27f..161862a7a 100644 --- a/fmriprep/interfaces/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motion_timeseries.tsv +++ b/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 From e8ac59b816c8915abba95e283a231ae7cddeffae Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 10 Feb 2025 14:38:29 -0500 Subject: [PATCH 3/6] test(rmsd): Verify reconstruction of FSL rmsd --- fmriprep/interfaces/tests/test_confounds.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/fmriprep/interfaces/tests/test_confounds.py b/fmriprep/interfaces/tests/test_confounds.py index 0af8c2e3e..fe14cc0cc 100644 --- a/fmriprep/interfaces/tests/test_confounds.py +++ b/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' From d157171947af8edab9744da246249003ada76a5d Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 10 Feb 2025 14:41:18 -0500 Subject: [PATCH 4/6] fix: Consistently use DataFrame.to_csv(..., na_rep='n/a') --- fmriprep/interfaces/confounds.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/fmriprep/interfaces/confounds.py b/fmriprep/interfaces/confounds.py index dae07939e..cc6152f09 100644 --- a/fmriprep/interfaces/confounds.py +++ b/fmriprep/interfaces/confounds.py @@ -191,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 @@ -224,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 @@ -252,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 @@ -315,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 From c28a6158e80bd054284533ff85f0ba5719606bb4 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Mon, 10 Feb 2025 14:46:35 -0500 Subject: [PATCH 5/6] rf: Calculate RMSD from motion transforms --- fmriprep/workflows/bold/base.py | 1 - fmriprep/workflows/bold/confounds.py | 18 ++++++------------ fmriprep/workflows/bold/fit.py | 11 +---------- fmriprep/workflows/bold/hmc.py | 14 ++------------ 4 files changed, 9 insertions(+), 35 deletions(-) diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index f7c8c4249..0130ebf65 100644 --- a/fmriprep/workflows/bold/base.py +++ b/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'), ]), diff --git a/fmriprep/workflows/bold/confounds.py b/fmriprep/workflows/bold/confounds.py index c21fbb47a..60aa88174 100644 --- a/fmriprep/workflows/bold/confounds.py +++ b/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')]), diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index c7bcbcbab..aa21ca72f 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/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: diff --git a/fmriprep/workflows/bold/hmc.py b/fmriprep/workflows/bold/hmc.py index 6e6b3448c..e85dacae5 100644 --- a/fmriprep/workflows/bold/hmc.py +++ b/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 From d1ebb0084050a79cbd557eee8ce506aa4aeda5d1 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 12 Feb 2025 13:35:28 -0500 Subject: [PATCH 6/6] Update fmriprep/interfaces/confounds.py Co-authored-by: Mathias Goncalves --- fmriprep/interfaces/confounds.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fmriprep/interfaces/confounds.py b/fmriprep/interfaces/confounds.py index cc6152f09..7be150bd9 100644 --- a/fmriprep/interfaces/confounds.py +++ b/fmriprep/interfaces/confounds.py @@ -91,8 +91,8 @@ def _run_interface(self, runtime): class _FSLRMSDeviationInputSpec(BaseInterfaceInputSpec): - xfm_file = File(exists=True, desc='Head motion transform file') - boldref_file = File(exists=True, desc='BOLD reference file') + 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):