Just iFOD1 and iFOD2 works, FACT not

Hi MRTrix community
I have a problem with the tckgen command for alogirthm different from iFOD. Here is the code:
Dabei ist alles in Python geschrieben und hier sind die Parameter:

opt_1 = ‘-number 200k -act 5TT.mif’
algorithm = “FACT”

os.system(’‘5ttgen fsl /home/ubuntu/E130-Daten/IngoHermann/%r/T1w/T1w_acpc_dc_restore_brain.nii.gz 5TT.mif -premasked’ % input_number)
os.system(‘mrconvert data.nii.gz %s.mif -fslgrad bvecs bvals -datatype float32 -stride 0,0,0,1’ % name_DWI)

os.system(‘dwibiascorrect -fsl %s.mif -’ % (name_DWI))

os.system('dwi2mask %s.mif %s_mask.mif' % (name_DWI, name_DWI))
os.system('dwi2response tournier %s.mif %s_response.txt' % (name_DWI, name_DWI))
os.system('dwiextract %s.mif %s_ex.mif -bzero' % (name_DWI, name_DWI))
os.system('dwiextract %s_ex.mif - -bzero | mrmath - mean meanb0.mif -axis 3' % name_DWI)
os.system('dwi2fod csd %s.mif %s_response.txt %s_fod.mif -mask %s_mask.mif' % (name_DWI_2, name_DWI, 
    os.system('tckgen %s_fod.mif %s_tckgen_%s_%s.tck -seed_image %s_mask.mif -mask %s_mask.mif -algorithm %s %s' % (name_DWI_2, name_DWI_2, algorithm, name, name_DWI, name_DWI, algorithm, opt_1))

And I get this image:
FACT (wrong labeled before as SD_Stream):


SD_Stream (wrong labeled before as FACT):

Why does it look so?

Hi @Ingo,

The issue here is that the main input file that needs to be provided to tckgen can change depending on the particular tracking algorithm being invoked; this is listed in the command’s help page.

So when you provide file %s_fod.mif to tckgen -algorithm fact, which is an image containing spherical harmonics, most likely of maximal harmonic order 8 (45 volumes), the FACT algorithm is mis-interpreting these data as 15 triplets of XYZ fibre directions in each voxel. This will inevitably lead to erroneous tracking (I’m guessing the images have been mis-labeled: the upper erroneous image is FACT, whereas the lower one is SD_Stream). tckgen is unable to detect the fact that you’ve provided the incorrect image since it actually satisfies the requirements for that algorithm (number of input volumes must be a multiple of 3).

To get tckgen -algorithm fact working, you can use either sh2peaks, or fod2fixel combined with fixel2voxel split_dir.

Cheers
Rob

Ok thank you, the labeling was swapped :D. With

fod2fixel %s_fod.mif -afd %s_fixel.msf
fixel2voxel %s_fixel.msf split_dir %s_voxel.mif

it worked quite good alike the iFOD algorithm. But why does SD_Stream look different and I also testet Nulldist and there one gets this noise image:
Nulldist2:


The new FACT:

Ok I fixed the SD_Stream problem. One has to use as option of the tckgen

-angle -77

or another number in this area because if you compare different angles one get:
Without any angle option:


With -angle 90:

With -angle 77:

OK, thanks for reporting those results; it opens up an interesting discussion.

For the sd_stream algorithm, the default step size is 1/10th the voxel size, and the angle of curvature constraint is applied at every step. Since the steps are so small, the default angle of curvature is correspondingly small: just 9 degrees per step (equivalent to 90 degrees per voxel, which is the same as that used in the default iFOD2 algorithm). By increasing the permitted angle to ~ 90 degrees per step, essentially no curvature constraint is applied at all, with the tracking algorithm instead relying entirely on FOD amplitude / mask constraints.

However with the default curvature constraint, in my quick experiment, over half of the streamlines generated are terminated due to the curvature constraint; and as you observed, there are ‘strips’ of white matter where there are no streamlines. This I suspect is due to the tracking algorithm performing interpolation directly on the FODs in their native spherical harmonic basis: If the angle between FOD peaks in successive voxels is relatively small, then the peak orientation will ‘rotate’ relatively smoothly as this SH interpolation is performed between the two voxels; however if the angle is larger, then the peak orientation will not rotate smoothly; interpolation will instead result in one peak diminishing until it disappears, and the other growing until the peak appears, as the position moves from one voxel centre to the other. Consequently, as a streamline tracks between those two voxels, it will experience a sudden change in orientation in one particular step, rather than a smooth transition in orientation across all steps, and this is more likely to exceed the curvature threshold.

One thing that can help in this situation is to up-sample the DWI volumes prior to spherical deconvolution, and then manually set the step size and angle parameters to what they would have been for the native resolution data. Interpolating the DWI data tends to produce smoother transitions in FOD peak orientation than interpolating the FODs themselves.

@jdtournier Any thoughts on adjusting the default curvature threshold for this algorithm?

Yes, @rsmith is spot on. This particular issue is discussed in detail in this paper (see figure 2), and is one of the reasons I don’t recommend its use in practice. The simplest ‘fix’ is to open up the angle of curvature to allow sharp turns, as you’ve already figured out. But that does allow unnatural ‘kinks’ in the trajectories, and increases the potential for artefacts such as streamlines making U-turns, etc. So it’s a delicate balancing act between allowing streamlines to progress along what must be the most likely direction of travel, and introducing too many spurious streamlines.

This issue isn’t going to be specific to this implementation by the way, any multi-fibre deterministic algorithm will have to deal with situations where the orientation being followed is no longer present at the next step (as illustrated in the paper above, or simply because the SNR is no longer sufficient to recover that fibre population, etc). The algorithm will then have to decide whether to terminate, or follow the next closest orientation. The question is then: how close is plausible? This is what the -angle option allows you to specify - but it is inherently a compromise.

What is undeniable though is that the default angle is not appropriate for this algorithm. In fact, in MRtrix 0.2, the default is to disable the angle constraint altogether (-curvature defaults to 0). I guess I’ve neglected to set this appropriately since I essentially never use that algorithm… I think the better default would be a hard value (i.e. not dependent on step size or voxel size, as per the current default), maybe between the 45° and 70° mark - any wider would allow complete changes in direction, which we should assume are implausible. @Ingo: how did you come up with the 77° value by the way? It seems very specific, was there a dramatic difference at that point…?

Sry for my late response but I was in holidays. The 77° was just an arbitrary number :smiley:

Can you still help me with the Nulldist2 algorithm as seen above? Why is that image so noisy?

Ah sorry, I forgot about that one in my response:

The nulldist / nulldist2 algorithms are supposed to be noisy; that’s their defining attribute (“Null Distribution”). I would suggest reading the relevant article that is cited in the tckgen command’s reference section.

This thread and the answers from @rsmith and @jdtournier have been very helpful in coming up to speed with different algorithms! So I’m getting the idea that the SD_Stream algorithm is essentially deprecated in favor of ifod1 and ifod2?

FWIW

I was getting ready to post a long question as to why SD_Stream produced tracts limited in space:

Compared to the ifod1 algorithm (with essentially the same parameters):

The I added the -angle 77 to the SD_Stream and got a better result:

Well, deprecated may be too strong a word here. But I personally don’t recommend it… Irrespective of the oddities of deterministic tracking on FODs, I think there’s just far too much uncertainty in the whole process (far beyond just imaging noise) to justify relying on a deterministic approach. The default is currently iFOD2, and that would be my recommendation for tractography in general. Of course, there may be cases where other approaches perform better (not sure they’d necessarily be more appropriate though), and bear in mind that others may also have different opinions on the matter (I know many other researchers in tractography are far more pro-deterministic than I am…).

And then there is also a whole spectrum of “more (or less) probabilistic/deterministic” of course. Using the -power option, one can “make” e.g. iFOD2 also behave more (or less) deterministic, all the way up to letting it behave very close to deterministic…

Hi @jdtournier,
I have a question about above figure: how to generate the 3d images ?

Thanks!

Best Wishes,
Jian

I’m not sure what you mean here. There’s an option to hide the main image in the view menu, which you can use in combination with the tractography tool’s ‘crop to slab’ option. Another possibility is the MRView ‘volume render’ mode, also in the view menu. This works best with clip planes, which you can add via the ‘view options’ tool. Otherwise, you’ll have to be a bit more specific as to what you’re after…

Hi @jdtournier,
Sorry, I didn’t fully clear.
I’ve got to the answer from your answer. Thank you! I can use the tractography tool’s ‘crop to slab’ option to display my results.
However, can we rotate with different angle to watch the 3d images in MRView?

Thanks!

Best Wishes,
Jian

This is the button you’re looking for:

Hi @jdtournier ,
Thank you for your detailed answer!

Best Wishes,
Jian

Hello,

I would like to use the FACT algorithm as well and tried the above methodology but to no success.

Here were the following steps:

dwi2fod csd eddy.mif wm_response.txt fod.mif -mask eddy_mask.mif
fod2fixel fod.mif fixel.msf

I had to ‘cd’ into the fixel directory since the fixel.msf was not accepted as an input for fixel2voxel.

fixel2voxel directions.mif none voxel.mif

tckgen voxel.mif dti.tck -algorithm FACT -seed_image eddy_mask.mif -mask eddy_mask.mif -select 5000000

Output:

tckgen: [ERROR] Number of volumes in FACT algorithm input image should be a multiple of 3

Any help would be appreciated.

Thanks,

Vinny

Hello,

I got the FACT to work by converting fod2fixel and then fixel2peaks. This produced an image that was a multiple of 3, and was used as input for tckgen with the FACT algorithm.

I’m assuming this is CSD followed by FACT.

Regards,

Vinny