How can I run a registration with 3D data? or what is the data type that I need and how do I input it? Thanks.
Preamble
Hi, firstly, I read the guidelines, and it mentioned something about “Reply as a Linked Topic”, but I could not find that option for these discussions:
- https://discourse.itk.org/t/image-registration-3d-volume/2278
- https://discourse.itk.org/t/deform-registration-3d-example-7-bspline/2344
- https://discourse.itk.org/t/3d-niigz-affine-registration/4007
- https://discourse.itk.org/t/3d-mri-image-registration/3144/4
- https://discourse.itk.org/t/3d-image-registration-with-dicom-series/3641/5
In some of these, the solution is based on SimpleITK or reference C++ code. The first link mentions some examples, which I tried to follow (MWE is based on https://itk.org/ITKExamples/src/Registration/Common/Perform2DTranslationRegistrationWithMeanSquares/Documentation.html):
Description of the issue
I am trying to perform 3D registration with the Python interface. I used Numpy to create a series of 3D arrays. I took the first and last “frames” of the series and converted them into ITK images. The number of dimensions reported by ITK is 3.
I wanted to create the data with Numpy, parse it with ITK and feed it into the registration function. When the program did not work, I thought of converting the data into VTK images and read them back with ITK.
When I use 2D images, the program works.
Expected behaviour
I want to get the registration parameters and displacements with ITK in Python.
Actual result
I get an error code about the type of data, which seems to be using 2D, but I need results for 3D.
Code (MWE, software version)
#!/usr/bin/env python
# This code is based on
# https://itk.org/ITKExamples/src/Registration/Common/Perform2DTranslationRegistrationWithMeanSquares/Documentation.html
# which holds the Apache 2.0 license, which can be found here:
# https://www.apache.org/licenses/LICENSE-2.0.html
#
# What differs from that work is licensed under the GPL 2.0
import sys
import itk
import numpy as np
from distutils.version import StrictVersion as VS
from benchmark3d.primitives import displ_series_obj3d
from skimage.morphology import disk
from skimage.transform import SimilarityTransform
from skimage.transform import warp
from scipy.ndimage import shift
# Draw a happy face
side_size = 100
obj3d = np.zeros((side_size, side_size, int(side_size/10)),
dtype=np.uint8)
transf_mat_1 = SimilarityTransform(
scale=0.8, translation=(5, (100 - 80) / 2)).inverse
disk1 = np.where(warp(disk(50, dtype=np.uint8),
transf_mat_1))
transf_mat_1 = SimilarityTransform(
scale=0.8, translation=(10, (100 - 80) / 2)).inverse
disk2 = np.where(warp(disk(50, dtype=np.uint8),
transf_mat_1))
obj3d[disk1] = 1
obj3d[disk2] = 0
obj3d[19:30, 49:60, :] = 1
obj3d[69:80, 49:60, :] = 1
# Happily create 4 frames
series_length = 4
holder = np.zeros((series_length, *obj3d.shape), dtype=np.uint8)
# Move the last frame
holder[-1] = shift(holder[-1],
((series_length - 1, 0, 0)))
# Load first and last frames as ITK data (images)
fixed_image = itk.image_from_array(holder[0])
moving_image = itk.image_from_array(holder[-1])
# Hack!!! ##################################################
# itk.imwrite(itk.image_from_array(holder[0, :, :, 0]),
# "itk_hello_registration.in1.tiff")
# itk.imwrite(itk.image_from_array(holder[-1, :, :, 1]),
# "itk_hello_registration.in4.tiff")
itk.imwrite(itk.image_from_array(holder[0]),
"itk_hello_registration.in1.vtk")
itk.imwrite(itk.image_from_array(holder[-1]),
"itk_hello_registration.in4.vtk")
del fixed_image, moving_image
PixelType = itk.ctype('float')
# fixed_image = itk.imread("itk_hello_registration.in1.tiff",
# PixelType)
# moving_image = itk.imread("itk_hello_registration.in4.tiff",
# PixelType)
fixed_image = itk.ImageFileReader.New(
FileName="itk_hello_registration.in1.vtk").GetOutput()
moving_image = itk.ImageFileReader.New(
FileName="itk_hello_registration.in4.vtk").GetOutput()
# ##########################################################
# Set outputs
output_image = "itk_hello_registration.out.vtk"
diff_image_after = "itk_hello_registration.after.vtk"
diff_image_before = "itk_hello_registration.before.vtk"
if VS(itk.Version.GetITKVersion()) < VS("4.9.0"):
print("ITK 4.9.0 is required.")
sys.exit(1)
# if len(sys.argv) != 6:
# print("Usage: " + sys.argv[0] + " <fixedImageFile> <movingImageFile> "
# "<outputImagefile> <differenceImageAfter> <differenceImageBefore>")
# sys.exit(1)
#
# fixedImageFile = sys.argv[1]
# movingImageFile = sys.argv[2]
# outputImageFile = sys.argv[3]
# diffImageAfter = sys.argv[4]
# diffImageBefore = sys.argv[5]
#
PixelType = itk.ctype('float')
#
# fixed_image = itk.imread(fixedImageFile, PixelType)
# moving_image = itk.imread(movingImageFile, PixelType)
Dimension = fixed_image.GetImageDimension()
FixedImageType = itk.Image[PixelType, Dimension]
MovingImageType = itk.Image[PixelType, Dimension]
TransformType = itk.TranslationTransform[itk.D, Dimension]
initialTransform = TransformType.New()
optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
LearningRate=4,
MinimumStepLength=0.001,
RelaxationFactor=0.5,
NumberOfIterations=200)
metric = itk.MeanSquaresImageToImageMetricv4[
FixedImageType, MovingImageType].New()
registration = itk.ImageRegistrationMethodv4.New(FixedImage=fixed_image,
MovingImage=moving_image,
Metric=metric,
Optimizer=optimizer,
InitialTransform=initialTransform)
movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()
initialParameters[0] = 0
initialParameters[1] = 0
movingInitialTransform.SetParameters(initialParameters)
registration.SetMovingInitialTransform(movingInitialTransform)
identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)
registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])
registration.Update()
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)
numberOfIterations = optimizer.GetCurrentIteration()
bestValue = optimizer.GetValue()
print("Result = ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations = " + str(numberOfIterations))
print(" Metric value = " + str(bestValue))
CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())
resampler = itk.ResampleImageFilter.New(Input=moving_image,
Transform=outputCompositeTransform,
UseReferenceImage=True,
ReferenceImage=fixed_image)
resampler.SetDefaultPixelValue(100)
OutputPixelType = itk.ctype('unsigned char')
OutputImageType = itk.Image[OutputPixelType, Dimension]
caster = itk.CastImageFilter[FixedImageType,
OutputImageType].New(Input=resampler)
writer = itk.ImageFileWriter.New(Input=caster, FileName=output_image)
writer.SetFileName(output_image)
writer.Update()
difference = itk.SubtractImageFilter.New(Input1=fixed_image,
Input2=resampler)
intensityRescaler = itk.RescaleIntensityImageFilter[FixedImageType,
OutputImageType].New(
Input=difference,
OutputMinimum=itk.NumericTraits[OutputPixelType].min(),
OutputMaximum=itk.NumericTraits[OutputPixelType].max())
resampler.SetDefaultPixelValue(1)
writer.SetInput(intensityRescaler.GetOutput())
writer.SetFileName(diff_image_after)
writer.Update()
resampler.SetTransform(identityTransform)
writer.SetFileName(diff_image_before)
writer.Update()
Error message
>>>
Traceback (most recent call last):
File "<string>", line 17, in __PYTHON_EL_eval
File "/tmp/pybEA2VC", line 116, in <module>
File "/usr/lib/python3.10/itk/support/template_class.py", line 734, in New
return self[list(keys)[0]].New(*args, **kwargs)
File "/usr/lib/python3.10/itk/itkImageRegistrationMethodv4Python.py", line 248, in New
template_class.New(obj, *args, **kargs)
File "/usr/lib/python3.10/itk/support/template_class.py", line 800, in New
itk.set_inputs(self, args, kargs)
File "/usr/lib/python3.10/itk/support/extras.py", line 1252, in set_inputs
attrib(itk.output(value))
TypeError: in method 'itkImageRegistrationMethodv4REGv4F2F2_SetFixedImage', argument 2 of type 'itkImageF2 const *'
Additional information:
Wrong number or type of arguments for overloaded function 'itkImageRegistrationMethodv4REGv4F2F2_SetFixedImage'.
Possible C/C++ prototypes are:
itkImageRegistrationMethodv4REGv4F2F2::SetFixedImage(itkImageF2 const *)
itkImageRegistrationMethodv4REGv4F2F2::SetFixedImage(unsigned long,itkImageF2 const *)
Software version
Table 1: Python and ITK versionsCompilation options of ITK
Python | 3.10.1 |
ITK | 5.2.1 |
Numpy | 1.21.5 |
Scipy | 1.7.3 |
Parabola(GNU/Linux) | 5.15.12-gnu-1 |
gcc | 11.1.0 |
safe_flags="-Wp,-D_FORTIFY_SOURCE=2,-D_GLIBCXX_ASSERTIONS"
safe_flags+=" -fcf-protection -fno-plt"
safe_flags+=" -fstack-clash-protection -Wformat"
safe_flags+=" -Werror=format-security"
generic_flags="-pipe -fno-plt -fPIC -fopenmp"
generic_flags+=" -march=native"
generic_flags+=" -mtune=native ${safe_flags}"
opt_flags="${generic_flags} -O3"
generic_flags="${generic_flags} -O2"
export COPTFLAGS="${opt_flags}"
export CXXOPTFLAGS="$COPTFLAGS"
export FOPTFLAGS="$COPTFLAGS"
export CPPFLAGS="$generic_flags"
export CXXFLAGS="$CPPFLAGS"
export CFLAGS="$generic_flags"
export FFLAGS="$generic_flags"
export FCFLAGS="$generic_flags"
export F90FLAGS="$generic_flags"
export F77FLAGS="$generic_flags"
export LANG=C
export OMPI_MCA_opal_cuda_support=0
export OMPI_MCA_mpi_oversubscribe=0
_usepython=true
# https://github.com/InsightSoftwareConsortium/ITK/search?q=typeinfo&type=issues
# https://insightsoftwareconsortium.atlassian.net/browse/ITK-3538?focusedCommentId=27765&page=com.atlassian.jira.plugin.system.issuetabpanels%3Acomment-tabpanel#comment-27765
# from here: https://github.com/InsightSoftwareConsortium/ITK/issues/1701#issuecomment-598311460
export LC_ALL=C LANG=C LANGUAGE=C
_confopts=(
-S "${srcdir}"/"${pkgname}"-"${pkgver}"
-B "${srcdir}"/"${pkgname}"-"${pkgver}"/build
-DCMAKE_BUILD_TYPE:STRING=Release
-DCMAKE_CXX_FLAGS:STRING="${CXXFLAGS} -std=c++98"
-DCMAKE_VERBOSE_MAKEFILE:BOOL=OFF
-DBUILD_TESTING:BOOL=OFF
-DBUILD_EXAMPLES:BOOL=OFF
-DBUILD_SHARED_LIBS:BOOL=ON
-DCMAKE_INSTALL_PREFIX:FILEPATH=/usr
-DModule_ITKReview:BOOL=ON
-DITK_USE_SYSTEM_JPEG:BOOL=ON
-DITK_USE_SYSTEM_PNG:BOOL=ON
-DITK_USE_SYSTEM_ZLIB:BOOL=ON
-DITK_USE_SYSTEM_TIFF:BOOL=ON
-DITK_USE_SYSTEM_GDCM:BOOL=ON
-DITK_LEGACY_SILENT:BOOL=ON
$( $_usepython && echo "-DITK_WRAP_PYTHON:BOOL=ON")
$( $_usepython && echo "-DModule_ITKReview:BOOL=OFF")
$( $_usepython && echo "-DITK_USE_SYSTEM_SWIG:BOOL=ON")
$( $_usepython && echo "-DITK_USE_SYSTEM_CASTXML:BOOL=ON")
-DITK_USE_SYSTEM_LIBRARIES:BOOL=ON
-DITK_USE_SYSTEM_EXPAT:BOOL=ON
-DITK_USE_SYSTEM_FFTW:BOOL=ON
-DITK_USE_SYSTEM_HDF5:BOOL=ON
-DITK_USE_64BITS_IDS:BOOL=ON
-DITK_LEGACY_REMOVE:BOOL=ON
-DModule_ITKIOMINC:BOOL=ON
-DModule_ITKIOTransformMINC:BOOL=ON
-DModule_SimpleITKFilters:BOOL=ON
-DITK_WRAP_IMAGE_DIMS:STRING="2;3;4"
)
cmake "${_confopts[@]}"
make -j1 -s -C "${srcdir}/${pkgname}-${pkgver}"/build