From f831d4def937d3ddbb082764b944cfa2b0804cc3 Mon Sep 17 00:00:00 2001 From: "Hans J. Johnson" Date: Thu, 27 Aug 2026 09:06:06 -0500 Subject: [PATCH] ENH: Add SimpleITK <-> ITK image conversion to the Python interface itk.image_from_simpleitk() and itk.simpleitk_from_image() convert between the two toolkits, and the filter decorator accepts a SimpleITK image wherever it accepts a NumPy array. Geometry uses the order-explicit spatial keys from #6710, falling back to the Get*() accessors. The bare 'spacing' key means (z,y,x) on an itk.Image and (x,y,z) on a SimpleITK Image, so reading it would reverse the spacing (#6706). SimpleITK images start at index 0, so the origin moves to the first stored voxel and the index is carried as ITK_original_index; the inverse restores both, keeping the pixels in the same physical location. Supersedes #6021. Co-Authored-By: Bradley Lowekamp <321061+blowekamp@users.noreply.github.com> Co-Authored-By: Matt McCormick <25432+thewtex@users.noreply.github.com> Co-Authored-By: Dzenan Zukic <1792121+dzenanz@users.noreply.github.com> --- .../Generators/Python/Tests/CMakeLists.txt | 5 + .../Python/Tests/simpleitk_protocol.py | 134 ++++++++++++++++++ .../Generators/Python/itk/support/extras.py | 134 ++++++++++++++++++ .../Generators/Python/itk/support/helpers.py | 22 ++- 4 files changed, 294 insertions(+), 1 deletion(-) create mode 100644 Wrapping/Generators/Python/Tests/simpleitk_protocol.py diff --git a/Wrapping/Generators/Python/Tests/CMakeLists.txt b/Wrapping/Generators/Python/Tests/CMakeLists.txt index 3fd4a97e55d..b7e4c3ec8d9 100644 --- a/Wrapping/Generators/Python/Tests/CMakeLists.txt +++ b/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 diff --git a/Wrapping/Generators/Python/Tests/simpleitk_protocol.py b/Wrapping/Generators/Python/Tests/simpleitk_protocol.py new file mode 100644 index 00000000000..f85678889ee --- /dev/null +++ b/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") + + +# -------------------------------------------------------------------------- +# A non-zero buffered region, readable in both key orders +# -------------------------------------------------------------------------- +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") diff --git a/Wrapping/Generators/Python/itk/support/extras.py b/Wrapping/Generators/Python/itk/support/extras.py index 92c28abc0b3..1ce7abc8a38 100644 --- a/Wrapping/Generators/Python/itk/support/extras.py +++ b/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 +# 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) + + 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. The origin is set to the physical + location of the first buffered voxel, so the pixels keep their position; the + ITK start index itself is not carried over. + + Raises + ------ + ValueError + If the buffered region differs from the largest possible region. A + SimpleITK image carries a single extent, so the distinction would be + lost without notice. + """ + import itk + import SimpleITK as sitk + + image = itk.output(image) + + # Updates the image, so the regions compared below are the ones just read. + array = itk.array_from_image(image) + + buffered_region = image.GetBufferedRegion() + largest_region = image.GetLargestPossibleRegion() + if buffered_region != largest_region: + raise ValueError( + "cannot convert an image whose buffered region differs from its " + "largest possible region: buffered index " + f"{list(buffered_region.GetIndex())} size {list(buffered_region.GetSize())}, " + f"largest index {list(largest_region.GetIndex())} size " + f"{list(largest_region.GetSize())}. A SimpleITK image holds one " + "extent, so the difference cannot be represented. Update the image " + "over its largest possible region before converting." + ) + + is_vector = image.GetNumberOfComponentsPerPixel() != 1 + sitk_image = sitk.GetImageFromArray(array, isVector=is_vector) + + start_index = tuple(buffered_region.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): + # SimpleITK stores strings only; entries that do not render are skipped. + pass + return sitk_image + + def dict_from_image(image: itkt.Image) -> dict: """Serialize a Python itk.Image object to a pickable Python dictionary.""" import itk diff --git a/Wrapping/Generators/Python/itk/support/helpers.py b/Wrapping/Generators/Python/itk/support/helpers.py index 3d271d141ee..2ccab7096ca 100644 --- a/Wrapping/Generators/Python/itk/support/helpers.py +++ b/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