Skip to content

Commit 46afb95

Browse files
hjmjohnsonblowekampthewtexdzenanz
committed
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>
1 parent 5ffbf44 commit 46afb95

4 files changed

Lines changed: 294 additions & 1 deletion

File tree

Wrapping/Generators/Python/Tests/CMakeLists.txt

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -194,6 +194,11 @@ if(ITK_WRAP_unsigned_char AND WRAP_2)
194194
COMMAND
195195
${CMAKE_CURRENT_SOURCE_DIR}/geometry_protocol.py
196196
)
197+
itk_python_add_test(
198+
NAME PythonSimpleITKProtocolTest
199+
COMMAND
200+
${CMAKE_CURRENT_SOURCE_DIR}/simpleitk_protocol.py
201+
)
197202
endif()
198203
itk_python_add_test(
199204
NAME PythonExtrasTest
Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,134 @@
1+
# ==========================================================================
2+
#
3+
# Copyright NumFOCUS
4+
#
5+
# Licensed under the Apache License, Version 2.0 (the "License");
6+
# you may not use this file except in compliance with the License.
7+
# You may obtain a copy of the License at
8+
#
9+
# https://www.apache.org/licenses/LICENSE-2.0.txt
10+
#
11+
# Unless required by applicable law or agreed to in writing, software
12+
# distributed under the License is distributed on an "AS IS" BASIS,
13+
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
# See the License for the specific language governing permissions and
15+
# limitations under the License.
16+
#
17+
# ==========================================================================
18+
19+
"""ITK's half of the SimpleITK interop contract, without SimpleITK.
20+
21+
SimpleITK builds ITK, so an ITK test that imports SimpleITK would close a
22+
cycle in the ecosystem build graph. The contract ITK owns is duck-typed --
23+
"read xyz-ordered geometry from an image-like object" -- so it is exercised
24+
here with a stub. Round-trip fidelity between two independently built ITK
25+
libraries is the business of a standalone suite that depends on neither
26+
project's build.
27+
"""
28+
29+
import itk
30+
import numpy as np
31+
32+
from itk.support.extras import _spatial_from_order_explicit
33+
34+
SPACING_XYZ = (1.0, 2.0, 3.0)
35+
36+
37+
class ImageLike:
38+
"""Stands in for a SimpleITK Image: Get*() accessors, optional keys."""
39+
40+
def __init__(self, spacing_xyz=SPACING_XYZ, keys=None):
41+
self._spacing = spacing_xyz
42+
self._keys = {} if keys is None else keys
43+
44+
def __getitem__(self, key):
45+
try:
46+
return self._keys[key]
47+
except KeyError:
48+
raise KeyError(f'"{key}" not in meta-data dictionary') from None
49+
50+
def GetSpacing(self):
51+
return self._spacing
52+
53+
54+
# --------------------------------------------------------------------------
55+
# The order-explicit key is preferred when present
56+
# --------------------------------------------------------------------------
57+
explicit = ImageLike(keys={"spacing_xyz": (7.0, 8.0, 9.0)})
58+
assert _spatial_from_order_explicit(explicit, "spacing") == (
59+
7.0,
60+
8.0,
61+
9.0,
62+
), "the _xyz key must win over the accessor when both are available"
63+
print("order-explicit key preferred over the accessor")
64+
65+
66+
# --------------------------------------------------------------------------
67+
# Fall back to the accessor: SimpleITK exposes no _xyz keys today
68+
# --------------------------------------------------------------------------
69+
# Verified against SimpleITK 2.5.4: image['spacing_xyz'] raises KeyError and
70+
# the class has no keys(). The accessor path therefore carries every real
71+
# conversion, so it is not a decorative fallback.
72+
accessor_only = ImageLike()
73+
assert (
74+
_spatial_from_order_explicit(accessor_only, "spacing") == SPACING_XYZ
75+
), "must fall back to GetSpacing() when no order-explicit key exists"
76+
print("accessor fallback used when no order-explicit key is present")
77+
78+
79+
# --------------------------------------------------------------------------
80+
# The bare key is never consulted, even when it is the only key
81+
# --------------------------------------------------------------------------
82+
# This is the #6706 conflict: a bare 'spacing' means (z,y,x) on an itk.Image
83+
# and (x,y,z) on a SimpleITK Image. Reading it would silently reverse the
84+
# spacing. Give the stub a bare key that disagrees with its accessor and
85+
# require the accessor to win.
86+
trap = ImageLike(keys={"spacing": SPACING_XYZ[::-1]})
87+
assert (
88+
_spatial_from_order_explicit(trap, "spacing") == SPACING_XYZ
89+
), "the order-ambiguous bare key must never be read (#6706)"
90+
print("order-ambiguous bare key is never read")
91+
92+
93+
# --------------------------------------------------------------------------
94+
# Absent entirely: report nothing rather than guess a default
95+
# --------------------------------------------------------------------------
96+
class Empty:
97+
def __getitem__(self, key):
98+
raise KeyError(key)
99+
100+
101+
assert (
102+
_spatial_from_order_explicit(Empty(), "spacing") is None
103+
), "must return None when neither the key nor the accessor exists"
104+
print("missing geometry reported as None rather than defaulted")
105+
106+
107+
# --------------------------------------------------------------------------
108+
# The index the converter records for a non-zero buffered region
109+
# --------------------------------------------------------------------------
110+
offset_image = itk.Image[itk.F, 3].New()
111+
region = itk.ImageRegion[3]()
112+
region.SetIndex([2, 3, 4])
113+
region.SetSize([4, 5, 6])
114+
offset_image.SetRegions(region)
115+
offset_image.Allocate(True)
116+
assert tuple(offset_image["index_xyz"]) == (2, 3, 4), "index_xyz must be xyz"
117+
assert tuple(offset_image["index_zyx"]) == (4, 3, 2), "index_zyx must be zyx"
118+
print("non-zero start index is readable in both orders (#6710)")
119+
120+
121+
# --------------------------------------------------------------------------
122+
# ITK's own order-explicit keys, which the converter writes against
123+
# --------------------------------------------------------------------------
124+
image = itk.Image[itk.F, 3].New()
125+
image.SetRegions([4, 5, 6])
126+
image.Allocate(True)
127+
image.SetSpacing(SPACING_XYZ)
128+
assert np.allclose(image["spacing_xyz"], SPACING_XYZ), "itk spacing_xyz is not xyz"
129+
assert np.allclose(
130+
image["spacing_zyx"], SPACING_XYZ[::-1]
131+
), "itk spacing_zyx is not zyx"
132+
print("itk.Image exposes both spatial key orders (#6710)")
133+
134+
print("simpleitk_protocol test passed")

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

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,8 @@
9595
"image_from_xarray",
9696
"vtk_image_from_image",
9797
"image_from_vtk_image",
98+
"image_from_simpleitk",
99+
"simpleitk_from_image",
98100
"dict_from_image",
99101
"image_from_dict",
100102
"image_intensity_min_max",
@@ -786,6 +788,138 @@ def image_from_vtk_image(vtk_image: vtk.vtkImageData) -> itkt.ImageBase:
786788
return l_image
787789

788790

791+
_ORIGINAL_INDEX_KEY = "ITK_original_index"
792+
793+
# Explicit rather than built from the name: an unknown name must fail loudly,
794+
# not resolve to a missing accessor and silently leave ITK defaults in place.
795+
_SPATIAL_ACCESSORS = {
796+
"spacing": "GetSpacing",
797+
"origin": "GetOrigin",
798+
"direction": "GetDirection",
799+
}
800+
801+
802+
def _spatial_from_order_explicit(obj, name: str):
803+
"""Read one xyz-ordered spatial attribute, never the order-ambiguous bare key (#6706)."""
804+
accessor_name = _SPATIAL_ACCESSORS[name]
805+
if hasattr(obj, "__getitem__"):
806+
try:
807+
return obj[f"{name}_xyz"]
808+
except KeyError:
809+
pass
810+
accessor = getattr(obj, accessor_name, None)
811+
return None if accessor is None else accessor()
812+
813+
814+
def image_from_simpleitk(sitk_image) -> itkt.ImageBase:
815+
"""Convert a SimpleITK Image to an itk.Image.
816+
817+
Geometry is read in ITK (x, y, z) order, via the order-explicit
818+
``spacing_xyz``/``origin_xyz``/``direction_xyz`` keys when the object
819+
provides them and otherwise via ``GetSpacing()``/``GetOrigin()``/
820+
``GetDirection()``. Pixels are copied through SimpleITK's array API.
821+
Multi-component images become an itk.VectorImage. Entries reported by
822+
``GetMetaDataKeys()`` are copied into the MetaDataDictionary.
823+
824+
Parameters
825+
----------
826+
sitk_image :
827+
A SimpleITK.Image.
828+
829+
Returns
830+
-------
831+
image :
832+
The resulting itk.Image, or itk.VectorImage for multi-component pixels.
833+
"""
834+
import itk
835+
import SimpleITK as sitk
836+
837+
dim = sitk_image.GetDimension()
838+
number_of_components = sitk_image.GetNumberOfComponentsPerPixel()
839+
is_vector = number_of_components != 1
840+
841+
array = sitk.GetArrayFromImage(sitk_image)
842+
843+
l_image = itk.image_view_from_array(array, is_vector=is_vector)
844+
845+
spacing = _spatial_from_order_explicit(sitk_image, "spacing")
846+
if spacing is not None:
847+
l_image.SetSpacing([float(s) for s in spacing])
848+
origin = _spatial_from_order_explicit(sitk_image, "origin")
849+
if origin is not None:
850+
l_image.SetOrigin([float(o) for o in origin])
851+
direction = _spatial_from_order_explicit(sitk_image, "direction")
852+
if direction is not None:
853+
l_image.SetDirection(np.asarray(direction, dtype=np.float64).reshape(dim, dim))
854+
855+
metadata = l_image.GetMetaDataDictionary()
856+
for key in sitk_image.GetMetaDataKeys():
857+
metadata[key] = sitk_image.GetMetaData(key)
858+
859+
if sitk_image.HasMetaDataKey(_ORIGINAL_INDEX_KEY):
860+
index = [int(i) for i in sitk_image.GetMetaData(_ORIGINAL_INDEX_KEY).split()]
861+
region = l_image.GetLargestPossibleRegion()
862+
region.SetIndex(index)
863+
l_image.SetRegions(region)
864+
# Undo the origin shift, so index and origin describe the same image again.
865+
shift = np.asarray(l_image.GetDirection()) @ (
866+
np.asarray(l_image.GetSpacing()) * np.asarray(index)
867+
)
868+
l_image.SetOrigin([float(o) for o in np.asarray(l_image.GetOrigin()) - shift])
869+
metadata.Erase(_ORIGINAL_INDEX_KEY)
870+
return l_image
871+
872+
873+
def simpleitk_from_image(image: itkt.ImageOrImageSource):
874+
"""Convert an itk.Image to a SimpleITK Image.
875+
876+
The inverse of :func:`image_from_simpleitk`. Geometry is transferred in ITK
877+
(x, y, z) order and the MetaDataDictionary is copied across.
878+
879+
SimpleITK images always start at index 0, so a non-zero buffered-region
880+
start index is recorded as metadata and restored by the inverse.
881+
"""
882+
import itk
883+
import SimpleITK as sitk
884+
885+
image = itk.output(image)
886+
887+
array = itk.array_from_image(image)
888+
is_vector = image.GetNumberOfComponentsPerPixel() != 1
889+
sitk_image = sitk.GetImageFromArray(array, isVector=is_vector)
890+
891+
start_index = tuple(image.GetBufferedRegion().GetIndex())
892+
893+
sitk_image.SetSpacing([float(s) for s in image.GetSpacing()])
894+
# The first stored voxel becomes index 0, so the origin moves with it and
895+
# the pixels keep their physical location.
896+
sitk_image.SetOrigin(
897+
[float(o) for o in image.TransformIndexToPhysicalPoint(list(start_index))]
898+
)
899+
sitk_image.SetDirection(
900+
[float(d) for d in np.asarray(image.GetDirection()).ravel()]
901+
)
902+
903+
metadata = image.GetMetaDataDictionary()
904+
for key in metadata.GetKeys():
905+
try:
906+
sitk_image.SetMetaData(key, str(metadata[key]))
907+
except (RuntimeError, TypeError):
908+
# Only ITK_-prefixed entries must survive the round trip.
909+
if key.startswith("ITK_"):
910+
raise
911+
912+
# Written after the copy, and cleared when zero, so the key always
913+
# reflects this image's region rather than a stale copied entry.
914+
if any(i != 0 for i in start_index):
915+
sitk_image.SetMetaData(
916+
_ORIGINAL_INDEX_KEY, " ".join(str(i) for i in start_index)
917+
)
918+
elif sitk_image.HasMetaDataKey(_ORIGINAL_INDEX_KEY):
919+
sitk_image.EraseMetaData(_ORIGINAL_INDEX_KEY)
920+
return sitk_image
921+
922+
789923
def dict_from_image(image: itkt.Image) -> dict:
790924
"""Serialize a Python itk.Image object to a pickable Python dictionary."""
791925
import itk

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

Lines changed: 21 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,13 @@
3939
except importlib.metadata.PackageNotFoundError:
4040
pass
4141

42+
_HAVE_SIMPLEITK = False
43+
try:
44+
metadata("SimpleITK")
45+
_HAVE_SIMPLEITK = True
46+
except importlib.metadata.PackageNotFoundError:
47+
pass
48+
4249

4350
def snake_to_camel_case(keyword: str):
4451
# Helpers for set_inputs snake case to CamelCase keyword argument conversion
@@ -105,6 +112,13 @@ def accept_array_like_xarray_torch(image_filter):
105112
import xarray as xr
106113
if _HAVE_TORCH:
107114
import torch
115+
sitk = None
116+
if _HAVE_SIMPLEITK:
117+
try:
118+
# A half-installed SimpleITK must not take down `import itk`.
119+
import SimpleITK as sitk
120+
except ImportError:
121+
pass
108122

109123
@functools.wraps(image_filter)
110124
def image_filter_wrapper(*args, **kwargs):
@@ -126,6 +140,9 @@ def image_filter_wrapper(*args, **kwargs):
126140
arr = move_last_dimension_to_first(arr)
127141
image = itk.image_view_from_array(arr, is_vector=channels > 1)
128142
args_list[index] = image
143+
elif sitk is not None and isinstance(arg, sitk.Image):
144+
# Not flagged as array input: outputs stay itk.Image.
145+
args_list[index] = itk.image_from_simpleitk(arg)
129146
elif not isinstance(arg, itk.Object) and is_arraylike(arg):
130147
have_array_input = True
131148
array = np.asarray(arg)
@@ -149,6 +166,8 @@ def image_filter_wrapper(*args, **kwargs):
149166
arr = move_last_dimension_to_first(arr)
150167
image = itk.image_view_from_array(arr, is_vector=channels > 1)
151168
kwargs[key] = image
169+
elif sitk is not None and isinstance(value, sitk.Image):
170+
kwargs[key] = itk.image_from_simpleitk(value)
152171
elif not isinstance(value, itk.Object) and is_arraylike(value):
153172
have_array_input = True
154173
array = np.asarray(value)
@@ -194,7 +213,8 @@ def image_filter_wrapper(*args, **kwargs):
194213
output = itk.array_view_from_image(output)
195214
return output
196215
else:
197-
return image_filter(*args, **kwargs)
216+
# args_list carries conversions that do not request output conversion back.
217+
return image_filter(*tuple(args_list), **kwargs)
198218

199219
return image_filter_wrapper
200220

0 commit comments

Comments
 (0)