Reading and writing of 5D images - a dimension disappears

Algorithms
,
1

Hi all,

We’re trying to enable 3D Slicer to read and write vector image sequences (displacement field or velocity field: 3D spatial + component + time = 5 dimensions) into a single file. The format of choice is NRRD, and we thought the most correct way would be to use ITK’s NrrdImageIO (the alternative would be to improve vtkTeemNRRDWriter/Reader, but ITK is a much more core dependency than teem).

I started testing the IO with a valid-looking 5D nrrd file that I created using pynrrd. The header of this file contains this:

  sizes: 2 20 30 40 4
  kinds: time domain domain domain RGBA-color

But when I read it with ITK, and write it back out, the header is changed to:

  sizes: 2 20 30 40 1
  kinds: time domain domain domain domain

I know that the first issue is at reading, because looking at the offsets, the last dimension is squashed to a single component, see
image

I tried different combinations of “kinds”, but the result was always similar (although the position of the vector component was pushed to the last place:

Original
  sizes: 2 4 40 30 20
  kinds: time RGBA-color domain domain domain
Write ITK
  sizes: 2 40 30 20 1
  kinds: time domain domain domain domain

Same thing with

  sizes: 4 2 20 30 40
  kinds: RGBA-color time domain domain domain

The code I used to read and write the 5D image was this: ITKReadWrite5DImage.cxx · GitHub

My question:

  • Is such 5D NRRD IO is supposed to be supported in the first place?
  • Is there is something I miss when reading the file? If not, is something broken in ITK?

Any other tips or pointers are appreciated as well. Thank you very much!

2

To get RGB interpreted correctly, using PixelType = unsigned short; should be changed to using PixelType = itk::RGBAPixel<unsigned short>;, along with #include "itkRGBAPixel.h".

Note: ITK uses a convention that channels (e.g. RGB/A) are the fastest varying component in the pixel array. Image formats which do not follow this, have their pixels rearranged accordingly during reading/writing.

2 Likes
3

Thank you very much for the quick answer! I’ll try this ASAP.

1 Like
4

It also good to look at the information the itk::ImageIO has read about the image. After calling ImageIO::ReadImageInformation, GetPixelType, GetNumberOfComponents, GetNumberOfDimensions and GetDimensions can be called to know what itk::Image template parameters would best match this file. When reading a file into itk::Image the template parameters are already specified and ITK does its best job to convert to the requested type.

Also SimpleITK can read images up to 5D, and auto magically detects if the input is a multi-component or not. It reads all multi-component images into a VectorImage and does not distinguish between RGB, covariant vector etc… So it could be useful for reading these images too.

2 Likes
5

Using the RGBA pixel type solved the reading issue, thanks a lot for the suggestion!

Just for the record (and for the sake of others reading this) I needed to reduce the number of dimensions to 4 because now the RGBA components are stored within the pixel.

Based on the offset table the dimensions correspond to the file that was read in (sizes: 4 40 30 20 2, kinds: RGBA-color domain domain domain time)

image

Now, when I write this image back to NRRD, the rest of the dimensions are reversed:

dimension: 5
space dimension: 4
sizes: 4 2 20 30 40
space directions: none (1,0,0,0) (0,1,0,0) (0,0,1,0) (0,0,0,1)
kinds: RGBA-color time domain domain domain

As I understand the first “size” is the fastest changing one, so now the NRRD contains the time component as the second fastest changing, then k, j, and i. To me it would make more sense to store it as it was: i, j, k, time.

There is no extra specification in the writer code, it directly writes out the read image (that is now correct, see above), see the gist (for convenience, here). What could be the reason for this reversing of the space dimensions? Is there a simple way to keep the order?

Thanks again for the great help!

6

I guess that C++ NRRD IO likes to normalize the order of dimensions to something it likes? I did not look at the above linked code, I just glanced at ITK’s wrapper for it.

7

Keep in mind the difference in indexing with ITK and numpy/python as you said that the original image was written with pynrrd. And ITK image is indexed as image(x,y,z)[c] or image(i,j,k)[c] while a numpy array is index like arr[k,j,i,c]. Here c is fastest followed by i etc… I am not sure what additional meta data is in the header that may indicate a specific ordering of axis.

8

Here is one complete header:

NRRD0004
# Complete NRRD file format specification at:
# http://teem.sourceforge.net/nrrd/format.html
type: unsigned char
dimension: 5
space: left-posterior-superior-time
sizes: 2 5 6 7 4
kinds: time domain domain domain RGBA-color
space directions: (1,0,0,0) (0,1,0,0) (0,0,1,0) (0,0,0,1) none
encoding: raw
space origin: (0,0,0,0)
9

Thank you very much for the answers!

Yes this is considered when writing out.

Do you know about some flexibility in the implementation that help influence the axis orders? Even if in memory the allocation is fine, it would be nice if the file we write out could contain the different time points in “blocks”. For example, when I read in the nrrd file that has the time as second axis to Slicer, I get something like an ice core (a narrow section of the images across time).
So this would mean kinds: RGBA-color domain domain domain time
Is there a way to achieve this without the need to recreate the buffer in memory? Thank you!

10

Would you be asking for essentially an inplace axis permutation filter aka in-place matrix transpose?

p.s. An alternative approach would be to create a streaming pipeline to read input slices then permute axis to assemble to the the desired axis order. Let me know if you want more details on how to implemented this efficiently.

11

Thanks Brad! Doing the reading pipeline is not an issue. We (with @lassoan) are just wondering if we could somehow have the writer write the axes with time being the slowest changing component. If it is not possible then we’ll do it on read time. Or if it is not readily usable but would be easy to add in the ITK writer, we can consider that too.

12

Looks at the NRRD code:

It looks like if you set the MetaDataDictionary of the NrrdImageIO you can control the details of the NRRD header. For example if you set “kinds[0]” to “time” that should be in the output header. Does that help?

13

Hi again,

I’ve made some more investigations about this recently.

If I create a dictionary and set it to the image that is the input to the writer (as well as call UseInputMetaDataDictionaryOn), then it is considered, but unfortunately not the kinds information.

After debugging into the IO, it seems that the first kind element is set from the pixel type, and the rest is always domain, see

I’d like to ask what is the vision of the ITK developers for writing multi-dimensional data such as what I propose (i.e. a sequence of 3D+component images, with the last kinds type time or list)? Is it in the plans to support such file format? Or if not a file format, then any other way to support this use case, such as streaming?

Thank you!

1 Like
14

NRRD has not been updated in ITK in a loooong time. And NRRD image IO is a part of teem so extracting it is non-trivial.

MetaImage is another mature format. It is being continually updated. If it is close to doing what you want, you might propose the changes you want via a PR.

Finally, there is work-in-progress OME-NGFF Zarr support. As neither the format nor its implementation are finalized, it is least hard to get this to support your use case, but it might require the most waiting and/or implementation effort.

@glk, @Stephen_Aylward and @matt.mccormick might want to comment, too.

15

ITK’s NrrdIO does not use teem anymore. Parts of teem related to NRRD were extracted into ITK.

15 years ago MetaImage was our main file format but had to switch to NRRD because MetaImage does not have a standard way to describe image axes (name, meaning, etc. of each axis) and there are controversies around how anatomical orientation information should be interpreted.

We could revisit this again, but considerable development work would be needed to add the missing axis description and axis reorientation/reordering features to MetaImageIO. It is not enough to just put it into the file format documentation that application developers must implement all these features: then the format will not be widely accepted (see why we dropped MetaImage) and inevitable inconsistency between documentation and implementations in various applications will make interpretation of the file format ambiguous (see Nifti’s struggle with image orientation in sform/qform).

Escaping forward to OME-NGFF is promising, we could then leave both nrrd and metaimage behind. Can ITK store 5D data (time sequence of color 3D images; time sequence of displacement field transforms) in OME-NGFF?

How 5D data is supposed to be stored in ITK? Is there a container that can store time sequences or ITK just stores one time point and application developers need to manage multiple time points? How does this work for IO classes (especially writing, where the writer class must be aware of the time sequence)?

Also SimpleITK can read images up to 5D, and auto magically detects if the input is a multi-component or not.

@blowekamp How SimpleITK represents 5D data (time sequence of 3D vector field or time sequence of 3D RGBA image) in memory?

1 Like
16

This enhancement request describes it:

As of right now, it is un-implemented.

1 Like
17

Hello @lassoan,

In memory, SimpleITK just represents a 5D image as an itk::ImageBase<5>. So it is up to the developer to determine the axis order. Most SimpleITK filters support 2D and 3D image processing. However, many image “grid” filters work with 4 and 5 dimensions images, to support working on a sub-volume. There have been recent optimizations to the index operator to support assignment of sub images and masked assignment to improve performance when working with a sub-dimension. Additional an sitk::Image::ToVector and sitk::Image::ToScale to support switching between a channel as a dimension and a channel as a VectorImage.

I have been working with OME-NGFF ZARR with large microscopy images. The images are stored as “TCZYX” to be compatible with the “bioformats” tool suit. Mostly dask is being used to read the ZARR in parallel at the require resolutions for the algorithms run in SimpleITK. Also, having the data stored in the correct chunk size and order is very important for algorithm and visualization purposes. If registration is being done on select channels ( 1 out of 8 ), it would be good if the image is not stored with the channels as the fast axis. However, if you are rendering all the channel then it may be desired for the channels to be in the fast axis, or at least all the channels being in the same chunk.

The last point goes to how the image should be stored in memory or on disk is really up to how the image will be used and the access patterns needed.

1 Like
18

Thanks @dzenanz for tagging me so that I knew about this conversation, albeit belatedly.

Yes, the NRRD format and the Nrrd library for its I/O originates in Teem, but as @lassoan pointed out, a separate NrrdIO library was created (by me) to have standalone self-contained C library to support NRRD file I/O, so that its use in ITK would not have any dependency on Teem. And yes, the NrrdIO part of ITK hasn’t been updated in a long time. I’d guess that that is due to: the stability of the format, the maturity of the code, the fact that I never used ITK and was never paid for ITK development, and that my work time doesn’t have as much for software development as I’d like.

The NRRD file format doesn’t know or care about any programming language or any software for image representation or array manipulation, and their associated conventions for syntactic axis ordering. What it really cares about about is whether an axis is faster vs slower, and how you move in world space as you increment the index along an axis. That is expressive enough to losslessly represent everything about DICOM, and lots of other array storage schemes. For the purposes of the format, you do have to pick an axis ordering (NRRD chooses fast-to-slow, mirroring the language of people who say “640 by 480”, if that still rings a bell for anyone), but the ability to label each axis, and document how the axis moves through space, is set up to hopefully facilitate having the on-disk data storage be exactly the same as whatever in-memory data ordering you’ve chosen, according to the locality properties of your algorithm.

Now, in ~2005, during my postdoc, I did help write some bridge code to go between ITK and NRRD. I wrote it at the direction of others who wanted to save out oriented images (like from DICOM) in a semantically lossless way, but I sensed then that the NRRD format, and the corresponding Nrrd struct in the NrrdIO library, was more flexible about axis ordering than ITK itself was. So - @dzenanz there may indeed be some axis ordering normalization in that bridge code; I forget now. And if it does axis ordering normalization, it may do so with assmptions that made sense in 2005 but not in 2024.

If it would be helpful to you all, and if you all can help me, I can work on reacquanting myself with the ITK-NrrdIO bridge code. I would like to do this (though it may proceed at a more leisurely rate than you’re used to) because my coding work on Teem has meanwhile continued, and I am working towards another release soon, which will have an updated NrrdIO library (not that anything has substantially changed in that code). Reconciling any changes that have accumulated with ITK’s NrrdIO with Teem’s Nrrd support is certainly for me (not you) to figure out, but I’d appreciate help/guidance on using whatever is the current mechanism for contributing and testing changes to ITK I/O code.

5 Likes
19

Contributing guide is a good start. If you have questions, ask.

20

Okay, thanks.

Can you remind me where (in the ITK source tree) is the ITK/NrrdIO bridge code?
And the NrrdIO library itself (ITK’s snapshot of it)?