ENH: New estimation API (featuring a TOPUP implementation!)

[oesteban]

Summary

Builds a new API that is uniform across implemented estimation methods:

• PEPOLAR: covering the default 3dQwarp, and adding FSL TOPUP
• Direct B₀ maps (with Spiral Echo Imaging - SEI)
• Phase-difference maps (with Gradient Recalled Echo - GRE sequences)
• Fieldmap-less (SyN-based)

CHANGES

• Start a new module sdcflows.interfaces.epi - There was a conceptual misunderstanding about the echo spacing / readout time. These inter-related timing parameters are relevant for all EPI datasets that are targetted for correction, as well as in PEPOLAR inputs because it scales the amount of distortion at each voxel. We had these interfaces within the fieldmap module (sdcflows.interfaces.fmap), which is not correct. You need the parameter fundamentally for unwarping, and it so happens that it is necessary when fieldmapping with PEPOLAR. A new module seems cleaner methodologically.
• Improved the FieldmapReportlet interface to allow displacement fields (i.e., 4D NIfTIs) to be plotted, just for the convenience of not having to extract the nonzero dimension.
• A new sdcflows.interfaces.utils.Flatten interface that unravels lists of 3D and 4D EPIs into a list of 3D EPIs with corresponding metadata expanded too. This is useful for PEPOLAR corrections and replaces hacky parts of the old preparation workflow. Potentially upstreamable to niworkflows
• The core of the PR: a new sdcflows.models submodule. This new module contains the estimation methods. The vision is to have a scikit-learn-like approach with some sort of fit() and predict() organization. sdcflows.models implements the fit() part of it. Inside we find:
  • sdcflows.models.fieldmap: contains both direct B₀ maps AND phase-difference maps. This was almost the case already, but some minimal changes have contributed to a more streamlined workflow. Phase-difference maps are massaged (phase-unwrapped, denoised, masked, converted to rad/s) and then passed in into the general direct B₀ workflow. If the input is a direct B₀, the phase-difference preprocessing is just skipped over. Both strategies share the exact same magnitude components processing.
  • sdcflows.models.pepolar: simplifies the existing workflow for 3dQwarp assuming the inputs are more massaged, in order to allow the FSL TOPUP implementation to have both the same interface.
  • sdcflows.models.syn: just transferred from the original sdcflows.workflows.syn
• sdcflows.workflows.base: updated to work with the new models submodule, and eliminated the many heuristics and validation checks that will now be done with the new features in ENH: New objects for better representation of fieldmap estimation #114.
• Updated tests in an attempt to maintain code coverage.

TODO (on a separate PR)

• A B-Spline fitting module:
  • Finishing up Improve/finish Phasediff SDC regularization with B-Splines #14
  • Regularizing the estimation of 3dQwarp (similarly to Improve/finish Phasediff SDC regularization with B-Splines #14).
  • Giving all estimation models a consistent output with two files (a la TOPUP): field coefficients and an actual field, scaled by the readout time.

[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-01 15:01:18 UTC

[oesteban]

Brace yourselves, this is going in - I want to test #119 ASAP.

[mattcieslak]

Question - were you able to figure out how to get the topup warp as an itk tranfsform out of the fsl outputs?

[oesteban]

Regarding the warp, @smoia had it almost finished in nipy/nitransforms#51 - I will try to push this to the end line shortly. But in principle, warpings are not of much concern.

The field coefficients are a different story. In order to make the output of all estimations consistent, all of them should have the same outputs of TOPUP, which means a fieldmap (I'm going to avoid calling it a warp, because that means you know the readout time of the target image) and the coefficients of a tensor B-Spline field supporting such field (#119).

I know this second part is going to be tricky because TOPUP dumps out a NIfTI where header information has not been adapted to the target image. So after finishing #119, I'll also have to write something that adapts FSL coefficients to a NIfTI with headers consistent with the original (interpolated) field. This shouldn't be super hard though.

With those two pieces, I think the new API will be pretty solid:

1. Estimate the field with as many fieldmap opportunities you encounter (and the BIDS PR that is in the works)
2. Dump the reference/magnitude, the field, and the coefficients to the fmap/ folder of the BIDS outputs (this is actually done by the downstream library, e.g., fMRIPrep, but it is good to mention here, and an "output" workflow that will be instantiated from the downstream library will also be written in this proper)
3. Align the reference/magnitude image to the target image to be corrected. Resample the BSpline coefficients on the target space (but keeping, of course, the control point separation).
4. Interpolate the field on the grid of the new image.
5. Scale the interpolated field with the actual readout-time/effective-echo-spacing of the target image
6. Convert to displacements field (amenable to nitransforms, and hence, to ITK)