I ran a subject (2x2x2 mm resolution, b0 vols + 60-dir b=700 AP+ 60-dir b=700 PA+ 60-dir b=3000 AP+60-dir b=3000 PA, generating after the dwiproproc recombination a b0+60 vol b=700+60 vol b=3000 multishell mif image) with either the dhollander or the msmt_5tt algorithm:
The top image is dhollander output, and the bottom msmt. It is apparent that CC/midline structures such as splenium are fatter in dhol than in msmt and overall dhol looks cleaner, however msmt ODFs seem to extend/branch out more laterally/peripherally. In other words, dhol appearance is more ‘centripetal’ vs msmt ‘centrifugal’. Is this the expected behavior?
Three more questions regarding this:
How is dhollander vs tournier on single shell data, any takes?
If one omits the -lmax flag with dwi2response, either ‘tournier’, ‘dhollander’, or ‘msmt_5tt’, is the default set based on the data (no of volumes, on teh set herein lmax would be 8) therefore one can omit it, or the default lmax is set on a fixed, lower value?
Since one explanation offered for the good dhollander performance vs msmt was the need for the latter for a good T1 registration, are there preferences among mrtrix experts for a verified, well-performing dwi-T1 registration method that we could implement (without guarantees or strings), or is it the wild west? We tried FLIRT/BBRegister and did not see big differences by eye-balling (both methods gererate areas with visible, if small, mismatches), but wanted to know you opinion for more ‘advanced’ (nonlinear) methods (eg. one recently published, regseg) ?
I’ve split off your original reply as a new topic, as it was a bit hidden in an existing topic (that was about another topic ). I hope you don’t mind.
I’ll get back to your questions tomorrow; I’ve got some good answers, but it may take some time to write them up decently… I’ve also got a figure or two to help illustrate some things for you.
In the mean time, a few quick remarks on the commands you’ve posted; as I’m spotting an error in there and I’ve got another piece of advise. It may be best to test (again) using these points, so as to make sure we’re still talking about the same findings:
The -lmax 8,0,0 in dwi2response (both dhollander and msmt_5tt) is wrong: the -lmax option for dwi2response specifies lmax’es for each b-value rather than for each tissue type. For these algorithms, the lmax’es that you specify only apply to the WM response, as the GM and CSF ones are inherently always isotropic (lmax=0) for all b-values. So in your case, this should have been something like -lmax 0,8,8 instead, reflecting the WM lmax for b=0, b=700 and b=3000. However, it’s even better to not supply the lmax option at all! The lmax is now set automatically to 10 for you for all non-b=0 shells of the WM response. Don’t worry about your 60 directions not being enough for that or something: the response is estimated from at least hundreds of voxels, so lmax=10 (and in principle even way beyond) is easily supported for the response functions. So long story short: ditch the lmax option to dwi2response.
I spotted a suspicious -mask mask_dilated.mif option to dwi2response there. The suspicious thing (for me) being the “dilated” there. This may make sense to dwi2fod, so you’re covering a safe extra margin when computing FODs, but for dwi2response this isn’t needed. In fact, this may even be detrimental there: as no useful response functions will ever be obtained from that margin beyond the brain and outer CSF. In extra-fact, dwi2response dhollander even on purpose erodes the brain mask: the definition/intuition of a “safety margin” for the purpose of response functions works quite the other way around: you’ve got loads of voxels in a whole brain mask, so you can be extra-picky about choosing only the best ones. It makes sense to avoid dangerous areas (like the very outer few voxels of a potentially not-entirely-accurate mask) at any cost. Your data looks like relatively normal developed brains, in which case dwi2response will happily compute a decent (enough) brain mask itself, and erode it for you. So long story short for this one: ditch the -mask option to dwi2response as well.
Your dwi2fod call is fine (and correct), but you can safely ditch the -lmax option here too. It’ll manage the lmax’es safely on its own for you.
So I would recommend re-running both dwi2response calls (and then dwi2fod of course as well) without the -mask and -lmax options, and let them do their thing that way.
Maybe a small additional request, so I can help you better (tomorrow when I work on a reply for you): could you zoom into a region (somewhere, anywhere around the cortex). The current screen shots are very hard to look at for those regions, as the resolution (of the screen shot) is quite low. It’s impossible to see what those FODs look like around these areas.
Finally, definitely also provide a copy-paste of the text (numbers) in each response function file (they’re just text files), e.g. as was done in this post. They typically explain much more about your response functions than the FODs (indirectly) do. The interactions in multi-tissue CSD can get quite complicated between the 3 tissue types at times, depending on certain properties of the actual response functions.
Ok, that was already quite some text. I hope my actual reply tomorrow won’t be much longer…
Here are two figures/insets of the left frontal cortex and wm, along with interhemispheric fissure (right margin) and part of ccal (right lower corner). Top figure is dhollander, bottom msmt (detail 3, scale 3, no interpolation, to accentuate differences) over the act_vis.mif structural.
Again, dhollander FODs seem to be better confined to structural wm (with exceptions, such as some portions rich in CSF), while enriched in CCallosum, in contrast, msmt ones are more ‘diffuse’, often with larger lobes at the gmwm junction or in the gm.
Please add these questions to my previous ones:
1/retake: If the default msmt_5tt and dhollander lmax for wm is 10, which is the default for tournier?
2. If we are not to worry for a default lmax of 10 when no or vols/grads would indicate an lmax of 8 for the RF, should this lack of worry extend to the case when the no of images is *much" lower (~30)? This is to cover most bases in case of a pipeline that may be used on image sets with widely different no of gradient directions.
3. when we should expect the dhollander paper out?
4. re dwi/structural registration: in addition to the larger discussion involving many aspects of this kind of registration, is there a tool other than FLIRT/bbregister that has successfully been used in your hands AFTER topup/eddy application (whether or not strictly for b=0/T1 transform)?
Thank you for taking the time to answer all these,
Great stuff, thanks for posting all that! It beautifully reproduces my findings up to now on a lot of datasets, as well as data from a few of other people around here who’ve also wondered about the differences. Let us take a look at how we can perfectly explain everything we see in your results. I did a quick plot of your reported response functions’ b-value dependent behaviour, which makes things a bit more visual (I always recommend taking a look at things like this):
WM-GM-CSF = blue-green-red (as in my abstracts, and also in Ben’s original MSMT-CSD paper), dhollander is full lines, msmt_5tt is dashed lines; horizontal axis is b-value (just 3 discrete points of course; your b-values); vertical axis is amplitude of each response (your numbers, but in general these are arbitrary units).
Ok, so, observations:
WM responses are as good as the same
CSF response function of msmt_5tt has a drastically lower b=0 intensity (compared to dhollander) and minor increase in b=700 intensity
GM response function of msmt_5tt has a significantly larger b=0 intensity (compared to dhollander) and minor increase in b=700 intensity.
What can have impacted msmt_5tt to explain the different CSF and GM responses?
Mis-registration of T1 vs DWI data: since CSF has the largest b=0 of all tissues by a wide margin, any (even subtle) mis-registrations will cause selection of “CSF” voxels that are non-CSF tissues, resulting in what can only become a severe drop in b=0 versus what it should have been. msmt_5tt has one additional form of protection, which is an additional FA threshold for GM and CSF (to stay below). This protects it a little bit against these problems where CSF interfaces with WM (high FA) specifically. However, it doesn’t protect against this issue at CSF-GM interfaces. And even at CSF-WM interfaces, it’s very limited in doing so, since “WM” FA drops very easily to quite low levels when even partial voluming only a bit with CSF. So we may still see (and get) voxels there with a low enough FA to make the cut, yet still significantly partial volumed by WM.
Mis-segmentation of tissues that aren’t 100% CSF as still very much close to “100% CSF”. Number one culprit to significantly cause this even in perfectly healthy brains: the choroid plexus structures in the ventricles. While these are sometimes barely decently visible in a T1 image, they are very real structures sitting there and changing the b=0 more than enough to be taken seriously. I note in your dataset, since you’ve also shown the 5tt image in the background (I presume) that in fact some of those got segmented as GM (which is also not what they are, see further on for GM response issues), which at least kept them partially away from the CSF selected voxels (one would hope).
Mis-registration of T1 vs DWI data: you effectively run the opposite risk compared to the CSF scenario here. You could first assume that this may “average” a bit out in a way: GM is mostly the cortex, which is mostly wedged in between WM and CSF (outside the brain) areas, so some random mis-registration may naturally land you with some GM-segmentation still overlapping (in DWI space) with actual GM, but some also with WM and some with CSF. Let’s naively assume 33% for each tissue type ending up in the “GM” segmentation or something like that… you could then naively assume that this saves the day a bit, since the b=0 of GM also sits in between that of WM and CSF. However, 2 important biases come to ruin the day here: first, the CSF has a much much higher b=0 whereas WM only has a moderately smaller b=0 intensity, so an “average kind of” mis-registration already biases towards the higher b=0. Second, there’s that FA threshold for selecting CSF at the end too, which prevents those accidental inclusions of WM voxels. So the bias will almost uniquely be towards the b=0 of CSF, which is indeed drastically higher than it should be (for GM).
Mis-segmentation of tissues that aren’t 100% GM, or not even GM at all: the cortex is not that wide, so the not quite 100% thing happens very easily. However, the “not even GM at all” is the biggest risk here: about any pathology of WM that lowers “fibre density”-like aspects, or even more so myelination, will result in a lower (than WM) T1w intensity. Think for instance WM hyperintensities in ageing or multiple sclerosis. We’ve got heaps of those over here (in our institute): 5ttgen classifies a lot of that stuff as (almost pure even) GM. However, in practice, it’s often a mix of tissue changes that tend towards a higher b=0 (well, it’s T2 hyperintensities anyway, right) then what you’d even see in healthy GM. Furthermore in your case: there’s the choroid plexus again. Shouldn’t be GM, and at your DWI resolution of 2mm isotropic, it’s going to have lots of partial voluming from CSF. However, an overwhelming amount of your GM voxels will still come from the actual cortex, so the choroid plexus isn’t going to be responsible for much problems in your case.
In practice, for your example, I reckon both “mis-registration” scenarios have played out, and to a (much) lesser extent the "mis-segmentation ones. I note for instance also a band of “GM” segmented at the interface of the thalamus and the CSF. I can see how that could’ve come from the T1, although it’s wider then what I’d expect: maybe you had an anisotropic T1w resolution? Or the T1 was blurry otherwise? I’d expect to see a band like that, just not that wide. Again, the CSF’s T2 will overwhelm there, and bias your GM response towards CSF properties.
If you do nail registration as good as realistically possible (combined with the best image quality imaginable), you may be able to reduce msmt_5tt's GM problems. Note they weren’t apparent in the data I presented in that abstract. I’ve also pushed this to high quality data, by using the preprocessed HCP data. I reason that’s quite up there in terms of data and preprocessing quality. See here the results of such a dataset:
Top row is msmt_5tt going through some steps (final column is the final voxel selection in diffusion space), bottom row is dhollander going through the same “main” steps as presented in the abstract. These are the responses that came out here:
Note the WM and GM are basically the same between both algorithms (even though the selected numbers of voxels for those 2 tissues were very different; this is good news in general, it does show that in good conditions the 2 very independent strategies go to a consensus… reassuring!). Note that the CSF still suffers from the problem described in the abstract (and reproduced on the forums here already a couple of times by people, and by your data now as well). If we assume very good registration quality in the HCP case, then segmentation issues it is: note the final CSF voxel selection for msmt_5tt's CSF is still quite generous, and indeed includes the choroid plexus and a lot of risky voxels at the edge of the mask. The dhollander algorithm is much more picky, and goes directly with the diffusion data’s properties (actually looking for a sharp decay as is expected for free diffusion). Look how it nicely avoids the choroid plexus and anything that even barely risks being partial volumed. But the dataset is so large, we can easily afford being this picky: we’ve still got 1824 CSF voxels (and lots of images per b-value) to work with. That should be more than enough to estimate just one average number for each b-value. After all, this is not about a segmentation to begin with: it’s about a very picky selection. (to see the full selection of voxels from dhollander for this HCP dataset in 3D, see the screen shot posted here)
Here’s the responses I got over here on a dataset (with a more realistic spatial resolution), yet with some lesions:
I can’t disclose much of the other results on these datasets at this point, but essentially lesions got selected as GM by 5ttgen, and hence the GM response function selection by msmt_5tt came mostly from those areas. Note the now very expected CSF problem also pops up (of course, as always). There we’re also slight mis-registrations in place here, due to partially uncorrected EPI distortions (motion / eddy currents were corrected though). Again, everything with the expected effects on the response functions.
Next, consequences for your FOD estimation when using your responses. I’ll talk in terms of biases caused by the biased GM and CSF responses that resulted from msmt_5tt. To understand this, it’s useful to consider tissue interfaces, because that’s where the contrast between 2 tissues will go wonky if either (or both) responses aren’t appropriate for what they’re supposed to represent. Here’s all the interfaces:
WM-CSF: troubles here caused by the underestimated CSF b=0, or in general underestimated CSF diffusivity. This scenario is mostly described in the original abstract. Essentially: more “CSF” is needed to represent the b=0 signal using the CSF response that has too little b=0. It’s actually more about b=0 versus the other b-values, and hence about the diffusivity, but the argument has the same format. The opposite holds for the WM response contribution here (even though the culprit is the CSF response, but it’s all about relative balance between both). Result: underestimation of your WM FOD, and it sometimes disappearing entirely for some levels of CSF partial voluming. Your responses and the resulting effect at the WM-CSF boundary is actually one of the most severe I’ve every seen. Don’t underestimate the impact for subsequent applications: e.g. the midbody of the corpus callosum (that is already quite a thin structure relative to DWI resolutions) is impacted severely, on both sides, so it grows significantly “thinner”. This impact quantification of FD, FC and FDC for this structure, and may cause strong biases in any connectome (as a significant number of the left-right connections of the brain go through this crucial bridge). Furthermore, in pathology (e.g. lesions), this may overestimate the free water content, and sometimes severely underestimate the amount of healthy intra-axonal space left. Pretty risky if e.g. surgery is informed by that.
WM-GM: troubles due to overestimated GM b=0. The (or a) simple(st) way to reason about this is to think of the GM response having become a bit more CSF-like. So the behaviour at the cortex becomes a little bit more like the scenario where we would have no GM response (i.e. 2-tissue CSD with WM and CSF). So essentially, some isotropic GM signal that would otherwise be captured by the GM response, is now not. And too much signal ends up in the WM compartment, which is your WM FOD. You can see this, as you’ve described, by larger FODs there… but note that in most cases in your more detailed zoomed screen shots, you can observe that these are messy FODs, or messy additions to FODs that were cleaner in the dhollander case. This is the “fingerprint” of isotropic signal being deconvolved by an anisotropic WM response, while using an SH basis that is cut off at a certain lmax. If you deconvolve perfectly isotropic signal that has just a tiny, tiny bit of noise on it with an anisotropic WM response, you don’t get a nice isotropic FOD, but a sphere full of random peaks. The noise is sharpened, and the SH result goes up and down like crazy. Not a useful addition to your WM FOD in those regions, that is… Again in case of pathology (again e.g. lesions), this may underestimate e.g. gliosis, and have some of its effects messing with your WM FODs. You also end up overestimating AFD then in these regions, and hence underestimate some of the damage in case of a genuine lesion.
GM-CSF: we can’t see this going with just your WM FOD; the GM and CSF outcome maps may show some effects here (in addition to those at the other interfaces described above). But both GM and CSF responses are biased towards each other… so you may see a (false) over-contrast of the GM-CSF spectrum. It’s as if you set the windowing on an image too tight in a way. The balance may also be disturbed dependent on which bias (of both responses) “wins out”. In pathology… well, this becomes a mess.
Finally, note that in general both of these GM and CSF response function biases strictly reduce the range of signals you can fit (even independent of meaning and interpretation). See for instance Figure 5 in this recent abstract: the triangle in the right hand side schematic will strictly shrink due to either of those biases. This schematic has actually already helped me a lot to translate formal observations about the response functions as well as 3-tissue CSD into visual interpretations, and more intuitive explanations. Using that schematic, you can perfectly obtain the above explanations at tissue boundaries under the presence of a response function bias as well. It’s a mighty handy tool (if I say so myself ) for white board discussion marathons (and my boss/supervisor has a big one in his office, so fun times guaranteed ). Also, both the inner workings of my SS3T-CSD as well as even the strategy of the dhollander algorithm can be made very visual on that actual schematic.
Alright, if I continue like this, I may end up with the actual paper almost written. I’ll briefly get to your questions, so I don’t overlook anything you were looking for:
This should now be clear from the above. Especially important for e.g. corpus callosum and fornix.
So from the above: be very careful about this, as a lot of this seems to be “false positive signal” that shouldn’t have ended up in the WM FOD, but rather in the GM compartment. Some of in in your zoomed screen shot is particularly messy, and also quite big compared to the actual nearby deeper WM. Not surprising: if this signal now ends up in the WM FOD, but has to be modelled by the WM response’s lower b=0 value, it’s going to be over-represented even more.
Even though both GM and CSF response biases of msmt_5tt work together in your case to give you this appearance (of dhollander versus the biased msmt_5tt), this isn’t a representation of both of them being more suited towards a certain region. Bigger FODs aren’t always better, as we’ve seen from our WM-GM boundary reasoning.
This is somewhat unrelated to all of the above. First things first: at the moment dhollanderusestournier internally for the final single-fibre WM selection, but dhollander already narrows things down for tournier to a safer selection of WM (both single-fibre as well as crossings and other geometries still). dhollander also determines a more relative number of voxels (compared to WM size and resolution) that scales better in case of pathologies or varying spatial resolutions, and then asks tournier for that. tournier out of the box (default) always goes for 300 voxels, independent of anything. You can set this number yourself of course as well if you run tournier. But in general, the core of the strategy to identify the single-fibre voxels is the same. In practice, for a single-shell response they also tend to give a very similar result. For multi-shell purposes though, tournier ignores lower (and b=0) b-values. dhollander has a sneaky outlier-rejection going on to avoid pesky CSF-partial volumed voxels and other unexpected T2 hyperintensities. This does help it to get slightly better multi-shell (or even b=0 and single-shell) response for WM. But this doesn’t matter for single-shell single-tissue CSD: for that one, only the shape of the response is the crucial thing. No drastic differences here. I’d advise to simply run tournier directly if you need a single-shell single-tissue CSD response. No need to overcomplicate things for that scenario. In the future, there is a chance that I change what dhollander does at this stage. It’ll always be inspired by the general strategy of tournier though, but there’s a few other things of concern going on, and an opportunity to slightly improve the multi-shell (or b=0 + single-shell) WM response; but I need to figure out if I can get it ultimately robust. (robustness is a property of the dhollander algorithm that I’d personally never want to sacrifice)
At the moment, the default for all algorithms will inherently end up being lmax=10. That is, for all WM responses (GM and CSF have lmax=0) and only for the non-b=0 shells of that WM response (again, b=0 has lmax=0). See some experiments and discussions and equally lengthy stand-alone posts on this: here, here and here are some relevant bits, but happily read that entire page of discussion (they’re all bits on the same page). If you’ve got any further questions on that front, don’t hesitate to ask for clarification. But the gist is: even with crappy 12 direction data (imagine that), and even with only 100 voxels for the single-fibre WM response (that’s a lot less than what you’d realistically always get), you’ve still got 1200 measurements, to estimate just 6 parameters (that’s the lmax=10 response model, i.e. lmax=10 for a zonal spherical harmonic basis). On top of that, we’ve got non-negativity constraints and monotonicity constrains for the WM response (if those 1200 measurements weren’t already enough to estimate just 6 parameters). And lmax=10 is certainly more than flexible enough to represent responses up to b=5000 in vivo. So it’s flexible enough, and mighty safe given the massive amounts of data we always have for this task.
As mentioned above, even a perfect registration doesn’t overcome all issues. But for it to at least be as perfect as possible, you’d need it to be sub-voxel accurate (relative to the DWI spatial resolution). And of course, up to that same level of accuracy, be able to overcome the many distortions that plague EPI data. And even then still, msmt_5tt will trip up over mis-segmentations… all while the actual information about the actual responses you’re after is actually right there in the DWI data, rather than in a T1. Note that the msmt_5tt algorithm was never meant to be the ultimate way of doing things; but at the time there was certainly need for at least some automated need to do this. It would otherwise have sounded very arbitrary in the paper(s) that we picked our responses “manually”… it also wouldn’t have been very convenient in practice to do this on cohorts of many tens or hundreds of subjects. msmt_5tt certainly did what it needed to do at the time. However, your question (about registration) remains important: for tasks like ACT, you still want that registration to be as good as it can get of course. It may be worthwhile asking that one specifically again in its own topic (with a well chosen title) to get some input from other peoples’ experience. I’ve seen a few different things here an there on the forum before too.
I should have read all question before I started typing above. So, as mentioned in the other answer: it’s also 10.
Yep, no worries, as mentioned in that other answer above. You’d even get away with 12 directions (or less!). However, the question you’ve got to ask for such a pipeline is what lmax for the CSD step is till appropriate, or even if the CSD step itself is still reasonable. But the response function estimation is (un)surprisingly safe. Also be careful with the results of such a pipeline on data with such differing qualities, as in: don’t compare such results with each other (i.e. based on different qualities relative to each other).
Looks like I almost wrote it in this post now. But other than that: it is my first and foremost priority at the moment, even before any SS3T-CSD shenanigans. Other projects (well, collaborations) do run in parallel, but those are now bootstrapped sufficiently (I hope). Time to focus on this one; lock myself up without internet and with a typewriter.
Again, very useful as a topic of its own; I expect different experiences in the wider MRtrix3 community. I’m always happy to learn what people do for these steps as well!