Probabilistic tracking on peaks


I would like to run probabilistic tracking on peaks (same format as the ones created from sh2peaks). I can do deterministic tracking using FACT, but I would like probabilistic. The peaks of course do not contain an orientation distribution which a probabilistic algorithm could sample from. But I would be fine with setting that manually for the entire image: e.g. sample from a normal distribution around each peak with a fixed standard deviation.
I guess one way would be to transform the peak image to an image of spherical harmonics and then I could run iFOD2. However, I do not really know how to do this transformation. Is there already some code in MRtrix that is doing something similar?


Unfortunately not, sorry… It would be trivial to code up though, if you feel like getting your hands dirty…? :grin:

You could start with one of our demo bits of code, figure out how to get it working, then throw in the relevant call to generate and store the SH coefficients for a delta function oriented along the desired direction. Should be a piece of cake. :cake:

Happy to assist if you get stuck… :+1:


Hi Jakob,

You might want to take a look at fixel2sh. The input needs to be in the fixel format (which is essentially vectors), but otherwise it does exactly what you’re after.



:blush: OK, that’s undeniably simpler…


Great. Thanks for your help.
When creating fixels I get a index.mif and directions.mif file. The index.mif file has 2 values for each voxel. I guess the first one is telling how many fixels there are and the second one is telling the index in the directions.mif file. If there is more than 1 fixel it is telling the index of the first fixel. For the index of the other fixels I just have to add +1/+2/+3… to the index?
Is this the correct interpretation of the format? If not, is there a place where it is documented?


I got it running and it is working as expected. Thanks a lot!

In case somebody else will need this in the future this is the code to transform a peak image (in this case with only 1 peak per voxel) to a fixel image. The fixel image can then be transformed to spherical harmonics.

from os.path import join
import nibabel as nib
import numpy as np

args = sys.argv[1:]
# Only works for peaks with one peak per voxel
peaks_in = args[0]  # expected shape: [x,y,z,3]
fixel_dir_out = args[1]

peaks_img = nib.load(peaks_in)
peaks = peaks_img.get_data()
s = peaks.shape

directions = []
index = np.zeros(list(s[:3]) + [2])
amplitudes = []

idx_ctr = 0
for x in range(s[0]):
    for y in range(s[1]):
        for z in range(s[2]):
            peak = peaks[x,y,z]
            peak_len = np.linalg.norm(peak)
            if peak_len > 0:
                peak_normalized = peak / (peak_len + 1e-20)
                index[x,y,z] = [1, idx_ctr]
                idx_ctr += 1, np.eye(4)), join(fixel_dir_out, "directions.nii.gz")), peaks_img.get_affine()), join(fixel_dir_out, "index.nii.gz")), np.eye(4)), join(fixel_dir_out, "amplitudes.nii.gz"))

Then run

fixel2sh amplitudes.nii.gz sh.nii.gz