Squeeze dicom ITK image from 4D to 3D

1

I have Dicom images in an unusual format. For example, usually, Dicom images are saved in a casepath/folder in multiple slices.
IM-0001-0240-0001.dcm
I use the following helper function to read the dicom images:

def dicomread(casePath, imgitk=None):
    """
    casePath: the folder location where all the dicom slices are in a subject/case
              e.g. r'E:\CAP Exam Data\Images\SPH-0706\66154284\20627700'
    return: 3D Image in numpy format and all image informations
    """
    if not os.path.isdir(casePath):
        print("The path doesn't exist.")
    reader = sitk.ImageSeriesReader()
    dicom_names = reader.GetGDCMSeriesFileNames(casePath)
    reader.SetFileNames(dicom_names)
    image = reader.Execute()
    image_array = sitk.GetArrayFromImage(image) # z, y, x
    origin = image.GetOrigin() # x, y, z
    spacing = image.GetSpacing() # x, y, z
    size = image.GetSize()
    direction = image.GetDirection() # image
    data_type = image.GetPixelIDTypeAsString()
    if imgitk == True:
        return image
    else:
        return image_array, origin, spacing, size, direction, data_type

For one image like I mentioned I get 4 dimensions instead of 3.

>> dcmItk = dicomread(r'/media/banikr/Elucid_data/dicoms', True)
>> dcmItk.GetSize()
(512, 512, 480, 1)

Is it possible to squeeze the itk image and use the 3D output in all itk operations?
For example, I use resample to align masks on 3D ITK images.

2

Hello @banikr,
This should do the trick inside your function (using SimpleITK>=2.1.0):

image = reader.Execute()
if image.GetDimension()==4 and image.GetSize()[3]==1:
  image = image[...,0]
1 Like
3

This actually worked. Thanks for the hint.

def dicomread(casePath, imgitk=None):
    """
    casePath: the folder location where all the dicom slices are in a subject/case
              e.g. r'E:\Elucid_data\CAP Exam Data\Images\SPH-0706\66154284\20627700'
    return: 3D Image in numpy format and all image informations
    """
    if not os.path.isdir(casePath):
        print("The path doesn't exist.")
    reader = sitk.ImageSeriesReader()
    dicom_names = reader.GetGDCMSeriesFileNames(casePath)
    reader.SetFileNames(dicom_names)
    image = reader.Execute() # type --> 'SimpleITK.SimpleITK.Image'
    if image.GetDimension() == 4 and image.GetSize()[3] == 1:
        image = sitk.GetArrayFromImage(image)[0,...]
        image = sitk.GetImageFromArray(image) # type --> 'SimpleITK.SimpleITK.Image'
    image_array = sitk.GetArrayFromImage(image) # z, y, x
    origin = image.GetOrigin() # x, y, z
    spacing = image.GetSpacing() # x, y, z
    size = image.GetSize()
    direction = image.GetDirection() # image
    data_type = image.GetPixelIDTypeAsString()
    if imgitk == True:
        return image
    else:
        return image_array, origin, spacing, size, direction, data_type
4

Hello @banikr,

Your new code has a bug, the image lost its spatial information (spacing/origin/direction matrix) in the transition to numpy and back to SimpleITK. If you don’t want to upgrade to SimpleITK 2.1.0 for the new syntax you can use the old one, both retain the spatial information:

image = reader.Execute()
if image.GetDimension()==4 and image.GetSize()[3]==1:
  image = image[:,:,:,0]
2 Likes
5

I didn’t notice that! Thanks for the catch. I think removing the sitk.Get.. modules works without any issues.

def dicomread(casePath, imgitk=None):
    """
    casePath: the folder location where all the dicom slices are in a subject/case
              e.g. r'E:\Elucid_data\CAP Exam Data\Images\SPH-0706\66154284\20627700'
              also works for casepath will only one .dcm extension with 4D size.
    return: 3D Image in numpy format and all image informations
    """
    if not os.path.isdir(casePath):
        print("The path doesn't exist.")
    reader = sitk.ImageSeriesReader()
    dicom_names = reader.GetGDCMSeriesFileNames(casePath)
    reader.SetFileNames(dicom_names)
    image = reader.Execute() # type --> 'SimpleITK.SimpleITK.Image'
    if image.GetDimension() == 4 and image.GetSize()[3] == 1:
#         image = sitk.GetArrayFromImage(image)[0,...]
        image = image[:,:,:,0]
#         image = sitk.GetImageFromArray(image) # type --> 'SimpleITK.SimpleITK.Image'
    if imgitk == True:
        return image
    else:
        image_array = sitk.GetArrayFromImage(image) # z, y, x
        origin = image.GetOrigin() # x, y, z
        spacing = image.GetSpacing() # x, y, z
        size = image.GetSize()
        direction = image.GetDirection() # image
        data_type = image.GetPixelIDTypeAsString()
        return image_array, origin, spacing, size, direction, data_type

The function results:

>> new_img = dicomread(dcmPath, True)
>> new_img.GetSize()

(512, 512, 240)

One other minor thing I noticed:
image = image[:,:,:,0] indexing works but image = image[...,0] generates IndexError: invalid index. Maybe due to data structure in ITK vs numpy or so.