ENH: Add SimpleITK <-> ITK image conversion to the Python interface (supersedes #6021)

Wrapping/Generators/Python/Tests/CMakeLists.txt

@@ -194,6 +194,11 @@ if(ITK_WRAP_unsigned_char AND WRAP_2)
         COMMAND
           ${CMAKE_CURRENT_SOURCE_DIR}/geometry_protocol.py
       )
+      itk_python_add_test(
+        NAME PythonSimpleITKProtocolTest
+        COMMAND
+          ${CMAKE_CURRENT_SOURCE_DIR}/simpleitk_protocol.py
+      )
     endif()
     itk_python_add_test(
       NAME PythonExtrasTest

Wrapping/Generators/Python/Tests/simpleitk_protocol.py

@@ -0,0 +1,134 @@
+# ==========================================================================
+#
+#   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.
+#
+# ==========================================================================
+
+"""ITK's half of the SimpleITK interop contract, without SimpleITK.
+
+SimpleITK builds ITK, so an ITK test that imports SimpleITK would close a
+cycle in the ecosystem build graph. The contract ITK owns is duck-typed --
+"read xyz-ordered geometry from an image-like object" -- so it is exercised
+here with a stub. Round-trip fidelity between two independently built ITK
+libraries is the business of a standalone suite that depends on neither
+project's build.
+"""
+
+import itk
+import numpy as np
+
+from itk.support.extras import _spatial_from_order_explicit
+
+SPACING_XYZ = (1.0, 2.0, 3.0)
+
+
+class ImageLike:
+    """Stands in for a SimpleITK Image: Get*() accessors, optional keys."""
+
+    def __init__(self, spacing_xyz=SPACING_XYZ, keys=None):
+        self._spacing = spacing_xyz
+        self._keys = {} if keys is None else keys
+
+    def __getitem__(self, key):
+        try:
+            return self._keys[key]
+        except KeyError:
+            raise KeyError(f'"{key}" not in meta-data dictionary') from None
+
+    def GetSpacing(self):
+        return self._spacing
+
+
+# --------------------------------------------------------------------------
+# The order-explicit key is preferred when present
+# --------------------------------------------------------------------------
+explicit = ImageLike(keys={"spacing_xyz": (7.0, 8.0, 9.0)})
+assert _spatial_from_order_explicit(explicit, "spacing") == (
+    7.0,
+    8.0,
+    9.0,
+), "the _xyz key must win over the accessor when both are available"
+print("order-explicit key preferred over the accessor")
+
+
+# --------------------------------------------------------------------------
+# Fall back to the accessor: SimpleITK exposes no _xyz keys today
+# --------------------------------------------------------------------------
+# Verified against SimpleITK 2.5.4: image['spacing_xyz'] raises KeyError and
+# the class has no keys(). The accessor path therefore carries every real
+# conversion, so it is not a decorative fallback.
+accessor_only = ImageLike()
+assert (
+    _spatial_from_order_explicit(accessor_only, "spacing") == SPACING_XYZ
+), "must fall back to GetSpacing() when no order-explicit key exists"
+print("accessor fallback used when no order-explicit key is present")
+
+
+# --------------------------------------------------------------------------
+# The bare key is never consulted, even when it is the only key
+# --------------------------------------------------------------------------
+# This is the #6706 conflict: a bare 'spacing' means (z,y,x) on an itk.Image
+# and (x,y,z) on a SimpleITK Image. Reading it would silently reverse the
+# spacing. Give the stub a bare key that disagrees with its accessor and
+# require the accessor to win.
+trap = ImageLike(keys={"spacing": SPACING_XYZ[::-1]})
+assert (
+    _spatial_from_order_explicit(trap, "spacing") == SPACING_XYZ
+), "the order-ambiguous bare key must never be read (#6706)"
+print("order-ambiguous bare key is never read")
+
+
+# --------------------------------------------------------------------------
+# Absent entirely: report nothing rather than guess a default
+# --------------------------------------------------------------------------
+class Empty:
+    def __getitem__(self, key):
+        raise KeyError(key)
+
+
+assert (
+    _spatial_from_order_explicit(Empty(), "spacing") is None
+), "must return None when neither the key nor the accessor exists"
+print("missing geometry reported as None rather than defaulted")
+
+
+# --------------------------------------------------------------------------
+# The index the converter records for a non-zero buffered region
+# --------------------------------------------------------------------------
+offset_image = itk.Image[itk.F, 3].New()
+region = itk.ImageRegion[3]()
+region.SetIndex([2, 3, 4])
+region.SetSize([4, 5, 6])
+offset_image.SetRegions(region)
+offset_image.Allocate(True)
+assert tuple(offset_image["index_xyz"]) == (2, 3, 4), "index_xyz must be xyz"
+assert tuple(offset_image["index_zyx"]) == (4, 3, 2), "index_zyx must be zyx"
+print("non-zero start index is readable in both orders (#6710)")
+
+
+# --------------------------------------------------------------------------
+# ITK's own order-explicit keys, which the converter writes against
+# --------------------------------------------------------------------------
+image = itk.Image[itk.F, 3].New()
+image.SetRegions([4, 5, 6])
+image.Allocate(True)
+image.SetSpacing(SPACING_XYZ)
+assert np.allclose(image["spacing_xyz"], SPACING_XYZ), "itk spacing_xyz is not xyz"
+assert np.allclose(
+    image["spacing_zyx"], SPACING_XYZ[::-1]
+), "itk spacing_zyx is not zyx"
+print("itk.Image exposes both spatial key orders (#6710)")
+
+print("simpleitk_protocol test passed")

Wrapping/Generators/Python/itk/support/extras.py

@@ -95,6 +95,8 @@
     "image_from_xarray",
     "vtk_image_from_image",
     "image_from_vtk_image",
+    "image_from_simpleitk",
+    "simpleitk_from_image",
     "dict_from_image",
     "image_from_dict",
     "image_intensity_min_max",
@@ -786,6 +788,138 @@ def image_from_vtk_image(vtk_image: vtk.vtkImageData) -> itkt.ImageBase:
     return l_image
 
 
+_ORIGINAL_INDEX_KEY = "ITK_original_index"
+
+# Explicit rather than built from the name: an unknown name must fail loudly,
+# not resolve to a missing accessor and silently leave ITK defaults in place.
+_SPATIAL_ACCESSORS = {
+    "spacing": "GetSpacing",
+    "origin": "GetOrigin",
+    "direction": "GetDirection",
+}
+
+
+def _spatial_from_order_explicit(obj, name: str):
+    """Read one xyz-ordered spatial attribute, never the order-ambiguous bare key (#6706)."""
+    accessor_name = _SPATIAL_ACCESSORS[name]
+    if hasattr(obj, "__getitem__"):
+        try:
+            return obj[f"{name}_xyz"]
+        except KeyError:
+            pass
+    accessor = getattr(obj, accessor_name, None)
+    return None if accessor is None else accessor()
+
+
+def image_from_simpleitk(sitk_image) -> itkt.ImageBase:
+    """Convert a SimpleITK Image to an itk.Image.
+
+    Geometry is read in ITK (x, y, z) order, via the order-explicit
+    ``spacing_xyz``/``origin_xyz``/``direction_xyz`` keys when the object
+    provides them and otherwise via ``GetSpacing()``/``GetOrigin()``/
+    ``GetDirection()``. Pixels are copied through SimpleITK's array API.
+    Multi-component images become an itk.VectorImage. Entries reported by
+    ``GetMetaDataKeys()`` are copied into the MetaDataDictionary.
+
+    Parameters
+    ----------
+    sitk_image :
+        A SimpleITK.Image.
+
+    Returns
+    -------
+    image :
+        The resulting itk.Image, or itk.VectorImage for multi-component pixels.
+    """
+    import itk
+    import SimpleITK as sitk
+
+    dim = sitk_image.GetDimension()
+    number_of_components = sitk_image.GetNumberOfComponentsPerPixel()
+    is_vector = number_of_components != 1
+
+    array = sitk.GetArrayFromImage(sitk_image)
+
+    l_image = itk.image_view_from_array(array, is_vector=is_vector)
+
+    spacing = _spatial_from_order_explicit(sitk_image, "spacing")
+    if spacing is not None:
+        l_image.SetSpacing([float(s) for s in spacing])
+    origin = _spatial_from_order_explicit(sitk_image, "origin")
+    if origin is not None:
+        l_image.SetOrigin([float(o) for o in origin])
+    direction = _spatial_from_order_explicit(sitk_image, "direction")
+    if direction is not None:
+        l_image.SetDirection(np.asarray(direction, dtype=np.float64).reshape(dim, dim))
+
+    metadata = l_image.GetMetaDataDictionary()
+    for key in sitk_image.GetMetaDataKeys():
+        metadata[key] = sitk_image.GetMetaData(key)
+
+    if sitk_image.HasMetaDataKey(_ORIGINAL_INDEX_KEY):
+        index = [int(i) for i in sitk_image.GetMetaData(_ORIGINAL_INDEX_KEY).split()]
+        region = l_image.GetLargestPossibleRegion()
+        region.SetIndex(index)
+        l_image.SetRegions(region)
+        # Undo the origin shift, so index and origin describe the same image again.
+        shift = np.asarray(l_image.GetDirection()) @ (
+            np.asarray(l_image.GetSpacing()) * np.asarray(index)
+        )
+        l_image.SetOrigin([float(o) for o in np.asarray(l_image.GetOrigin()) - shift])
+        metadata.Erase(_ORIGINAL_INDEX_KEY)
+    return l_image
+
+
+def simpleitk_from_image(image: itkt.ImageOrImageSource):
+    """Convert an itk.Image to a SimpleITK Image.
+
+    The inverse of :func:`image_from_simpleitk`. Geometry is transferred in ITK
+    (x, y, z) order and the MetaDataDictionary is copied across.
+
+    SimpleITK images always start at index 0, so a non-zero buffered-region
+    start index is recorded as metadata and restored by the inverse.
+    """
+    import itk
+    import SimpleITK as sitk
+
+    image = itk.output(image)
+
+    array = itk.array_from_image(image)
+    is_vector = image.GetNumberOfComponentsPerPixel() != 1
+    sitk_image = sitk.GetImageFromArray(array, isVector=is_vector)
+
+    start_index = tuple(image.GetBufferedRegion().GetIndex())
+
+    sitk_image.SetSpacing([float(s) for s in image.GetSpacing()])
+    # The first stored voxel becomes index 0, so the origin moves with it and
+    # the pixels keep their physical location.
+    sitk_image.SetOrigin(
+        [float(o) for o in image.TransformIndexToPhysicalPoint(list(start_index))]
+    )
+    sitk_image.SetDirection(
+        [float(d) for d in np.asarray(image.GetDirection()).ravel()]
+    )
+
+    metadata = image.GetMetaDataDictionary()
+    for key in metadata.GetKeys():
+        try:
+            sitk_image.SetMetaData(key, str(metadata[key]))
+        except (RuntimeError, TypeError):
+            # Only ITK_-prefixed entries must survive the round trip.
+            if key.startswith("ITK_"):
+                raise
+
+    # Written after the copy, and cleared when zero, so the key always
+    # reflects this image's region rather than a stale copied entry.
+    if any(i != 0 for i in start_index):
+        sitk_image.SetMetaData(
+            _ORIGINAL_INDEX_KEY, " ".join(str(i) for i in start_index)
+        )
+    elif sitk_image.HasMetaDataKey(_ORIGINAL_INDEX_KEY):
+        sitk_image.EraseMetaData(_ORIGINAL_INDEX_KEY)
+    return sitk_image
+
+
 def dict_from_image(image: itkt.Image) -> dict:
     """Serialize a Python itk.Image object to a pickable Python dictionary."""
     import itk

Wrapping/Generators/Python/itk/support/helpers.py

@@ -39,6 +39,13 @@
 except importlib.metadata.PackageNotFoundError:
     pass
 
+_HAVE_SIMPLEITK = False
+try:
+    metadata("SimpleITK")
+    _HAVE_SIMPLEITK = True
+except importlib.metadata.PackageNotFoundError:
+    pass
+
 
 def snake_to_camel_case(keyword: str):
     # Helpers for set_inputs snake case to CamelCase keyword argument conversion
@@ -105,6 +112,13 @@ def accept_array_like_xarray_torch(image_filter):
         import xarray as xr
     if _HAVE_TORCH:
         import torch
+    sitk = None
+    if _HAVE_SIMPLEITK:
+        try:
+            # A half-installed SimpleITK must not take down `import itk`.
+            import SimpleITK as sitk
+        except ImportError:
+            pass
 
     @functools.wraps(image_filter)
     def image_filter_wrapper(*args, **kwargs):
@@ -126,6 +140,9 @@ def image_filter_wrapper(*args, **kwargs):
                     arr = move_last_dimension_to_first(arr)
                 image = itk.image_view_from_array(arr, is_vector=channels > 1)
                 args_list[index] = image
+            elif sitk is not None and isinstance(arg, sitk.Image):
+                # Not flagged as array input: outputs stay itk.Image.
+                args_list[index] = itk.image_from_simpleitk(arg)
             elif not isinstance(arg, itk.Object) and is_arraylike(arg):
                 have_array_input = True
                 array = np.asarray(arg)
@@ -149,6 +166,8 @@ def image_filter_wrapper(*args, **kwargs):
                         arr = move_last_dimension_to_first(arr)
                     image = itk.image_view_from_array(arr, is_vector=channels > 1)
                     kwargs[key] = image
+                elif sitk is not None and isinstance(value, sitk.Image):
+                    kwargs[key] = itk.image_from_simpleitk(value)
                 elif not isinstance(value, itk.Object) and is_arraylike(value):
                     have_array_input = True
                     array = np.asarray(value)
@@ -194,7 +213,8 @@ def image_filter_wrapper(*args, **kwargs):
                         output = itk.array_view_from_image(output)
                 return output
         else:
-            return image_filter(*args, **kwargs)
+            # args_list carries conversions that do not request output conversion back.
+            return image_filter(*tuple(args_list), **kwargs)
 
     return image_filter_wrapper