Probabilistic tracking on peaks

Hi,
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.

Cheers,
Thijs

: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)
                directions.append(peak_normalized)
                amplitudes.append(peak_len)
                index[x,y,z] = [1, idx_ctr]
                idx_ctr += 1

nib.save(nib.Nifti2Image(np.array(directions), np.eye(4)), join(fixel_dir_out, "directions.nii.gz"))
nib.save(nib.Nifti2Image(index, peaks_img.get_affine()), join(fixel_dir_out, "index.nii.gz"))
nib.save(nib.Nifti2Image(np.array(amplitudes), np.eye(4)), join(fixel_dir_out, "amplitudes.nii.gz"))

Then run

fixel2sh amplitudes.nii.gz sh.nii.gz

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.

The answers above relate to the second option (i.e. generating an SH image), but the former is perfectly feasible also. The tckgen code is designed to be relatively easy to insert additional tracking algorithms. Essentially all you would need is a copy of the FACT algorithm, which upon entering a new voxel, randomly perturbs the streamline orientation after selecting the nearest fixel, and then stores the sampled orientation for that particular fixel for that streamline (so that a new randomly perturbed orientation is not generated if the next streamline vertex is still in the same voxel).

Of course this would require C++ chops to pull off. It would also not be possible to perform sub-voxel interpolation with such an algorithm, whereas iFOD1/2 interpolate the FOD SH coefficient directly.