diff --git a/eitprocessing/datahandling/__init__.py b/eitprocessing/datahandling/__init__.py
index 190182068..74025328d 100644
--- a/eitprocessing/datahandling/__init__.py
+++ b/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):
diff --git a/eitprocessing/datahandling/continuousdata.py b/eitprocessing/datahandling/continuousdata.py
index 8ca919465..d8d0f80d1 100644
--- a/eitprocessing/datahandling/continuousdata.py
+++ b/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:
diff --git a/eitprocessing/datahandling/datacollection.py b/eitprocessing/datahandling/datacollection.py
index 25bd8874a..db256b53d 100644
--- a/eitprocessing/datahandling/datacollection.py
+++ b/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
diff --git a/eitprocessing/datahandling/eitdata.py b/eitprocessing/datahandling/eitdata.py
index ea50d9bf0..77a5d509b 100644
--- a/eitprocessing/datahandling/eitdata.py
+++ b/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:
diff --git a/eitprocessing/datahandling/intervaldata.py b/eitprocessing/datahandling/intervaldata.py
index 2ae95826a..cdb9a1ddc 100644
--- a/eitprocessing/datahandling/intervaldata.py
+++ b/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)
diff --git a/eitprocessing/datahandling/mixins/arrays.py b/eitprocessing/datahandling/mixins/arrays.py
new file mode 100644
index 000000000..1a86b72c9
--- /dev/null
+++ b/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__}`)")
diff --git a/eitprocessing/datahandling/mixins/slicing.py b/eitprocessing/datahandling/mixins/slicing.py
index e48e620fe..45fe728b2 100644
--- a/eitprocessing/datahandling/mixins/slicing.py
+++ b/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
diff --git a/eitprocessing/datahandling/pixelmap.py b/eitprocessing/datahandling/pixelmap.py
index 87f441f02..90def0039 100644
--- a/eitprocessing/datahandling/pixelmap.py
+++ b/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,21 +430,60 @@ 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)
@@ -450,17 +491,17 @@ def __add__(self, other: npt.ArrayLike | float | PixelMap) -> PixelMap:
 
     __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)
@@ -468,7 +509,7 @@ def __mul__(self, other: npt.ArrayLike | float | PixelMap) -> PixelMap:
 
     __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)
diff --git a/eitprocessing/datahandling/sequence.py b/eitprocessing/datahandling/sequence.py
index f4a4b768f..623aa523b 100644
--- a/eitprocessing/datahandling/sequence.py
+++ b/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
diff --git a/eitprocessing/datahandling/sparsedata.py b/eitprocessing/datahandling/sparsedata.py
index 2f0bd8365..813f61c2f 100644
--- a/eitprocessing/datahandling/sparsedata.py
+++ b/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}')"
diff --git a/eitprocessing/roi/__init__.py b/eitprocessing/roi/__init__.py
index 1d01472f1..592dafda6 100644
--- a/eitprocessing/roi/__init__.py
+++ b/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)
diff --git a/eitprocessing/roi/pixelmaskcollection.py b/eitprocessing/roi/pixelmaskcollection.py
index 28168e2d4..69c9c01cc 100644
--- a/eitprocessing/roi/pixelmaskcollection.py
+++ b/eitprocessing/roi/pixelmaskcollection.py
@@ -9,6 +9,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
 from eitprocessing.roi import PixelMask
 
@@ -16,7 +17,7 @@
 
 
 @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
diff --git a/eitprocessing/roi/watershed.py b/eitprocessing/roi/watershed.py
index 3820bbed5..62636b56e 100644
--- a/eitprocessing/roi/watershed.py
+++ b/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)
diff --git a/tests/mixins/conftest.py b/tests/mixins/conftest.py
new file mode 100644
index 000000000..ce2adac7c
--- /dev/null
+++ b/tests/mixins/conftest.py
@@ -0,0 +1,91 @@
+import numpy as np
+import pytest
+
+from eitprocessing.datahandling.continuousdata import ContinuousData
+from eitprocessing.datahandling.datacollection import DataCollection
+from eitprocessing.datahandling.eitdata import EITData, Vendor
+from eitprocessing.datahandling.intervaldata import IntervalData
+from eitprocessing.datahandling.pixelmap import PixelMap
+from eitprocessing.datahandling.sequence import Sequence
+from eitprocessing.datahandling.sparsedata import SparseData
+from eitprocessing.roi import PixelMask
+from eitprocessing.roi.pixelmaskcollection import PixelMaskCollection
+
+
+@pytest.fixture
+def data_object(request: pytest.FixtureRequest) -> object:
+    """Return the object named by the parameter, for `indirect` parametrization over several types of object."""
+    return request.getfixturevalue(request.param)
+
+
+@pytest.fixture
+def continuous_data() -> ContinuousData:
+    """Return a ContinuousData fixture with a short ramp as values."""
+    return ContinuousData(
+        label="cd",
+        name="cd",
+        unit="a.u.",
+        category="impedance",
+        time=np.arange(10.0),
+        values=np.arange(10.0),
+        sample_frequency=1.0,
+    )
+
+
+@pytest.fixture
+def eit_data() -> EITData:
+    """Return an EITData fixture with four frames of 2x2 pixels."""
+    return EITData(
+        path="somewhere",
+        nframes=4,
+        time=np.arange(4.0),
+        sample_frequency=20.0,
+        vendor=Vendor.DRAEGER,
+        pixel_impedance=np.ones((4, 2, 2)),
+        suppress_simulated_warning=True,
+    )
+
+
+@pytest.fixture
+def sparse_data() -> SparseData:
+    """Return a SparseData fixture with three values."""
+    return SparseData(label="sd", name="sd", unit=None, category="breath", time=np.arange(3.0), values=[1, 2, 3])
+
+
+@pytest.fixture
+def interval_data() -> IntervalData:
+    """Return an IntervalData fixture with a single interval."""
+    return IntervalData(label="id", name="id", unit=None, category="breath", intervals=[(0.0, 1.0)], values=[1])
+
+
+@pytest.fixture
+def empty_sequence() -> Sequence:
+    """Return an empty Sequence.
+
+    Named `empty_sequence` to avoid shadowing the `sequence` fixture in `tests/conftest.py`, which loads test data.
+    """
+    return Sequence()
+
+
+@pytest.fixture
+def data_collection() -> DataCollection:
+    """Return an empty DataCollection."""
+    return DataCollection(EITData)
+
+
+@pytest.fixture
+def pixel_map() -> PixelMap:
+    """Return a PixelMap fixture with 2x2 pixels."""
+    return PixelMap([[1.0, 2.0], [3.0, 4.0]])
+
+
+@pytest.fixture
+def pixel_mask() -> PixelMask:
+    """Return a PixelMask fixture with 2x2 pixels."""
+    return PixelMask(np.ones((2, 2)))
+
+
+@pytest.fixture
+def pixel_mask_collection() -> PixelMaskCollection:
+    """Return an empty PixelMaskCollection."""
+    return PixelMaskCollection()
diff --git a/tests/mixins/test_arrays.py b/tests/mixins/test_arrays.py
new file mode 100644
index 000000000..cdbebe73c
--- /dev/null
+++ b/tests/mixins/test_arrays.py
@@ -0,0 +1,87 @@
+import numpy as np
+import pytest
+import scipy.signal
+
+from eitprocessing.datahandling.continuousdata import ContinuousData
+from eitprocessing.datahandling.mixins.arrays import NotAnArray
+
+
+@pytest.mark.parametrize(
+    ("data_object", "attribute"),
+    [
+        ("pixel_map", "values"),
+        ("pixel_mask", "mask"),
+        ("continuous_data", "values"),
+        ("eit_data", "pixel_impedance"),
+        ("sparse_data", "values"),
+        ("interval_data", "values"),
+        ("empty_sequence", None),
+        ("data_collection", None),
+        ("pixel_mask_collection", None),
+    ],
+    indirect=["data_object"],
+)
+def test_error_names_the_data_attribute(data_object: NotAnArray, attribute: str | None):
+    """Every object using the mixin refuses conversion and points at its own data attribute, if it has one."""
+    with pytest.raises(
+        TypeError, match=f"`{type(data_object).__name__}` objects can not be used as an array"
+    ) as excinfo:
+        np.asarray(data_object)
+
+    if attribute:
+        assert f"Pass the `{attribute}` attribute instead." in str(excinfo.value)
+    else:
+        assert "Pass the" not in str(excinfo.value)
+
+
+def test_conversion_to_array_is_refused(continuous_data: ContinuousData):
+    """`__array__`: `numpy.asarray()` and friends."""
+    for convert in (np.asarray, np.array, lambda obj: np.array(obj, dtype=float)):
+        with pytest.raises(TypeError, match="can not be used as an array"):
+            convert(continuous_data)
+
+
+def test_ufuncs_are_refused(continuous_data: ContinuousData):
+    """`__array_ufunc__`: `numpy.sin()`, `numpy.add()`, `array + object`, ..."""
+    with pytest.raises(TypeError, match=r"can not be used as an array \(attempted `numpy.sin`\)"):
+        np.sin(continuous_data)
+
+    with pytest.raises(TypeError, match=r"can not be used as an array \(attempted `numpy.add`\)"):
+        np.add(np.ones(10), continuous_data)
+
+    with pytest.raises(TypeError, match=r"can not be used as an array \(attempted `numpy.add`\)"):
+        _ = np.ones(10) + continuous_data
+
+
+def test_array_functions_are_refused(continuous_data: ContinuousData):
+    """`__array_function__`: the rest of the numpy API, which dispatches before conversion."""
+    with pytest.raises(TypeError, match=r"can not be used as an array \(attempted `numpy.mean`\)"):
+        np.mean(continuous_data)
+
+    with pytest.raises(TypeError, match=r"can not be used as an array \(attempted `numpy.concatenate`\)"):
+        np.concatenate([continuous_data, continuous_data])
+
+    with pytest.raises(TypeError, match=r"can not be used as an array \(attempted `numpy.stack`\)"):
+        np.stack([continuous_data])
+
+
+def test_scipy_functions_are_refused(continuous_data: ContinuousData):
+    """Scipy converts its input to an array before doing anything else, so `__array__` covers it."""
+    with pytest.raises(TypeError, match="can not be used as an array"):
+        scipy.signal.detrend(continuous_data)
+
+    with pytest.raises(TypeError, match="can not be used as an array"):
+        scipy.signal.filtfilt([1.0], [1.0], continuous_data)
+
+
+def test_slicing_is_unaffected(continuous_data: ContinuousData):
+    """The mixin should only block numpy; regular Python behaviour must keep working."""
+    assert len(continuous_data) == 10
+    assert len(continuous_data[2:5]) == 3
+    assert np.array_equal(continuous_data[2:5].values, np.arange(2.0, 5.0))
+    assert continuous_data.t[2.0:5.0] == continuous_data[2:5]
+
+
+def test_underlying_data_can_still_be_used(continuous_data: ContinuousData):
+    """The suggestion in the error message should actually work."""
+    assert np.mean(continuous_data.values) == pytest.approx(4.5)
diff --git a/tests/mixins/test_slicing.py b/tests/mixins/test_slicing.py
index e32e0a3dd..0598cfc2d 100644
--- a/tests/mixins/test_slicing.py
+++ b/tests/mixins/test_slicing.py
@@ -1,9 +1,79 @@
 import pytest
 
+from eitprocessing.datahandling.continuousdata import ContinuousData
 from eitprocessing.datahandling.loading import load_eit_data
+from eitprocessing.datahandling.mixins.slicing import NotIterable, SelectByIndex
 from eitprocessing.datahandling.sequence import Sequence
 
 
+def test_select_by_index_is_not_iterable():
+    """The opt-out is inherited by every sliceable class, so a new one can not forget it."""
+    assert issubclass(SelectByIndex, NotIterable)
+
+
+@pytest.mark.parametrize(
+    ("data_object", "attribute"),
+    [
+        ("continuous_data", "values"),
+        ("eit_data", "pixel_impedance"),
+        ("sparse_data", "values"),
+        ("interval_data", "values"),
+        ("empty_sequence", None),
+    ],
+    indirect=["data_object"],
+)
+def test_iteration_is_refused(data_object: NotIterable, attribute: str | None):
+    """Without `__iter__`, iteration falls back to `__getitem__` and never terminates.
+
+    Each of these raises rather than hanging, which is what the test is really checking.
+    """
+    with pytest.raises(TypeError, match=f"`{type(data_object).__name__}` objects are not iterable") as excinfo:
+        list(data_object)
+
+    if attribute:
+        assert f"Iterate over the `{attribute}` attribute instead." in str(excinfo.value)
+    else:
+        assert "Iterate over" not in str(excinfo.value)
+
+
+def test_iteration_is_refused_in_every_form(continuous_data: ContinuousData):
+    with pytest.raises(TypeError, match="`ContinuousData` objects are not iterable"):
+        iter(continuous_data)
+
+    with pytest.raises(TypeError, match="`ContinuousData` objects are not iterable"):
+        for _item in continuous_data:
+            pass
+
+    with pytest.raises(TypeError, match="`ContinuousData` objects are not iterable"):
+        _first, _second = continuous_data
+
+    with pytest.raises(TypeError, match="`ContinuousData` objects are not iterable"):
+        sum(continuous_data)
+
+
+def test_slicing_still_works_without_iteration(continuous_data: ContinuousData):
+    assert len(continuous_data) == 10
+    assert len(continuous_data[2:5]) == 3
+    assert continuous_data[2:5] == continuous_data.select_by_index(2, 5)
+    assert list(continuous_data.values) == list(range(10))
+
+
+def test_data_collection_still_iterates_over_keys(draeger_20hz_healthy_volunteer: Sequence):
+    """`DataCollection` is a `UserDict` and does not slice by index, so it stays iterable."""
+    collection = draeger_20hz_healthy_volunteer.eit_data
+
+    assert not isinstance(collection, NotIterable)
+    assert list(collection) == list(collection.keys())
+
+
+def test_sequence_data_still_iterates_over_labels(draeger_20hz_healthy_volunteer: Sequence):
+    """`Sequence` is not iterable, but its `data` accessor is."""
+    assert list(draeger_20hz_healthy_volunteer.data) == list(draeger_20hz_healthy_volunteer.data.keys())
+
+    with pytest.raises(TypeError, match="`Sequence` objects are not iterable"):
+        list(draeger_20hz_healthy_volunteer)
+
+
 @pytest.mark.parametrize(
     "sequence",
     ["draeger_20hz_healthy_volunteer", "draeger_20hz_healthy_volunteer_pressure_pod"],
diff --git a/tests/test_parameter_tiv.py b/tests/test_parameter_tiv.py
index 3410e4eda..c5b0fd8fb 100644
--- a/tests/test_parameter_tiv.py
+++ b/tests/test_parameter_tiv.py
@@ -74,7 +74,7 @@ def mock_continuous_data():
         category="relative impedance",
         description="Global impedance created for testing TIV parameter",
         parameters={},
-        derived_from="mock_eit_data",
+        derived_from=["mock_eit_data"],
         time=np.linspace(0, 18, (18 * 1000), endpoint=False),
         values=mock_global_impedance(),
         sample_frequency=1000,
diff --git a/tests/test_pixel_breath.py b/tests/test_pixel_breath.py
index 43ee5791b..d78f618b8 100644
--- a/tests/test_pixel_breath.py
+++ b/tests/test_pixel_breath.py
@@ -63,7 +63,7 @@ def mock_continuous_data():
         category="relative impedance",
         description="Global impedance created for testing pixel breath feature",
         parameters={},
-        derived_from="mock_eit_data",
+        derived_from=["mock_eit_data"],
         time=np.linspace(0, 2 * np.pi, 400),
         values=mock_global_impedance(),
         sample_frequency=399 / 2 * np.pi,
diff --git a/tests/test_pixelmap.py b/tests/test_pixelmap.py
index ebc3f9c15..162a8c843 100644
--- a/tests/test_pixelmap.py
+++ b/tests/test_pixelmap.py
@@ -1,6 +1,8 @@
 import copy
+import operator as op
 import sys
 import warnings
+from collections.abc import Callable
 from dataclasses import FrozenInstanceError
 
 import frozendict
@@ -9,6 +11,7 @@
 from matplotlib import pyplot as plt
 from matplotlib.colors import CenteredNorm, Colormap, Normalize
 from matplotlib.ticker import PercentFormatter, ScalarFormatter
+from numpy import typing as npt
 from numpy.exceptions import ComplexWarning
 
 from eitprocessing.datahandling.pixelmap import (
@@ -24,6 +27,7 @@
 from eitprocessing.plotting import _PLOT_CONFIG_REGISTRY, reset_plot_config, set_plot_config_parameters
 from eitprocessing.plotting.helpers import AbsolutePercentFormatter, AbsoluteScalarFormatter
 from eitprocessing.plotting.pixelmap import PixelMapPlotConfig
+from eitprocessing.roi import PixelMask
 
 
 def test_init_values():
@@ -506,6 +510,80 @@ def test_div():
     assert np.array_equal(pm2_div_pm1.values, [[np.nan, 2, 3 / 2, 4 / 3]], equal_nan=True)
 
 
+MATCH_NOT_COMBINABLE = "Pixel maps can only be combined with other pixel maps or with scalars"
+
+
+@pytest.mark.parametrize(
+    "other",
+    [
+        pytest.param(np.array([[1.0, 2.0, 3.0, 4.0]]), id="array"),
+        pytest.param(np.array([[1.0, 2.0]]), id="array of different shape"),
+        pytest.param([[1, 2, 3, 4]], id="nested list"),
+        pytest.param(np.arange(4.0), id="1D array"),
+    ],
+)
+@pytest.mark.parametrize("operator", [op.add, op.sub, op.mul, op.truediv])
+def test_operators_reject_arrays(operator: Callable, other: npt.ArrayLike):
+    """Pixel maps combine with pixel maps and scalars, but never with raw arrays."""
+    pm = PerfusionMap([[0.0, 1.0, 2.0, 3.0]])
+
+    with pytest.raises(TypeError, match=MATCH_NOT_COMBINABLE):
+        _ = operator(pm, other)
+
+    with pytest.raises(TypeError, match=MATCH_NOT_COMBINABLE):
+        _ = operator(other, pm)
+
+
+@pytest.mark.parametrize("operator", [op.add, op.sub, op.mul, op.truediv])
+def test_operators_reject_other_types(operator: Callable):
+    pm = PerfusionMap([[1.0, 2.0]])
+
+    with pytest.raises(TypeError, match=MATCH_NOT_COMBINABLE):
+        _ = operator(pm, "not a pixel map")
+
+    with pytest.raises(TypeError, match=MATCH_NOT_COMBINABLE):
+        _ = operator(pm, PixelMask(np.ones((1, 2))))
+
+
+@pytest.mark.parametrize("scalar_type", [int, float, np.int64, np.float64])
+@pytest.mark.parametrize(
+    ("operator", "expected"),
+    [
+        (op.add, [[3.0, 5.0]]),
+        (op.mul, [[2.0, 6.0]]),
+        (op.sub, [[-1.0, 1.0]]),
+        (op.truediv, [[0.5, 1.5]]),
+    ],
+)
+def test_operators_accept_scalars_on_either_side(operator: Callable, expected: list, scalar_type: type):
+    """Numpy scalars do not defer to Python's reflected operators, so `__array_ufunc__` handles them."""
+    pm = PerfusionMap([[1.0, 3.0]])
+    scalar = scalar_type(2)
+
+    result = operator(pm, scalar)
+    assert isinstance(result, PerfusionMap)
+    assert np.array_equal(result.values, expected)
+
+    # The reflected operators compute `scalar <operator> pixel_map`, which differs for non-commutative operators.
+    reflected = operator(scalar, pm)
+    assert isinstance(reflected, PerfusionMap)
+    assert np.array_equal(reflected.values, operator(np.array([[2.0, 2.0]]), pm.values))
+
+
+@pytest.mark.parametrize(
+    "function",
+    [np.mean, np.sin, np.asarray, np.nanmax, lambda pm: np.stack([pm]), lambda pm: np.multiply(pm, pm)],
+)
+def test_numpy_functions_are_refused(function: Callable):
+    """A pixel map is not an array; numpy should refuse it rather than build an object array."""
+    with pytest.raises(TypeError, match="`PixelMap` objects can not be used as an array"):
+        function(PixelMap([[1.0, 2.0]]))
+
+
+def test_numpy_functions_work_on_values():
+    assert np.mean(PixelMap([[1.0, 3.0]]).values) == pytest.approx(2.0)
+
+
 def test_nan():
     # preventing 0-values, because the will lead to extra nan-values when dividing
     pm1 = PendelluftMap(np.reshape(np.arange(1, 101), (10, 10)))