Ok, so having reviewed the code here, as I understand it:

1. Problem is motion between RPE pair (or other separate SE EPI input) and DWI data to be corrected for EPI distortions. We should be quite careful what we do/change with/about this scenario, as I reckon this is probably the most common use case nowadays (i.e. acquiring an RPE pair for this purpose).
2. Proposed implemented strategy is looking for the first b=0 image in the DWI data, and add this to the start of the topup input image. As (I think) I understand it, topup also deals with motion itself, and expresses its output relative to the first input image. So that should then make the estimated EPI distortion appropriate relative to that first b=0 volume that came from the DWI data. That in and of itself sounds good to me. Only uncertainty I've got here: does topup deal correctly with the different/changed balance between PE images? What I mean is: say I've got a simple RPE pair, as an AP and PA image and the DWI dataset is AP. The strategy will add an AP image from the diffusion data to the pair, so there's now two AP images and one PA image. Is this a problem? I suppose this is dealt with (correctly) somewhere?
3. Concerns/limitations:
   a. Scenarios that aren't eligible for correction (any more): different echo times between SE EPI images and DWI dataset, and DWI datasets without b=0 images. More generally, the scenario where the SE EPI images don't have an image in the DWI dataset with corresponding contrast, because one of the latter has to be added to the former set. Does this happen a lot in practice? And if it does, shouldn't for instance the old strategy be applied, where the user reasons motion is negligible or something? That's what is the current strategy any way, so at the moment this is already assumed (but it may not per se be good of course).
   b. The strategy may still not be perfect, because the first b=0 volume in the DWI dataset may not be the first volume of the DWI dataset, so not per se the best reference volume, as motion is dealt with in eddy relative to the first volume in the DWI dataset (right?).

So well, I'm not entirely sure about this one. I'm not against the strategy; it does make sense to deal with it like this. But my biggest worry would be how common those non-handled cases are (a question to which I don't really know the anwer...). Also: how big is the impact of (a bit of) motion between the SE EPI images relative to the DWI data on the EPI distortion correction? EPI distortions are spatially still "relatively" smooth, right? Everything is relative of course, but I guess I'm trying to get a grip on how big the original problem really is versus the (potential) concerns that the proposed solution may bring. This is not a negative assessment of the proposed solution: I'm just unsure about a few things that are key to making a good assessment of the proposed solution.

@jdtournier and @dchristiaens , I've added you as (proposed) reviewers. Not to put time pressure on you guys (I understand if you're busy), but just to hear what are your ideas on this, and whether my semi-assessment makes sense (or to add nuance to it where needed). If too many people are still undecided about this and we detect it may need more discussion (entirely justified if the case), then I would personally hold off packing it in RC2. If you both are confident this solution is safe to proceed with already, then by all means, allow it to be merged in.

I'm not sure I'm the best person to review this - I don't get to process data very often these days, and when I do, it's typically neonatal, so heavily non-standard...

That said, a few comments:

does topup deal correctly with the different/changed balance between PE images? What I mean is: say I've got a simple RPE pair, as an AP and PA image and the DWI dataset is AP. The strategy will add an AP image from the diffusion data to the pair, so there's now two AP images and one PA image.

This has come up in discussion before, and I'm not sure what the answer is here. I'd expect it might introduce some bias since the cost function is a simple least-squares on the difference, but I'm not sure I have any particular insights into how much of an impact this might have. My best guess is: not a huge amount - but that's based on gut feeling rather than hard facts...

Scenarios that aren't eligible for correction (any more): different echo times between SE EPI images and DWI dataset, and DWI datasets without b=0 images.

Were these ever eligible for correction? I don't see that our scripts ever allowed these kinds of cases...

The strategy may still not be perfect, because the first b=0 volume in the DWI dataset may not be the first volume of the DWI dataset, so not per se the best reference volume, as motion is dealt with in eddy relative to the first volume in the DWI dataset (right?).

I guess this is what this comment is about... I agree it would make sense to at least verify that the first volume is the same b=0 image as was used in the TOPUP stage. I assume this is in the works (?).

Only uncertainty I've got here: does topup deal correctly with the different/changed balance between PE images?

I deal with this: code from here. If the SE-EPI data are "phase-encode balanced", i.e. the sum of the phase encode directions of all volumes is a vector [0, 0, 0], then I erase a volume from the SE-EPI image as I insert the first b=0 volume, in order to preserve that "balanced-ness". If this is not the case in the input image, then I just say "bugger it" and insert the DWI b=0 volume as-is.

Exactly how topup behaves in the latter scenario, I honestly don't know, will depend on the particulars of the implementation. Could grab all the b=0 images from a DWI series and a single reversed phase-encode volume and feed them in and see what happens. But I'm going off the premise that "balanced" = "good", and that failing to preserve such would probably attract criticism.

... different echo times between SE EPI images ...

Not aware of anybody doing this presently, but it's actually how I advocated that it should be done at the time, and wanted to do an ISMRM abstract on it, so would have liked to preserve compatibility with such. But nowadays there's other alternatives that would be more worthwhile pursuing.

... and DWI datasets without b=0 images.

Again, haven't encountered such myself. The proposed code will error out on this. One possibility would be to add a volume from the SE-EPI series to the start of the DWIs before running eddy, then erase it afterwards, on the assumption that eddy will appropriately deal with alignment between shells; but given that the method is called "Post-eddy alignment of shells", it's entirely possible that this wouldn't solve the problem, depending on if the realignment re-applies the inhomogeneity field. So I'm not sure I want to further complicate the script to support this, it's getting pretty hard to follow as it is.

And if it does, shouldn't for instance the old strategy be applied, where the user reasons motion is negligible or something?

Maybe. At least with a warning. Though I worry that warnings are going unnoticed; should fix that.

The strategy may still not be perfect, because the first b=0 volume in the DWI dataset may not be the first volume of the DWI dataset, so not per se the best reference volume, as motion is dealt with in eddy relative to the first volume in the DWI dataset (right?).

I guess this is what this comment is about... I agree it would make sense to at least verify that the first volume is the same b=0 image as was used in the TOPUP stage. I assume this is in the works (?).

My interpretation is that topup defines the estimated field with respect to the first volume in its input, and eddy assumes that this is "intrinsically" aligned to the first volume in its input. So yes, if there's motion between the first volume and the first b=0 volume, you're screwed; hence the floating TODO. Again not an acquisition I've ever seen...

Solution is similar to a previous comment: re-arrange the volumes as they are input to eddy, then re-arrange again afterwards. It's just, again, more complexity on a script that's already doing my head in.

Also: how big is the impact of (a bit of) motion between the SE EPI images relative to the DWI data on the EPI distortion correction?

It was enough to cause a bit of concern on the FSL mailing list, and enough for me to get contacted and harassed about it. Certainly MCRI were concerned since kids have a lot of motion, and inter-protocol motion can often be bigger than within-protocol (i.e. "The scanners's stopped, I can get off this pressure point now"), even with correct instructions from radiographers.

EPI distortions are spatially still "relatively" smooth, right?

This is still unclear. With the method I worked on at the time, I was picking up some high-frequency stuff at fluid-tissue interfaces, but it's not clear whether or not it was over-fitting. Looks like topup's resolution maxes out at 4mm, but there can be sharp transitions between those points.

Scenario where first volume of input DWIs is not a b=0 is now handled & tested.

Also permitted processing where the input DWI contains no b=0 volumes to proceed with a warning rather than erroring out.

I've also tested a case where the DWI and SE-EPI images were defined on a different voxel grid, and that completed without issue.

... and DWI datasets without b=0 images.

This crashes eddy.

Quickly settling the previous concerns that I had about quite exotic scenarios: all good, we probably shouldn't worry too much about those. But I'm starting to get cold feet about the actual main change of this pull request (replacing an SE EPI image by one from the DWI data):

Ok, I've done a few quick tests with topup to try and (de)mystify some things. The result of this is some concern... Here goes a quick overview:

Data from a collaborator, so probably can't share the raw data, but the results still make a point. The data that I got came with a separate pair of AP and PA images, and the DWI data was AP. The DWI data itself started with a b=0 image, and 3 others (so total of 4 b=0 images) were spread through the acquisition of a single b=2800 shell. This is quite representative for what I think is the most common scenario nowadays. Furthermore, there was only very little motion in this case; a pretty "safe" case so as to say. I ran topup with:

1. Just the AP and PA pair, as one would "usually" do, both inherently via dwipreproc as on master, as well as via the example scenario in the FSL documentation
2. Replicated the AP image 5 times (i.e. exact copies), followed by the PA image. Goal: test the balance question we had directly. I.e., "balance due to numbers of volumes"
3. Multiplied the AP image intensities by a factor 5 (but just 1 copy as in the usual scenario), followed by the PA image. Goal: test a difference balance question, i.e., "balance due to intensity (difference)" and just see how dependent topup is on normalised intensities
4. The AP image, followed by the 4 AP images from the DWI data, followed by the PA image. That's the same balance (due to numbers of volumes) as experiment 2, but then with actual different b=0 images. This experiment (and what it was supposed to check; i.e. just add some different images and a tiny bit of motion versus experiment 2) was flawed in its setup though, since there was another difference that I hadn't anticipated when setting it up. See below.

So, results:

1. Field map looks normal (i.e. as expected), motion parameters file shows no motion for the first volume (which is just hard coded, the first volume isn't actually checked for motion relative to itself) and a tiny bit for the second (PA) volume. All makes sense.
2. Produces very close to the same field map as experiment 1. Motion parameters shows no motion for the first volume, but (very) little motion for all others, even though the next 4 volumes are exact copies of the first volume. A bit funny, but not problematic. This confirms that topup deals with "balance due to numbers of volumes" correctly. I was first a bit more naive (in my head) about what topup does in practice, but reading a bit about it, this actually makes sense. This could lead us to feel safe about bringing in all b=0 volumes from the DWI data, so we don't need to throw an SE EPI one out and we can use all data. But we're not there yet...
3. Same thing: produces very close to the same field map as experiment 1. I was a bit surprised here, because I (too) would've thought topup just used a least squares metric. But so: it deals well with gross global intensity differences. Motion parameters in line with experiment 1 too.
4. So given the previous experiments' findings, this one should've worked well too. But it didn't: the field map was definitely obviously different. Then I noticed that the b=0 images from the DWI data we're somewhere around a factor 1.8 ~ 2 different in intensity range compared to the SE EPI pair. But even that shouldn't be an issue, if we combine our successful conclusions from experiments 2 and 3. But then I noticed the contrast was also not exactly the same. I don't have all the acquisition parameters (nor the expertise to judge this correctly), but: even though the contrast was different, the distortions definitely were not different: all those AP images from the DWI data looked very similar in shape (due to distortions) to the AP image from "the pair". So the pair does, as far as I can tell, a very good job at estimating an appropriate field map for the DWI data, but mixing in b=0 image(s) from the DWI data to compute the field map becomes a disaster...

Just to at least give some visual insight and the confidence that the distortions were the same (and motion was very limited), yet the contrast and intensity level wasn't, here's those input images:

AP from "the pair":

4 b=0 AP images from DWI data:

PA from "the pair":

All the above are shown using the same intensity range / windowing, bottom of the range being 0. Note the DWI b=0's have lower maximum intensity (e.g. in CSF), but actually slightly higher minimum intensities (e.g. in WM). So the DWI b=0's show less contrast. Note that (and I've checked this from many different angles) the distortions between the AP from "the pair" versus the APs from the DWI data are very (very, very) similar.

So this is where I'm getting cold feet... as feeding this into the proposed dwipreproc in this pull request would go quite wrong (I haven't tested it explicitly, but it follows almost directly from the above experiments now). This opens up maybe a "bigger" question: how much (and what) should we aspire to automate in dwipreproc, and of course by doing so, take quite some responsibility for. At the moment (master) the value of dwipreproc is mostly to make interfacing with FSL's tools much easier, and probably most importantly (for us) in line with MRtrix-style of interfacing. We essentially wrap FSL's tools here to not break a user's MRtrix-style environment and experience (as well as concatenate a few of those tools, and offer more direct ways for the most common "normal" scenarios). All of these things are of great value, and very (very) much justify the existence of dwipreproc. But what this pull request is proposing does go fundamentally beyond: it automatically makes certain ("design") kinds of decisions for the user. Of course we can document all that, but it would still happen automatically. I'm starting to wonder if that's not (far) too dangerous...?

To turn this the other way around: if a user would want to do what this pull request implements and offers directly (e.g. let's consider the example of replacing the AP image of their pair by the first b=0 image from the DWI data), it's pretty simple for them to do this (consciously) by mrconvert ... -coord 3 0 piping into mrcat ... with their PA image, right? And at the moment for even just their AP-PA pair out of the box, they need to do that mrcat step already anyway... so concatenating something else and creating this kind of customised setup is pretty simple with the existing dwipreproc. Wouldn't it be safer to leave this design decision with them? What we're proposing here now goes even beyond the FSL manual's example scenarios (so isn't even explicitly "advised" there, even though I admit it can make good sense).

So well... I don't know, but I'm starting to feel that it might be a better idea to stick with a (very) convenient wrapper-style dwipreproc as it is now, and not automate too many things that could be (too) debatable and should maybe be left to the user to decide upon given their own scenario at hand.

Just before I leave the office (where I'm starting to live almost... erm... 🆘), something that just came to mind: a difference that I think may have been there between the acquisition of the AP-PA pair versus the DWI data may have been multi-band. The DWI data was certainly acquired multiband, but I'm not sure if the AP-PA pair also was... I'm far from an expert in how this may affect the image contrast, so I could be talking complete nonsense here in implicitly suggesting that this explains the contrast difference. But it's certainly not an uncommon case nowadays, since many labs are starting to implement multiband for their DWI acquisition.

I was going to say that there'd need to be something fundamentally different in the acquisition; ourselves we literally use the diffusion sequence with one b-value of 0, so it's even twice-refocused. But multi-band would indeed have a TR effect; single volume acquisitions would almost certainly not be multi-band. The distortions might be different due to differing bandwidth (sounds like the recommendation is to disable GRAPPA with multi-band); this should be captured in the header entries if things are being imported & working as they should, but certainly won't work if using the older -rpe_* options.

Could activate this feature only when a command-line option is provided; that'd mean it could be used by those for which it's applicable, but this particular solution may not remain in place down the track, so far from ideal.

We need to carefully consider dwipreproc more generally (again): there's all sorts of other things I could capture in header entries and use, and I could try the alternative solution to #874, but the "convenience" of being able to use the older -rpe_* options may be opening the door for bad pre-processing. And trying to handle more and more use cases may end up requiring more time and effort than I'm willing to put into wrapping somebody else's methods, but removing functionality that was once provided never goes down well.

Going to have to sleep on this one. Good catch...

My instinct is to close this and omit from #1052. Any discussion regarding this code / solution specifically can happen here; discussion of alternative possible solution(s) / scope of dwipreproc in general to happen in #874.

My instinct is to close this and omit from #1052.

Agreed.

The core of the "issue" (motion between DWI series and e.g. RPE single volume acquisitions) has definitely been brought to the attention of the FSL developers as well... so while discussion can continue of course, I'd also sit back a bit and see what they come up with; both as part of their software as well as recommendations. As to our own users: for those that face a scenario where this issue is relevant, I'd say it's still pretty easy for them to take the necessary steps (even via the current dwipreproc) to deal with it if they wish (and are themselves certain that their data allows for it correctly).

Closing this off without merge to clean up the PR list. Not deleting the branch in case there's some parts of the implementation that may be useful when re-attempting #874.