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.