Making a 3D sequence from 4D volume

Hi,

I have loaded in a series of 2D sequences (3 dimensions) and stitched them together into a 4D image. Now I would like to save the image as a 3D sequence. I think that to do this I need an image that has 3 dimensions and “NumberOfComponentsPerPixel” is the size of the 4th dimension (“domain, domain, domain, list”), but I haven’t been able to figure out how to do that. The reason I think this, is because if I have a 4D numpy array and use GetImageFromArray, I get an image with 3 dimensions and “NumberOfComponentsPerPixel” is greater than 1, and when I save it as a seq.nrrd I am able to load it as a sequence volume in slicer.

Is that the correct way of doing it, or is there an easier way?

Thanks,

Juan

Hello @jrojasUNC,

I was trying to follow the question but can’t figure exactly what is going on. Could you post some pseudocode with what you expect the image dimensions and pixel types to be (vector/scalar). Note that the GetImageFromArray has a isVector flag to tell SimpleITK whether the numpy array should be interpreted as regular dimensions or as vector pixels.

I think your code works for you and you are asking if this is the correct approach? I wasn’t sure about this, hence the text above.

In any case, if your code works then that is the approach you should be using :slight_smile:

I guess the question is how to make ITK write a .nrrd file with “kinds: list domain domain domain” or “kinds: domain domain domain list”.

It seems that the only way to do write a 4D file is to create a 3D volume with multiple scalar components (one scalar component for each time point). This is stored in memory so that all time points are stored for the first voxel, then for the second voxel, etc. The problem with this memory organization, it is computationally quite expensive to add, remove, or extract additional 3D volumes.

Can ITK be configured to store a 4D volume as a vector of 3D volumes in memory (so that in memory all voxels of t=0 are stored, followed by all voxels of t=1, then all voxels of t=2…)?

Can ITK store 4D color volumes (in total 5D: 3 spatial + color + time) in memory and read/write files?

I have so far no problems with saving 4D images in MetaImage format (scalar or RGB). Just create 4D image and appropriate writer.

ObjectType = Image
NDims = 4
BinaryData = True
BinaryDataByteOrderMSB = False
CompressedData = False
TransformMatrix = -1 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 1
Offset = 99.5 -301.5 -149 0
CenterOfRotation = 0 0 0 0
ElementSpacing = 0.38867200000000002 0.38867200000000002 10 1
DimSize = 512 512 2 3
AnatomicalOrientation = ???? 
ElementNumberOfChannels = 3
ElementType = MET_UCHAR
ElementDataFile = LOCAL

MetaImage RGB8 4D

MetaImage RGB16 4D

To generate them from set of 3D images wrote small function just using iterators.

Edit: Saved to Nrrd

Nrrd RGB8 4D

Looks like it is OK, isn’t it?

NRRD0004
# Complete NRRD file format specification at:
# http://teem.sourceforge.net/nrrd/format.html
type: unsigned char
dimension: 5
space dimension: 4
sizes: 3 512 512 2 3
space directions: none (-0.38867200000000002,0,0,0) (0,0.38867200000000002,0,0) (0,0,-10,0) (0,0,0,1)
kinds: RGB-color domain domain domain domain
encoding: raw
space origin: (99.5,-301.5,-149,0)

Edit 2: Source is DICOM Supplemental Palette series 3 x 3D multi-frame files (D Clunie), no chance to open with colors with ITK DICOM support, may be Slicer can open in multi-volume mode as scalar, not sure

DICOM

Edit 3: sorry, formatting was bad and typos, corrected

Thank you very much for looking into this. This is almost correct, but the issue is that there are 4 spatial dimensions. Can you create image with “kinds: RGB-color domain domain domain list”? How would you declare this image type in ITK? Can you load such an image? When you load it, is the image type preserved or it is read as a multi-component scalar image?

We used metaimage for many years but in recent years we had to switch to nrrd because metaimage does not have a stabdard way to specify axis kinds, so you cannot reliably store data of 3 or more dimensions. E.g., if you have a 3x10x20x3 image then you don’t know if it is a 3D RGB volume, or a time sequence of a 3D scalar volume, and which axis correspond to which dimension.

2 Likes

To save i only passed image to writer, for example something like this, nothing special.

typedef itk::Image<unsigned char, 4> T4d;
typedef itk::ImageFileWriter<T4d> WriterType;
typename WriterType::Pointer writer = WriterType::New();
writer->SetInput(image);
writer->SetFileName(f);
try
{
  writer->Update();
}
catch (itk::ExceptionObject & ex)
{
  return QString(ex.GetDescription());
}

Works with .mha/.mhd or .nrrd extensions.
I didn’t look closer into, didn’t ever know there are different possibilities (list or domain for 4th dimension). To load i decompose 4D into set of 3D, but it is application specific stuff, IMHO.

Edit: corrected typos, sorry

IMHO, number of dimensions returned by itk::ImageIOBase is correct, but may be i only had luck and didn’t run into problems. All this is rather rarely used and difficult to validate.

#include "itkImageFileReader.h" // required!!
#include "itkImageIOBase.h"
#include "itkImageIOFactory.h"
#include <iostream>

int main(int argc, char ** argv)
{
  if (argc < 2)
  {
    std::cout << "File name is required" << std::endl;
    return 0;
  }
  unsigned int dim = 0;
  itk::ImageIOBase::IOPixelType pixel_type =
    itk::ImageIOBase::UNKNOWNPIXELTYPE;
  itk::ImageIOBase::IOComponentType component_type =
    itk::ImageIOBase::UNKNOWNCOMPONENTTYPE;
  unsigned int component_number = 0;
  itk::ImageIOBase::Pointer imageIO;
  try
  {
    imageIO = itk::ImageIOFactory::CreateImageIO(
      argv[1],
      itk::ImageIOFactory::ReadMode);
  }
  catch (itk::ExceptionObject & ex)
  {
    std::cout << ex.GetDescription() << std::endl;
    return 0;
  }
  if (imageIO.IsNull())
  {
    std::cout << "Can not create ImageIO for "
      << argv[1] << std::endl;
    return 0;
  }
  //
  imageIO->SetFileName(std::string(argv[1]));
  try
  {
    imageIO->ReadImageInformation();
  }
  catch (itk::ExceptionObject & ex)
  {
    std::cout << ex.GetDescription() << std::endl;
    return 0;
  }
  catch (...)
  {
    std::cout << "Unknown exception" << std::endl;
    return 0;
  }
  component_type   = imageIO->GetComponentType();
  dim              = imageIO->GetNumberOfDimensions();
  component_number = imageIO->GetNumberOfComponents();
  pixel_type       = imageIO->GetPixelType();
  //
  std::cout << std::endl << imageIO->GetNameOfClass()
    << " " << std::endl;
  std::cout << "Image dimensions: " << dim << std::endl;
  std::cout << "Pixel type: ";
  //
  switch(pixel_type)
  {
  case itk::ImageIOBase::SCALAR:
    std::cout << "scalar";
    break;
  case itk::ImageIOBase::RGB:
    std::cout << "RGB";
    break;
  case itk::ImageIOBase::RGBA:
    std::cout << "RGBA";
    break;
  case itk::ImageIOBase::OFFSET:
    std::cout << "offset";
    break;
  case itk::ImageIOBase::VECTOR:
    std::cout << "vector";
    break;
  case itk::ImageIOBase::POINT:
    std::cout << "point";
    break;
  case itk::ImageIOBase::COVARIANTVECTOR:
    std::cout << "covariant vector";
    break;
  case itk::ImageIOBase::SYMMETRICSECONDRANKTENSOR:
    std::cout << "symmetric second rank tensor";
    break;
  case itk::ImageIOBase::DIFFUSIONTENSOR3D:
    std::cout << "diffusion tensor 3d";
    break;
  case itk::ImageIOBase::COMPLEX:
    std::cout << "complex";
    break;
  case itk::ImageIOBase::FIXEDARRAY:
    std::cout << "fixed array";
    break;
  case itk::ImageIOBase::MATRIX:
    std::cout << "matrix";
    break;
  default:
    std::cout << "unknown pixel type";
    break;
  }
  std::cout << std::endl;
  //
  std::cout << "Component type: ";
  switch(component_type)
  {
  case itk::ImageIOBase::UCHAR:
    std::cout << "unsigned char";
    break;
  case itk::ImageIOBase::CHAR:
    std::cout << "char";
    break;
  case itk::ImageIOBase::USHORT:
    std::cout << "unsigned short";
    break;
  case itk::ImageIOBase::SHORT:
    std::cout << "short";
    break;
  case itk::ImageIOBase::UINT:
    std::cout << "unsigned int";
    break;
  case itk::ImageIOBase::INT:
    std::cout << "int";
    break;
  case itk::ImageIOBase::FLOAT:
    std::cout << "float";
    break;
  case itk::ImageIOBase::DOUBLE:
    std::cout << "double";
    break;
  case itk::ImageIOBase::ULONG:
    std::cout << "unsigned long";
    break;
  case itk::ImageIOBase::LONG:
    std::cout << "long";
    break;
#if ((ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR > 12) || ITK_VERSION_MAJOR >= 5)
  case itk::ImageIOBase::ULONGLONG:
    std::cout << "unsigned long long";
    break;
  case itk::ImageIOBase::LONGLONG:
    std::cout << "long long";
    break;
#endif
  default:
    std::cout << "unknown component";
    break;
  }
  std::cout << " (" << component_number << ")" << std::endl;
  return 0;
}
1 Like

@lassoan ITK doesn’t have special treatment of time dimension. That needs to be enforced by the application (e.g. Slicer). Current NRRD writer in ITK can only have non-domain kind for the first dimension.

What you seem to want is the ability to write using ImageRGB3DplusTime = std::vector<itk::Image<itk::RGBPixel<unsigned char>, 3>>; into a single .nrrd file, possibly with other containers besides vector. This would require significant changes to ITK’s IO plumbing.

2 Likes

Yes, I was interested in exactly this. I was hoping that with so much improvements in ITK in the last few years, there is stronger support for 3D+t data, too, but it seems that this still mostly has to be managed at application level. We do this already in 3D Slicer, I was just wondering if we could leverage some more ITK infrastructure for this.

Thank you @dzenanz and @mihail.isakov for the answers.

2 Likes

Hi all,

Thank you for all the comments.

If I have a numpy array that is 100x512x128x200 where the last dimension is the time dimension, I use SimpleITK.GetImageFromArray(), which produces an image that has 3 space dimensions (GetSize() = (128, 512, 100)) and 200 components per pixel (GetNumberOfComponentsPerPixel() = 200).If I use SimpleITK.ImageFileWriter() to write it as a “seq.nrrd”, and I can load it into slicer as a 3D + t sequence. @lassoan, it seems that it is doing what you were suggesting for making a 4D image- 3D volume with multiple scalar components.

However, if I have an SimpleITK image where GetSize() = (128, 512, 200, 100) and I save that as a seq.nrrd, this image doesn’t load as a 3D sequence.

GetImageFromArray() is able to format the data to have 3 space dimensions and one time dimension (3D volume with multiple scalar components). Can you “tap into” that functionality and format a 4D (4 space dimensions) simpleITK image into a 3D volume with multiple components so that it can be saved as a seq.nrrd and loaded as a sequence in slicer?

Thanks,

Juan

What command do you use to write the file?

Could you copy-paste here the header of the nrrd file that Slicer cannot load as a volume sequence (and maybe also upload the entire file somewhere and copy the link here)?

It is not clear to me what exactly you are wanting from SimpleITK and ITK in producing a NRRD file to be compatible with 3D Slicer time sequence. I am happy to help however I can.

I have been using SimpleITK with the Bioformats microscopy ImageIO plugin, which commonly presents image files a 5D volumes of XYZTC. The extra dimensions are time then a dimension of number of channel. I am also saving images as XYZC, where the number of channels is approaching 100. For this application it is important that the channels are stored as a dimension and not the “fast” multi-component memory arrangement so that, a single channel can be easily read, displayed and operated on. When I load a 4D or 5D image in ImageJ it prompts me how to interpret the dimensions XYZT, XYZC, etc…

SimpleITK has complete support for 3D images. It’s support for 4D images is the Image class itself, IO, and the grid operations to “glue” and extract volumes together. I have been pondering extending this grid support up to 5D Image. This is independent of support for multi-component or “Vector” images.

Bellow are a couple recent patches in the SimpleITK master branch which might be related:


Please let me know if there is some additional support needed in SimpleITK. There are many options available for manipulation of these many dimensional images. Also I am interested in getting these multi-dimensional images loaded into Slicer as well.

1 Like

How do you represent XYZTC files in SimpleITK in memory now? How do you plan to improve it?

Currently, in Slicer in memory we store images in vtkImageData (CXYZ) and Sequences extension manages a vector of any data nodes (resulting in CXYZ+T volumes, where each CXYZ array is contiguous memory block, but memory for each time point is managed separately).

Me too! 3D Slicer supports 4D and 5D volumes and time sequence of any other data types (transforms, polygon and volumetric meshes, annotations, camera position, display settings, etc.) and many operations on them (sequence registration, segmentation, cropping, visualization - including volume rendering, etc.). See a recent example here:

Loading/saving of 3D+t nrrd file works well (we use it all the time), so I would expect to work out of the box with (Simple)ITK-written files. We should test other 4D and 5D configurations and maybe also update the direct memory transfer utility classes in Slicer.

@lassoan for writing of all images I use the following code:

    writer = sitk.ImageFileWriter()
    writer.SetFileName(output_filename)
    writer.Execute(final_image) 

The header looks like:

NRRD0004
# Complete NRRD file format specification at:

type: unsigned char
dimension: 4
space dimension: 4
sizes: 128 512 194 10
space directions: (1,0,0,0) (0,1,0,0) (0,0,1,0) (0,0,0,1)
kinds: domain domain domain domain
encoding: raw
space origin: (0,0,0,0)

The file can be found here https://www.dropbox.com/s/10ym46uim95frwq/RawCineVolume.seq.nrrd?dl=0

since kinds: domain domain domain domain, I can’t load it as a volume sequence.

When I use GetImageFromArray() to convert a 4D numpy array into an SimpleITK image and write as a seq.nrrd the header is:

NRRD0004
# Complete NRRD file format specification at:

type: unsigned char
dimension: 4
space: left-posterior-superior
sizes: 194 128 512 10
space directions: none (1,0,0) (0,1,0) (0,0,1)
kinds: vector domain domain domain
encoding: raw
space origin: (0,0,0)

the file is here https://www.dropbox.com/s/fe6fsvizqlyp3mm/RawCineVolume_good.seq.nrrd?dl=0, and I can load it as a volume sequence into Slicer.

@blowekamp Thank you for sharing the current patches. In GetImageFromArray() “4D arrays are automatically considered 3D vector images”, so the function is doing what I need. Is there a way to use the functionality of GetImageFromArray() to change a 4D simpleITK image into a 3D vector (ie kinds: vector domain domain domain)?

May be 128 512 10 194 were better?
Can open that file, but looks somewhat strange

@mihail.isakov Can you load it into slicer as a “volume sequence”? that is what I’m trying to accomplish. I can load it as a volume, but then it is limited to 3 dimensions.

Sorry, not in Slicer

I can make Slicer’s Sequences extension to load any kind of 4D or 5D sequence but we would need to have an agreement about what configurations are valid and needs to be supported.

For example, I would be reluctant to allow loading anything from a sequence of kinds: domain domain domain domain with space dimension: 4 as it does not make any sense. Space is 3 dimensional (or may be restricted to 2), but 4 … We also would not know which axis we should use for time.

As you confirmed, kinds: domain domain domain list and kinds: list domain domain domain already work (the specification is clear).

space dimension: 4 as it does not make any sense

why not, 4th can be time (3D+t), or 3rd can be time (2D+t)

kinds: domain domain domain list

it is Nrrd terminology, MetaImage doesn’t have such wording and works well for me.
But it is, of course, your decision how to do this in Slicer.