Join GitHub today
GitHub is home to over 31 million developers working together to host and review code, manage projects, and build software together.Sign up
image_data.rst - typo in strides for PIR orientation #856
I am trying to understand the definition of data strides in documentation. Quite confused me that strides of PIR data storage are specified as -2,-3,1. According to my understanding it should be 3,-1,-2. I have checked this also by running mrconvert and inspecting the result in freesurfer's mri_info which provides the info about the direction of image axes.
Also, maybe it should be more clearly described (ideally with some example) how the strides specification is translated into data ordering. There is only provided RAS example which is trivial.
My understanding for conversion is as follows:
From data ordering to strides: for each direction (RAS) which defines x,y,z axes of NIfTI coordinate system, find its order and direction in data ordering to be translated. For example PIR ordering: Take first letter from RAS, R (corresponding to x NIfTI axis), it is in 3rd position of PIR and correct direction, therefore first stride is 3. Take second letter, A, it is in 1st position of PIR and reverse direction, therefore second stride is -1. Take S, it is in 2nd position of PIR and reverse direction, therefore third stride is -2.
The understanding of symbolic strides could be also be inferred to imagine how the nested for cycles in hypothetical data reading code are organized. I.e. the index which increments fastest has stride 1 etc.
Is my understanding correct?
Yes, this one is never easy to communicate, maybe I should try to explain that better. I think you're right, the PIR example in the docs is wrong, it should read
So the way I interpret strides, taking your PIR example, is by saying:
giving the expected strides
But really, the only real way to understand this is by thinking in terms of the strides needed to navigate the data properly - i.e. the mapping from voxel coordinate to index into the linear data array. If you imagine a 2×2×2 image, its actual strides in RAS would be
But this is for a 2×2×2 image. If the image had been 5×5×5, its actual strides for a PIR image would have been
As you can see, the actual strides are intimately dependent on the image dimensions. So to make the strides specification independent of the image dimensions, we express them as symbolic strides, where all that matters is their relative absolute magnitudes. In this way, we interpret
The above treatment extends naturally to arbitrary numbers of dimensions, which is probably the reason why we settled on this terminology. There are no additional labels for the 4th, 5th, etc dimensions in NIfTI. Besides, NIfTI can't store data that arbitrarily anyway, we can only modify the first 3 strides, but never have volume-contiguous storage (although I think there's a way to do it in the specs, but I don't think many software packages support that).
Does this help clarify things...? If so, I might have to think about how to update the docs with something like this. Feel free to tell me whether anything I've said actually helped or not, or suggest ways of improving what's in the docs - always appreciated.
added a commit
Dec 8, 2016
thank you for the feedback and clarification. I was not used to think in strides, therefore it requires some time to get used to that. I had to think about your examples with differing dimension across axes, since it was not clear how to interpret to which axis the particular size belongs. If I say that image dimension for PIR image is 3x4x5, it means, that
If I print the image info by mrinfo -dim, or fslinfo (FSL), or mri_info (FreeSurfer) what is the interpretation of values? It seems that for fslinfo and mri_info the interpretation is 2) (i.e. the order of numbers changes when the strides via mrconvert -strides are changed), whereas for mrinfo the interpretation is 1) if the -norealign option was NOT provided. If the -norealign IS provided, then the interpretation is also 2), if I understand it correctly.
However, it is so easy to get confused with this...
Definitely... I'll try to clarify how I think about this stuff, hopefully I'll answer your actual question somewhere in there - just bear with me.
So the main thing to understand first off is that there is an interaction between the transformation matrix and the strides. The way I think about this (and is how these issues are represented within the MRtrix3 codebase itself), is via 3 different coordinate systems (note that this isn't necessarily how everyone conceptualises these things, it's just my way):
So the issue is that there is some redundancy between these two things: if I was to switch the strides so that we now store voxels contiguously along y, then x in voxel space instead of the standard x then y, that is essentially equivalent to using a modified transform where the image x & y axis are exchanged - without explicitly changing the strides as such. In both cases, the voxel values are stored in the same order on file, and the location of any given voxel is the same in scanner space. This is where things get a bit messy, since different packages will deal with this differently (well, to be honest I think it would be more appropriate to say that MRtrix3 deals with this differently from everything else...).
If we focus on NIfTI in particular, things are relatively simple since that format assumes fixed, standard strides: voxel intensities are always stored contiguously along the x-axis, then each row is stored contiguously along the y-axis, etc. (what we would refer to in MRtrix as strides
there's a bit of confusion here with what is meant by the 'NIfTI axis'. Using your example, the size along the voxel space x-axis is indeed 3, but whether that corresponds with scanner space L->R depends on the transform. For a PIR image, that would actually map to A->P, which is a closer match to your interpretation 2... So, I think I can see where the confusion comes from, but before I go into that, I just want to clarify a few things about how MRtrix3 handles things:
However, when you don't use the
So having clarified this, I think the confusion probably arose due to my implicit assumption about what is meant by an image with dimensions 3×4×5. Since in MRtrix3 we always have these dimensions reported after the transform has been modified, that means that the dimension along the x-axis genuinely is that along L->R, etc. So when I talked about an image with dimensions 3×4×5, and then storing that image as RAS or PIR, those dimensions assumed your interpretation 1 - and if you think about it, it has to, anything else would mean the resulting images would no longer be equivalent. However, technically what is stored in the NIfTI header is quite different now, with a very different transform matrix, different dimensions along each axis, and voxel values stored in a different order on disk. Most NIfTI-handling tools would report exactly what was stated in the headers. In contrast, in MRtrix3, NIfTI is just one of a number of possible image formats, and is therefore treated in ways that simplify concurrent processing of these different types of images. In particular, MRtrix3 will modify the information to make it more consistent with any other images that might be processed alongside it. For instance, if you were to create a small RAS image with dimensions 3×4×5, e.g.:
then convert that to PIR:
and inspect with
This is as it should be: the information content is identical between the two images. All we've done is order the voxels differently on file... So any MRtrix3 application that needs to process these two images together will work as expected. For example, subtracting the two images voxel-wise works exactly as expected, returning an image that contains zero values everywhere - even though the voxel values in the two images are ordered differently on file:
However, it you look at the information in the headers without allowing MRtrix3 to modify anything, i.e. using
Clearly, if you were to subtract that from the original RAS without the modifications we normally do, the images would simply not match, and the application would either bail out with a message reporting the mismatch, or it would just crash out. This is why MRtrix3 performs these modifications.
Hopefully that clarifies what's going on here...?
Here is a dump of issue 74 on GitLab referenced above, just for completeness:
Should MRtrix reshuffle the transform on load by default?
J-Donald Tournier @jdtournier commented 2 years ago
The current default for MRtrix is to massage the image transform into as near an axial orientation as possible by switching columns and/or inverting them. This is all done in a self-consistent manner, of course: any column switch will be accompanied by the corresponding changes to the data strides, and any inversion of an image axis will be accompanied by a change in the origin component of the image transform, leaving the data themselves untouched.
I've always thought this was a nice feature since it allows images to be displayed in a much consistent manner in the viewer, regardless of the way they were acquired and/or stored. It also allows images that have been stored using different transforms / strides to be processed together sensibly - most of the time. Issues occur when operating on images, one of which has had its transform modified such that it has been rotated from its original orientation by an amount sufficient to trigger a transform reshuffle (or at least, to trigger a different transform reshuffle). The problem then is that trivial operations such as a straight copy will end up looping over different axes from what might have been intended - for example what was x in the original image might now be considered to be y after the reshuffle - which obviously isn't what the developer would have intended. For now, this issue has been fixed in
So the question is: should that behaviour be removed altogether, and perhaps limited to the viewer alone? It would be trivial to shift that code out of the default path in the
Is this something worth worrying about, or is the current behaviour the safest way to proceed, with the option of disabling it in the few cases where it is not desired...? If the latter, suggestions are also welcome for a cleaner way of preventing the transform reshuffle...
Dave Raffelt @draffelt commented 2 years ago
Just a couple of questions to make sure I understand (I've always wanted to get my head around this):
If I've understood correctly, then I'd probably vote for what you are doing now, just disabling it where it is a problem.
J-Donald Tournier @jdtournier commented 2 years ago
I must admit, as I was writing the original issue, I went from not too sure the reshuffle was necessary to pretty convinced it was absolutely required. Especially since MRtrix allows applications to specify the output strides, so there is no guarantee that an image processed by MRtrix will end up with the same strides as the original. This is in my opinion a good thing, but it does mean we have to take special precautions to ensure it doesn't produce weirdness down the track.
One other option that I hadn't thought of till now, is to only reshuffle the transform on image load, and leave it as-is on image create. This would ensure that any images read will align properly, but prevent situations like the one I've outlined above when creating a new rotated image. In 99% of cases, the input transform is used as-is to create the output image, which would imply that it's already been reshuffled - and any app creating the output transform from scratch will probably create a near-axial one anyway. That would prevent the issue in
Dave Raffelt @draffelt commented 2 years ago
Thanks for explaining it all, makes good sense now. Sounds like a good solution to me. Maybe we could also add an option to
J-Donald Tournier @jdtournier commented 2 years ago
Good idea re
Introducing image space in context of mrtrix brought a confusion to me since it was not clear to me, how do you define this space. What direction the particular x,y,z axes of image space really point at? The data are usually 3D array and can be viewed in whatever orientation, so the definition of image space cannot be related to way the image is displayed. Do I grasp it correctly according to your statement
that the direction of x,y,z axes of voxel/image space in your understanding is defined by their corresponding columns in qform matrix (i.e. first column defines x axis, second column defines y axis, etc)?
In contrast, my previous understanding of axes in voxel and image space was primarily related to way how the voxel indexes in file are arranged and their position is mapped by qform, i.e. x axis defines most contiguous dimension in file (with stride 1), y axis is second etc.
That it is also in conformance with NIfTI standard, where
The concept of "variable strides" introduces one more level of complexity above that. Maybe it would help to provide in the documentation extension of the original nifti1 formula (or point to particular piece of source code) how precisely the index from file data array is mapped to the mrtrix scanner space RAS coordinate.
By NIfTI axis I meant the axis in scanner space, which is defined as x pointing from left to right, y from posterior to anterior, z from inferior to superior (RAS).
On the other hand, I understand that by introducing concept of strides you can assure that orientation of x,y,z axes in coordinate system you use in mrtrix is constant, identical to NIfTI standard
This is in contrast to image space coordinate system (as I understand its definition) where each axis can point to whatever direction.
This approach has big advantages, for example (and here I am replicating your previous appended post):
But, it seems to me, that to grasp this concept it requires to abandon image space and think only in mrtrix/scanner space.
I have one more think to discuss which is only partially related to this topics (and maybe deserves stand-alone chapter): I found that when qform_code is set o 0 in NIfTI dataset, fslview does not display orientation labels. In contrast, mrview produces
And, this generally possibly brings complications to all image handling routines which make use of mrtrix coordinate system (which implies that image transform has to be valid). What is general behaviour in this situation? Does mrtrix allow to process these files?
It seems the more I try to explain this, the more confusion ensues...
OK, the first thing is to cleanly separate the NIfTI standard from MRtrix3's implementation and handling of images in general. Bear in mind that MRtrix started life before NIfTI came along, and supports a number of other image formats besides NIfTI - in fact NIfTI support came relatively late to MRtrix, well after its initial release (bearing in mind that it had been used in-house for quite a few years before then). So a lot of the image handling is designed to be agnostic to the actual format of the images, and to provide a coherent view of the data regardless of the format these data might be stored in.
Back to NIfTI: it is clearly described in the file you linked to (which is actually included wholesale in MRtrix3), and that provides the basis for interpreting these images. The NIfTI standard has no concept of (variable) strides: it states unambiguously (as you mentioned) that:
It also states that the axes along which i,j,k point in scanner space are provided by the qform/sform (although it's not clear which should be used when both are present). So that's pretty much all we need to know to unambiguously map a voxel value in the file to its position in real space.
Everything else after that is how MRtrix3 handles these images - and this applies not just to NIfTI format images, but any other format too, be it DICOM, mif/mih, mgh/mgz, Analyse, etc. But first, I'll try to clarify exactly what space is what, where strides come into it, and how that all relates to NIfTI.
So all this says that when the strides of the image are 1,2,3, then everything is as you expect. This is what
The bit that is causing confusion here is probably what happens next, the additional transform reshuffle, which happens later in a code path common to all image formats. This is nothing to do with NIfTI per se, it's what MRtrix3 will do to all images regardless of their original format. The reason it can happen at all is that MRtrix3 has a very flexible concept of strides: voxel values can be organised in any order, the data are just treated as a N-dimensional array, and there is no particular expectation that values along the first axis should be contiguous. In fact, for performance reasons MRtrix3 will often rearrange the values so that they are contiguous along the volume axis (i.e. axis 3, indexed from 0) - if the output image format supports it (although I think the only format that fully supports arbitrary strides is MRtrix3's own mif/mih format). To navigate this array, all we need to know are the strides. I've already explained how these are handled, and how a symbolic specification (e.g. 2,1,-3) can be expanded to the actual strides given the dimensions of the array, hopefully that was clear enough then.
It's important to realise at this stage that the strides are conceptually entirely independent of the image transform. The latter only relates to the first 3 (spatial) axes anyway. Now I can take a NIfTI image interpreted as per the standard (i.e. strides 1,2,3), with a potentially large rotation in the transform, and that is already a fully valid image. The reason for the further modifications that MRtrix3 makes is as explained previously: to allow sensible processing of the data alongside other images that might have been processed differently, and been regridded from e.g. PIR to RAS, etc. even though the voxel centres lie on the same exact grid. The only way to allow images regridded like this to be processed as expected alongside their original non-regridded version, in a way that is transparent to the user (or indeed the app developer) is to do the modifications I mentioned previously. We have the flexibility in MRtrix3 to modify the strides, so what we can do is take the axes defined by the incoming transform (i.e. qform/sform for a NIfTI image), and permute/flip them around so that they point close to their corresponding anatomical axes (i.e. scanner space), and adjust the strides accordingly. Note that this now means that image & voxel space no longer correspond with the original qform/sform axes in the NIfTI header, but with the transform matrix created by this realignment process. The voxel space (in MRtrix3 speak - not NIfTI) should now have its axes pointing close to the anatomical L->R, P->A, I->S axes, and you now need to take account of the strides to figure out the offset into the data array of a voxel at index i,j,k.
So a voxel at i,j,k in a PIR NIfTI image that has been imported via this process is actually at -j,-k,i in the original image (note that this means you need a non-zero base offset here, otherwise you will end up with negative offsets). What this translate to as far as MRtrix3 is concerned, is into a different formula for computing the offset. Assuming the image dimensions were 3×4×5 (as per the example above), where these dimensions relate to the scanner frame (i.e. the dimension along the axis closest to L->R is 3, etc), if the image was stored RAS, the actual strides would be 1,3,12, and the offset for a voxel at x,y,z would be x+3y+12z. If the image was PIR, its dimensions as stored in the NIfTI header would actually be 4×5×3 (as per one of my previous posts), its unmodified strides would be 1,4,20, and the offset for a voxel at i,j,k would be i+4j+20k - but note that in this case, i,j,k would not correspond to the x,y,z used in the RAS image. After the transform reshuffle, the new strides would be 20,-1,-4, and the offset for a voxel at x,y,z would be 19+20x-y-4z (note the additional non-zero base offset here), and in this case, the x,y,z indices used here are a direct match to the RAS case. You can verify all this by running e.g.
If you don't do these kinds of modifications, you can end up with issues such as derived images (that may have been regridded) not matching with their original non-regridded image in the display (many reports of that on the various mailing lists and forums in relation to
That's a sensible suggestion, but would be very problematic to do in practice, I think. First off, some image formats simply don't provide orientation information - Analyze is a case in point. More to the point, as you allude to, we need to allow processing of these images in general, and in many cases we just have to put something in there, otherwise a lot of the MRtrix3 commands will fall over. This is mostly a consequence of the explicit design decision I made over a decade ago that processing should happen in the scanner frame where possible. For example, even something as simple as importing the bvecs/bvals can't be done in MRtrix3 if the image doesn't have a valid transform, since MRtrix3 will internally convert the DW gradient table to scanner space. As least by resetting the transform to unity, the bvecs/bvals can be imported for that image. Same for computing tensors and FODs: the coefficients stored in both cases are provided with respect to the scanner frame, not the image axes. Tractography also won't work, since the input FODs or tensors need to be provided with respect to scanner space, and the streamlines are to be stored with respect to scanner space also. Also, any image registration needs a valid transform. So I think the simplest approach is simply to warn the user about this, but nonetheless insert some vaguely sensible default in there to allow subsequent commands to proceed.
One possible approach is to detect these cases in MRView, and not display the orientation labels in this case. But I'm not sure that's such a good idea. The point is, in any other command, the image will be interpreted as if it were pure axial (the user will get the warning, but that's it). So it makes sense to have the viewer also show that interpretation, if only so users can see that this is what would happen if they were to process their images with any other commands. Another issue is that the menu entries to change orientation are labelled axial / sagittal / coronal, so the very fact that we display the images and you can chose one of those orientations implies that we've made these very assumptions about the transform. I can see why
You also have to take into account that MRView indeed needs this information (since stuff is displayed in scanner space); so hiding the labels doesn't mean they aren't used and relied upon even on that very moment. Apart from the (maybe more formal/abstract) same arguments @jdtournier provides about this in the last paragraph above, maybe consider this simple example: assume you open such an image (where the information is missing) as an overlay to an image where it isn't missing. This alone already shows you need to assume something somehow; how else will you have a definition of how to overlay the problematic image onto the other one? Also, from a GUI point of view, the labels will already be present anyway (because the main image didn't have any issues). We don't show the labels for separate images (or overlays): we just show them for the space, and the images are shown in this space. In a way, even the main image already is "an overlay": the real/scanner/... space (labels and all) "is already there", and anything that is loaded is displayed in it.
Dear @jdtournier ,
Thank you really for the patience, I think that finally you clarified all the aspects of the problem very well :-)
As for the issue how to handle the missing image orientation in