May be related to Creating a displacement field image from a numpy array - #11 by Alex_Reimelt
# create_displ.py creates a displacement field with ITK
# warping does not work as expected
#
# Copyright (C) 2023 eDgar
#
# This program is free software: you can redistribute it
# and/or modify it under the terms of the GNU General Public
# License as published by the Free Software Foundation,
# version 3 of the License.
#
# This program is distributed in the hope that it will be
# useful, but WITHOUT ANY WARRANTY; without even the implied
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
# PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public
# License along with this program. If not, see
# <http://www.gnu.org/licenses/>.
import itk
import numpy as np
ndim = 3
nz = 4
ny = 30
nx = 20
img_shape = (nz, ny, nx)
img = np.zeros(img_shape)
img[:,
(ny // 3):(ny * 2 // 3),
(nx // 4):(nx * 3 // 4)] = 1
# import matplotlib.pyplot as plt
# plt.imshow(img[0])
# plt.show()
fixed_image = itk.image_from_array(
img.astype(np.float32))
VoxelType = itk.F
ImageType = itk.Image[VoxelType, ndim]
VectorVoxelType = itk.F
VectorDimension = ndim
VectorTypeDim = itk.Vector[VectorVoxelType, VectorDimension]
ImageDimension = ndim
VectorImageType = itk.Image[VectorTypeDim, ImageDimension]
# Unitary displacement field (translate)
DisplacementFieldNumpy = np.ones(
(nz, ny, nx, ndim), dtype=np.float32)
# Old-school
PyBuffer = itk.PyBuffer[VectorImageType]
# DisplacementField = PyBuffer.GetImageFromArray(
# DisplacementFieldNumpy.astype(np.float32))
DisplacementField = PyBuffer.GetImageViewFromArray(
DisplacementFieldNumpy, is_vector=True)
DisplacementField.shape
DisplacementField.SetSpacing(fixed_image.GetSpacing())
DisplacementField.SetOrigin(fixed_image.GetOrigin())
# itk.array_from_image(DisplacementField).shape
DisplacementField = itk.image_view_from_array(
DisplacementFieldNumpy,
# is_vector=True
ttype=VectorImageType
)
DisplacementField.SetSpacing(fixed_image.GetSpacing())
DisplacementField.SetOrigin(fixed_image.GetOrigin())
interpolator = itk.LinearInterpolateImageFunction[ImageType, itk.D].New()
warpFilter = itk.WarpImageFilter[ImageType, ImageType, VectorImageType].New(
Interpolator=interpolator,
OutputSpacing=fixed_image.GetSpacing(),
OutputOrigin=fixed_image.GetOrigin(),
DisplacementField=DisplacementField,
Input=fixed_image.astype(itk.F))
moving_image = warpFilter.GetOutput()
itk.imwrite(DisplacementField, "/tmp/displ_field_by_itk.vtk")
itk.imwrite(moving_image, "/tmp/warped_by_itk.vtk")
itk.imwrite(fixed_image, "/tmp/fixed_image_by_itk.vtk")