Need to transpose array obtained from itk instead of reversing of itk.size() in ITKReader

[take5v]

Description

In ITKReader, itk.size(img) should not be reversed since it outputs correct size of an image. Instead, itk.array_view_from_image(img, keep_axes=False) provides the transposed array and such numpy array must be transposed because of the differences of storing the data in numpy and itk.

This PR is wrong: #1045. See also InsightSoftwareConsortium/ITK#2006

Expected behavior

1. Need to remove the reverse for itk.size.
2. Transpose array obtained from itk.array_view_from_image.

[aylward]

Note that work has been done to now automate the conversion of Pytorch tensors to ITK images. See this commit that happened today:

InsightSoftwareConsortium/ITK@c42436f

Additionally, ITK Python now supports 4D images even when those 4D images are actually representing vectors of 3D images or 3D images of vectors.

Much testing is needed, and the MONAI code will perhaps need updating, but I think fixing ITK was the right first step.

@thewtex - Thanks!

[thewtex]

@take5v we need to avoid itk.array_view_from_image(img, keep_axes=False) because it results in data storage that is different from NumPy and PyTorch. We need to avoid unnecessary related performance penalties.

#1045 adds the correct reversal in itk.size to correspond to the index ordering of a C-order NumPy array, described in InsightSoftwareConsortium/ITK#2006. ITK uses IJK ordering, which C-order NumPy arrays are KJI.

[Nic-Ma]

Hi @thewtex and @aylward ,

Thanks for your detailed explanation here.
Could you please help submit a MR to MONAI to update something once ITK fixed the order?
Let's keep this ticket open to track it.

Thanks in advance.

[take5v]

Hi @aylward, @thewtex, @Nic-Ma!

I have made a simple test in order to check ITKReader currently implemented in MONAI. Beforehand, I converted nifty anatomical.nii (availble in MONAI test data) to DICOM.

TEST_CASE = [DICOM_PATH, NIFTY_PATH]

class TestITKLoader(unittest.TestCase):
    @parameterized.expand([TEST_CASE])
    def test_array_shape(self, dicom_path, nifty_path):
        itk_reader = ITKReader()
        itk_img = itk_reader.read(dicom_path)
        itk_array, itk_meta_data = itk_reader.get_data(itk_img)

        nifty_reader = NibabelReader()
        nifty_img = nifty_reader.read(nifty_path)
        nifty_array, nifty_metadata = nifty_reader.get_data(nifty_img)

        self.assertTupleEqual(itk_array.shape, nifty_array.shape)
        self.assertTupleEqual(tuple(
            itk_meta_data['spatial_shape']), tuple(nifty_metadata['spatial_shape']))

Currently this test is failed:
itk_array shape is (25, 41, 33) whereas nifty_array shape is (33, 41, 25)
itk_meta_data['spatial_shape'] is (25, 41, 33) whereas nifty_metadata['spatial_shape'] is (33, 41, 25)

That was a reason why I suggested to remove shape reversing and to add transposing of itk.array_view_from_image which is, at least, fixes a bug in ITKReader. As @thewtex and @aylward pointed out, it may be a better solution for this.

Please find attached DICOM files of the test image (anatomical.nii):
anatomical_DICOM.zip

Additionally, this issue should be marked as a bug.

[thewtex]

@Nic-Ma @take5v just a note that I will take some time to follow up on this over the next few days.

[thewtex]

To follow-up, I made some progress on this, but I will come back to it after I return in the new year. Happy pre-2021!

[Nic-Ma]

Hi @thewtex ,

Have you got a chance to make some progress on this ticket?

Thanks.

[thewtex]

@Nic-Ma I know there is a lot of confusion in the community around indexing, why the order is the way it is, how indexing order relates to convension, memory layout, and performance. And why the order is not the way we sometimes want it to be.

In summary, the order provided by ITKReader is correct. It is consistent with NumPy index conventions, the memory layout, and optimal performance. The current NibabelReader returns a reversed shape that some folks in the neuroimaging community are used to, and this can work when working within the scope of nipy tools. However, it conflicts with most other image processing tools in the ecosystem, e.g. scikit-image, itk, imagej. And, while it may not make a big difference for neuroimaging datasets, which are often nearly anisotropically sampled, it will have unnecessary performance penalties for other imaging datasets, where in-plane resolution is often much higher than out-of-plane resolution, and 2D models are applied.

I am working on documentation and examples that explain this, but it will take some time.

[Nic-Ma]

Hi @thewtex ,

Really appreciate your help here!

Thanks in advance.

[brad-t-moore]

I'm not sure if this is related to a bug I'm currently running into (I haven't submitted an issue yet).

print(img.GetSpacing())
print(img.GetLargestPossibleRegion().GetSize())
type(img)
(55, 512, 512)
itkVectorD3 ([0.976562, 0.976562, 5])
itkSize3 ([512, 512, 55])
itk.itkImagePython.itkImageD3

And then when I'm trying to apply a Spacingd to that image (as loaded through MONAI) it's misapplying. The output is from print statements in the Compose apply transform loop (type(input_), type(transform), input['image'].shape):

val_transforms = Compose(
    [
        LoadImaged(keys=["image", "label"], reader='ITKReader'),
        AddChanneld(keys=["image", "label"]),
        Orientationd(keys=["image", "label"], axcodes="RAS"),
        Spacingd(keys=["image", "label"], pixdim=(1.0, 1.0, 1.0), mode=("bilinear", "nearest")),
        ScaleIntensityRanged(
            keys=["image"], a_min=-57, a_max=164, b_min=0.0, b_max=1.0, clip=True,
        ),
        CropForegroundd(keys=["image", "label"], source_key="image"),
        ToTensord(keys=["image", "label"]),
    ]
)
<class 'dict'>
<class 'monai.transforms.io.dictionary.LoadImaged'>
<class 'dict'>
<class 'monai.transforms.utility.dictionary.AddChanneld'>
image.shape:  (55, 512, 512)
<class 'dict'>
<class 'monai.transforms.spatial.dictionary.Orientationd'>
image.shape:  (1, 55, 512, 512)
<class 'dict'>
<class 'monai.transforms.spatial.dictionary.Spacingd'>
image.shape:  (1, 55, 512, 512)
<class 'dict'>
<class 'monai.transforms.intensity.dictionary.ScaleIntensityRanged'>
image.shape:  (1, 54, 500, 2556)

It's taking that last index and multiplying it by 5 due to the "z-axis" being a originally around 5mm in spacing. I'm calling this a bug because somewhere it's correctly calculating the spacing ratio but then applying it to the wrong index.

This is the state of the input_ in the Compose transform before applying Spacingd:

ipdb>  input_['image'].shape
(1, 55, 512, 512)
ipdb>  input_['image_meta_dict']
{'0008|0016': '1.2.840.10008.5.1.4.1.1.7.3', '0008|0018': '1.2.826.0.1.3680043.2.1125.1.90455018032395752885544164701193819', '0008|0020': '20210128', '0008|0030': '164956.536161 ', '0008|0050': '', '0008|0060': 'OT', '0008|0090': '', '0010|0010': '', '0010|0020': '', '0010|0030': '', '0010|0040': '', '0020|000d': '1.2.826.0.1.3680043.2.1125.1.28882890755346393261644168155624487', '0020|000e': '1.2.826.0.1.3680043.2.1125.1.68989042725981558895947946072747819', '0020|0010': '', '0020|0011': '', '0020|0013': '', '0020|0052': '1.2.826.0.1.3680043.2.1125.1.19359276313777369875433246900418117', '0028|0002': '1', '0028|0004': 'MONOCHROME2 ', '0028|0008': '55', '0028|0009': '(5200,9230)', '0028|0010': '512', '0028|0011': '512', '0028|0100': '16', '0028|0101': '16', '0028|0102': '15', '0028|0103': '0', '0028|1052': '-16383.75 ', '0028|1053': '15.999755859375 ', '0028|1054': 'US', 'origin': array([499.02319336, 499.02319336,   0.        ]), 'spacing': array([0.97656202, 0.97656202, 5.        ]), 'direction': array([[-1.,  0.,  0.],
       [ 0., -1.,  0.],
       [ 0.,  0.,  1.]]), 'original_affine': array([[ -0.97656202,   0.        ,   0.        , 499.02319336],
       [  0.        ,  -0.97656202,   0.        , 499.02319336],
       [  0.        ,   0.        ,   5.        ,   0.        ],
       [  0.        ,   0.        ,   0.        ,   1.        ]]), 'affine': array([[ 9.76562023e-01,  0.00000000e+00,  0.00000000e+00,
         4.46288844e+02],
       [ 0.00000000e+00,  9.76562023e-01,  0.00000000e+00,
        -4.76836760e-07],
       [ 0.00000000e+00,  0.00000000e+00,  5.00000000e+00,
         0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         1.00000000e+00]]), 'spatial_shape': array([ 55, 512, 512]), 'filename_or_obj': './PVM-Girder-TestData/spleen_10.dcm'}

And here it is afterward:

ipdb>  input_['image'].shape
(1, 54, 500, 2556)
ipdb>  input_['image_meta_dict']
{'0008|0016': '1.2.840.10008.5.1.4.1.1.7.3', '0008|0018': '1.2.826.0.1.3680043.2.1125.1.90455018032395752885544164701193819', '0008|0020': '20210128', '0008|0030': '164956.536161 ', '0008|0050': '', '0008|0060': 'OT', '0008|0090': '', '0010|0010': '', '0010|0020': '', '0010|0030': '', '0010|0040': '', '0020|000d': '1.2.826.0.1.3680043.2.1125.1.28882890755346393261644168155624487', '0020|000e': '1.2.826.0.1.3680043.2.1125.1.68989042725981558895947946072747819', '0020|0010': '', '0020|0011': '', '0020|0013': '', '0020|0052': '1.2.826.0.1.3680043.2.1125.1.19359276313777369875433246900418117', '0028|0002': '1', '0028|0004': 'MONOCHROME2 ', '0028|0008': '55', '0028|0009': '(5200,9230)', '0028|0010': '512', '0028|0011': '512', '0028|0100': '16', '0028|0101': '16', '0028|0102': '15', '0028|0103': '0', '0028|1052': '-16383.75 ', '0028|1053': '15.999755859375 ', '0028|1054': 'US', 'origin': array([499.02319336, 499.02319336,   0.        ]), 'spacing': array([0.97656202, 0.97656202, 5.        ]), 'direction': array([[-1.,  0.,  0.],
       [ 0., -1.,  0.],
       [ 0.,  0.,  1.]]), 'original_affine': array([[ -0.97656202,   0.        ,   0.        , 499.02319336],
       [  0.        ,  -0.97656202,   0.        , 499.02319336],
       [  0.        ,   0.        ,   5.        ,   0.        ],
       [  0.        ,   0.        ,   0.        ,   1.        ]]), 'affine': array([[ 1.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         4.46288844e+02],
       [ 0.00000000e+00,  1.00000000e+00,  0.00000000e+00,
        -4.76836760e-07],
       [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00,
         0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         1.00000000e+00]]), 'spatial_shape': array([ 55, 512, 512]), 'filename_or_obj': './PVM-Girder-TestData/spleen_10.dcm'}

@thewtex

EDIT: this is using the release packages via pip (itk 5.1.2, monai 0.4.0)

[thewtex]

@brad-t-moore original_affine is not correct in this case.

print(img.GetSpacing())
itkVectorD3 ([0.976562, 0.976562, 5])

This is ITK, Fortran order.

image.shape: (55, 512, 512)

This is NumPy, C order.

But the affine matrix is not in NumPy C order.

'original_affine': array([[ -0.97656202,   0.        ,   0.        , 499.02319336],
       [  0.        ,  -0.97656202,   0.        , 499.02319336],
       [  0.        ,   0.        ,   5.        ,   0.        ],
       [  0.        ,   0.        ,   0.        ,   1.        ]]

[thewtex]

this is using the release packages via pip (itk 5.1.2, monai 0.4.0)

We will soon release ITK 5.2, RC 2, which will make it possible to simplify the MONAI code. I will up update MONAI and look into this original_affine issue.