SimpleITK Resize images

Question

I have a set o 3D volumes that I am reading with SimpleITK

import SimpleITK as sitk
for filename in filenames:
    image = sitk.ReadImage(filename)

Each of the volumes has different size, spacing, origin and direction. This code yields different values for different images:

print(image.GetSize())
print(image.GetOrigin())
print(image.GetSpacing())
print(image.GetDirection())

My question is: how do I transform the images to have the same size and spacing so that they all have the same resolution and size when converted to numpy arrays. Something like:

import SimpleITK as sitk
for filename in filenames:
    image = sitk.ReadImage(filename)
    image = transform(image, fixed_size, fixed_spacing)
    array = sitk.GetArrayFromImage(image)

Wicklund

Answer 1

The way to do this is to use the Resample function with fixed/arbitrary size and spacing. Below is a code snippet showing construction of this "reference_image" space:

reference_origin = np.zeros(dimension)
reference_direction = np.identity(dimension).flatten()
reference_size = [128]*dimension # Arbitrary sizes, smallest size that yields desired results. 
reference_spacing = [ phys_sz/(sz-1) for sz,phys_sz in zip(reference_size, reference_physical_size) ]

reference_image = sitk.Image(reference_size, data[0].GetPixelIDValue())
reference_image.SetOrigin(reference_origin)
reference_image.SetSpacing(reference_spacing)
reference_image.SetDirection(reference_direction)

For a turnkey solution have a look at this Jupyter notebook which illustrates how to do data augmentation with variable sized images in SimpleITK (code above is from the notebook). You may find the other notebooks from the SimpleITK notebook repository of use too.

Lovesick

What transformations should be used in the Resample() function in the case of simply wanting to resize? Also could please comment on the interpolator and default value parameters? Thank you – Wicklund

I created a github gist which just shows how to do the resampling. For most situations the linear interpolator is good enough, but for label interpolation you need the nearest neighbor. The default_intensity_value is the value to set if mapped outside the image bounds. – Lovesick

sorry to bother you again. When resampling to a given resolution and size and then back there are bound to be losses. With regular images I found these losses to be admissible but with labeled images I some time get losses of 20% of the segmented region (even when using nearest neighbour interpolation). Any suggestion on how to address this issue? – Wicklund

The only reason for using nearest neighbor is that it doesn't introduce new labels but it will make aliasing artifacts worse. Possibly treat the segmentation image as a real valued image instead of binary and use any other interpolator followed by thresholding to get a binary image back. Also, smooth the image when going to coarser resolution. This is just mitigating the issue, not solving it. For future discussions please use discourse.itk.org (simpleitk tag), the stackoverflow character limit on comments is restrictive. – Lovesick

Answer 2

According to SimpleITK's documentation, the process of image resampling involves 4 steps:

1. Image - the image we resample, given in the coordinate system;
2. Resampling grid - a regular grid of points given in a coordinate system which will be mapped to the coordinate system;
3. Transformation - maps points from the coordinate system to coordinate system;
4. Interpolator - a method for obtaining the intensity values at arbitrary points in the coordinate system from the values of the points defined by the Image

The following snippet is for downsampling the image preserving its coordinate system properties:

def downsamplePatient(patient_CT, resize_factor):

    original_CT = sitk.ReadImage(patient_CT,sitk.sitkInt32)
    dimension = original_CT.GetDimension()
    reference_physical_size = np.zeros(original_CT.GetDimension())
    reference_physical_size[:] = [(sz-1)*spc if sz*spc>mx  else mx for sz,spc,mx in zip(original_CT.GetSize(), original_CT.GetSpacing(), reference_physical_size)]
    
    reference_origin = original_CT.GetOrigin()
    reference_direction = original_CT.GetDirection()

    reference_size = [round(sz/resize_factor) for sz in original_CT.GetSize()] 
    reference_spacing = [ phys_sz/(sz-1) for sz,phys_sz in zip(reference_size, reference_physical_size) ]

    reference_image = sitk.Image(reference_size, original_CT.GetPixelIDValue())
    reference_image.SetOrigin(reference_origin)
    reference_image.SetSpacing(reference_spacing)
    reference_image.SetDirection(reference_direction)

    reference_center = np.array(reference_image.TransformContinuousIndexToPhysicalPoint(np.array(reference_image.GetSize())/2.0))
    
    transform = sitk.AffineTransform(dimension)
    transform.SetMatrix(original_CT.GetDirection())

    transform.SetTranslation(np.array(original_CT.GetOrigin()) - reference_origin)
  
    centering_transform = sitk.TranslationTransform(dimension)
    img_center = np.array(original_CT.TransformContinuousIndexToPhysicalPoint(np.array(original_CT.GetSize())/2.0))
    centering_transform.SetOffset(np.array(transform.GetInverse().TransformPoint(img_center) - reference_center))
    centered_transform = sitk.Transform(transform)
    centered_transform.AddTransform(centering_transform)

    # sitk.Show(sitk.Resample(original_CT, reference_image, centered_transform, sitk.sitkLinear, 0.0))
    
    return sitk.Resample(original_CT, reference_image, centered_transform, sitk.sitkLinear, 0.0)

Using the snippet above in a brain CT scan we get:

Svend

centered_transform = sitk.Transform(transform) should be: centered_transform = sitk.CompositeTransform(transform) – Discompose