# Probabilistic tracking on peaks
**URL:** https://community.mrtrix.org/t/probabilistic-tracking-on-peaks/2030
**Category:** Uncategorized
**Created:** [November 7, 2018, 10:43am UTC](https://community.mrtrix.org/t/probabilistic-tracking-on-peaks/2030 "2018-11-07T10:43:27Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@wasserth](https://community.mrtrix.org/u/wasserth)
#### Post date: [November 7, 2018, 10:43am UTC](https://community.mrtrix.org/t/probabilistic-tracking-on-peaks/2030/1 "2018-11-07T10:43:27Z")
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?
---
### Author:  [@jdtournier](https://community.mrtrix.org/u/jdtournier)
#### Post date: [November 13, 2018, 11:52am UTC](https://community.mrtrix.org/t/probabilistic-tracking-on-peaks/2030/2 "2018-11-13T11:52:50Z")
> [@wasserth](#):
>
> 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…? 😁
You could start with one of our [demo bits of code](http://www.mrtrix.org/developer-documentation/example_per_voxel_multithreaded_4_d_processing.html), figure out how to [get it working](http://www.mrtrix.org/developer-documentation/build_page.html), then throw in the [relevant call](http://www.mrtrix.org/developer-documentation/group __spherical__ harmonics.html#gae9a8ddabd463387e10a641f4fcd32858) to generate and store the SH coefficients for a delta function oriented along the desired direction. Should be a piece of cake. 🍰™
Happy to assist if you get stuck… 👍
---
### Author:  [@ThijsDhollander](https://community.mrtrix.org/u/ThijsDhollander)
#### Post date: [November 13, 2018, 11:24pm UTC](https://community.mrtrix.org/t/probabilistic-tracking-on-peaks/2030/3 "2018-11-13T23:24:36Z")
Hi Jakob,
You might want to take a look at [fixel2sh](https://mrtrix.readthedocs.io/en/latest/reference/commands/fixel2sh.html). The input needs to be in the fixel format (which is essentially vectors), but otherwise it does exactly what you’re after.
Cheers,
Thijs
---
### Author:  [@jdtournier](https://community.mrtrix.org/u/jdtournier)
#### Post date: [November 14, 2018, 8:24am UTC](https://community.mrtrix.org/t/probabilistic-tracking-on-peaks/2030/4 "2018-11-14T08:24:36Z")
> [@ThijsDhollander](#):
>
> You might want to take a look at [fixel2sh](https://mrtrix.readthedocs.io/en/latest/reference/commands/fixel2sh.html).
😊 OK, that’s undeniably simpler…
---
### Author:  [@wasserth](https://community.mrtrix.org/u/wasserth)
#### Post date: [November 14, 2018, 9:39am UTC](https://community.mrtrix.org/t/probabilistic-tracking-on-peaks/2030/5 "2018-11-14T09:39:22Z")
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?
---
### Author:  [@wasserth](https://community.mrtrix.org/u/wasserth)
#### Post date: [November 14, 2018, 11:57am UTC](https://community.mrtrix.org/t/probabilistic-tracking-on-peaks/2030/6 "2018-11-14T11:57:09Z")
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.
```auto
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
```auto
fixel2sh amplitudes.nii.gz sh.nii.gz
```
---
### Author:  [@rsmith](https://community.mrtrix.org/u/rsmith)
#### Post date: [November 25, 2018, 6:58am UTC](https://community.mrtrix.org/t/probabilistic-tracking-on-peaks/2030/7 "2018-11-25T06:58:02Z")
> 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](https://github.com/MRtrix3/mrtrix3/blob/master/src/dwi/tractography/algorithms/fact.h), 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.