Thresholding fixels using voxelwise threshold

1

Dear MRtrixers,

I am trying to replicate some tractography results based on deterministic tractography using FOD data (I know not many here are fans of using a deterministic approach in this case). I am running tckgen with the FACT algorithm, giving as input an image created with fod2fixel, combined with fixel2voxel (and the split_dir option).

However, before I’m interested in cropping the fixels whose amplitude is smaller (<5%) compared with the maximal fixel amplitude of the same voxel. Essentially, I’m want to follow this work by Dell’Acqua.

However, I’m not sure how this can be implemented. Is there a simple way to apply a voxel-wise threshold for all fixels?

Many thanks,
Best,

Roey

2

Indeed, not what we’d typically recommend…

However, given your data are in fixel format, I think you might be able to do what you’re after with a combination of:

  • fixel2voxel using the max operation to generate a voxel-wise image of the maximum fixel intensity
  • voxel2fixel to map that back onto fixels
  • and finally, a call to mrcalc like this:
    mrcalc fixel_amplitudes.mif fixel_max.mif 0.05 -mult -gt fixel_amplitudes.mif 0 -if thresholded_fixels.mif
    to set amplitudes lower than 5% of the max to zero.

This is untested, by the way, but hopefully enough to give you a feel for how you might achieve this.

3

Thank you so much, Donald!
Indeed, I was missing the voxel2fixel step after calculating the maximum fixel intensity. I will test according to your suggestion.
Many thanks.

4

Dear all,

First of all, here’s an example of the code I’m using now to crop the fixels I give as input to the FACT tractography. Names of files are written in between <>. I now create a binary mask of fixels to keep (based on some voxel-wise relative amplitude threshold). I later apply this mask to crop the fixels (step 6).

% (1) Segment FOD to fixels (create fixelDir with a fixels-peak-amplitudes file) fod2fixel -peak

% (2) Split the fixels to 3 volumes per fixel (fixelsNum x 3 file), later used for FACT tractography (create fixelPeakDirsFile file)
fixelFile = fullfile(fixelDir,‘directions.mif’);
fixel2voxel split_dir

% (3) Find the maximal amplitude of each voxel (create fixelPeaksMaxFile file)
fixelPeaksAmpsFile = fullfile(fixelDir,fixelPeaksAmpsFile);
fixel2voxel max

% (4) Map this maximal amplitude back to fixel space (fixelsNum x 1 file) (create fixelPeaksMaxFixelSpaceFile file)
voxel2fixel

% (5) Create a binary mask of the fixels, at 10% of the maximal FOD in each (create thresholdFixelsMaskFile file)
mrcalc <fullfile(fixelDir,fixelPeaksMaxFixelSpaceFile)> 0.1 -mult -gt 1 0 -if

% (6) Crop the fixels using the fixels mask (this creates a new fixelDir: croppedFixelDir)
fixelcrop

% (7) Extract the cropped fixels to a file with 3 volumes per fixel (fixelsNum x 3 file) (create croppedFixelPeakDirsFile file)
croppedFixelFile = fullfile(croppedFixelDir,‘directions.mif’);
croppedFixelPeakDirsFile = fullfile(croppedFixelDir,‘fixel_dirs3vols.mif’);
fixel2voxel split_dir

% (8) Run FACT tractography (create tckFile)
tckgen -algorithm FACT -seed_image -mask -select 500000 croppedFixelPeakDirsFile>

Similarly, I also create an absolute threshold mask of 0.2 amplitude.

My problem now is that the tractography results seem to be contaminated by false positives that come from spurious fixels peaks. For example, in the corpus callosum many streamlines have anterior-posterior orientation (showing here dummy runs with a total of 50k streamlines):

1RelativeThresh0000
1RelativeThresh00001535×1130 442 KB

If I increase the cutoff in the tractography to 0.2 (instead of the default 0.05), I get similar results:

20000
200001526×1130 456 KB

However, if I first crop the fixels using my absolute binary mask of 0.2 amplitude, the results are much cleaner:

20000
200001526×1130 432 KB

Is this possibly related to the following post?

My question is: Do you have an idea why the cutoff option doesn’t do what I expect it to do?

In addition, judging by the cube-like structure shown in the tractography image above, I understand that perhaps I should have upsampled the diffusion FOD data before continuing to the next steps.

Finally, another short question. In many works (e.g., this one) I’ve seen that people also use an absolute threshold on the peak FOD, which is defined as 3 “multiples of the amplitude of a spherical FOD obtained from an isotropic voxel of gray matter”. Could you please tell me how you interpret this? My first thought was taking the l=0 volume of the FOD file, and average across gray-matter voxels (using the 5tt file, for example). However, this doesn’t seem right, and perhaps I need to convert this to an FOD amplitude somehow.

Sorry for the lengthy descriptions.

Many thanks!
Roey

5

I now create a binary mask of fixels to keep …

(emphasis mine)

It’s worth noting that the behaviour of the tracking may vary slightly depending on whether you remove fixels that are below some threshold, versus modulating the amplitudes of fixels to ensure that they fall below the global amplitude threshold used within the FACT algorithm in tckgen. Consider a streamline whose orientation is in between two fixels, slightly closer to the smaller one, but the smaller fixel is below the voxel-wise volume fraction threshold. In the former case, the streamline will be attributed to the larger fixel, since it’s the fixel whose orientation is closest to the incoming tangent, and tracking will continue. In the latter, the streamline will be attributed to the smaller fixel, but its amplitude will be detected as sub-threshold, and hence the streamline will be terminated there.

Don’t know how particular you want to be on such details, so just thought I’d raise.

Do you have an idea why the cutoff option doesn’t do what I expect it to do?

It’s hard to know from screenshots alone; a cursory look at the code suggests that the fixel amplitude threshold should be doing as intended. We’ll need to try to duplicate the steps you’ve done here to see if there’s an issue.

I’ve seen that people also use an absolute threshold on the peak FOD …

Firstly, it’s important to be specific about what measure is being used, whether a peak FOD amplitude or a fixel integral. The former is obtained with the -peak option in fod2fixel; the latter is what we prefer to use for FBA, and is what fod2fixel will give you with the -afd option.

This is important because if you want to apply a threshold based on your description here:

3 multiples of the amplitude of a spherical FOD obtained from an isotropic voxel of gray matter

, if you were to derive this amplitude, it would only make sense to apply such a threshold to FOD peak amplitudes, not fixel integrals, and there’s a chance that that may not be consistent with how you’ve been applying fixel amplitude thresholds previously.

If you have a spherical harmonic expansion where the l=0 term is some fixed value (e.g. as averaged from GM voxels) and all other SH terms are zero (i.e. isotropic FOD), then converting that l=0 coefficient to an amplitude requires division by sqrt(4pi).

Rob

Edit: Whoops, that’ll teach me for addressing posts from bottom to top… Already addressed elsewhere.

6

Dear Rob,
Thank you for these insights. It’s (quite literally) not an angle I had thought of.
I will have to think about the FOD integral vs peak issue.
Thanks again,
Roey