FSL's topup implementation (fmriprep >=21) significantly worse than AFNI's Qwarp (fmriprep <20)

[Gilles86]

What would you like to see added in fMRIPrep?

I have noticed that since fmriprep version 21, the PEPOLAR SDC workflow using AFNI's qwarp was replaced by the topup-program of FSL.

During my time working on ultra-high (<1mm) resolution 7T data, I've used both FSL and AFNI for PEPOLAR SDC correction, and I remember getting much better results using AFNI's implementation. I also remember Chris Gorgolewski pointing this out in an earlier version of fmriprep, when he implemented versions with both.

It now has become a concrete problem for me. I just collected a 3T dataset at a modest 2.5mm resolution with LR-encoding. The distortions are pretty bad and fmriprep 21, using fsl's topup, underestimates them severly. The result of this is that the distortions of my LR runs are very different from my RL runs. With fmriprep 20 this problem was way smaller.

Note how in fmriprep 21, the gray matter in the EPI image really extends outside of the GM/CSF border in right frontal cortex. The issue is much less-pronounced in fmriprep 20. (note that fmriprep 21 does report having applied the pepolar SDC; also see the github repo below)

fmriprep 20

fmriprep 21

So:

1. I am just curious: what was the decision to switch to fsl's topup based on?
2. More importantly, I think it would be great if there would be a command line switch to apply AFNI's PEPOLAR SDC workflow rather than FSL's. I noticed that in sdcflows, both are currently implemented. So this shouldn't be too complex:
   https://github.com/nipreps/sdcflows/blob/master/sdcflows/workflows/fit/pepolar.py
3. I personally would even argue that AFNI's implementation of pepolar SDC should be the default, but I suspect there were good reasons to choose FSL's implementation as the new default.
4. As an aside, I also noticed that the visual representation of the SDC in the fmriprep html-report now only has a contour of the outline of the brain mask. I personally think this is much less informative from the GM/WM and GM/CSF boundaries, which were used fmriprep 20 and below.

Do you have any interest in helping implement the feature?

Yes, but I would need guidance

Additional information / screenshots

You can find a bunch of examples of the problem in my github repo here:
https://github.com/Gilles86/topupvsafni

[effigies]

Hi @Gilles86. Thanks for this report.

1. I am just curious: what was the decision to switch to fsl's topup based on?

We found a number of datasets where 3dqwarp performed significantly worse than TOPUP. @oesteban may remember more clearly, but I believe the issue was that TOPUP performed better on high spatial-frequency distortion.

More importantly, I think it would be great if there would be a command line switch to apply AFNI's PEPOLAR SDC workflow rather than FSL's. I noticed that in sdcflows, both are currently implemented. So this shouldn't be too complex:
https://github.com/nipreps/sdcflows/blob/master/sdcflows/workflows/fit/pepolar.py

It seems like this is probably correct. Have you manually tried TOPUP and 3dqwarp on your dataset to verify that it's the quality of the tools on your data, and not some issue with how we're preparing the data for the tools?

[oesteban]

The reasons to abandon 3dQwarp:

• inability to use images with different total readout times
• inability to estimate when using orthogonal directions
• lack of exposure of BSpline basis
• lack of specific regularization (which is related to all above).

Personally, I'd love to avoid FSL in general (for the commercial limitations) and topup in particular (it's a rather opaque tool).

Unfortunately, field map estimation and SDC are pretty hairy problems. My lack of time is making SDCFlows accumulate a number of use cases like yours where it fails dramatically.

I'll have a deep look into this tomorrow CEST. THANKS much for your patience

[oesteban]

2. More importantly, I think it would be great if there would be a command line switch to apply AFNI's PEPOLAR SDC workflow rather than FSL's. I noticed that in sdcflows, both are currently implemented. So this shouldn't be too complex:
https://github.com/nipreps/sdcflows/blob/master/sdcflows/workflows/fit/pepolar.py

This would go against the philosophy of fMRIPrep. The first necessary step here is to determine where SDCflows is failing. At the most abstract level there are two main possibilities:

1. The field is not properly estimated (and hence, it doesn't really matter what goes afterward)
2. The field is okay, but using it to unwarp the data results in problematic results.

Then, within option (1) there could be several reasons for it not to work:

• Misconfiguration: not passing the appropriate settings to SDCFlows, and/or automated preparation steps within SDCFlows fail.
• Poor performance of topup
• Bug in the generalization of topup results to reusable field representation

Within option (2) there are other problems:

• Misconfiguration: the effective echo spacing is wrong (and hence the fieldmap is wrongly scaled)
• Registration between the field reference and the actual EPI dataset to be corrected fails
• Bug in the mapping of the BSpline coefficients into the target EPI
• Bug in the interpolation of the field in EPI space
• Bug in the application of the corresponding displacements field

Why the former implementation with 3dQwarp was less brittle? Beyond the limitations on my previous comment, the main reason is that the susceptibility distortion workflows have exploded in their complexity, and, at the same time, have been pulled out from fMRIPrep to work more generally (e.g., with diffusion data).

For example, the implementation with AFNI does not give any flexibility on the acquisition parameters, and the effective echo spacing / total readout time parameter had to be the same on all the EPIs involved.

The new implementation allows:

• Mismatched EES/TRT between EPI images used in estimation (they can even have orthogonal PE directions).
• Mismatched EES/TRT on any target EPI

In addition, the new implementation intends to make the fieldmap more robust across sessions: previously, using a fieldmap estimated on one session to correct EPIs from other sessions was risky. Now, since only the coefficients are mapped without interpolations, there's guarantee that the fieldmap is going to be very very consistent across all EPIs.

I will prioritize SDCFlows so that it gets more stable ASAP. But, unless 3dQwarp addresses its limitations, we have no other option than to stick with TOPUP.

(at least that's my view on the problem)

[oesteban]

In more succinct words: this is not a regression, this is a new bug and addressing the problem comparatively to fMRIPrep < 20 is IMHO the wrong approach.

[Gilles86]

Okay. It seems like you (obviously) gave this more thought than I did!

One comment though: for the pepolar-approach, if you're using identical EPIs but with opposite phase encodings, it doesn't really matter what ees you provide, right?

And where do you propose we go next? Would it be helpful to have some of my data?

[oesteban]

One comment though: for the pepolar-approach, if you're using identical EPIs but with opposite phase encodings, it doesn't really matter what ees you provide, right?

Partially correct for (1) but incorrect for (2).

For topup, the actual TRT values don't matter. What matters is the ratios between them.

In the case of (2) the TRT is critical with a decent implementation of susceptibility distortion correction. The reason is that you first estimate the fieldmap in Hz. To convert that to a displacements field you need to scale by the TRT and voxel size.

And where do you propose we go next? Would it be helpful to have some of my data?

This would definitely help in two ways. First, because you've hit an instance of a problem with very obvious consequences (easier to debug). And second, because if we manage to spot the problem, you will quickly benefit that the fix is rolled out :)

So, if you don't have any problem with sharing some data privately, I will allocate time next week to debug this. SDCFlows is currently a can of worms and I'm the main blocker to take action.

[Gilles86]

OK, @oesteban . I sent you an email. :)

[oesteban]

Okay, I bear good and bad news.

The good news is that, as in #2835, the fieldmap estimation looks good:

• The field looks reasonable (the red/transparent/blue hue)
• The EPI shown as background does seem to have an acceptable correction (although the temporal lobes are too close to the edge of the FoV and the correction there is harder to assess. (Note please consider a more generous FoV for these blip-up/down EPIs going forward.)

The bad news is I just started testing, so I have no clue of how long this may take. One thing that is suggested by the data is that either the fieldmap estimation has the wrong sign or the phase encoding of your EPIs have. Indeed, your SDC'd images look worse than before correction.

Interestingly, this is exactly where I got with @mgxd's data one year ago. I hope that so much accumulated evidence allows me to catch the problem.

[Gilles86]

ok. thx for the update and working on this. Hope you catch the potential bug!

[oesteban]

Okay, I narrowed it down to a field flip when generalizing topup's coefficients output. I still don't know if this flip on the PE direction is caused by the fact that the first PE that goes into topup is i-, the fact that your images are LAS-oriented (i.e., fMRIPrep may be disregarding the polarity flip of the LR axis), or some interaction between both elements.

I have also looked into @mgxd's data and I believe there could be additional issues because his data are LAS too, but the first PE direction fed into topup is j.

[oesteban]

I just confirmed the flip. I'm not sure though what should be patched - the generalization of topup coefficients, or the generation of an actual displacements field that unwarps the dataset.

If it were the latter, then we should see this problem on GRE and SE fieldmaps too, from time to time. I'll start with topup and generalization to see where I can get.

[oesteban]

@Gilles86 - as discussed in nipreps/sdcflows#278, I believe I've gotten to a point where topup is executed and processed more reliably.

However, I have the impression that either your images in the fmap/ folder or the BOLD runs under func/ have the PE polarity inversed.

If you open the fieldmap EPI with polarity i and one BOLD run of the same PE direction you'll see that the distortions are opposed. If you then open the fieldmap EPI with PE i-, then you see it matches the distortion of the fieldmap.

[Gilles86]

@oesteban – thank you so much for your efforts!

Is there an easy way that I can run a development version of fmriprep to test your work on my data myself?

If you open the fieldmap EPI with polarity i and one BOLD run of the same PE direction you'll see that the distortions are opposed. If you then open the fieldmap EPI with PE i-, then you see it matches the distortion of the fieldmap.

This is not my experience. If I open sub-34_ses-1_task-task_run-1_bold.nii next to sub-34_ses-1_dir-LR_run-1_epi.nii, or sub-34_ses-1_task-task_run-2_bold.nii next to sub-34_ses-1_dir-RL_run-2_epi.nii, they –to me– seem to have clearly opposite distortions.

Maybe you are looking at different files or I'm misunderstanding you somehow? Are you saying that I have used i-encoding where I should have used i--encoding? How would one even see that? Shouldn't the warp field remain the same? It's just the field map that then would have a different sign, no?

[oesteban]

Are you saying that I have used i-encoding where I should have used i--encoding?

Yes, this is what I suspect happened. If the two files have clearly opposite distortions, then their metadata should have opposite PhaseEncodingDirection values.

Shouldn't the warp field remain the same? It's just the field map that then would have a different sign, no?

Actually the opposite: the warp field (as in, displacements field that unwarps the distortion) will have reversed directions between two opposed PE runs. The field map (as in, the deviation from nominal $B_0$ in Hz units) is the same in both cases. The PhaseEncodingDirection metadata just anticipates in which direction the voxels will be displaced along the PE axis.

[Gilles86]

Yes, this is what I suspect happened. If the two files have clearly opposite distortions, then their metadata should have opposite PhaseEncodingDirection values.

Sure. But I that's what I see in my data? Could you give me a very concrete example of where this goes wrong?

[oesteban]

Could you give me a very concrete example of where this goes wrong?

subject 39, that's the one I was testing on and inferred all of them would have the PE confused. However, you are right that for subject 34, if I take the two first examples you gave ( sub-34_ses-1_task-task_run-1_bold.nii and sub-34_ses-1_dir-LR_run-1_epi.nii) they seem to be correct.

could you confirm that subject 39 shows this issue? Meanwhile, I'm going to extend my tests to the rest of subjects you gave me (this far I only used 38 and 39).

[dlevitas]

I've also noticed similar issues in the SDC workflows between fmriprep versions 20 and 21. I've attached two screenshots, one for each subject's run 1 across two fmriprep versions. Version 20.2.0 appears fine, while version 21.0.2 performs quite poorly.

I have spin-echo fieldmaps with 180 degree flipped PEDs: PA (j) and AP (j-). For functional, I have BOLD and a corresponding SBRef, each with PED of (j-). @oesteban hypothesizes in nipreps/sdcflows#278 that "the fieldmap estimated by topup is only correct in polarity if the input EPIs are RA{S,I} or AR{S,I} oriented." My fmap, bold, and sbref data are all in RPI orientation, according to AFNI's 3dinfo command.

I'd be happy to share my data and fmriprep html reports if that would be helpful.

[Gilles86]

Should this issue now be fixed in the 22.1.1 version of fmriprep?

[effigies]

This should be resolved in the 23.2.0a1 release.