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..54e89e06287 --- /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") + + +# -------------------------------------------------------------------------- +# 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") diff --git a/Wrapping/Generators/Python/itk/support/extras.py b/Wrapping/Generators/Python/itk/support/extras.py index 92c28abc0b3..3af878ad53e 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 +_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 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