ENH: A 3D tensor B-Spline approximator and extrapolator

[oesteban]

This PR finally adds an implementation for B-Spline smoothing and extrapolation of fieldmaps. It also makes the outputs of all estimation methods more consistent, as the outputs of TOPUP are two - the field and the B-Spline coefficients supporting the field.

Improving over TOPUP, this B-Spline modeling allows for multi-level B-Splines. Typically, you will want to use two grids of B-Splines. One with few control points to capture the global, low-spatial-frequency bias that makes the outline of the brain align better with the anatomy after correction, and with more dense coverage of the FoV to capture the more local high-frequency warping typically found at temporal lobes and the ventromedial prefrontal cortex.

This base implementation does not include a B-Spline interpolator interface that is necessary if you move the control points to the target EPI you want to undistort (which is the better way of moving the field, rather than the dense field representation). To make this possible, we will need to "un-project" the coefficients field and bring them into alignment with the original fieldmap when these are not "plumb" (i.e., scanner coordinate system is not aligned with the data array axes). (EDITED 11/25 - it does include the interpolator now)

References: #71, #22.
Resolves: #72.
Resolves: #14.

---

EDITED - Text below was included on 11/25/2020, when the PR was nearing a conclusion.

CHANGES

• Add the BSplineApprox interface: this is the nuts-and-bolts of this PR. This interface implements an approximation method to fit a grid of tensor-product cubic B-Spline basis, regularized with ridge regression. The output of this interface are coefficients maps (one for each resolution level of the B-Spline basis) embedded in NIfTI files. These NIfTI files have an affine aligned with the original fieldmap, but with the voxel sizes of the selected B-Spline basis.
• Add the Coefficients2Warp interface, which takes in coefficients NIfTI files (like those generated by BSplineApprox, and when a conversion interface is set in place, those generated by TOPUP) and interpolates/extrapolates the field. It also converts from Hz to mm. To do so, the full affine of the target EPI being corrected is applied - meaning that it will perform correction along the PE axis considering that the affine can be "oblique" (see https://github.com/nipreps/sdcflows/pull/119/files#diff-c4cb2b592bcb50f08366cd3726e2ea4ff006fe1e43db8ff4b2e20cc6dd14161dR282). In other words, it does not naively multiply by the pixel spacing assuming that the image is "plumb" (using AFNI's nomenclature).
  BUGFIX - While working on this interface, an important bug calculating the Total Readout Time from the EffectiveEchoSpacing metadata has been solved - if the echo spacing info is "effective" it means that parallel acquisition and other acceleration factors have been already accounted for in the given figure (see https://github.com/nipreps/sdcflows/pull/119/files#diff-f9d322bceffe55f3a5c8b80d0c03c1a923ee838531b19ca12ad2d9435ceb4df9R107-R115).
• [Workflows] Adapted init_fmap_wf so that now, all types of fieldmaps (phasediff, phase1/2, and fieldmap) are smoothed with BSplineApprox. That has also allowed dropping the old niflow and a bunch of "artisanal" settings from the good-old, mitical epiunwarp.fsl script of MGH. Now, fieldmaps are more consistent with TOPUP outputs, as both options have the same relevant outputs:
  • a fieldmap reference (magnitude file / EPI file)
  • a preprocessed fieldmap (in Hz), given in alignment with the fieldmap reference
  • coefficients file(s), in Hz, given in alignment with the preprocessed fieldmap (for the case of BSplineApprox, only) - this is where we need to work on TOPUP's coefficients file, to store the BSpline grid separation and rotation w.r.t. canonical on the NIfTI's header.
• [Workflows] Revised the init_coeff2epi_wf workflow so that coefficients are moved into the target EPI space without resampling them. To do so, the affine of the NIfTI file is rotated and recentered following ANTs registration outputs (see https://github.com/nipreps/sdcflows/pull/119/files#diff-7121558e9d85fc90428b6d6b2f87d9586e7ff8682789782d847a831c12e80b1eR130). IMHO this is one of the coolest tiny details of this new implementation that will make it killer altogether. One test runs an SBRef-to-magnitude registration, see, e.g., this report: https://1436-189236569-gh.circle-artifacts.com/0/tmp/tests/sdcflows/sub-HCP101006/figures/sub-HCP101006_space-sbref_fieldmap.svg
• [Workflows] Implemented a draft of init_unwarp_wf where only the coefficients moved (without resampling, remember) on to the target's space are used to reconstruct the fieldmap. This is a real breakthrough w.r.t. existing implementation of fieldmaps, esp. non-PEPOLAR:
  • The fieldmap resampled using the projected coefficients looks pretty good (see https://1436-189236569-gh.circle-artifacts.com/0/tmp/tests/sdcflows/sub-HCP101006/figures/sub-HCP101006_task-rest_dir-LR_desc-fieldmap_bold.svg). For reference, the original coefficients file was generated by another test - see https://1436-189236569-gh.circle-artifacts.com/0/tmp/tests/sdcflows/sub-HCP101006/figures/sub-HCP101006_desc-phasediff_fieldmap.svg - the apparent differences in intensity are caused by the more liberal mask of the SBRef, compared to the magnitude.
  • The correction does look good: https://1436-189236569-gh.circle-artifacts.com/0/tmp/tests/sdcflows/sub-HCP101006/figures/sub-HCP101006_space-sbref_fieldmap.svg - don't get distracted by those funny artifacts outside the brain, that's because the extrapolated field there is very strong and signal is being picked up from far away. However that should be masked out. The important issue here is having the brain match the anatomical image, and that is clearly achieved with a pretty low-resolution BSpline grid. I would expect better results when running in production where we can apply a two-levels BSplines pyramid.
• Tests - as codecov registers, the diff is 89% covered (which is pretty nice). I'm sure some additional tests could be added, but I think for the time being this is already great.
• Documentation - revised existing documentation and improved some formulas. Attempted to make these formulas more descriptive of our implementation. Added lots of documentation about the B-Splines implementation.

NB. AFAIK, the implementation of B-Spline in numpy/scipy fall short for what we needed here. For the 2D case I could have reused very specific spline evaluation functions but nothing exists for the 3D case that we could have leveraged.

[pep8speaks]

Hello @oesteban, Thank you for updating!

Cheers! There are no style issues detected in this Pull Request. 🍻 To test for issues locally, pip install flake8 and then run flake8 sdcflows.

Comment last updated at 2020-11-25 19:50:24 UTC

[oesteban]

Okay, this is finally working. Please check the reports after this smoothing:

• Phasediff 1: before (master) and after (this PR)
• Phasediff 2: before (master) and after (this PR)
• Phasediff HCP: before (master) and after (this PR)

I think the results are pretty cool. Please note that, in debug mode the fieldmap is not extrapolated outside the brain mask, and that the B-Spline grid is single and with the control points pretty far apart (that means the generated field will filter out high-frequency components of the original data).

There's perhaps one missing check to do - there might be a scaling component because the field is centered about zero, it seems that the dips and peaks could be fitted a bit more closely.

[effigies]

I haven't worked through the workflow, but here are some comments on the interfaces. Please ping me again if I don't finish this tomorrow.

[oesteban]

Rebased on top of #123 - oesteban/sdcflows@enh/new-unwarp-api...oesteban:enh/bspline-regularization

[codecov]

Codecov Report

Merging #119 (daa8930) into dev/1.4.0 (f137b47) will increase coverage by 5.49%.
The diff coverage is 88.78%.

@@              Coverage Diff              @@
##           dev/1.4.0     #119      +/-   ##
=============================================
+ Coverage      71.30%   76.79%   +5.49%     
=============================================
  Files             14       16       +2     
  Lines            784      935     +151     
  Branches          96      110      +14     
=============================================
+ Hits             559      718     +159     
+ Misses           200      187      -13     
- Partials          25       30       +5     

Flag	Coverage Δ
travis	63.85% <79.51%> (+5.55%)	⬆️
unittests	76.71% <88.78%> (+5.52%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
sdcflows/interfaces/fmap.py	67.39% <ø> (ø)
sdcflows/workflows/fit/pepolar.py	63.51% <ø> (ø)
sdcflows/workflows/apply/registration.py	56.75% <6.66%> (-25.86%)	⬇️
sdcflows/interfaces/epi.py	82.53% <62.50%> (-3.67%)	⬇️
sdcflows/interfaces/bspline.py	96.15% <96.15%> (ø)
sdcflows/interfaces/reportlets.py	81.35% <100.00%> (+1.35%)	⬆️
sdcflows/workflows/apply/correction.py	100.00% <100.00%> (ø)
sdcflows/workflows/fit/fieldmap.py	98.07% <100.00%> (+31.41%)	⬆️
sdcflows/workflows/outputs.py	67.24% <0.00%> (-3.45%)	⬇️
... and 2 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update f137b47...daa8930. Read the comment docs.

[oesteban]

@mgxd after the latest force-push to remove a little update of docs/conf.py that I have moved directly into master, I think this is ready for you to have a look. I have edited the PR's message to include some information that I think will help you navigate through this monster.

Comments from @effigies and @mattcieslak (or anyone else) are also welcome of course (although I may defer to the moment we merge dev/1.4.0 into master for the feedback if I feel that merging this will allow me to get further, and I know Chris is not available these days, for a start).

[mgxd]

Some comments / suggestions from my first pass. Overall I think this looks great.