mrcat DWI images with PE scheme on axis 4 does not work

[maxpietsch]

Currently mrcat concatenates the pe_schemes and throws an error as the output has only half the number of volumes as the created scheme suggests. mrcat: [ERROR] Number of volumes in image "dwi.mif" (300) does not match that in phase encoding table (600).

We probably should check whether both schemes are identical and use only one, warn and strip it otherwise?

[thijsdhollander]

I'm slightly confused... doesn't this mean your inputs then also had half the number of volumes compared to their schemes...? The PE schemes handling is analogous to the gradient scheme handling, right?

[maxpietsch]

$ mrinfo dwi.mif
************************************************
Image:               "dwi.mif"
************************************************
  Dimensions:        100 x 100 x 64 x 300
  Voxel size:        1.51515 x 1.51515 x 1.5 x 1
  Data strides:      [ -2 3 4 1 ]
  Format:            MRtrix
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1          -0           0      -78.41
                               -0           1           0      -54.45
                               -0          -0           1      -39.96
  dw_scheme:         0.019866,0.411523,0.911183,0
  [300 entries]      -0.385352,-0.305637,0.870684,2600
                     ...
                     -0.305057,-0.445414,-0.841752,2600
                     -0.19586,0.857924,0.47498,1000
  mrtrix_version:    3.0_RC3-266-gebbabc8d-dirty
  pe_scheme:         1,0,0,0.1
  [300 entries]      -1,0,0,0.1
                     ...
                     0,1,0,0.1
                     0,-1,0,0.1

Concatenate twice:

 $ ~/mrtrix3_dev [dev]$ mrcat dwi.mif dwi.mif -axis 4 - | mrinfo -
mrcat: [ERROR] Number of volumes in image "dwi.mif" (300) does not match that in phase encoding table (600)

thrice:

$ mrcat dwi.mif dwi.mif dwi.mif -axis 4 - | mrinfo -
mrcat: [ERROR] Number of volumes in image "dwi.mif" (300) does not match that in phase encoding table (900)

mrcat does not abort when the pe_scheme is removed in one of the files

mrcat dwi.mif dwi.mif $(mrconvert dwi.mif -clear_property pe_scheme - ) dwi.mif -axis 4 - | mrinfo -

but then the output does not have a pe_scheme header entry and the diffusion encoding scheme is now 1200 entries long (the image has 300 volumes, 4 volume groups):

************************************************
Image:               "/home/mp14/tmp13/mrtrix-tmp-82MOey.mif"
************************************************
  Dimensions:        100 x 100 x 64 x 300 x 4
  Voxel size:        1.51515 x 1.51515 x 1.5 x 1 x ?
  Data strides:      [ -2 3 4 1 5 ]
  Format:            Internal pipe
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1          -0           0      -78.41
                               -0           1           0      -54.45
                               -0          -0           1      -39.96
  command_history:   mrcat "dwi.mif" "/home/mp14/tmp13/mrtrix-tmp-3SJ78r.mif" "dwi.mif" "dwi.mif" "-axis" "4" "-"  (version=3.0_RC3-51-g52a2540d)
  dw_scheme:         0.019866,0.411523,0.911183,0
  [1200 entries]     -0.385352,-0.305637,0.870684,2600
                     ...
                     -0.305057,-0.445414,-0.841752,2600
                     -0.19586,0.857924,0.47498,1000
  mrtrix_version:    3.0_RC3-51-g52a2540d

[thijsdhollander]

Oh, ok, now I see what I overlooked: axis 4. Interesting to see how dw_scheme does concatenate in that case... so there is indeed a discrepancy between how both end up being handled in this scenario. But (I think) this raises the more general question of what should happen with either dw_scheme as well as pe_scheme when concatenating along any dimension that is not axis 3 (above, I was wrongly reasoning you were talking about axis 3 actually; hence my confusion). We're now up to having about four categories of dimension(s) in MRtrix... in principle, for each operation (and mrcat in particular) we should probably start to think through each case more specifically, and see if it's dealt with well:

• Axes 0-1-2 : spatial
• Axis 3 : volumes
• Axis 4 : volume groups (this one is maybe not generalised enough yet throughout MRtrix, even though there's dedicated keys to switch through volume groups in mrview)
• Axis 5-... : "all the rest" category; still, there should probably be a consistent strategy to deal with those (e.g. in mrcat)

Opinions, @MRtrix3/mrtrix3-devs ?

[Lestropie]

This really needs more detailed information from anyone actually using > 4D data as to what they're actually representing within axis 4 onwards, and therefore whether there is an appropriate strategy for handling dw_scheme and pe_scheme that generalises. mrcat may behave differently to other commands due to testing for axis > 2 rather than axis == 3; but the desired behaviour needs to be explicit before such a 'bug' can be 'fixed'.

[maxpietsch]

I use mrcat -axis 4 to visually compare DWI files (usually with identical dw_ and pe_scheme) in mrview. Hence, the desired behaviour for my use case is stated above: check whether both schemes are identical and use only one, warn and strip it otherwise. The bug label refers to mrcat throwing an exception with a misleading error message.

[Lestropie]

That's valid enough given your particular use case. My concern is that there may be other equally-valid use cases where the expectation, and hence desired behaviour, differs.

For instance, it would seem at least somewhat reasonable to concatenate two 4D images along axis 4, where the diffusion sensitisation is identical between the two volume groups but the phase encoding of the second volume group is reversed with respect to the first. In this case, proper representation of the phase encoding information would require either turning pe_scheme into a 3D matrix, or making an assumption about how pe_scheme may or may not include a concatenation across volume groups, and what the order of concatenation is with respect to axis 3. The alternative is to discard such information; which I suppose for the sake of removing a misleading error message is fine, until such hypothetical time that someone tries to actually use it in such a way.

Just checking whether or not we're in need of a higher-order discussion about how to handle header key-value data in 5D images, as opposed to just fixing the issue as it appears - which should be a matter of both changing this line from: if (axis > 2) { to: if (axis == 3) { , and then having a separate code branch for axis > 3 to catch and erase any key-value data that's not precisely consistent across images.

[jdtournier]

My 2 cents for what it's worth:

There are no relevant commands currently that interpret these fields and are capable of processing 5D data - they all expect 4D data (dwi2tensor, dwi2fod, ...). So the only requirement at this stage is to make sure that whatever we do does not introduce ambiguity in these type of data. On that basis, I support @maxpietsch's suggestion: check they're identical and keep one, or warn and discard both (applies to both dw_scheme and pe_scheme). At least this means that if someone subsequently extracts a 4D subset from the 5D image with a call to mrconvert -coord 4 n, the information is either correct, or not present - never ambiguous.

The practical reality is that this is rarely used functionality, primarily used by us as a means of visualisation. I don't anticipate this will impact anyone's workflow for the foreseeable future. But this issue currently prevents this kind of casual usage (typically to pipe directly into mrview), so we just need a simple fix that will allow this operation without introducing further problems. We can revisit if & when we come up with more complex processing tools that need to operate on 5D data.

Also, the specific uses cases where it might make sense to concatenate along the 5th dimension despite differences in DW or PE schemes are I reckon few & far between. It's just as easy to concatenate along the 4th dimension, as indeed is what all of the relevant tools would expect. I'm not saying there's never going to be a good reason to do this, only pointing out that there are already accepted ways of merging these types of data sensibly using the existing framework.

[Lestropie]

If the problematic use case encountered is not opening Pandora's Box w.r.t. 5D data, then the simple fix is the way to go. But I wanted to be sure that that was indeed the situation. There's probably a large number of other cases where things could potentially go awry (e.g. of the top of my head, I'm guessing dwiextract on a 5D input would likely not do what one would intuitively expect); but if you're not actually trying to handle such data over and above an mrcat | mrview, and are not expecting others to do the same, such hypothetical limitations are not worth frying brains over.

[thijsdhollander]

Happy with @jdtournier's arguments; they make full sense to me. Whatever we do now, or in the future, I reckon, should always honour

applies to both dw_scheme and pe_scheme

They're a kind of property that applies on a same level, as one entry in each one of those schemes goes together with one acquired volume.

Actually, doing things as @jdtournier suggest, results in simply treating the 5th dimension (axis number 4) as the first 3 dimensions (axes 0, 1 and 2), in the context of the DW and PE schemes. I now realise this makes things a whole lot easier to understand (or discuss, document, ...): the schemes simply go together with the 4th dimension (axis 3), without any further complexities. A valid scheme has the same number of entries (lines) as the length of axis 3.

So all we've got to do in this context, is treat axis 4 in the same way as axes 0, 1 and 2. Does e.g. mrcat along either axis 0, 1 or 2 already work correctly? It seems it should be performing the same operation or processing of the dw_scheme and pe_scheme that is being suggested for axis 4.

[Lestropie]

So all we've got to do in this context, is treat axis 4 in the same way as axes 0, 1 and 2. Does e.g. mrcat along either axis 0, 1 or 2 already work correctly?

Not in terms of what's being discussed here.

Specifically for mrcat, the proposed behaviour for concatenation along axis 4 (keep DW / PE schemes if identical across inputs, discard entirely otherwise) should in fact be applied when concatenation occurs along any axis other than 3. As soon as, for a fixed position along axis 3, such sidecar information is not homogeneously applicable to all locations along any other axis (spatial or otherwise), it needs to be discarded.

[thijsdhollander]

Yep, exactly.

There's another little challenge in there though:

keep DW / PE schemes if identical across inputs, discard entirely otherwise

Challenge being mostly for DW schemes probably: to what precision do they have to be identical to be regarded as... "identical"?

[jdtournier]

By the way, this reminds me of a conversation we had a while back about formalising 'per-volume' header entries - I can't find the discussion now, if you do find the relevant issue, please link to it here.

But this clearly relates to this one. To recap, at the time, we'd discussed having a format for entries that flagged them up for treatment as per-volume parameters, which would be handled appropriately by mrcat, mrconvert -coord, etc. Problem being that it would likely require new names for these entries (or at least some different formatting), breaking compatibility with the existing 'standard'. So that would mean instead of having the dw_scheme entries as:

dw_scheme: 0,0,0,0

We might have them as:

dw_scheme@3: 0,0,0,0

which would signal that each entry is expected to correspond to sequential locations along axis 3. There's other ways to do this too, this is just an example to illustrate how it might be done. This would allow custom entries to advertise themselves as per-volume, and the MRtrix3 backend would treat them appropriately.

Not sure how much appetite there might be for this, just thoughts I'd raise it as it might be relevant in the context of how such entries should be handled in a 5D context...

[Lestropie]

Challenge being mostly for DW schemes probably: to what precision do they have to be identical to be regarded as... "identical"?

This is an issue for pe_scheme as well. e.g. See dwipreproc code here for dw_scheme and here for pe_scheme.

... formalising 'per-volume' header entries ...

You could potentially provide such a capability for user-specified header entries without breaking backwards compatibility: Existing keys would simply be added to a "whitelist" of sorts, that would be handled along with any keys that conform to some specified format. But that'd shift the focus toward user-specified fields rather than these, in which case you'd need examples of such to justify the effort. There's also the issue of precision as mentioned above.

[thijsdhollander]

... formalising 'per-volume' header entries ...

While not a bad idea (at all), unless indeed someone has the time/motivation to put into it, I'd not prioritise it at the moment.

So well, trying to go to a conclusion here: I've got the feeling we're all happy with the dw_scheme and pe_scheme going with axis 3 exclusively, and treating any other axes (0, 1, 2, 4 and beyond) in a similar manner? Correct me if I'm wrong of course.

Only issue remaining is the precision thing to compare schemes. Not sure on that one myself... or am I just overthinking it? :thinking:

[jdtournier]

Agree this is low priority, just thought it might be worth pointing out given that we're discussing how we should handle volume-wise entries generically...

This would also suggest that it might be an idea to set up a common code path to handle such entries, which we could then use to handle both dw_scheme and pe_scheme, and later on open up for other entries.

[maxpietsch]

Another argument for automatic handling of axis-specific (user-specified) header entries or at least stripping non-matching entries:

mrinfo a.mif -property image_scale
1.0
mrinfo b.mif -property image_scale
2.0

mrcat a.mif b.mif - -axis 3 -quiet | mrinfo - -property image_scale
1.0

[jdtournier]

😕 what is image_scale supposed to represent? I hope it's not the data's intensity scaling entries, we only support a unique scale + offset per whole image...

[maxpietsch]

I store volume-specific metadata in the header, image_scale is just an exemplary custom-defined header entry. On concatenation, the same would happen to our dwi_norm_scale_factor, lognorm_scale, ...

[jdtournier]

OK 😅

This is going to be tricky to get right though, if I understand correctly. For example: if you concatenate 2 4D images with different numbers of volumes, each with its own lognorm_scale value, what happens? Ideally, it would replicate the entry for all volumes in the first image, and replicate it with the different correct value for the second. So in the one case, we have a single entry, in the latter we have one entry per volume. So this argues that either we should store these values per volume anyway (for consistency), or we adopt the convention that if all values are identical, we only store it once for the whole image (my preference). This is going to require some thought...

[maxpietsch]

Yes, we'd need a mechanism to differentiate between global and volume-wise values. In the case of mrcat * -axis 3 out, it's not too hard though. If a subsequent command expects a single entry but the header has a list of these, then I recon, it should fail. mrcat on the other hand should either duplicate values on the axis specified or at least discard them if they are not the same.

if you concatenate 2 4D images with different numbers of volumes

this is currently not possible as the image dimensions do not match. One of the reasons why I wanted to have amrpad in.mif pad -as reference.mif -nd out.mif. For the above to work, mrpad would need to append empty entries to the padded volumes.

[jdtournier]

this is currently not possible as the image dimensions do not match

Was talking about concatenating along the volume axis - but yours is a good point too...

[maxpietsch]

Was talking about concatenating along the volume axis - but yours is a good point too...

:nerd_face: Turns out my brain is zero-based.