Creating a displacement field image from a numpy array

Beginner Questions
1

Hello,

in my PhD project I analyze 3D microCT datasets of lung tissue samples to find out more about the alveolar shape. One topic is the simulation of an atelectasis by warping the image. In order to do that (using the WarpImageFilter or the ResampleImageFilter) I have to create a displacement vector field (with ITK Python). Therefore, I have to convert a 3D numpy array into an itk image using the GetImageFromArray function.

Here´s my code:

array1 = []

for i in range (-5,5):
    for j in range(-5,5):
        for k in range(-5,5):
            if i == 0 and j == 0 and k == 0:
                array1.append([0, 0, 0])
            else:
                x = (float(i)/float(i**2 + j**2 + k**2))
                y = (float(j)/float(i**2 + j**2 + k**2))
                z = (float(k)/float(i**2 + j**2 + k**2))
                array1.append([x, y, z])


displacementFieldFileName = itk.image_from_array(np.reshape(array1, (10,10,10,3)), is_vector = True)

In the last line, I receive the following error message:

Traceback (most recent call last):
File “Test_Displacement.py”, line 39, in
displacementFieldFileName = itk.image_from_array(np.reshape(array1, (10,10,10,3)), is_vector = True)
File “/home/reimelta/.local/lib/python2.7/site-packages/itkExtras.py”, line 297, in GetImageFromArray
return _GetImageFromArray(arr, “GetImageFromArray”, is_vector)
File “/home/reimelta/.local/lib/python2.7/site-packages/itkExtras.py”, line 291, in _GetImageFromArray
templatedFunction = getattr(itk.PyBuffer[ImageType], function)
File “/home/reimelta/.local/lib/python2.7/site-packages/itkTemplate.py”, line 340, in getitem
raise TemplateTypeError(self, tuple(cleanParameters))
itkTemplate.TemplateTypeError: itk.PyBuffer is not wrapped for input type itk.Image[itk.Vector[itk.D,3],3].

A similar topic can be found here:

Any help will be highly appreciated!
Alex

1 Like
2

You should make your NumPy array 32-bit floats, by throwing in .astype(np.float32 or specifying somethere dtype=np.float32. Or build ITK yourself and make sure to enable wrapping for 3-dimensional Vector images and double scalar type.

3

Thank you for your reply! I tried out the first two suggestions but am getting the following error with both of them:

Traceback (most recent call last):
  File "Test_Displacement.py", line 61, in <module>
    fieldReader.SetFileName(displacementFieldFileName)
TypeError: in method 'itkImageFileReaderIF3_SetFileName', argument 2 of type 'std::string const &'

Here is the code in case I did not use the suggestions correctly:

array1 = []

for i in range (-5,5):
    for j in range(-5,5):
        for k in range(-5,5):
            if i == 0 and j == 0 and k == 0:
                array1.append([0, 0, 0])
            else:
                x = ((float(i)/float(i**2 + j**2 + k**2)))
                y = ((float(j)/float(i**2 + j**2 + k**2)))
                z = ((float(k)/float(i**2 + j**2 + k**2)))
                
                array1.append([x, y, z])

array2 = np.array(array1, dtype = np.float32)

displacementFieldFileName = itk.image_from_array(np.reshape(array2, (10,10,10,3)), is_vector = True)

or:

array1 = []

for i in range (-5,5):
    for j in range(-5,5):
        for k in range(-5,5):
            if i == 0 and j == 0 and k == 0:
                array1.append([0, 0, 0])
            else:
                x = ((float(i)/float(i**2 + j**2 + k**2)))
                y = ((float(j)/float(i**2 + j**2 + k**2)))
                z = ((float(k)/float(i**2 + j**2 + k**2)))
                
                array1.append([x, y, z])

displacementFieldFileName = itk.image_from_array(np.reshape(array1, (10,10,10,3)).astype(np.float32), is_vector = True)

What am I doing wrong?

Thanks in advance
Alex

4

Where is line 61? I don’t see ImageFileReader anywhere in your code. Some other part of your code might be giving you the issue now.

5

Not sure if that error report is complete, check that your displacementFieldFileName is actually a python string. '....'

6

I´m sorry, the the last line was missing, here you can see the whole script:

#!/usr/bin/env python
import itk
import numpy as np
import sys

if len(sys.argv) != 3:
    print('Usage: ' + sys.argv[0] +
          ' <InputFileName> <OutputFileName>')
    sys.exit(1)

inputFileName = sys.argv[1]
outputFileName = sys.argv[2]

Dimension = 3

VectorComponentType = itk.F
VectorPixelType = itk.Vector[VectorComponentType, Dimension]

DisplacementFieldType = itk.Image[itk.F, Dimension]

PixelType = itk.UC
ImageType = itk.Image[PixelType, Dimension]


array1 = []

for i in range (-5,5):
    for j in range(-5,5):
        for k in range(-5,5):
            if i == 0 and j == 0 and k == 0:
                array1.append([0, 0, 0])
            else:
                x = ((float(i)/float(i**2 + j**2 + k**2)))
                y = ((float(j)/float(i**2 + j**2 + k**2)))
                z = ((float(k)/float(i**2 + j**2 + k**2)))
                
                array1.append([x, y, z])

displacementFieldFileName = itk.image_from_array(np.reshape(array1, (10,10,10,3)).astype(np.float32), is_vector = True)


reader = itk.ImageFileReader[ImageType].New()
reader.SetFileName(inputFileName)

fieldReader = itk.ImageFileReader[DisplacementFieldType].New()
fieldReader.SetFileName(displacementFieldFileName)
fieldReader.Update()

deformationField = fieldReader.GetOutput()

warpFilter = \
        itk.WarpImageFilter[ImageType, ImageType, DisplacementFieldType].New()

interpolator = itk.LinearInterpolateImageFunction[ImageType, itk.D].New()

warpFilter.SetInterpolator(interpolator)

warpFilter.SetOutputSpacing(deformationField.GetSpacing())
warpFilter.SetOutputOrigin(deformationField.GetOrigin())
warpFilter.SetOutputDirection(deformationField.GetDirection())

warpFilter.SetDisplacementField(deformationField)

warpFilter.SetInput(reader.GetOutput())

writer = itk.ImageFileWriter[ImageType].New()
writer.SetInput(warpFilter.GetOutput())
writer.SetFileName(outputFileName)

writer.Update()
7

In your script, displacementFieldFileName is not a FileName, but an ITK image, change its name directly to deformationField. You can remove all the fieldReader part.

8

Thanks for the advice. I will do so!

9

The script now delivers an output (and no error message), but only for small arrays. If the given array is 10 by 10 by 10 or larger, the output file is all black. The dataset that I would like to process has 500 voxels in each direction. Could it be that the image origin and the origin of the deformation field do not fit to each other or something similar?

Here´s how the code with the correct size would look like:

#!/usr/bin/env python
import itk
import numpy as np
import sys

if len(sys.argv) != 3:
    print('Usage: ' + sys.argv[0] +
          ' <InputFileName> <OutputFileName>')
    sys.exit(1)

inputFileName = sys.argv[1]
outputFileName = sys.argv[2]

Dimension = 3

VectorComponentType = itk.F
VectorPixelType = itk.Vector[VectorComponentType, Dimension]

DisplacementFieldType = itk.Image[VectorPixelType, Dimension]

PixelType = itk.UC
ImageType = itk.Image[PixelType, Dimension]


array1 = []

for i in range (-250,250):
    for j in range(-250,250):
        for k in range(-250,250):
            if i == 0 and j == 0 and k == 0:
                array1.append([0, 0, 0])
            else:
                x = (float(i)/float(i**2 + j**2 + k**2))
                y = (float(j)/float(i**2 + j**2 + k**2))
                z = (float(k)/float(i**2 + j**2 + k**2))
                
                array1.append([x, y, z])

deformationField = itk.image_from_array(np.reshape(array1, (500,500,500,3)).astype(np.float32), is_vector = True)

reader = itk.ImageFileReader[ImageType].New()
reader.SetFileName(inputFileName)

warpFilter = \
        itk.WarpImageFilter[ImageType, ImageType, DisplacementFieldType].New()

interpolator = itk.LinearInterpolateImageFunction[ImageType, itk.D].New()

warpFilter.SetInterpolator(interpolator)

warpFilter.SetOutputSpacing(deformationField.GetSpacing())
warpFilter.SetOutputOrigin(deformationField.GetOrigin())
warpFilter.SetOutputDirection(deformationField.GetDirection())

warpFilter.SetDisplacementField(deformationField)

warpFilter.SetInput(reader.GetOutput())

writer = itk.ImageFileWriter[ImageType].New()
writer.SetInput(warpFilter.GetOutput())
writer.SetFileName(outputFileName)

writer.Update()

Thank you again for your help.

1 Like
10

First, I fail to find a reader.Update() after reader.SetFileName, so reader->GetOutput() is not populated.
Second, I am not familiar with the WarpImageFilter, but I am sure this is not doing what you expect:

You haven’t set the Spacing, Origin and Direction of your deformationField image, so those are the default. If you want to copy that metadata from your input image:

reader = itk.ImageFileReader[ImageType].New()
reader.SetFileName(inputFileName)
reader.Update()
deformationField = itk.image_from_array(np.reshape(array1, (500,500,500,3)).astype(np.float32), is_vector = True)
deformationField.SetSpacing(reader.GetOutput().GetSpacing())
deformationField.SetOrigin(reader.GetOutput().GetOrigin())
deformationField.SetDirection(reader.GetOutput().GetDirection())

Hope it helps!

1 Like
11

Thank you! I will try it out.

12

It worked out kind of how I wanted it to be. The only smaller issue is that on the first and last 5 to 6 images of my stack there are empty spaces. But they can be cut out or I will try to use the ResampleImage Filter.

Thank you both!

3 Likes