How to `itk.resample_image_filter` with a `BSplineTransform` with 3D images?

Beginner Questions
, , ,
1

1. How to itk.resample_image_filter with a BSplineTransform with 3D images?

Hello, everyone!

I am trying to use parts of the DeformableRegistration15.cxx in Python, but I don’t know what to do when ResampleImageFilter is combined with GetCoefficientImages (from a BSplineTransform).

1.1. What I have tried

1.1.1. Inputs (data)

Listing 1: Block to print the input to ResampleImageFilter and resample_image_filter.
print("moving_image_type:\t", moving_image_type)
print("coord_type:\t", moving_image_type)
print("identity_transform:\t", identity_transform)
print("bspline_transform_coarse\t", bspline_transform_coarse)
print("bspline_transform_fine\t", bspline_transform_fine)
interpolator_fine = itk.BSplineResampleImageFunction[
    moving_image_type, coord_type].New()

moving_image_type:	 <class 'itk.itkImagePython.itkImageF3'>
coord_type:	 <class 'itk.itkImagePython.itkImageF3'>
identity_transform:	 IdentityTransform (0x55a08d43e540)
  RTTI typeinfo:   itk::IdentityTransform<double, 3u>
  Reference Count: 1
  Modified Time: 127
  Debug: Off
  Object Name:
  Observers:
    none

bspline_transform_coarse	 BSplineTransform (0x55a08c2b42b0)
  RTTI typeinfo:   itk::BSplineTransform<double, 3u, 3u>
  Reference Count: 3
  Modified Time: 14797
  Debug: Off
  Object Name:
  Observers:
    none
  CoefficientImage: [ 0x55a08cd32630, 0x55a08cd328f0, 0x55a08cade1e0 ]
  TransformDomainOrigin: [0, 0, 0]
  TransformDomainPhysicalDimensions: [99, 99, 9]
  TransformDomainDirection: 1 0 0
0 1 0
0 0 1

  TransformDomainMeshSize: [2, 2, 2]
  GridSize: [5, 5, 5]
  GridOrigin: [-49.5, -49.5, -4.5]
  GridSpacing: [49.5, 49.5, 4.5]
  GridDirection: 1 0 0
0 1 0
0 0 1


bspline_transform_fine	 BSplineTransform (0x55a08f8fa480)
  RTTI typeinfo:   itk::BSplineTransform<double, 3u, 3u>
  Reference Count: 1
  Modified Time: 14848
  Debug: Off
  Object Name:
  Observers:
    none
  CoefficientImage: [ 0x55a08f198a00, 0x55a08f198cc0, 0x55a08f198f80 ]
  TransformDomainOrigin: [0, 0, 0]
  TransformDomainPhysicalDimensions: [99, 99, 9]
  TransformDomainDirection: 1 0 0
0 1 0
0 0 1

  TransformDomainMeshSize: [2, 2, 2]
  GridSize: [5, 5, 5]
  GridOrigin: [-49.5, -49.5, -4.5]
  GridSpacing: [49.5, 49.5, 4.5]
  GridDirection: 1 0 0
0 1 0
0 0 1

1.1.2. Failed attempt with resample_image_filter

Listing 2: This block fails probably because I don’t know how to provide the GetCoefficientImages.
# https://examples.itk.org/src/filtering/imagegrid/resampleavectorimage/documentation
upsampler = itk.resample_image_filter(
    bspline_transform_coarse.GetCoefficientImages(),
    interpolator=interpolator_fine,
    transform=identity_transform,
    size=bspline_transform_fine.GetCoefficientImages())

swig/python detected a memory leak of type 'itk::FixedArray< itk::SmartPointer< itk::Image< double,3 > >,3 > *', no destructor found.
swig/python detected a memory leak of type 'itk::FixedArray< itk::SmartPointer< itk::Image< double,3 > >,3 > *', no destructor found.
Traceback (most recent call last):
  File "<string>", line 17, in __PYTHON_EL_eval
  File "<string>", line 2, in <module>
  File "/usr/lib/python3.10/site-packages/itk/support/helpers.py", line 176, in image_filter_wrapper
    return image_filter(*args, **kwargs)
  File "/usr/lib/python3.10/site-packages/itk/itkResampleImageFilterPython.py", line 5542, in resample_image_filter
    instance = itk.ResampleImageFilter.New(*args, **kwargs)
  File "/usr/lib/python3.10/site-packages/itk/support/template_class.py", line 734, in New
    raise itk.TemplateTypeError(self, input_type)
itk.support.extras.TemplateTypeError: itk.ResampleImageFilter is not wrapped for input type `None`.

To limit the size of the package, only a limited number of
types are available in ITK Python. To print the supported
types, run the following command in your python environment:

    itk.ResampleImageFilter.GetTypes()

Possible solutions:
* If you are an application user:
** Convert your input image into a supported format (see below).
** Contact developer to report the issue.
* If you are an application developer, force input images to be
loaded in a supported pixel type.

    e.g.: instance = itk.ResampleImageFilter[itk.Image[itk.SC,2], itk.Image[itk.SC,2]].New(my_input)

* (Advanced) If you are an application developer, build ITK Python yourself and
turned to `ON` the corresponding CMake option to wrap the pixel type or image
dimension you need. When configuring ITK with CMake, you can set
`ITK_WRAP_${type}` (replace ${type} with appropriate pixel type such as
`double`). If you need to support images with 4 or 5 dimensions, you can add
these dimensions to the list of dimensions in the CMake variable
`ITK_WRAP_IMAGE_DIMS`.

Supported input types:

itk.Image[itk.SC,2]
itk.Image[itk.SC,3]
itk.Image[itk.SC,4]
itk.Image[itk.SS,2]
itk.Image[itk.SS,3]
itk.Image[itk.SS,4]
itk.Image[itk.UC,2]
itk.Image[itk.UC,3]
itk.Image[itk.UC,4]
itk.Image[itk.US,2]
itk.Image[itk.US,3]
itk.Image[itk.US,4]
itk.Image[itk.F,2]
itk.Image[itk.F,3]
itk.Image[itk.F,4]
itk.Image[itk.D,2]
itk.Image[itk.D,3]
itk.Image[itk.D,4]
itk.Image[itk.Vector[itk.D,2],2]
itk.Image[itk.Vector[itk.D,2],3]
itk.Image[itk.Vector[itk.D,2],4]
itk.Image[itk.Vector[itk.D,3],2]
itk.Image[itk.Vector[itk.D,3],3]
itk.Image[itk.Vector[itk.D,3],4]
itk.Image[itk.Vector[itk.D,4],2]
itk.Image[itk.Vector[itk.D,4],3]
itk.Image[itk.Vector[itk.D,4],4]
itk.Image[itk.Vector[itk.F,2],2]
itk.Image[itk.Vector[itk.F,2],3]
itk.Image[itk.Vector[itk.F,2],4]
itk.Image[itk.Vector[itk.F,3],2]
itk.Image[itk.Vector[itk.F,3],3]
itk.Image[itk.Vector[itk.F,3],4]
itk.Image[itk.Vector[itk.F,4],2]
itk.Image[itk.Vector[itk.F,4],3]
itk.Image[itk.Vector[itk.F,4],4]
itk.Image[itk.RGBPixel[itk.UC],2]
itk.Image[itk.RGBPixel[itk.UC],3]
itk.Image[itk.RGBPixel[itk.UC],4]
itk.Image[itk.RGBPixel[itk.US],2]
itk.Image[itk.RGBPixel[itk.US],3]
itk.Image[itk.RGBPixel[itk.US],4]
itk.Image[itk.RGBAPixel[itk.UC],2]
itk.Image[itk.RGBAPixel[itk.UC],3]
itk.Image[itk.RGBAPixel[itk.UC],4]
itk.Image[itk.RGBAPixel[itk.US],2]
itk.Image[itk.RGBAPixel[itk.US],3]
itk.Image[itk.RGBAPixel[itk.US],4]
itk.VectorImage[itk.UC,2]
itk.VectorImage[itk.SC,2]
itk.VectorImage[itk.SS,2]
itk.VectorImage[itk.US,2]
itk.VectorImage[itk.F,2]
itk.VectorImage[itk.D,2]
itk.VectorImage[itk.UC,3]
itk.VectorImage[itk.SC,3]
itk.VectorImage[itk.SS,3]
itk.VectorImage[itk.US,3]
itk.VectorImage[itk.F,3]
itk.VectorImage[itk.D,3]
itk.VectorImage[itk.UC,4]
itk.VectorImage[itk.SC,4]
itk.VectorImage[itk.SS,4]
itk.VectorImage[itk.US,4]
itk.VectorImage[itk.F,4]
itk.VectorImage[itk.D,4]

1.1.3. Useless result with ResampleImageFilter

Listing 3: This block works.
# https://examples.itk.org/src/filtering/imagegrid/resampleavectorimage/documentation
upsampler = itk.ResampleImageFilter[moving_image_type, moving_image_type].New()
print(upsampler)

ResampleImageFilter (0x55a08bde5270)
  RTTI typeinfo:   itk::ResampleImageFilter<itk::Image<float, 3u>, itk::Image<float, 3u>, double, double>
  Reference Count: 1
  Modified Time: 15126
  Debug: Off
  Object Name:
  Observers:
    none
  Inputs:
    Primary: (0) *
    ReferenceImage: (0)
    Transform: (0x55a08ffa49c0) *
  Indexed Inputs:
    0: Primary (0)
    1: ReferenceImage (0)
  Required Input Names: Primary, Transform
  NumberOfRequiredInputs: 1
  Outputs:
    Primary: (0x55a08ee2f840)
  Indexed Outputs:
    0: Primary (0x55a08ee2f840)
  NumberOfRequiredOutputs: 1
  Number Of Work Units: 32
  ReleaseDataFlag: Off
  ReleaseDataBeforeUpdateFlag: Off
  AbortGenerateData: Off
  Progress: 0
  Multithreader:
    RTTI typeinfo:   itk::PoolMultiThreader
    Reference Count: 1
    Modified Time: 15106
    Debug: Off
    Object Name:
    Observers:
      none
    Number of Work Units: 32
    Number of Threads: 8
    Global Maximum Number Of Threads: 128
    Global Default Number Of Threads: 8
    Global Default Threader Type: itk::MultiThreaderBaseEnums::Threader::Pool
    SingleMethod: 0
    SingleData: 0
  DynamicMultiThreading: On
  CoordinateTolerance: 1e-06
  DirectionTolerance: 1e-06
  DefaultPixelValue: 0
  Size: [0, 0, 0]
  OutputStartIndex: [0, 0, 0]
  OutputSpacing: [1, 1, 1]
  OutputOrigin: [0, 0, 0]
  OutputDirection: 1 0 0
0 1 0
0 0 1

  Transform: 0x55a08f6d3800
  Interpolator: 0x55a08fd5e4b0
  Extrapolator: 0
  UseReferenceImage: Off

1.1.4. Error with ResampleImageFilter and GetCoefficientImages()

Listing 4: This block will not work, because of my wrong call to GetCoefficientImages.
# https://examples.itk.org/src/filtering/imagegrid/resampleavectorimage/documentation
upsampler = itk.ResampleImageFilter[moving_image_type, moving_image_type].New(
    Input=bspline_transform_coarse.GetCoefficientImages(),
    Interpolator=interpolator_fine,
    Transform=identity_transform,
    Size=bspline_transform_fine.GetCoefficientImages())

Traceback (most recent call last):
  File "<string>", line 17, in __PYTHON_EL_eval
  File "<string>", line 2, in <module>
  File "/usr/lib/python3.10/site-packages/itk/itkResampleImageFilterPython.py", line 516, in New
    template_class.New(obj, *args, **kargs)
  File "/usr/lib/python3.10/site-packages/itk/support/template_class.py", line 800, in New
    itk.set_inputs(self, args, kargs)
  File "/usr/lib/python3.10/site-packages/itk/support/extras.py", line 1543, in set_inputs
    attrib(itk.output(value))
TypeError: Expecting argument of type itkImageF3 or itkImageSourceIF3.
Additional information:
Wrong number or type of arguments for overloaded function 'itkImageToImageFilterIF3IF3_SetInput'.
  Possible C/C++ prototypes are:
    itkImageToImageFilterIF3IF3::SetInput(itkImageF3 const *)
    itkImageToImageFilterIF3IF3::SetInput(unsigned int,itkImageF3 const *)

1.1.5. Error with ResampleImageFilter and GetCoefficientImages()[0]

Listing 5: This block will not work, because of my wrong call to GetCoefficientImages (or trying to access an inexisting element).
# https://examples.itk.org/src/filtering/imagegrid/resampleavectorimage/documentation
upsampler = itk.ResampleImageFilter[moving_image_type, moving_image_type].New(
    Input=bspline_transform_coarse.GetCoefficientImages()[0],
    Interpolator=interpolator_fine,
    Transform=identity_transform,
    Size=bspline_transform_fine.GetCoefficientImages())

swig/python detected a memory leak of type 'itk::FixedArray< itk::SmartPointer< itk::Image< double,3 > >,3 > *', no destructor found.
Traceback (most recent call last):
  File "<string>", line 17, in __PYTHON_EL_eval
  File "<string>", line 3, in <module>
TypeError: 'SwigPyObject' object is not subscriptable

1.2. Software

import itk
import sys
print("Self-compiled version of ITK:\t", itk.__version__)
print("Python version:\t", sys.version)

Self-compiled version of ITK:	 5.3.0
Python version:	 3.10.4 (main, Mar 23 2022, 23:05:40) [GCC 11.2.0]
2

Hi @edgar ,

To resample an image with itk.resample_image_filter and a b-spline transform, pass the transform in as the transform argument.

For the input, use the moving image. The sampling grid can be specified with a reference image, e.g. the fixed image, or by the parameters of the sampling gride, including the ‘size’.

HTH,
Matt

3

Hi and thanks, @matt.mccormick. Sorry for the wall of text. In DeformableRegistration15.cxx, GetCoefficientImages()[k] is used within a loop to provide the input image. How can I access that image in Python? I get an error if I try (see above). Thanks!

4

Hi, I’m sure everyone is busy, and it’s alright. If someone has the time, I would appreciate the help. Thanks.

5

Can you please check if the data type is correct ?
It expects float format as per the error message.

6

Also please post a minimum end-to-end Python code to reproduce the error.

7

About the code: sure, here it is.

#!/usr/bin/env python
# This code is partially based on
# DeformableRegistration15.cxx
# which holds the Apache 2.0 license (see below).
#
# Copyright 2022 edgar
# License:
#
#    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
#    <https://www.gnu.org/licenses/>.
#
# Copyright Insight Software Consortium
#
# Licensed under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0.txt
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on
# an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

import sys
import itk
import numpy as np
from distutils.version import StrictVersion as VS
from scipy import ndimage as ndi
from skimage.morphology import disk
from scipy.ndimage import shift


if VS(itk.Version.GetITKVersion()) < VS("5.3.0"):
    print("ITK 5.3.0 is required.")
    sys.exit(1)

# 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))

fixed_image_arr = holder[0].transpose((2, 1, 0))
moving_image_arr = holder[-1].transpose((2, 1, 0))
fixed_image = itk.image_from_array(
    fixed_image_arr.astype(np.float32))
moving_image = itk.image_from_array(
    moving_image_arr.astype(np.float32))

spatial_dim = fixed_image.GetImageDimension()
fixed_image_type = type(fixed_image)
moving_image_type = type(moving_image)

origin = [0] * spatial_dim
length_per_px_ratios = [1] * spatial_dim
[dict(img).update({
    # Set origin ([0, 0, 0]); fixed_image.GetOrigin()
    "origin": origin,
    # Axial scale (length unit / pixel for each axis)
    "spacing": length_per_px_ratios,
    # Set orientation of axes; fixed_image.GetDirection()
    "direction": np.eye(spatial_dim)})
 for img in (fixed_image, moving_image)]

coord_type = itk.D
versor_type = coord_type
versor_transform_class = itk.VersorRigid3DTransform[versor_type]
rigid_transform = versor_transform_class.New()
trnfm_init_class = itk.CenteredTransformInitializer[
    # versor_transform_class,
    itk.VersorRigid3DTransform[itk.D],
    fixed_image_type,
    moving_image_type]
optimiser_class = itk.RegularStepGradientDescentOptimizer
optimiser = optimiser_class.New()
metric = itk.MattesMutualInformationImageToImageMetric[
    fixed_image_type, moving_image_type].New(
        NumberOfHistogramBins=50)
metric.ReinitializeSeed(76926294)
metric.SetUseCachingOfBSplineWeights(True)
metric.SetUseExplicitPDFDerivatives(True)
interpolator = itk.LinearInterpolateImageFunction[
    moving_image_type, coord_type].New()
registrant = itk.ImageRegistrationMethod[
    fixed_image_type, moving_image_type].New(
        Metric=metric,
        Optimizer=optimiser,
        Interpolator=interpolator,
        FixedImage=fixed_image,
        MovingImage=moving_image)

identity_transform = itk.IdentityTransform[
    coord_type, spatial_dim].New()

initialiser = trnfm_init_class.New(
    Transform=rigid_transform,
    FixedImage=fixed_image,
    MovingImage=moving_image)
initialiser.MomentsOn()
initialiser.InitializeTransform()
fixed_region = fixed_image.GetBufferedRegion()
number_of_pixels = fixed_region.GetNumberOfPixels()
registrant.SetFixedImageRegion(fixed_region)
registrant.SetInitialTransformParameters(
    rigid_transform.GetParameters())
registrant.SetTransform(rigid_transform)

translation_scale = 1 / 1000
optimiser_scales = [1.0, 1.0, 1.0,
                    translation_scale,
                    translation_scale,
                    translation_scale]
optimiser.SetScales(optimiser_scales)
optimiser.SetMaximumStepLength(0.2000)
optimiser.SetMinimumStepLength(0.0001)
optimiser.SetNumberOfIterations(200)

metric.SetNumberOfSpatialSamples(int(10000))
registrant.Update()
rigid_transform.SetParameters(registrant.GetLastTransformParameters())

affine_transform_class = itk.AffineTransform[
    coord_type, spatial_dim]
affine_transform = affine_transform_class.New(
    Center=rigid_transform.GetCenter(),
    Translation=rigid_transform.GetTranslation(),
    Matrix=rigid_transform.GetMatrix())
registrant.SetTransform(affine_transform)
registrant.SetInitialTransformParameters(
    affine_transform.GetParameters())
optimiser.SetScales([1] * 9 + [translation_scale] * 3)
optimiser.SetMaximumStepLength(0.2)
optimiser.SetMinimumStepLength(1e-4)
optimiser.SetNumberOfIterations(200)
metric.SetNumberOfSpatialSamples(int(50000))

registrant.Update()

affine_transform.SetParameters(registrant.GetLastTransformParameters())

spline_order = 3
deformable_transform_class = itk.BSplineTransform[
    coord_type, spatial_dim, spline_order]
bspline_transform_coarse = deformable_transform_class.New()
grid_nodes_1_dim_coarse = int(5)
fixed_origin = fixed_image.GetOrigin()
fixed_physical_dimensions = [fixed_image.GetSpacing()[idx] * (
    fixed_image.GetLargestPossibleRegion().GetSize()[idx] - 1)
                             for idx in range(spatial_dim)]
bspline_transform_coarse.SetTransformDomainOrigin(fixed_origin)
bspline_transform_coarse.SetTransformDomainPhysicalDimensions(
    fixed_physical_dimensions)
mesh_size = \
    bspline_transform_coarse.GetTransformDomainMeshSize()
mesh_size.Fill(grid_nodes_1_dim_coarse - spline_order)
bspline_transform_coarse.SetTransformDomainMeshSize(mesh_size)
bspline_transform_coarse.SetTransformDomainDirection(
    fixed_image.GetDirection())
number_of_bspline_param = \
    bspline_transform_coarse.GetNumberOfParameters()
optimiser.SetScales(
    [1.0 for idx in range(number_of_bspline_param)])
composite_transform_class = itk.CompositeTransform[
    coord_type, spatial_dim]
composite_transform = composite_transform_class.New()
composite_transform.AddTransform(affine_transform)
composite_transform.AddTransform(bspline_transform_coarse)
composite_transform.SetOnlyMostRecentTransformToOptimizeOn()
bspline_transform_coarse.GetParameters().Fill(1.0)
registrant.SetInitialTransformParameters(
    bspline_transform_coarse.GetParameters())
registrant.SetTransform(composite_transform)

optimiser.SetMaximumStepLength(10.0)
optimiser.SetMinimumStepLength(0.01)
optimiser.SetRelaxationFactor(0.7)
optimiser.SetNumberOfIterations(50)

metric.SetNumberOfSpatialSamples(number_of_bspline_param * 100)
registrant.Update()
final_parameters = registrant.GetLastTransformParameters()
bspline_transform_coarse.SetParameters(final_parameters)


bspline_transform_fine = deformable_transform_class.New(
    TransformDomainOrigin=fixed_origin,
    TransformDomainPhysicalDimensions=fixed_physical_dimensions,
    TransformDomainMeshSize=mesh_size,
    TransformDomainDirection=fixed_image.GetDirection())
grid_nodes_1_dim_fine = int(20)
mesh_size.Fill(grid_nodes_1_dim_fine - spline_order)
number_of_bspline_param = \
    bspline_transform_fine.GetNumberOfParameters()

interpolator_fine = itk.BSplineResampleImageFunction[
    moving_image_type, coord_type].New()

upsampler = itk.ResampleImageFilter[moving_image_type, moving_image_type].New(
    Input=bspline_transform_coarse.GetCoefficientImages(),
    Interpolator=interpolator_fine,
    Transform=identity_transform,
    Size=bspline_transform_fine.GetCoefficientImages())

upsampler = itk.resample_image_filter[moving_image_type, moving_image_type](
    Input=bspline_transform_coarse.GetCoefficientImages()[0],
    interpolator=interpolator_fine,
    transform=identity_transform,
    size=bspline_transform_fine.GetCoefficientImages())

About the data, the idea is to use GetCoefficientImages as the input (just as in the example), but let me know if this is the information that you need:

moving_image_type:	 <class 'itk.itkImagePython.itkImageF3'>
fixed_image_type:	 <class 'itk.itkImagePython.itkImageF3'>

Thank you!

8

If I print the object using the following code,

t1 = bspline_transform_coarse.GetCoefficientImages()
print(t1)

I get:

<Swig Object of type 'itk::FixedArray< itk::SmartPointer< itk::Image< double,3 > >,3 > *' at 0x7f0e70fa7480>

which is a fixed array of image with double type. So you can’t directly pass this as input.

9

Thanks. In https://itk.org/Doxygen/html/Examples_2RegistrationITKv4_2DeformableRegistration15_8cxx-example.html, they use
bspline_transform_coarse.GetCoefficientImages()[k]. Is there a way to do the same with Python? How to access the k^{th} element? Is it even possible?

10

You can first get the deformation field and then upsample it and then decompose it to BSpline.

convertFilter = itk.TransformToDisplacementFieldFilter.IVF33D.New()
convertFilter.SetTransform(bsplineTransform)
convertFilter.UseReferenceImageOn()
convertFilter.SetReferenceImage(fixedImage)
convertFilter.Update()
field = convertFilter.GetOutput()
field = np.array(field)
np.save('displacement_field.npy', field)

Here, upsample the field and then pass it as input to the itk::BSplineDecompositionImageFilter to create a new B-spline transform for next iteration.

I will have to check if there is a way to image from GetCoefficientImages method in Python.

1 Like
11

Ok, thanks.

12

Will there be an implementation?