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.
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')
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:
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.