ENH fod2fixel: Output segmented FODs

[fionaEyoung]

When I want to run an analysis on individual fibre populations of an FOD image represented as SH distributions, I currently go through a very longwinded sequence: fod2fixel -> fixel2peaks -> separate peaks (e.g. mrconvert afd_scaled_peaks.mif -coord 3 0:2 f1.mif) -> peaks2fixel -> fixel2sh Additionally, the resulting "fixel ODFs" are essentially delta functions (scaled by AFD or peak amp or whatever) at the fixel peak directions, without any of the detail (e.g. dispersion) present in the original FOD lobe (like in the white boxed voxels in the images below).

Left: Original WM FOD, Right: First fixel converted to SH

It would be super useful to be able to directly acces the results of the FMLS FOD segmenter, before only the peak directions, AFD etc are retained and written out to fixel.

I've taken a look at the code to see if I have any hope of implementing this myself. Where I get lost is in fmls.cpp - I can't tell if the amplitude values of each direction sample are retained anywhere, which could be used to (approximately) reconstruct the individual lobes with something like amp2sh.

I appreciate the wishlist is very long 😉 I'd be game to try this myself but I'll probably need some *pointers, as my C++ knowledge is not advanced enough to handle the amount of abstraction in MRtrix's code.

[Lestropie]

Hi Fiona,

I would need to know exactly what you're trying to do with the output of that sequence of steps to know whether that particular sequence makes sense or whether there is some alternative. But my intentions for development in this context separately to your proposal may help to carve a path forward.

The key information that is yielded by the FMLS FOD segmentation process, but that is not currently exported to the file system, is the lookup table. This maps each of the 1,281 directions used for segmentation to a fixel index; and similarly, for each fixel there is a binary mask specifying which of those 1,281 directions was used in the construction of that fixel and the calculation of its properties. If those data were written to the fixel directory for subsequent use, then commands like tcksift and tcksift2 would correspondingly be modified to receive those data as input rather than the FOD image: the only reason they currently do the latter is so that they can access the lookup table information as it is useful for better informing streamline-fixel assignment. This has been on my internal todo list for many years, and would logically go with a bunch of other changes I have been chipping away at over a long period of time, but I have not commenced writing any corresponding code.

I suspect that this is a more direct expression of the limitation that you're encountering, and so directly addressing that limitation may better facilitate whatever it is you're intending to do downstream of this. If you're interested, then having fod2fixel optionally export those data would probably not be too difficult to implement, and you could then do with that information as you wish; modifying things like tcksift and tcksift2 to use that information is harder, but wouldn't be necessary for your own use case.

Rob

[fionaEyoung]

Hi Rob

Thanks for clarifying what's going on within the FMLS segmentation, I understand it much better now.* In that case, I believe that either the LUT or masks would indeed provide the information I'm looking for. I did some very crude experimentation with a fake directions mask (fixel_1_mask.mif) and achieved more or less what I'm aiming for by doing

sh2amp fod_sh.mif dirs.txt fod_amp.mif
mrcalc fixel_1_mask.mif fod_amp.mif 0 -if fixel_1_amp.mif
amp2sh fixel_1_amp.mif fixel_1_sh.mif

The trouble is with the high dimensionality (1000+ volume 4D image per fixel) things get a bit slow... Maybe the lookup table would be better on that front.

As for what I'm trying to do with the output, well that's partly a mystery to me too, such are the joys of research! But broadly I'd like to use it just like an FOD image, e.g. for tractography, running calculations on the SH coefficients, etc. It's all experimentation at the moment! I like the concept of fixels, but I worry about reducing each fibre population to a single peak direction plus scalar metrics, when the original FOD can capture a bit more nuance.

---

(*except this bit:

commands like tcksift and tcksift2 would correspondingly be modified to receive those data as input rather than the FOD image: the only reason they currently do the latter is so that they can access the lookup table information

How does the FOD image convey the lookup table information?)

[Lestropie]

amp2sh fixel_1_amp.mif fixel_1_sh.mif

The FOD exported here won't be of precisely the same shape as the lobe that you are attempting to mask, due to SH truncation. Have you considered the prospects of operating directly on the dixel representation?

Maybe the lookup table would be better on that front.

In the simple case, you just export an X x Y x Z x 1281 image, of 8-bit unsigned integer type, where each value indexes the relevant fixel offset within that voxel, and the 1281 cartesian directions get embedded in the image header. This is IMO quite achievable if you know a little C++ and don't mind wrapping your head around the basics of the MRtrix3 API. More than happy to provide feedback if you're willing to give it a try.

I like the concept of fixels, but I worry about reducing each fibre population to a single peak direction plus scalar metrics, when the original FOD can capture a bit more nuance.

Agreed; you might find this abstract interesting. Have not touched that code for many years though, and its optimisation convergence is far from guaranteed unfortunately.

How does the FOD image convey the lookup table information?

It's not "conveyed" in the FOD image. Those commands receive as input the FOD image and themselves run the very same FOD segmentation that is invoked by fod2fixel. They do that in order to retain and utilise the lookup table information that is not included in the output of fod2fixel. But it means that if a pipeline is invoking both fod2fixel and tcksift(2), FOD FMLS segmentation is actually being executed twice. Would be preferable to instead run fod2fixel, have its output include the lookup table information, and then feed that as the input to tcksift(2) instead of the FOD image.

[fionaEyoung]

This is IMO quite achievable

Finally gave this a whirl, turns out it was indeed achievable. Using the changes in #2621 I can achieve the behaviour I described above:

fod2fixel fod_sh.mif fixels -lut lut.mif
sh2amp fod_sh.mif fixels/lut.mif fod_amp.mif # modified sh2amp, see #2621 
f=0 # 1st fixel
mrcalc fixels/lut.mif $f -eq A1_FOD_wm_dixel.mif 0 -if fixel_${f}_amp.mif
amp2sh fixel_${f}_amp.mif fixel_${f}_sh.mif

fod_sh.mif:       fixel_0_amp.mif:     fixel_0_sh.mif: