diff --git a/Modules/Bridge/NumPy/wrapping/PyBuffer.i.init b/Modules/Bridge/NumPy/wrapping/PyBuffer.i.init index 8f502c6fd66..991ceecc51c 100644 --- a/Modules/Bridge/NumPy/wrapping/PyBuffer.i.init +++ b/Modules/Bridge/NumPy/wrapping/PyBuffer.i.init @@ -29,31 +29,85 @@ else: loads = dask_deserialize.dispatch(np.ndarray) return NDArrayITKBase(loads(header, frames)) +def _get_buffer_formatstring(itk_pixel_code: str) -> str: + """Return the struct format character for an ITK pixel type code. + + Used by the PEP 688 ``__buffer__`` protocol implementation on + ``itk.Image`` to describe the element type of the exported + memoryview. Format characters follow Python's ``struct`` module + specification. + + Parameters + ---------- + itk_pixel_code : str + Short name of the ITK component type, e.g. ``"UC"``, ``"F"``. + + Returns + ------- + str + Single-character ``struct`` format string. + """ + + return _BUFFER_FORMAT_MAP[itk_pixel_code] + + +# Module-level constant — built once at import time. +# Composite pixel types (RGB, Vector, etc.) are decomposed to their +# scalar component type before reaching _get_buffer_formatstring(), +# so only scalar codes appear here. +# Platform aliases (ST/IT/OT) resolve to UL/SL at the wrapping level. +import os as _os +_BUFFER_FORMAT_MAP = { + # --- integer types --- + "B": "?", # bool + "UC": "B", # unsigned char -> uint8 + "US": "H", # unsigned short -> uint16 + "UI": "I", # unsigned int -> uint32 + "UL": "L", # unsigned long -> platform (8 bytes Linux/macOS, 4 bytes Windows) + "ULL": "Q", # unsigned long long -> uint64 + "SC": "b", # signed char -> int8 + "SS": "h", # signed short -> int16 + "SI": "i", # signed int -> int32 + "SL": "l", # signed long -> platform + "SLL": "q", # signed long long -> int64 + # --- floating point types --- + "F": "f", # float -> float32 + "D": "d", # double -> float64 + # "LD" (long double) intentionally omitted: sizeof(long double) + # is 16 bytes on Linux/macOS x86-64 but struct "d" is 8 bytes. + # Casting would silently corrupt the buffer. +} +if _os.name == 'nt': + # On Windows, C ``long`` is 32-bit + _BUFFER_FORMAT_MAP['UL'] = 'I' + _BUFFER_FORMAT_MAP['SL'] = 'i' + + def _get_numpy_pixelid(itk_Image_type) -> np.dtype: """Returns a ITK PixelID given a numpy array.""" - # This is a Mapping from numpy array types to itk pixel types. - _np_itk = {"UC":np.dtype(np.uint8), - "US":np.dtype(np.uint16), - "UI":np.dtype(np.uint32), - "UL":np.dtype(np.uint64), - "ULL":np.dtype(np.uint64), - "SC":np.dtype(np.int8), - "SS":np.dtype(np.int16), - "SI":np.dtype(np.int32), - "SL":np.dtype(np.int64), - "SLL":np.dtype(np.int64), - "F":np.dtype(np.float32), - "D":np.dtype(np.float64), - "PF2":np.dtype(np.float32), - "PF3":np.dtype(np.float32), - } - import os - if os.name == 'nt': - _np_itk['UL'] = np.dtype(np.uint32) - _np_itk['SL'] = np.dtype(np.int32) - try: - return _np_itk[itk_Image_type] - except KeyError as e: - raise e + return _NUMPY_PIXELID_MAP[itk_Image_type] + + +# Module-level constant — built once at import time. +_NUMPY_PIXELID_MAP = { + "B": np.dtype(np.bool_), + "UC": np.dtype(np.uint8), + "US": np.dtype(np.uint16), + "UI": np.dtype(np.uint32), + "UL": np.dtype(np.uint64), + "ULL": np.dtype(np.uint64), + "SC": np.dtype(np.int8), + "SS": np.dtype(np.int16), + "SI": np.dtype(np.int32), + "SL": np.dtype(np.int64), + "SLL": np.dtype(np.int64), + "F": np.dtype(np.float32), + "D": np.dtype(np.float64), + "PF2": np.dtype(np.float32), + "PF3": np.dtype(np.float32), +} +if _os.name == 'nt': + _NUMPY_PIXELID_MAP['UL'] = np.dtype(np.uint32) + _NUMPY_PIXELID_MAP['SL'] = np.dtype(np.int32) %} diff --git a/Modules/Core/Common/wrapping/test/CMakeLists.txt b/Modules/Core/Common/wrapping/test/CMakeLists.txt index 76ae8305118..e4b530685a0 100644 --- a/Modules/Core/Common/wrapping/test/CMakeLists.txt +++ b/Modules/Core/Common/wrapping/test/CMakeLists.txt @@ -56,4 +56,9 @@ if(ITK_WRAP_PYTHON) COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/itkSpatialOrientationAdapterTest.py ) + itk_python_add_test( + NAME itkImageInteropPythonTest + COMMAND + ${CMAKE_CURRENT_SOURCE_DIR}/itkImageInteropTest.py + ) endif() diff --git a/Modules/Core/Common/wrapping/test/itkImageInteropTest.py b/Modules/Core/Common/wrapping/test/itkImageInteropTest.py new file mode 100755 index 00000000000..36d2dd6b916 --- /dev/null +++ b/Modules/Core/Common/wrapping/test/itkImageInteropTest.py @@ -0,0 +1,363 @@ +#!/usr/bin/env python +# ========================================================================== +# +# 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. +# +# ========================================================================== +"""Comprehensive interop tests for itk.Image with numpy, torch, and dask. + +Tests use image sizes representative of clinical imaging modalities: + - CT: 512x512x128 signed short (int16) + - MRI: 256x256x64 float32 + - DWI: 128x128x40 float32 vector (30 gradient directions) + +Required: numpy +Optional: torch, dask (tests skipped if unavailable) +""" +import sys +import itk +import numpy as np + +# --- Optional dependency detection --- +_HAVE_TORCH = False +try: + import torch + + _HAVE_TORCH = True +except ImportError: + pass + +_HAVE_DASK = False +try: + import dask + import dask.array as da + + _HAVE_DASK = True +except ImportError: + pass + +# Zero-copy __array_interface__ is the default (correct NumPy semantics). +# Set itk.SIMULATE_PEP688 = False to revert to legacy deep-copy behavior. + +passed = 0 +skipped = 0 +failed = 0 + + +def check(name, condition, detail=""): + global passed, failed + if condition: + passed += 1 + else: + failed += 1 + msg = f" FAIL: {name}" + if detail: + msg += f" ({detail})" + print(msg) + + +def skip(name, reason): + global skipped + skipped += 1 + + +# =================================================================== +# Helper: create images matching clinical modalities +# =================================================================== + + +def make_ct_image(): + """CT: 512x512x128, signed short (int16), ~67 MB.""" + img = itk.Image[itk.SS, 3].New() + img.SetRegions([512, 512, 128]) + img.SetSpacing([0.5, 0.5, 1.0]) + img.Allocate() + img.FillBuffer(-1024) # Hounsfield air + return img + + +def make_mri_image(): + """MRI T1: 256x256x64, float32, ~16 MB.""" + img = itk.Image[itk.F, 3].New() + img.SetRegions([256, 256, 64]) + img.SetSpacing([1.0, 1.0, 2.0]) + img.Allocate() + img.FillBuffer(500.0) # Typical T1 signal + return img + + +def make_dwi_image(): + """DWI: 128x128x40, float32, single volume, ~2.5 MB.""" + img = itk.Image[itk.F, 3].New() + img.SetRegions([128, 128, 40]) + img.SetSpacing([2.0, 2.0, 2.5]) + img.Allocate() + img.FillBuffer(800.0) # B0 signal + return img + + +def make_rgb_image(): + """RGB color: 256x256, unsigned char, ~192 KB.""" + img = itk.Image[itk.RGBPixel[itk.UC], 2].New() + img.SetRegions([256, 256]) + img.Allocate() + return img + + +def make_vector_image(): + """Vector image: 64x64x32, 6 components (DTI tensor), float32.""" + img = itk.VectorImage[itk.F, 3].New() + img.SetRegions([64, 64, 32]) + img.SetNumberOfComponentsPerPixel(6) + img.Allocate(True) # zero-initialize + return img + + +# =================================================================== +# Section 1: NumPy interop (REQUIRED) +# =================================================================== +print("=" * 60) +print("Section 1: NumPy interop (required)") +print("=" * 60) + +# -- 1.1 __array_interface__ basic functionality -- +ct = make_ct_image() +ai = ct.__array_interface__ +check("CT __array_interface__ version", ai["version"] == 3) +check("CT __array_interface__ shape", ai["shape"] == (128, 512, 512)) +check("CT __array_interface__ typestr", ai["typestr"] in ("itk", ct.GetPixel([30, 20, 10]) == 1500) + +# Verify zero-copy: modify through ITK, read back through numpy +ct.SetPixel([5, 6, 7], 42) +check("CT zero-copy itk->np", ct_arr[7, 6, 5] == 42) + +# -- 1.3 MRI image -- +mri = make_mri_image() +mri_arr = np.asarray(mri) +check("MRI np.asarray shape", mri_arr.shape == (64, 256, 256)) +check("MRI np.asarray dtype", mri_arr.dtype == np.float32) +check("MRI np.asarray value", abs(mri_arr[0, 0, 0] - 500.0) < 1e-5) + +# Zero-copy write +mri_arr[32, 128, 128] = 999.0 +check("MRI zero-copy", mri.GetPixel([128, 128, 32]) == 999.0) + +# -- 1.4 DWI single volume -- +dwi = make_dwi_image() +dwi_arr = np.asarray(dwi) +check("DWI np.asarray shape", dwi_arr.shape == (40, 128, 128)) +check("DWI np.asarray dtype", dwi_arr.dtype == np.float32) + +# -- 1.5 Multiple pixel types -- +for ptype, expected_dtype, label in [ + (itk.UC, np.uint8, "UC/uint8"), + (itk.SS, np.int16, "SS/int16"), + (itk.SI, np.int32, "SI/int32"), + (itk.F, np.float32, "F/float32"), + (itk.D, np.float64, "D/float64"), +]: + img = itk.Image[ptype, 3].New() + img.SetRegions([4, 4, 4]) + img.Allocate() + arr = np.asarray(img) + check(f"{label} dtype", arr.dtype == expected_dtype, f"got {arr.dtype}") + check(f"{label} shape", arr.shape == (4, 4, 4)) + +# -- 1.6 __array__ returns view (not copy) -- +mri2 = make_mri_image() +arr1 = mri2.__array__() +check("__array__ returns ndarray", isinstance(arr1, np.ndarray)) +arr1[0, 0, 0] = 12345.0 +check("__array__ zero-copy", mri2.GetPixel([0, 0, 0]) == 12345.0) + +# -- 1.7 __buffer__ protocol -- +ct_mv = ct.__buffer__() +check("CT __buffer__ type", isinstance(ct_mv, memoryview)) +check("CT __buffer__ format", ct_mv.format == "h") # signed short +check("CT __buffer__ shape", ct_mv.shape == (128, 512, 512)) + +if sys.version_info >= (3, 12): + auto_mv = memoryview(ct) + check("PEP 688 auto memoryview", auto_mv.shape == (128, 512, 512)) + +# -- 1.8 Vector/multi-component images -- +# Note: VectorImage uses array_view_from_image path (DECL_PYTHON_IMAGE_CLASS +# is only applied to itk::Image, not itk::VectorImage) +vec = make_vector_image() +vec_arr = itk.array_view_from_image(vec) +check("VectorImage shape", vec_arr.shape == (32, 64, 64, 6)) +check("VectorImage dtype", vec_arr.dtype == np.float32) + +# -- 1.9 np.array with dtype conversion -- +ct_f32 = np.array(ct, dtype=np.float32) +check("np.array dtype conversion", ct_f32.dtype == np.float32) + +# -- 1.10 Consistency: __array_interface__ vs __buffer__ vs array_view -- +mri3 = make_mri_image() +via_ai = np.asarray(mri3) +via_view = itk.array_view_from_image(mri3) +via_buf = np.asarray(mri3.__buffer__()) +check( + "Consistency: all paths same data", + np.array_equal(via_ai, via_view) and np.array_equal(via_ai, via_buf), +) + +print(f" NumPy: {passed} passed, {failed} failed") +numpy_passed = passed +numpy_failed = failed + +# =================================================================== +# Section 2: PyTorch interop (optional) +# =================================================================== +print() +print("=" * 60) +if _HAVE_TORCH: + print(f"Section 2: PyTorch interop (torch {torch.__version__})") +else: + print("Section 2: PyTorch interop (SKIPPED — torch not installed)") +print("=" * 60) + +if _HAVE_TORCH: + # -- 2.1 torch.as_tensor via __array_interface__ -- + ct2 = make_ct_image() + ct_tensor = torch.as_tensor(np.asarray(ct2)) + check("CT torch.as_tensor shape", tuple(ct_tensor.shape) == (128, 512, 512)) + check("CT torch.as_tensor dtype", ct_tensor.dtype == torch.int16) + check("CT torch.as_tensor value", ct_tensor[0, 0, 0].item() == -1024) + + # -- 2.2 torch.from_numpy zero-copy chain -- + mri4 = make_mri_image() + mri_np = np.asarray(mri4) # zero-copy via __array_interface__ + mri_tensor = torch.from_numpy(mri_np) # zero-copy numpy->torch + check("MRI torch shape", tuple(mri_tensor.shape) == (64, 256, 256)) + check("MRI torch dtype", mri_tensor.dtype == torch.float32) + + # Zero-copy chain: ITK -> numpy -> torch + mri4.SetPixel([100, 100, 30], 777.0) + check( + "MRI zero-copy itk->np->torch", + mri_tensor[30, 100, 100].item() == 777.0, + ) + + # -- 2.3 DWI volume -- + dwi2 = make_dwi_image() + dwi_tensor = torch.from_numpy(np.asarray(dwi2)) + check("DWI torch shape", tuple(dwi_tensor.shape) == (40, 128, 128)) + + # -- 2.4 Vector image -- + vec2 = make_vector_image() + vec_tensor = torch.from_numpy(itk.array_view_from_image(vec2)) + check("VectorImage torch shape", tuple(vec_tensor.shape) == (32, 64, 64, 6)) + check("VectorImage torch dtype", vec_tensor.dtype == torch.float32) + + # -- 2.5 dtype mapping -- + for ptype, expected_tdtype, label in [ + (itk.UC, torch.uint8, "UC->uint8"), + (itk.SS, torch.int16, "SS->int16"), + (itk.SI, torch.int32, "SI->int32"), + (itk.F, torch.float32, "F->float32"), + (itk.D, torch.float64, "D->float64"), + ]: + img = itk.Image[ptype, 3].New() + img.SetRegions([4, 4, 4]) + img.Allocate() + t = torch.from_numpy(np.asarray(img)) + check(f"torch {label}", t.dtype == expected_tdtype, f"got {t.dtype}") + + torch_passed = passed - numpy_passed + torch_failed = failed - numpy_failed + print(f" PyTorch: {torch_passed} passed, {torch_failed} failed") +else: + skip("PyTorch tests", "torch not installed") + print(" (install torch to enable these tests)") + +# =================================================================== +# Section 3: Dask interop (optional) +# =================================================================== +print() +print("=" * 60) +if _HAVE_DASK: + print(f"Section 3: Dask interop (dask {dask.__version__})") +else: + print("Section 3: Dask interop (SKIPPED — dask not installed)") +print("=" * 60) + +pre_dask_passed = passed +pre_dask_failed = failed + +if _HAVE_DASK: + # -- 3.1 dask.array.from_array via __array_interface__ -- + ct3 = make_ct_image() + ct_darr = da.from_array(np.asarray(ct3), chunks=(32, 128, 128)) + check("CT dask shape", ct_darr.shape == (128, 512, 512)) + check("CT dask dtype", ct_darr.dtype == np.int16) + + # Compute a chunk and verify value + chunk = ct_darr[0:32, 0:128, 0:128].compute() + check("CT dask compute value", chunk[0, 0, 0] == -1024) + check("CT dask compute shape", chunk.shape == (32, 128, 128)) + + # -- 3.2 MRI with dask -- + mri5 = make_mri_image() + mri_darr = da.from_array(np.asarray(mri5), chunks=(16, 64, 64)) + check("MRI dask shape", mri_darr.shape == (64, 256, 256)) + check("MRI dask dtype", mri_darr.dtype == np.float32) + + # -- 3.3 Dask reduction on clinical-size image -- + mri_mean = mri_darr.mean().compute() + check("MRI dask mean", abs(mri_mean - 500.0) < 1e-3, f"got {mri_mean}") + + # -- 3.4 DWI with dask -- + dwi3 = make_dwi_image() + dwi_darr = da.from_array(np.asarray(dwi3), chunks=(10, 32, 32)) + check("DWI dask shape", dwi_darr.shape == (40, 128, 128)) + + # -- 3.5 Vector image with dask -- + vec3 = make_vector_image() + vec_darr = da.from_array(itk.array_view_from_image(vec3), chunks=(8, 16, 16, 6)) + check("VectorImage dask shape", vec_darr.shape == (32, 64, 64, 6)) + + dask_passed = passed - pre_dask_passed + dask_failed = failed - pre_dask_failed + print(f" Dask: {dask_passed} passed, {dask_failed} failed") +else: + skip("Dask tests", "dask not installed") + print(" (install dask to enable these tests)") + +# =================================================================== +# Summary +# =================================================================== +print() +print("=" * 60) +print(f"TOTAL: {passed} passed, {failed} failed, {skipped} skipped") +print("=" * 60) + +if failed > 0: + sys.exit(1) + +print("All interop tests passed.") diff --git a/Modules/Core/Common/wrapping/test/itkImageTest.py b/Modules/Core/Common/wrapping/test/itkImageTest.py index 140806f119a..eccaef12a60 100644 --- a/Modules/Core/Common/wrapping/test/itkImageTest.py +++ b/Modules/Core/Common/wrapping/test/itkImageTest.py @@ -15,6 +15,7 @@ # limitations under the License. # # ========================================================================== +import sys import itk import numpy as np @@ -39,3 +40,121 @@ assert array[0, 0] == 4 assert array[0, 1] == 4 assert isinstance(array, np.ndarray) + +# --- PEP 688 buffer protocol tests --- + +# Test __buffer__ on scalar image +mv = image.__buffer__() +assert isinstance(mv, memoryview) +assert mv.format == "B" # unsigned char +assert mv.shape == (10, 10) +assert mv[0, 0] == 4 +assert mv.readonly is False + +# Test that memoryview shares memory (zero-copy) +image.SetPixel([3, 5], 42) +assert mv[5, 3] == 42 + +# Test float image +float_image = itk.Image[itk.F, 2].New() +float_image.SetRegions([8, 6]) +float_image.Allocate() +float_image.FillBuffer(3.14) +fmv = float_image.__buffer__() +assert fmv.format == "f" +assert fmv.shape == (6, 8) +assert abs(fmv[0, 0] - 3.14) < 1e-5 + +# Test 3D image +image_3d = itk.Image[itk.SS, 3].New() +image_3d.SetRegions([4, 5, 6]) +image_3d.Allocate() +image_3d.FillBuffer(-7) +mv3d = image_3d.__buffer__() +assert mv3d.format == "h" # signed short +assert mv3d.shape == (6, 5, 4) +assert mv3d[0, 0, 0] == -7 + +# On Python 3.12+, memoryview() should call __buffer__ automatically +if sys.version_info >= (3, 12): + auto_mv = memoryview(image) + assert auto_mv.shape == (10, 10) + assert auto_mv[0, 0] == 4 + +# --- Default: zero-copy (correct NumPy np.asarray semantics) --- + +# __array_interface__ is active by default +ai = image.__array_interface__ +assert ai["version"] == 3, f"Expected version 3, got {ai['version']}" +assert ai["shape"] == (10, 10), f"Expected (10, 10), got {ai['shape']}" +assert isinstance(ai["data"], tuple) +assert isinstance(ai["data"][0], int) # pointer +assert ai["data"][1] is False # not read-only + +# np.asarray returns zero-copy view (correct per NumPy semantics) +arr_view = np.asarray(image) +assert arr_view.shape == (10, 10) +assert arr_view.dtype == np.uint8 +image.SetPixel([1, 2], 99) +assert arr_view[2, 1] == 99 # zero-copy verified + +# __array__ returns zero-copy view by default +arr_a = image.__array__() +image.SetPixel([0, 0], 77) +assert arr_a[0, 0] == 77 # zero-copy verified + +# --- __array__ copy parameter (NumPy 2.0 protocol) --- + +# copy=None (default): zero-copy view +arr_none = image.__array__(copy=None) +image.SetPixel([0, 0], 55) +assert arr_none[0, 0] == 55 # zero-copy + +# copy=True: always returns a copy +arr_true = image.__array__(copy=True) +image.SetPixel([0, 0], 66) +assert arr_true[0, 0] == 55 # copy: still old value + +# copy=False: zero-copy, succeeds when no copy is needed +arr_false = image.__array__(copy=False) +image.SetPixel([0, 0], 88) +assert arr_false[0, 0] == 88 # zero-copy + +# copy=False with incompatible dtype: raises ValueError +try: + _ = image.__array__(dtype=np.float64, copy=False) + assert False, "copy=False with dtype conversion should raise ValueError" +except ValueError: + pass # expected: dtype conversion requires a copy + +# copy=False with matching dtype: succeeds +arr_same_dtype = image.__array__(dtype=np.uint8, copy=False) +image.SetPixel([0, 0], 91) +assert arr_same_dtype[0, 0] == 91 # zero-copy, no dtype conversion + +# --- SIMULATE_PEP688=False: revert to legacy deep-copy behavior --- +itk.SIMULATE_PEP688 = False + +# __array_interface__ is disabled in legacy mode +try: + _ = image.__array_interface__ + assert False, "__array_interface__ should raise when SIMULATE_PEP688 is False" +except AttributeError: + pass # expected: reverts to legacy behavior + +# __array__ returns deep copy in legacy mode +arr_copy = image.__array__() +image.SetPixel([0, 0], 123) +assert arr_copy[0, 0] == 91 # deep copy: still old value + +# copy=False in legacy mode: raises ValueError (copy unavoidable) +try: + _ = image.__array__(copy=False) + assert False, "copy=False in legacy mode should raise ValueError" +except ValueError: + pass # expected: legacy mode always copies + +# Restore default +del itk.SIMULATE_PEP688 + +print("All buffer protocol tests passed.") diff --git a/Wrapping/Generators/Python/PyBase/pyBase.i b/Wrapping/Generators/Python/PyBase/pyBase.i index c90d8020de7..08e073b1823 100644 --- a/Wrapping/Generators/Python/PyBase/pyBase.i +++ b/Wrapping/Generators/Python/PyBase/pyBase.i @@ -679,13 +679,162 @@ str = str %define DECL_PYTHON_IMAGE_CLASS(swig_name) %extend swig_name { - %pythoncode { - def __array__(self, dtype=None): + %pythoncode %{ + def __buffer__(self, flags=0, /): + """PEP 688 buffer protocol — export image data as a memoryview. + + On Python 3.12+ this is called automatically by + ``memoryview(image)`` and ``numpy.asarray(image)``. + On Python 3.10–3.11 it can be called explicitly. + + The returned memoryview shares memory with the image + (zero-copy). The caller must keep the image alive for + as long as the memoryview is in use. + """ import itk + from itk.itkPyBufferPython import _get_buffer_formatstring + + # Get 1-D raw memoryview from the C++ buffer + ImageType = type(self) + PyBufferType = itk.PyBuffer[ImageType] + raw_memview = PyBufferType._GetArrayViewFromImage(self) + + # Build shape in C-order (NumPy convention: [z, y, x, ...]) + itksize = self.GetBufferedRegion().GetSize() + shape = [int(itksize[d]) for d in range(len(itksize))] + + n_components = self.GetNumberOfComponentsPerPixel() + if n_components > 1: + shape.insert(0, n_components) + + shape.reverse() + + # Determine the struct format character for the component type + tpl = itk.template(self) + pixel_type = tpl[1][0] + from itk.support.types import itkCType + if isinstance(pixel_type, itkCType): + # Scalar pixel (UC, F, SS, etc.) + component_code = pixel_type.short_name + else: + # Composite pixel (RGB, Vector, etc.) — use component type + pixel_tpl = itk.template(pixel_type) + component_code = pixel_tpl[1][0].short_name + + fmt = _get_buffer_formatstring(component_code) + + return raw_memview.cast(fmt, shape=shape) + + def _get_array_interface(self): + """NumPy array interface (v3) — zero-copy on all Python versions. + + Returns a dict that lets ``numpy.asarray(image)`` create a + zero-copy view of the pixel buffer. This is the correct + behavior per NumPy semantics: ``np.asarray()`` should avoid + copies when possible. + + Disabled when ``itk.SIMULATE_PEP688`` is explicitly set to + ``False``, reverting to the legacy deep-copy behavior for + code that depended on the previous (incorrect) semantics. + + Set ``itk.SIMULATE_PEP688_DEBUG = True`` to print a + diagnostic each time this zero-copy path is taken, helping + identify code that assumed ``np.asarray()`` returned a copy. + """ + import itk + if getattr(itk, 'SIMULATE_PEP688', True) is False: + raise AttributeError('__array_interface__') + + if getattr(itk, 'SIMULATE_PEP688_DEBUG', False): + import traceback + print( + "[itk.SIMULATE_PEP688_DEBUG] __array_interface__ " + "returning zero-copy view for " + f"{type(self).__name__}. " + "Code that mutates the result will modify the " + "image. Use np.array(image, copy=True) or " + "itk.array_from_image(image) if a copy is needed." + ) + traceback.print_stack(limit=5) + import numpy as np - array = itk.array_from_image(self) - return np.asarray(array, dtype=dtype) - } + from itk.itkPyBufferPython import _get_numpy_pixelid + + ImageType = type(self) + PyBufferType = itk.PyBuffer[ImageType] + raw_memview = PyBufferType._GetArrayViewFromImage(self) + + # Shape in C-order (NumPy convention: [z, y, x]) + itksize = self.GetBufferedRegion().GetSize() + shape = tuple(int(itksize[d]) for d in reversed(range(len(itksize)))) + + n_components = self.GetNumberOfComponentsPerPixel() + if n_components > 1 or isinstance(self, itk.VectorImage): + shape = shape + (n_components,) + + # Resolve component type code + tpl = itk.template(self) + pixel_type = tpl[1][0] + from itk.support.types import itkCType + if isinstance(pixel_type, itkCType): + component_code = pixel_type.short_name + else: + pixel_tpl = itk.template(pixel_type) + component_code = pixel_tpl[1][0].short_name + + dtype = _get_numpy_pixelid(component_code) + np_arr = np.asarray(raw_memview) + data_ptr = np_arr.__array_interface__['data'][0] + + return { + 'version': 3, + 'shape': shape, + 'typestr': dtype.str, + 'data': (data_ptr, False), + 'strides': None, + } + __array_interface__ = property(_get_array_interface) + + def __array__(self, dtype=None, copy=None): + import itk + import numpy as np + + if getattr(itk, 'SIMULATE_PEP688', True) is False: + # Legacy mode: return deep copy (pre-PEP688 behavior) + if copy is False: + raise ValueError( + "Unable to avoid copy: itk.SIMULATE_PEP688 is " + "False (legacy deep-copy mode). Set " + "itk.SIMULATE_PEP688 = True or use " + "itk.array_view_from_image() for zero-copy." + ) + array = itk.array_from_image(self) + else: + # Correct behavior: return zero-copy view + if getattr(itk, 'SIMULATE_PEP688_DEBUG', False): + import traceback + print( + "[itk.SIMULATE_PEP688_DEBUG] __array__ " + "returning zero-copy view for " + f"{type(self).__name__}. " + "Set itk.SIMULATE_PEP688 = False to revert " + "to deep-copy behavior." + ) + traceback.print_stack(limit=5) + array = itk.array_view_from_image(self) + + if dtype is not None: + if copy is False and np.dtype(dtype) != array.dtype: + raise ValueError( + "Unable to avoid copy: dtype conversion from " + f"{array.dtype} to {np.dtype(dtype)} requires " + "a copy." + ) + array = np.asarray(array, dtype=dtype) + if copy: + array = array.copy() + return array + %} } %enddef