Incorrect resampling result with SimpleITK.sitkLinear Interpolator option

[makarovNick]

SimpleITK.sitkLinear affects result
Applying ResampleImageFilter on image produces incorrect result with interpolator = sitk.sitkLinear

To Reproduce
Steps to reproduce the behavior:

1. Operating system, version, and architecture

• OS: [Ubuntu 18.4]
  Architecture: x86_64
  CPU(s): 16

2. Programming language: Python 3.9.14]

3. Version of SimpleITK:
   tested on 2.3.1

4. How was SimpleITK installed?

• binary distribution [e.g. python -m pip install SimpleITK]

6. A minimal working example which causes the error.

import numpy as np
import SimpleITK as sitk

img_size = np.array([100, 100, 100])
img_spacing = np.array([0.15, 0.15, 0.15])
arr = np.random.randint(0, 32000, img_size)
img = sitk.GetImageFromArray(arr)
img.SetSpacing(img_spacing.tolist())

resampler = sitk.ResampleImageFilter()
resampler.SetInterpolator(sitk.sitkLinear)
resampler.SetOutputSpacing((0.3, 0.3, 0.3))
resampler.SetSize((img_size * 2).tolist())

result = resampler.Execute(img)

for eps in range(4, 18):
    new_spacing = img_spacing + (0, 0, 1 / (10**eps))
    img.SetSpacing(new_spacing)
    new_result = resampler.Execute(img)
    print(img.GetSpacing()[2], np.all(sitk.GetArrayFromImage(result) == sitk.GetArrayFromImage(new_result)))

Got:

0.15009999999999998 False
0.15001 False
0.150001 False
0.1500001 False
0.15000001 False
0.150000001 False
0.1500000001 False
0.15000000001 False
0.150000000001 False
0.1500000000001 False
0.15000000000001 False
0.150000000000001 False
0.1500000000000001 False
0.15 True

Expected behavior
For example output with sitk.sitkNearestNeighbor:

0.15009999999999998 True
0.15001 True
0.150001 True
0.1500001 True
0.15000001 True
0.150000001 True
0.1500000001 True
0.15000000001 True
0.150000000001 True
0.1500000000001 True
0.15000000000001 True
0.150000000000001 True
0.1500000000000001 True
0.15 True

Additional context
I understand possible reasons why result can vary, though I guess it will be correct if output image will be the same
in case of output_size - size * (output_spacing / input_spacing) < 0.5

[blowekamp]

Your assumption on expected behavior is incorrect.

There are likely large differences between these random pixel values which make the changes in spacing significant enough for the linear interpolation of the faction of a pixel to be off by more than an half integer. Also there is the boundary condition and interpolation occurring.

You may want to investigate the behavior of sitk::Image::EvaluateAtContinuousIndex

[makarovNick]

@blowekamp just checked with image of size 3x3x3 and values from 0 to 1, sometimes it fails to

img_size = np.array([3, 3, 3])
img_spacing = np.array([0.15, 0.15, 0.15])
arr = np.random.randint(0, 2, img_size)

[makarovNick]

something weird is happening here:

arr = np.array([[[1., 0., 1.],
                [1., 0., 0.],
                [0., 1., 1.]],

               [[0., 1., 0.],
                [1., 1., 1.],
                [0., 1., 1.]],

               [[0., 1., 1.],
                [1., 1., 1.],
                [0., 1., 0.]]])

from skimage.transform import rescale
print(rescale(arr.astype(float), (2, 2, 2)).round())

array([[[1., 1., 0., 0., 1., 1.],
        [1., 1., 0., 0., 1., 1.],
        [1., 1., 0., 0., 0., 0.],
        [1., 1., 1., 0., 0., 0.],
        [0., 0., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.]],

       [[1., 1., 0., 0., 1., 1.],
        [1., 1., 0., 0., 1., 1.],
        [1., 1., 0., 0., 0., 0.],
        [1., 1., 1., 0., 0., 0.],
        [0., 0., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.]],

       [[1., 1., 1., 1., 0., 0.],
        [1., 1., 1., 1., 0., 0.],
        [1., 1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.]],

       [[0., 0., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.]],

       [[0., 0., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.]],

       [[0., 0., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.],
        [0., 0., 1., 1., 1., 1.]]])

print(sitk.GetArrayViewFromImage(result))

array([[[1., 1., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.]],

       [[0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.]],

       [[0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0.]]], dtype=float32)

what did I do wrong?...

[makarovNick]

So I found out a way to make results for integer more stable:

resampler.SetOutputPixelType(sitk.sitkFloat32)
result = sitk.Cast(sitk.Round(resample.Execute(image)), image.GetPixelID())

Result:

0.15009999999999998 False
0.15001 False
0.150001 False
0.1500001 False
0.15000001 True
0.150000001 True
0.1500000001 True
0.15000000001 True
0.150000000001 True
0.1500000000001 True
0.15000000000001 True
0.150000000000001 True
0.1500000000000001 True
0.15 True