Fix #475

eitprocessing/datahandling/__init__.py

@@ -3,11 +3,12 @@
 
 from typing_extensions import Self
 
+from eitprocessing.datahandling.mixins.arrays import NotAnArray
 from eitprocessing.datahandling.mixins.equality import Equivalence
 
 
 @dataclass(eq=False)
-class DataContainer(Equivalence):
+class DataContainer(Equivalence, NotAnArray):
     """Base class for data container classes."""
 
     def __bool__(self):

eitprocessing/datahandling/continuousdata.py

@@ -45,7 +45,7 @@ class ContinuousData(DataContainer, SelectByTime):
     parameters: dict[str, Any] = field(default_factory=dict, repr=False, metadata={"check_equivalence": True})
     derived_from: Any | list[Any] = field(default_factory=list, repr=False, compare=False)
     time: np.ndarray = field(kw_only=True, repr=False)
-    values: np.ndarray = field(kw_only=True, repr=False)
+    values: np.ndarray = field(kw_only=True, repr=False, metadata={"array_attribute": True})
     sample_frequency: float | None = field(kw_only=True, repr=False, metadata={"check_equivalence": True}, default=None)
 
     def __post_init__(self) -> None:

eitprocessing/datahandling/datacollection.py

@@ -6,6 +6,7 @@
 from eitprocessing.datahandling.continuousdata import ContinuousData
 from eitprocessing.datahandling.eitdata import EITData
 from eitprocessing.datahandling.intervaldata import IntervalData
+from eitprocessing.datahandling.mixins.arrays import NotAnArray
 from eitprocessing.datahandling.mixins.equality import Equivalence
 from eitprocessing.datahandling.mixins.slicing import HasTimeIndexer
 from eitprocessing.datahandling.sparsedata import SparseData
@@ -18,7 +19,7 @@
 V_classes = V.__constraints__
 
 
-class DataCollection(Equivalence, UserDict, HasTimeIndexer, Generic[V]):
+class DataCollection(Equivalence, UserDict, HasTimeIndexer, NotAnArray, Generic[V]):
     """A collection of a single type of data with unique labels.
 
     A DataCollection functions largely as a dictionary, but requires a data_type argument, which must be one of the data

eitprocessing/datahandling/eitdata.py

@@ -48,7 +48,7 @@ class is meant to hold data from (part of) a singular continuous measurement.
     label: str | None = field(default=None, compare=False, metadata={"check_equivalence": True})
     description: str = field(default="", compare=False, repr=False)
     name: str | None = field(default=None, compare=False, repr=False)
-    pixel_impedance: np.ndarray = field(repr=False, kw_only=True)
+    pixel_impedance: np.ndarray = field(repr=False, kw_only=True, metadata={"array_attribute": True})
     suppress_simulated_warning: InitVar[bool] = False
 
     def __post_init__(self, suppress_simulated_warning: bool) -> None:

eitprocessing/datahandling/intervaldata.py

@@ -54,7 +54,7 @@ class IntervalData(DataContainer, SelectByIndex, HasTimeIndexer):
     unit: str | None = field(metadata={"check_equivalence": True}, repr=False)
     category: str = field(metadata={"check_equivalence": True}, repr=False)
     intervals: list[Interval | tuple[float, float]] = field(repr=False)
-    values: list[Any] | None = field(repr=False, default=None)
+    values: list[Any] | None = field(repr=False, default=None, metadata={"array_attribute": True})
     parameters: dict[str, Any] = field(default_factory=dict, metadata={"check_equivalence": True}, repr=False)
     derived_from: list[Any] = field(default_factory=list, compare=False, repr=False)
     description: str = field(compare=False, default="", repr=False)

eitprocessing/datahandling/mixins/arrays.py

@@ -0,0 +1,60 @@
+from __future__ import annotations
+
+from dataclasses import fields, is_dataclass
+from typing import TYPE_CHECKING, Any
+
+if TYPE_CHECKING:
+    from collections.abc import Callable
+
+    import numpy as np
+    from typing_extensions import Never
+
+
+class NotAnArray:
+    """Mixin class that prevents numpy and scipy from converting an object to an array.
+
+    Objects in this package keep their numeric data in an attribute (e.g. `values`), rather than being arrays
+    themselves. Passing such an object directly to a numpy or scipy function is virtually always a mistake, but numpy
+    does not treat it as one: it falls back to the sequence protocol (`__len__`/`__getitem__`) or wraps the object in a
+    0-dimensional object array. For sliceable objects that returns a copy of the object for every index, which is
+    prohibitively slow; for other objects it silently produces an object array that does not contain the data.
+
+    This mixin closes the three routes numpy uses, so any attempt raises a `TypeError` explaining what to pass instead:
+
+    - `__array__` is used by `numpy.asarray()`/`numpy.array()`, and therefore by most of scipy, which converts its
+      input before doing anything else;
+    - `__array_ufunc__` is used by ufuncs, e.g. `numpy.sin()`, `numpy.add()` and `array + object`;
+    - `__array_function__` is used by the rest of the numpy API, e.g. `numpy.mean()` and `numpy.concatenate()`, which
+      dispatches before any conversion happens.
+
+    Regular Python behaviour (slicing, `len()`, iteration, comparison) is unaffected.
+
+    The error message points at the field holding the data, if there is one. Mark that field with
+    `field(metadata={"array_attribute": True})` to have it named.
+    """
+
+    @property
+    def _array_attribute(self) -> str | None:
+        """Name of the field holding the data, or None if no field is marked as such.
+
+        Only looked up when an error is raised, so the cost of scanning the fields does not matter.
+        """
+        if not is_dataclass(self):
+            return None
+        return next((field.name for field in fields(self) if field.metadata.get("array_attribute")), None)
+
+    def _refuse_array_conversion(self, attempted: str = "") -> Never:
+        """Raise a `TypeError` explaining that this object can not be used as an array."""
+        msg = f"`{type(self).__name__}` objects can not be used as an array{attempted}."
+        if self._array_attribute:
+            msg += f" Pass the `{self._array_attribute}` attribute instead."
+        raise TypeError(msg)
+
+    def __array__(self, dtype: np.dtype | None = None, copy: bool | None = None) -> Never:
+        self._refuse_array_conversion()
+
+    def __array_ufunc__(self, ufunc: np.ufunc, method: str, *inputs: Any, **kwargs: Any) -> Never:  # noqa: ANN401
+        self._refuse_array_conversion(f" (attempted `numpy.{ufunc.__name__}`)")
+
+    def __array_function__(self, func: Callable, types: Any, args: Any, kwargs: Any) -> Never:  # noqa: ANN401
+        self._refuse_array_conversion(f" (attempted `{func.__module__}.{func.__name__}`)")

eitprocessing/datahandling/mixins/slicing.py

@@ -10,16 +10,40 @@
 import numpy as np
 
 if TYPE_CHECKING:
-    from typing_extensions import Self
+    from typing_extensions import Never, Self
 
 
-class SelectByIndex(ABC):
+class NotIterable:
+    """Mixin class that prevents iteration through the legacy sequence protocol.
+
+    A class that implements `__getitem__` but not `__iter__` is iterable whether it wants to be or not: `iter()` falls
+    back to calling `__getitem__` with 0, 1, 2, ... until `IndexError` is raised. For a class whose `__getitem__`
+    slices rather than returning single elements, that is never what the caller meant. Each call returns a copy of the
+    whole object holding a single sample, which is quadratic, and an out-of-range index returns an empty copy rather
+    than raising `IndexError`, so `list(obj)` and `for item in obj` never terminate at all.
+
+    Defining `__iter__` opts out of that fallback, so iteration raises a `TypeError` explaining what to iterate over
+    instead.
+    """
+
+    def __iter__(self) -> Never:
+        """Refuse to iterate, rather than falling back to `__getitem__`."""
+        msg = f"`{type(self).__name__}` objects are not iterable."
+        if attribute := getattr(self, "_array_attribute", None):
+            msg += f" Iterate over the `{attribute}` attribute instead."
+        raise TypeError(msg)
+
+
+class SelectByIndex(NotIterable, ABC):
     """Adds slicing functionality to subclass by implementing `__getitem__`.
 
     Subclasses must implement a `_sliced_copy` function that defines what should
     happen when the object is sliced. This class ensures that when calling a
     slice between square brackets (as e.g. done for lists) then return the
     expected sliced object.
+
+    Subscripting means slicing here, not element access, so subclasses are not sequences. `NotIterable` makes that
+    explicit by blocking the iteration Python would otherwise infer from `__getitem__`.
     """
 
     label: str

eitprocessing/datahandling/pixelmap.py

@@ -23,8 +23,9 @@
 only a single plotting configuration).
 
 Mathematical Operations:
-    `PixelMap` instances support basic mathematical operations (+, -, *, /) with other `PixelMap` instances, arrays, or
-    scalar values. The operations are applied element-wise to the underlying values.
+    `PixelMap` instances support basic mathematical operations (+, -, *, /) with other `PixelMap` instances or scalar
+    values. The operations are applied element-wise to the underlying values. Arrays are not supported: convert them
+    with `PixelMap(values)` first, so that shape and type are validated explicitly.
 
     - Addition (+): Returns a `PixelMap` with values added element-wise.
 
@@ -36,8 +37,8 @@
       warning.
 
     When operating with another `PixelMap` of any type, operations typically return the base PixelMap type, except for
-    subtraction which returns a DifferenceMap. When operating with scalars or arrays, operations return the same type as
-    the original `PixelMap`.
+    subtraction which returns a DifferenceMap. When operating with scalars, operations return the same type as the
+    original `PixelMap`.
 
     Note: Some `PixelMap` subclasses (like `TIVMap` and `PerfusionMap`) do not allow negative values. Operations that
     might produce negative values with these maps will display appropriate warnings.
@@ -54,6 +55,7 @@
 from numpy import typing as npt
 from typing_extensions import Self
 
+from eitprocessing.datahandling.mixins.arrays import NotAnArray
 from eitprocessing.utils import make_capture
 
 if TYPE_CHECKING:
@@ -67,7 +69,7 @@
 
 
 @dataclass(frozen=True)
-class PixelMap:
+class PixelMap(NotAnArray):
     """Map representing a single value for each pixel.
 
     At initialization, values are conveted to a 2D numpy array of floats. The values are immutable after initialization,
@@ -88,7 +90,7 @@ class PixelMap:
         allow_negative_values (bool): Whether negative values are allowed in the pixel map.
     """
 
-    values: np.ndarray
+    values: np.ndarray = field(metadata={"array_attribute": True})
     _: KW_ONLY
     label: str | None = None
     plot_config: InitVar[PixelMapPlotConfig]
@@ -428,47 +430,86 @@ def __replace__(self, /, **changes) -> Self:
 
     update = __replace__
 
-    def _validate_other(self, other: npt.ArrayLike | float | PixelMap) -> np.ndarray | float:
-        other_values = other.values if isinstance(other, PixelMap) else other
-
-        if isinstance(other, (float, int)):
+    #: Ufuncs backing an arithmetic operator, mapped to the reflected operator handling `other <operator> pixel_map`.
+    _reflected_operators: ClassVar[dict[str, str]] = {
+        "add": "__radd__",
+        "subtract": "__rsub__",
+        "multiply": "__rmul__",
+        "divide": "__rtruediv__",
+    }
+
+    def __array_ufunc__(self, ufunc: np.ufunc, method: str, *inputs, **kwargs) -> PixelMap:  # noqa: RET503
+        """Handle `other <operator> pixel_map`, and refuse all other numpy operations.
+
+        `PixelMap` is not an array (see `NotAnArray`), so numpy operations on it are refused. The exception is an
+        arithmetic operator with the pixel map on the right hand side of a non-pixel map, e.g.
+        `np.float64(3) * pixel_map`. Numpy routes those here rather than deferring to Python's reflected operator
+        protocol, so they are forwarded to the reflected operator, which accepts scalars and rejects arrays with a
+        helpful message. Calling a ufunc on two pixel maps, e.g. `np.multiply(map_a, map_b)`, is refused: use the
+        operator instead.
+        """
+        if (
+            method == "__call__"
+            and not kwargs
+            and len(inputs) == 2  # noqa: PLR2004, ignore hardcoded value
+            and inputs[1] is self
+            and not isinstance(inputs[0], PixelMap)
+            and (reflected := self._reflected_operators.get(ufunc.__name__)) is not None
+        ):
+            return getattr(self, reflected)(inputs[0])
+
+        self._refuse_array_conversion(f" (attempted `numpy.{ufunc.__name__}`)")
+
+    def _validate_other(self, other: float | PixelMap) -> np.ndarray | float:
+        """Check that `other` can be combined with this pixel map, and return the values to operate on.
+
+        Pixel maps can only be combined with other pixel maps of the same shape, or with scalars. Arrays are
+        deliberately rejected: an array carries no label or plotting configuration, and silently treating it as a
+        pixel map hides shape and unit mistakes. Convert it explicitly with `PixelMap(array)` instead.
+        """
+        if isinstance(other, (int, float, np.number)):
             return other
 
-        other_values = np.array(other_values)
+        if not isinstance(other, PixelMap):
+            msg = (
+                f"Can't combine `{type(self).__name__}` with `{type(other).__name__}`. Pixel maps can only be "
+                "combined with other pixel maps or with scalars. Convert an array with `PixelMap(values)` first."
+            )
+            raise TypeError(msg)
 
-        if (os := other_values.shape) != (ss := self.values.shape):
+        if (os := other.values.shape) != (ss := self.values.shape):
             msg = f"Shape of PixelMaps (self: {ss}, other: {os}) do not match."
             raise ValueError(msg)
 
-        return other_values
+        return other.values
 
-    def __add__(self, other: npt.ArrayLike | float | PixelMap) -> PixelMap:
+    def __add__(self, other: float | PixelMap) -> PixelMap:
         new_values = self.values + self._validate_other(other)
         if isinstance(other, PixelMap):
             return PixelMap(new_values)
         return self.update(values=new_values, label=None)
 
     __radd__ = __add__
 
-    def __sub__(self, other: npt.ArrayLike | float | PixelMap) -> PixelMap:
+    def __sub__(self, other: float | PixelMap) -> PixelMap:
         new_values = self.values - self._validate_other(other)
         if isinstance(other, PixelMap):
             return DifferenceMap(new_values)
         return self.update(values=new_values, label=None)
 
-    def __rsub__(self, other: npt.ArrayLike | float | PixelMap) -> PixelMap:
+    def __rsub__(self, other: float | PixelMap) -> PixelMap:
         new_values = -self.values + self._validate_other(other)
         return self.update(values=new_values, label=None)
 
-    def __mul__(self, other: npt.ArrayLike | float | PixelMap) -> PixelMap:
+    def __mul__(self, other: float | PixelMap) -> PixelMap:
         new_values = self.values * self._validate_other(other)
         if isinstance(other, PixelMap):
             return PixelMap(new_values)
         return self.update(values=new_values, label=None)
 
     __rmul__ = __mul__
 
-    def __truediv__(self, other: npt.ArrayLike | float | PixelMap) -> PixelMap:
+    def __truediv__(self, other: float | PixelMap) -> PixelMap:
         other_values = self._validate_other(other)
         if isinstance(other_values, np.ndarray) and 0 in other_values:
             warnings.warn("Dividing by 0 will result in `np.nan` value.", UserWarning, stacklevel=2)
@@ -480,7 +521,7 @@ def __truediv__(self, other: npt.ArrayLike | float | PixelMap) -> PixelMap:
             return PixelMap(new_values)
         return self.update(values=new_values, label=None)
 
-    def __rtruediv__(self, other: npt.ArrayLike | float | PixelMap) -> PixelMap:
+    def __rtruediv__(self, other: float | PixelMap) -> PixelMap:
         other_values = self._validate_other(other)
         if 0 in self.values:
             warnings.warn("Dividing by 0 will result in `np.nan` value.", UserWarning, stacklevel=2)

eitprocessing/datahandling/sequence.py

@@ -9,6 +9,7 @@
 from eitprocessing.datahandling.datacollection import DataCollection
 from eitprocessing.datahandling.eitdata import EITData
 from eitprocessing.datahandling.intervaldata import IntervalData
+from eitprocessing.datahandling.mixins.arrays import NotAnArray
 from eitprocessing.datahandling.mixins.equality import Equivalence
 from eitprocessing.datahandling.mixins.slicing import SelectByTime
 from eitprocessing.datahandling.sparsedata import SparseData
@@ -25,7 +26,7 @@
 
 
 @dataclass(eq=False)
-class Sequence(Equivalence, SelectByTime):
+class Sequence(Equivalence, SelectByTime, NotAnArray):
     """Sequence of timepoints containing respiratory data.
 
     A Sequence object is a representation of data points over time. These data can consist of any combination of EIT

eitprocessing/datahandling/sparsedata.py

@@ -43,7 +43,7 @@ class SparseData(DataContainer, SelectByTime):
     description: str = field(compare=False, default="", repr=False)
     parameters: dict[str, Any] = field(default_factory=dict, metadata={"check_equivalence": True}, repr=False)
     derived_from: list[Any] = field(default_factory=list, compare=False, repr=False)
-    values: Any | None = None
+    values: Any | None = field(default=None, metadata={"array_attribute": True})
 
     def __repr__(self) -> str:
         return f"{self.__class__.__name__}('{self.label}')"

eitprocessing/roi/__init__.py

@@ -34,6 +34,7 @@
 from typing_extensions import Self
 
 from eitprocessing.datahandling.eitdata import EITData
+from eitprocessing.datahandling.mixins.arrays import NotAnArray
 from eitprocessing.datahandling.pixelmap import PixelMap
 
 if TYPE_CHECKING:
@@ -43,7 +44,7 @@
 
 
 @dataclass(frozen=True)
-class PixelMask:  # noqa: PLW1641
+class PixelMask(NotAnArray):  # noqa: PLW1641
     """Mask pixels by selecting or weighing them individually.
 
     A mask is a 2D array with a value for each pixel. Most often, this value is NaN (`np.nan`, 'not a number') or 1, and
@@ -92,7 +93,7 @@ class PixelMask:  # noqa: PLW1641
 
     """
 
-    mask: np.ndarray
+    mask: np.ndarray = field(metadata={"array_attribute": True})
     plot_config: InitVar[PixelMapPlotConfig]
     label: str | None = None
     keep_zeros: InitVar[bool] = field(default=False, kw_only=True)

eitprocessing/roi/pixelmaskcollection.py

@@ -9,14 +9,15 @@
 from typing_extensions import Self
 
 from eitprocessing.datahandling.eitdata import EITData
+from eitprocessing.datahandling.mixins.arrays import NotAnArray
 from eitprocessing.datahandling.pixelmap import PixelMap
 from eitprocessing.roi import PixelMask
 
 T = TypeVar("T", bound=EITData | PixelMap)
 
 
 @dataclass(frozen=True)
-class PixelMaskCollection:
+class PixelMaskCollection(NotAnArray):
     """A collection of pixel masks, each representing a specific region of interest (ROI) in the EIT data.
 
     This class allows for the application of multiple masks to numpy arrays, EIT data or pixel maps, enabling the

eitprocessing/roi/watershed.py

@@ -170,7 +170,7 @@ def apply(  # noqa: PLR0915
             included_peaks = np.argwhere(included_marker_indices)
             excluded_peaks = np.argwhere(peaks_loc_bool & ~included_marker_indices)
 
-            included_watershed_regions = np.where(included_region, watershed_regions, np.nan)
+            included_watershed_regions = np.where(included_region.mask, watershed_regions, np.nan)
 
             capture("included peaks", included_peaks)
             capture("excluded peaks", excluded_peaks)