diff --git a/Documentation/docs/migration_guides/index.md b/Documentation/docs/migration_guides/index.md index 6d3d78fdab6..a2f75135e8c 100644 --- a/Documentation/docs/migration_guides/index.md +++ b/Documentation/docs/migration_guides/index.md @@ -36,5 +36,6 @@ advances with every commit. itk_6_migration_guide itk_5_migration_guide +python_spatial_key_migration joint_histogram_mutual_information_metric_correction ``` diff --git a/Documentation/docs/migration_guides/python_spatial_key_migration.md b/Documentation/docs/migration_guides/python_spatial_key_migration.md new file mode 100644 index 00000000000..0aace05aa1c --- /dev/null +++ b/Documentation/docs/migration_guides/python_spatial_key_migration.md @@ -0,0 +1,53 @@ +# Python spatial key migration + +The bare `image['origin']`, `image['spacing']`, and `image['direction']` +string keys on `itk.Image` are deprecated in favor of order-explicit keys. +The bare keys are order-ambiguous across toolkits: ITK returns NumPy +`(z, y, x)` order while SimpleITK returns `(x, y, z)` order for the same key +names, so generic code consuming these keys silently corrupts geometry when +handed the other toolkit's image +([issue #6706](https://github.com/InsightSoftwareConsortium/ITK/issues/6706)). + +## The order-explicit keys + +| Key | Order | Value | +|---|---|---| +| `origin_xyz`, `spacing_xyz` | ITK `(x, y, z)` | `GetOrigin()` / `GetSpacing()` | +| `direction_xyz` | ITK | `GetDirection()` as a NumPy matrix | +| `index_xyz`, `size_xyz` | ITK | LargestPossibleRegion index / size | +| `origin_zyx`, `spacing_zyx`, `index_zyx`, `size_zyx` | NumPy `(z, y, x)` | `np.flip` of the `_xyz` value | +| `direction_zyx` | NumPy | `np.flip(direction_xyz, axis=None)` (= P·D·P) | + +The `_zyx` frame is a complete coordinate frame matching `np.array(image)` +indexing, with physical coordinates also listed in `(z, y, x)` order: + +```python +# forward, equals np.flip(TransformContinuousIndexToPhysicalPoint(i)): +p = image['origin_zyx'] + image['direction_zyx'] @ (image['spacing_zyx'] * i) +# inverse, equals np.flip(TransformPhysicalPointToContinuousIndex(p)): +i = np.linalg.inv(image['direction_zyx']) @ (p - image['origin_zyx']) / image['spacing_zyx'] +# np.array(image)[0, 0, 0] sits at continuous index image['index_zyx'] +``` + +Writing `index_*` or `size_*` calls `SetRegions()`; a size change requires a +subsequent `Allocate()`. + +## Migrating + +The bare keys returned `(z, y, x)` order, so the drop-in replacement is the +`_zyx` key: + +```python +spacing = image['spacing'] # deprecated +spacing = image['spacing_zyx'] # identical values +spacing = image['spacing_xyz'] # GetSpacing() order, matches SimpleITK +``` + +## Warning behavior + +- By default, bare-key access emits a `DeprecationWarning` (hidden outside + `__main__` unless enabled with `-W` or `warnings.simplefilter`). +- When ITK is configured with `ITK_FUTURE_LEGACY_REMOVE=ON`, the warning is a + `FutureWarning`, which Python displays everywhere by default. +- The `ITK_PYTHON_FUTURE_LEGACY_REMOVE` environment variable overrides the + configured default at runtime (`itkConfig.FutureLegacyRemove`). diff --git a/Wrapping/Generators/Python/CMakeLists.txt b/Wrapping/Generators/Python/CMakeLists.txt index c77412c9633..9d03572d479 100644 --- a/Wrapping/Generators/Python/CMakeLists.txt +++ b/Wrapping/Generators/Python/CMakeLists.txt @@ -165,6 +165,12 @@ if(NOT EXTERNAL_WRAP_ITK_PROJECT) ${ITK_PYTHON_PACKAGE_DIR} ) + if(ITK_FUTURE_LEGACY_REMOVE) + set(ITK_WRAP_PYTHON_FUTURE_LEGACY_REMOVE "True") + else() + set(ITK_WRAP_PYTHON_FUTURE_LEGACY_REMOVE "False") + endif() + configure_file( "${CMAKE_CURRENT_SOURCE_DIR}/itkConfig.template.in.py" "${ITK_WRAP_PYTHON_ROOT_BINARY_DIR}/itkConfig.py" diff --git a/Wrapping/Generators/Python/PyBase/pyBase.i b/Wrapping/Generators/Python/PyBase/pyBase.i index 6a94404b3bf..df92df7822c 100644 --- a/Wrapping/Generators/Python/PyBase/pyBase.i +++ b/Wrapping/Generators/Python/PyBase/pyBase.i @@ -599,18 +599,44 @@ str = str These keys are used in the dictionary resulting from dict(image). These keys include MetaDataDictionary keys along with - 'origin', 'spacing', and 'direction' keys, which - correspond to the image's Origin, Spacing, and Direction. However, - they are in (z, y, x) order as opposed to (x, y, z) order to - correspond to the indexing of the shape of the pixel buffer - array resulting from np.array(image). + order-explicit spatial keys: 'origin_xyz', 'spacing_xyz', + 'direction_xyz', 'index_xyz', 'size_xyz' in ITK (x, y, z) + order, and 'origin_zyx', 'spacing_zyx', 'direction_zyx', + 'index_zyx', 'size_zyx' in NumPy (z, y, x) order matching + the indexing of np.array(image). 'index' and 'size' describe + the LargestPossibleRegion; writing them calls SetRegions(), + and a size change requires a subsequent Allocate(). + 'direction_zyx' is + np.flip(direction_xyz, axis=None), the change-of-basis with + physical coordinates also listed in (z, y, x) order. + + The bare 'origin', 'spacing', and 'direction' keys are + deprecated aliases of the '_zyx' keys; the same bare names in + SimpleITK follow (x, y, z) order, so order-ambiguous bare + keys are being retired (see issue #6706). """ meta_keys = self.GetMetaDataDictionary().GetKeys() # Ignore deprecated, legacy members that cause issues result = list(filter(lambda k: not k.startswith('ITK_original'), meta_keys)) result.extend(['origin', 'spacing', 'direction']) + result.extend(['origin_xyz', 'spacing_xyz', 'direction_xyz', 'index_xyz', 'size_xyz']) + result.extend(['origin_zyx', 'spacing_zyx', 'direction_zyx', 'index_zyx', 'size_zyx']) return result + def _warn_bare_spatial_key(self, key): + import warnings + import itkConfig + category = FutureWarning if itkConfig.FutureLegacyRemove else DeprecationWarning + warnings.warn( + f"itk.Image['{key}'] is deprecated because the bare key is " + f"order-ambiguous across toolkits: ITK returns NumPy (z,y,x) " + f"order while SimpleITK returns (x,y,z) order for the same " + f"key. Use '{key}_zyx' (same values as this key) or " + f"'{key}_xyz' (GetOrigin()/GetSpacing()/GetDirection() " + f"order). See the Python spatial key migration guide and " + f"issue #6706.", + category, stacklevel=3) + def __getitem__(self, key): """Access metadata keys, see help(image.keys), for string keys, otherwise provide NumPy indexing to the pixel buffer @@ -619,14 +645,24 @@ str = str import itk if isinstance(key, str): import numpy as np - if key == 'origin': - return np.flip(np.asarray(self.GetOrigin()), axis=None) - elif key == 'spacing': - return np.flip(np.asarray(self.GetSpacing()), axis=None) - elif key == 'direction': - return np.flip(itk.array_from_matrix(self.GetDirection()), axis=None) - else: - return self.GetMetaDataDictionary()[key] + region = self.GetLargestPossibleRegion() + xyz_getters = { + 'origin_xyz': lambda: np.asarray(self.GetOrigin()), + 'spacing_xyz': lambda: np.asarray(self.GetSpacing()), + 'direction_xyz': lambda: itk.array_from_matrix(self.GetDirection()), + 'index_xyz': lambda: np.asarray(region.GetIndex()), + 'size_xyz': lambda: np.asarray(region.GetSize()), + } + if key in xyz_getters: + return xyz_getters[key]() + if key.endswith('_zyx') and key[:-4] + '_xyz' in xyz_getters: + return np.flip(xyz_getters[key[:-4] + '_xyz'](), axis=None) + if key in ('origin', 'spacing', 'direction'): + self._warn_bare_spatial_key(key) + if key == 'direction': + return np.flip(itk.array_from_matrix(self.GetDirection()), axis=None) + return np.flip(np.asarray(getattr(self, 'Get' + key.capitalize())()), axis=None) + return self.GetMetaDataDictionary()[key] else: return itk.array_view_from_image(self).__getitem__(key) @@ -637,12 +673,36 @@ str = str order, i.e. [z, y, x] versus [x, y, z].""" if isinstance(key, str): import numpy as np - if key == 'origin': - self.SetOrigin(np.flip(value, axis=None)) - elif key == 'spacing': - self.SetSpacing(np.flip(value, axis=None)) - elif key == 'direction': - self.SetDirection(np.flip(value, axis=None)) + if key == 'origin_xyz': + self.SetOrigin(np.asarray(value, dtype=float)) + elif key == 'spacing_xyz': + self.SetSpacing(np.asarray(value, dtype=float)) + elif key == 'direction_xyz': + self.SetDirection(np.asarray(value, dtype=float)) + elif key == 'origin_zyx': + self.SetOrigin(np.flip(np.asarray(value, dtype=float))) + elif key == 'spacing_zyx': + self.SetSpacing(np.flip(np.asarray(value, dtype=float))) + elif key == 'direction_zyx': + self.SetDirection(np.flip(np.asarray(value, dtype=float), axis=None)) + elif key in ('index_xyz', 'size_xyz', 'index_zyx', 'size_zyx'): + arr = np.asarray(value).astype(int) + if key.endswith('_zyx'): + arr = np.flip(arr) + region = self.GetLargestPossibleRegion() + if key.startswith('index'): + region.SetIndex([int(v) for v in arr]) + else: + region.SetSize([int(v) for v in arr]) + self.SetRegions(region) + elif key in ('origin', 'spacing', 'direction'): + self._warn_bare_spatial_key(key) + if key == 'origin': + self.SetOrigin(np.flip(value, axis=None)) + elif key == 'spacing': + self.SetSpacing(np.flip(value, axis=None)) + else: + self.SetDirection(np.flip(value, axis=None)) else: self.GetMetaDataDictionary()[key] = value else: diff --git a/Wrapping/Generators/Python/Tests/CMakeLists.txt b/Wrapping/Generators/Python/Tests/CMakeLists.txt index c7a19b26b35..6434dcd4019 100644 --- a/Wrapping/Generators/Python/Tests/CMakeLists.txt +++ b/Wrapping/Generators/Python/Tests/CMakeLists.txt @@ -186,6 +186,11 @@ if(ITK_WRAP_unsigned_char AND WRAP_2) DATA{${WrapITK_SOURCE_DIR}/images/cthead1.png} 5 ) + itk_python_add_test( + NAME PythonGeometryProtocolTest + COMMAND + ${CMAKE_CURRENT_SOURCE_DIR}/geometry_protocol.py + ) itk_python_add_test( NAME PythonExtrasTest COMMAND diff --git a/Wrapping/Generators/Python/Tests/geometry_protocol.py b/Wrapping/Generators/Python/Tests/geometry_protocol.py new file mode 100644 index 00000000000..e1f87319fbc --- /dev/null +++ b/Wrapping/Generators/Python/Tests/geometry_protocol.py @@ -0,0 +1,197 @@ +# ========================================================================== +# +# 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. +# +# ========================================================================== + +"""Order-explicit spatial key protocol: the NumPy (z,y,x) frame is a full +isomorphic image of the ITK (x,y,z) frame under phi = flip (P.M.P for +matrices), verified against the C++ geometry oracles in both domains.""" + +import warnings + +import numpy as np + +import itk + +rng = np.random.default_rng(20260726) + + +def rot(axis, deg): + c, s = np.cos(np.radians(deg)), np.sin(np.radians(deg)) + if axis == "z": + return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]]) + if axis == "y": + return np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]]) + return np.array([[1, 0, 0], [0, c, -s], [0, s, c]]) + + +def phi_v(v): + return np.flip(np.asarray(v, dtype=float)) + + +def phi_m(m): + return np.flip(np.asarray(m, dtype=float), axis=None) + + +D_XYZ = rot("z", 30) @ rot("x", 20) +E_XYZ = rot("y", -40) @ rot("z", 10) +TOL = 1e-12 + + +def make_image(start_index): + image = itk.Image[itk.F, 3].New() + region = itk.ImageRegion[3]() + region.SetIndex(start_index) + region.SetSize([7, 9, 11]) + image.SetRegions(region) + image.Allocate() + image.SetSpacing([1.0, 2.0, 3.0]) + image.SetOrigin([5.0, 10.0, 15.0]) + image.SetDirection(itk.matrix_from_array(D_XYZ)) + return image + + +def continuous_index_oracle(image, ci_xyz): + ci = itk.ContinuousIndex[itk.D, 3]() + for k in range(3): + ci.SetElement(k, float(ci_xyz[k])) + return np.asarray(image.TransformContinuousIndexToPhysicalPoint(ci)) + + +for start_index in ([0, 0, 0], [4, -3, 10]): + image = make_image(start_index) + + o_xyz = image["origin_xyz"] + s_xyz = image["spacing_xyz"] + d_xyz = image["direction_xyz"] + i_xyz = image["index_xyz"] + n_xyz = image["size_xyz"] + o_zyx = image["origin_zyx"] + s_zyx = image["spacing_zyx"] + d_zyx = image["direction_zyx"] + i_zyx = image["index_zyx"] + n_zyx = image["size_zyx"] + + # xyz keys match member functions; zyx keys are phi of xyz + assert np.allclose(o_xyz, np.asarray(image.GetOrigin())) + assert np.allclose(s_xyz, np.asarray(image.GetSpacing())) + assert np.allclose(d_xyz, itk.array_from_matrix(image.GetDirection())) + assert np.array_equal( + i_xyz, np.asarray(image.GetLargestPossibleRegion().GetIndex()) + ) + assert np.array_equal(n_xyz, np.asarray(image.GetLargestPossibleRegion().GetSize())) + for a, b in ((o_zyx, o_xyz), (s_zyx, s_xyz), (i_zyx, i_xyz), (n_zyx, n_xyz)): + assert np.allclose(a, phi_v(b)) + assert np.allclose(d_zyx, phi_m(d_xyz)) + assert np.asarray(itk.array_view_from_image(image)).shape == tuple( + int(v) for v in n_zyx + ) + + # phi is a homomorphism: composition and inversion close in either domain + assert np.allclose(phi_m(D_XYZ @ E_XYZ), phi_m(D_XYZ) @ phi_m(E_XYZ)) + assert np.allclose(phi_m(np.linalg.inv(D_XYZ)), np.linalg.inv(phi_m(D_XYZ))) + + for _ in range(100): + ci = rng.uniform(-5, 15, 3) + p_oracle = continuous_index_oracle(image, ci) + # forward computed in each domain vs oracle + p_itk = o_xyz + d_xyz @ (s_xyz * ci) + p_np = o_zyx + d_zyx @ (s_zyx * phi_v(ci)) + assert np.abs(p_itk - p_oracle).max() < TOL + assert np.abs(p_np - phi_v(p_oracle)).max() < TOL + + p = rng.uniform(-50, 80, 3) + ci_oracle = np.asarray( + image.TransformPhysicalPointToContinuousIndex([float(v) for v in p]) + ) + # inverse computed in each domain vs oracle + ci_itk = np.linalg.inv(d_xyz) @ (p - o_xyz) / s_xyz + ci_np = np.linalg.inv(d_zyx) @ (phi_v(p) - o_zyx) / s_zyx + assert np.abs(ci_itk - ci_oracle).max() < TOL + assert np.abs(ci_np - phi_v(ci_oracle)).max() < TOL + + # region completeness: physical location of array voxel [0, 0, 0] + p0_oracle = np.asarray(image.TransformIndexToPhysicalPoint(start_index)) + assert np.allclose(o_xyz + d_xyz @ (s_xyz * i_xyz), p0_oracle) + assert np.allclose(phi_v(o_zyx + d_zyx @ (s_zyx * i_zyx)), p0_oracle) + + # setters round-trip in both conventions + other = make_image(start_index) + other["origin_zyx"] = o_zyx + other["spacing_zyx"] = s_zyx + other["direction_zyx"] = d_zyx + assert np.allclose(np.asarray(other.GetOrigin()), o_xyz) + assert np.allclose(np.asarray(other.GetSpacing()), s_xyz) + assert np.allclose(itk.array_from_matrix(other.GetDirection()), d_xyz) + other["origin_xyz"] = o_xyz + 1.0 + other["direction_xyz"] = E_XYZ + assert np.allclose(np.asarray(other.GetOrigin()), o_xyz + 1.0) + assert np.allclose(itk.array_from_matrix(other.GetDirection()), E_XYZ) + + # index/size writes go through SetRegions + resized = make_image(start_index) + resized["index_xyz"] = [1, 2, 3] + assert np.array_equal( + np.asarray(resized.GetLargestPossibleRegion().GetIndex()), [1, 2, 3] + ) + resized["index_zyx"] = phi_v([1, 2, 3]).astype(int) + assert np.array_equal(resized["index_xyz"], [1, 2, 3]) + resized["size_zyx"] = [13, 9, 7] + resized.Allocate() + assert np.array_equal( + np.asarray(resized.GetLargestPossibleRegion().GetSize()), [7, 9, 13] + ) + assert np.asarray(itk.array_view_from_image(resized)).shape == (13, 9, 7) + + # bare keys warn but keep their historical (z,y,x) behavior + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + assert np.allclose(image["origin"], o_zyx) + assert np.allclose(image["spacing"], s_zyx) + assert np.allclose(image["direction"], d_zyx) + image["origin"] = o_zyx + deprecations = [w for w in caught if issubclass(w.category, DeprecationWarning)] + assert len(deprecations) == 4 + assert "6706" in str(deprecations[0].message) + + # ITK_FUTURE_LEGACY_REMOVE escalates bare-key access to FutureWarning + import itkConfig + + saved = itkConfig.FutureLegacyRemove + try: + itkConfig.FutureLegacyRemove = True + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + image["origin"] + assert len(caught) == 1 and issubclass(caught[0].category, FutureWarning) + finally: + itkConfig.FutureLegacyRemove = saved + + # explicit keys never warn + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + image["origin_zyx"] + image["direction_xyz"] + image["origin_zyx"] = o_zyx + assert not [w for w in caught if issubclass(w.category, DeprecationWarning)] + + # all protocol keys are advertised + advertised = set(image.keys()) + for suffix in ("xyz", "zyx"): + for base in ("origin", "spacing", "direction", "index", "size"): + assert f"{base}_{suffix}" in advertised + +print("Test finished.") diff --git a/Wrapping/Generators/Python/itkConfig.template.in.py b/Wrapping/Generators/Python/itkConfig.template.in.py index ad83217162f..b6cd5d98331 100644 --- a/Wrapping/Generators/Python/itkConfig.template.in.py +++ b/Wrapping/Generators/Python/itkConfig.template.in.py @@ -78,6 +78,9 @@ def _strtobool(val: str) -> bool: "ITK_PYTHON_DEFAULTFACTORYLOADING", "True" ) NotInPlace: bool = _get_environment_boolean("ITK_PYTHON_NOTINPLACE", "False") +FutureLegacyRemove: bool = _get_environment_boolean( + "ITK_PYTHON_FUTURE_LEGACY_REMOVE", "@ITK_WRAP_PYTHON_FUTURE_LEGACY_REMOVE@" +) from os import environ as _environ