Hi @lukevano,
To answer this (and your other related post): yes, that approach can work, but in contrast to using the tensor eigenvectors (with the -modulate none
option), the direction vectors produced by sh2peaks
are scaled to the amplitude of the fODF along that direction – they are not unit vectors, so they need to be normalised to unit length before computing the angle using the acos
function.
This ought to do the trick:
mrconvert data.nii.gz dwi.mif -fslgrad bvecs bvals
dwi2response tournier dwi.mif response.txt
dwi2fod csd dwi.mif response.txt fod.mif
sh2peaks -num 1 fod.mif first_peak.mif
# addition to normalise peak vector:
mrmath first_peak.mif norm -axis 3 first_peak_amplitude.mif
mrcalc first_peak.mif first_peak_amplitude.mif -div first_peak_norm.mif
# then use normalised peak for subsequent angle calculation:
mrconvert first_peak_norm.mif -coord 3 2 z-component.mif
mrcalc z-component.mif -acos angle.mif
mrconvert angle.mif angle.nii.gz
I also note you can simplify all this and avoid the need to store so many intermediate files if you want. This relies on the use of Unix pipes and the fact that most MRtrix commands natively support NIfTI anyway:
dwi2response tournier -fslgrad bvecs bvals data.nii.gz response.txt
dwi2fod csd dwi.mif response.txt - | sh2peaks - -num 1 peak1.mif
mrmath peak1.mif norm -axis 3 - | mrcalc peak1.mif - -div - | mrconvert - -coord 3 2 - | mrcalc - -acos angle.nii.gz
Finally, if you want the answer in degrees rather than radians, you can add this to the final mrcalc
call:
... | mrcalc - -acos 180 pi -div -mult angle.nii.gz
Happy to discuss all this face to face if you’re ever in St Thomas, by the way!
All the best,
Donald