Mrtransform does not output same data strides as -template volume

preprocessing

#1

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


#2

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…


#3

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 :fire: :rofl:


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)

#4

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

But I have reservations about this:

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?


#5

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.

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 !

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 :blush:


#6

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.

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

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.

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…

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… :stuck_out_tongue_winking_eye: