@thijsdhollander, Interesting to hear about the issues in the current approach. I must admit I haven't had time to have a thorough look into it... I tried running dwi2response on the data I was using to test my estimate_response script, and came across the dreaded no voxel left in mask error... I just assumed it struggled to converge sometimes, but clearly it's a bigger issue than that. I guess this justifies us doing something about it now, then.

So as to how to sort this out, I don't think a simple renaming will do. Clearly the issue is more complex than that, and even changing the commands to output the mask rather than the actual response doesn't really help much, if at all. It still leaves all the issues that we've seen come up time and again on the forum, and the more troubling issues mentioned by @thijsdhollander above. Even if the intention had been to use the current approach to identify a good single-fibre mask that we could then use to estimate the various responses from, it seems clear that it's not up to the task.

One of the problems I can see here is that dwi2response actually needs to implement all of the 4 steps I listed above - there's simply no way to refine the mask without an updated version of the response. It wouldn't make sense for it to not output the response. Yet as implemented there is no freedom to change any of the steps in it. I'd argued before (#189) that we could script the whole thing as-is, since the only piece missing was a sh2fixel command as the orientation estimator. However, I would argue that sh2peaks is actually a better tool for the job, since the amplitude of the peak is lowered by dispersion - this is an effect we try to minimise for AFD analyses, but actually works in our favour here: it biases the results towards low dispersion regions, resulting in a sharper response, as @dchristiaens pointed out.

The script I've proposed here is actually much closer to the original Tax et al. approach. The only differences are that it uses a much more sensible (and stable!) initialisation, and use the top 300 voxels rather than a hard user-defined threshold. With these two changes we end up with something that is much faster and robust, using tools we've already written. And as I've stated above, I think operating on peak ratios rather than segmented FOD lobe volume ratios is actually the right thing to do.

So what I'm clearly coming to here is that I think the current dwi2response command needs to be replaced... I do think it needs to remain available, at least for a while, to ease the transition, but we need to provide a solution that works a lot better than that.

Clearly before we do something like this we need to come to a consensus about the need to do this, and whether the approach I'm proposing is the right way to go. Personally, I think even if this script turns out not to be completely optimal, it has the desirable property of being fairly robust (I can't think of any way it would fail to run and provide a useful set of response coefficients). I've tried it with very low anisotropy, very low b-value neonatal data (b=500, 47 DW directions - really very little anisotropy in these data), and it works fine there (although I did have to remove the mask refinement step, which made it run slower - I'm looking into that right now). It clearly works just fine in standard adult data.

So, as always, opinions welcome...

However, I would argue that sh2peaks is actually a better tool for the job, since the amplitude of the peak is lowered by dispersion - this is an effect we try to minimise for AFD analyses, but actually works in our favour here: it biases the results towards low dispersion regions, resulting in a sharper response

dwi2response does try to reject 'single-fibre' voxels with too much dispersion; that's one addition I made on top of the publication. Whether it's practically better or worse than just using peak amplitude ratio as a surrogate for both density and dispersion: I'm sure I tested it at the time and chose the full segmentation approach, but I don't remember the numbers.

I also have a sneaking suspicion that when using the sh2peaks approach, it's not actually voxels with a peak amplitude ratio above some threshold that are being selected, but those voxels where a second peak is simply not found based on the 60 starting directions (from memory Chantal was using a ratio of like 0.01). So it's probably less stable when doing a per-voxel selection/rejection approach compared to the FMLS segmenter. In the context of top-300-voxels, you may end up with an undefined solution if you have > 300 voxels with no second peak identified - this will depend on whether your 'single-fibre-ness' metric takes into account first peak amplitude in the absence of a second peak.

The script I've proposed here is actually much closer to the original Tax et al. approach. The only differences are that it uses a much more sensible (and stable!) initialisation, and use the top 300 voxels rather than a hard user-defined threshold.

I'd disagree: she explicitly justified starting with a too-fat response rather than a too-thin response, and uses a per-voxel data-driven threshold rather than a # of voxels threshold. So the current dwi2response is definitely closer to that publication than this script. Not that that matters in the slightest: the whole point of this process is going with what works, not what was published.

Just to comment on a few of the issues raised by @Lestropie:

• Happy to switch to a fixel-based approach if we can verify that it performs better. I'm not married to the sh2peaks approach, it just seems simple and satisfies the requirements for the job.
• Good point about what happens when there is no second peak. There's no threshold in sh2peaks by default, it'll just spit out all the peaks it finds. But in low lmax reconstructions (i.e. the first pass), the voxels we're interested in will typically be perfectly unimodal, so there simply isn't a second peak. The default behaviour in sh2peaks is to fill the orientations for these with NaNs... I've just added an option to sh2peaks to specify the value for no peaks, so we can write a zero in those cases. The actual peak ratio is subsequently computed as |peak1|² / ( |peak2|² + 1), so there is no divide by zero possible. Removing NaNs from the output fixes the issues I was seeing in the neonatal data, since the initial SF mask basically excluded all the best voxels...
• As you see from the equation for the peak ratio, if the second peak is zero, it favours high-amplitude peaks, which I think is the right thing to do here.

So I've made some changes, which I've just committed. Just to highlight the main differences:

• Mark non-existent peaks as zero (as discussed above)
• initial FOD estimation with 'sharp', but lmax=4 response (i.e. [1 -1 1]), reverting to full lmax (default or as requested) for subsequent iterations. This is more stable, but also has the added bonus of speeding up the first iteration substantially.
• initial mask is refined after first iteration by selecting the top 3000 voxels (i.e. 10 times the number needed) and dilating by one voxel to include neighbours. This speeds up subsequent iterations drastically, and pretty much all the voxels that would otherwise have been included in the single-fibre mask are present in this mask, so the bias introduced by this should be minimal, if any. In my data, it was only for the b=500 neonatal data where there were any differences, and even then, the differences were tiny.
• I've also added a few rudimentary options. This is just for convenience during testing, of course we'll switch to Python eventually...

Let me know when/if it falls over...

OK, just figured there is an issue with the peak ratio calculation... I'll have a look into it later.

OK, I've had a go at improving the criterion by which fibres are assessed for inclusion in the single-fibre mask. I've settled on this formula:

  sqrt(|peak1|) * (1 - |peak2| / |peak1|)^2

Rationale being that we want to remain scale invariant: here scaling all the values by some multiplier will still give a map that looks identical (just scaled by some global scale factor), so that thresholding by number of voxels gives the same mask. We also want to promote voxels with higher main peak amplitude, but still allow lower amplitude peaks to be included if they turn out to be highly suitable: that's the sqrt(|peak1|) term - I decided to use a sqrt() here since it softens the bias towards the highest peak amplitudes, but also maintains the scaling-invariance of the map, so no fragile threshold values required (I also tried a log(|peak1|) there, but that's too soft, still allowed a lot of nonsense noisy voxels to have high values). We also clearly want to promote voxels where peak1 is much larger than peak2, but without creating problems when peak2 is zero or near-zero, as Rob pointed out. That's the second term: (1 - |peak2| / |peak1|)^2, which is 1 if peak2 is zero, and rapidly reduces to zero as peak2 approaches peak1. So with that, things seem pretty stable, and I get a decent-looking response within ~15s on a routine dataset. Also seems to work fine at various b-values on neonatal data, which is more challenging.

Everyone, feel free to have a go at breaking it...

Small update since you're playing with this:

I'm currently working on #419 and #422. Part of the changes will be a single dwi2response script, that includes this approach (reimplemented in the Python libraries), the Tax approach (reimplemented in Python), and others. So this branch will likely be deleted...

Sounds like a great plan! Go for it. We could also merge msdwi2fod into it eventually...?

I'll assume you mean msdwi2response rather than msdwi2fod...

Yes, that's the plan. I want to duplicate the current msdwi2response functionality as one algorithm, then have another that determines tissue masks based on the ACT 5TT image (like how Ben did it in the MSMT paper). Eventually when his matrix factorization is done that can be added also.

Even ye olde FA threshold approach gets its own algorithm in the script :-P

😁 yes, I did mean msdwi2response...

Sounds good!

Closing since implemented properly by @Lestropie in commit fece64a...