Clarify and correct b-vector format and orientation

[jeremie-fouquet]

Description
It would be good to:

1. Clearly identify the output format that brkraw is using for bvec and bval. For now, it seems to be the FSL format, based for instance on this page: https://mrtrix.readthedocs.io/en/latest/concepts/dw_scheme.html#fsl-format.
2. If the FSL format is chosen, make sure the vectors in the bvec file are in the image space.
3. Update the documentation to reflect this.

Suggested solution
If the FSL format is chosen, the PVM_DwGradRead, PVM_DwGradPhase, and PVM_DwGradSlice parameters specify the b-vectors in the read/phase/slice directions. In theory, we would only need to transform them into the radiological convention, according to https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/FAQ#What_conventions_do_the_bvecs_use.3F. This page also provides a bit more context: https://www.mrtrix.org/2016/04/15/bug-in-fsl-bvecs-handling/.

Alternative solutions

• The PVM_DwGradVec parameter gives the b-vector in the scanner space (according to PV6.0.1 and PV360 documentation, in the "x,y,z coordinate system"). I think applying the inverse image space->scanner space transform to these b-vectors should in theory be sufficient to get the b-vectors into the image space. We would then also need to transform them into the radiological convention.
• Decide to go for a different format, like the .b format from mrtrix3 (https://mrtrix.readthedocs.io/en/latest/concepts/dw_scheme.html#mrtrix-format)

Additional information

• This seems to also be an open question for the BIDS format:
  -[ENH] Optional alternative to bval/bvec bids-standard/bids-specification#349
  -Incomplete description of DWI bvecs format bids-standard/bids-specification#243

[Nerdynien]

Description It would be good to:

1. Clearly identify the output format that brkraw is using for bvec and bval. For now, it seems to be the FSL format, based for instance on this page: https://mrtrix.readthedocs.io/en/latest/concepts/dw_scheme.html#fsl-format.
2. If the FSL format is chosen, make sure the vectors in the bvec file are in the image space.
3. Update the documentation to reflect this.

Suggested solution If the FSL format is chosen, the PVM_DwGradRead, PVM_DwGradPhase, and PVM_DwGradSlice parameters specify the b-vectors in the read/phase/slice directions. In theory, we would only need to transform them into the radiological convention, according to https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/FAQ#What_conventions_do_the_bvecs_use.3F. This page also provides a bit more context: https://www.mrtrix.org/2016/04/15/bug-in-fsl-bvecs-handling/.

Alternative solutions

• The PVM_DwGradVec parameter gives the b-vector in the scanner space (according to PV6.0.1 and PV360 documentation, in the "x,y,z coordinate system"). I think applying the inverse image space->scanner space transform to these b-vectors should in theory be sufficient to get the b-vectors into the image space. We would then also need to transform them into the radiological convention.
• Decide to go for a different format, like the .b format from mrtrix3 (https://mrtrix.readthedocs.io/en/latest/concepts/dw_scheme.html#mrtrix-format)

Additional information

• This seems to also be an open question for the BIDS format:
  -[ENH] Optional alternative to bval/bvec bids-standard/bids-specification#349
  -Incomplete description of DWI bvecs format bids-standard/bids-specification#243

@jeremie-fouquet do you have a current workaround/solution to generate the appropriate bval/bvec files in the correct space to continue an analysis in FSL? this issue has been haunting me for a while!

[jeremie-fouquet]

My current workaround is a homemade script that reorders the x-y-z components of the .bvec file to match my acquisition orientation (plus I don't rotate the acquisition FOV so I don't have to apply extra rotations). I confirmed the reordering was right by checking the orientation of the tractography fibers in a white matter tract.

This script is specific to our acquisition orientation. A general solution would be much better...

[araikes]

@jeremie-fouquet et al,

I'm testing the edits you made in #96 and notice that the b-vectors that are produced are not unit norm vectors. As indicated here (https://dicomifier.readthedocs.io/en/latest/diffusion/bruker/index.html), the DwGradVec field is the gradient amplitudes relative to the maximum rather than unit vectors.

The consequence of this is that like MRtrix will attempt to scale the b-values by the squared amplitude of the vectors (https://mrtrix.readthedocs.io/en/latest/concepts/dw_scheme.html#b-value-scaling). I know that FSL's DTIFIT will also normalize the vectors and I'm not sure how it treats the b-values internally when doing so.

It would probably be beneficial produce unit norm vectors as the default output so that no unexpected b-value scaling happens behind the scenes.

Note: I used dicomifier (see link above) to convert some data that I'm working on right now as well. Their method produces bvecs that are unit norm and correctly ordered relative to the image space and oriented (confirmed with tractography) without any additional manipulation. However, the bvecs are slightly different (+/- 0.05) and the NIFTI has a different orientation compared to brkraw's output. Additionally, dicomifier doesn't produce BIDS-compliant JSONs, so brkraw is preferable.

Currently, I use MRtrix's dwigradcheck plus a visual assessment to reorient the bvecs and ensure radiologic convention.

[jeremie-fouquet]

Good point, @araikes. Since pull request #109, b-vectors should be normalized as a temporary fix until we develop a coherent scheme to output them in a space we agree on.

[araikes]

@jeremie-fouquet and @dvm-shlee,

Please note relevant discussions and specification issues:

1. https://bids-specification.readthedocs.io/en/stable/appendices/coordinate-systems.html#introduction
2. https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/01-magnetic-resonance-imaging-data.html#required-gradient-orientation-information
3. https://neurostars.org/t/rotate-bvecs-when-warping-dwi-to-t1w/26136/11
4. https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FDT/FAQ#What_conventions_do_the_bvecs_use.3F

[dvm-shlee]

@araikes @jeremie-fouquet Thank you both for your outstanding contributions. To further enhance our handling of diffusion-specific issues, I've implemented a few changes:

Separated b-Vector Handling: I've isolated the handling of the b-vector format into a separate class to facilitate better troubleshooting and precision. Diffusion helper module

Orientation Helper Update: The Orientation Helper has also been separated to simplify future updates. Orientation Helper to simplify updates. Assuming the b-vector adheres to the Bruker scanner convention (readout: x, phase: y, slice: -z), applying an affine transformation using VisuCoreOrient should prove effective, however,

For greater converted output consistency across PV version or sequence, I would like to request anyone's support to collecting dataset to test and validate Bvec data as below

1. Ground Truth b-Vector and FA Image: We should create a manually corrected, ideal ground truth b-vector, accompanied by an RGB-mapped FA image calculated from this vector.
2. Sample Datasets for Testing: I plan to provide sample datasets for each version (5.1, 6.1, 7, 360, if feasible) for thorough testing. Contribution guidelines for these datasets will be forthcoming at the BrkRaw dataset repository. Additionally, the next update will include a CI routine to ensure all data shared there are converted correctly, with priority given to correcting any detected issues.

Thank you!

[araikes]

Hi @dvm-shlee,

We've been working toward a ground-truth diffusion phantom and I'm following up to see where that project is at. If I'm able to share the data I will, otherwise I'll put together a reasonably verbose set of outputs for evaluation from PV360.

I will say this: The ordering of the b-vectors is extremely dependent on the output orientation of the NIFTI. For example, from our mouse data, I've determined the following (in PV360, brkraw 0.3.11 converted datasets):

import argparse
import nibabel as nb
import numpy as np
import json
from dmriprep.utils.vectors import DiffusionGradientTable


# Load the gradients and image together using DMRIPrep. This uses DIPY notation and formatting of the gradient table.
gradients = DiffusionGradientTable(dwi_file = args.input, 
                                   bvecs = args.bvec,
                                   bvals = args.bval, 
                                   b_scale = False,
                                   bvec_norm_epsilon = 0.1)

# Reorder vectors to match image as output by brkraw 
# Ordering is based on image orientation and pulse sequence (for SPR+ as the default)
# Use QSIPrep notation
if orient == "LPI+":
    transform = np.array(([1,1],[0,1],[2,1]))
elif orient == "IPL+":
    transform = np.array(([0,1], [1, -1], [2,1]))
elif orient == "RPS+":
    transform = np.array(([1,1], [0, 1], [2,1]))
elif (orient == "SPR+") and ('uaDtiEpi' in pulseseq):
    transform =  np.array(([0,1], [1,1], [2,-1]))
else:
    transform = np.array(([0,-1], [2,1], [1,1]))
    
output_bvecs = np.zeros_like(gradients.bvecs.T)
for this_axnum, (axnum, flip) in enumerate(transform):
    output_bvecs[this_axnum] = gradients.bvecs.T[int(axnum)] * flip

Beyond that, if I'm then reorienting to a common output orientation (in this case LPS):

new_axcodes = ("L", "P", "S")

if not input_axcodes == new_axcodes:
    # Re-orient
    input_orientation = nb.orientations.axcodes2ornt(input_axcodes)
    desired_orientation = nb.orientations.axcodes2ornt(new_axcodes)
    transform_orientation = nb.orientations.ornt_transform(input_orientation, desired_orientation)
    reoriented_img = input_img.as_reoriented(transform_orientation)
    output_bvecs_lps = np.zeros_like(output_bvecs)
    if (orient != "LPI+"):
        for ax, (axnum, flip) in enumerate(transform_orientation):
            output_bvecs_lps[ax] = output_bvecs[int(axnum)] * flip
    else:
        transform = np.array(([0,1], [1,-1], [2,1]))
        for ax, (axnum, flip) in enumerate(transform):
            output_bvecs_lps[ax] = output_bvecs[int(axnum)] * flip
else:
    print("Already correctly oriented.")

So there's a pretty wide range of transforms required that don't inherently align with the way things are currently being written. However, I do understand that there's updates on the immediate horizon that may address this heterogeneity as well.