Input Orientation and Axis Indexing in PyRadiomics (SimpleITK vs NumPy)

Community Radiomics
1

Is there a specific orientation (axis order) expected by Pyradiomics when processing input images? I’m asking in the context of the opposite axes indexing between numpy/nibabel and SimpleITK :

SimpleITK and numpy indexing access is in opposite order!

SimpleITK: image[x,y,z]
numpy: image_numpy_array[z,y,x]

Since Pyradiomics uses SimpleITK internally, should the input axis order be interpreted as x, y, z?

If so, does this imply the following mapping for the force2Ddimension?

  • 0 == x == axial
  • 1 == y == coronal
  • 2 == z == sagittal
force2Ddimension details 
force2Ddimension: 0  # axial slices, for coronal slices, use dimension 1 and for sagittal, dimension 2

Does Pyradiomics automatically reorder axes internally to match SimpleITK’s conventions (x, y, z), or is it the user’s responsibility to ensure the input image has the correct orientation?

2

I don’t recall exactly what pyradiomics does, but I believe it resamples the segmentation to match the geometry of the source image (or maybe the other way around).

Keep in mind that medical images can be acquired in many different orientations, so you should not make assumptions about the relationship between world space and index space. In Slicer that information is in the ijkToRAS matrix and related ones. See this info: https://slicer.readthedocs.io/en/latest/user_guide/coordinate_systems.html

3

You’are right! This is exactly why I want to reorient all my images to the same orientation before calling pyradiomics. Please note that I’m calling pyradiomics from the CLI, not via Slicer.

I should have mentioned in my first comment what I meant by axial, coronal, and sagittal, so clarifying it now:

  • 0 == x == axial == Superior-Inferior
  • 1 == y == coronal == Anterior-Posterior
  • 2 == z == sagittal == Right-Left

(source for numbers 0, 1, 2 and planes axial, coronal, saggital; source for x, y, z)

Which effectively translates into SAR (Superior-Anterior-Right) orientation, right?

This is however contradicted by this source:

ITK uses the LPS convention.

4

Just to restate what I said earlier, the 0, 1, 2 and x, y, z in image space that you are referring to do not correspond to any particular world space direction or slicing plane. That mapping is defined by an affine matrix.

You’ll need to study the pyradiomics code, but as I recall it takes care of the resampling for you.

But if you want to make sure all your test data is anatomically aligned, you need to be sure to understand and apply the transformations I’ve described.

1 Like
5

Thank you!

Just to confirm I understand it correctly, let me illustrate what I mean by an example:

  • I have two types of 3D MRI images: T2w axial images and T2w sagittal images.
  • T2w axial images have a high in-plane resolution in the axial plane and thick slices in the superior-inferior (S-I) axis.
  • T2w sagittal images have a high in-plane resolution in the sagittal plane and thick slices in the right-left (R-L) axis.
  • I reoriented both images to the RAS coordinate system

When inspecting the reoriented images using SimpleITK, I get the following:

>>> import SimpleITK as sitk
>>> image_ax = sitk.ReadImage('sub-001_acq-ax_T2w_RAS.nii.gz')
>>> image_ax.GetSpacing()
(0.5624999403953552, 0.5625000596046448, 3.7499988079071045)


>>> import SimpleITK as sitk
>>> image_sag = sitk.ReadImage('sub-001_acq-sag_T2w_RAS.nii.gz')
>>> image_sag.GetSpacing()
(3.4499998092651367, 0.48828125, 0.4882812798023224)
  • the axial image (sub-001_acq-ax_T2w_RAS.nii.gz) has correctly the highest spacing (i.e., thick slices) in the last dimension (S-I)
  • the sagittal image (sub-001_acq-sag_T2w_RAS.nii.gz) has correctly the highest spacing in the first dimension (R-L)

Then it matches the pyradiomics example, which states:

force2Ddimension: 0 # axial slices, for coronal slices, use dimension 1 and for sagittal, dimension 2.

6

Just popping in to try to clarify what @valosekj is referring to by the above quote.

We’re coming from the world of NIFTI1 images, nibabel, and FSL (specifically tools like the FSLeyes image viewer and fslhd). As far as I understand:

  • For nibabel and FSL, affine matrices are assumed to always map the voxel array to RAS+ world space.
  • However, in voxel space, the array can be “reoriented” (i.e. axes transposition and inversion). For FSL, this can be done with fslswapdim. In order to preserve the “always maps to RAS+” assumption, if the data array is changed, then the affine will be changed in kind.
  • The tool fslhd reports the header information in NIFTI1 image files. Part of the output information is the axis ordering corresponding to the qform/sform affine matrices. For example, on a sample image I’ve picked out, the command reports:
    sto_xyz:1       0.000000 0.000000 1.000000 -32.219284 
sto_xyz:2       -1.000000 0.000000 0.000000 31.786606 
sto_xyz:3       0.000000 1.000000 0.000000 -19.321671 
sto_xyz:4       0.000000 0.000000 0.000000 1.000000 
sform_xorient   Anterior-to-Posterior
sform_yorient   Inferior-to-Superior
sform_zorient   Left-to-Right
  • Because of this information, our lab will refer to such images in shorthand as “AIL-/PSR+ images”. In other words, for this example image, the voxel-space array has an orientation associated with it (0 = x = AP, 1 = y = IS, 2 = z = LR), but only because of the inherent assumption that the affine will always map to RAS+ physical space.
    • Side note: This can be intuited via the above affine. The -1 in the first column flips the so-called “PSR+ image” to be ASR+, and the non-diagonal matrix reorients the axes from ASR+ to RAS+. This kind of intuition is part of why we find it helpful to describe voxel space in physical terms.
  • ((For reference, the reason why we describe images in “voxel space terms” is because we do a lot of Python-based image processing where we operate directly on the array using tools such as numpy, skimage, scipy, dipy, etc. Since we process the array using tools that ignore the context of the image header, it is important for us to know which array axes map to which physical axes (which then get universally transformed into RAS+ physical coordinates by the affine matrix).))

So, I suppose what @valosekj is truly asking is:

  • Assuming the images have an affine that universally maps the array to RAS+ physical coordinates, does axis order matter in voxel space?
  • In other words, will “reorienting” the image in voxel space (i.e. modifying the data array + header in tandem) make any difference for pyradiomics, given the assumption that the affine will always map the array to RAS+?

Kind regards,
Joshua

7

OK, I think we have found the source of our confusion:

In this code snippet, the comment implies that 0 = axial, 1 = coronal, and 2 = sagittal. This was confusing for us because we are used to RAS+ orientation (see my previous post).

However, I then found a documentation page on force2Ddimension which says:

force2Ddimension [0]: int, range 0-2. Specifies the ‘slice’ dimension for a by-slice feature extraction. A value of 0 represents the native acquisition plane for the images (usually axial for CT and axial, coronal or sagittal for MRI).

In other words, it seems that the code snippet’s comments were written from a CT-specific perspective (0 = axial), whereas we’re coming from an MRI perspective (0 = sagittal = LR, because our affines are RAS+ based).

I think the confusion has been cleared up? In our case, it seems that we should (A) make sure our voxel-space images are in RAS+ orientation (possibly irrelevant), (B) that the affine transforms the image into RAS+ world space (more relevant?), and (C) select 2 to specify the SI axis (i.e. axial slices). (Note: using 2 for axial slices is what we are used to in our processing.)

Kind regards,
Joshua

8

@JoostJM or @fedorov can you jump in here?