Support for b-tensor encoding

[bjeurissen]

With this issue, I want to kick of the discussion about how to store b-tensor encoding information.

The b-tensor as a natural extension of the b-vector/value combination

The most natural format seems to be to simply store the 6 unique elements of the b-tensor. This has the benefit that it allows to represent arbitrary b-tensor shapes, whilst still providing relatively straightforward access to the directional information (through the eigenvectors), b-value (trace), and b-shape (proportions of eigenvalues). The shell-clustering can operate directly on the vector of sorted b-tensor eigenvalues instead of simply the b-value as we do now.

This approach would add a lot of flexibility, without too much added complexity as much of the tensor math that we have in place for DTI can simply be reused. Moreover, given that there are already DICOM tags in place to store the full b-matrix (https://dicom.innolitics.com/ciods/mr-spectroscopy/mr-spectroscopy-multi-frame-functional-groups/52009230/00189117/00189601), this is also how manufacturers could start providing this information once b-tensor encoding makes it into product sequences.

Some open questions/issues:

Tensor parametrization

Diffusion tensors are currently stored in MRtrix as [XX YY ZZ XY XZ YZ]. I assume we would stick to this order and scaling? The b-tensor community (e.g. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4916005/#APP1, https://github.com/markus-nilsson/md-dmri, https://github.com/bjeurissen/ciwlls) seems to favour a variation of Mandel Notation [XX YY ZZ XY*sqrt(2) XZ*sqrt(2) YZ*sqrt(2)], which simplifies some of the tensor mathematics (e.g. inner and outer products of tensors can be performed by inner and outer products of their 6x1 representations). Looking to DICOM, they seem to suggest [XX XY XZ YY YZ ZZ], although each entry has its own tag, so it might not suggest and order after all (https://dicom.innolitics.com/ciods/mr-spectroscopy/mr-spectroscopy-multi-frame-functional-groups/52009230/00189117/00189601). Note that regardless of which parametrization we use, we should probably use it across the board (i.e. also for diffusion tensors).

Conversion to and coexistence with dw_scheme

How do we handle precedence and conflicts? Unless we amend dw_scheme with e.g. b-tensor anisotropy, only LTE schemes can be converted to dw_scheme. What about b-tensors with small secondary and tertiary eigenvalues?

How does this fit into the bigger picture of advanced contrast mechanisms?

E.g., varying TE, small, big delta, funky waveforms, per slice encoding, ... This is an interesting source of inspiration: https://arxiv.org/pdf/2103.14485.pdf, but it could be taking things a bit too far at this point... EDIT: See #2323

Gradient polarity is lost

When using the b-tensor, the polarity of the gradient can no longer be recovered, which does not play well with EDDY.

[Lestropie]

Prior initial listing in #2323.

1. I was never a fan of storing only the b-tensor anisotropy as a fourth / fifth column rather than the full b-tensor.
2. In other contexts we have chosen to stick with the DICOM standard rather than alternatives, and as you say it would mean consistency with existing tensor handling code, so I would advocate for that rather than Mandel.
3. I think I would prefer to use a different header key rather than augmenting dw_scheme. The current interpretation of dw_scheme is that the first three columns are spatial 3-vector, fourth column is b-value, then optional additional columns are "any other information you'd like to store along with each volume". Trouble is, there would be no way to discriminate between "an extra three columns in order to represent the b-tensor" and "three additional columns of unknown content". So processes like shelling would need to look first for the new key (eg. btensor_scheme), and then for dw_scheme, and also throw an error if both are present and they're incompatible.
4. Variable TE can be encoded in EchoTime as a vector of equal length to the number of volumes. While this isn't part of the BIDS specification AFAIK (for eg. relaxometry they encode each echo as a separate image with its own JSON), it wouldn't surprise me if it becomes the case at some point. This is already supported in MRtrix3 at DICOM import for the multi-TE non-DWI data I've encountered.
5. For more advanced DWI contrasts, I think that there would need to be a mature plan in place to be able to do something meaningful with that additional information before investing the effort in supporting it in the back-end. And if we were to pursue point 3 above and then encounter this need later, it could exist as a third unique key, with precedence at read time applied accordingly.
6. Is the loss of gradient polarity at the point of DICOM generation? I'd have thought that the tensors could be antipodally flipped.
7. I would have thought that there'd be a lot more about these data that would not play well with eddy w.r.t. such data over and above just gradient polarity? I wonder if it's worth pinging FMRIB to see if they have any specific plans about how to start supporting their data, since their bvecs / bvals representation breaks down in the same way as our own.

[bjeurissen]

6. Is the loss of gradient polarity at the point of DICOM generation? I'd have thought that the tensors could be antipodally flipped.

What I meant was that the b-matrix representation itself (regardless of whether we read it from DICOM or create it ourselves) cannot completely replace dw_scheme as it lacks the ability to store the gradient polarity due to its [x; y; z]*[x y z] nature. This is fine for modelling, where tensors can indeed be flipped, but not for eddy as the polarity affects the nature of the EC distortions.

7. I would have thought that there'd be a lot more about these data that would not play well with eddy w.r.t. such data over and above just gradient polarity? I wonder if it's worth pinging FMRIB to see if they have any specific plans about how to start supporting their data, since their bvecs / bvals representation breaks down in the same way as our own.

Indeed, eddy does not support non-LTE encoding at the moment. My comment about eddy was related to performing eddy on LTE data when only the b-matrix is available. This breaks down as you no longer have access to the gradient polarity.

[bjeurissen]

Based on the above, would like to propose:

1. Keep dw_scheme key as is
2. Add new btensor_scheme key that takes precedence over dw_scheme and use this for encodings we cannot describe at the moment, i.e. isotropic diffusion encoding and b-tensor encoding
3. Add 3D shell clustering for btensor_scheme based on the 3 btensors eigenvalues
4. Add reorientation logic for btensor_scheme (as we have for dw_scheme)
5. For images that only have the btensor_scheme and not the dw_scheme, prevent using eddy in dwifslpreproc as it is not supported correctly anyway.

Some additional points:

1. On the DICOM side, ideally we would use https://dicom.innolitics.com/ciods/enhanced-mr-color-image/enhanced-mr-color-image-multi-frame-functional-groups/52009229/00189117/00189075 and create a correct btensor_scheme for cases where directionality is ISOTROPIC or BMATRIX. If diffusion directionality is DIRECTIONAL, do we also create btensor_scheme or only dw_scheme?
2. To represent the btensor, do we follow the order implied in DICOM standard (XX,XY,XZ,YY,YZ,ZZ) or use the MRtrix convention for diffusion tensors (XX, YY, ZZ, XY, XZ, YZ)? And if we go for the former, do we adjust how diffusion tensors are represented as well?

[bjeurissen]

Additionally, it would be good to add a check for axial symmetry of the btensor_scheme. Not having axial symmetry is fine in general and will work with commands such as dwi2tensor, but it will cause problems for commands that inherently assume axial symmetry of the encoding such as dwi2response, dwi2fod, ...

[dchristiaens]

I'm also very interested to have support for multidimensional encoding - happy to help out.

Based on the above, would like to propose ...

Sounds like a good plan to me. Reorientation of the B-tensor scheme may prove tricky to get for encoding tensors that are not axially-symmetric, I expect. Regarding preprocessing, once we have shell clustering for B-tensor encoding you can use SHARD for motion correction - I do this regularly by manually "hacking" dw_scheme.

Re your questions:

On the DICOM side, ... If diffusion directionality is DIRECTIONAL, do we also create btensor_scheme or only dw_scheme?

Only dw_scheme, I'd say. Duplication risks inconsistency.

To represent the btensor, do we follow the order implied in DICOM standard (XX,XY,XZ,YY,YZ,ZZ) or use the MRtrix convention for diffusion tensors (XX, YY, ZZ, XY, XZ, YZ)? And if we go for the former, do we adjust how diffusion tensors are represented as well?

I would store the B-tensors in accordance with the MRtrix convention for DTI. Changing the DTI representation now is too big a change for users, comparable to the introduction of the orthonormal SH basis in MRtrix3. And having different tensor representations for the B-tensor and the diffusion tensor is just ugly...

[Lestropie]

I'm happy to go with MRtrix ordering of diffusion tensor components. Just make sure it's documented.

There will be all sorts of weird edge cases that you'll need to consider to make the handling robust; eg. headers that have somehow obtained both btensor_scheme and dw_scheme (ideally they should be mutually exclusive); performing a volume concatenation of two images where one has btensor_scheme and the other has dw_scheme. I can probably come up with various ways to try to break an implementation, but there might not be too much benefit in thinking about it until after something functional is implemented.

Regarding dwifslpreproc: While I've not read anything specific anywhere, I've a suspicion that it may be possible to "hack" eddy by just filling bvals with unique identifiers for each group of volumes, whether that grouping is influenced by b-value / TE / tensor shape. On the topic of polarity, is it definitely the case that for sensitisation gradients that are not STE or LTE, polarity of that sensitisation is lost in the DICOM encoding? I'd have thought that that would be preserved, such that generating a hack bvecs for eddy could contain something that would give some reasonable information about relationships of eddy currents between volumes.

[bjeurissen]

I'm happy to go with MRtrix ordering of diffusion tensor components. Just make sure it's documented.

Based on the feedback of @dchristiaens, this is what I have started implementing. As you say, many of the edge cases and their consequences will probably only become totally clear once we have something functional implemented.

On the topic of polarity, I think the b-tensor representation simply cannot encode it.

As a simple example, both

g = [1 0 0]
b = 1
B = g'*g*b

and

g = [-1 0 0]
b = 1
B = g'*g*b

will produce the same b-tensor:

B =
     1     0     0
     0     0     0
     0     0     0

On the topic of eddy: I think it is just not equipped to deal with anything other than LTE. It might be possible to hack bvals and bvecs to get reasonable results on b-tensor encoding or data with multiple TE's but in my experience this does not even work well for something as relatively simple as LTE with multiple TE's. For me it is unclear exactly how eddy uses the bval/bvec pair to predict eddy current distortions. Does it only rely on the bvec to relate distortions within a single shell or does it also use bval to relate distortions between shells? Was anyone able to find the actual formulas used? In any case, I think for the time being, the easiest solution is to only allow dwifslpreproc to be run on data with dw_scheme.

[bjeurissen]

On the topic of DICOM: did anyone come across DICOM data with the "Diffusion Directionality Attribute" (0018,9075) set to BMATRIX? I think I have only seen it in Bruker generated data...

[bjeurissen]

Some concerns about how this will be supported by different commands:

• Spherical b-tensors only:

  • dwi2adc --> currently we don't properly encode trace weighted images, we could do this with btensor_scheme

• Linear b-tensors only (or at least close to linear):

  • dwifslpreproc

• Axially symmetric b-tensors:

  • dwi2fod
  • dwi2response --> EDIT suggested by @jdtournier: might need to be reworked to be modular (separate voxel selection and direction selection from the actual calculation of the response using amp2response)
  • amp2response --> EDIT: monotonicity constraint has to take into account b-delta (use sign of b-delta)
  • amp2sh
  • sh2amp
  • sh2response
  • tckgen nontensor
  • tckglobal

• Any b-tensors:

  • dwi2mask --> only has to be able to extract shells?
  • dwi2tensor --> naturally works with b-tensors
  • dwibiascorrect --> only has to be able to extract b=0
  • dwicat --> has to be able to extract b=0, btensor_scheme has to be updated
  • dwidenoise --> does not use gradient info
  • dwiextract --> requires specification of 3 b-tensor eigenvalues per shell, instead of a single b-value
  • dwigradcheck --> requires b-tensor version of sign and axis flips
  • dwinormalise --> only has to be able to extract b=0?
  • dwishellmath --> only has to be able to cluster shells
  • tckgen tensor --> naturally works with b-tensors
  • mrinfo --> btensor_scheme support including shell clustering
  • mrcat --> btensor_scheme has to be updated
  • mrconvert --> btensor_scheme has to be updated
  • mrdegibbs --> not affected
  • mrtransform --> has to reorient btensor_scheme, not hard to do

Any command that allows shell specification will have to support specifying the three b-tensor eigenvalues in addition to just a single b-value --> will get a bit messy

[Lestropie]

• How common are non axially symmetric b-tensors? If they're almost never used, there's an argument for having a command-line interface that selects volume groups based on the scalar shape anisotropy measure.

• There are some complications that arise just from the nature of volume groups having more than one attribute by which they could be filtered, ie not anything specific to the nature of b-tensors. Say for instance you have four volume groups, say LTE and STE at two different b-values. You need to be able to select both LTEs, or both low b-values, or just the low b-value LTE group, or all but the low b-value LTE group. One way to do this is to have dwiextract apply the user's selection criteria as an OR operation. Ie. A volume group only has to satisfy either requested b-value or requested shape to be included in the output. An AND operation can be achieved by piping one dwiextract call into another. Another corollary here is that, if the command-line interface for selection of specific volume groups becomes too cumbersome, then having that interface duplicated across all DWI-related commands may be clumsy. It may instead be preferable to defer that complexity exclusively to dwiextract, and rely more heavily on command-line piping to feed data subsets to other commands. Not necessarily advocating for it up front, just establishing that it's a possibility.

• dirstat currently missing; that will be important for demonstrating the nature of the groupings.

• Indeed these enhancements may expose that the command-line options for mrinfo for accessing DWI grouping info may be in the wrong place.

• RE dwifslpreproc: My understanding right now is that, contrary to what is described in the relevant journal article, eddy treats each "shell" (/ volume group) independently, ie. there's no relationship imposed on the signal relationship between shells. This is why they have the "Post-Eddy Alignment of Shells (PEAS)" mechanism. However the exact nature of the model relating bvecs to eddy current distortions I am not sure.

[bjeurissen]

examples of the suggested syntax to specify shells in the (axially symmetric) b-tensor universum by @jdtournier:

• all b-zeros, all PTE b=1000, all LTE b=1000:

  -shells 0,1000@-0.5,1000@1

• all LTE shells:

  -shells @1

[bjeurissen]

After discussion with @daljit46, I think it would be nice to try to approach this in a more generic way.

It would be nice to have support for generic per-volume attributes of at least these forms:

• scalar: e.g. TE, TR, Flip angle, big delta, small delta, b-value, b-shape, …
• vector: e.g. gradient direction, phase encoding direction, b-tensor elements
  • needs a flag to indicate whether and how to update on transform
  • needs a distance measure that can be used for clustering into shells

Together with correct generic handling of this in mrconvert and mrcat. And with a generic mechanism to apply updates on spatial transform and a generic mechanism to allow selection and clustering on an arbitrary list of attributes.

Since most of the difficulties with supporting btensor_scheme would be covered like that, it would avoid having to repeat this exercise for variable TE, TR, ...

[Lestropie]

We've had at least one long discussion in the past on the topic of "generic axis-linked metadata". But I think it's all undocumented.

1. If the desire is to have multiple metadata fields that can vary per volume, that alone shouldn't be too difficult to expand. There's already been centralisation of the responsible functions, just need to increase the list of known fields that need to be manipulated in that way, and ensure that all commands (not just mrconvert / mrcat) are invoking those centralised functions rather than custom code.

2. Where it gets trickier is if you want any metadata field to be able to also nominate the image axis along which it varies. The code for manipulating such data isn't actually too difficult as long as it's properly centralised; but there needs to be a standardised way to encode such alongside the metadata.

3. One of the prospects I had raised when this came up previously was having an updated MIF/MIH format where the header information is a JSON rather than a std::map<std::string, std::string>. In addition to image axis indices being typed, it would also allow nested dictionaries, so eg. a metadata field could have a longform name, units, associated axis index, etc.. Technically not requisite for the exercise, but it would solve a couple of projected problems.

4. How crucial is it for your current work to have a totally generic solution in place? If everything is axis 3, and you just need proper handling of header merge, concatenation, and volume subset selection, for a relatively short list of fields (some of which already exist), going all-out might be an unnecessarily high barrier. I would at least start with a separate Issue listing where we can discuss that separately to the nuances of b-tensor encoding metadata.

[bjeurissen]

If the desire is to have multiple metadata fields that can vary per volume, that alone shouldn't be too difficult to expand. There's already been centralisation of the responsible functions, just need to increase the list of known fields that need to be manipulated in that way, and ensure that all commands (not just mrconvert / mrcat) are invoking those centralised functions rather than custom code.

I think this would get us very far already for a lot of use cases.

Where it gets trickier is if you want any metadata field to be able to also nominate the image axis along which it varies. The code for manipulating such data isn't actually too difficult as long as it's properly centralised; but there needs to be a standardised way to encode such alongside the metadata.

That would be nice. Already fixing the previous point, might be a good warmup for that.

How crucial is it for your current work to have a totally generic solution in place? If everything is axis 3, and you just need proper handling of header merge, concatenation, and volume subset selection, for a relatively short list of fields (some of which already exist), going all-out might be an unnecessarily high barrier. I would at least start with a separate Issue listing where we can discuss that separately to the nuances of b-tensor encoding metadata.

I think this would indeed cover already the types of data that people can come up with using stock sequences.

[Lestropie]

The intention of the current code base is to manage the key header operations as follows:

• "Concatenation": https://github.com/MRtrix3/mrtrix3/blob/3.0.4/core/header.cpp#L746

  • Special handling for if the metadata field varies along the same axis as the axis of concatenation (works for DWI & phase encoding along axis 3; may not be currently working for slice timing along slice encoding axis)
  • Any other metadata needs to be "merged" (see below), including potentially having some tolerance on some numerical fields (the slice timing vector has been a problem here previously, but conceptually some other fields may require some tolerance also)

• "Merge": https://github.com/MRtrix3/mrtrix3/blob/3.0.4/core/header.cpp#L106 Any time that an output image is being produced from more than one input image, for any operation other than a concatenation.

• "Slicing": Ideally the header manipulation functionality should not be occurring internally within eg. dwiextract / mrconvert. Instead, those commands should be using their own bespoke logic to determine the set of indices along each axis that are to be included in the output image, and providing that information to Adaptor::Subset, which should (doesn't currently) be making the requisite modifications to header metadata: https://github.com/MRtrix3/mrtrix3/blob/3.0.4/core/adapter/subset.h