A new workflow for calculating EPI references

[oesteban]

With the integration of new modalities (dMRIPrep) and specimens (babies, rodents) we are realizing of two problems with how fMRIPrep generates "BOLDRefs":

• The workflow does not generalize well to DWIs
• The workflow has a masking step that is atlas-based, which impedes generalization to new specimens (babies, rodents) and may not be reliable with DWIs of human adults.
• The workflow is currently unaware of SDC

Here, I propose to discuss the following workflow to be implemented in niworkflows, to be used by all downstream nipreps:

• The workflow takes away the masking part:
  • as this will depend on the specimen and the contrast, it is better not to address it at this point, and
  • because the mapping to the anatomical reference is estimated, a pretty precise mask is generated by projecting it from anatomical space.
• The workflow generalizes across modalities, and in particular, depends on Make EstimateReference modality-agnostic #600 for the manipulation of potential references.
• The workflow uses the new API of SDCFlows to apply unwarping in one step. This unwarping is by-passed if there's no fieldmap to be applied (e.g., you don't have fieldmaps or you indeed are estimating an EPI reference for some fieldmap-less estimation). NOTE: the first prototype of the workflow should be designed without this SDC connection just not to make it overly complicated for a proof-of-concept.
• We will need to come up with a robust heuristic for the final reference grid. Probably just the highest resolution of each of the dimensions (reusing our EstimateSamplingGrid interface). Although not common, the runs could theoretically have different sampling grids (probably more likely to happen with DWIs).

Assigning Mathias & @eilidhmacnicol.

Requesting feedback from @effigies, @mattcieslak, @josephmje, @dPys, @tsalo.

[oesteban]

Actually, I'm thinking that the last (BBR) / ANTs could be better addressed outside this particular workflow. For instance, if we want to eventually implement an anatomical-free pipeline (for legacy datasets), then the final alignment to anatomical is not necessary.

[effigies]

If I'm understanding correctly, this unifies essentially all steps up to T1w registration, which we currently have as estimate reference, run hmc and apply sdc. A few notes:

1. Reference finalization might want to be more customizable than selecting the highest SNR (I'll call it maxSNR). This would be a good place to implement mean EPI (Use subject specific BOLD template for T1w coregistration fmriprep#450), which I think would basically be this, but instead of using maxSNR in its orig space, we would use mri_robust_template to find a central location and resample all runs to that location.
2. Do we want to resample all runs into a reference grid before estimating EPI-T1w? What's the advantage over just pulling forward the transformation list?
3. What's the advantage of estimating HMC in this workflow? If I want to parallelize fMRIPrep invocations over BOLD runs, each job would need to calculate the EPI reference (which takes all runs) but the motion estimates would be wasted duplicate effort.
4. I agree with not attempting EPI-T1w registration here. EPInorm (EPInorm to replace EPI->T1w->MNI registration fmriprep#620) would be much easier to implement if we keep this just to EPI.
5. I believe the outputs of this workflow should be:
   • EPI reference (1)
   • STC (or reference to original) EPI series (per-run)
   • HMC transforms (per-run) if calculated
   • SDC transforms (per-run)
   • EPI-to-reference transforms (per-run)

[oesteban]

this unifies essentially all steps up to T1w registration, which we currently have as estimate reference, run hmc and apply sdc

Yes, this seems like a fair description.

1. Reference finalization might want to be more customizable than selecting the highest SNR (I'll call it maxSNR). This would be a good place to implement mean EPI (nipreps/fmriprep#450), which I think would basically be this, but instead of using maxSNR in its orig space, we would use mri_robust_template to find a central location and resample all runs to that location.

Seems very reasonable. Let's start with the maxSNR idea but being mindful that it should be easily replaceable by alternatives.

2. Do we want to resample all runs into a reference grid before estimating EPI-T1w? What's the advantage over just pulling forward the transformation list?

Yes, the idea is to run single-channel (T1w template) to single-channel (the overall EPI average) registration - and hence, run antsRegistration/bbregister just once.

With this question, I think you are pointing at that node with the label "Resample displacements field into each reference and volume", right?
If so, the idea is to be able to correct each volume selected across runs with their particular (and unique) realization of the distortion. In other words, I'd like to account for the fact that each volume has suffered its own different distortion. Definitely not crucial for humans (and even less in rats, which are more restrained and move less) but it is theoretically the right thing to do (and can give the necessary accuracy on very limited FoV and populations with substantial motion).

3. What's the advantage of estimating HMC in this workflow? If I want to parallelize fMRIPrep invocations over BOLD runs, each job would need to calculate the EPI reference (which takes all runs) but the motion estimates would be wasted duplicate effort.

General case: You can skip HMC only under the assumption that, because you are selecting contiguous volumes (hence, little chance for large movements) in the beginning of the scan (hence, the subject is not yet so uncomfortable, less motion), the selected volumes will be aligned (or that HMC is going to actually introduce motion rather than correcting). However, that is not the case in DWI (b=0 references are typically scattered along the sequence), and may not be the case for BOLD - if you don't select just a few volumes. I agree this is something we might want to consider more carefully. That said, the first HMC pass (i.e., intra-run) should be optional and therefore we could set a skip for references selected with the nonsteady state detector code (for instance, to keep it closer to current fMRIPrep's behavior).

Directly answering the proposed use-case: I guess I have little to counter-argue. One of the assumptions under the proposed workflow is that you will obtain a session-wise reference of the EPI (one longstanding request from a good bunch of fMRIPrep users), and it seems that use case would advise against parallelizing runs (because you need to ensure all those identical HMC runs will be repeatable to the last bit between processes). Other than that, I guess that running 3dVolReg on the same 5 volumes 5 times is not computationally the end of the world or a big waste.

If still interested in the parallelization and worried about replicability across parallel processes, we could recommend those users to use bids-filters for the parallelization so that the EPI references generated are single-run.

4. I agree with not attempting EPI-T1w registration here. EPInorm (nipreps/fmriprep#620) would be much easier to implement if we keep this just to EPI.

Yup, forgot to mention EPInorm here but I definitely had it in mind when sketching out the workflow.

5. I believe the outputs of this workflow should be:

   • EPI reference (1)
   • STC (or reference to original) EPI series (per-run)
   • HMC transforms (per-run) if calculated
   • SDC transforms (per-run)
   • EPI-to-reference transforms (per-run)

I would add a trivial output that is the same as the list of inputs, but dropping their paths. All the multiple outputs will have the same length so that downstream workflows can select a set of outputs using the KeySelect interface (very much like we are selecting standard spaces generated by the anatomical workflow). I think you were going on a similar direction with the ("reference to original EPI series (per-run)"). It could come handy to also have an output with booleans indicating which one is the reference (or one output with the key = original file name without path of the reference)

Is STC = slice timing correction? If so, this would be the only one I'd take out of the list. AFAIK, DWIs are not corrected for slice timing in general (cc/ @arokem, @garikoitz, @josephmje, @mattcieslak). I might be missing the point of having it here, so please elaborate more. Forwarding a reference to the original EPI series seems like something reasonable that may also help Nipype internal signaling.

[oesteban]

(thanks a lot for the feedback!)

[effigies]

Reordering for clarity.

2. Do we want to resample all runs into a reference grid before estimating EPI-T1w? What's the advantage over just pulling forward the transformation list?

Yes, the idea is to run single-channel (T1w template) to single-channel (the overall EPI average) registration - and hence, run antsRegistration/bbregister just once.

With this question, I think you are pointing at that node with the label "Resample displacements field into each reference and volume", right?

No, I'm talking about "Resample all runs into reference grid". In general, fMRIPrep has resampled the references only until a final single-shot interpolation could be done. If this approach requires moving to multiple interpolation steps, I'm not sure it's worth the generalizability.

But I don't think it's necessary. As long as we prepare the EPI-T1w registration reference, and keep track of the displacements from the per-run references to that one, then that seems like it's enough. We then concatenate (HMC, SDC, EPI-to-reference, reference-to-T1w, T1w-to-template) as usual.

Is STC = slice timing correction? If so, this would be the only one I'd take out of the list. AFAIK, DWIs are not corrected for slice timing in general [...]. I might be missing the point of having it here, so please elaborate more.

The point is that fMRIPrep currently does head motion estimation -> slice timing correction -> head motion correction, which is the best we can do short of 4D realignment to avoiding letting either step interfere with the other. I agree that it doesn't seem applicable to DWI, but if there are any corrections or estimations that need to be handled similarly, we could make this pluggable as well.

3. ...

...

It appears that what you're describing is running HMC on reference scans only. If that is the case, and we will calculate intra-run HMC elsewhere, then I misunderstood.

That will also resolve the HMC/STC issue just discussed, since those will be handled outside this workflow.

To try to restate my position clearly: We should run the minimum required HMC in order to effectively calculate the reference volumes. Full-run HMC should either not take place within this workflow (my preference), or we should at least code in such a way that it would not be difficult to add a feature where we perform full-run HMC on only a subset of runs.

[oesteban]

No, I'm talking about "Resample all runs into reference grid". In general, fMRIPrep has resampled the references only until a final single-shot interpolation could be done. If this approach requires moving to multiple interpolation steps, I'm not sure it's worth the generalizability.

Oh, right. The plan is definitely not using this reference as a source to accumulate transformations, this is fundamentally to have the best reference possible for the alignment to the anatomical image and HM estimation of full-runs (see towards the end of this comment).

The point is that fMRIPrep currently does head motion estimation -> slice timing correction -> head motion correction, which is the best we can do short of 4D realignment to avoiding letting either step interfere with the other. I agree that it doesn't seem applicable to DWI, but if there are any corrections or estimations that need to be handled similarly, we could make this pluggable as well.

Yup, but please note that, as proposed, the workflow does not really resample full timeseries, just the selected volumes to build the reference ("template" if you wish to use a nomenclature consistent to sMRIPrep's). I'm not very sure that missing STC of the selected reference volumes will make a big difference... but it is not undoable, that is true.

Regarding the STC problem itself, I'm pretty convinced that allowing 4D transformations in NiTransforms is the best route - especially because the nipy code is a bit out of shape. But even if you resurrect it, you will want to be able to also apply other transforms at the same time.

It appears that what you're describing is running HMC on reference scans only. If that is the case, and we will calculate intra-run HMC elsewhere, then I misunderstood.

Yes, here HMC is proposed for only the reference volumes that end up building the "session-template". This "session-template" is then the target for HM estimation of each individual run (we could actually generate a reference without SDC to make this HM estimation easier, otherwise they will be potentially mismatched and there's no need for that).

Full-run HMC should either not take place within this workflow (my preference)

Yup, confirming this is the actual proposal :) - it is not done in this workflow.

[eilidhmacnicol]

This has been quite a high-level discussion, so apologies for taking so long to comment; I've been getting my head around the ideas suggested from across the framework.

If I understand correctly, once #602 is merged, then we need to implement it in a workflow that first arranges inputs by phase encoding direction (what is t_RO? read out time?) into a list of lists, which are fed into RobustAverage, then N4/clipping (via init_enhance_and_skullstrip_bold_wf) and then EstimateReference (which still needs updating since it was taken out of RobustAverage, as I understand it?). This leaves us with implementing rigid registration as the last part of this workflow.

Then, in a higher-order workflow, the outputs of the rigid registration are passed forward for SDC and registration (either to an anatomical image, or to an EPItemplate in the future case of EPInorm).

Have I got that right?

The only thing I have to note, is that if we are still using init_enhance_and_skullstrip_bold_wf for the N4 and clipping is that the workflow contains BET, so is not rodent/babies friendly in its current form.

[oesteban]

@mgxd can you port your prelim work over here (https://github.com/nipreps/nibabies/blob/582d07f530e0c9b33151677f98db9c84c9da4da0/nibabies/workflows/bold/reference.py) and start using #609 for integration testing?

I'd do it myself, but I don't want to take credit for your commit - I guess I could override the Authored-by meta in my commit message?

[oesteban]

Okay, done! I'll give it a crack.

[oesteban]

A baseline workflow for this has been merged (#610). Please follow up with its offspring of issues (#611, #612, #613, #614, #615)