# Creating a displacement field image from a numpy array

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

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.

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>
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?

Alex

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.

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

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 +
' <InputFileName> <OutputFileName>')
sys.exit(1)

inputFileName = sys.argv
outputFileName = sys.argv

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)

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)

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

writer.Update()

``````

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.

Thanks for the advice. I will do so!

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 +
' <InputFileName> <OutputFileName>')
sys.exit(1)

inputFileName = sys.argv
outputFileName = sys.argv

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)

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)

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

writer.Update()

``````

Thank you again for your help.

1 Like

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()
deformationField = itk.image_from_array(np.reshape(array1, (500,500,500,3)).astype(np.float32), is_vector = True)