Nipype should be a wrapper, not a scaffold

[oesteban]

A quick and lazy glimpse into the code gave me the impression that the meat of the different algorithms is stuffed into nipype's _run_interface method implementations. Then, the different components are glued together with nipype workflows.

I think that will become problematic in two ways:

• Non-nipype programmers will find a tall entry barrier on this.
• How do you plan to migrate to Pydra when the time comes?

I think that developing a clear API that can easily wrapped inside Nipype interfaces, or integrated with any other code, would be beneficial. For instance, if you want to have a command line interface it seems to me that Nipype is too heavy of a dependency that helps nothing to that.

WDYT?

[oesteban]

If you agree, I think we should open a Wiki page on this repo and sketch out this API. This document must define:

• Inputs: a list of bare inputs that are necessary to run this. By bare I mean that no massaging should happen. For example:

  • Gradient table in RASB format, normalized.
  • Flag whether the acquisition is shelled, non-shelled or DSI.
  • Reference b=0 - e.g., like the one dMRIPrep generates that is pretty sophisticated (head motion, drift-regression, etc).
  • DWIs, e.g., without b=0?,
  • Model specification

• Outputs:

  • List of composite (HMC+EMC) affines, one per DWI
  • List of decomposed HMC rigid-body affines, one per DWI.
  • List of residual affines (HMC+EMC - HMC), one per DWI (somewhat encoding only-eddy distortions).
  • List of jacobian maps?
  • Any other metrics we can think could be useful as an output

Something like this?

[oesteban]

Then, for the pipelining I also think we should give up on nipype. It feels to me that a smart utilization of python's parallelization (any appetite to test asyncio?) would be better.

We could directly do this with pydra as well (which will use asyncio underneath), I don't know how @dPys would feel about this. Pretty sure @satra will encourage us to do so :D.

[dPys]

Thanks for fleshing out a deeper discussion on this @oesteban .

My immediate first thoughts are that this probably does still need to be housed in a workflow engine of some kind. It is a workflow fundamentally, so If not nipype, then probably Pydra if it's ready for it. The idea for asyncio is interesting because there might be ways still to send some of the heavy memory labor to disk. As for how to go about those kinds of sweeping refactorings more practically-speaking, it's a totally separate question since as you know my bandwidth is still pretty limited at the moment

The only required inputs for now are a dwi series with any number of shells, and at least one B0, bval, and bvec (see the jupyter demo). A RASB input as an alternative option to bval/bvec is low-hanging fruit and so should probably be added asap. The pipeline detects shelled-ness automatically so probably don't need to specify it manually. A reference B0 is actually currently generated using dmriprep's init_b0_reference_wf (see https://github.com/dPys/EddyMotionCorrection/blob/main/emc/workflows/correct.py#L703), but another optional cli input could be a custom B0 specified by the user or pulled from elsewhere? Although the eddy correction does require a B0 reference of some kind, theoretically the LOO HMC should still succeed without it. That said, any particular reason to support the no-B0 case? A flag for model specification is definitely needed, but should be easy to include. And on that note, an easy-to-use cli should be a priority. Another important input for us to include will be any SDC transforms created beforehand, since these could be applied to improve the quality of the registrations if they are available. Also, support for multiple dwi runs...

I agree with all of the proposed outputs listed.

Let's create a wiki page!

[oesteban]

A RASB input as an alternative option to bval/bvec is low-hanging fruit and so should probably be added asap. The pipeline detects shelled-ness automatically so probably don't need to specify it manually. A reference B0 is actually currently generated using dmriprep's init_b0_reference_wf

Yes, the idea is to remove as many responsibilities as possible away from this package. The package should be good at just one thing.

• The RASB format is probably the best way to avoid internal manipulations. Also it is easy to implement by others.
• If we are using the B0 workflow from dMRIPrep, let's not copy the code here. Assume you get a B0 with same resolution and matrix as the inputs (or whatever other requirements you may come up with if you think it through).
• Shelled-ness - again, it's better not bothering here. All these "addons" really add brittleness to the code and distract away from the real value.
• Workflow: it seems to me that we are not really running any compute graph, but rather iterating over just a few sequential steps. Not sure a workflow engine is not an overkill here. Disgregating the method also will pose problems if you want to try to use third-party optimizers (say we want to try with some Bayesian optimization).

I can have a crack at refactoring all of this, if you like.

[dPys]

A RASB input as an alternative option to bval/bvec is low-hanging fruit and so should probably be added asap. The pipeline detects shelled-ness automatically so probably don't need to specify it manually. A reference B0 is actually currently generated using dmriprep's init_b0_reference_wf. And in general, there shouldn't be any redundant code with what is already in dmriprep or niworkflows, but I could be wrong.

Yes, the idea is to remove as many responsibilities as possible away from this package. The package should be good at just one thing. Yep, agreed.

• The RASB format is probably the best way to avoid internal manipulations. Also it is easy to implement by others. Agreed, though wouldn't we still want to be sure that users who haven't yet conformed to the still new RASB standard can also use the tool as well?
• If we are using the B0 workflow from dMRIPrep, let's not copy the code here. Assume you get a B0 with same resolution and matrix as the inputs (or whatever other requirements you may come up with if you think it through).

So the B0 wf is actually not copied here, it's imported directly from dmriprep. But I 100% agree that eventually that portion should be removed and left up to the user! I had just included it initially to ease development.

On that note as well, I should mention that my vision for this package is to free itself entirely from all non-pypi 3rd party dependencies (e.g. FSL, ANTs, freesurfer, AFNI). Without the dmriprep B0 wf imported, we are nearly there. Currently, all registrations are done in-house: https://github.com/dPys/EddyMotionCorrection/blob/main/emc/utils/register.py and Dipy is used whenever possible.

• Shelled-ness - again, it's better not bothering here. All these "addons" really add brittleness to the code and distract away from the real value. So I might have miscommunicated this-- I just meant counting the number of unique bvals. It's checked only in one place from the gtab to determine which diffusion model is appropriate for the HMC. Of course, we could also leave this up to the user and assume the user knows in advance which models are compatible with their data or not, and then raise an error if incompatibility is detected? :-)

• Workflow: it seems to me that we are not really running any compute graph, but rather iterating over just a few sequential steps. Not sure a workflow engine is not an overkill here. Disgregating the method also will pose problems if you want to try to use third-party optimizers (say we want to try with some Bayesian optimization).

I can have a crack at refactoring all of this, if you like.

There are actually at least 10 nested workflows involved at the moment, which make use of 5 MapNodes that seem like they would be difficult to emulate with everyday python, but it may be possible. Might be worth taking a deep-dive into https://github.com/dPys/EddyMotionCorrection/blob/main/emc/workflows/correct.py. That said, if you can think of a good way to condense this without nipype, I'm totally game since as you mentioned before it would reduce barriers to contributing for other future developers. Feel free to PR a WIP?

My initial goal was just to get a prototype EMC tool other than FSL's Eddy that works. After that, sky is the limit. We do still need to run some benchmarking to see how well the corrections are actually doing, however. From pure visual inspection, it appears to be doing a very good job, even when using just the tensor model, but exactly how many iterations of the LOO are most optimal and using which specific models to do the HMC prediction is still unclear. If the package performs as well or better than Eddy, for example, then it would almost definitely be worth the effort to do formal code cleanup and get this out of alpha status. If not, then maybe not. So a benchmarking study is needed, badly.

p.s. I love the idea of incorporating Bayesian optimization. What exactly did you have in mind?

[oesteban]

Cool, seems we are heading in the same direction :).

Okay, before I proceed to suggest changes (a PR), let me double-check I'm getting this correctly. Say we have a DWI with 8 b0 + 64 high-b directions, 2-shells for the following questions:

1. The SignalPrediction MapNode: https://github.com/dPys/EddyMotionCorrection/blob/ea730f4de9e0ceeeec315e9cb9d169e6963dce72/emc/workflows/correct.py#L268-L275 This MapNode generates the corresponding 64 directions simulated in a LOO loop, where the remaining 63 are used in training a particular model. Is this correct?

   Nitpick: Why you don't want to train with q-samples that are close to the predicted? Is over-fitting a problem here (I'd say no, right)? https://github.com/dPys/EddyMotionCorrection/blob/ea730f4de9e0ceeeec315e9cb9d169e6963dce72/emc/interfaces/model.py#L76-L80

2. Register the 64 original DWIs to the predicted ones: https://github.com/dPys/EddyMotionCorrection/blob/ea730f4de9e0ceeeec315e9cb9d169e6963dce72/emc/workflows/correct.py#L284-L291 I'm interpreting that synchronize = True ensures we are doing the dot-product and not the tensor-product between the two lists. I.e., DWI with index 3 is registered with the simulated index 3 exclusively. Is that correct?
3. The rest of the code: 3.1. Initializes data and settings (e.g., that Patch2Self in disguise - is denoising a must-have for this idea to work?) 3.2. Composes iterations over the workflow 3.3. Massages outputs (e.g., that CombineMotions node?)

Did I get it at least close?

[dPys]

Cool, seems we are heading in the same direction :).

Okay, before I proceed to suggest changes (a PR), let me double-check I'm getting this correctly. Say we have a DWI with 8 b0 + 64 high-b directions, 2-shells for the following questions:

1. The SignalPrediction MapNode: https://github.com/dPys/EddyMotionCorrection/blob/ea730f4de9e0ceeeec315e9cb9d169e6963dce72/emc/workflows/correct.py#L268-L275

   This MapNode generates the corresponding 64 directions simulated in a LOO loop, where the remaining 63 are used in training a particular model. Is this correct?

Correct. And I'm currently initializing this procedure using a tensor model on the first LOO pass, followed by subsequent iterations with the model of the user's choice. It is precisely this stage where more optimization can conceivably be done (e.g. stacking multiple diffusion models, bootstrapped aggregation, etc.).

Nitpick: Why you don't want to train with q-samples that are close to the predicted? Is over-fitting a problem here (I'd say no, right)? https://github.com/dPys/EddyMotionCorrection/blob/ea730f4de9e0ceeeec315e9cb9d169e6963dce72/emc/interfaces/model.py#L76-L80

Overfitting is a problem here. If the training samples are too close along the sphere to the sample to predict, it will bias the prediction for particular directions, particularly if the vectors are not equally distributed along the sphere.

On a slightly related note one thing that could be interesting w.r.t. a Bayesian approach might be to try drawing additional, simulated vectors from the posterior distribution of vectors in order to get more ("pseudo"-) training/samples for prediction. Perhaps priors could be set based on the outputs of the first LOO iteration with the tensor model...

2. Register the 64 original DWIs to the predicted ones: https://github.com/dPys/EddyMotionCorrection/blob/ea730f4de9e0ceeeec315e9cb9d169e6963dce72/emc/workflows/correct.py#L284-L291

   I'm interpreting that synchronize = True ensures we are doing the dot-product and not the tensor-product between the two lists. I.e., DWI with index 3 is registered with the simulated index 3 exclusively. Is that correct?

Correct, it ensures that we are synchronizing the indices across the raw DWI with those in the predicted DWI, because we want the dot-product.

3. The rest of the code: 3.1. Initializes data and settings (e.g., that Patch2Self in disguise - is denoising a must-have for this idea to work?)

Denoising is not a must-have, but as Matt once noted (and I found also while developing this), it does seem to greatly improve both the quality of the registrations and the quality of the predictions. The nice thing about this approach is that the final motion-corrected outputs are actually not directly denoised, which would risk changing various assumptions about the data that some users might not like. Nevertheless, they still reap all the benefits of the learning gains from the denoising, because that learning is transferred indirectly -- it improves the quality of the transforms that are ultimately applied to the raw dwi series at the end of the workflow.

3.2. Composes iterations over the workflow 3.3. Massages outputs (e.g., that CombineMotions node?)

Yep, and yep. And generate reports on the outputs, which should be expanded upon because they are still very bare-bones.

Did I get it at least close? You're right on the money. The only thing I might add is that the eddy correction is performed before the HMC to ensure that, like with the denoising, the LOO prediction is drawing upon training examples that have already scrubbed of eddy currents that would otherwise greatly detract from the quality of the diffusion model fit/predict methods in HMC. Optional SDC and Gibbs de-ringing beforehand should conceivably yield similar gains...

[oesteban]

Denoising is not a must-have, but as Matt once noted (and I found also while developing this), it does seem to greatly improve both the quality of the registrations and the quality of the predictions. The nice thing about this approach is that the final motion-corrected outputs are actually not directly denoised, which would risk changing various assumptions about the data that some users might not like. Nevertheless, they still reap all the benefits of the learning gains from the denoising, because that learning is transferred indirectly -- it improves the quality of the transforms that are ultimately applied to the raw dwi series at the end of the workflow.

Consistently with the above, I'm going to recommend taking it out. In general, to follow the principle of simplicity - this library should do just one thing. What if a much better algorithm comes out next month? Moving away from Patch2Self will then take some maintenance burden that is not justified. But, in particular, I'm concerned about a potential double-dipping with Patch2Self - it shares exactly the same LOO scheme to apply regression.

In other words, I would not be surprised that head motion and eddy currents are partially corrected by Patch2Self. So, your simulation is indeed much better, but your estimation of the artifacts is not. I still need to learn more about to Patch2Self, but following today's DIPY tutorial, it seems to me that data MUST be eddy-corrected and realigned for it to meet the noise independence assumption. WDYT @arokem ?

[oesteban]

Overfitting is a problem here. If the training samples are too close along the sphere to the sample to predict, it will bias the prediction for particular directions, particularly if the vectors are not equally distributed along the sphere.

Wait, this is a different problem. One problem is that samples do not uniformly cover the sphere (or half of it) and a different one is discarding nearby samples.

In general, this method will not work well with heavily nonuniformly sampled data. But same works for diffusion models down the line, so there's not much to say here. Users should not expect magic so that garbage-in, garbage-out does not apply to your method. I don't think this is surprising.

However, if the dataset is sufficiently uniformly sampled, then having more density around the left-out sample cannot hurt. Overfitting is good here, because it should give you a better reference for registration (especially if the neighbors are not neighboring in acquisition time).

[arokem]

On Mon, Mar 15, 2021 at 9:37 AM Oscar Esteban @.***> wrote:

Denoising is not a must-have, but as Matt once noted (and I found also while developing this), it does seem to greatly improve both the quality of the registrations and the quality of the predictions. The nice thing about this approach is that the final motion-corrected outputs are actually not directly denoised, which would risk changing various assumptions about the data that some users might not like. Nevertheless, they still reap all the benefits of the learning gains from the denoising, because that learning is transferred indirectly -- it improves the quality of the transforms that are ultimately applied to the raw dwi series at the end of the workflow.

Consistently with the above, I'm going to recommend taking it out. In general, to follow the principle of simplicity - this library should do just one thing. What if a much better algorithm comes out next month? Moving away from Patch2Self will then take some maintenance burden that is not justified. But, in particular, I'm concerned about a potential double-dipping with Patch2Self - it shares exactly the same LOO scheme to apply regression.

In other words, I would not be surprised that head motion and eddy currents are partially corrected by Patch2Self. So, your simulation is indeed much better, but your estimation of the artifacts is not. I still need to learn more about to Patch2Self, but following today's DIPY tutorial, it seems to me that data MUST be eddy-corrected and realigned for it to meet the noise independence assumption. WDYT @arokem https://github.com/arokem ?

Eleftherios and Shreyas actually just answered that at the Q&A for the workshop. You might expect better algorithm performance if emc is done first, but it's not necessary. That said, I agree that maybe we decouple these things: denoising and emc.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/dPys/EddyMotionCorrection/issues/12#issuecomment-799564813, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAA46NUSVFMZPOP4YMHOK2LTDYZUPANCNFSM4Y5N2ODQ .

[dPys]

Denoising is not a must-have, but as Matt once noted (and I found also while developing this), it does seem to greatly improve both the quality of the registrations and the quality of the predictions. The nice thing about this approach is that the final motion-corrected outputs are actually not directly denoised, which would risk changing various assumptions about the data that some users might not like. Nevertheless, they still reap all the benefits of the learning gains from the denoising, because that learning is transferred indirectly -- it improves the quality of the transforms that are ultimately applied to the raw dwi series at the end of the workflow.

Consistently with the above, I'm going to recommend taking it out. In general, to follow the principle of simplicity - this library should do just one thing. What if a much better algorithm comes out next month? Moving away from Patch2Self will then take some maintenance burden that is not justified. But, in particular, I'm concerned about a potential double-dipping with Patch2Self - it shares exactly the same LOO scheme to apply regression.

In other words, I would not be surprised that head motion and eddy currents are partially corrected by Patch2Self. So, your simulation is indeed much better, but your estimation of the artifacts is not. I still need to learn more about to Patch2Self, but following today's DIPY tutorial, it seems to me that data MUST be eddy-corrected and realigned for it to meet the noise independence assumption. WDYT @arokem ?

The primary motivation for including preliminary denoising here was to improve the quality of the registrations, which it does indeed do. I do see, however, a good argument that it could be like double-dipping for the LOO prediction, particularly if the same regularization is used. Therefore, I'm happy adding it to the chopping-block, at least for now, along with the B0 template + mask generation (since those were intended as temporary anyway). And come to think of it, requiring that the user supply a premade B0 template and brain mask seems pretty reasonable given that's what FSL's EDDY stipulates. It could also afford some nice flexibility in the sense that users could then experiment with group or even population B0-templates for the eddy correction piece...

It's worth reiterating that the quality of prediction is sensitive to the quality of the input data. To me, this means any of three things w.r.t. deciding where, in the broader scheme of things, EMC should commence: 1) the data will already have been maximally preprocessed (i.e. SDC-corrected, denoised, etc.) before the correction and thus EMC happens near the end; 2) the input data is assumed to be basically raw, which (most likely) will result in high prediction / inter-volume registration error; or 3) EMC itself is allowed to be good at more than one thing -- e.g. it preprocesses a training dataset enough to get the most bang-for-our-buck out of prediction, but once the LOO predicted transforms are obtained, they are simply applied to the raw, unpreprocessed data to then resume with other preprocessing steps of the user's choice. I've personally been leaning towards (3), but am very open to discussing more. I'd say there is good reason why FSL packaged eddy correction, SDC, and head motion correction into a single function. The quality of the predictive model will only go so far. Just how far? We don't actually know yet...

Relatedly, it does not surprise me at all that Patch2Self's performance would be more optimal following EMC-- if the input features will be better, so will be the predictions. That said, we can't have things both ways, or at least not efficiently. What we have to figure out w.r.t. order-of-operations is which preprocessing stage consistently has the largest global influence on final data quality (e.g. Eddy correction, HMC, denoising, or something else), and then ensure (i.e. in or out of the context of this package) that the steps leading up to that stage have done everything possible to scrub the input data of all other types of artifact that could detract from the performance of that stage. Luckily, I don't think we humans need to make these decisions at all. I fundamentally agree with @jelleveraart 's original point that our best bet is to let the data speak for itself to any extent we can.

At this stage, I say we focus our efforts on benchmarking first and foremost. From there, we can proceed to conduct a sensitivity analysis to directly quantify how each benchmark varies as a consequence of preliminary preprocessing choices and their order-of-operation. With this approach, we can also better suss out any potential interaction indices with respect to how the data was acquired (i.e. # shells, # directions, q-space sampling scheme). Once those issues are settled statistically, I think that the optimal scope of EMC as a package will crystallize on its own. What do you think?

[oesteban]

@dPys I fully agree with the essence of what you said above. But we also need to be practical here, and I think we diverge in our interpretation or our plan of converting theories and principles into lines of code. Let's try to find agreement there:

At this stage, I say we focus our efforts on benchmarking first and foremost.

I think we first need a basic validation of the idea, before delving into sophisticated analyses. Basically, we need to be able to answer yes to these two questions:

• Can dMRIPrep (or DIPY, or QSIPREP) readily use this EMC implementation (this also includes a clean and easy API that is minimalistic enough to cater to any tool), and it runs in a reasonable amount of time?
• Results look sound on a few test datasets to the eyes of our experts (Ariel, Jelle, yourself)?

After those two "yes", we need to also say yes to the question that originated this effort:

• Is it worth replacing FSL with this package in general?

Anything below the bar of "both perform similarly, but our EMC is Apache 2.0" should make us think.

It's worth reiterating that the quality of prediction is sensitive to the quality of the input data.

Yup, but the interest here is not prediction, is HM and EC estimation. The particular risk of Patch2Self here is that you can get predictions particularly gorgeous but masking head motion. Your prediction will have been applied an hm/ec affine inherently, ruling it out as target for registration in the attempt to estimate the same effect.

Also, you want your method to work even with not-great data. If your method only takes VIP data then it's not probably worth the researcher's time:

• HCP and folks with access to extra-cool Connectome scanners will keep using their current procedures because their building is already pretty tall at this point. They can't afford replacing one foundation.
• The rest of researchers with not-so-great data will not be able to get anything good from your method and not use it.

This is all to say that we don't need to solve everything in one shot. Don't let the perfect be enemies of the good and let's get a valuable scientific MVP quick. We have all the ingredients:

• A well defined, important problem.
• A good idea for improving the current state-of-art, thanks to Matt and their work.
• A good infrastructure to make it work nicely, thanks to DIPY
• A good infrastructure to test integration, thanks to dMRIPrep
• Two brilliant minds at the implementation, thanks to you and Ariel.

So, I see all that benchmarking and integration in dMRIPrep's proper. Let's work on this tool as if it were to make a PR into DIPY's codebase.

WDYT?

[dPys]

@dPys I fully agree with the essence of what you said above. But we also need to be practical here, and I think we diverge in our interpretation or our plan of converting theories and principles into lines of code. Let's try to find agreement there:

At this stage, I say we focus our efforts on benchmarking first and foremost.

I think we first need a basic validation of the idea, before delving into sophisticated analyses. Basically, we need to be able to answer yes to these two questions:

• Can dMRIPrep (or DIPY, or QSIPREP) readily use this EMC implementation (this also includes a clean and easy API that is minimalistic enough to cater to any tool), and it runs in a reasonable amount of time?
• Results look sound on a few test datasets to the eyes of our experts (Ariel, Jelle, yourself)?

I do agree that the pragmatics of the implementation are important, even vital. So I think I can swing meeting halfway on question 1 :-) And as for the second question, I'd say that just based on pure visual inspection, we may even already be there. The corrections are looking pretty darn good as a starting point, even just with the tensor model.

After those two "yes", we need to also say yes to the question that originated this effort:

• Is it worth replacing FSL with this package in general?

Anything below the bar of "both perform similarly, but our EMC is Apache 2.0" should make us think.

Haha, agreed. I think that regardless of how well EMC performs, competing with EDDY will be a steep climb-- they have had nearly a decade to optimize what it's doing, extensive documentation and testing, plus they have the speed of C++ under the hood. On that note, we should seriously consider cythonizing here...

It's worth reiterating that the quality of prediction is sensitive to the quality of the input data.

Yup, but the interest here is not prediction, is HM and EC estimation. The particular risk of Patch2Self here is that you can get predictions particularly gorgeous but masking head motion. Your prediction will have been applied an hm/ec affine inherently, ruling it out as target for registration in the attempt to estimate the same effect.

Let's return to this point later after we've dug into things a bit little more?

Also, you want your method to work even with not-great data. If your method only takes VIP data then it's not probably worth the researcher's time:

• HCP and folks with access to extra-cool Connectome scanners will keep using their current procedures because their building is already pretty tall at this point. They can't afford replacing one foundation.
• The rest of researchers with not-so-great data will not be able to get anything good from your method and not use it.

This is all to say that we don't need to solve everything in one shot. Don't let the perfect be enemies of the good and let's get a valuable scientific MVP quick. We have all the ingredients:

• A well defined, important problem.
• A good idea for improving the current state-of-art, thanks to Matt and their work.
• A good infrastructure to make it work nicely, thanks to DIPY
• A good infrastructure to test integration, thanks to dMRIPrep
• Two brilliant minds at the implementation, thanks to you and Ariel.

So, I see all that benchmarking and integration in dMRIPrep's proper. Let's work on this tool as if it were to make a PR into DIPY's codebase.

WDYT?

I think sounds like a plan. Now in terms of where we go from here, I'm leaving the current working implementation as the main branch until we're ready to replace it. In the meantime, I've hacked together a new WIP in a feature branch called "bare_bones", which takes out denoising and B0 template + mask generation. The latter two images have been relegated to the list of inputs required by the user. Feel free to get started making changes!

[oesteban]

wow 4 additions / 488 deletions (https://github.com/dPys/EddyMotionCorrection/compare/bare_bones)

That's a cleanup!

[dPys]

Quick and easy!

[oesteban]

@arokem @dPys, I'm realizing of a very low-hanging fruit optimization (if not actually an algorithmic improvement):

is there a reason for generating the predictions in bulk? After all, we need to fit the model and generate a prediction for every gradient orientation right? Why not leveraging those orientations that have been already aligned ? This, in addition to randomly choose the next direction to be registered could really improve registration IMHO.

WDYT?

[oesteban]

Okay, this is a sketch of the interface I'm considering:

"""A SHORELine-like algorithm for the realignment of dMRI data."""
import attr
import numpy as np
from pathlib import Path

@attr.s(slots=True)
class EddyMotionEstimation:
    """Estimate rigid-body head-motion and distortions derived from eddy-currents."""

    datapath = attr.ib(default=None)
    """Path to the input DWI dataset."""
    gradients = attr.ib(default=None)
    """A 2D numpy array of the gradient table in RAS+B format."""
    bzero = attr.ib(default=None)
    """A *b=0* reference map, preferably obtained by some smart averaging."""
    model = attr.ib(default="SFM", type="str")
    """The prediction model - options: SFM, SHORE, Tensor, DKI."""
    transforms = attr.ib(default=None)
    """Current transform parameters (list of affine matrices)."""

    def fit(
        *,
        X=None,
        init=None,
        **kwargs,
    ):
        """
        Run the algorithm.

        Parameters
        ----------
        X : :obj:`~nibabel.SpatialImage`
          A DWI image, as a nibabel image object.
        init : :obj:`~numpy.ndarray`
          An Nx4x4 array of initialization transforms.

        """
        loo_index = np.random.shuffle(range(len(X)))

        self.transforms = [
            _emc(X, index=i)
            for i in loo_index
        ]
        raise NotImplementedError

    def predict():
        """Generate a corrected NiBabel SpatialImage object."""

def _emc(x, index, model="SFM"):
    return NotImplementedError

The _emc() function predicts one volume (that pointed by X) and runs registration. Returns the transform parameters. Should accept kwargs for both the model and the registration.

Finally, this example shows the sequential option (therefore the shuffling is not really necessary). If we are going to apply the transforms as they become available, then the shuffling will be important.

comments?

[dPys]

I really like the idea of incorporating the fit-predict methods for this, if anything just for its simplicity and user-friendliness. However, we should be careful not to oversimplify this even in its prototypical form.

A few comments if we are indeed moving in this direction: 1) With each transformation generated in the HMC, don't forget that we'll also need to reorient the vectors accordingly, and that the LOO gets re-implemented across multiple iterations. That is, just one iteration is generally not enough, so we'll want to at least include a num_iters parameter.

2) I can't believe I'm saying this, but now I'm thinking it best that we break up the eddy correction and the HMC a bit more. Currently they are contained in separate workflows. Although these functionalities definitely belong in the same package, they are probably too much for the same class, though maybe not a superclass. One idea would be to do: EddyCurrentCorrect as one class with fit and transform methods. HeadMotionCorrect as one class with fit, predict, and transform methods. This class would also require the transforms generated by EddyCurrentCorrect as init. EddyMotionCorrect, a superclass version of what you posted above that runs the recommended fit, transform, fit, predict, transform methods of the former classes in sequence.

Thoughts on this?

3) Following the above logic, we'll want to be extra sensitive to resource-handling. If we abandon MapNodes altogether, we'll want to replace it with something (e.g. joblib), even from the get-go, to support some type of parallelism. Non-parallelized SFM prediction would take days to run even a single LOO iteration on a 72-direction dwi series. Similarly w.r.t memory, it simply isn't possible to do all of this efficiently in python will all objects loaded into memory. As the package currently does, we'll want to write out at least some intermediate files to disk. Although it might be overkill, right now I'm writing out each and every dwi volume -- predictor and predicted, moving and transformed -- to disk via runtime.cwd, which only just barely helps to keep consumption low enough to run this in LOO-style.

[oesteban]

2. I can't believe I'm saying this, but now I'm thinking it best that we break up the eddy correction and the HMC a bit more.

Open-source/open-science rocks!

support some type of parallelism.

Certainly. One question to @arokem is what parallelism is supported from DIPY itself (both DWI modeling and registration). Maybe we just want to enable that.

Another source of optimization would be to only predict those samples that will really be used in registration (typically, registration will randomly sample the fixed image, so we could restrict prediction to that random sample).

Similarly w.r.t memory, it simply isn't possible to do all of this efficiently in python will all objects loaded into memory.

This is right. I'm thinking of actually having a data-representation object that encapsulates:

• The data array
• NIfTI headers information: orientation (affine, q/s-forms), something else needed?
• B0 reference data array, with voxel coordinates matched with the data array
• A slot for the fieldmap (a displacements field), with voxel coordinates matched with the data array
• Gradient information
• EM transforms

This object will have a to_filename convenience function that writes an HDF5 file - which will allow us to very easily access data without using too much memory.

WDYT?

[oesteban]

You might expect better algorithm performance if emc is done first, but it's not necessary.

This is exactly what I would have expected, and the reason why I think Patch2Self would indeed correct for some head motion and eddy currents (and that's why you "can" run it without emc first).

The problem (see Lindquist et al. 2019 for the case of fMRI), as I mentioned above, is that running EMC after Patch2Self would reintroduce motion and eddy!

[oesteban]

About the data object - what do you think about this - https://github.com/nipreps/nipreps-book/blame/enh/hmc-unit/docs/dmriprep/hmc.md#L97-L131 ?

[dPys]

Love it. And I agree that HDF5 is a great place to start more cleanly grappling with the memory requirements. Have you tried the DWI class yet? I'm tinkering with it now, but it might be helpful to see an example of how you envision we'd implement it here.