Fibre orientation wrt B0

[lottevdvoort]

@ThijsDhollander

Good bits to check is the “end” bits of bundles towards the cortex and even in the cortical gray matter: those are the areas where not only false positive lobes and fixels might start to appear, but more importantly, where they then quickly appear in larger numbers. You might for instance see voxels there that suddenly have 4, 5, or even many more fixels. Those are certainly false positives, so if you see many of those, you might want to try to increase thresholds.

I do indeed see a sudden increase of the number of fixels at the very end of the corticospinal tract and the interface with gray matter. However, I am using msmt csd so I run fod2fixel only on the wm fods. Do you mean that the fixels you’re talking about in gray matter would still be visible for the wm fods as false positives?

@ThijsDhollander

So as mentioned just above here: you can first extract the relevant metric using either -afd or -peak when performing fod2fixel . Then you can use fixel2voxel : provide it with the relevant data file as input (first argument), and choose the relevant operation as second argument. In this case, that would be max : the maximum of any fixel-wise metric across all fixels in each voxel.

Thanks! That seems to work although the maximum afd in this voxel image is different form the one in my entire fixel image when I colour code the latter by afd.mif. Aren’t those supposed to be the same since max extracts the largest afd for every voxel?

@ThijsDhollander

split_data on the other hand would give you the value for all fixels. That might then actually make it harder again to get the one (“largest”) you’d want to show.

The reason why I initially wanted to use the split_data operation is because I want to make a similar afd map for the second and third largest fixel as well. I thought this would be possible through

fixel2voxel fixels/afd.mif split_data -number 1 peak1_afd.mif

and then by changing -number to 2 and extracting -coord 3 1 I could get the afd for the second largest fixel. The output from both commands (the one you suggested with max operation and the one I sugegsted with split_data) look similar but differ for some voxels. Could you tell me what the exact difference is between these two commands? Both however show the same maximum afd value which differs form my “total” afd value as mentioned above.

Hope this makes sense and thanks for your help!

[ThijsDhollander]

@lottevdvoort

I do indeed see a sudden increase of the number of fixels at the very end of the corticospinal tract and the interface with gray matter. However, I am using msmt csd so I run fod2fixel only on the wm fods. Do you mean that the fixels you’re talking about in gray matter would still be visible for the wm fods as false positives?

It’s worded a bit confusingly (or my brain is a bit clouded at this late hour )… but so basically, if you see many (3+) fixels in the cortical ribbon, I’d argue several of those are false positives. Increasing the threshold to segment the fixels there, that should be reduced to 2 or even just 1 in many places. However, the possibility to achieve this without in turn losing a lot of true positive fixels in other places, strongly depends on a 3-tissue CSD model, i.e. crucially on the presence of a GM compartment in the model. This compartment then “filters out” the GM signal from the WM FOD, or putting it differently: it “cleans up” the WM FOD, so it doesn’t have so many false positive peaks.

Having said that, this depends also on your data quality and how you went about modelling it of course. I think I didn’t ask that bit yet above though: what is/are your b-value(s), and how many gradient directions do you have for each b-value (and how many b=0 images do you have)?

@lottevdvoort

Thanks! That seems to work although the maximum afd in this voxel image is different form the one in my entire fixel image when I colour code the latter by afd.mif. Aren’t those supposed to be the same since max extracts the largest afd for every voxel?

If you show the relevant “data file” (e.g. afd.mif or other) that is stored in your fixel folder (where there’s also index.mif and directions.mif) in the fixel plot tool in mrview, the values you see there (e.g. via the colour scale) should match for the maximum fixel value in each voxel versus the value you extract via max. Make sure you’re each time talking about / looking at the same fixel “data file” though; i.e. the file you provide to fixel2voxel ..... max ..... should be the same as the one you’re looking at in the fixel plot tool, in terms of “things that match”.

@lottevdvoort

The output from both commands (the one you suggested with max operation and the one I sugegsted with split_data ) look similar but differ for some voxels.

That’s strange; unless we’re somewhere miscommunicating… Is the “some voxels” just a few voxels, or a lot of voxels…? Might be worthwhile to show us both images.

@lottevdvoort

Both however show the same maximum afd value which differs form my “total” afd value as mentioned above.

Ok, I’m a bit confused. So they do show the same maximum AFD value? (above I seemed to have understood they didn’t, for a few voxels?). Here you’re mentioning “total” AFD though: do you mean the sum/integral of all AFD in a voxel? That’s something else of course; that’s essentially summing the AFD of all fixels, rather than picking just the maximum one.

So… I’m a bit confused about what you’re seeing that matches, and what you’re seeing that differs between images…

[lottevdvoort]

@ThijsDhollander

Having said that, this depends also on your data quality and how you went about modelling it of course. I think I didn’t ask that bit yet above though: what is/are your b-value(s), and how many gradient directions do you have for each b-value (and how many b=0 images do you have)?

I have a NODDI data set with 97 directions: 8 b0, 27 b1000 and 62 b2500 and I used a reversed PE b0 for distortion correction.

@ThijsDhollander

If you show the relevant “data file” (e.g. afd.mif or other) that is stored in your fixel folder (where there’s also index.mif and directions.mif) in the fixel plot tool in mrview, the values you see there (e.g. via the colour scale) should match for the maximum fixel value in each voxel versus the value you extract via max .

Yes, that’s what I was expecting but somehow the maximum value in the extracted max is lower. I have used a threshold of 0.25 when segmenting fixels. With this threshold, the amount of fixels at the ends of the WM tracts, near the cortical gray matter, were mostly reduced to 1 or 2 or sometimes even 0 without sudden termination of visible ‘fixel’ tracts in the deep WM, like you said. Below are the image, the first one is the afd map shown calculated by your suggested max method and the second onen by using split_data and -coord 3 0. Both images are overlaid with the same fixels as used for calculating the afd maps.

I first compared both methods without thresholding the fixels and then only a few voxels differed; they showed different afd values. When using split_data somehow the lowest afd is 0 although this should be 0.25. I am certain that I’m using the same fixels/afd every time. Does mrview somehow scale the images?

@ThijsDhollander

Here you’re mentioning “total” AFD though: do you mean the sum/integral of all AFD in a voxel? That’s something else of course; that’s essentially summing the AFD of all fixels, rather than picking just the maximum one.

Ah now that makes sense! I thought that the maximum afd (as shown by the colour bar) belongs to one certain fixel in a voxel but if I understand correctly it’s the sum of all fixels in a voxel. Yep that would be more logical of course . I guess you can then ignore my earlier comments on this. Although I still do not understand why the split_data method shows an afd of 0 as the lowest value.

[ThijsDhollander]

@lottevdvoort

I have a NODDI data set with 97 directions: 8 b0, 27 b1000 and 62 b2500 and I used a reversed PE b0 for distortion correction.

Ok, that should be alright. We’ve recently started to discover and got to the bottom of MSMT-CSD even over-cleaning the WM FOD in some scenarios, including parts of cortical GM. But the bias we’re coming across should (in a weird sense) work “in your favour” here.

So now it’s “just” about picking a good threshold. Putting your earlier statement into this context:

@lottevdvoort

I do indeed see a sudden increase of the number of fixels at the very end of the corticospinal tract and the interface with gray matter. However, I am using msmt csd so I run fod2fixel only on the wm fods. Do you mean that the fixels you’re talking about in gray matter would still be visible for the wm fods as false positives?

So if you’re seeing that indeed, but with such data and MSMT-CSD, I think you can probably still safely increase the threshold a bit. It’s always a balance, so when false positives disappear, some true positives might go with it still; but I’m guessing that the balance is still in your favour here: there’s an opportunity to get rid of (potentially) loads of false positives near/in the cortex. For your application, as I’ve mentioned somewhere above, I reckon you’d also favour a “clean” fixel segmentation more than an inclusive but potentially messy one.

@lottevdvoort

Yes, that’s what I was expecting but somehow the maximum value in the extracted max is lower. I have used a threshold of 0.25 when segmenting fixels. With this threshold, the amount of fixels at the ends of the WM tracts, near the cortical gray matter, were mostly reduced to 1 or 2 or sometimes even 0 without sudden termination of visible ‘fixel’ tracts in the deep WM, like you said. Below are the image, the first one is the afd map shown calculated by your suggested max method and the second onen by using split_data and -coord 3 0 . Both images are overlaid with the same fixels as used for calculating the afd maps.

Ok, not sure this is the correct explanation, but based on your screenshots, I reckon it might have to do with how you “look” at the images. Note that each source (both the fixels as well as the first and the second voxel-wise intensity image) are shown with a different colour scale, i.e. windowing of that scale. A few things of note to get a good grip on this:

• When fixel images are opened just like you did (load the index and then colour by the data file; that’s correct!), it sets the windowing (i.e. the min and max limits of the colour bar) to the actual minimum and maximum of the values in the entire fixel image. You can change those manually in the tool’s relevant fields (bottom right in your screen shot), but that’s how they are set by default.

• For the main image you load however (i.e. the voxel wise intensity image), it sets the min-max windowing by default only based on the first slice that’s shown when you load it.

• The 2 described approaches (either using max operator, or your approach with split_data and the -coord selection) might differ in what they do for voxels that have no fixels in them. Sensible outcomes here would be either ending up with a zero in those voxels, or alternatively with a NaN (not a number) value. However, this will then influence the default windowing, as that doesn’t take into account NaNs, but it does take into account zeros. You can see this happening in your screenshots: the first one windows from a number larger than zero, whereas the second one windows from exactly zero. Also, when interpolation is on (which it is by default, hence the “smooth” quality of the images you see), you can kind of “see” NaN values, as they’re not interpolated. Those are the crisp looking “sharp” edges between the black background and the non-black intensities in the first screenshot.

So putting all that together, I think your perception of the difference between either image might be due to how they’re windowed? You can manually set the windowing of the main loaded image by opening the “view options” tool.

@lottevdvoort

I first compared both methods without thresholding the fixels and then only a few voxels differed; they showed different afd values. When using split_data somehow the lowest afd is 0 although this should be 0.25. I am certain that I’m using the same fixels/afd every time. Does mrview somehow scale the images?

…and this further aligns with that. That “lowest AFD” of zero is all those voxels without fixels (after segmentation). Don’t trust the viewer’s default windowing for this purpose. So manually set the windowing of either of those images to match the other: the images should look quite similar then. Apart from the NaN’s not being interpolated though. I reckon that aspect is the “few different voxels” you were spotting locally. I always prefer to switch interpolation off in any case, because it’s often misleading in several ways about the data you’re looking at. You can switch interpolation of by pressing the i button on your keyboard, or alternatively via the “colour map” menu.

@lottevdvoort

Ah now that makes sense! I thought that the maximum afd (as shown by the colour bar) belongs to one certain fixel in a voxel but if I understand correctly it’s the sum of all fixels in a voxel. Yep that would be more logical of course . I guess you can then ignore my earlier comments on this. Although I still do not understand why the split_data method shows an afd of 0 as the lowest value.

No, we’re getting ourselves into another confusion here. You should instead ignore my comment here: it’s now clear from your explanation that you definitely were not looking at the sum or “total” in any way. That was just a guess, based on your wording of “total …”; but it’s something you would’ve had to do on purpose, not something the viewer would bring up just like that. So to be clear: this is not it. I’m now pretty sure it’s the windowing in the viewer; the key message being: don’t use it as a guide for what is the minimum and maximum in your images. For most loaded “main” images in the viewer, its even only based on the first loaded slice, and not the entire volume.

So I think we’ve got it sorted out now; if all the above makes sense to you of course!

[lottevdvoort]

Thank you for the explanation about the viewer, makes things a lot more clear!

@ThijsDhollander

For the main image you load however (i.e. the voxel wise intensity image), it sets the min-max windowing by default only based on the first slice that’s shown when you load it.

Yep, I see now that there are voxels in the afd map that have values that go beyond the maximum intensity shown in the bar.

@ThijsDhollander

I think you can probably still safely increase the threshold a bit.

You mean going even higher than 0.25? A threshold of 0.20 already took care of the false positives near the cortical GM but I increased it to 0.25 just to be sure and because I indeed rather have some false negatives than false positives.

On last question; there are two volumes in the 4th dimension of the max output. The second volume is filled with zero’s but what exactly is this?

[ThijsDhollander]

@lottevdvoort

Yep, I see now that there are voxels in the afd map that have values that go beyond the maximum intensity shown in the bar.

Alright, good to hear that clarified things!

@lottevdvoort

You mean going even higher than 0.25? A threshold of 0.20 already took care of the false positives near the cortical GM but I increased it to 0.25 just to be sure and because I indeed rather have some false negatives than false positives.

No, 0.25 looks good! I hadn’t yet read that part when I was replying to the earlier bits. Your images definitely look sufficiently conservative: no false positives in sight near the cortex at all.

@lottevdvoort

On last question; there are two volumes in the 4th dimension of the max output. The second volume is filled with zero’s but what exactly is this?

That’s odd. So you only performed fixel2voxel with the max operator, and no other options? In that case, it sounds like a bug in MRtrix. I’d simply ignore the zero-filled volume. I’ll alert the developers, so they can look into it.

[ThijsDhollander]

Here you go: an issue for the developers to look into. Wait and see!

But in any case, it sounds like you can safely ignore the zero-filled volume; because of all the testing you’ve gone through, it looks like the first volume matches the alternative strategy, as well as the fixels’ data themselves. So I reckon you’re safe!

Update: looks like it’s just waiting for a release. But it’s confirmed to be a bug at least.

[lottevdvoort]

@ThijsDhollander

Update: looks like it’s just waiting for a release. But it’s confirmed to be a bug at least.

Alright!

Thanks again for all your help!

[lottevdvoort]

Hi @ThijsDhollander,

I still had a question about how the ‘size’ of fixels are determined by split_data? I was expecting this would be based on the input fixel file you provide to fixel2voxel but apparently it’s not. Meaning; if I provide an afd file to fixel2voxel and extract the afdmaps of the three ‘largest’ fixels via …

fixel2voxel wmfixels/afd.mif split_data - | mrconvert - -coord 3 0 peak1_afd.mif
fixel2voxel wmfixels/afd.mif split_data - | mrconvert - -coord 3 1 peak2_afd.mif
fixel2voxel wmfixels/afd.mif split_data - | mrconvert - -coord 3 2 peak3_afd.mif

… then for the same voxel in all those afdmaps (considering a voxel that has 3 fixels) the afd is not necessarily higher in the peak2 map than the peak3 map. What exactly is it that split_data stores in the 4th dimension? I thought these were the fixels in every voxels based on the afd value…

[ThijsDhollander]

@lottevdvoort

the afd is not necessarily higher in the peak2 map than the peak3 map. What exactly is it that split_data stores in the 4th dimension? I thought these were the fixels in every voxels based on the afd value…

Ok, just double-checked myself to be certain: you’re right, they’re not at all times guaranteed to be sorted from large to small (I wasn’t aware of that myself when I typed up any of the above!). Other than that though, the extracted values do match the AFD values from the original fixels. But basically, to get what you’re after in the example you provide, you can’t rely on the ordering in which they are stored.

Looking further into the options, I was then (now) thinking you could get this nonetheless by adding the -weighted option to fixel2voxel, and providing that option with the same afd.mif file. The idea here would be that the value you extract (the main argument to the command) is separate from the data that defines the “size” and that the latter would be used for sorting. However, trying that myself, I encountered a warning:

fixel2voxel: [WARNING] Option -weighted has no meaningful interpretation for the operation specified; ignoring

This is strange as there’s also the -number option to limit the extraction to a given (maximum) number of fixels. The documentation of this option reads:

  -number N
     use only the largest N fixels in calculation of the voxel-wise statistic;
     in the case of "split_data" and "split_dir", output only the largest N
     fixels, padding where necessary.

Suggesting at least a notion of sorting (otherwise there’d be no notion of the “largest N fixels”). But while that option does limit the number of extracted fixels, it still doesn’t sort them; probably also implying that its promise of extracting the largest N fixels might be broken (though I haven’t tested that). Finally, I tried to combine -number and -weighted, so as to explicitly introduce the notion of size again, hoping that -number would use it. But alas, the same warning popped up, and the output didn’t sort; so I’m again also left wondering whether the largest N promise holds here…

This is hard to check because most of the time the fixels do come out “almost” sorted. I suppose they just come out in the order they’re stored in the fixel data. That fixel data itself in turn probably gets its ordering just from how the fixels are segmented, but because the algorithm has a few rules about merging or separating individual peaks into single fixels based on relative amplitudes of peaks and valleys, I can imagine the final fixels don’t come out per se in the order the original separate peaks were encountered (e.g. when peaks are merged, I suppose). Lots of supposing and guessing here on my part in this little paragraph… so don’t take it for granted (but it does seem sensible to me in this way). @rsmith can probably help you out with the exact behaviour or hidden/unexpected quirks of fod2fixel followed by fixel2voxel.

However, that aside, the thing I did encounter with the -number not sorting makes me strongly wonder/question whether the “largest N fixels” promise there actually holds in practice, especially when the -weighted is sticking to the warning of not being meaningful for split_data… Even if this does still check out; the documentation could certainly be improved a bit, e.g. to at least mention when sorting can (not) be expected or relied upon. I’ll create another issue to sort out the documentation, or bug if there is one.

Cheers,
Thijs

[lottevdvoort]

Thanks for your quick response!

I tried the same thing with the -weighted option but surprisingly got the same warning

One thing I know for sure is that the way split_data and the -number option separate different fixels based on their ‘size’ is the same so there must be some logic behind it. It is definitely partially based on the AFD cause in most cases fixel1 shows a higher AFD than fixel2 and 3, just not in all…

[ThijsDhollander]

@lottevdvoort

I tried the same thing with the -weighted option but surprisingly got the same warning

A beautiful example of convergent thinking.

@lottevdvoort

It is definitely partially based on the AFD cause in most cases fixel1 shows a higher AFD than fixel2 and 3, just not in all…

Yes, this is what I kind of described (based on some guessing on what feels sensible on my part) in one of those paragraphs above: if my guessing there is right, then this is “just” the order in which they were already stored before fixel2voxel. That order itself, I’m guessing, is naturally quite likely to be not too different from the AFD values, due to how the fod2fixel algorithm operates. But also, due to how it operates, it’s not guaranteed to result in the kind of sorting you’re after (“it just often comes close”). This is, somewhat ironically, the very reason why this should probably be better documented: you have to look pretty closely to release it’s not exactly sorted. If it was 100% random with little to no pattern to it, you’d probably be less likely to assume it might be sorted. Either way, it’s obviously something that can benefit from more explicit documentation / warning.

I’ve created an issue for the developers over here: fixel2voxel split_data: -weighted and -number confusion · Issue #1790 · MRtrix3/mrtrix3 · GitHub . Thanks once more for your input / highlighting this issue here!

[kate3]

@jdtournier

spli

Is there an alternative option to split_dir with the new fixel2voxel command in matrix?

[kate3]

I try and use that option however I am given the error “unexpected value supplied for operation, valid choices are mean, sum, product, min, max, absmax, magmax, count, complexity, sf, dec_unit, dec_scaled, none.”

[jdtournier]

I have a feeling this functionality has been moved over into a new command fixel2peaks – rationale for the change was discussed on this GitHub issue, for those who might be interested.

[kate3]

The angle map I get is not ideal, and so I am trying to improve on this. I think what may be the problem in using fixel2peaks, is that the image is scaled based on the afd, and so I need to find a way to ‘descale’ this image, as I propose the scaling has affected the angle obtained?? Im not sure if that would be why or if the issue is something else

[rsmith]

I think what may be the problem in using fixel2peaks, is that the image is scaled based on the afd, and so I need to find a way to ‘descale’ this image, as I propose the scaling has affected the angle obtained?

This can be done manually:

mrmath peak.mif norm - | mrcalc peak.mif - -div peak_dir.mif

This assumes a volume with 3 volumes and just a single XYZ orientation per voxel; if encoding more than this you’d need to use mrconvert -coord 3 to extract the volumes corresponding to each peak, normalise each individually, then stitch them back together with mrcat.

[kate3]

Ok so I have a 122x122x68x3 image, so essentially 3 volumes so I won’t need to carry out mrconvert? Is that correct

[kate3]

Is the split_data functionality of fixel2voxel also possible? Reading this: https://github.com/MRtrix3/mrtrix3/pull/1874
it says “so it makes sense to replace split_data with none and ditch split_dir in favour of a dedicated command…” so I’m assuming the split_data functionality is still possible with the new fixel2voxel command, but by specifying ‘none’ as the operation? Thank you

[rsmith]

Ok so I have a 122x122x68x3 image, so essentially 3 volumes so I won’t need to carry out mrconvert? Is that correct

Correct: if there’s 3 volumes in the 4D image, then there’s only one fibre direction per voxel, and so you don’t need to do any trickery in order to separately normalise multiple fibre directions in each voxel.

I’m assuming the split_data functionality is still possible with the new fixel2voxel command, but by specifying ‘none’ as the operation?

Correct: that is a simple renaming of the operation, the functionality remains the same.