Search code examples
pythonnumpyimagingsimpleitkmedical-imaging

SimpleITK - Coronal/sagittal views problems with size

I'm trying to extract the all three views (Axial, Sagittal and Coronal) from a CTA in DICOM format, using the SimpleItk library.

I can correctly read the series from a given directory:

    ...
    import SimpleITK as sitk
    ...
    reader = sitk.ImageSeriesReader()
    dicom_names = reader.GetGDCMSeriesFileNames(input_dir)
    reader.SetFileNames(dicom_names)
    # Execute the reader
    image = reader.Execute()
    ...

and then, using numpy arrays as stated in this questions, I'm able to extract and save the 3 views.

    ...
    image_array = sitk.GetArrayFromImage(image)
    ...
    for i in range(image_array.shape[0]):
        output_file_name = axial_out_dir + 'axial_' + str(i) + '.png'
        logging.debug('Saving image to ' + output_file_name)
        imageio.imwrite(output_file_name, convert_img(image_array[i, :, :], axial_min, axial_max), format='png')
    ...

The other 2 are made by saving image_array[:, i, :] and image_array[:, :, i], while convert_img(..) is a function that only converts the data type, so it does not alter any shape.

However, the coronal and sagittal views are stretched, rotated and with wide black band (in some slice they are very wide).

Here's the screenshot from Slicer3d:

Sclicer3D screenshot

while this is the output of my code:

Axial

Axial

Sagittal

Sagittal

Coronal

Coronal

Image shape is 512x512x1723, which result in axial pngs being 512x512 pixel, coronal and sagittal being 512x1723, thus this seems correct.

Should I try using PermuteAxes filter? The problem is that I was not able to find any documentation regarding its use in python (neither in other language due to 404 in documentation page)

There is also a way to improve the contrast? I have used the AdaptiveHistogramEqualization filter from simpleitk but it's way worse than Slicer3D visualization, other than being very slow.

Any help is appreciated, thank you!

Solution

  • Answering my own question if someone need it.

    Given the fact that I need to use the slices in a deep learning framework and for data augmentation, I need them to be resampled in a new spacing which is (1.0, 1.0, 1.0).

    Solved it by using this function:

    def resample_image(itk_image, out_spacing=(1.0, 1.0, 1.0)):
    """
    Resample itk_image to new out_spacing
    :param itk_image: the input image
    :param out_spacing: the desired spacing
    :return: the resampled image
    """
    # get original spacing and size
    original_spacing = itk_image.GetSpacing()
    original_size = itk_image.GetSize()
    # calculate new size
    out_size = [
        int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0]))),
        int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1]))),
        int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])))
    ]
    # instantiate resample filter with properties and execute it
    resample = sitk.ResampleImageFilter()
    resample.SetOutputSpacing(out_spacing)
    resample.SetSize(out_size)
    resample.SetOutputDirection(itk_image.GetDirection())
    resample.SetOutputOrigin(itk_image.GetOrigin())
    resample.SetTransform(sitk.Transform())
    resample.SetDefaultPixelValue(itk_image.GetPixelIDValue())
    resample.SetInterpolator(sitk.sitkNearestNeighbor)
    return resample.Execute(itk_image)

    and then saving by using numpy arrays as stated in the original question.