Reputation: 223
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)
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)
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)
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)
Upvotes: 1
Views: 1621
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