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?
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
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?
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
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.
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.
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;
}
@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.
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.
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?
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.
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.
@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)?
@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.
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).
