Upsampling with SimpleITK

clarkbab · January 4, 2024, 11:27pm

Hi there,

I’m running into an issue when upsampling an image with SimpleITK. I’ve created a simple grid to test my upsampling code, and my SimpleITK implementation seems to introduce some shift. Any ideas what I’m doing wrong? I’ve attached a reference (scipy) implementation and my SimpleITK version below.

Scipy

import numpy as np
import scipy

# Create data.
spacing = 1
data = np.array([
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.]])

# Upsample.
def resample(data, spacing, output_spacing):    
    scale = spacing / output_spacing
    resampled = scipy.ndimage.zoom(data, scale, order=1)
    return resampled

n_upsample = 5
_, axs = plt.subplots(1, n_upsample, figsize=(12, 6))
for ax in axs:
    ax.set_title(f"spacing={spacing}\nsize={data.shape[0]}")
    ax.imshow(np.transpose(data), origin='lower')
    output_spacing = spacing / 2
    data = resample(data, spacing, output_spacing)
    spacing = output_spacing

Output:

SimpleITK

import matplotlib.pyplot as plt
import numpy as np
import SimpleITK as sitk

# Create data.
spacing = (1, 1, 1)
data = np.array([
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.]])

# Upsample.
def resample(data, spacing, output_spacing):    
    image = sitk.GetImageFromArray(data)
    image.SetSpacing(spacing)
    size = image.GetSize()
    
    size_scaling = np.array(spacing) / np.array(output_spacing)
    output_size = tuple(int(s * sc) for s, sc in zip(size, size_scaling))
    output_origin = image.GetOrigin()
    
    resampled = sitk.Resample(
        image,
        output_size,
        outputOrigin=output_origin,
        outputSpacing=output_spacing,
    )
    
    return sitk.GetArrayFromImage(resampled)

n_upsample = 5
_, axs = plt.subplots(1, n_upsample, figsize=(12, 6))
for ax in axs:
    ax.set_title(f"spacing={spacing[0]}\nsize={data.shape[0]}")
    ax.imshow(np.transpose(data), origin='lower')
    output_spacing = tuple(np.array(spacing) / 2)
    data = resample(data, spacing, output_spacing)
    spacing = output_spacing

Output:

Thanks,
Brett

blowekamp · January 5, 2024, 2:11pm

Hello,

The has to do with the ITK Image being a geometric object with an origin defined at the center of the 0,0 pixel. Details of the ITK Image geometry can be reviewed here:
https://simpleitk.readthedocs.io/en/master/fundamentalConcepts.html

You can also compare the results with the ExpandImageFilter.

Also look into a SimpleITKUtilities function resize which parameters operate on pixels, and the spacing is computed with various options and features.

clarkbab · January 27, 2024, 3:01am

Thanks Bradley,

I was aware of the ITK geometry but wasn’t able to find the correct resampling parameters. In the end, adjusting the ‘output_origin’ worked as shown in the following code. The formula for calculating the output origin works for both upsampling and downsampling.

This was implemented to upsample low-resolution saliency maps from a convolutional neural network back to the input space (reversing the effect of max pooling layers) and worked well

import matplotlib.pyplot as plt
import numpy as np
import SimpleITK as sitk

# Create data.
spacing = (1, 1, 1)
data = np.array([
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.]])

# Upsample.
def resample(data, spacing, output_spacing):    
    image = sitk.GetImageFromArray(data)
    image.SetSpacing(spacing)
    size = image.GetSize()
    
    size_scaling = np.array(spacing) / np.array(output_spacing)
    output_size = tuple(int(s * sc) for s, sc in zip(size, size_scaling))
    spacing_scaling = np.array(output_spacing) / np.array(spacing)
    output_origin = tuple(sp * (s - 1) / 2 for s, sp in zip(spacing_scaling, spacing))
    
    resampled = sitk.Resample(
        image,
        output_size,
        outputOrigin=output_origin,
        outputSpacing=output_spacing,
    )
    
    return sitk.GetArrayFromImage(resampled)

n_upsample = 5
_, axs = plt.subplots(1, n_upsample, figsize=(12, 6))
for ax in axs:
    ax.set_title(f"spacing={spacing[0]}\nsize={data.shape[0]}")
    ax.imshow(np.transpose(data), origin='lower')
    output_spacing = tuple(np.array(spacing) / 2)
    data = resample(data, spacing, output_spacing)
    spacing = output_spacing

Thanks,
Brett

blowekamp · January 29, 2024, 1:46pm

Hello,

Your resampling results almost OK. Please be careful converting back and from between numpy and SimpleITK descriptive arrays as they are in different order. The SimpleITK.Utilities.resize method contains this functionality. Consider reusing this method:


data = np.array([
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 1., 1., 0., 1., 1., 0., 1.],
    [0., 0., 0., 0., 0., 0., 0., 0.],
    [0., 1., 1., 0., 1., 1., 0., 1.]])

n_upsample = 5
_, axs = plt.subplots(1, n_upsample, figsize=(12, 6))
img = sitk.GetImageFromArray(data)
for ax in axs:
    ax.set_title(f"spacing={img.GetSpacing()[0]}\nsize={data.shape[0]}")
    ax.imshow(np.transpose(data), origin='lower')
    output_size = tuple(int(s*2) for s in data.shape[::-1])
    img = sitkUtils.resize(img, output_size, sitk.sitkLinear)
    print(img.GetSpacing())
    data = sitk.GetArrayFromImage(img)

I suspect the difference between the outputs is that the SimpleITK.utilities.resize method correctly changes the origin of the image so that the extent of the image is the same. This has to do with the origin being defined as the center of the voxel not the corner.