Simultaneously segmenting FODs to fixels while warping to template

[draffelt]

Computing fixel apparent fibre density (FD) and fixel correspondence in the current fixel-based analysis (FBA) pipeline is performed by:

1. warp FOD to template with FOD reorientation and optional modulation - using mrtransform
2. FOD segmentation to get AFD per fixel - using fod2fixel
3. for each fixel in template fixel mask, find the corresponding subject fixel based on direction - this happens as a first step within fixelcfestats

I'm planning to write a single command that combines the first two steps. Basically it will accept a subject FOD image, a warp, and a template image. For each voxel in template space, the warp will be used to find the corresponding FOD in subject space. The (interpolated) FOD will be then segmented into fixels (in subject space), and the fixel data (orientations and FD) saved in template space (with reorentation of the fixel direction). Step 3 above will be performed using the dedicated fixelcorrespondence command, and this step will be removed from fixelcfestats.

I can see several benefits to this approach:

1. FODs will be segmented in native space. This may be beneficial if severe warps are required to map the subject to template (i.e. extreme Jacobians), which can distort the FODs during reorientation. While using many apodised PSFs mitigates this, it's still possible to get changes to the FOD lobe for really extreme Jacobians (e.g. single lobe splitting into two).
2. Once fixelcorrespondence is performed, modulation can be applied in template space with a simple fixel-wise multiplication with the 'fibre cross-section (FC)' measure (from our upcoming fixel-based morphology manuscript). Hopefully this will make clearer what is happening from the users perspective (i.e. the combination of fibre density and cross-section to get third measure called "fibre density and cross-section, FDC" (what we used to call modulated AFD)). It will simplify the number of steps in a FBA, and speed it up since warping FODs (with and without modulation) with large numbers of aPSFs @ lmax = 8 is slow, and fixel segmentation happens once only (pre modulation).
3. It will be possible in the future to substitute this single fod segmentation/warping command for another to perform FBA of measures derived from other diffusion MRI models. For example, instead of warping and interpolating FODs then segmenting to get fixels, one could warp and interpolate DWIs then compute fixels measures via CHARMED etc. This idea of keeping FBA in MRtrix generic fits in nicely with the FBA flow diagram in my current manuscript.
4. If fixel angular correspondence is a separate step performed before (and not within) fixelcfestats, then it means we could input fixel data for all subjects as a single folder using the proposed fixel fomat. Basically the index and directions file would be the same for all subjects, with only the fixel data being different. We could just label each subject's data as 1,2,3 etc within the fixel folder to assign each subject to rows of the data/design matrices. A whole bunch of the current loading and indexing code in fixelcfestats would be simplified. I could also remove the template fixel mask input from fixelcfestats too. Usage would then be a trivial fixelcfestats input_fixel_folder design contrast tractogram output_stats_fixel_folder

Aside from being a brain dump, I'm posting as an issue mainly to get your thoughts on if you think this warrants a separate command, or if it could be buried in fod2fixel as a -warp option. In terms of the fod2fixel code, @Lestropie it would be a completely different branch from the current code but that doesn't matter from a user's perspective.

[Lestropie]

Adding some comments I made in discussion into here.

1. For cases such as having a single lobe split into two, I'd still like something like fixelcorrespondence to explicitly detect such cases (e.g. two subject fixels mapping to a single template fixel) and handle appropriately (e.g. sum FD), rather than only permitting a one-to-one fixel mapping between subjects and template, which could potentially be introducing unwanted variance. I'd also like to see for myself the effects that are occurring to cause this sort of thing, just to confirm in my own mind that we're solving the right problem the right way - the lobe splitting may be a symptom of point 2 below.
2. This is making the implicit assumption that all fibres within the fixel are coherently oriented. If there's any orientation dispersion, then technically different directions within the lobe will experience slightly different modulations, which will result in modification of the shape of the FOD. Whether one approach is 'better' than the other given the limitations of each requires experience in working with the data that I simply don't have; but this would need to be acknowledged and discussed nonetheless.
3. Re-gridding images (DWIs or FODs) to template space then subsequently fitting a model will still probably be more familiar to readers / users, and should still be possible to do. Having a single command to warp DWIs and fit a model is hypothetically possible, but we shouldn't sell it as a prerequisite. This should probably also be likened to QSDR, which fits the diffusion model in template space taking the subject-to-template warp field into account.

Last point: I'd probably vote for a new command. If this were a fod2fixel option, you'd then have to explicitly match fixels between the FD and FC images (wherever the latter were to come from) in order to enable explicit modulation. As its own command, you could output data for both of these two, ensuring that the target fixel image that the data refer to is the same, so a straight element-wise multiplication could be done to derive FDC.

[draffelt]

For cases such as having a single lobe split into two, I'd still like something like fixelcorrespondence to explicitly detect such cases

Agreed that fixelcorrespondence can be improved in the future. I have a few ideas myself. A benefit of it being independent from fixelcfestats is that a user can try/swap in another method if he/she desires. Keeping in mind that not all DWI fixel models have information on spread etc.

This is making the implicit assumption that all fibres within the fixel are coherently oriented. If there's any orientation dispersion, then technically different directions within the lobe will experience slightly different modulations, which will result in modification of the shape of the FOD.

Yes, fixel vs dixel modulation may be subtly different, however I'd be surprised if it had any effect on the end results. I suppose one may introduce more variance than the other, but I should be able to get a feel for this with our clinical studies. I do like the fact that with fixel modulation we have a FBA method (and tutorial) that can be applied to any fixel-based FD measure by only changing a single step. As also noted below, FC is computed currently from a template fixel mask and each subject's warp. I.e. it is already strictly fixel based... I don't estimate FC by first breaking up the template FOD into dixels, computing FC per dixel, then somehow combining all dixel FC for each lobe/fixel. So personally I think it's cleaner to have the entire FD, FC, FDC analysis being fixel-based. Modulating FODs via dixels to get FDC, then computing FC via fixels kinda breaks the whole FDC = FD * FC concept (even though I admit the result is likely to be the same).

Re-gridding images (DWIs or FODs) to template space then subsequently fitting a model will still probably be more familiar to readers / users, and should still be possible to do. Having a single command to warp DWIs and fit a model is hypothetically possible, but we shouldn't sell it as a prerequisite.

Just to be clear, this is not what I was proposing for analysing another measures (such as CHARMED). Warping DWI (qspace) data to template space still requires reorientation with some sort of fibre basis functions, and therefore may still introduce issues/spreading due to distortions. I was proposing that the warp defines where to interpolate the DWI, the model is fit in subject space on the interpolated DWI data, and then fixel orientations are reoriented before saving in template space.

Last point: I'd probably vote for a new command. If this were a fod2fixel option, you'd then have to explicitly match fixels between the FD and FC images (wherever the latter were to come from) in order to enable explicit modulation.

The FC images are computed only from the template fixel mask (and warps from each subject). So only FD measures would still require fixelcorrespondence to be in 'template space'. However I think performing fixelcorrespondence separately is the way to go, since it can be applied to any measure and can be substituted for different methods if the users wants.

As its own command, you could output data for both of these two, ensuring that the target fixel image that the data refer to is the same, so a straight element-wise multiplication could be done to derive FDC.

I still need to commit a new command to compute FC from the new warp format. Though I'm not sure combining it into a single command that also computes FD from FODs make sense. FC is DWI model independent, so I think it would be better to be a stand alone command. On that note my old FC command was called fixelmodulate. Can anyone think of a better name?

I'm not yet convinced that the proposed command should be part of fod2fixel, but if it's a stand alone command, I'm struggling to think of what name I would use. It would be something like fod2fixeltemplate.

[draffelt]

Ok, So I've come up with a simpler (and more unix-like) solution to this:

1. A FOD image is transformed to template space without reorientation.
2. The FODs are segmented with fod2fixel (i.e. with no FOD lobe distortions from reorientation).
3. Fixel directions are reoriented with the new fixelreorient (still in the registration branch in Gitlab).
4. Angular correspondence is then obtained with the current fixelcorrespondence.

I have tested the difference between this approach and the old method (warping FODs with reorientation...then segmenting to fixels). You might be glad to know there is very little difference between the two in both fixel number, orientation and AFD integral, suggesting apodised PSFs in FOD reorientation are doing their job. I could see where the odd FOD lobe was split into two (and vice versa), however in general the result is very similar. I should probably note that it was only a quick and dirty registration with jacobian determinants, min: 0.469422, max: 1.79063. I will test it out with hi-res data and more severe warps to be sure.

So the benefit of this proposed method is that any fixel-based measure can now be analysed in MRtrix by:

1. Warping the raw DWI data without reorientation (presumably with FOD computed warps)
2. computing the DWI model (e.g. CHARMED)
3. then using fixelreorient and fixelcorrespondence as per FOD analysis above.

So, the question now is (@jdtournier), how can we best modify piping so that we can also pipe fixel images and combine some of the above steps?

Also, are there any suggestions on the naming of fixelcorrespondence? Fixelangularcorrespondence is probably more appropriate but too long.

[thijsdhollander]

Small side note: those minor differences you might be seeing, might not even be due to the apodised PSF reorientation doing its job "slightly less well". Even when it's doing it's job perfectly, the difference in order of the reorientation and segmentation steps may in certain cases slightly change the result, due to the segmentation's criteria to separate (or not separate) lobes: it uses the amplitude of peaks relative to the amplitude of the lowest valley in between them to decide on this, but exactly this amplitude ratio can change due to reorientation (as in: even when it's done theoretically correct, without discretisation issues due to apodised PSFs for a finite number of orientations).

...but all in all, it's indeed reassuring that this only has minor impact on typical realistic results. :+1:

[jdtournier]

Sounds fine to me - I must admit I'm feeling a little out of the loop on this stuff at the moment... But just on the topic of piping fixel images:

The best thing to do is something similar to what is used currently for regular images. In contrast to tractography, which is typically processed in a stream fashion (i.e. no random access required), images typically need to be processed to completion by the first command before processing can start with the next. So streaming the data through a pipeline doesn't work so well - it can work, but it introduces quite a bit of latency as the data needs to be copied through into the address space of the next command. The approach used for images is to write them to temporary storage, pass the filename through to the next command once processing is complete, and the next command will access the image via memory mapping. Basically this ensures a zero-copy pipeline (although this not strictly true any more now that we have delayed write-back - maybe we can still use direct memory mapping for temporary images...? Might start another issue on the topic).

The difficulty with fixel images is that they may eventually become folders...? However, as things stand we should be able to pipe them through as-is, given that they ride on top of the regular Image framework. I haven't looked into it, but I wouldn't be surprised if this happened to just work... Have you tried it?

[Lestropie]

However, as things stand we should be able to pipe them through as-is, given that they ride on top of the regular Image framework. I haven't looked into it, but I wouldn't be surprised if this happened to just work... Have you tried it?

I'd expect it to not work out of the box. Piping and sparse images use two separate handlers. You'd need to implement a piped sparse image handler; then, because the standard pipe handler would load the file name from standard in but its sanity check would fail, you'd need to write that temp file name to a string variable somewhere so that further down the handler list the piped sparse image handler could see the "-", load that file name, and wrap the sparse image handling, deleting the file on close.

[jdtournier]

Ah, yes. I'd forgotten this was implemented using a dedicated handler... Yes, it would require a bit more work, but nothing impossible.

[draffelt]

Will close this for now since the main issue has been resolved. Will have a think about piping fixel images and reopen another issue if required.