rf: Calculate RMSD from motion transforms

[effigies]

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

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 72.92%. Comparing base (f53185e) to head (d1ebb00).
Report is 7 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #3427      +/-   ##
==========================================
+ Coverage   72.70%   72.92%   +0.21%     
==========================================
  Files          57       57              
  Lines        4327     4358      +31     
  Branches      546      547       +1     
==========================================
+ Hits         3146     3178      +32     
+ Misses       1065     1064       -1     
  Partials      116      116              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[effigies]

Comparing the ds005 output from this (right) and master (left):

[effigies]

@oesteban What do you think about this?

[oesteban]

Were we not generated rmsd? (looking at fmriprep/interfaces/tests/data/sub-01_task-mixedgamblestask_run-01_desc-motion_timeseries.tsv)

It does look to replicate both our and their work well, which is great. Do I understand correctly if I say this increases our dependency on MCFLIRT instead of reducing it? Perhaps we should report both (our traditional rmsd and your fsl-like calculation)?

Also, cross-comparison with https://github.com/nipreps/nifreeze/blob/237a9cb0c4a3ca9681ac6adda3f3330337f5c4b7/src/nifreeze/registration/utils.py#L82-L114 would be very interesting :)

[effigies]

Sorry, I realized that I said more live to Mathias than I put into this PR. The basic problem is when passing:

fmriprep $BIDS $OUT participant --derivatives fmriprep=$MINIMAL

The HMC pipeline was generating ITK transforms, motion parameters, and RMSD. If we were getting pre-computed transforms, then we were not running HMC and so not getting motion parameters or RMSD. #3424 dropped the motion parameters from the HMC pipeline and reconstructed them in the confounds pipeline, and this aims to do the same thing for RMSD.

I don't think this increases our dependence on FSL; in fact it means we can generate FSL-equivalent motion parameters / RMSD, given rigid transform matrices from any HMC estimation tool.

[oesteban]

That's very cool. I can try to have a deeper look into the code tomorrow, if you can wait.

[mgxd]

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

[effigies]

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.

[effigies]

Going to merge. @oesteban feel free to do a post-review merge, and we can back it out if you aren't comfortable with it.