I have a 3d dicom nii.gz file with shape [900,900,500]. And I read it by using itk::ImageReader. Then I want to split it into 6480 with shapes [50,50,25] by moving ImageBufferPointer. But I don’t know if the order in memory is row-major(C order) or column-major(F order). And I can’t find it in the API documentation or the User Guide.
Besides, that nii.gz file is generated by itk::ImageSeriesFileReader. Does the shape means [w,h,d], or [h,w,d]?
Any suggestions are greatly appreciated.
1 Like
You should take a look at Data Representation chapter, especially section 4.1.7 Importing Image Data from a Buffer.
In ITK, shapes mean IJK, which usually (identity orientation) means WHD. And ordering is C-style, not Fortran-style.
1 Like
I did look at this chapter before, but didn’t find the relevant information. Interestingly, I used itkImageSeriesReader to read a series of DCM files and produce a nii.gz file. Since I didn’t find itk’s API about C-order or F-order, I read it with python’s nibabel library. it shows F-order.
Here is the python code and the output:
import nibabel as nib
data = nib.load("./my_data.nii.gz").get_fdata()
print(data.flags)
C_CONTIGUOUS : False
F_CONTIGUOUS : True
OWNDATA : True
WRITEABLE : True
ALIGNED : True
WRITEBACKIFCOPY : False
UPDATEIFCOPY : False
Maybe this nii.gz file has been modified. I create a nii.gz file with random pixel.Here is my c++ code:
/*Some unimportant details are omitted*/
using PixelType = signed short;
const unsigned short Dimension = 3;
using ImageType = itk::Image<PixelType, Dimension>;
using ImageReaderType = itk::ImageFileReader<ImageType>;
using ImageWriterType = itk::ImageFileWriter<ImageType>;
ImageType::SizeType imageSize = imageIn->GetLargestPossibleRegion().GetSize();
ImageType::IndexType imageIndex = imageIn->GetLargestPossibleRegion().GetIndex();
ImageType::Pointer imageOut = ImageType::New();
ImageType::RegionType imageOutRegion(imageIndex, imageSize);
imageOut->SetRegions(imageOutRegion);
imageOut->Allocate();
ImageWriterType::Pointer imageWriter = ImageWriterType::New();
imageWriter->SetFileName(fileName.substr(0, 16) + s.str());
imageWriter->SetInput(imageOut);
try
{
imageWriter->Update();
}
catch (const itk::ExceptionObject &e)
{
std::cerr << e.what() << '\n';
return EXIT_FAILURE;
}
Then I read the file with python, it also shows the F-order.
More importantly, if I default the ImageType is F-order, then the calculated strides are correct and can be correctly divided into 6480 blocks; if the default ImageType is C-order, it is incorrect
Hi,
nibabel reads in nifti into F-order arrays.
ITK reads in nifti (and other images) into C-order. For comparison with nibabel:
import itk
import numpy as np
image = itk.imread('image.nii.gz')
array = np.asarray(image)
print(array.shape)
print(array.flags)
ITK’s API indexes in IJK order, though. NumPy indexes in KJI order.
1 Like
Hi, thank you for your reply., and I tried it again:
I generate the nii.gz file with c++ itk library using itkImageSeriesReader, the code is totally the same as the Data Representation chapter, section 1.12.3 Reading a 2D DICOM Series and Writing a Volume.
Here is my code:
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include <vector>
#include <string>
int main(int argc, char **argv)
{
if(argc < 2)
{
std::cerr << "Usage: \n"
<< argv[0] << " "
<< "DicomFilesDir" << " "
<< "StoreDciomFileName.dcm" << " "
<< "{seriesIdentifier}" << std::endl;
}
using PixelType = signed short;
const unsigned int Dimension = 3;
using ImageType = itk::Image<PixelType, Dimension>;
using ReaderType = itk::ImageSeriesReader<ImageType>;
ReaderType::Pointer reader = ReaderType::New();
using ImageIOType = itk::GDCMImageIO;
ImageIOType::Pointer dicomIO = ImageIOType::New();
reader->SetImageIO(dicomIO);
using NamesGeneratorType = itk::GDCMSeriesFileNames;
NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
nameGenerator->SetUseSeriesDetails(true);
nameGenerator->SetDirectory(argv[1]);
using SeriesIdContainer = std::vector<std::string>;
const SeriesIdContainer &seriesUID = nameGenerator->GetSeriesUIDs();
SeriesIdContainer::const_iterator seriesItr = seriesUID.begin();
SeriesIdContainer::const_iterator seriesEnd = seriesUID.end();
while(seriesItr != seriesEnd)
{
std::cout << seriesItr->c_str() << std::endl;
++seriesItr;
}
std::string seriesIdentifier;
if(argc > 3)
{
seriesIdentifier = argv[3];
}
else
{
seriesIdentifier = seriesUID.begin()->c_str();
}
using FileNamesContainer = std::vector<std::string>;
FileNamesContainer fileNames;
fileNames = nameGenerator->GetFileNames(seriesIdentifier);
reader->SetFileNames(fileNames);
try
{
reader->Update();
}
catch(const itk::ExceptionObject& exp)
{
std::cerr << exp.what() << '\n';
return EXIT_FAILURE;
}
using WriterType = itk::ImageFileWriter<ImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(argv[2]);
writer->SetInput(reader->GetOutput());
writer->Update();
return EXIT_SUCCESS;
}
Then I test my stride code which split a volume with shape [900,900, 500] into 6480 small volumes with shape [50,50,25], it works with the strides which are calculated by F-order. for the higher speed, I use a raw pointer to get the 6480 small volume. It works a bit like NumPy.libs.stride_tricks.as_strided(src, dst_shape, stride), the dst_shape is [18,18,20,50,50,25]:
/*c++ code*/
small_volume[I,j,k,l,m,n] = *(image->GetBufferPointer() + i* stride[0] + j *stride[1] + k * stride[2] + l *stride[3] + m * stride[4] + m * stride[5])
in C-order, the stride is [50*900*500, 50*500,25,900*500,500,1];
in F-order, the stride is [50, 900*50, 900*900*25,1,900,900*900]
Here is a simple code to verify it.
data = np.random.randint(low=0,high=500, size=(900,900,500), dtype=np.int16) #default C-order
samll_volume = np.lib.stride_tricks.as_strided(data, (18,18,20,50,50,25),
(900*500*50*2,
500*50*2,
25*2,
900*500*2,
500*2,
1*2))
counts = 0
for i in range(18):
for j in range(18):
for k in range(20):
if np.all(samll_volume[i,j,k,...] == data[i*50:(i+1) * 50, j*50:(j+1)*50, k*25:(k+1)*25]):
counts += 1
print(counts) #6480
But when I use the C-order stride in c++ code, it does not work, the F-order stride works.
What you seem to be looking for is called RegionOfInterestImageFilter in ITK. That way you don’t have to calculate the strides by hand. If you want to apply a filter only to a region of an image, you could take a look at this example.
Hi， thanks for your suggestion.
I am a green hand with itk. When I googled the solution to the question， some websites told me the RegionOfInterestImageFilter can solve it，I tried this. But my requirement is to need all the small volumes， not just one of them. Maybe I can change the RegionType::IndexType to Iteratively generate all small volumes. But I think it could be slower than using the raw pointer. I haven’t done relevant experiments yet, it’s just my inference. I will have a try!
Thanks again for your suggestion
Take a look at this example. If you replace median filter with RegionOfInterest you can do arbitrary processing with the resulting chunk.