Resample produces blurry results when just cropping

Theoretialperson · October 5, 2021, 8:37am

I found out I was getting slightly blurry results when using Resample just with cropping:

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

def myshow(img, title=None, margin=0.05, dpi=80):
    nda = sitk.GetArrayViewFromImage(img)[:,:,62]
    spacing = img.GetSpacing()
        
    ysize = nda.shape[0]
    xsize = nda.shape[1]
      
    figsize = ((1 + margin) * ysize / dpi, (1 + margin) * xsize / dpi)

    fig = plt.figure(title, figsize=figsize, dpi=dpi)
    ax = fig.add_axes([margin, margin, 1 - 2*margin, 1 - 2*margin])
    
    extent = (0, xsize*spacing[1], 0, ysize*spacing[0])

    t = ax.imshow(nda,
            extent=extent,
            interpolation='none',
            cmap='gray',
            origin='lower',
            vmin=0.0, vmax=1.0)
    
    if(title):
        plt.title(title)

    plt.show()

if __name__ == '__main__':
    scale = 1.0
    spacing = 0.2
    image_size = 250
    patch_size = 125
    sampling = sitk.sitkBSplineResamplerOrder3
    #sampling = sitk.sitkNearestNeighbor

    grid = ((np.indices(np.repeat(image_size, 3)) // 4).sum(axis=0) % 2).astype(float)
    grid = sitk.GetImageFromArray(grid)
    grid.SetSpacing(np.ones(3) * spacing)
    grid.SetOrigin((30.9, 10.2, 28.4))

    #myshow(grid, 'Grid Input')

    outputSpacing   = scale * np.array(grid.GetSpacing())
    outputDirection = np.diag([1.0, 1.0, 1.0]) @ np.array(grid.GetDirection()).reshape((3,3))

    imageCenter = 0.5 * (image_size - 1) * np.array(grid.GetSpacing())
    patchCenter = outputDirection.dot(0.5 * (patch_size - 1) * outputSpacing)
    translation = imageCenter - patchCenter
    print(translation)

    outputOrigin = np.array(grid.GetOrigin())# + translation

    resampled1 = sitk.Resample(
        grid, (np.repeat(patch_size, 3)).tolist(),
        sitk.Transform(3, sitk.sitkIdentity), sampling,
        defaultPixelValue=255,
        outputSpacing=outputSpacing, outputOrigin=outputOrigin, outputDirection=outputDirection.flatten()
    )
    myshow(resampled1, 'Resampled')

Even if you disable cropping by making patch_size = 250, the result is blurry, as can be seen in the close-up below.

Why is this happening? Is it due to numerical accuracy? Can I alleviate this problem?

zivy · October 5, 2021, 1:31pm

Hello @Theoretialperson,

This isn’t a bug it is a feature .

When resampling a signal with high frequency content like the grid created here you are highlighting several aspects of resampling:

The effect of the interpolating function, in this case sitkBSplineResamplerOrder3 uses values from a larger neighborhood and thus you get the blurring / gradual change in value. Switching to sitkLinear results in a sharp grid.
The effect of the resampling rate a combination of the transformation and the choice of output grid size, spacing, and direction. See the difference between running the code with the + translation commented out vs. enabling it (using sitkLinear).

The quality of the resampled signal is always dependent on the interpolator the original signal content, frequency wise and the resampling rate (Nyquist theorem).

To crop an image/signal it is best to directly do it (grid[start_x:end_x, start_y:end_y, start_z:end_z]).

Theoretialperson · October 5, 2021, 2:17pm

@zivy You are right, I was testing with sitkNearestNeighbor but with the translation, which introduced aliasing issues, but this does not happen without the translation part. The problem originates mostly from the origin of the grid, which can not be represented exactly.

I asked the question originally because I’m trying to find the best way to build a patch extractor/augmenter for medical images (e.g. CT scans). One way is to extract a slightly larger patch and apply any transformations to it and crop to the final patch size. The other way is to just use Resample with outputDirection, outputSpacing, etc. This approach is faster, but the problem is that can result in more blurry images when no transformation is applied. I would have to implement logic to use normal cropping when no transforms are applied.

lassoan · October 6, 2021, 1:33am

Have a look at torchio and monai. These libraries may already provide the patch extraction and data augmentation features that you need and they are specifically developed for medical imaging. They already use ITK and the developers are open for collaboration, so if you miss something then it is probably better to contribute to these packages than redevelop everything from scratch.

Theoretialperson · October 14, 2021, 8:55am

@zivy

I still thought it’s weird you’re getting a blurred result when not applying any transform so I made a quick test:

import SimpleITK as sitk
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import map_coordinates

if __name__ == '__main__':
    image_size = np.array([64, 64])

    grid_array = ((np.indices(image_size) // np.array([4, 4])[:,None,None]).sum(axis=0) % 2).astype(float)

    plt.subplot(131)
    plt.gca().set_title('original')
    plt.imshow(grid_array, cmap='gray')

    # sitk Resample
    image = sitk.GetImageFromArray(grid_array)
    resampled_image = sitk.Resample(image, sitk.Transform(2, sitk.sitkIdentity), sitk.sitkBSplineResamplerOrder3, 1.0)
    resampled_grid = sitk.GetArrayViewFromImage(resampled_image)

    plt.subplot(132)
    plt.gca().set_title('sitk cubic spline')
    plt.imshow(resampled_grid, cmap='gray')
    # sitk Resample

    # scipy map_coordinates
    si, sj = grid_array.shape
    coords = np.mgrid[0:si, 0:sj]
    resampled_grid = map_coordinates(grid_array, coords, order=3, cval=1.0, mode='constant')

    plt.subplot(133)
    plt.gca().set_title('scipy cubic spline')
    plt.imshow(resampled_grid, cmap='gray')
    # scipy map_coordinates

    plt.show()

And this is the result:

As you can see, although scipy also does cubic interpolation, without any transform, the result is identical to the original image. Do you still think this is a feature and not a bug? If yes, is there a way to achieve the same results as scipy’s implementation?

Theoretialperson · October 14, 2021, 9:45am

I also tried using itk for the resampling:

    image = itk.GetImageFromArray(grid_array)

    interpolator = itk.BSplineInterpolateImageFunction.New(image)
    interpolator.SetSplineOrder(3)

    transform = itk.IdentityTransform.New()
    resampled_image = itk.resample_image_filter(
        image, interpolator=interpolator, transform=transform,
        size=itk.size(image), output_spacing=itk.spacing(image), output_origin=itk.origin(image)
    )
    resampled_grid = itk.GetArrayViewFromImage(resampled_image)

    plt.subplot(132)
    plt.gca().set_title('itk cubic spline')
    plt.imshow(resampled_grid, cmap='gray')

Which also produced the same result as the original image:

zivy · October 14, 2021, 12:48pm

Hello @Theoretialperson,

Well, this time diving into the code I found the specific issue, and it is the choice of interpolator. The interpolator enumerations in SimpleITK, include several that are specifically targeted towards interpolating BSpline coefficient images (sitkBSplineResampler, sitkBSplineResamplerOrder3 etc.) These are not intended for use with a standard image. For a standard image the interpolator type is sitkBSpline and the spline order is three, no choice here.

Will update the documentation to clarify this.

import SimpleITK as sitk
import numpy as np

if __name__ == '__main__':
    image_size = np.array([64, 64])

    grid_array = ((np.indices(image_size) // np.array([4, 4])[:,None,None]).sum(axis=0) % 2).astype(float)

    image = sitk.GetImageFromArray(grid_array)
    resampled_image = sitk.Resample(image, sitk.Transform(2, sitk.sitkIdentity), sitk.sitkBSpline, 1.0)

    diff_arr = sitk.GetArrayFromImage(image-resampled_image)
    print(f' min and max intensity differences: {diff_arr.min()}, {diff_arr.max()}')
    
    sitk.Show(image, 'original')
    sitk.Show(resampled_image, 'resampled')

Theoretialperson · October 14, 2021, 5:22pm

Thanks! I was scratching my head about why my AI was performing worse when using 3rd order spline interpolation instead of linear. After changing the interpolator to sitkBSpline, it fixed the issue of blurriness and the interpolated images are clearly sharper than with just linear interpolation, as it should.

Any reason why the spline order cannot be specified?

zivy · October 14, 2021, 5:35pm

Hello @Theoretialperson,

The limitation on specifying the spline order is part of the “Simple” philosophy, exposing the ITK interface only for the common use cases. 3’rd order BSpline is usually what we end up using anyways when working with BSpline. Linear is often the sweet spot as a compromise between speed and accuracy.

If you still feel the need for more sophisticated interpolation take a look at the windowed-sinc family (listed in the link above).