Dear all,
My name is Luke and I am psychiatrist and neuroimaging researcher working for King’s College London. I am trying to emulate aspects from the processing pipeline detailed in the following paper:
“In vivo assessment of anisotropy of apparent magnetic susceptibility in white matter from a single orientation acquisition” https://doi.org/10.1016/j.neuroimage.2021.118442
I have DWI data (data.nii.gz) which I have generated through a typical FSL pipeline. In line with the pipeline presented in the paper, I need to work out "“The predominant fibre orientation was determined as the angle between the first peak of the orientation distribution function (ODF) and the B0 direction using constrained spherical deconvolution (Tournier et al. 2007: 10.1016/j.neuroimage.2007.02.016).”
I am using the steps below to generate data from the first peak.
mrconvert data.nii.gz dwi.mif -fslgrad bvecs bvals
dwi2response tournier dwi.mif response.txt
dwi2fod csd dwi.mif response.txt fod.mif
sh2peaks fod.mif peaks.mif
mrconvert peaks.mif peaks.nii.gz
sh2peaks -num 1 fod.mif first_peak.mif
mrconvert first_peak.mif first_peak.nii.gz
# Load the peaks NIfTI file
peaks_img = nib.load('/data/project/CARDS/working/007/DTI/peaks.nii.gz')
peaks_data = peaks_img.get_fdata()
# Define the B0 field direction (assuming it's along the z-axis, i.e., [0, 0, 1])
B0 = np.array([0, 0, 1])
# Calculate the angle between each peak and B0 at each voxel
def compute_angle(peak, B0):
dot_product = np.dot(peak, B0)
peak_magnitude = np.linalg.norm(peak)
B0_magnitude = np.linalg.norm(B0)
cos_angle = dot_product / (peak_magnitude * B0_magnitude)
return np.arccos(np.clip(cos_angle, -1, 1)) # Ensure valid range for arccos
# Initialize an array to store the angle values (alpha) for each voxel
alpha = np.full(peaks_data.shape[:3], np.nan)
# Loop through each voxel and calculate the angle for the primary peak
for x in range(peaks_data.shape[0]):
for y in range(peaks_data.shape[1]):
for z in range(peaks_data.shape[2]):
# Extract the peak vector for the voxel (only 3 components in this case)
peak_vector = peaks_data[x, y, z, :]
# Check if the peak vector is not all zeros (to avoid invalid calculations)
if np.linalg.norm(peak_vector) > 0:
# Compute the angle between the peak vector and B0
alpha[x, y, z] = compute_angle(peak_vector, B0)
Could I check if this appears to be the correct way to work out the fibre-to-field angle which represents the orientation of the fibre with respect to B0? Many thanks.
The angle file does not look correct.
My first_peak.mif image looks great (attached below). This is a 4D image with values between -1.5 and 1.7
Is there a simple way that I can generate a map where the values are the angle between peak direction and B0 (the z-axis)?
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.
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:
Dear Donald,
I did not know you are at St. Thomas! I’m finishing up my PhD and working in East London but if we set up any new DTI projects I definitely discuss with you from the outset.
Thanks for explaining the need to normalize the peaks. It looks to me as though there might be an error with the normalization using the code you’ve outlined as the values are not between -1 and 1. The subsequent angle image has quite a few NaN values. Could you double check the code and see if there is an error?
