Merge pull request #78 from inhuszar/fix/numpy2x

NumPy 2.x compatibility fixes for image interpolation and I/O

Author: jnolan14

setup.py

@@ -1,5 +1,6 @@
 #!/usr/bin/env python
 
+import os
 import re
 import pathlib
 import platform
@@ -16,12 +17,16 @@ class bdist_wheel_abi3(bdist_wheel):  # noqa: D101
     def get_tag(self):  # noqa: D102
         python, abi, plat = super().get_tag()
 
+        # Only meaningful when we are *actually* building an abi3 wheel.
+        # We target cp311 as the minimum CPython that supports the features
+        # this project relies on for limited API mode.
         if python.startswith("cp"):
             return "cp311", "abi3", plat
 
         return python, abi, plat
 
 
+
 requirements = [
     'numpy',
     'scipy',
@@ -50,8 +55,10 @@ def get_tag(self):  # noqa: D102
 )
 macros = []
 setup_opts = {}
-if sys.version_info.minor >= 11 and platform.python_implementation() == "CPython":
-    # Can create an abi3 wheel (typed memoryviews first available in 3.11)!
+
+ABI3 = (os.environ.get("SURFA_ABI3", "0") == "1")
+
+if ABI3 and sys.version_info >= (3, 11) and platform.python_implementation() == "CPython":
     ext_opts["define_macros"].append(("Py_LIMITED_API", "0x030B0000"))
     ext_opts["py_limited_api"] = True
     setup_opts["cmdclass"] = {"bdist_wheel": bdist_wheel_abi3}
@@ -63,6 +70,7 @@ def get_tag(self):  # noqa: D102
 # since we interface the c stuff with numpy, it's another hard
 # requirement at build-time
 import numpy as np
+
 include_dirs = [np.get_include()]
 
 # extract the current version
@@ -97,7 +105,7 @@ def get_tag(self):  # noqa: D102
     packages=packages,
     ext_modules=extensions,
     include_dirs=include_dirs,
-    package_data={'': ['*.pyx'], '': ['*.h']},
+    package_data={'': ['*.pyx', '*.h']},
     install_requires=requirements,
     classifiers=[
         'Development Status :: 3 - Alpha',

surfa/image/framed.py

@@ -509,8 +509,18 @@ def reorient(self, orientation, copy=True, inplace=False):
         for i in range(self.basedim):
             if world_axes_src[i] != world_axes_trg[i]:
                 data = np.swapaxes(data, world_axes_src[i], world_axes_trg[i])
-                swapped_axis_idx = np.where(world_axes_src == world_axes_trg[i])
-                world_axes_src[swapped_axis_idx], world_axes_src[i] = world_axes_src[i], world_axes_src[swapped_axis_idx]
+                # NumPy 2.x no longer supports mixed scalar/sequence tuple-assignment
+                # with np.where tuple indexing here, so resolve a scalar index first.
+                swapped_axis_idx = np.where(world_axes_src == world_axes_trg[i])[0]
+                if swapped_axis_idx.size == 0:
+                    raise RuntimeError(
+                        "Failed to resolve source axis during reorientation."
+                    )
+                swapped_axis = int(swapped_axis_idx[0])
+                world_axes_src[swapped_axis], world_axes_src[i] = (
+                    world_axes_src[i],
+                    world_axes_src[swapped_axis],
+                )
 
         # align directions
         dot_products = np.sum(affine[:3, :3] * trg_matrix[:3, :3], axis=0)

surfa/image/interp.pyx

@@ -93,9 +93,10 @@ def interpolate(source, target_shape, method, affine=None, disp=None, fill=0):
     # ensure correct byteorder
     # TODO maybe this should be done at read-time?
     swap_byteorder = sys.byteorder == 'little' and '>' or '<'
-    if (source.dtype.byteorder == swap_byteorder):
-       source = source.byteswap().view(source.dtype.newbyteorder())
-    
+    if source.dtype.byteorder == swap_byteorder:
+        # NumPy 2 removed ndarray.newbyteorder(); use a dtype view instead.
+        source = source.byteswap().view(source.dtype.newbyteorder('='))
+
     # a few types aren't supported, so let's just convert to float and convert back if necessary
     unsupported_dtype = None
     if source.dtype in (np.bool_,):

surfa/io/framed.py

@@ -2,6 +2,13 @@
 import gzip
 import numpy as np
 
+# NumPy compatibility: prefer the numpy boolean scalar dtype object.
+# `np.bool_` exists in both NumPy 1.x and 2.x; keep a fallback for safety.
+if hasattr(np, 'bool_'):
+    _np_bool_dtype = np.bool_  # NumPy 1.x
+else:
+    _np_bool_dtype = np.bool   # NumPy 2.x
+
 from surfa import Volume
 from surfa import Slice
 from surfa import Overlay
@@ -407,7 +414,7 @@ def save(self, arr, filename, intent=FramedArrayIntents.mri):
 
             # determine supported dtype to save as (order here is very important)
             type_map = {
-                np.bool_: 0,
+                _np_bool_dtype: 0,
                 np.uint8: 0,
                 np.int32: 1,
                 np.floating: 3,
@@ -670,7 +677,7 @@ def save(self, arr, filename, intent=FramedArrayIntents.mri):
 
         # convert to a valid output type (for now this is only bool but there are probably more)
         type_map = {
-            np.bool_: np.uint8,
+            _np_bool_dtype: np.uint8,
         }
         dtype_id = next((i for dt, i in type_map.items() if np.issubdtype(arr.dtype, dt)), None)
         data = arr.data if dtype_id is None else arr.data.astype(dtype_id)

surfa/io/fsnifti1extension.py

@@ -1,5 +1,13 @@
 import numpy as np
 
+# NumPy compatibility: handle np.int_ removal in NumPy 2.0
+# In NumPy 1.x, np.int_ is an alias for the platform's native int type
+# In NumPy 2.x, np.int_ was removed; use np.integer for type checking
+if hasattr(np, 'int_'):
+    _np_int_types = (np.int_, np.integer)  # NumPy 1.x: check both for safety
+else:
+    _np_int_types = (np.integer,)          # NumPy 2.x: only np.integer exists
+
 from surfa.io import fsio
 from surfa.io import utils as iou
 from surfa.transform.geometry import ImageGeometry
@@ -566,7 +574,7 @@ def _from_framedimage(self, image):
             self.version = 1
             self.intent = image.metadata.get('intent', FramedArrayIntents.mri)
             self.dof = 1
-            if isinstance(self.intent, np.int_):
+            if isinstance(self.intent, _np_int_types):
                 self.intent = self.intent.item()  # convert numpy int to python int
 
             if self.intent == FramedArrayIntents.warpmap:

test/test_numpy_compat.py

@@ -0,0 +1,46 @@
+import sys
+
+import numpy as np
+import pytest
+
+
+sf = pytest.importorskip('surfa')
+
+
+def _non_native_float32(data):
+    byteorder = '>' if sys.byteorder == 'little' else '<'
+    return np.asarray(data, dtype=np.float32).astype(np.dtype(np.float32).newbyteorder(byteorder), copy=True)
+
+
+def _assert_machine_precision_equal(a, b):
+    assert a.shape == b.shape
+    assert a.dtype == b.dtype
+    assert np.issubdtype(a.dtype, np.floating)
+
+    scale = max(
+        float(np.max(np.abs(a))),
+        float(np.max(np.abs(b))),
+        1.0,
+    )
+    atol = np.finfo(a.dtype).eps * scale
+    assert np.allclose(a, b, rtol=0.0, atol=atol), f'max abs diff={np.max(np.abs(a - b))}, atol={atol}'
+
+
+@pytest.mark.parametrize('method', ['linear', 'nearest'])
+def test_resize_non_native_endian_matches_native(method):
+    rng = np.random.default_rng(123456)
+    data_native = rng.normal(loc=0.1, scale=1.7, size=(9, 7, 5, 2)).astype(np.float32)
+    data_non_native = _non_native_float32(data_native)
+
+    vol_native = sf.Volume(data_native)
+    vol_non_native = sf.Volume(data_non_native)
+
+    assert vol_non_native.framed_data.dtype.byteorder not in ('=', '|')
+
+    voxsize = np.array([1.3, 0.9, 1.7], dtype=np.float32)
+    resized_native = vol_native.resize(voxsize, method=method)
+    resized_non_native = vol_non_native.resize(voxsize, method=method)
+
+    assert resized_native.framed_data.shape == resized_non_native.framed_data.shape
+    assert resized_native.framed_data.dtype == resized_non_native.framed_data.dtype
+    _assert_machine_precision_equal(resized_native.framed_data, resized_non_native.framed_data)