Resampling to volume to smaller size and smaller voxel spacing

Beginner Questions
1

Hi ITK community,

I have recently started learning SimpleITK and have been stuck on a specific type of Resampling that I have been trying to implement. I have PET/CT Nifti images. My typical PET and CT volumes have the following characteristics:

PET Size: (168, 168, 324)
PET Spacing: (4.07283, 4.07283, 3.0)

CT Size: (512, 512, 486)
CT Spacing: (0.9765625, 0.9765625, 2.0)

I want to downsample my CT images to match their corresponding PET Size, but I also want all my voxel spacing for both PET and CT to be isotropic and equal to (1.0, 1.0, 1.0). So far, I have tried codes by people who have already asked similar questions, but most of them ended up cropping my final image. Is it possible to perform this resampling?

I am attaching the latest version of my code where I am trying to resample a PET image to a volume where I keep the size the same as the original PET image, but change the voxel spacing to (1.0, 1.0, 1.0). This code results in getting the final resampled image cropped. Please let me know what I am doing wrong.

Thanks in advance for your help. :slight_smile:

def reSample_fixedSize_fixedSpacing (img_sitk,interpolation, new_spacing):
   
    dimension = img_sitk.GetDimension()
    reference_size = img_sitk.GetSize()#[new_size,new_size,new_size]
    reference_physical_size = np.zeros(img_sitk.GetDimension())
    reference_physical_size[:] = [(sz-1)*spc if sz*spc>mx  else mx for sz,spc,mx in zip(img_sitk.GetSize(), img_sitk.GetSpacing(), reference_physical_size)]
    reference_direction = img_sitk.GetDirection()
    reference_origin = img_sitk.GetOrigin()
    
    reference_image = sitk.Image(reference_size, img_sitk.GetPixelIDValue())
    reference_image.SetOrigin(reference_origin)
    reference_image.SetSpacing(new_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(img_sitk.GetDirection())
    transform.SetTranslation(np.array(img_sitk.GetOrigin()) - reference_origin)
  
    centering_transform = sitk.TranslationTransform(dimension)
    img_center = np.array(img_sitk.TransformContinuousIndexToPhysicalPoint(np.array(img_sitk.GetSize())/2.0))
    centering_transform.SetOffset(np.array(transform.GetInverse().TransformPoint(img_center) - reference_center))
    centered_transform = sitk.CompositeTransform([centering_transform,transform])
    new_img = sitk.Resample(img_sitk, reference_image, centered_transform, interpolation, 0.0)
    return new_img


ptpath = '14695318_20170802_pt.nii.gz'
ptimg = sitk.ReadImage(ptpath, imageIO="NiftiImageIO")

resamp_ptimg = reSample_fixedSize_fixedSpacing(ptimg, sitk.sitkBSpline, [1.0, 1.0, 1.0])

info = os.path.basename(ptpath).split('_')
filename = info[0] + '_' + info[1]
sitk.WriteImage(resamp_ptimg,filename+'_pt_resamp.nii')
2

Hello @shadab,

Welcome to SimpleITK!

What you are missing is the understanding that in ITK/SimpleITK an image is a physical/spatial object with a metric size, see Fundamental Concepts. In ITK/SimpleITK nomenclature the term Size means number of pixels/voxels and the term Spacing refers to the distance between consecutive indexes and is a metric unit (usually mm). Defining a resampling grid is done in physical space.

“I want to downsample my CT images to match their corresponding PET Size, but I also want all my voxel spacing for both PET and CT to be isotropic and equal to (1.0, 1.0, 1.0).”

If you want the PET image to be isotropic (1.0,1.0,1.0), while retaining the original number of voxels, you are defining a grid that is smaller than the original PET image because the new spacing is smaller than the original one and the Size is the same as the original:

pet_original_size = (168, 168, 324)
pet_original_spacing = (4.07283, 4.07283, 3.0)
pet_original_physical_size = [(sz-1)*spc for sz,spc in zip(pet_original_size,pet_original_spacing)]

pet_new_spacing = [1.0,1.0,1.0]
pet_physical_size_using_new_spacing_original_size = [(sz-1)*spc for sz,spc in zip(pet_original_size, pet_new_spacing)]

print(pet_original_physical_size)
print(pet_physical_size_using_new_spacing_original_size)

If you want both images to have the same spacing and size you need to:

  1. Resample PET using given new_spacing and compute the required new_size to maintain the original physical size.
  2. Resample CT using new_spacing and the PET’s new_size.
2 Likes
3

Hi @zivy

Thank you so much for your detailed explanation. I understand the difference between Size and “Physical size” of the image in the world of SimpleITK. But, my aim is to do exactly what you mentioned in your reply. I want to change the physical size of the image, but I want to keep my Spacing and Size to whatever I specify them to be. So, that would mean downsampling CT to corresponding PET sizes while retaining the number of voxels for PET after resampling and setting the spacing for both PET and CT to (1.0, 1.0, 1.0) after resampling.

I was just wondering if this is possible, as no one really does this (as far as I have noticed)!

4

Hello @shadab ,

Sorry, but not really sure what it is you are attempting to do. Based on your requirements you are intentionally cropping the PET image, so not sure why you are surprised that the image is cropped. Be aware that this is not a symmetric cropping around the center of the image.

Possibly if you provide the expected input and output we could better help. Is the following what you are thinking?

Input:

PET Size: (168, 168, 324)
PET Spacing: (4.07283, 4.07283, 3.0)
PET Origin: ?

CT Size: (512, 512, 486)
CT Spacing: (0.9765625, 0.9765625, 2.0)
CT Origin: ?

Output:

PET Size: (168, 168, 324)
PET Spacing: (1.0, 1.0, 1.0)
PET Origin: ?

CT Size: (168, 168, 324)
CT Spacing: (1.0, 1.0, 1.0)
CT Origin: ?
5

HI @zivy

I want the whole 3D image to be contained (maybe “resized” is the correct terminology here?) in a smaller physical volume after resampling while having the final Spacing and final Size to whatever I specify.

As for the inputs and outputs, you are correct.

Input:

PET Size: (168, 168, 324)
PET Spacing: (4.07283, 4.07283, 3.0)
PET Origin: (-342.45751953125, -491.60284423828125, -1466.5)
PET Direction: (1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)

CT Size: (512, 512, 486)
CT Spacing: (0.9765625, 0.9765625, 2.0)
CT Origin: (-249.51171875, -399.51171875, -1467.5)
CT Direction: (1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)

Output:

PET Size: (168, 168, 324)
PET Spacing: (1.0, 1.0, 1.0)
PET Origin: ?
PET Direction: ?

CT Size: (168, 168, 324)
CT Spacing: (1.0, 1.0, 1.0)
CT Origin: ?
CT Direction: ?

For the Output Origin and Direction for PET and CT, as long as the final resampled CT image is well-aligned/registered to the final resampled PET image, it will be great.

6

Hello @shadab,

I suspect the code below is what you want. It means that you are willing to accept that the objects will change their physical size and that the change is non-isotropic.

import SimpleITK as sitk

def scale_resample(image, interpolator, new_size, new_spacing):
    original_physical_size = [(sz-1)*spc for sz,spc in zip(image.GetSize(), image.GetSpacing())]
    new_physical_size = [(sz-1)*spc for sz,spc in zip(new_size, new_spacing)]

    tx = sitk.ScaleTransform(image.GetDimension(), [ops/nps for nps, ops in zip(new_physical_size, original_physical_size)])
    return sitk.Resample(image,
                         size = new_size,
                         transform = tx,
                         interpolator = interpolator,
                         outputSpacing = new_spacing)

file_name =   # set the file name
new_size = (168, 168, 324)
new_spacing = (1.0, 1.0, 1.0)
image = sitk.ReadImage(file_name)

sitk.Show(scale_resample(image, sitk.sitkLinear, new_size, new_spacing))
7

Hi @zivy

Thanks again for your reply. I used the code you gave above and now it returns an empty 3D array (all elements are zeros).

The sitk.Show() command didn’t really work for me, and I got the following error:

RuntimeError: Exception thrown in SimpleITK Show: D:\a\1\sitk\Code\IO\src\sitkImageViewer.cxx:620:
sitk::ERROR: No ImageJ/Fiji application found.

I just ended up getting the array from the generated image and used this:

def scale_resample(image, interpolator, new_size, new_spacing):
    original_physical_size = [(sz-1)*spc for sz,spc in zip(image.GetSize(), image.GetSpacing())]
    new_physical_size = [(sz-1)*spc for sz,spc in zip(new_size, new_spacing)]

    tx = sitk.ScaleTransform(image.GetDimension(), [ops/nps for nps, ops in zip(new_physical_size, original_physical_size)])
    return sitk.Resample(image,
                         size = new_size,
                         transform = tx,
                         interpolator = interpolator,
                         outputSpacing = new_spacing)


#%%
ptpath = '14695318_20170802_pt.nii.gz'
ptimg = sitk.ReadImage(ptpath)
new_size = ptimg.GetSize()
new_spacing = [1.0, 1.0, 1.0]

#%%
resampled_ptimg = scale_resample(ptimg, sitk.sitkLinear, new_size, new_spacing)
resampled_ptarray = sitk.GetArrayFromImage(resampled_ptimg)

After that, I tried to print any element that is non-zero but didn’t get any output:

flatarray = resampled_ptarray.flatten()
for i in flatarray:
    if i != 0:
        print(i)
8

Hello @shadab,

Results are surprising, can you share a sample image where the code doesn’t work?

With respect to the Show function see the FAQ.

10

Hello @shadab,

Thanks for sharing the data. The problem with the code above was that I assumed the image origin is at [0,0,0] and direction cosine matrix was the identity. The code below works for the general setting. Hopefully this addresses your question.

def scale_resample(image, interpolator, new_size, new_spacing):
    original_physical_size = [(sz-1)*spc for sz,spc in zip(image.GetSize(), image.GetSpacing())]
    new_physical_size = [(sz-1)*spc for sz,spc in zip(new_size, new_spacing)]

    scale_tx = sitk.ScaleTransform(image.GetDimension(), [ops/nps for nps, ops in zip(new_physical_size, original_physical_size)])
    rigid_tx = sitk.Euler3DTransform()
    rigid_tx.SetMatrix(image.GetDirection())
    rigid_tx.SetTranslation(image.GetOrigin())
    #accomodate for non identity direction cosine matrix, ensure that scaling is
    #along the image axes and not the canonical ones.
    tx = sitk.CompositeTransform([rigid_tx, scale_tx, rigid_tx.GetInverse()])

    #output image has the same origin and direction as the input, but is a
    #scaled version of the content
    return sitk.Resample(image,
                         size = new_size,
                         transform = tx,
                         interpolator = interpolator,
                         outputOrigin = image.GetOrigin(),
                         outputDirection = image.GetDirection(),
                         outputSpacing = new_spacing)
11

Hi @zivy

I realized something about my image only now, and now my goals have changed. Now I want to keep my PET images as it is, and only resample CT images to PET size and spacing. But the origin of PET and their corresponding CT images are not the same, and the Direction cosines may also need not be the same. The physical size of the CT image is also slightly smaller than the corresponding PET image. But in the current image pairs, the PET and CT images are well aligned/registered. I am using the following code to resample CT to PET size and spacing now:

import SimpleITK as sitk

def resample_ct2pt(ctimg, ptimg, interpolator):#, interpolator, new_size, new_spacing):
    resampler = sitk.ResampleImageFilter()
    resampler.SetOutputOrigin(ctimg.GetOrigin())
    resampler.SetOutputDirection(ctimg.GetDirection())
    resampler.SetSize(ptimg.GetSize())
    resampler.SetOutputSpacing(ptimg.GetSpacing())
    resampler.SetInterpolator(interpolator)
    resampled_ctimg = resampler.Execute(ctimg)
    
    return resampled_ctimg 

ctpath = # ct image here
ptpath = # pet image here
info = os.path.basename(ptpath).split('_')
filename = info[0] + '_' + info[1]

ctimg = sitk.ReadImage(ctpath, imageIO="NiftiImageIO")
ptimg = sitk.ReadImage(ptpath, imageIO="NiftiImageIO")

resampled_ctimg = resample_ct2pt(ctimg, ptimg, sitk.sitkLinear)
sitk.WriteImage(resampled_ctimg,filename+'_ct_new.nii')

I have the following two screenshots: The top image is the PET/CT overlapped in the Lifex software before resampling, and the bottom image is after resampling. I have a feeling that they might not be as well aligned as before (because of the new grey region that appears on the right side of the image in the bottom images). How can I ensure they are aligned in the same way after resampling as they were before?

Thank you for your help so far.

12

Hello @shadab,

  1. The grey region is due to the defaultPixelValue, which has a default value of 0. When working with CT a more appropriate default value is -1000 (HU value of air).
  2. A quick way to check spatial correspondence between images is to use ITK-SNAP with its linked cursor. Open two instances of the software with the two images, place the cursor in one of them on an easily identifiable point. If the images are aligned, the cursor in the second instance will move to the corresponding point. This is similiar to the linked cursor approach used in this notebook.
  3. When resampling a high resolution CT onto a low resolution grid (PET is usually much lower resolution) you will likely need to blur the CT before resampling, for details see this notebook section titled Before we start, a word of caution.

Finally, no need for a resample_ct2pt function as this is a one liner:

resampled_ctimg = sitk.Resample(sitk.SmoothingRecursiveGaussian(ctimg, sigma), ptimg, defaultPixelValue=-1000)
2 Likes