python - 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)

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.

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: