Compose image from different modality with different number of slices

Algorithms
,
1

Hi,

I want to combine images from different modality into one composed image and then perform segmentation on it. My q is does itk have provision to compose two n dimensional images where number of slices in both differ? For example I want to merge PET scan and CT Scan image but CT scan has ~800 dicom images whereas PET scan has ~300 dicom images

Thanks,
Jiten

2

Hello @Jigabytes,

This can be done in two steps:

  1. Resample PET onto the CT image grid.
  2. Compose the resampled PET with the CT to get a combined two channel image.

Hopefully this is what you wanted (not 100% sure I understood your request).

1 Like
3

Thanks much Zivy! Yeah I think you are right on the money here ! I just need to understand overall what is resampling concept but yeah I need to get PET image synched up into CT and then compose them once they are of same dimensions, types etc…

I will give it a try and will let you know the feedback.

Thanks,
Jiten

4

Hi Zivy,

Any information from CT scan that I need to pass while resampling pet scan? Like which transform, interpolator I should be using?

Thanks,
Jiten

5

Hi Zivy,

I was able to apply resample and compose on single slice of pet and ct scan each. Just to check the results. For resampling of pet image I just used pet image and ct image as argument, I did not give any transform or interpolator values, I observed that the centre of pet image gets shifted. I am not sure which transform should I apply here Any suggestions?

Thanks,
Jiten

6

Here is what I did
Pet image was 192 X 192
CT image was 512 x 512
so after doing simple math I came up with following translation transform
point = (1.0, 1.0)
identity_transform = sitk.Transform(2, sitk.sitkIdentity)
translation = sitk.TranslationTransform(2)
translation.SetParameters((160.0, 160.0))
transform_point(translation, point)

After applying above inverse of above translation transform I was able to see resampled pet image as before without resampling.

Let me know if I am doing thing in right direction?

Thanks,
Jiten

7

Since you have a CT and PET image from different frame of references I would look into a full registration first. This might be a bit tricky depending on what the PET image looks like but it could give you a good result.

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
12

Hello @mattias and @zivy thanks much for your responses.

I tried your code on images I have Below is more details on output of my program on images that I have. Do you still think I need to perform image registration?
As per my knowledge below scanner is the modern one and registration won’t be required.

I will try to plot these images to make sure I see what I expect from the composition. In addition I was using pydicom for reading images but itk for segmentation. I changed reading of dicom file to read from ITK and now results are much better!

Thanks,
Jiten

Details of output on my pet ct images.

PET information:
image size: (192, 192, 1)
image spacing: (3.6458332538605, 3.6458332538605, 1.0)
pixel type: 64-bit float
number of channels: 1
CT information:
image size: (512, 512, 1)
image spacing: (0.976562, 0.976562, 1.0)
pixel type: 32-bit signed integer
number of channels: 1
Resampled PETinformation:
image size: (512, 512, 1)
image spacing: (0.976562, 0.976562, 1.0)
pixel type: 64-bit float
number of channels: 1
Combined PET-CT information:
image size: (512, 512, 1)
image spacing: (0.976562, 0.976562, 1.0)
pixel type: vector of 64-bit float
number of channels: 2

13

If your images are taken in a combined PET/CT machine and are from the same session you will not need a registration since they are taken in the same frame of reference (usually), you can check this by comparing the “Frame of Reference UID” (0020, 0052) dicom tag of the two. If they are the same all you need to do is to re-sample the PET to the CT like @zivy suggested.

If you want to complicate things you might still need some registration/processing to compensate for things like patient movement during the session.

1 Like
14

Hi @mattias Thanks for the tag. I confirmed and they were same so no registration required! I was able to compose PET CT successfully. Thanks a ton to you and @zivy for assistance here. On to next problem :slight_smile:

15

Hi @zivy, while doing resampling of images do we always have to resample smaller dimension (192 x 192 ) image to a bigger dimension (512 x 512) image or we can resample bigger dimension (512 x 512) image to smaller dimension (192 x 192 ) as well? Thanks Jiten

16

Hello @Jigabytes,

No, resampling to a smaller sized grid will work just fine.

Please note that if your original image has high frequencies you may encounter aliasing artifacts in the smaller image because you have fewer samples. The common way to deal with this is to smooth (low pass filter) the original image before resampling.

17

Thanks @zivy for the info. So does Curvature flow filter qualify as low pass filter from filters available in simpleITK? How to find what all low pass filter are available in ITK APIs to smooth the original image before resampling. In my use case I am using Curvature flow API to first get smooth images and then resampling it. Thanks, Jiten

19

pet_image_resampled = sitk.Resample(pet_image, ct_image)
powerful function, very good, had some doubt, but was wrong, tested and deleted above post, sorry

1 Like