# C order or F order？

**URL:** https://discourse.itk.org/t/c-order-or-f-order/5287
**Category:** Beginner Questions
**Tags:** dicom
**Created:** [August 23, 2022, 2:35am UTC](https://discourse.itk.org/t/c-order-or-f-order/5287 "2022-08-23T02:35:47Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![liguangops](https://discourse.itk.org/user_avatar/discourse.itk.org/liguangops/32/2900_2.png) [@liguangops](https://discourse.itk.org/u/liguangops)
#### Post date: [August 23, 2022, 2:35am UTC](https://discourse.itk.org/t/c-order-or-f-order/5287/1 "2022-08-23T02:35:47Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [August 23, 2022, 1:58pm UTC](https://discourse.itk.org/t/c-order-or-f-order/5287/2 "2022-08-23T13:58:01Z")

</div>

You should take a look at [Data Representation](https://itk.org/ITKSoftwareGuide/html/Book1/ITKSoftwareGuide-Book1ch4.html) 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.

---

<div class="post-metadata">

### Author: ![liguangops](https://discourse.itk.org/user_avatar/discourse.itk.org/liguangops/32/2900_2.png) [@liguangops](https://discourse.itk.org/u/liguangops)
#### Post date: [August 25, 2022, 1:43am UTC](https://discourse.itk.org/t/c-order-or-f-order/5287/3 "2022-08-25T01:43:10Z")

</div>

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:

```python
import nibabel as nib

data = nib.load("./my_data.nii.gz").get_fdata()
print(data.flags)

```

```nohighlight
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:

```nohighlight
/*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

---

<div class="post-metadata">

### Author: ![matt.mccormick](https://discourse.itk.org/user_avatar/discourse.itk.org/matt.mccormick/32/7_2.png) [@matt.mccormick](https://discourse.itk.org/u/matt.mccormick)
#### Post date: [August 25, 2022, 1:57am UTC](https://discourse.itk.org/t/c-order-or-f-order/5287/4 "2022-08-25T01:57:50Z")

</div>

Hi,

nibabel reads in nifti into F-order arrays.

ITK reads in nifti (and other images) into C-order. For comparison with nibabel:

```auto
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.

---

<div class="post-metadata">

### Author: ![liguangops](https://discourse.itk.org/user_avatar/discourse.itk.org/liguangops/32/2900_2.png) [@liguangops](https://discourse.itk.org/u/liguangops)
#### Post date: [August 25, 2022, 9:05am UTC](https://discourse.itk.org/t/c-order-or-f-order/5287/5 "2022-08-25T09:05:34Z")

</div>

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:

```nohighlight
#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]:

```nohighlight
/*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.

```python
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.

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [August 25, 2022, 2:24pm UTC](https://discourse.itk.org/t/c-order-or-f-order/5287/6 "2022-08-25T14:24:11Z")

</div>

What you seem to be looking for is called [RegionOfInterestImageFilter](https://itk.org/Doxygen/html/classitk_1_1RegionOfInterestImageFilter.html) 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](https://examples.itk.org/src/core/common/applyafilteronlytoaspecifiedregionofanimage/documentation).

---

<div class="post-metadata">

### Author: ![liguangops](https://discourse.itk.org/user_avatar/discourse.itk.org/liguangops/32/2900_2.png) [@liguangops](https://discourse.itk.org/u/liguangops)
#### Post date: [August 26, 2022, 1:46am UTC](https://discourse.itk.org/t/c-order-or-f-order/5287/7 "2022-08-26T01:46:24Z")

</div>

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

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [August 26, 2022, 1:09pm UTC](https://discourse.itk.org/t/c-order-or-f-order/5287/8 "2022-08-26T13:09:21Z")

</div>

Take a look at this [example](https://github.com/OpenImaging/demos/blob/master/ProcessImageChunks/Code.py). If you replace median filter with RegionOfInterest you can do arbitrary processing with the resulting chunk.