Resampling difference between SimpleITK and scipy's zoom

Beginner Questions
,
1

I am checking a simple resampling operation to double the size of input array. SimpleITK’s resample output seems to be shifted to left and top by one voxel. Is there any other argument such as origin of the grid that I need to change to match scipy’s zoom output?

Scipy Zoom:

Sample script:

import numpy as np
from scipy.ndimage import zoom
import SimpleITK as sitk

def sitk_resample(arr, new_shape, new_spacing):

    arr_t = np.transpose(arr, (2, 1, 0))
    sitk_arr = sitk.GetImageFromArray(arr_t)

    resampled_arr = sitk.Resample(sitk_arr,
                                  new_shape,
                                  sitk.Transform(),
                                  sitk.sitkNearestNeighbor,
                                  sitk_arr.GetOrigin(),
                                  new_spacing,
                                  sitk_arr.GetDirection(),
                                  0.0,
                                  sitk_arr.GetPixelIDValue()
                                 )

    resampled_arr = sitk.GetArrayFromImage(resampled_arr)
    resampled_arr = np.transpose(resampled_arr, (2, 1, 0))

    return resampled_arr


def scipy_zoom(arr, new_shape, new_spacing):
    zoomed_arr = [np.zeros(new_shape)]
    for i in range(1, 9):
        tmp = arr.copy()
        tmp[tmp!=i] = 0
        tmp[tmp>0] = 1
        zoom_arr = zoom(tmp, zoom=new_spacing, order=0)
        zoom_arr[zoom_arr<0.5] = 0
        zoomed_arr.append(zoom_arr)

    zoomed_arr = np.argmax(zoomed_arr, axis=0)
    
    return zoomed_arr

def main():
    arr = np.zeros((4, 4, 4), dtype=np.uint8)
    arr[1, 1, 1] = 1
    arr[1, 2, 1] = 2
    arr[2, 1, 1] = 3
    arr[2, 2, 1] = 4
    arr[1, 1, 2] = 5
    arr[1, 2, 2] = 6
    arr[2, 1, 2] = 7
    arr[2, 2, 2] = 8
    print("Input_arr\n", arr)

    new_shape = [8, 8, 8]
    new_spacing = [.5, .5, .5]

    sitk_out = sitk_resample(arr, new_shape, new_spacing)

    zoom_out = scipy_zoom(arr, new_shape, 1./np.array(new_spacing))

    print("\n\tsitk_out\t\tzoom_out\n")

    for i in range(new_shape[0]):
        for j in range(new_shape[0]):
            print(sitk_out[i,j,:], "\t", zoom_out[i,j,:])
        print()


if __name__ == "__main__":
    main()

Sample Output:

Input_arr
 [[[0 0 0 0]
  [0 0 0 0]
  [0 0 0 0]
  [0 0 0 0]]

 [[0 0 0 0]
  [0 1 5 0]
  [0 2 6 0]
  [0 0 0 0]]

 [[0 0 0 0]
  [0 3 7 0]
  [0 4 8 0]
  [0 0 0 0]]

 [[0 0 0 0]
  [0 0 0 0]
  [0 0 0 0]
  [0 0 0 0]]]

	sitk_out		      zoom_out

[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 1 1 5 5 0 0 0] 	 [0 0 0 0 0 0 0 0]
[0 1 1 5 5 0 0 0] 	 [0 0 0 0 0 0 0 0]
[0 2 2 6 6 0 0 0] 	 [0 0 0 0 0 0 0 0]
[0 2 2 6 6 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 1 1 5 5 0 0 0] 	 [0 0 0 0 0 0 0 0]
[0 1 1 5 5 0 0 0] 	 [0 0 1 1 5 5 0 0]
[0 2 2 6 6 0 0 0] 	 [0 0 1 1 5 5 0 0]
[0 2 2 6 6 0 0 0] 	 [0 0 2 2 6 6 0 0]
[0 0 0 0 0 0 0 0] 	 [0 0 2 2 6 6 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 3 3 7 7 0 0 0] 	 [0 0 0 0 0 0 0 0]
[0 3 3 7 7 0 0 0] 	 [0 0 1 1 5 5 0 0]
[0 4 4 8 8 0 0 0] 	 [0 0 1 1 5 5 0 0]
[0 4 4 8 8 0 0 0] 	 [0 0 2 2 6 6 0 0]
[0 0 0 0 0 0 0 0] 	 [0 0 2 2 6 6 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 3 3 7 7 0 0 0] 	 [0 0 0 0 0 0 0 0]
[0 3 3 7 7 0 0 0] 	 [0 0 3 3 7 7 0 0]
[0 4 4 8 8 0 0 0] 	 [0 0 3 3 7 7 0 0]
[0 4 4 8 8 0 0 0] 	 [0 0 4 4 8 8 0 0]
[0 0 0 0 0 0 0 0] 	 [0 0 4 4 8 8 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 3 3 7 7 0 0]
[0 0 0 0 0 0 0 0] 	 [0 0 3 3 7 7 0 0]
[0 0 0 0 0 0 0 0] 	 [0 0 4 4 8 8 0 0]
[0 0 0 0 0 0 0 0] 	 [0 0 4 4 8 8 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]
[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]
2

Hello @couzulli,

The voxel/pixel location in ITK/SimpleITK is not at the corner of the grid, it is the center of the element. I suspect this is the reason for the discrepancy you are experiencing. Please see the SimpleITK Fundamental Concepts page for details.

Also, note that the original image and the resampled image in the code do not occupy the same physical region. The resampled image extends 0.5mm out beyond the bounds of the original. For the resampled image to occupy the same physical region as the original and have a spacing of 0.5, its size would be 7 voxels.

3

@zivy Thank you for the clarification!

Per your comment and the information from the link which states that “image’s physical extent starts half a voxel before the origin”, I updated origin argument in the above code to [-.5, -.5, -.5] and the results match now.