Mrtransform does not output same data strides as -template volume

[Romain_Valabregue]

Hello

My appologies if this has alreay been discuss (but I could not find it)

I notice a confusing behaviour of mrtransform when using with -template option.
if the template has a data stride of -1 the output will have a data stride of 1

This induce bug (I had hard time to figure out) : when you reslice a volume to a template reference to expect to have the exact same voxel order, so that you can use your dummy matlab nifti reader will do the right job ( with only taking the voxel value and reoder them in 3D with volume dimension)
If the data stride are different the 3D matrix will no more be coherent. I agree that a good nifit reader should take it into account, but on the other hand after a reslice you expect to do no have to care about it …

So in short I think it will be better (an more consistent) if mrtransform with -template option would output the exact same voxel order.

Many thanks for providing this great soft

Romain

[jdtournier]

I can see the argument for that… Currently, the strides of the output image are based on the strides of the input image. It would be trivial to set them according to the strides of the template, we’d just need to add something like:

Stride::set (output_header, template_header);

at line 254 or so in cmd/mrtransform.cpp.

But I’m not convinced that would always be a good idea. The specific use case I have in mind is specifying a 3D image as template when transforming a 4D dataset, when the input dataset was laid out in volume-contiguous order (i.e. strides like 2,3,4,1): the approach above would immediately remove the volume-contiguous ordering, which will in general not be great (performance implications, for instance). This is unlikely to be an issue for NIfTI since it doesn’t support volume-contiguous storage, but it is an issue for the .mif format. Handling this well is possible, but would require a bit of thought…

Also, there is an argument to be made for preserving the strides of the input image. I’m not convinced I’m 100% sure what the right thing to do is here, when there are two valid sources of information for how to set the strides.

Maybe we could add a -stride option to mrtransform, would that help? Another possibility would be to allow the -stride option to accept an image as argument, and we could then take the strides directly from its header, a bit like how we handle the -include and -exclude options in tckgen? That might be a nice addition in general, I mind consider that in general everywhere the -stride is accepted…

[rsmith]

The specific use case I have in mind is specifying a 3D image as template when transforming a 4D dataset …

The -template option is specifically for defining the voxel grid; any additional dimensionality in this image is irrelevant in the context of its use with this option. So this could potentially be handled something like follows:

• Relative strides between the three spatial axes are determined by image provided via -template option.

• Strides of additional axes relative to one another, and relative to the 3 spatial axes as a group, are determined by input image.

• If the input image has a non-spatial axis with a stride that ‘bisects’ the strides of the spatial axes (e.g. 1,3,4,2) … burn it all down

---

Of course, you can always avoid such ambiguity by explicitly performing a trick like:

mrtransform input.mif -template template.mif - | mrconvert - output.mif -stride $(mrinfo template.mif -strides)

[jdtournier]

Yes, that’s more or less what I had in mind when I said:

@jdtournier

Handling this well is possible, but would require a bit of thought…

But I have reservations about this:

@jdtournier

Also, there is an argument to be made for preserving the strides of the input image. I’m not convinced I’m 100% sure what the right thing to do is here, when there are two valid sources of information for how to set the strides.

And right now, I think the better way to handle this is to leave the strides as-is by default (and document this in the help page), add the -strides option to mrtransform, and allow that option to take an image as input – pretty much exactly what we now have in this pull request. This would allow this kind of usage:

mrtransform in.mif -template template.mif -strides template.mif output.mif

I guess that would be a decent solution?

[Romain_Valabregue]

Hi

I am not sure to fully understand your point, I guess this is because I am not use to the mif format (and volume contiguous order)

My use of mrtransform -template is to reslice a volume in the template space in order to be do further calculation (mrcalc marstat or with fslmaths ect …)

For me this is the principale use of mrtransform (but I may be wrong). In this case it make sens to choose the default option to write the same stride as the template. I would add an other option like " -keep-stides " to specify if you want to preserve the input strides.

@jdtournier

I guess that would be a decent solution?

Yes this is a decent solution (although not the best one in my point of view). but it will do the job, so I am fine with it !

@rsmith

The -template option is specifically for defining the voxel grid; any additional dimensionality in this image is irrelevant in the context of its use with this option. So this could potentially be handled something like follows:

I do not get your point there, why this only define the voxel grid ? the voxel order should be part of it. (and for the reslice to be correct it also have to take into account voxel dimension and volume orientation ).

Following is just examples, to show what problem it causes.

For most of the software it is important to have the exact same voxel order in order to compute 2 images.
I check and this is not the case for mrstat (and I guess for other mrcalc utility). which mean you do not have to take care if the stride is the same, mrtrix will take it into account to keep the 2 volumes coherent and all computations will be correct. Unfortunately the correct use of the nifti header is very particular to mrtrix .

For instance this is not the case with fsl utilities. I do not understand why they handle the first dimension differently. Actually if you have one image with 1 2 3 and the other with -1 2 3, fslmaths with take it into account (and the result is correct) but not if you have 1 2 3 and 1 -2 3 (the bad new it that fslmaths (or fslstat) does not complain, just the results is bad). if you have 1 2 3 and 2 1 3 then fslmaths will complain " image must be the same size" ) → very border line choices …

So in short I thing to have nifti files with different voxel ordering is not a good idea, (unless you only use the best software ( I mean mrtrix of course !). But I still have old habbit with fsl utilities

…

[jdtournier]

@Romain_Valabregue

I am not sure to fully understand your point, I guess this is because I am not use to the mif format (and volume contiguous order)

Yes, this is a feature that is specific to the MRtrix3 mif format, and not found in NIfTI. What volume-contiguous storage allows is to have for example all the DWI values for a given voxel co-located on disk and/or in RAM, which makes it much easier for the system to stream all these data in one go from the disk into RAM and then into the CPU - no need to seek to different locations, and the CPU cache can be used more efficiently. So this is particularly advantageous in applications where 4D datasets need to be processed voxel-wise. For this reason, you’ll note that applications like tckgen will explicitly pre-load the data into RAM if the data are not already stored in volume-contiguous format (it does this to re-order the data to make it volume-contiguous): during tracking, every step requires all the SH coefficients for the voxels near the relevant spatial location (those used in the interpolator), so helping the system to load them efficiently makes a fair difference to performance.

You can avoid this preload if the data happen to already be in the right order (and in 32-bit floating point format) – the application can access the data directly via memory-mapping in this case. This is why for instance dwi2fod will produce its output in 32-bit floating-point format in volume-contiguous order: it’s immediately in the right format for tckgen to use.

@Romain_Valabregue

My use of mrtransform -template is to reslice a volume in the template space in order to be do further calculation (mrcalc marstat or with fslmaths ect …)

Yes, I appreciate that. However, this is not a problem for mrcalc, or mrstats, or any other MRtrix3 command: the whole MRtrix3 backend is designed to handle this transparently. See here for a description of what happens.

But it will indeed be a problem for fslmath…

@Romain_Valabregue

For me this is the principale use of mrtransform (but I may be wrong). In this case it make sens to choose the default option to write the same stride as the template. I would add an other option like " -keep-stides " to specify if you want to preserve the input strides.

Well, if you think about it, this -keep_strides idea is actually the default already: we keep the strides of the input image… What you mean by ‘keep strides’ here is actually ‘take the strides from the template image’ (note that -template is an option, mrtransform can and will operate without it). And you’re right that one of the uses of mrtransform is to regrid one image onto another, but it’s by no means the only use.

@Romain_Valabregue

I do not get your point there, why this only define the voxel grid ? the voxel order should be part of it. (and for the reslice to be correct it also have to take into account voxel dimension and volume orientation ).

In NIfTI, you’re right that the voxel grid is also defined by the order of the axes. In MRtrix3 however, that is kept logically distinct: you can have images with different strides (and hence different ordering of the axes; or equivalently, rotated 90° in various ways; or equivalently, stored e.g. RAS, LPI, PSL, …) that nonetheless map onto the same exact voxel grid in real space. And MRtrix3 is designed to behave as expected in these situations: if the images map onto each other voxel-for-voxel in the viewer, they will also map onto each other in other computations. This is actually an issue that has been debated quite a bit amongst ourselves, and dates back from the very early days of MRtrix (well before its initial release back in 2008…). Have a look here for some of the background if you’re really interested…

@Romain_Valabregue

So in short I thing to have nifti files with different voxel ordering is not a good idea, (unless you only use the best software ( I mean mrtrix of course !). But I still have old habbit with fsl utilities

Yes, this is an issue that you have raised in the past, and we made quite a few changes to ensure these issues were kept to a minimum. So converting between NIfTI & mif format should preserve the strides, and even computations that do involve changing the stride to volume contiguous storage should preserve the strides if converted to NIfTI. You’ve hit on one of those cases where there’s an ambiguity, since there’s two legitimate sources for that information – and I can see the rationale for expecting the axis ordering to follow the template, especially if coming from a NIfTI perspective. But like I said, I’m not convinced either one is more ‘correct’, without knowledge of the specific use case that the user has in mind. So given all of these considerations, I reckon we’re going to go with the approach suggested… although it’s not black & white, I’m still open to different opinions on that one!

Also, when you say that having NIfTI files with different voxel ordering is not a good idea, I agree that users should keep an eye on this and make sure all their data maintain consistent strides, for the reasons you mention. But this is not something we can easily enforce at the software level: MRtrix used to force NIfTI-standard 0,1,2 ordering in its output, irrespective of the input strides, but as I’m sure you’ll recall from this old conversation, this caused very similar issues to those you’ve raised here…

[Romain_Valabregue]

4 years later, I forget about that point and get the same bug (but with mrgrid). Using mrgrid with the -template option will not preserve the strides. For that reason, I still think that it is not a “strict” reslicing operation.

Reading this thread again I fully understand the internal mrview point of view (and the optimization issue) but my point is that it prevents the basic reslicing tool (mrgrid input.nii -template …) to be used outside mrtrix.

There are plenty alternative to perform the reslicing operation (and I pretty sure they all preserve the strides). But I like very much the piping options of mrtrix command (for instance with combination with mrcalc), it is very useful to spare a few “writting to files” steps

So why not keep the same default behavior and add an option to preserve the stride ?
as @jdtournier suggested :

mrtransform in.mif -template template.mif -strides template.mif output.mif

or just an option like
–preserve-stride

[maxpietsch]

Not sure I understand, but does
mrgrid in.mif regrid -template ref.mif -stride ref.mif out.mif
not work as expected?

Maybe we should add a warning to mrgrid and mrtransform if the order of the spatial strides of template and output do not match and no stride option was provided?

[jdtournier]

I’ve just double-checked just to make: the use of the -stride template.mif option works as expected in mrtransform and mrgrid, so that is indeed already possible.

But I think there is a legitimate discussion to be had around whether the use of the -template option should mean that the output strides follow the template rather than the original input (in the absence of an additional -stride option). Currently the strides match those of the input regardless of the presence of the -template option, and that makes sense in MRtrix-land where the strides and the image transform are considered as separate entities. With NIfTI however, the only way to achieve different strides is purely via 90⁰ rotations in the image transform, and many 3rd-party software packages will make stronger assumptions about the transforms matching (i.e. they will make no allowances for the potential of a 90⁰ rotation via the image transform being spatially equivalent to a regridding without a rotation – if that makes sense). So I can see the appeal of defaulting to the strides of the template if one is provided, especially if there is no overriding rationale to use one over the other.

The one case I’d be concerned about is how to handle situations where the input is volume-contiguous (e.g. a dMRI dataset), but the template image is not (e.g. an anatomical reference). I can see cases where losing the volume-contiguity would be detrimental to performance of subsequent processing steps, and it wouldn’t be obvious to the user. This wouldn’t impact NIfTI at all since it doesn’t support volume-contiguous storage anyway. So my concern would be that by catering for NIfTI, we could compromise performance for data stored using our own MRtrix storage. I’m not sure how to handle that transparently… Open to suggestions though!

[maxpietsch]

One can argue either way. I favour consistent behaviour between mrtranform and mrgrid and defaulting to preserving properties unless explicitly requested. The issue with interoperability with other packages can be addressed if not solved via warnings if the template and output strides do not match (mrgrid and mrtransform: check spatial stride order by maxpietsch · Pull Request #2460 · MRtrix3/mrtrix3 · GitHub). This also increases awareness of these stride issues which unfortunately seems necessary.

[jdtournier]

@maxpietsch

This also increases awareness of these stride issues which unfortunately seems necessary.

Yes, that’s a definite take-home message from this discussion…

Let’s carry on the discussion on that GitHub pull request of yours.

[Romain_Valabregue]

@jdtournier

I’ve just double-checked just to make: the use of the -stride template.mif option works as expected in mrtransform and mrgrid, so that is indeed already possible.

Great I did not realize it was already there ! (and it works)
I’ll add this by defautl in my scripts.
many thanks !

@maxpietsch

One can argue either way.

either way depending on the objective ?
But for the case of nifti images and for interoperability purpose, there is only one way … regrid operation should follow the stride of the template.

but agreed a warning is a good compromise

[Romain_Valabregue]

Hi
sorry to come back here on this interoperability issue, but I just realize, that even with this -stride
option there are still some (minor) changes :

It adds an extra 4th dimension (equal to 1) even my input and my template are just 3D volume

for instance

Dimensions: 136 x 206 x 153 x 1

(Luckily this does not affect mri_segstats for which I need this reslice) but still why this extra singleton dimension ?

[jdtournier]

Hi Romain,

I can’t reproduce this on my system. Can you take me through the exact steps so I can figure out where the issue is? For reference, this is what I tried on my system:

# input image:
$ mrinfo in.nii 
************************************************
Image name:          "in.nii"
************************************************
  Dimensions:        84 x 84 x 48
  Voxel size:        3.04762 x 3.04762 x 3
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9997    -0.01396    -0.01968      -118.9
                          0.01234      0.9967    -0.08022        -109
                          0.02074     0.07996      0.9966      -95.82
  comments:          MCI017^
  mrtrix_version:    3.0.3-98-gf586a656

# template image:
$ mrinfo template.nii 
************************************************
Image name:          "template.nii"
************************************************
  Dimensions:        192 x 254 x 256
  Voxel size:        0.9 x 0.898438 x 0.898438
  Data strides:      [ 3 -1 2 ]
  Format:            NIfTI-1.1
  Data type:         unsigned 16 bit integer (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9992     0.04071   0.0004318      -89.22
                         -0.04071      0.9992     0.00267      -82.66
                       -0.0003228   -0.002686           1      -129.6
  comments:          TOURNIER DONALD (
  mrtrix_version:    3.0.3-98-gf586a656

# apply mrtransform without -stride:
$ mrtransform in.nii -template template.nii out.nii 
mrtransform: [100%] reslicing "in.nii"

# check output:
$ mrinfo out.nii 
************************************************
Image name:          "out.nii"
************************************************
  Dimensions:        192 x 254 x 256
  Voxel size:        0.9 x 0.898438 x 0.898438
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9992     0.04071   0.0004318      -89.22
                         -0.04071      0.9992     0.00267      -82.66
                       -0.0003228   -0.002686           1      -129.6
  comments:          MCI017^
  mrtrix_version:    3.0.3-98-gf586a656

# apply mrtransform with -stride:
$ mrtransform in.nii -template template.nii -stride template.nii out.nii -force
mrtransform: [WARNING] existing output files will be overwritten
mrtransform: [100%] reslicing "in.nii"

# check output:
$ mrinfo out.nii 
************************************************
Image name:          "out.nii"
************************************************
  Dimensions:        192 x 254 x 256
  Voxel size:        0.9 x 0.898438 x 0.898438
  Data strides:      [ 3 -1 2 ]
  Format:            NIfTI-1.1
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9992     0.04071   0.0004318      -89.22
                         -0.04071      0.9992     0.00267      -82.66
                       -0.0003228   -0.002686           1      -129.6
  comments:          MCI017^
  mrtrix_version:    3.0.3-98-gf586a656