Interpolation error(?) with zero translation in Affine

[ntatsisk]

Hello,

Doing an Affine transform with translate_params=0 leads some kind of interpolation error between the input and the output tensor, which in theory should be identical. When setting mode="nearest" there is no such error as expected.

Code to reproduce:

import torch
from monai.transforms import Affine

def calculate_max_error(tensor_shape, iterations=1000):
    # Zero translation
    affine_transform = Affine(translate_params=0, padding_mode="zeros")

    # Calculate max absolute error between input_tensor and output_tensor
    errors = []
    for i in range(iterations):
        input_tensor = 1000*torch.rand(tensor_shape)
        output_tensor, output_affine = affine_transform(input_tensor)
        diff = output_tensor - input_tensor
        errors.append(diff.abs().max())

    return torch.tensor(errors).max().item()

def print_error_for_shape(tensor_shape):
    print("Shape: {} - error: {}".format(tensor_shape, calculate_max_error(tensor_shape)))

print_error_for_shape((1, 2, 2, 2))
print_error_for_shape((1, 1, 1, 3))
print_error_for_shape((1, 3, 3, 3))
print_error_for_shape((1, 1, 1, 4))
print_error_for_shape((1, 4, 4, 4))
print_error_for_shape((1, 8, 8, 8))
print_error_for_shape((1, 9, 9, 9))
print_error_for_shape((1, 25, 25, 25))
print_error_for_shape((1, 27, 27, 27))
print_error_for_shape((1, 50, 50, 50))

which gives the following output:

Shape: (1, 2, 2, 2) - error: 0.0
Shape: (1, 1, 1, 3) - error: 0.0
Shape: (1, 3, 3, 3) - error: 0.0
Shape: (1, 1, 1, 4) - error: 6.103515625e-05
Shape: (1, 4, 4, 4) - error: 0.00018310546875
Shape: (1, 8, 8, 8) - error: 0.00030517578125
Shape: (1, 9, 9, 9) - error: 0.0
Shape: (1, 25, 25, 25) - error: 0.0013427734375
Shape: (1, 27, 27, 27) - error: 0.005126953125
Shape: (1, 50, 50, 50) - error: 0.01068115234375

Is this normal? Why is this interpolation artifact created in the first place? Equally puzzling: why some of the cases give zero error?

I am using:

MONAI version: 1.1.dev2248
Numpy version: 1.23.5
Pytorch version: 1.13.0+cu117
MONAI flags: HAS_EXT = False, USE_COMPILED = False, USE_META_DICT = False
MONAI rev id: 3400bd91422ccba9ccc3aa2ffe7fecd4eb5596bf

Thank you!

--
xref for our convenience: InsightSoftwareConsortium/ITKElastix#126

[wyli]

thanks for reporting, with float64 it'll be more precise

affine_transform = Affine(translate_params=0, padding_mode="zeros", dtype=torch.float64)

Shape: (1, 2, 2, 2) - error: 0.0
Shape: (1, 1, 1, 3) - error: 0.0
Shape: (1, 3, 3, 3) - error: 0.0
Shape: (1, 1, 1, 4) - error: 0.0
Shape: (1, 4, 4, 4) - error: 0.0
Shape: (1, 8, 8, 8) - error: 0.0
Shape: (1, 9, 9, 9) - error: 0.0
Shape: (1, 25, 25, 25) - error: 5.949706330908577e-13
Shape: (1, 27, 27, 27) - error: 0.0
Shape: (1, 50, 50, 50) - error: 6.532264746478278e-12

[ntatsisk]

Thanks @wyli for the reply! The solution that you provided indeed addresses the problem for zero translation. However, when I eventually start using non-zero values I notice higher interpolation error compared to e.g. ITK resampling. I might be wrong though. Is it fair to say that when normalized=False (default value) then Affine does actually two consecutive interpolations? One for the "coordinate normalization" and one for requested transformation/resampling?

Any thought or resource/link will help me a lot because I am working on an affine matrix bridge between ITK and MONAI.

[wyli]

that's an interesting question, coordinate normalization is remapping it to -1 and 1 so that the pytorch resampling backend can understand it

I made some tests comparing monai and itkpython

MONAI/tests/test_meta_affine.py

Lines 163 to 175 in e270741

	def test_linear_consistent_dict(self, xform_cls, input_dict, atol):
	"""xform cls testing itk consistency"""
	img = LoadImaged(keys, image_only=True, simple_keys=True)({keys[0]: FILE_PATH, keys[1]: FILE_PATH_1})
	img = EnsureChannelFirstd(keys)(img)
	ref_1 = {k: _create_itk_obj(img[k][0], img[k].affine) for k in keys}
	output = self.run_transform(img, xform_cls, input_dict)
	ref_2 = {k: _create_itk_obj(output[k][0], output[k].affine) for k in keys}
	expected = {k: _resample_to_affine(ref_1[k], ref_2[k]) for k in keys}
	# compare ref_2 and expected results from itk
	diff = {k: np.abs(itk.GetArrayFromImage(ref_2[k]) - itk.GetArrayFromImage(expected[k])) for k in keys}
	avg_diff = {k: np.mean(diff[k]) for k in keys}
	for k in keys:
	self.assertTrue(avg_diff[k] < atol, f"{xform_cls} avg_diff: {avg_diff}, tol: {atol}")

another related unit tests is

MONAI/tests/test_warp.py

Line 110 in 99b9dd7

def test_itk_benchmark(self):

please let me know if you have specific code examples for your questions, I'm happy to look into them.