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?