Compose image from different modality with different number of slices

Algorithms
,
8

Hello @Jigabytes,

Per @mattias advice, if your PET/CT images are not aligned then you will first have to look into registration (see the SimpleITK examples or notebooks for details on that) .

Given that images acquired with modern machines are intrinsically registered you may not need to perform registration.

The main thing to note is that an image occupies a region in physical space, defined by image origin, size (number of pixels), spacing, and direction cosine matrix. ITK/SimpleITK deal with that under the hood. See the SimpleITK fundamental concepts documentation.

The code below does what you want, tested on PET/CT data downloaded from here.

import SimpleITK as sitk

image_dirs = ['PET', 'CT']

def read_first_series_in_dir(dir):
    series_reader = sitk.ImageSeriesReader()
    series_IDs = sitk.ImageSeriesReader.GetGDCMSeriesIDs(dir)
    series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(dir, series_IDs[0])
    series_reader.SetFileNames(series_file_names)    
    return series_reader.Execute()

images = [read_first_series_in_dir(image_dir) for image_dir in image_dirs]

pet_image = images[0]
ct_image = images[1]


print('PET information:')
print('\timage size: {0}'.format(pet_image.GetSize()))
print('\timage spacing: {0}'.format(pet_image.GetSpacing()))
print('\tpixel type: ' + pet_image.GetPixelIDTypeAsString())
print('\tnumber of channels: ' + str(pet_image.GetNumberOfComponentsPerPixel()))

print('CT information:')
print('\timage size: {0}'.format(ct_image.GetSize()))
print('\timage spacing: {0}'.format(ct_image.GetSpacing()))
print('\tpixel type: ' + ct_image.GetPixelIDTypeAsString())
print('\tnumber of channels: ' + str(ct_image.GetNumberOfComponentsPerPixel()))

# Resample PET onto CT grid using default interpolator and identity transformation.
pet_image_resampled = sitk.Resample(pet_image, ct_image)

# Compose the PET and CT image into a single two channel image. 
# The pixel types of all channels need to match, so we upcast the CT from 
# 32bit signed int to the PET pixel type of 64bit float.
pet_ct_combined = sitk.Compose(pet_image_resampled, sitk.Cast(ct_image, pet_image_resampled.GetPixelID()))

print('Combined PET-CT information:')
print('\timage size: {0}'.format(pet_ct_combined.GetSize()))
print('\timage spacing: {0}'.format(pet_ct_combined.GetSpacing()))
print('\tpixel type: ' + pet_ct_combined.GetPixelIDTypeAsString())
print('\tnumber of channels: ' + str(pet_ct_combined.GetNumberOfComponentsPerPixel()))
1 Like