Slice 3D volume along a specific axis using python ITK

I found Read DICOM Series and Write 3D Image — v5.3.0 and it works, but I can only get slices from one axis (along z in my case). Is there a way to specify which axis I’m slicing along (I need x and y as well).

I saw a thread (Extract an oblique 2D slice from a 3D volume) that suggested that I can use resample image filter, but I couldn’t find any examples for it.

ps: I used to do it by reading the DICOM files using SimpleITK, converting to VTK, and then doing the slicing in VTK. But now I need to use ITK-specific functions on the sliced images, so it would be better (at least, I think) if the only thing that VTK was doing is to display the image, with everything else being done by ITK.

1 Like

You could take a look at this thread:

1 Like

With SimpleITK you can extract the slices like this:

import SimpleITK as sitk

vol = sitk.Image(100, 100, 100, sitk.sitkUInt8)

# extract the 42nd X, Y and Z slices
xslice = vol[42, :, :]
yslice = vol[:, 42, :]
zslice = vol[:, :, 42]

You can also try FREDtools with uses SimpleITK. An example code:

import fredtools as ft
import matplotlib.pyplot as plt

# read 3D image from mhd file (or create/read any other format as a SimpleITK object)
img=ft.readMHD('file_name.mhd', displayInfo=True)

# Get 2D SimpleITK from 3D SimpleITK
# Any slice, for instance 'XY', 'ZX', or even something like 'X-Z' can be calculated
# The slice will be interpolated from 3D, including 'linear' or 'nearest' interpolation methods
# The 'point' describes (x,y,z) coordinates in the image units (e.g. mm) to get the slice through

slice=ft.getSlice(img, point=[0,4,5], plane='XY', interpolation='linear')

# display the 2D slice or use ft.showSlice(...) method
plt.imshow(ft.arr(slice), extent=ft.getExtMpl(slice), cmap='jet')

check https://fredtools.ifj.edu.pl/ manual and tutorials on github for more details and find more functionalities of the FREDtools.

Disclaimer: I am one of the authors.

Thanks, but I specifically need to use ITK to do the slicing, not SimpleITK, because I need to use the sliced images in ITK Elastix. As I understand, the only way to go from SimpleITK to ITK is by converting it to a numpy array, which is not ideal.

As far as I know, the SimpleITK and ITK image objects are very similar. FREDtools provides two routines:

fredtools.SITK2ITK(imgSITK)
and
fredtools.ITK2SITK(imgITK)

to convert between them.

The approach proposed by @dchen is also very useful but allows getting a slice only through the centre of a voxel (equivalent to the nearest interpolation method). Moreover, you have to calculate the voxel from the real-world coordinates.

You could use the PermuteAxesImageFilter to create different versions of your volume. Map the X axis to Z, and then slice that volume to get X slices, and similarly for Y to Z.

Maybe apply this suggestion:

1 Like