RegionOfInterestImageFilter to extract image slice from volume

Hello,

I’d like to extract a 2D slice from a MRI NIfTI file. I believe the RegionOfInterestImageFilter function in itk will do the job:

https://itk.org/ITKExamples/src/Filtering/ImageGrid/ExtractRegionOfInterestInOneImage/Documentation.html

However, I’m having difficulties using this function in Python. I’ve tried writing the code this way, without success:

image_slice11 = itk.region_of_interest_filter(image, [512,512,1], [0,0,11])

Can you provide advice on how that should be implemented?

Best regards!

It feels like you are confusing Region of Interest with ExtractImageFilter. What happens if you replace RoI by Extract in your code?

1 Like

Many thanks for the suggestion, @dzenanz!

When I replace the function in the line of code, I get the error below.

image_slice11 = itk.extract_image_filter(fixed, [512,512,1], [0,0,11])

TypeError: Expecting argument of type itkImageF3 or itkImageSourceIF3.
Additional information:
Wrong number or type of arguments for overloaded function 'itkImageToImageFilterIF3IF3_SetInput'.
  Possible C/C++ prototypes are:
    itkImageToImageFilterIF3IF3::SetInput(itkImageF3 const *)
    itkImageToImageFilterIF3IF3::SetInput(unsigned int,itkImageF3 const *)

I don’t understand what’s the problem, since I’ve checked and the input image ‘fixed’ is of type ‘itkImageF3’. Could it be that I’m setting the other inputs wrong?

I appreciate your help. Best!

What you seem to want is a subset of this example. The example uses the more verbose functional syntax, but you could try updating it to procedural syntax.

1 Like

Thank you, for sharing this example!

Another option to be aware of: a) use itk.xarray_from_image, b) do the slicing with pandas/numpy syntax, c) then itk.image_from_xarray. This may be more convenient while preserving spatial metadata.

1 Like

Thank you @matt.mccormick, for suggesting this convenient alternative!

Dear @matt.mccormick,

I’ve implemented your suggestion to extract a single 2D slice, using the following lines of code:

fixed_array = itk.xarray_from_image(fixed)

fixed_array_slice = fixed_array[11,:,:]

fixed_slice = itk.image_from_xarray(fixed_array_slice)

Unfortunately, I got the error:

ValueError: Array does not have the expected shape

Could it be because fixed_array_slice has dimensions (512, 512) and itk.image_from_xarray is expecting a 3D array? How can I solve this?

Does xarray have expand_dims? Or another way to do the equivalent of fixed_array_slice = np.expand_dims(fixed_array_slice, 0)? This assumes that xarray metadata is for 3D image, while fixed_array_slice has 2D shape.

Thank you @dzenanz! I’ve tried expand_dims and it worked correctly this way:

fixed_array_slice.expand_dims(dim = ["z"], axis = 0)

However, I want the image created after image_from_xarray to be a 2D object of dimensions (512,512), instead of the resulting (1,512,512). Is there a way to remove the z dimension from an image, converting from 3D to 2D?

itk.extract_image_filter reduces dimension :smiley:

What happens if you change
image_slice11 = itk.extract_image_filter(fixed, [512,512,1], [0,0,11])
into
image_slice11 = itk.extract_image_filter(fixed, [512,512,0], [0,0,11])
Size of the region along Z axis might need to be 0 instead of 1.

Scratch that. Looking at the example, you will need to provide ExtractionRegion, not size+index. Take a look at CreateAnImageRegion.

I have adapted the examples you’ve shared, in order to extract slice 11 from the 3D image and then set the dimensions of the image to 2D (512,512) :

indx_slice = 11
extractFilter = itk.ExtractImageFilter.New(fixed)
extractFilter.SetDirectionCollapseToSubmatrix()

# set up the extraction region [one slice]
inputRegion = fixed.GetBufferedRegion()
size = inputRegion.GetSize()
size[2] = 1  # we extract along z direction
print(size)
start = inputRegion.GetIndex()
print(start)
sliceNumber = indx_slice
start[2] = sliceNumber

Dimension = 2

size = itk.Size[Dimension]()
size.Fill(512)

index = itk.Index[Dimension]()
index.Fill(1)

RegionType = itk.ImageRegion[Dimension]
desiredRegion = RegionType()
desiredRegion.SetIndex(index)
desiredRegion.SetIndex(size)

extractFilter.SetExtractionRegion(desiredRegion)
extractFilter.Update()
fixed_slice = extractFilter.GetOutput()

The error I’m getting is:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/var/folders/zt/fmlplr_j0054qnlh5slzkxbm0000gn/T/ipykernel_4255/3506052614.py in <module>
     32 #desiredRegion.SetIndex(start)
     33 
---> 34 extractFilter.SetExtractionRegion(desiredRegion)
     35 extractFilter.Update()
     36 fixed_slice = extractFilter.GetOutput()

TypeError: in method 'itkExtractImageFilterIF3IF3_SetExtractionRegion', argument 2 of type 'itkImageRegion3'

Can you help me identify what I’m doing wrong?

Hi @SofiaCV, little bug there, use SetSize(size) instead of SetIndex(size).

The problem you are having is that your fixed image is 3D but your region is 2D, and the filter cannot accept a region of different dimension.

Try something like:

Dimension = fixed.GetImageDimension()
size = itk.Size[Dimension]()
size.Fill(512)
size[2] = 1
index = itk.Index[Dimension]()
index.Fill(0)
index[2] = sliceNumber

desiredRegion = itk.ImageRegion[Dimension]
desiredRegion.SetSize(size)
desiredRegion.SetIndex(index)
extractFilter.SetExtractionRegion(desiredRegion)
...

Not sure now if size[2] needs to be 1 or 0, but try it out as @dzenan suggested.

3 Likes

Hi @phcerdan, thank you for noticing it! Here is the code after the suggestions:

Dimension = fixed.GetImageDimension()
indx_slice = 11
extractFilter = itk.ExtractImageFilter.New(fixed)
extractFilter.SetDirectionCollapseToIdentity()

inputRegion = fixed.GetBufferedRegion()
size = inputRegion.GetSize()
size[2] = 1  # we extract along z direction
start = inputRegion.GetIndex()
sliceNumber = indx_slice
start[2] = sliceNumber

RegionType = itk.ImageRegion[Dimension]
desiredRegion = RegionType()
desiredRegion.SetIndex(start)
desiredRegion.SetSize(size)

extractFilter.SetExtractionRegion(desiredRegion)
extractFilter.Update()
fixed_slice = extractFilter.GetOutput()

It worked! Thank you all for your input.

3 Likes