3D image registration with Python and without SimpleITK

Beginner Questions
,
1

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:

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 versions

Compilation 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
2

The main problem is that fixed_image = itk.ImageFileReader.New(FileName="itk_hello_registration.in1.vtk").GetOutput() created a UC3 image, instead of F3 image which you seem to expect.

Additionally, itk.ImageRegistrationMethodv4.New() creates registration instantiated for F2 → F2 image types by default. The trick is to specify F3.

Attached is a program which runs without error on my computer. I took a look at the files on disk, but they are all all-zeroes.

edgar_mwe.py (6.9 KB)

2 Likes
3

Hi, @dzenanz. Thank you. It also works for me (see results below). I am going to open another question regarding how to do this without creating the extra VTK file.

Result =
 Translation X = 0.0
 Translation Y = 0.0
 Translation Z = 4.0
 Iterations    = 2
 Metric value  = 0.0

You were getting all zeros, because I forgot to add a line of code. This may help:

# Draw a 3D disk
side_size = 10
obj2d = disk(side_size / 2, dtype=np.uint8)
obj3d = np.asarray([disk(side_size / 2, dtype=np.uint8)
                    for idx in range(5)])
# Create 4 frames
series_length = 4
holder = np.asarray([obj3d for idx in range(series_length)])
# Move the last frame
holder[-1] = shift(holder[-1],
                   (series_length - 1, 0, 0))

NB: To get the translation in the last coordinate, I just added

translationAlongZ = finalParameters.GetElement(2)

and

print(" Translation Z = " + str(translationAlongZ))
2 Likes