-
Notifications
You must be signed in to change notification settings - Fork 300
rf: Calculate RMSD from motion transforms #3427
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.
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
Conversation
05e25c8
to
c28a615
Compare
Codecov ReportAll modified and coverable lines are covered by tests ✅
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. |
@oesteban What do you think about this? |
Were we not generated rmsd? (looking at 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 :) |
Sorry, I realized that I said more live to Mathias than I put into this PR. The basic problem is when passing:
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. |
That's very cool. I can try to have a deeper look into the code tomorrow, if you can wait. |
np.diag(t @ t.T) | ||
+ np.trace(M.transpose(0, 2, 1) @ M, axis1=1, axis2=2) * Rmax**2 / 5 |
There was a problem hiding this comment.
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
?
There was a problem hiding this comment.
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:
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.
Co-authored-by: Mathias Goncalves <goncalves.mathias@gmail.com>
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. |
This builds on #3427 to add a couple final things to the transform-only workflow: ``` fmriprep $BIDS $OUT participant --derivatives fmriprep=$MINIMAL ``` BOLD masks were not getting generated from the boldref and dummy scans were not being estimated. This PR factors out the validation and dummy scan detection from `init_raw_boldref_wf` and uses the full workflow in fit mode and the subworkflow in transform mode. Will rebase after #3427 is merged. Can compare now with: c28a615...6c78e54
25.0.0 (March 25, 2025) New feature release in the 25.0.x series. This release substantially improves support for pre-computed derivatives. Previous releases would miss some derivatives and rerun the computations. Note that derivatives from previous versions will be accepted, so it should not be necessary to recompute derivatives from previous versions. The recommended command line is:: fmriprep BIDS_DIR OUT_DIR participant --derivatives fmriprep=PRECOMP_DIR Note that multiple derivatives can be specified, for example:: fmriprep BIDS_DIR OUT_DIR participant \ anat=PRECOMPUTED_ANATOMICAL_DIR \ func=PRECOMPUTED_FUNCTIONAL_DIR When the same file is found in multiple derivatives, the last one found takes precedence. Additionally, `--force-*` flags have been consolidated into a single `--force` flag that can take multiple, space-separated arguments. Structural processing changes ----------------------------- We now output white, pial and midthickness fsLR meshes on the subject surface. Look for `sub-<subject>_hemi-<L|R>_space-fsLR_*_<surf>.surf.gii` files. Brain extraction has been modified slightly to more closely match the `antsBrainExtraction.sh` workflow distributed by ANTs. The impact should be minimal, but in rare cases this fixes a crash. Fieldmap processing changes --------------------------- SyN-SDC fieldmap filtering is now single-level, following the improvements for gradient-echo fieldmaps in 24.1. Jacobian-weighting during fieldmap unwarping is now on by default *only* for PEPolar fieldmaps. To enable for other fieldmap types, use `--force fieldmap-jacobian`. All merged pull requests ------------------------ * FIX: Detect and apply precomputed fieldmaps (#3439) * FIX: Calculate bold mask and dummy scans in transform-only runs (#3428) * FIX: Use consistent skull-stripping pre- and post- SDC (#3415) * FIX: Use removeprefix instead of lstrip or ternary operator (#3409) * FIX: Listify sessions when generating reports (#3408) * FIX: Ensure fieldmap is resampled correctly in report (#3387) * FIX: Stop excluding FS minc_modify_header used during fallback registration (#3372) * FIX: Repair and test query for precalculated baseline/boldref files (#3370) * FIX: Repair search for precomputed transforms (#3369) * ENH: Enable Jacobians only for PEPOLAR by default, allow forcing (#3443) * ENH: Create `--force` flag that accepts a list, replacing individual `--force-*` flags (#3442) * ENH: Output fsLR meshes on subject surfaces (#3411) * ENH: Flexibilize "sophisticated" pepolar to allow monomodal execution (#3393) * ENH: Update FSL packages for reported bug fixes (#3374) * RF: Calculate RMSD from motion transforms (#3427) * RF: Reconstruct motion confounds from minimal derivatives (#3424) * RF: Replace deprecated pkgutil.find_loader (#3384) * RF: Upgrade nitransforms and remove workarounds (#3378) * DOC: Fix xfm extension in the outputs docs (#3435) * DOC: Mention fMRIPost-AROMA in parser documentation (#3356) * MNT: Remove CLI flags with expired deprecation periods (#3445) * MNT: Update pinned environment (#3440) * MNT: Bump pins, update RTD config (#3425) * MNT: Declare linux/amd64 platform during Docker build (#3422) * MNT: Bump astral-sh/setup-uv from 4 to 5 (#3417) * MNT: Test support for Python 3.13 (#3416) * MNT: Install Workbench CLI via conda (#3410) * MNT: Update minimum dependencies, test with tox-uv (#3412) * MNT: Install c3d through conda (#3382) * CI: Fetch tags and 200 commits to support describe (#3381) * CI: Build docker images in GHA, store cache inline and push to GHCR (#3380)
Reconstructed from: