Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

rf: Calculate RMSD from motion transforms #3427

Merged
merged 6 commits into from
Feb 12, 2025
Merged

rf: Calculate RMSD from motion transforms #3427

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
d3c1a41
feat: Add implementation of FSL RMS Deviation (rmsd)
effigies Feb 10, 2025
50e4b92
test: Update test timeseries.tsv
effigies Feb 10, 2025
e8ac59b
test(rmsd): Verify reconstruction of FSL rmsd
effigies Feb 10, 2025
d157171
fix: Consistently use DataFrame.to_csv(..., na_rep='n/a')
effigies Feb 10, 2025
c28a615
rf: Calculate RMSD from motion transforms
effigies Feb 10, 2025
d1ebb00
Update fmriprep/interfaces/confounds.py
effigies Feb 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 61 additions & 5 deletions fmriprep/interfaces/confounds.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Comment on lines +133 to +134
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

for my understanding - why the difference in transposing t and M?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

t has shape (n, 3) so t @ t.T is (n, n). We then take the diagonal, so we get a vector of shape (n,). If we swapped the order, we'd get (3, 3), and thus a 3-vector.

M has shape (n, 3, 3) so M.transpose(0,2, 1) @ M has shape (n, 3, 3). The trace(..., axis1=1, axis2=2) is the sum of the diagonals of the (3, 3) submatrices, so (n,). The trace(M) == trace(M.T), so we could just as easily do:

Suggested change
np.diag(t @ t.T)
+ np.trace(M.transpose(0, 2, 1) @ M, axis1=1, axis2=2) * Rmax**2 / 5
np.diag(t @ t.T)
+ np.trace(M @ M.transpose(0, 2, 1), axis1=1, axis2=2) * Rmax**2 / 5

It's all just a way of getting sums of squares, and definitely not the most efficient, but it pretty closely follows the FSL code, while using vectorized operations.

),
]
)

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')
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down
34 changes: 17 additions & 17 deletions ...prep/interfaces/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motion_timeseries.tsv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
20 changes: 20 additions & 0 deletions fmriprep/interfaces/tests/test_confounds.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand Down
1 change: 0 additions & 1 deletion fmriprep/workflows/bold/base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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'),
]),
Expand Down
18 changes: 6 additions & 12 deletions fmriprep/workflows/bold/confounds.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@
FMRISummary,
FramewiseDisplacement,
FSLMotionParams,
FSLRMSDeviation,
GatherConfounds,
RenameACompCor,
)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -225,7 +224,6 @@ def init_bold_confs_wf(
'bold_mask',
'hmc_boldref',
'motion_xfm',
'rmsd_file',
'skip_vols',
't1w_mask',
't1w_tpms',
Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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'),
Expand Down Expand Up @@ -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')]),
Expand All @@ -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')]),

Expand Down
11 changes: 1 addition & 10 deletions fmriprep/workflows/bold/fit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -285,7 +283,6 @@ def init_bold_fit_wf(
'motion_xfm',
'boldref2anat_xfm',
'boldref2fmap_xfm',
'rmsd_file',
],
),
name='outputnode',
Expand All @@ -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'
)
Expand Down Expand Up @@ -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'),
Expand Down Expand Up @@ -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:
Expand Down
14 changes: 2 additions & 12 deletions fmriprep/workflows/bold/hmc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down
Loading