ENH: Expose displacement field transform adaptors in Python

Modules/Filtering/DisplacementField/wrapping/itkDisplacementFieldTransform.wrap

@@ -1,8 +1,3 @@
 itk_wrap_class("itk::DisplacementFieldTransform" POINTER)
-foreach(d ${ITK_WRAP_IMAGE_DIMS})
-  itk_wrap_template("${ITKM_D}${d}" "${ITKT_D},${d}")
-  if(ITK_WRAP_vector_float)
-    itk_wrap_template("${ITKM_F}${d}" "${ITKT_F},${d}")
-  endif()
-endforeach()
+set(WRAPPER_TEMPLATES "${itk_Wrap_DisplacementFieldTransform}")
 itk_end_wrap_class()

Modules/Registration/Common/wrapping/CMakeLists.txt

@@ -1,2 +1,10 @@
 itk_wrap_module(ITKRegistrationCommon)
+set(displacement_field_adaptor_templates)
+foreach(d ${ITK_WRAP_IMAGE_DIMS})
+  list(
+    APPEND
+    displacement_field_adaptor_templates
+    "${ITKM_DFTD${d}} # ${ITKT_DFTD${d}}"
+  )
+endforeach()
 itk_auto_load_and_end_wrap_submodules()

Modules/Registration/Common/wrapping/itkDisplacementFieldTransformParametersAdaptor.wrap

@@ -0,0 +1,5 @@
+itk_wrap_include("itkDisplacementFieldTransform.h")
+
+itk_wrap_class("itk::DisplacementFieldTransformParametersAdaptor" POINTER)
+set(WRAPPER_TEMPLATES "${displacement_field_adaptor_templates}")
+itk_end_wrap_class()

Modules/Registration/Common/wrapping/itkTransformParametersAdaptor.wrap

@@ -0,0 +1,5 @@
+itk_wrap_include("itkDisplacementFieldTransform.h")
+
+itk_wrap_class("itk::TransformParametersAdaptor" POINTER)
+set(WRAPPER_TEMPLATES "${displacement_field_adaptor_templates}")
+itk_end_wrap_class()

Modules/Registration/Common/wrapping/itkTransformParametersAdaptorBase.wrap

@@ -1,9 +1,17 @@
 itk_wrap_include("itkTransform.h")
 itk_wrap_class("itk::TransformParametersAdaptorBase" POINTER)
 unique(real_types "${WRAP_ITK_REAL};D")
-foreach(d ${ITK_WRAP_DIMS})
+foreach(d ${ITK_WRAP_IMAGE_DIMS})
   foreach(t ${real_types})
     itk_wrap_template("${ITKM_${t}}${d}" "itk::Transform< ${ITKT_${t}},${d},${d} >")
+    set(adaptor_name "itkTransformParametersAdaptorBase${ITKM_${t}}${d}")
+    string(APPEND ITK_WRAP_PYTHON_SWIG_EXT
+      "DECL_PYTHON_STD_VEC_RAW_TO_SMARTPTR_TYPEMAP(${adaptor_name}, ${adaptor_name}_Pointer)\n"
+      "%template(vector${adaptor_name}) std::vector< ${adaptor_name}_Pointer >;\n"
+    )
+    add_python_config_template("vector" "std::vector" "vector${adaptor_name}"
+      "itk::TransformParametersAdaptorBase< itk::Transform< ${ITKT_${t}},${d},${d} > >"
+    )
   endforeach()
 endforeach()
 itk_end_wrap_class()

Modules/Registration/RegistrationMethodsv4/wrapping/test/CMakeLists.txt

@@ -1,4 +1,24 @@
 if(ITK_WRAP_PYTHON)
+  foreach(d 2 3)
+    list(FIND ITK_WRAP_IMAGE_DIMS ${d} wrap_dimension_index)
+    if(
+      ITKBridgeNumPy_ENABLED
+      AND
+        ITK_WRAP_float
+      AND
+        wrap_dimension_index
+          GREATER
+          -1
+    )
+      itk_python_add_test(
+        NAME itkSyNImageRegistrationMethodMultiresolution${d}DPythonTest
+        COMMAND
+          ${CMAKE_CURRENT_SOURCE_DIR}/itkSyNImageRegistrationMethodTest.py
+          ${d}
+      )
+    endif()
+  endforeach()
+
   itk_python_expression_add_test(
     NAME itkBSplineSyNImageRegistrationMethodPythonTest
     EXPRESSION "instance = itk.BSplineSyNImageRegistrationMethod.New()"

Modules/Registration/RegistrationMethodsv4/wrapping/test/itkSyNImageRegistrationMethodTest.py

@@ -0,0 +1,130 @@
+# ==========================================================================
+#
+#   Copyright NumFOCUS
+#
+#   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
+#
+#          https://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
+
+
+dimension = int(sys.argv[1])
+image_size = 32
+blob_sigma = 4.0
+shift = 2.0
+smoothing_sigmas = [1, 0]
+iterations_per_level = [20, 20]
+geometry_tolerance = 1e-6
+center = np.full(dimension, (image_size - 1) / 2)
+moving_center = center.copy()
+moving_center[0] += shift
+coordinates = np.indices((image_size,) * dimension, dtype=np.float32)
+fixed_squared_radius = sum(
+    (coordinates[dimension - 1 - axis] - center[axis]) ** 2 for axis in range(dimension)
+)
+moving_squared_radius = sum(
+    (coordinates[dimension - 1 - axis] - moving_center[axis]) ** 2
+    for axis in range(dimension)
+)
+image = itk.image_from_array(
+    np.exp(-fixed_squared_radius / (2 * blob_sigma**2)).astype(np.float32)
+)
+moving_image = itk.image_from_array(
+    np.exp(-moving_squared_radius / (2 * blob_sigma**2)).astype(np.float32)
+)
+Image = type(image)
+Transform = itk.DisplacementFieldTransform[itk.D, dimension]
+Adaptor = itk.DisplacementFieldTransformParametersAdaptor[Transform]
+AdaptorBase = itk.TransformParametersAdaptorBase[
+    itk.Transform[itk.D, dimension, dimension]
+]
+
+for shrink_factors in ([2, 1], [4, 2], [1, 1]):
+    registration = itk.SyNImageRegistrationMethod[Image, Image].New(
+        FixedImage=image,
+        MovingImage=moving_image,
+        Metric=itk.MeanSquaresImageToImageMetricv4[Image, Image].New(),
+    )
+    Field = itk.Image[itk.Vector[itk.D, dimension], dimension]
+    fields = []
+    for _ in range(2):
+        field = Field.New()
+        field.CopyInformation(image)
+        field.SetRegions(image.GetBufferedRegion())
+        field.AllocateInitialized()
+        fields.append(field)
+    initial_transform = Transform.New(
+        DisplacementField=fields[0], InverseDisplacementField=fields[1]
+    )
+    registration.SetInitialTransform(initial_transform)
+    registration.SetNumberOfLevels(len(shrink_factors))
+    registration.SetShrinkFactorsPerLevel(shrink_factors)
+    registration.SetSmoothingSigmasPerLevel(smoothing_sigmas)
+    registration.SetNumberOfIterationsPerLevel(iterations_per_level)
+
+    adaptors = itk.vector[AdaptorBase]()
+    for factor in shrink_factors:
+        shrink = itk.ShrinkImageFilter[Image, Image].New(
+            Input=image, ShrinkFactors=factor
+        )
+        shrink.Update()
+        level_image = shrink.GetOutput()
+        adaptor = Adaptor.New()
+        adaptor.SetRequiredSpacing(level_image.GetSpacing())
+        adaptor.SetRequiredSize(level_image.GetBufferedRegion().GetSize())
+        adaptor.SetRequiredDirection(level_image.GetDirection())
+        adaptor.SetRequiredOrigin(level_image.GetOrigin())
+        adaptors.push_back(adaptor)
+
+    try:
+        registration.SetTransformParametersAdaptorsPerLevel(itk.vector[AdaptorBase]())
+    except RuntimeError:
+        pass
+    else:
+        raise AssertionError("An adaptor is required for each level")
+
+    registration.SetTransformParametersAdaptorsPerLevel(adaptors)
+    assert len(registration.GetTransformParametersAdaptorsPerLevel()) == len(
+        shrink_factors
+    )
+    registration.Update()
+
+    transform = registration.GetModifiableTransform()
+    mapped_center = np.array(transform.TransformPoint(center))
+    assert np.linalg.norm(mapped_center - moving_center) < shift / 2
+    inverse = transform.GetInverseTransform()
+    mapped_moving_center = np.array(inverse.TransformPoint(moving_center))
+    assert np.linalg.norm(mapped_moving_center - center) < shift / 2
+
+    for field in (
+        transform.GetDisplacementField(),
+        transform.GetInverseDisplacementField(),
+    ):
+        assert tuple(field.GetLargestPossibleRegion().GetSize()) == tuple(
+            level_image.GetLargestPossibleRegion().GetSize()
+        )
+        np.testing.assert_allclose(
+            field.GetOrigin(), level_image.GetOrigin(), atol=geometry_tolerance
+        )
+        np.testing.assert_allclose(
+            field.GetSpacing(), level_image.GetSpacing(), atol=geometry_tolerance
+        )
+        np.testing.assert_allclose(
+            itk.array_from_matrix(field.GetDirection()),
+            itk.array_from_matrix(level_image.GetDirection()),
+            atol=geometry_tolerance,
+        )

Wrapping/WrapITKTypes.cmake

@@ -82,6 +82,16 @@ endforeach()
 end_wrap_type()
 set(itk_Wrap_Vector ${WRAPPER_TEMPLATES})
 
+wrap_type("itk::DisplacementFieldTransform" "DFT")
+foreach(d ${ITK_WRAP_IMAGE_DIMS})
+  add_template("${ITKM_D}${d}" "${ITKT_D},${d}")
+  if(ITK_WRAP_vector_float)
+    add_template("${ITKM_F}${d}" "${ITKT_F},${d}")
+  endif()
+endforeach()
+end_wrap_type()
+set(itk_Wrap_DisplacementFieldTransform ${WRAPPER_TEMPLATES})
+
 wrap_type("itk::CovariantVector" "CV" "itkCovariantVector.h")
 foreach(vector_dim ${ITK_WRAP_VECTOR_COMPONENTS_INCREMENTED})
   add_template("${ITKM_F}${vector_dim}" "${ITKT_F},${vector_dim}")