CodingAddict
CodingAddict

Reputation: 223

Transpose 3-D MRI image axis to different image view

How can I transpose 3D MRI image with SimpleITK to new image view, but keep aspect ration as it was in the initial view? Here is my code:

def show_n_slices(array, start_from=0, step=10, columns=6, figsize=(18,10)):
    """ Plot N slices of a 3D image.
    :param array: N-dimensional numpy array (i.e. 3D image)
    :param start_from: slice index to start from
    :param step: step to take when moving to the next slice
    :param columns: number of columns in the plot
    :param figsize: figure size in inches
    """
    array = np.swapaxes(array, 0, 2)
    fig, ax = plt.subplots(1, columns, figsize=figsize)
    slice_num = start_from
    for n in range(columns):
        ax[n].imshow(array[:, :, slice_num], 'gray')
        ax[n].set_xticks([])
        ax[n].set_yticks([])
        ax[n].set_title('Slice number: {}'.format(slice_num), color='r')
        slice_num += step

    fig.subplots_adjust(wspace=0, hspace=0)
    plt.show()

def print_sitk_info(itk_image):
    """ Prints SimpleITK image information
    :param itk_image: SimpleITK image object.
    """
    print(f"[INFO]: Shape - {itk_image.GetSize()}")
    print(f"[INFO]: Spacing - {itk_image.GetSpacing()}")
    print(f"[INFO]: Origin - {itk_image.GetOrigin()}")
    print(f"[INFO]: Direction - {itk_image.GetDirection()}\n")

print("[INFO]: Image Info before Resampling to Isotropic Resolution:")
print_sitk_info(itk_image)
show_n_slices(sitk.GetArrayFromImage(itk_image), start_from=20, step=5)

enter image description here

And here is the code where I try to transpose image axes in order to obtain different image view:

array2 = array.transpose(2, 0 ,1)
itk_image2 = sitk.GetImageFromArray(array2)
itk_image2.SetOrigin(itk_image.GetOrigin())
itk_image2.SetSpacing(itk_image.GetSpacing())
itk_image2.SetDirection(itk_image.GetDirection())
array2 = sitk.GetArrayFromImage(itk_image2)

print_sitk_info(itk_image2)
show_n_slices(array2, start_from=30, step=5)

enter image description here

What am I doing wrong here? Do I have to resample image size next?

EDIT [Based on Dave Chens answer]:

Result of adding spacing in the correct order:

array2 = array.transpose(2, 0 ,1)
spacing = itk_image.GetSpacing()
itk_image2 = sitk.GetImageFromArray(array2)
itk_image2.SetOrigin(itk_image.GetOrigin())
itk_image2.SetSpacing([spacing[2], spacing[0], spacing[1]])
itk_image2.SetDirection(itk_image.GetDirection())
array2 = sitk.GetArrayFromImage(itk_image2)

print_sitk_info(itk_image2)
show_n_slices(array2, start_from=30, step=5)

enter image description here

Result of SimpleITK PermuteAxesImageFilter:

pa = sitk.PermuteAxesImageFilter()
pa.SetOrder([2, 0 , 1])
i = pa.Execute(itk_image)
print_sitk_info(i)
show_n_slices(sitk.GetArrayFromImage(i), start_from=30, step=5)

enter image description here

Upvotes: 1

Views: 1621

Answers (1)

Dave Chen
Dave Chen

Reputation: 2085

The issue is that you're using the pixel spacing from the original image. You need to transpose the pixel spacings the same way you transpose the pixels. So you'd do something like this:

spacing = itk_image.GetSpacing()
spacing2 = [spacing[2], spacing[0], spacing[1]]
itk_image2.SetSpacing(spacing2)

Also, note that SimpleITK has a PermuteAxes function which can do the same thing and will preserve the image meta-data, so you wouldn't have to copy it over.

UPDATE: Ok, below is the code that I used. I run it in a Jupyter Notebook and use ITKWidgets to render the 3d image in the notebook. Note that each of the view calls at the bottom should be in its own notebook cell.

import SimpleITK as sitk
import itkwidgets

fnames = []
for i in range(1,33):
    name = "Image-"+str(i)+".dcm"
    fnames.append(name)
reader = sitk.ImageSeriesReader()
reader.SetFileNames(fnames)
img = reader.Execute()
print(img)

pa = sitk.PermuteAxesImageFilter()
pa.SetOrder([2,0,1])
img2 = pa.Execute(img)
print(img2)

# Render cell #1
itkwidgets.view(img)

# Render cell #2
itkwidgets.view(img2)

Upvotes: 2

Related Questions