Merge branch 'main' into testing/cuda_cf_api

Three files conflicted:

- ClusterFinderMT::collect(): both sides rewrote it. Keep this branch's
  version (f4bbbc9), which reads the stop flags before sweeping, over
  main's (afce91b), which still samples queue emptiness mid-condition.
  The staged splice had also lost a brace, swallowing popFront().
- python/aare/ClusterFinder.py: restore the CUDA API deleted by the
  resolution (_cuda_available, _CUDA_ALGORITHMS, ClusterFinderCUDA,
  find_cluster_views_batched_iter, ClusterFinderCUDAGraph), and main's
  min_pedestal_samples forwarding, dropped as an adjacent hunk.
- bind_ClusterVector.hpp: keep py::module_local() (133cedf) alongside
  main's new docstrings; without it _aare and _aare_cuda cannot both
  register Cluster/ClusterVector under pybind11 3.

Does not build on its own: ClusterFinderFrozen and the CUDA drivers
still call Pedestal APIs main removed. Fixed in the follow-up.

python/CMakeLists.txt

@@ -69,6 +69,7 @@ set(PYTHON_FILES
     aare/ClusterVector.py
     aare/Cluster.py
     aare/FastPedestal.py
+    aare/Pedestal.py
     aare/calibration.py
     aare/factory.py
     aare/experimental.py

python/aare/ClusterFinder.py

@@ -18,17 +18,33 @@ def _get_class(name, cluster_size, dtype):
 
 
 
-def ClusterFinder(image_size, cluster_size=(3,3), n_sigma=5, dtype = np.int32, capacity = 1024):
+def ClusterFinder(image_size, *, cluster_size=(3,3), n_sigma=5, dtype = np.int32, capacity = 1024, min_pedestal_samples = 1000):
     """
     Factory function to create a ClusterFinder object. Provides a cleaner syntax for 
     the templated ClusterFinder in C++.
+
+    Parameters
+    ----------
+    image_size : tuple
+        The size of the image as a tuple (height, width).
+    cluster_size : tuple, optional
+        The size of the cluster to find as a tuple (height, width). Default is (3,3).
+    n_sigma : int, optional
+        Multiplier of the standard deviation used as a threshold to identify potential photon pixels. Default is 5.
+    dtype : data-type, optional
+        The data type of the image. Default is np.int32.
+    capacity : int, optional
+        The maximum number of clusters than can be stored before reallocating. Default is 1024.
+    min_pedestal_samples : int, optional
+        The minimum number of pedestal samples to accumulate before using the pedestal. Default is 1000.
     """
     cls = _get_class("ClusterFinder", cluster_size, dtype)
-    return cls(image_size, n_sigma=n_sigma, capacity=capacity)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity, min_pedestal_samples=min_pedestal_samples)
 
 
 
-def ClusterFinderFrozen(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32, capacity=1024):
+def ClusterFinderFrozen(image_size, *, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
+                        capacity=1024, min_pedestal_samples=1000):
     """
     Factory function to create a ClusterFinderFrozen object.
 
@@ -37,19 +53,40 @@ def ClusterFinderFrozen(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int3
     pedestal updates are deferred to the end of the frame). This mirrors the CUDA
     kernel's per-frame update model, so running it against ClusterFinderCUDA
     isolates pedestal-update timing as the sole variable.
+
+    min_pedestal_samples matches ClusterFinder so the twins warm up identically;
+    a comparison run must pass the same value to both.
     """
     cls = _get_class("ClusterFinderFrozen", cluster_size, dtype)
-    return cls(image_size, n_sigma=n_sigma, capacity=capacity)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity,
+               min_pedestal_samples=min_pedestal_samples)
 
 
-def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3):
+def ClusterFinderMT(image_size, *, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3, min_pedestal_samples = 1000): 
     """ 
     Factory function to create a ClusterFinderMT object. Provides a cleaner syntax for 
     the templated ClusterFinderMT in C++.
+
+    Parameters
+    ----------
+    image_size : tuple
+        The size of the image as a tuple (height, width).
+    cluster_size : tuple, optional
+        The size of the cluster to find as a tuple (height, width). Default is (3,3).
+    n_sigma : int, optional
+        Multiplier of the standard deviation used as a threshold to identify potential photon pixels. Default is 5.
+    dtype : data-type, optional
+        The data type of the image. Default is np.int32.
+    capacity : int, optional
+        The maximum number of clusters than can be stored before reallocating. Default is 1024.
+    n_threads : int, optional
+        The number of threads to use for processing. Default is 3.
+    min_pedestal_samples : int, optional
+        The minimum number of pedestal samples to accumulate before using the pedestal. Default is 1000.
     """
 
     cls = _get_class("ClusterFinderMT", cluster_size, dtype)
-    return cls(image_size, n_sigma=n_sigma, capacity=capacity, n_threads=n_threads)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity, min_pedestal_samples=min_pedestal_samples, n_threads=n_threads)
 
 
 def _cuda_available():
@@ -63,9 +100,10 @@ def _cuda_available():
 _CUDA_ALGORITHMS = {"fixed_window": ""}
 
 
-def ClusterFinderCUDA(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
+def ClusterFinderCUDA(image_size, *, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
                       max_clusters_per_frame=2048, n_streams=4,
-                      time_kernels=False, algorithm="fixed_window"):
+                      min_pedestal_samples=1000, algorithm="fixed_window",
+                      time_kernels=False):
     """
     Factory function to create a ClusterFinderCUDA object. Provides a cleaner
     syntax for the templated ClusterFinderCUDA in C++. API mirrors
@@ -87,6 +125,14 @@ def ClusterFinderCUDA(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
         but as tight as possible to minimize PCIe traffic. Default 2048.
     n_streams : int, optional
         Number of CUDA streams for H2D/kernel/D2H pipelining. Default 4.
+    min_pedestal_samples : int, optional
+        Frames the host pedestal accumulates before it is ready, and the
+        reciprocal of the steady-state update weight. Matches ClusterFinder so
+        a CPU/GPU comparison warms up identically. Default 1000.
+    algorithm : str, optional
+        GPU algorithm. "fixed_window" (default): a cluster_size window around
+        each pixel with local-maximum suppression. Memory and streams are
+        managed internally whichever algorithm is chosen.
     time_kernels : bool, optional
         Enable per-frame CUDA-event kernel timing, exposed via
         avg_kernel_time_ms(). Off by default because it adds two event records
@@ -95,10 +141,6 @@ def ClusterFinderCUDA(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
         contention the events measure queue wait, not execution, and over-read
         by up to ~3.5x. Use Nsight Systems for exclusive kernel times. When
         disabled, avg_kernel_time_ms() returns NaN.
-    algorithm : str, optional
-        GPU algorithm. "fixed_window" (default): a cluster_size window around
-        each pixel with local-maximum suppression. Memory and streams are
-        managed internally whichever algorithm is chosen.
 
     Example
     -------
@@ -137,6 +179,7 @@ def ClusterFinderCUDA(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
                n_sigma=n_sigma,
                max_clusters_per_frame=max_clusters_per_frame,
                n_streams=n_streams,
+               min_pedestal_samples=min_pedestal_samples,
                time_kernels=time_kernels)
 
 def find_cluster_views_batched_iter(cf, frames, first_frame=0, chunk=None):
@@ -202,8 +245,9 @@ def find_cluster_views_batched_iter(cf, frames, first_frame=0, chunk=None):
         view.release()
 
 
-def ClusterFinderCUDAGraph(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
-                           max_clusters_per_frame=2048, n_streams=4):
+def ClusterFinderCUDAGraph(image_size, *, cluster_size=(3,3), n_sigma=5, dtype=np.int32,
+                           max_clusters_per_frame=2048, n_streams=4,
+                           min_pedestal_samples=1000):
     """
     Factory function to create a ClusterFinderCUDAGraph object. Uses pre-recorded
     CUDA Graphs to reduce per-frame CPU API overhead (~23 µs vs ~31 µs for the
@@ -224,6 +268,8 @@ def ClusterFinderCUDAGraph(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.i
         Hard upper bound on clusters per frame. Default 2048.
     n_streams : int, optional
         Number of CUDA streams (one graph per stream). Default 4.
+    min_pedestal_samples : int, optional
+        Frames the host pedestal accumulates before it is ready. Default 1000.
 
     Note
     ----
@@ -240,7 +286,8 @@ def ClusterFinderCUDAGraph(image_size, cluster_size=(3,3), n_sigma=5, dtype=np.i
     return cls(image_size,
                n_sigma=n_sigma,
                max_clusters_per_frame=max_clusters_per_frame,
-               n_streams=n_streams)
+               n_streams=n_streams,
+               min_pedestal_samples=min_pedestal_samples)
 
 
 def ClusterCollector(clusterfindermt, dtype=np.int32):
@@ -263,19 +310,55 @@ def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32):
 
 
 def ClusterFile(fname, cluster_size=(3,3), dtype=np.int32, chunk_size = 1000, mode = "r"):
-    """
-    Factory function to create a ClusterFile object. Provides a cleaner syntax for
-    the templated ClusterFile in C++.
+    """Create a reader or writer for a binary cluster file.
+
+    Parameters
+    ----------
+    fname : path-like
+        Cluster file to open.
+    cluster_size : tuple[int, int], default=(3, 3)
+        Cluster dimensions stored in the file.
+    dtype : numpy dtype, default=numpy.int32
+        Data type of the cluster values stored in the file.
+    chunk_size : int, default=1000
+        Maximum number of selected clusters returned by ``chunks()`` and
+        default iteration. Must be positive when iterating over chunks.
+    mode : {"r", "w", "a"}, default="r"
+        Open for reading, truncate and write, or append, respectively.
+
+    Returns
+    -------
+    ClusterFile
+        The compiled ClusterFile specialization matching ``cluster_size`` and
+        ``dtype``.
+
+    Notes
+    -----
+    The file format contains no cluster shape or data-type metadata. Supplying
+    values that do not match the file causes its bytes to be interpreted
+    incorrectly. Use ``frames()`` to iterate over complete frames, including
+    empty or fully filtered frames with their stored frame numbers. Use
+    ``chunks()`` or ``chunks(chunk_size)`` to iterate over selected clusters
+    in batches. Chunks may combine frames, so their frame number is not
+    reliable per-cluster metadata.
+
+    Iterators consume the current file position without rewinding; use one
+    traversal at a time. Each result owns its storage and remains valid after
+    advancing the iterator or closing the file.
+
+    Examples
+    --------
 
     .. code-block:: python
 
         from aare import ClusterFile
-        
-        with ClusterFile("clusters.clust", cluster_size=(3,3), dtype=np.int32) as cf:
-            # cf is now a ClusterFile_Cluster3x3i object but you don't need to know that.
-            for clusters in cf:
-                # Loop over clusters in chunks of 1000 
-                # The type of clusters will be a ClusterVector_Cluster3x3i in this case
+
+        with ClusterFile(
+            "clusters.clust", cluster_size=(3, 3), dtype=np.int32
+        ) as cf:
+            for clusters in cf.chunks():
+                # Process clusters in chunks of at most 1000.
+                ...
 
     """
 

python/aare/ClusterVector.py

@@ -7,16 +7,33 @@ from .ClusterFinder import _get_class
 
 def ClusterVector(cluster_size=(3,3), dtype = np.int32):
     """
-    Factory function to create a ClusterVector object. Provides a cleaner syntax for
-    the templated ClusterVector in C++.
+    Create an empty ClusterVector for a supported cluster size and pixel dtype.
 
-    .. code-block:: python
+    Parameters
+    ----------
+    cluster_size : tuple[int, int], default=(3, 3)
+        Cluster dimensions in the x and y directions.
+    dtype : numpy.dtype, default=numpy.int32
+        Pixel storage type. The cluster size and dtype combination must have a
+        compiled binding.
 
-        from aare import ClusterVector
-        
-        ClusterVector(cluster_size=(3,3), dtype=np.float64)
+    Returns
+    -------
+    ClusterVector
+        An empty vector with frame number 0 and space reserved for at least
+        1024 clusters.
+
+    Raises
+    ------
+    ValueError
+        If the requested size and dtype combination is unavailable.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aare import ClusterVector
+    >>> clusters = ClusterVector(cluster_size=(3, 3), dtype=np.float64)
     """
 
     cls = _get_class("ClusterVector", cluster_size, dtype)
     return cls()
-

python/aare/FastPedestal.py

@@ -22,7 +22,7 @@ def FastPedestal(rows, cols, n_samples=1000, dtype=np.float64):
         cols: Number of image columns.
         n_samples: Initialization frame count and steady-state update-weight
             denominator.
-        dtype: Output dtype for the mean, variance, and standard deviation.
+        dtype: Output dtype for the mean and standard deviation.
             Supported values are ``np.float64``, ``np.float32``, and
             ``np.int16``.
     """
@@ -41,7 +41,7 @@ def from_file(filename, n_samples=1000, skip_first=0, dtype=np.float64):
         filename: Input image file.
         n_samples: Number of frames used for initialization.
         skip_first: Number of leading frames to ignore.
-        dtype: Output dtype for the mean, variance, and standard deviation.
+        dtype: Output dtype for the mean and standard deviation.
 
     Raises:
         RuntimeError: If fewer than ``n_samples`` frames remain after

python/aare/Pedestal.py

@@ -0,0 +1,26 @@
+# SPDX-License-Identifier: MPL-2.0
+
+import numpy as np
+
+from .factory import _get_typed_class
+
+
+def Pedestal(rows, cols, n_samples=1000, dtype=np.float64):
+    """Create an empty per-pixel running pedestal.
+
+    This factory hides the dtype suffix used by the templated C++ bindings.
+    Call ``push()`` to update the statistics and cached mean for each frame.
+    Statistics are available during initialization and are zero for empty pixels.
+    Internal moments and variance always use double precision.
+
+    Args:
+        rows: Number of image rows.
+        cols: Number of image columns.
+        n_samples: Number of samples accumulated before switching to
+            steady-state updates with weight ``1 / n_samples``.
+        dtype: Output dtype for the mean and standard deviation.
+            Supported values are ``np.float64``, ``np.float32``, and
+            ``np.int16``.
+    """
+    cls = _get_typed_class("Pedestal", dtype)
+    return cls(rows, cols, n_samples)

python/aare/RawFile.py

@@ -43,7 +43,7 @@ class RawFile(_aare.RawFile):
         Returns:
             int: Number of frames in file.
         """
-        return super().frames_in_file
+        return self.total_frames
 
     def __enter__(self):
         return self

python/aare/__init__.py

@@ -39,15 +39,18 @@ from ._aare import (
     RawSubFile,
 )
 from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder
-from ._aare import DetectorType, ReadoutMode 
+from ._aare import DetectorType, FrameDiscardPolicy, ReadoutMode
 from ._aare import hitmap
 from ._aare import ROI
 from ._aare import corner 
 
+from ._aare import UDPPortPosition
+
 # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
 
 from ._version import __version__
 from .FastPedestal import FastPedestal
+from .Pedestal import Pedestal
 from .ClusterFinder import (ClusterFinder, ClusterFinderFrozen, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile,
                             ClusterFinderCUDA, ClusterFinderCUDAGraph, _cuda_available, find_cluster_views_batched_iter)
 from .ClusterVector import ClusterVector

python/aare/calibration.py

@@ -5,9 +5,12 @@ def load_calibration(fname, hg0=False):
     """
     Load calibration data from a file.
     
-    Parameters:
-    fname (str): Path to the calibration file.
-    hg0 (bool): If True, load HG0 calibration data instead of G0.
+    Parameters
+    ----------
+    fname : str
+        Path to the calibration file.
+    hg0 : bool
+        If True, load HG0 calibration data instead of G0.
 
     """
     gains = 3

python/src/bind_Cluster.hpp

@@ -37,7 +37,6 @@ void define_Cluster(py::module &m, const std::string &typestr) {
             return cluster;
         }))
 
-        // TODO! Review if to keep or not
         .def_property_readonly(
             "data",
             [](Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> &c)
@@ -60,6 +59,26 @@ void define_Cluster(py::module &m, const std::string &typestr) {
         .def_readonly("y",
                       &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::y)
 
+        .def("__repr__",
+             [class_name](const Cluster<Type, ClusterSizeX, ClusterSizeY,
+                                        CoordType> &self) {
+                 fmt::memory_buffer data;
+                 fmt::format_to(std::back_inserter(data), "[[");
+                 for (size_t i = 0; i < self.data.size(); ++i) {
+                     if (i != 0) {
+                         fmt::format_to(std::back_inserter(data),
+                                        i % ClusterSizeX == 0 ? "], [" : ", ");
+                     }
+                     fmt::format_to(std::back_inserter(data), "{}",
+                                    self.data[i]);
+                 }
+                 fmt::format_to(std::back_inserter(data), "]]");
+
+                 return fmt::format("{}(x={}, y={}, data={})", class_name,
+                                    self.x, self.y,
+                                    fmt::string_view(data.data(), data.size()));
+             })
+
         .def(
             "max_sum_2x2",
             [](Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> &self) {
@@ -67,8 +86,16 @@ void define_Cluster(py::module &m, const std::string &typestr) {
                 return py::make_tuple(max_sum.sum,
                                       static_cast<int>(max_sum.index));
             },
-            R"(calculates sum of 2x2 subcluster with highest energy and index relative to cluster center 0: top_left, 1: top_right, 2: bottom_left, 3: bottom_right
-            )");
+            R"doc(
+            Return the highest-sum center-adjacent 2x2 subcluster.
+
+            Returns
+            -------
+            tuple
+                ``(sum, index)``, where index is 0 for top-left, 1 for
+                top-right, 2 for bottom-left, or 3 for bottom-right relative to
+                the cluster center.
+            )doc");
 }
 
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
@@ -80,7 +107,14 @@ void reduce_to_3x3(py::module &m) {
         [](const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
             return reduce_to_3x3(cl);
         },
-        py::return_value_policy::move, R"(Reduce cluster to 3x3 subcluster)");
+        py::return_value_policy::move, py::arg("cluster"), R"doc(
+        Return the 3x3 block around the cluster's center index.
+
+        Both input dimensions must be at least 3, and at least one must be
+        greater than 3.
+        The input coordinates are preserved and output data is stored in
+        row-major order.
+        )doc");
 }
 
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
@@ -92,16 +126,12 @@ void reduce_to_2x2(py::module &m) {
         [](const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
             return reduce_to_2x2(cl);
         },
-        py::return_value_policy::move,
-        R"(
-        Reduce cluster to 2x2 subcluster by taking the 2x2 subcluster with 
-        the highest photon energy.
+        py::return_value_policy::move, py::arg("cluster"), R"doc(
+        Return the highest-sum center-adjacent 2x2 block.
 
-        RETURN: 
-        
-        reduced cluster (cluster is filled in row major ordering starting at the top left. Thus for a max subcluster in the top left corner the photon hit is at the fourth position.)
-        
-        )");
+        The input coordinates are preserved and output data is stored in
+        row-major order.
+        )doc");
 }
 
-#pragma GCC diagnostic pop
+#pragma GCC diagnostic pop

python/src/bind_ClusterFile.hpp

@@ -11,6 +11,7 @@
 #include <pybind11/stl.h>
 #include <pybind11/stl/filesystem.h>
 #include <string>
+#include <utility>
 
 // Disable warnings for unused parameters, as we ignore some
 // in the __exit__ method
@@ -20,18 +21,90 @@
 namespace py = pybind11;
 using namespace ::aare;
 
+template <typename Range> class ClusterFileIterator {
+    Range m_range;
+    std::optional<typename Range::Iterator> m_iterator;
+    bool m_done{false};
+
+  public:
+    explicit ClusterFileIterator(Range range) : m_range(std::move(range)) {}
+
+    typename Range::Iterator::value_type next() {
+        if (m_done) {
+            throw py::stop_iteration();
+        }
+        // Advance before returning the next result, never after yielding it.
+        if (m_iterator) {
+            ++*m_iterator;
+        } else {
+            m_iterator.emplace(m_range.begin());
+        }
+        if (*m_iterator == m_range.end()) {
+            m_done = true;
+            throw py::stop_iteration();
+        }
+        return std::move(**m_iterator);
+    }
+};
+
+template <typename Iterator>
+void define_cluster_file_iterator(py::module &m, const std::string &name) {
+    py::class_<Iterator>(m, name.c_str())
+        .def(
+            "__iter__", [](Iterator &self) -> Iterator & { return self; },
+            py::return_value_policy::reference_internal)
+        .def("__next__", &Iterator::next);
+}
+
 template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
           typename CoordType = uint16_t>
 void define_ClusterFile(py::module &m, const std::string &typestr) {
 
     using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
+    using File = ClusterFile<ClusterType>;
+    using FrameIterator = ClusterFileIterator<typename File::FrameRange>;
+    using ChunkIterator = ClusterFileIterator<typename File::ChunkRange>;
 
     auto class_name = fmt::format("ClusterFile_{}", typestr);
+    define_cluster_file_iterator<FrameIterator>(m,
+                                                class_name + "_FrameIterator");
+    define_cluster_file_iterator<ChunkIterator>(m,
+                                                class_name + "_ChunkIterator");
 
-    py::class_<ClusterFile<ClusterType>>(m, class_name.c_str())
+    py::class_<ClusterFile<ClusterType>>(
+        m, class_name.c_str(),
+        "Read and write legacy binary cluster files. The format contains no "
+        "cluster type or shape metadata, so this class must match the file.")
         .def(py::init<const std::filesystem::path &, size_t,
                       const std::string &>(),
-             py::arg(), py::arg("chunk_size") = 1000, py::arg("mode") = "r")
+             py::arg("fname"), py::arg("chunk_size") = 1000,
+             py::arg("mode") = "r",
+             "Open a cluster file. Mode must be 'r' to read, 'w' to truncate "
+             "and write, or 'a' to append.")
+        .def(
+            "frames", [](File &self) { return FrameIterator(self.frames()); },
+            py::keep_alive<0, 1>(),
+            "Iterate over complete frames from the current position, including "
+            "empty and fully filtered frames with their stored frame numbers. "
+            "Results own their storage. The iterator keeps the file alive but "
+            "shares its cursor; use only one traversal at a time.")
+        .def(
+            "chunks", [](File &self) { return ChunkIterator(self.chunks()); },
+            py::keep_alive<0, 1>(),
+            "Iterate from the current position using the constructor's chunk "
+            "size, which must be positive. Results own their storage. The "
+            "iterator keeps the file alive but shares its cursor; use only "
+            "one traversal at a time.")
+        .def(
+            "chunks",
+            [](File &self, size_t chunk_size) {
+                return ChunkIterator(self.chunks(chunk_size));
+            },
+            py::arg("chunk_size"), py::keep_alive<0, 1>(),
+            "Iterate over chunks of up to chunk_size selected clusters. The "
+            "size must be positive and does not change the constructor's "
+            "default. Chunks may span frames; their frame numbers are not "
+            "per-cluster metadata. Only the final chunk can be short.")
         .def(
             "read_clusters",
             [](ClusterFile<ClusterType> &self, size_t n_clusters) {
@@ -39,32 +112,60 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
                     self.read_clusters(n_clusters));
                 return v;
             },
-            py::return_value_policy::take_ownership, py::arg("n_clusters"))
-        .def("read_frame",
-             [](ClusterFile<ClusterType> &self) {
-                 auto v = new ClusterVector<ClusterType>(self.read_frame());
-                 return v;
-             })
-        .def("set_roi", &ClusterFile<ClusterType>::set_roi, py::arg("roi"))
-        .def("tell", &ClusterFile<ClusterType>::tell)
+            py::return_value_policy::take_ownership, py::arg("n_clusters"),
+            "Read up to n_clusters without preserving frame boundaries. The "
+            "result may combine frames, so its frame number is not reliable "
+            "per-cluster metadata.")
+        .def(
+            "read_frame",
+            [](ClusterFile<ClusterType> &self) -> py::object {
+                auto clusters = self.read_frame();
+                if (!clusters) {
+                    return py::none();
+                }
+                return py::cast(
+                    new ClusterVector<ClusterType>(std::move(*clusters)),
+                    py::return_value_policy::take_ownership);
+            },
+            "Read and return the next complete frame with its frame number, "
+            "or None at end of file.")
+        .def("set_roi", &ClusterFile<ClusterType>::set_roi, py::arg("roi"),
+             "Select clusters whose centers lie within the half-open ROI.")
+        .def("tell", &ClusterFile<ClusterType>::tell,
+             "Return the current byte position in the file.")
         .def("estimate_n_clusters",
-             &ClusterFile<ClusterType>::estimate_n_clusters)
+             &ClusterFile<ClusterType>::estimate_n_clusters,
+             "Estimate the number of clusters from the file size. Frame "
+             "headers can make this larger than the actual count.")
         .def(
             "set_noise_map",
-            [](ClusterFile<ClusterType> &self, py::array_t<int32_t> noise_map) {
+            [](ClusterFile<ClusterType> &self,
+               py::array_t<int32_t, py::array::c_style> noise_map) {
                 auto view = make_view_2d(noise_map);
                 self.set_noise_map(view);
             },
-            py::arg("noise_map"))
+            py::arg("noise_map"),
+            "Set a two-dimensional, C-contiguous int32 noise map indexed as "
+            "[y, x]. The map must cover every cluster center coordinate.")
 
-        .def("set_gain_map",
-             [](ClusterFile<ClusterType> &self, py::array_t<double> gain_map) {
-                 auto view = make_view_2d(gain_map);
-                 self.set_gain_map(view);
-             })
+        .def(
+            "set_gain_map",
+            [](ClusterFile<ClusterType> &self,
+               py::array_t<double, py::array::c_style> gain_map) {
+                auto view = make_view_2d(gain_map);
+                self.set_gain_map(view);
+            },
+            py::arg("gain_map"),
+            "Set a two-dimensional, C-contiguous float64 gain map in "
+            "ADU/energy, indexed as [y, x]. Clusters whose complete footprint "
+            "extends beyond the map are retained with all data values set to "
+            "zero.")
 
-        .def("close", &ClusterFile<ClusterType>::close)
-        .def("write_frame", &ClusterFile<ClusterType>::write_frame)
+        .def("close", &ClusterFile<ClusterType>::close,
+             "Close the file. Calling close more than once is safe.")
+        .def("write_frame", &ClusterFile<ClusterType>::write_frame,
+             py::arg("clusters"),
+             "Write one ClusterVector, including its frame number.")
         .def("__enter__", [](ClusterFile<ClusterType> &self) { return &self; })
         .def("__exit__",
              [](ClusterFile<ClusterType> &self,
@@ -74,14 +175,8 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
                  self.close();
              })
         .def("__iter__", [](ClusterFile<ClusterType> &self) { return &self; })
-        .def("__next__", [](ClusterFile<ClusterType> &self) {
-            auto v = new ClusterVector<ClusterType>(
-                self.read_clusters(self.chunk_size()));
-            if (v->size() == 0) {
-                throw py::stop_iteration();
-            }
-            return v;
-        });
+        .def("__next__",
+             [](File &self) { return ChunkIterator(self.chunks()).next(); });
 }
 
 #pragma GCC diagnostic pop

python/src/bind_ClusterFinder.hpp

@@ -32,8 +32,10 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
 
     py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>(
         m, class_name.c_str())
-        .def(py::init<Shape<2>, pd_type, size_t>(), py::arg("image_size"),
-             py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
+        .def(py::init<Shape<2>, pd_type, size_t, size_t>(),
+             py::arg("image_size"), py::arg("n_sigma") = 5.0,
+             py::arg("capacity") = 1'000'000,
+             py::arg("min_pedestal_samples") = 1000)
 
         .def_property(
             "nSigma",

python/src/bind_ClusterFinderMT.hpp

@@ -32,10 +32,11 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
 
     py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>(
         m, class_name.c_str())
-        .def(py::init<Shape<2>, pd_type, size_t, size_t, size_t>(),
+        .def(py::init<Shape<2>, pd_type, size_t, size_t, size_t, size_t>(),
              py::arg("image_size"), py::arg("n_sigma") = 5.0,
              py::arg("capacity") = 2048, py::arg("n_threads") = 3,
-             py::arg("queue_depth") = 16)
+             py::arg("queue_depth") = 16,
+             py::arg("min_pedestal_samples") = 1000)
         .def("push_pedestal_frame",
              [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
                 py::array_t<uint16_t> frame) {

python/src/bind_ClusterVector.hpp

@@ -29,37 +29,73 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
 
     py::class_<ClusterVector<
         Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>, void>>(
-        m, class_name.c_str(), py::buffer_protocol(),
-        py::module_local())
+        m, class_name.c_str(), R"doc(
+        A contiguous, move-only container of fixed-size clusters.
 
-        .def(py::init()) // TODO change!!!
+        The class supports the Python buffer protocol, so ``numpy.array`` can
+        either copy its data or create a zero-copy view. A zero-copy view is
+        valid only while the ClusterVector's underlying allocation and size
+        remain unchanged.
+        )doc",
+        py::buffer_protocol(), py::module_local())
+
+        .def(py::init(), R"doc(
+        Create an empty ClusterVector with frame number 0 and space reserved
+        for at least 1024 clusters.
+        )doc")
 
         .def(
             "__call__",
-            [](ClusterVector<ClusterType> &self, py::array_t<bool> mask) {
+            [](ClusterVector<ClusterType> &self,
+               py::array_t<bool, py::array::c_style> mask) {
                 return self(make_view_1d(mask));
             },
-            py::arg("mask"), R"(
-            Create a copy of the clustervector and apply a boolean mask to the ClusterVector.
+            py::arg("mask").noconvert(), R"doc(
+            Return a filtered copy of this ClusterVector.
 
             Parameters
             ----------
+            mask : numpy.ndarray
+                One-dimensional, writable, C-contiguous array with dtype
+                ``numpy.bool_`` and one element per cluster.
 
-            mask : 1d boolean numpy array
-                Mask to apply to the ClusterVector. Must be the same length as the number of clusters in the ClusterVector. 
+            Returns
+            -------
+            ClusterVector
+                Selected clusters in their original order. The frame number is
+                preserved.
+            )doc")
 
-            )")
+        .def(
+            "push_back",
+            [](ClusterVector<ClusterType> &self, const ClusterType &cluster) {
+                self.push_back(cluster);
+            },
+            py::arg("cluster"), R"doc(
+            Append one cluster.
 
-        .def("push_back",
-             [](ClusterVector<ClusterType> &self, const ClusterType &cluster) {
-                 self.push_back(cluster);
-             })
+            Notes
+            -----
+            Do not call this method while a zero-copy NumPy view of the
+            ClusterVector exists. Reallocation invalidates the view, while an
+            append without reallocation leaves its shape unchanged.
+            )doc")
 
-        .def("sum",
-             [](ClusterVector<ClusterType> &self) {
-                 auto *vec = new std::vector<Type>(self.sum());
-                 return return_vector(vec);
-             })
+        .def(
+            "sum",
+            [](ClusterVector<ClusterType> &self) {
+                auto *vec = new std::vector<Type>(self.sum());
+                return return_vector(vec);
+            },
+            R"doc(
+            Return the sum of all pixels in each cluster.
+
+            Returns
+            -------
+            numpy.ndarray
+                One value per cluster, in container order and with the cluster
+                pixel dtype.
+            )doc")
         .def(
             "sum_2x2",
             [](ClusterVector<ClusterType> &self) {
@@ -68,24 +104,45 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
 
                 return return_vector(vec);
             },
-            R"(calculates sum of 2x2 subcluster with highest energy and index relative to cluster center 0: top_left, 1: top_right, 2: bottom_left, 3: bottom_right
-            )")
-        .def_property_readonly("size", &ClusterVector<ClusterType>::size)
-        .def("empty", &ClusterVector<ClusterType>::empty)
-        .def("item_size", &ClusterVector<ClusterType>::item_size)
-        .def_property_readonly("fmt",
-                               [typestr](ClusterVector<ClusterType> &self) {
-                                   return fmt_format<ClusterType>;
-                               })
+            R"doc(
+            Return the highest-sum center-adjacent 2x2 subcluster for each
+            cluster.
+
+            Returns
+            -------
+            numpy.ndarray
+                Structured array with ``sum`` and ``index`` fields. Indices are
+                0 for top-left, 1 for top-right, 2 for bottom-left, and 3 for
+                bottom-right, relative to the cluster center.
+            )doc")
+        .def_property_readonly("size", &ClusterVector<ClusterType>::size,
+                               "Number of stored clusters.")
+        .def("empty", &ClusterVector<ClusterType>::empty,
+             "Return True when no clusters are stored.")
+        .def("item_size", &ClusterVector<ClusterType>::item_size,
+             "Return the size in bytes of one stored cluster, including "
+             "padding.")
+        .def_property_readonly(
+            "fmt",
+            [typestr](ClusterVector<ClusterType> &self) {
+                return fmt_format<ClusterType>;
+            },
+            "PEP 3118 format string for one stored cluster.")
 
         .def_property_readonly("cluster_size_x",
-                               &ClusterVector<ClusterType>::cluster_size_x)
+                               &ClusterVector<ClusterType>::cluster_size_x,
+                               "Cluster size in the x dimension.")
         .def_property_readonly("cluster_size_y",
-                               &ClusterVector<ClusterType>::cluster_size_y)
+                               &ClusterVector<ClusterType>::cluster_size_y,
+                               "Cluster size in the y dimension.")
         .def_property_readonly("capacity",
-                               &ClusterVector<ClusterType>::capacity)
+                               &ClusterVector<ClusterType>::capacity,
+                               "Number of clusters that fit without "
+                               "reallocation.")
         .def_property("frame_number", &ClusterVector<ClusterType>::frame_number,
-                      &ClusterVector<ClusterType>::set_frame_number)
+                      &ClusterVector<ClusterType>::set_frame_number,
+                      "Signed 32-bit frame number; 0 can indicate clusters "
+                      "from multiple frames.")
         .def_buffer(
             [typestr](ClusterVector<ClusterType> &self) -> py::buffer_info {
                 return py::buffer_info(
@@ -99,31 +156,38 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
             });
 
     // Free functions using ClusterVector
-    m.def("hitmap",
-          [](std::array<size_t, 2> image_size, ClusterVector<ClusterType> &cv) {
-              // Create a numpy array to hold the hitmap
-              // The shape of the array is (image_size[0], image_size[1])
-              // note that the python array is passed as [row, col] which
-              // is the opposite of the clusters [x,y]
-              py::array_t<int32_t> hitmap(image_size);
-              auto r = hitmap.mutable_unchecked<2>();
+    m.def(
+        "hitmap",
+        [](std::array<size_t, 2> image_size, ClusterVector<ClusterType> &cv) {
+            // Create a numpy array to hold the hitmap
+            // The shape of the array is (image_size[0], image_size[1])
+            // note that the python array is passed as [row, col] which
+            // is the opposite of the clusters [x,y]
+            py::array_t<int32_t> hitmap(image_size);
+            auto r = hitmap.mutable_unchecked<2>();
 
-              // Initialize hitmap to 0
-              for (py::ssize_t i = 0; i < r.shape(0); i++)
-                  for (py::ssize_t j = 0; j < r.shape(1); j++)
-                      r(i, j) = 0;
+            // Initialize hitmap to 0
+            for (py::ssize_t i = 0; i < r.shape(0); i++)
+                for (py::ssize_t j = 0; j < r.shape(1); j++)
+                    r(i, j) = 0;
 
-              // Loop over the clusters and increment the hitmap
-              // Skip out of bound clusters
-              for (const auto &cluster : cv) {
-                  auto x = cluster.x;
-                  auto y = cluster.y;
-                  if (x < image_size[1] && y < image_size[0])
-                      r(cluster.y, cluster.x) += 1;
-              }
+            // Loop over the clusters and increment the hitmap
+            // Skip out of bound clusters
+            for (const auto &cluster : cv) {
+                auto x = cluster.x;
+                auto y = cluster.y;
+                if (x < image_size[1] && y < image_size[0])
+                    r(cluster.y, cluster.x) += 1;
+            }
 
-              return hitmap;
-          });
+            return hitmap;
+        },
+        py::arg("image_size"), py::arg("clusters"), R"doc(
+        Count cluster centers into an ``int32`` image whose shape is given by
+        ``image_size`` as ``(rows, columns)``. Element ``[y, x]`` contains the
+        number of cluster centers at that coordinate. Out-of-bounds centers are
+        ignored.
+        )doc");
 }
 
 template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
@@ -136,15 +200,12 @@ void define_2x2_reduction(py::module &m) {
             return new ClusterVector<Cluster<Type, 2, 2, CoordType>>(
                 reduce_to_2x2(cv));
         },
-        R"(
-        Reduce cluster to 2x2 subcluster by taking the 2x2 subcluster with 
-        the highest photon energy.
-
-        Parameters
-
-        cv : ClusterVector (clusters are filled in row-major ordering starting at the top left. Thus for a max subcluster in the top left corner the photon hit is at the fourth position.)
-
-        )",
+        R"doc(
+        Reduce every cluster to its highest-sum center-adjacent 2x2 block.
+        Returns a new ClusterVector; cluster order, coordinates, and frame
+        number are preserved. Input pixel data is interpreted in row-major
+        order.
+        )doc",
         py::arg("clustervector"));
 }
 
@@ -159,13 +220,13 @@ void define_3x3_reduction(py::module &m) {
             return new ClusterVector<Cluster<Type, 3, 3, CoordType>>(
                 reduce_to_3x3(cv));
         },
-        R"(
-        Reduce cluster to 3x3 subcluster 
-
-        Parameters
-        
-        cv : ClusterVector   
-        )",
+        R"doc(
+        Reduce every cluster to the 3x3 block around its center index.
+        Both input dimensions must be at least 3, and at least one must be
+        greater than 3.
+        Returns a new ClusterVector; cluster order, coordinates, and frame
+        number are preserved.
+        )doc",
         py::arg("clustervector"));
 }
 

python/src/bind_Defs.hpp

@@ -27,4 +27,77 @@ void define_defs_bindings(py::module &m) {
     moench05.attr("nRows") = Moench05::nRows;
     moench05.attr("nCols") = Moench05::nCols;
     moench05.attr("adcNumbers") = Moench05::adcNumbers;
+
+    py::class_<ROI>(m, "ROI")
+        .def(py::init<>())
+        .def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"),
+             py::arg("xmax"), py::arg("ymin"), py::arg("ymax"))
+        .def_readwrite("xmin", &ROI::xmin)
+        .def_readwrite("xmax", &ROI::xmax)
+        .def_readwrite("ymin", &ROI::ymin)
+        .def_readwrite("ymax", &ROI::ymax)
+
+        .def(
+            "slice",
+            [](const ROI &self) {
+                return std::make_tuple(py::slice(self.ymin, self.ymax, 1),
+                                       py::slice(self.xmin, self.xmax, 1));
+            },
+            R"doc(
+            The slice can be used to index into a 2D array to extract the region of interest.
+
+            Returns
+            -------
+            tuple of slice
+                A tuple containing two slice objects for the x and y dimensions.
+        )doc")
+
+        .def("__str__",
+             [](const ROI &self) {
+                 return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}",
+                                    self.xmin, self.xmax, self.ymin, self.ymax);
+             })
+        .def("__repr__",
+             [](const ROI &self) {
+                 return fmt::format(
+                     "<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin,
+                     self.xmax, self.ymin, self.ymax);
+             })
+        .def("__iter__",
+             [](const ROI &self) {
+                 return py::make_iterator(&self.xmin, &self.ymax + 1); // NOLINT
+             })
+
+        .def("__eq__", [](const ROI &self, const ROI &other) {
+            return self.xmin == other.xmin && self.xmax == other.xmax &&
+                   self.ymin == other.ymin && self.ymax == other.ymax;
+        });
+
+    py::enum_<DetectorType>(m, "DetectorType")
+        .value("Jungfrau", DetectorType::Jungfrau)
+        .value("Eiger", DetectorType::Eiger)
+        .value("Mythen3", DetectorType::Mythen3)
+        .value("Moench", DetectorType::Moench)
+        .value("Moench03", DetectorType::Moench03)
+        .value("Moench03_old", DetectorType::Moench03_old)
+        .value("ChipTestBoard", DetectorType::ChipTestBoard)
+        .value("Unknown", DetectorType::Unknown);
+
+    py::enum_<UDPPortPosition>(m, "UDPPortPosition")
+        .value("LEFT", UDPPortPosition::LEFT)
+        .value("RIGHT", UDPPortPosition::RIGHT)
+        .value("TOP", UDPPortPosition::TOP)
+        .value("BOTTOM", UDPPortPosition::BOTTOM);
+
+    py::enum_<TimingMode>(m, "TimingMode")
+        .value("AUTO_TIMING", TimingMode::AUTO_TIMING)
+        .value("TRIGGER_EXPOSURE", TimingMode::TRIGGER_EXPOSURE)
+        .value("GATED", TimingMode::GATED)
+        .value("BURST_TRIGGER", TimingMode::BURST_TRIGGER)
+        .value("TRIGGER_GATED", TimingMode::TRIGGER_GATED);
+
+    py::enum_<FrameDiscardPolicy>(m, "FrameDiscardPolicy")
+        .value("NoDiscard", FrameDiscardPolicy::NoDiscard)
+        .value("Discard", FrameDiscardPolicy::Discard)
+        .value("DiscardPartial", FrameDiscardPolicy::DiscardPartial);
 }

python/src/bind_Eta.hpp

@@ -1,6 +1,8 @@
 #include "aare/CalculateEta.hpp"
+#include "aare/to_string.hpp"
 
 #include <cstdint>
+#include <fmt/format.h>
 // #include <pybind11/native_enum.h> only for version 3
 #include <pybind11/pybind11.h>
 
@@ -18,7 +20,11 @@ void define_eta(py::module &m, const std::string &typestr) {
         .def_readwrite("c", &Eta2<T>::c,
                        "eta corner value cTopLeft, cTopRight, "
                        "cBottomLeft, cBottomRight")
-        .def_readwrite("sum", &Eta2<T>::sum, "photon energy of cluster");
+        .def_readwrite("sum", &Eta2<T>::sum, "photon energy of cluster")
+        .def("__repr__", [class_name](const Eta2<T> &self) {
+            return fmt::format("{}(x={}, y={}, c={}, sum={})", class_name,
+                               self.x, self.y, to_string(self.c), self.sum);
+        });
 }
 
 void define_corner_enum(py::module &m) {
@@ -101,4 +107,4 @@ void register_calculate_3x3eta(py::module &m) {
         },
         R"(calculates eta3x3 taking into account cross pixels in cluster)",
         py::arg("cluster"));
-}
+}

python/src/bind_FastPedestal.hpp

@@ -15,8 +15,7 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) {
 
     py::class_<FastPedestal<SUM_TYPE>>(
         m, name.c_str(),
-        "Maintain a per-pixel running mean, population variance, and "
-        "standard deviation.",
+        "Maintain a per-pixel running mean and population standard deviation.",
         py::buffer_protocol())
         .def(py::init<uint32_t, uint32_t, uint32_t>(), py::arg("rows"),
              py::arg("cols"), py::arg("n_samples"),
@@ -33,15 +32,6 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) {
                 return return_image_data(mean);
             },
             "Return a copy of the cached mean. The pedestal must be ready.")
-        .def(
-            "var",
-            [](FastPedestal<SUM_TYPE> &self) {
-                auto variance = new NDArray<SUM_TYPE, 2>{};
-                *variance = self.variance();
-                return return_image_data(variance);
-            },
-            "Return the population variance, normalized by n_samples, as a "
-            "NumPy array. The pedestal must be ready.")
         .def(
             "std",
             [](FastPedestal<SUM_TYPE> &self) {

python/src/bind_Pedestal.hpp

@@ -0,0 +1,132 @@
+// SPDX-License-Identifier: MPL-2.0
+
+#include "aare/Pedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+template <typename SUM_TYPE>
+void define_pedestal_bindings(py::module &m, const std::string &name) {
+
+    py::class_<Pedestal<SUM_TYPE>>(
+        m, name.c_str(),
+        "Maintain a per-pixel running mean and population standard deviation. "
+        "Statistics are available during initialization and are zero for empty "
+        "pixels.",
+        py::buffer_protocol())
+        .def(py::init<uint32_t, uint32_t, uint32_t>(), py::arg("rows"),
+             py::arg("cols"), py::arg("n_samples"),
+             "Construct an empty pedestal. Each pixel accumulates n_samples "
+             "values before switching to exponential updates.")
+        .def(py::init<uint32_t, uint32_t>(), py::arg("rows"), py::arg("cols"),
+             "Construct an empty pedestal with n_samples=1000.")
+        .def(
+            "mean",
+            [](Pedestal<SUM_TYPE> &self) {
+                auto mea = new NDArray<SUM_TYPE, 2>{};
+                *mea = self.mean();
+                return return_image_data(mea);
+            },
+            "Return a copy of the cached mean. Empty pixels return zero.")
+        .def(
+            "view",
+            [](py::object self_py) {
+                return py::module_::import("numpy").attr("asarray")(self_py);
+            },
+            "Return a non-owning, non-writable NumPy view of the cached mean.")
+        .def(
+            "std",
+            [](Pedestal<SUM_TYPE> &self) {
+                auto std = new NDArray<SUM_TYPE, 2>{};
+                *std = self.std();
+                return return_image_data(std);
+            },
+            "Return the population standard deviation as a NumPy array. "
+            "Empty pixels return zero.")
+        .def(
+            "__array_ufunc__",
+            [](py::object self, py::object ufunc, const std::string &method,
+               py::args inputs, py::kwargs kwargs) -> py::object {
+                if (method != "__call__" || inputs.size() != 2 ||
+                    inputs[1].ptr() != self.ptr() ||
+                    py::cast<std::string>(ufunc.attr("__name__")) !=
+                        "subtract") {
+                    return py::reinterpret_borrow<py::object>(
+                        Py_NotImplemented);
+                }
+
+                auto mean =
+                    py::module_::import("builtins").attr("memoryview")(self);
+                return ufunc(inputs[0], mean, **kwargs);
+            },
+            "Support subtracting a Pedestal from a NumPy array.")
+        .def("clear", py::overload_cast<>(&Pedestal<SUM_TYPE>::clear),
+             "Reset all statistics and per-pixel sample counts to zero.")
+        .def_property_readonly("rows", &Pedestal<SUM_TYPE>::rows,
+                               "Number of image rows.")
+        .def_property_readonly("cols", &Pedestal<SUM_TYPE>::cols,
+                               "Number of image columns.")
+        .def_property_readonly(
+            "n_samples", &Pedestal<SUM_TYPE>::n_samples,
+            "Initialization sample count per pixel and steady-state "
+            "update-weight denominator.")
+        .def(
+            "clone",
+            [&](Pedestal<SUM_TYPE> &pedestal) {
+                return Pedestal<SUM_TYPE>(pedestal);
+            },
+            "Return an independent copy of the pedestal and its state.")
+        // TODO! add push for other data types
+        .def(
+            "push",
+            [](Pedestal<SUM_TYPE> &pedestal,
+               py::array_t<uint16_t, py::array::c_style> &f) {
+                if (f.ndim() != 2) {
+                    throw py::value_error("Frame must be 2-dimensional");
+                }
+                auto v = make_view_2d(f);
+                pedestal.push(v);
+            },
+            py::arg("frame").noconvert(),
+            "Accumulate or exponentially update every pixel from a "
+            "C-contiguous uint16 frame matching the pedestal shape. After "
+            "n_samples values per pixel, new values have weight 1 / n_samples.")
+        .def(
+            "push_with_threshold",
+            [](Pedestal<SUM_TYPE> &pedestal,
+               py::array_t<uint16_t, py::array::c_style> &f,
+               py::array_t<SUM_TYPE, py::array::c_style> &threshold) {
+                if (f.ndim() != 2) {
+                    throw py::value_error("Frame must be 2-dimensional");
+                }
+                if (threshold.ndim() != 2) {
+                    throw py::value_error("Threshold must be 2-dimensional");
+                }
+                auto frame_view = make_view_2d(f);
+                auto threshold_view = make_view_2d(threshold);
+                pedestal.push_with_threshold(frame_view, threshold_view);
+            },
+            py::arg("frame").noconvert(), py::arg("threshold").noconvert(),
+            "Push only pixels where abs(frame - mean) is strictly less than "
+            "threshold. Both arrays must be C-contiguous with the pedestal "
+            "shape; frame must be uint16 and threshold must use the output "
+            "dtype. Rejected pixels keep their statistics and sample counts.")
+        .def_buffer([](Pedestal<SUM_TYPE> &self) {
+            auto mean = self.view();
+            return py::buffer_info(
+                const_cast<SUM_TYPE *>(mean.data()), sizeof(SUM_TYPE),
+                py::format_descriptor<SUM_TYPE>::format(), 2,
+                {static_cast<py::ssize_t>(mean.shape(0)),
+                 static_cast<py::ssize_t>(mean.shape(1))},
+                {static_cast<py::ssize_t>(mean.strides()[0] * sizeof(SUM_TYPE)),
+                 static_cast<py::ssize_t>(mean.strides()[1] *
+                                          sizeof(SUM_TYPE))},
+                true);
+        });
+}

python/src/bind_RawFile.hpp

@@ -11,6 +11,7 @@
 
 #include <cstdint>
 #include <filesystem>
+#include <fmt/format.h>
 #include <pybind11/iostream.h>
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
@@ -21,14 +22,23 @@
 namespace py = pybind11;
 using namespace ::aare;
 
+namespace {
+std::runtime_error raw_file_read_error(RawFile &file,
+                                       const std::string &message) {
+    return std::runtime_error(fmt::format(
+        "Error reading frame index {} from file '{}': {}", file.tell(),
+        file.master().master_fname().string(), message));
+}
+} // namespace
+
 void define_raw_file_io_bindings(py::module &m) {
     py::class_<RawFile>(m, "RawFile")
         .def(py::init<const std::filesystem::path &>())
         .def("read_frame",
              [](RawFile &self) {
                  if (self.n_modules_in_roi().size() > 1) {
-                     throw std::runtime_error(
-                         "File contains multiple ROIs - use read_ROIs()");
+                     throw raw_file_read_error(
+                         self, "File contains multiple ROIs - use read_ROIs()");
                  }
 
                  std::vector<size_t> shape;
@@ -51,16 +61,19 @@ void define_raw_file_io_bindings(py::module &m) {
             "read_n",
             [](RawFile &self, size_t n_frames) {
                 if (self.n_modules_in_roi().size() > 1) {
-                    throw std::runtime_error(
+                    throw raw_file_read_error(
+                        self,
                         "File contains multiple ROIs - use read_n_ROIs() to "
                         "read a specific ROI or use read_ROIs and "
                         "read one frame at a time.");
                 }
                 // adjust for actual frames left in the file
                 n_frames =
-                    std::min(n_frames, self.total_frames() - self.tell());
+                    std::min(n_frames, self.tell() < self.total_frames()
+                                           ? self.total_frames() - self.tell()
+                                           : size_t{0});
                 if (n_frames == 0) {
-                    throw std::runtime_error("No frames left in file");
+                    throw raw_file_read_error(self, "No frames left in file");
                 }
                 std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
 
@@ -89,12 +102,13 @@ void define_raw_file_io_bindings(py::module &m) {
             "read_roi",
             [](RawFile &self, const size_t roi_index) {
                 if (self.num_rois() == 0) {
-                    throw std::runtime_error(LOCATION + "No ROIs defined.");
+                    throw raw_file_read_error(self,
+                                              LOCATION + "No ROIs defined.");
                 }
 
                 if (roi_index >= self.num_rois()) {
-                    throw std::runtime_error(LOCATION +
-                                             "ROI index out of range.");
+                    throw raw_file_read_error(
+                        self, LOCATION + "ROI index out of range.");
                 }
 
                 // return headers from all subfiles
@@ -141,7 +155,8 @@ void define_raw_file_io_bindings(py::module &m) {
             "read_rois",
             [](RawFile &self) {
                 if (self.num_rois() == 0) {
-                    throw std::runtime_error(LOCATION + "No ROIs defined.");
+                    throw raw_file_read_error(self,
+                                              LOCATION + "No ROIs defined.");
                 }
 
                 size_t number_of_ROIs = self.num_rois();
@@ -196,19 +211,22 @@ void define_raw_file_io_bindings(py::module &m) {
             "read_n_with_roi",
             [](RawFile &self, const size_t num_frames, const size_t roi_index) {
                 if (self.num_rois() == 0) {
-                    throw std::runtime_error(LOCATION + "No ROIs defined.");
+                    throw raw_file_read_error(self,
+                                              LOCATION + "No ROIs defined.");
                 }
 
                 if (roi_index >= self.num_rois()) {
-                    throw std::runtime_error(LOCATION +
-                                             "ROI index out of range.");
+                    throw raw_file_read_error(
+                        self, LOCATION + "ROI index out of range.");
                 }
 
                 // adjust for actual frames left in the file
                 size_t n_frames =
-                    std::min(num_frames, self.total_frames() - self.tell());
+                    std::min(num_frames, self.tell() < self.total_frames()
+                                             ? self.total_frames() - self.tell()
+                                             : size_t{0});
                 if (n_frames == 0) {
-                    throw std::runtime_error("No frames left in file");
+                    throw raw_file_read_error(self, "No frames left in file");
                 }
                 std::vector<size_t> shape{
                     n_frames, self.roi_geometries(roi_index).pixels_y(),
@@ -306,4 +324,4 @@ void define_raw_file_io_bindings(py::module &m) {
         .def_property_readonly("n_modules", &RawFile::n_modules)
         .def_property_readonly("n_modules_in_roi", &RawFile::n_modules_in_roi)
         .def_property_readonly("num_rois", &RawFile::num_rois);
-}
+}

python/src/bind_calibration.hpp

@@ -73,10 +73,27 @@ py::array_t<T> pybind_calculate_pedestal_g0(
 }
 
 void bind_calibration(py::module &m) {
+    py::options options;
+    options.disable_function_signatures();
+
     m.def("apply_calibration", &pybind_apply_calibration<double>,
           py::arg("raw_data").noconvert(), py::kw_only(),
           py::arg("pd").noconvert(), py::arg("cal").noconvert(),
-          py::arg("n_threads") = 4);
+          py::arg("n_threads") = 4,
+          R"(
+        Apply the calibration to the raw data and return the result as a 3D array of doubles.
+        
+        Parameters
+        ----------
+        raw_data : array_like
+            3D array of shape ``(frames, rows, cols)``.
+        pd : array_like
+            Pedestal array.
+        cal : array_like
+            Calibration array.
+        n_threads : int, optional
+            Number of threads to use. Defaults to 4.
+        )");
 
     m.def("apply_calibration", &pybind_apply_calibration<float>,
           py::arg("raw_data").noconvert(), py::kw_only(),

python/src/bind_testing.hpp

@@ -0,0 +1,65 @@
+// SPDX-License-Identifier: MPL-2.0
+// Probes used by the python test-suite to exercise the helpers in
+// np_helper.hpp. Internal, no API guarantees.
+#pragma once
+
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "np_helper.hpp"
+
+#include <algorithm>
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <type_traits>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+// Arrays are taken with the default flags (no c_style) so that nothing is
+// made contiguous by pybind11 before it reaches make_view
+template <typename T, ssize_t Ndim>
+void define_make_view_probes(py::module &m, const std::string &suffix) {
+
+    // Copy what C++ sees through the view back to python
+    m.def(("view_roundtrip_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_view<Ndim>(arr);
+              return return_image_data(new NDArray<T, Ndim>(view));
+          },
+          py::arg("arr").noconvert());
+
+    // Same as above but through a read only view
+    m.def(("const_view_roundtrip_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_const_view<Ndim>(arr);
+              static_assert(
+                  std::is_same_v<decltype(view), NDView<const T, Ndim>>);
+              auto out = new NDArray<T, Ndim>(view.shape());
+              std::copy(view.begin(), view.end(), out->begin());
+              return return_image_data(out);
+          },
+          py::arg("arr").noconvert());
+
+    // Write through the view
+    m.def(("view_fill_" + suffix).c_str(),
+          [](py::array_t<T> arr, T value) {
+              auto view = make_view<Ndim>(arr);
+              view = value;
+          },
+          py::arg("arr").noconvert(), py::arg("value"));
+}
+
+void define_testing_bindings(py::module &m) {
+    auto testing = m.def_submodule(
+        "testing", "Internal helpers for the test-suite, no API guarantees");
+
+    // noconvert + overloads gives dispatch on dtype
+    define_make_view_probes<double, 1>(testing, "1d");
+    define_make_view_probes<double, 2>(testing, "2d");
+    define_make_view_probes<double, 3>(testing, "3d");
+    define_make_view_probes<uint16_t, 1>(testing, "1d");
+    define_make_view_probes<uint16_t, 2>(testing, "2d");
+    define_make_view_probes<uint16_t, 3>(testing, "3d");
+}

python/src/file.hpp

@@ -2,6 +2,7 @@
 #include "aare/CtbRawFile.hpp"
 #include "aare/File.hpp"
 #include "aare/Frame.hpp"
+#include "aare/ROI.hpp"
 #include "aare/RawFile.hpp"
 #include "aare/RawMasterFile.hpp"
 #include "aare/RawSubFile.hpp"
@@ -27,16 +28,6 @@ using namespace ::aare;
 
 void define_file_io_bindings(py::module &m) {
 
-    py::enum_<DetectorType>(m, "DetectorType")
-        .value("Jungfrau", DetectorType::Jungfrau)
-        .value("Eiger", DetectorType::Eiger)
-        .value("Mythen3", DetectorType::Mythen3)
-        .value("Moench", DetectorType::Moench)
-        .value("Moench03", DetectorType::Moench03)
-        .value("Moench03_old", DetectorType::Moench03_old)
-        .value("ChipTestBoard", DetectorType::ChipTestBoard)
-        .value("Unknown", DetectorType::Unknown);
-
     PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber,
                          bunchId, timestamp, modId, row, column, reserved,
                          debug, roundRNumber, detType, version, packetMask);
@@ -170,28 +161,5 @@ void define_file_io_bindings(py::module &m) {
         .def_property_readonly("stop", &ScanParameters::stop)
         .def_property_readonly("step", &ScanParameters::step);
 
-    py::class_<ROI>(m, "ROI")
-        .def(py::init<>())
-        .def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"),
-             py::arg("xmax"), py::arg("ymin"), py::arg("ymax"))
-        .def_readwrite("xmin", &ROI::xmin)
-        .def_readwrite("xmax", &ROI::xmax)
-        .def_readwrite("ymin", &ROI::ymin)
-        .def_readwrite("ymax", &ROI::ymax)
-        .def("__str__",
-             [](const ROI &self) {
-                 return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}",
-                                    self.xmin, self.xmax, self.ymin, self.ymax);
-             })
-        .def("__repr__",
-             [](const ROI &self) {
-                 return fmt::format(
-                     "<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin,
-                     self.xmax, self.ymin, self.ymax);
-             })
-        .def("__iter__", [](const ROI &self) {
-            return py::make_iterator(&self.xmin, &self.ymax + 1); // NOLINT
-        });
-
 #pragma GCC diagnostic pop
 }

python/src/module.cpp

@@ -17,18 +17,19 @@
 #include "bind_FastPedestal.hpp"
 #include "bind_Interpolator.hpp"
 #include "bind_MultiThreadedFileReader.hpp"
+#include "bind_Pedestal.hpp"
 #include "bind_PedestalTrackingPixelHistogram.hpp"
 #include "bind_PixelHistogram.hpp"
 #include "bind_PixelMap.hpp"
 #include "bind_RawFile.hpp"
 #include "bind_calibration.hpp"
+#include "bind_testing.hpp"
 
 // TODO! migrate the other names
 #include "ctb_raw_file.hpp"
 #include "file.hpp"
 #include "fit.hpp"
 #include "jungfrau_data_file.hpp"
-#include "pedestal.hpp"
 #include "raw_master_file.hpp"
 #include "raw_sub_file.hpp"
 #include "var_cluster.hpp"
@@ -90,6 +91,7 @@ PYBIND11_MODULE(_aare, m) {
     define_jungfrau_data_file_io_bindings(m);
 
     bind_calibration(m);
+    define_testing_bindings(m);
 
     DEFINE_CLUSTER_BINDINGS(int, 3, 3, uint16_t, i);
     DEFINE_CLUSTER_BINDINGS(double, 3, 3, uint16_t, d);
@@ -130,9 +132,6 @@ PYBIND11_MODULE(_aare, m) {
     DEFINE_BINDINGS_CLUSTERFINDER(double, 9, 9, uint16_t, d);
     DEFINE_BINDINGS_CLUSTERFINDER(float, 9, 9, uint16_t, f);
 
-    define_3x3_reduction<int, 3, 3, uint16_t>(m);
-    define_3x3_reduction<double, 3, 3, uint16_t>(m);
-    define_3x3_reduction<float, 3, 3, uint16_t>(m);
     define_3x3_reduction<int, 5, 5, uint16_t>(m);
     define_3x3_reduction<double, 5, 5, uint16_t>(m);
     define_3x3_reduction<float, 5, 5, uint16_t>(m);
@@ -143,9 +142,6 @@ PYBIND11_MODULE(_aare, m) {
     define_3x3_reduction<double, 9, 9, uint16_t>(m);
     define_3x3_reduction<float, 9, 9, uint16_t>(m);
 
-    reduce_to_3x3<int, 3, 3, uint16_t>(m);
-    reduce_to_3x3<double, 3, 3, uint16_t>(m);
-    reduce_to_3x3<float, 3, 3, uint16_t>(m);
     reduce_to_3x3<int, 5, 5, uint16_t>(m);
     reduce_to_3x3<double, 5, 5, uint16_t>(m);
     reduce_to_3x3<float, 5, 5, uint16_t>(m);

python/src/np_helper.hpp

@@ -41,31 +41,60 @@ template <typename T> py::array return_vector(std::vector<T> *vec) {
                           free_when_done); // numpy array references this parent
 }
 
-// todo rewrite generic
-template <class T, int Flags>
-auto get_shape_3d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<3>{arr.shape(0), arr.shape(1), arr.shape(2)};
+// Create a NDView of a numpy array. NDView assumes C-order strides so we
+// reject anything that is not Ndim dimensional and C-contiguous instead of
+// silently reading the data in the wrong order.
+template <ssize_t Ndim, class T, int Flags>
+auto get_checked_shape(const py::array_t<T, Flags> &arr) {
+    if (arr.ndim() != Ndim) {
+        throw py::value_error(
+            fmt::format("Expected {}D array, got {}D", Ndim, arr.ndim()));
+    }
+    if (!(arr.flags() & py::array::c_style)) {
+        throw py::value_error(
+            "Array is not C-contiguous. Use np.ascontiguousarray(arr)");
+    }
+    aare::Shape<Ndim> shape{};
+    for (ssize_t i = 0; i < Ndim; ++i) {
+        shape[i] = arr.shape(i);
+    }
+    return shape;
+}
+
+template <ssize_t Ndim, class T, int Flags>
+auto make_view(py::array_t<T, Flags> &arr) {
+    auto shape = get_checked_shape<Ndim>(arr);
+    return aare::NDView<T, Ndim>(arr.mutable_data(), shape);
+}
+
+// Read only view, also works for numpy arrays that are not writeable
+template <ssize_t Ndim, class T, int Flags>
+auto make_const_view(const py::array_t<T, Flags> &arr) {
+    auto shape = get_checked_shape<Ndim>(arr);
+    return aare::NDView<const T, Ndim>(arr.data(), shape);
 }
 
 template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 3>(arr.mutable_data(), get_shape_3d<T, Flags>(arr));
+    return make_view<3>(arr);
 }
-
-template <class T, int Flags>
-auto get_shape_2d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<2>{arr.shape(0), arr.shape(1)};
-}
-
-template <class T, int Flags>
-auto get_shape_1d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<1>{arr.shape(0)};
-}
-
 template <class T, int Flags> auto make_view_2d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 2>(arr.mutable_data(), get_shape_2d<T, Flags>(arr));
+    return make_view<2>(arr);
 }
 template <class T, int Flags> auto make_view_1d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 1>(arr.mutable_data(), get_shape_1d<T, Flags>(arr));
+    return make_view<1>(arr);
+}
+
+template <class T, int Flags>
+auto make_const_view_3d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<3>(arr);
+}
+template <class T, int Flags>
+auto make_const_view_2d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<2>(arr);
+}
+template <class T, int Flags>
+auto make_const_view_1d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<1>(arr);
 }
 
 template <typename ClusterType> struct fmt_format_trait; // forward declaration

python/src/pedestal.hpp

@@ -1,125 +0,0 @@
-// SPDX-License-Identifier: MPL-2.0
-
-#include "aare/Pedestal.hpp"
-#include "np_helper.hpp"
-
-#include <cstdint>
-#include <filesystem>
-#include <pybind11/numpy.h>
-#include <pybind11/pybind11.h>
-#include <pybind11/stl.h>
-
-namespace py = pybind11;
-
-template <typename SUM_TYPE>
-void define_pedestal_bindings(py::module &m, const std::string &name) {
-
-    py::class_<Pedestal<SUM_TYPE>>(m, name.c_str(), py::buffer_protocol())
-        .def(py::init<int, int, int>())
-        .def(py::init<int, int>())
-        .def("mean",
-             [](Pedestal<SUM_TYPE> &self) {
-                 auto mea = new NDArray<SUM_TYPE, 2>{};
-                 *mea = self.mean();
-                 return return_image_data(mea);
-             })
-        .def("view",
-             [](py::object self_py) {
-                 auto &self = self_py.cast<Pedestal<SUM_TYPE> &>();
-                 auto v = self.view();
-                 std::array<py::ssize_t, 2> shape{
-                     static_cast<py::ssize_t>(v.shape(0)),
-                     static_cast<py::ssize_t>(v.shape(1))};
-                 std::array<py::ssize_t, 2> byte_strides{
-                     static_cast<py::ssize_t>(v.strides()[0]) *
-                         static_cast<py::ssize_t>(sizeof(SUM_TYPE)),
-                     static_cast<py::ssize_t>(v.strides()[1]) *
-                         static_cast<py::ssize_t>(sizeof(SUM_TYPE))};
-                 auto arr = py::array_t<SUM_TYPE>(shape, byte_strides, v.data(),
-                                                  self_py);
-                 arr.attr("setflags")(py::arg("write") = false);
-                 return arr;
-             })
-        .def("variance",
-             [](Pedestal<SUM_TYPE> &self) {
-                 auto var = new NDArray<SUM_TYPE, 2>{};
-                 *var = self.variance();
-                 return return_image_data(var);
-             })
-        .def("std",
-             [](Pedestal<SUM_TYPE> &self) {
-                 auto std = new NDArray<SUM_TYPE, 2>{};
-                 *std = self.std();
-                 return return_image_data(std);
-             })
-        .def("cached_std",
-             [](Pedestal<SUM_TYPE> &self) {
-                 auto standard_deviation = new NDArray<SUM_TYPE, 2>{};
-                 *standard_deviation = self.cached_std();
-                 return return_image_data(standard_deviation);
-             })
-        .def(
-            "__array_ufunc__",
-            [](py::object self, py::object ufunc, const std::string &method,
-               py::args inputs, py::kwargs kwargs) -> py::object {
-                if (method != "__call__" || inputs.size() != 2 ||
-                    inputs[1].ptr() != self.ptr() ||
-                    py::cast<std::string>(ufunc.attr("__name__")) !=
-                        "subtract") {
-                    return py::reinterpret_borrow<py::object>(
-                        Py_NotImplemented);
-                }
-
-                auto mean =
-                    py::module_::import("builtins").attr("memoryview")(self);
-                return ufunc(inputs[0], mean, **kwargs);
-            },
-            "Support subtracting a Pedestal from a NumPy array.")
-        .def("clear", py::overload_cast<>(&Pedestal<SUM_TYPE>::clear))
-        .def_property_readonly("rows", &Pedestal<SUM_TYPE>::rows)
-        .def_property_readonly("cols", &Pedestal<SUM_TYPE>::cols)
-        .def_property_readonly("n_samples", &Pedestal<SUM_TYPE>::n_samples)
-        .def_property_readonly("sum", &Pedestal<SUM_TYPE>::get_sum)
-        .def_property_readonly("sum2", &Pedestal<SUM_TYPE>::get_sum2)
-        .def("clone",
-             [&](Pedestal<SUM_TYPE> &pedestal) {
-                 return Pedestal<SUM_TYPE>(pedestal);
-             })
-        // TODO! add push for other data types
-        .def("push",
-             [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> &f) {
-                 auto v = make_view_2d(f);
-                 pedestal.push(v);
-             })
-        .def(
-            "push_with_threshold",
-            [](Pedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &f,
-               py::array_t<SUM_TYPE, py::array::c_style> &threshold) {
-                auto frame_view = make_view_2d(f);
-                auto threshold_view = make_view_2d(threshold);
-                pedestal.push_with_threshold(frame_view, threshold_view);
-            },
-            py::arg("frame").noconvert(), py::arg("threshold").noconvert())
-        .def(
-            "push_no_update",
-            [](Pedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &f) {
-                auto v = make_view_2d(f);
-                pedestal.push_no_update(v);
-            },
-            py::arg().noconvert())
-        .def("update_mean", &Pedestal<SUM_TYPE>::update_mean)
-        .def_buffer([](Pedestal<SUM_TYPE> &self) {
-            auto mean = self.view();
-            return py::buffer_info(
-                const_cast<SUM_TYPE *>(mean.data()), sizeof(SUM_TYPE),
-                py::format_descriptor<SUM_TYPE>::format(), 2,
-                {static_cast<py::ssize_t>(mean.shape(0)),
-                 static_cast<py::ssize_t>(mean.shape(1))},
-                {static_cast<py::ssize_t>(mean.strides()[0] * sizeof(SUM_TYPE)),
-                 static_cast<py::ssize_t>(mean.strides()[1] *
-                                          sizeof(SUM_TYPE))},
-                true);
-        });
-}

python/src/raw_master_file.hpp

@@ -66,7 +66,8 @@ void define_raw_master_file_bindings(py::module &m) {
 
         .def_property_readonly("total_frames_expected",
                                &RawMasterFile::total_frames_expected)
-        .def_property_readonly("geometry", &RawMasterFile::geometry)
+        .def_property_readonly("detector_layout",
+                               &RawMasterFile::detector_layout)
         .def_property_readonly("udp_interfaces_per_module",
                                &RawMasterFile::udp_interfaces_per_module)
         .def_property_readonly("analog_samples", &RawMasterFile::analog_samples,
@@ -107,9 +108,46 @@ void define_raw_master_file_bindings(py::module &m) {
                     return std::nullopt;
                 }
             })
-        .def_property_readonly("period", [](RawMasterFile &self) {
-            double seconds =
-                std::chrono::duration<double>(self.period()).count();
-            return seconds;
-        });
+        .def_property_readonly("rois", &RawMasterFile::rois, R"(
+            Get the ROIs defined in the master file
+
+            Returns
+            ----------
+                List[ROI]
+                    List of ROIs (default complete ROI)
+            )")
+        .def_property_readonly("udp_port_types", &RawMasterFile::udp_port_types,
+                               R"(
+            Get the types of UDP ports
+
+            Returns
+            ----------
+                Optional[List[UDPPortPosition]]
+                    Optional vector of UDP port types as strings (only present for
+                    masterfile version >= 8.1)
+            )")
+        .def_property_readonly("disabled_udp_ports",
+                               &RawMasterFile::disabled_udp_ports, R"(
+            Get the indices of disabled UDP ports   
+
+            Returns
+            ----------
+                List[int]
+                    Vector of disabled UDP port indices relative to UDP port types (empty if none are disabled)
+            )")
+        .def_property_readonly(
+            "period",
+            [](RawMasterFile &self) {
+                double seconds =
+                    std::chrono::duration<double>(self.period()).count();
+                return seconds;
+            })
+
+        .def("__enter__",
+             [](RawMasterFile &self) -> RawMasterFile & { return self; })
+
+        .def("__exit__",
+             []([[maybe_unused]] RawMasterFile &self,
+                [[maybe_unused]] py::object, [[maybe_unused]] py::object,
+                [[maybe_unused]] py::object) { return; });
 }

python/tests/test_Cluster.py

@@ -3,7 +3,7 @@ import pytest
 import numpy as np
 
 from aare import _aare #import the C++ module
-from aare import corner 
+from aare import ROI, corner 
 from conftest import test_data_path
 
 
@@ -91,6 +91,22 @@ def test_calculate_eta():
     assert eta2.sum == 4
 
 
+def test_cluster_representations():
+    """Test Cluster representation"""
+
+    cluster = _aare.Cluster2x2i(5, 6, np.array([1, 2, 3, 4], dtype=np.int32))
+    assert repr(cluster) == "Cluster2x2i(x=5, y=6, data=[[1, 2], [3, 4]])"
+
+def test_eta_representations():
+    """Test Eta representation"""
+    eta = _aare.Etai()
+    eta.x = 0.25
+    eta.y = 0.75
+    eta.c = corner.cBottomRight
+    eta.sum = 42
+    assert repr(eta) == "Etai(x=0.25, y=0.75, c=BottomRight, sum=42)"
+
+
 def test_max_sum(): 
     """Max 2x2 Sum"""
     cluster = _aare.Cluster3x3i(5,5,np.array([1, 1, 1, 2, 3, 1, 2, 2, 1], dtype=np.int32))
@@ -146,4 +162,3 @@ def test_3x3_reduction():
 
 
 
-

python/tests/test_ClusterFile.py

@@ -6,10 +6,279 @@ import boost_histogram as bh
 import time
 from pathlib import Path
 import pickle
+import struct
+import gc
+import weakref
 
-from aare import ClusterFile
+from aare import ClusterFile, ROI
 from conftest import test_data_path
 
+
+@pytest.fixture
+def iterator_file(tmp_path):
+    fname = tmp_path / "iteration.clust"
+    records = [
+        (-11, []),
+        (42, [(5, 6, 10), (7, 8, 20)]),
+        (103, []),
+        (104, [(5, 6, 30), (7, 8, 40)]),
+        (105, []),
+    ]
+    fname.write_bytes(b"".join(
+        struct.pack("@iI", number, len(clusters))
+        + b"".join(struct.pack("@HH9i", x, y, *([value] * 9))
+                   for x, y, value in clusters)
+        for number, clusters in records
+    ))
+    return fname
+
+
+def test_frames_preserve_empty_frames_and_frame_numbers(iterator_file):
+    with ClusterFile(iterator_file) as reader:
+        frames = reader.frames()
+        assert iter(frames) is frames
+        assert reader.tell() == 0
+        results = list(frames)
+        assert [frame.frame_number for frame in results] == [-11, 42, 103, 104, 105]
+        assert [frame.size for frame in results] == [0, 2, 0, 2, 0]
+        for _ in range(2):
+            with pytest.raises(StopIteration):
+                next(frames)
+        assert list(reader.frames()) == []
+
+
+@pytest.mark.parametrize("chunk_size, sizes", [(1, [1, 1, 1, 1]), (2, [2, 2]),
+                                             (3, [3, 1]), (10, [4])])
+@pytest.mark.parametrize("explicit", [False, True])
+def test_chunks_cross_frames(iterator_file, chunk_size, sizes, explicit):
+    with ClusterFile(iterator_file, chunk_size=chunk_size) as reader:
+        chunks = reader.chunks(chunk_size=chunk_size) if explicit else reader.chunks()
+        assert iter(chunks) is chunks
+        assert reader.tell() == 0
+        results = list(chunks)
+        assert [chunk.size for chunk in results] == sizes
+        data = np.concatenate([np.asarray(chunk) for chunk in results])
+        np.testing.assert_array_equal(data["x"], [5, 7, 5, 7])
+        np.testing.assert_array_equal(data["data"][:, 0, 0], [10, 20, 30, 40])
+        for _ in range(2):
+            with pytest.raises(StopIteration):
+                next(chunks)
+
+
+def test_chunk_override_preserves_default_iteration(iterator_file):
+    with ClusterFile(iterator_file, chunk_size=2) as reader:
+        assert iter(reader) is reader
+        assert next(reader.chunks(1)).size == 1
+        assert [chunk.size for chunk in reader] == [2, 1]
+
+
+@pytest.mark.parametrize("use_roi", [False, True])
+@pytest.mark.parametrize("use_noise", [False, True])
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_iterators_apply_filters_and_gain(iterator_file, use_roi, use_noise, method):
+    with ClusterFile(iterator_file, chunk_size=3) as reader:
+        if use_roi:
+            reader.set_roi(ROI(0, 6, 0, 12))
+        if use_noise:
+            reader.set_noise_map(np.full((12, 12), 15, dtype=np.int32))
+        reader.set_gain_map(np.full((12, 12), 2.0))
+        results = list(getattr(reader, method)())
+        if method == "frames":
+            assert [frame.frame_number for frame in results] == [-11, 42, 103, 104, 105]
+        expected = [value / 2 for x, value in [(5, 10), (7, 20), (5, 30), (7, 40)]
+                    if (not use_roi or x == 5) and (not use_noise or value > 15)]
+        data = np.concatenate([np.asarray(result) for result in results])
+        np.testing.assert_array_equal(data["data"][:, 0, 0], expected)
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_iterators_handle_all_rejected_clusters(iterator_file, method):
+    with ClusterFile(iterator_file) as reader:
+        reader.set_roi(ROI(0, 1, 0, 1))
+        results = list(getattr(reader, method)())
+        assert len(results) == (5 if method == "frames" else 0)
+        assert all(result.size == 0 for result in results)
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_iterator_keeps_file_alive(iterator_file, method):
+    reader = ClusterFile(iterator_file)
+    owner = weakref.ref(reader)
+    iterator = getattr(reader, method)()
+    del reader
+    gc.collect()
+    assert owner() is not None
+    results = list(iterator)
+    assert sum(result.size for result in results) == 4
+    del iterator
+    gc.collect()
+    assert owner() is None
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_retained_results_and_numpy_views_own_storage(iterator_file, method):
+    with ClusterFile(iterator_file, chunk_size=2) as reader:
+        reader.read_frame()  # Skip the initial empty frame.
+        iterator = getattr(reader, method)()
+        first = next(iterator)
+        array = np.asarray(first)
+        expected = array.copy()
+        results = list(iterator)
+        np.testing.assert_array_equal(np.asarray(first), expected)
+    del first, results, iterator, reader
+    gc.collect()
+    np.testing.assert_array_equal(array, expected)
+
+
+def test_frame_iteration_resumes_without_reading_ahead(iterator_file):
+    with ClusterFile(iterator_file) as reader:
+        reader.read_frame()
+        for frame in reader.frames():
+            assert frame.frame_number == 42
+            break
+        assert reader.read_frame().frame_number == 103
+        assert [frame.frame_number for frame in reader.frames()] == [104, 105]
+
+
+def test_partial_chunk_must_be_completed_before_frame_iteration(iterator_file):
+    with ClusterFile(iterator_file) as reader:
+        for chunk in reader.chunks(1):
+            assert chunk.size == 1
+            break
+        with pytest.raises(RuntimeError, match="clusters left"):
+            next(reader.frames())
+        assert reader.read_clusters(1).size == 1
+        assert [frame.frame_number for frame in reader.frames()] == [103, 104, 105]
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_empty_iterators(tmp_path, method):
+    fname = tmp_path / "empty.clust"
+    fname.touch()
+    with ClusterFile(fname) as reader:
+        iterator = getattr(reader, method)()
+        for _ in range(2):
+            with pytest.raises(StopIteration):
+                next(iterator)
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+@pytest.mark.parametrize("started", [False, True])
+def test_iterators_reject_closed_files(iterator_file, method, started):
+    reader = ClusterFile(iterator_file, chunk_size=1)
+    iterator = getattr(reader, method)()
+    if started:
+        next(iterator)
+    reader.close()
+    with pytest.raises(RuntimeError, match="not opened for reading"):
+        next(iterator)
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+@pytest.mark.parametrize("mode", ["w", "a"])
+def test_iterators_reject_writing_modes(tmp_path, method, mode):
+    with ClusterFile(tmp_path / "output.clust", mode=mode) as reader:
+        with pytest.raises(RuntimeError, match="not opened for reading"):
+            next(getattr(reader, method)())
+
+
+def test_chunk_sizes_must_be_positive(iterator_file):
+    with ClusterFile(iterator_file, chunk_size=0) as reader:
+        with pytest.raises(ValueError, match="greater than zero"):
+            reader.chunks()
+        with pytest.raises(ValueError, match="greater than zero"):
+            reader.chunks(0)
+        with pytest.raises(ValueError, match="greater than zero"):
+            next(reader)
+        with pytest.raises(TypeError):
+            reader.chunks(-1)
+        assert reader.tell() == 0
+        assert reader.read_clusters(0).size == 0
+        assert len(list(reader.frames())) == 5
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_iterators_defer_read_errors_until_next_result(tmp_path, method):
+    fname = tmp_path / "truncated.clust"
+    cluster = struct.pack("@HH9i", 5, 6, *range(9))
+    frame = struct.pack("@iI", 42, 1) + cluster
+    fname.write_bytes(frame + frame[:-1])
+    with ClusterFile(fname, chunk_size=1) as reader:
+        iterator = getattr(reader, method)()
+        assert next(iterator).size == 1
+        assert reader.tell() == len(frame)
+        with pytest.raises(RuntimeError):
+            next(iterator)
+
+
+@pytest.mark.parametrize(
+    "shape, dtype",
+    [(shape, dtype)
+     for shape in [(2, 2), (3, 3), (5, 5), (7, 7), (9, 9)]
+     for dtype in [np.int32, np.float32, np.float64]]
+    + [((3, 3), np.int16)],
+)
+def test_iterators_are_bound_for_each_cluster_type(tmp_path, shape, dtype):
+    fname = tmp_path / "typed_frames.clust"
+    values = list(range(shape[0] * shape[1]))
+    record = struct.pack("@HH" + np.dtype(dtype).char * len(values), 5, 6, *values)
+    fname.write_bytes((struct.pack("@iI", -42, 1) + record) * 3)
+    with ClusterFile(fname, cluster_size=shape, dtype=dtype) as reader:
+        frame = next(reader.frames())
+        assert frame.frame_number == -42
+        np.testing.assert_array_equal(np.asarray(frame)["data"].reshape(-1), values)
+        chunks = list(reader.chunks())
+        assert len(chunks) == 1
+        assert chunks[0].size == 2
+        np.testing.assert_array_equal(np.asarray(chunks[0])["data"].reshape(-1), values * 2)
+
+
+def test_read_frame_returns_none_at_eof(tmp_path):
+    fname = tmp_path / "empty.clust"
+    fname.touch()
+
+    with ClusterFile(fname) as f:
+        assert f.read_frame() is None
+
+
+def test_read_frame_raises_for_malformed_file(tmp_path):
+    fname = tmp_path / "malformed.clust"
+    fname.write_bytes(b"\x00")
+
+    with ClusterFile(fname) as f, pytest.raises(RuntimeError):
+        f.read_frame()
+
+
+@pytest.mark.parametrize("use_roi", [False, True])
+@pytest.mark.parametrize("method", ["read_clusters", "default", "frames", "chunks"])
+@pytest.mark.parametrize(
+    "size",
+    [
+        pytest.param(3, id="partial-frame-number"),
+        pytest.param(4, id="missing-cluster-count"),
+        pytest.param(7, id="partial-cluster-count"),
+        pytest.param(48, id="missing-cluster-record"),
+        pytest.param(87, id="partial-cluster-record"),
+    ],
+)
+def test_chunk_reads_reject_incomplete_frames(tmp_path, use_roi, method, size):
+    fname = tmp_path / "incomplete.clust"
+    cluster = struct.pack("@HH9i", 5, 6, *range(9))
+    frame = struct.pack("@iI", 42, 2) + cluster * 2
+    fname.write_bytes(frame[:size])
+
+    with ClusterFile(fname) as reader:
+        if use_roi:
+            reader.set_roi(ROI(0, 10, 0, 10))
+        with pytest.raises(RuntimeError):
+            if method == "default":
+                next(reader)
+            elif method == "read_clusters":
+                reader.read_clusters(10)
+            else:
+                next(getattr(reader, method)())
+
+
 @pytest.mark.withdata
 def test_cluster_file(test_data_path): 
     """Test ClusterFile""" 

python/tests/test_ClusterVector.py

@@ -83,12 +83,15 @@ def test_make_a_hitmap_from_cluster_vector():
 
 def test_2x2_reduction(): 
     cv = ClusterVector((3,3))
+    cv.frame_number = -135
 
     cv.push_back(_aare.Cluster3x3i(5, 5, np.array([1, 1, 1, 2, 3, 1, 2, 2, 1], dtype=np.int32)))
     cv.push_back(_aare.Cluster3x3i(5, 5, np.array([2, 2, 1, 2, 3, 1, 1, 1, 1], dtype=np.int32)))
 
-    reduced_cv = np.array(_aare.reduce_to_2x2(cv), copy=False) 
+    reduced = _aare.reduce_to_2x2(cv)
+    reduced_cv = np.array(reduced, copy=False)
 
+    assert reduced.frame_number == cv.frame_number
     assert reduced_cv.size == 2
     assert reduced_cv[0]["x"] == 5
     assert reduced_cv[0]["y"] == 5
@@ -100,14 +103,17 @@ def test_2x2_reduction():
     
 def test_3x3_reduction(): 
     cv = _aare.ClusterVector_Cluster5x5d()
+    cv.frame_number = 246
     
     cv.push_back(_aare.Cluster5x5d(5,5,np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0, 2.0, 3.0,
                                    1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], dtype=np.double)))
     cv.push_back(_aare.Cluster5x5d(5,5,np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0, 2.0, 3.0,
                                    1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], dtype=np.double)))
     
-    reduced_cv = np.array(_aare.reduce_to_3x3(cv), copy=False)  
+    reduced = _aare.reduce_to_3x3(cv)
+    reduced_cv = np.array(reduced, copy=False)
 
+    assert reduced.frame_number == cv.frame_number
     assert reduced_cv.size == 2
     assert reduced_cv[0]["x"] == 5
     assert reduced_cv[0]["y"] == 5
@@ -130,3 +136,26 @@ def test_masking():
     assert cv_masked_array[0]["x"] == 1
     assert cv_masked_array[0]["y"] == 2
     assert (cv_masked_array[0]["data"] == np.ones((3,3),dtype=np.int32)).all()
+
+
+def test_masking_requires_c_contiguous_array():
+    cv = _aare.ClusterVector_Cluster3x3i()
+    cv.push_back(_aare.Cluster3x3i(1, 2, np.ones(9, dtype=np.int32)))
+    cv.push_back(_aare.Cluster3x3i(3, 4, np.ones(9, dtype=np.int32)))
+
+    mask = np.array([True, False, True, False], dtype=bool)[::2]
+    assert not mask.flags.c_contiguous
+
+    with pytest.raises(TypeError):
+        cv(mask)
+
+
+def test_masking_requires_one_dimension():
+    cv = _aare.ClusterVector_Cluster3x3i()
+    cv.push_back(_aare.Cluster3x3i(1, 2, np.ones(9, dtype=np.int32)))
+    cv.push_back(_aare.Cluster3x3i(3, 4, np.ones(9, dtype=np.int32)))
+
+    mask = np.array([[True, False]], dtype=bool)
+
+    with pytest.raises(ValueError):
+        cv(mask)

python/tests/test_CtbRawFile.py

@@ -0,0 +1,58 @@
+
+
+import pytest
+import numpy as np
+from aare import CtbRawFile
+
+
+@pytest.mark.withdata
+def test_basic_properites(test_data_path):
+    f = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    assert f.total_frames == 47
+    assert f.image_size_in_bytes == 640
+    assert f.master.analog_samples == 10
+    assert f.master.digital_samples is None
+    assert f.master.transceiver_samples is None
+
+
+@pytest.mark.withdata
+def test_seek_and_read(test_data_path):
+    f = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    f.seek(0)
+    header, raw_data = f.read_frame()
+    assert header['frameNumber'] == 48 #we know that the first frame in this file is 48
+
+    #Seek a multiple of frames in file
+    f.seek(9)
+    header, raw_data = f.read_frame()
+    assert header['frameNumber'] == 57 
+
+@pytest.mark.withdata
+def test_read_all_frames_direct_and_iterate(test_data_path):
+    """
+    We want to make sure that iterating and reading give the same result.
+    """
+    f0 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    f1 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+
+    for h0, d0 in f0:
+        h1, d1 = f1.read_frame()
+        assert h0 == h1
+        assert np.all(d0 == d1)
+
+@pytest.mark.withdata
+def test_read_all_frames_with_index_and_iterate(test_data_path):
+    """
+    We want to make sure that iterating and reading give the same result.
+    """
+    f0 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    f1 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+
+    frames_read = 0
+    for i,(h0, d0) in enumerate(f0):
+        h1, d1 = f1.read_frame(i)
+        assert h0 == h1
+        assert np.all(d0 == d1)
+        frames_read += 1
+
+    assert frames_read == f0.total_frames

python/tests/test_FastPedestal.py

@@ -143,3 +143,17 @@ def test_fast_pedestal_rejects_wrong_shape():
 
     with pytest.raises(RuntimeError, match="shape"):
         pedestal.add_init_frame(np.zeros((2, 2), dtype=np.uint16))
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+def test_fast_pedestal_std_uses_double_variance(dtype):
+    pedestal = FastPedestal(1, 1, n_samples=2, dtype=dtype)
+    pedestal.add_init_frame(np.array([[0]], dtype=np.uint16))
+    pedestal.add_init_frame(np.array([[1000]], dtype=np.uint16))
+
+    assert pedestal.std().dtype == dtype
+    np.testing.assert_array_equal(pedestal.std(), [[500]])
+
+    pedestal.push_ema(np.array([[1000]], dtype=np.uint16))
+    expected = np.array([[np.sqrt(187500.0)]], dtype=dtype)
+    np.testing.assert_array_equal(pedestal.std(), expected)

python/tests/test_Pedestal.py

@@ -1,7 +1,106 @@
 import numpy as np
 import pytest
 
-from aare import Pedestal_d, Pedestal_f
+from aare import Pedestal, Pedestal_d, Pedestal_f, Pedestal_i16
+
+
+@pytest.mark.parametrize(
+    ("dtype", "pedestal_type"),
+    [(np.float64, Pedestal_d), (np.float32, Pedestal_f), (np.int16, Pedestal_i16)],
+)
+def test_pedestal_factory(dtype, pedestal_type):
+    pedestal = Pedestal(rows=2, cols=3, n_samples=4, dtype=dtype)
+
+    assert isinstance(pedestal, pedestal_type)
+    assert pedestal.rows == 2
+    assert pedestal.cols == 3
+    assert pedestal.n_samples == 4
+
+
+def test_pedestal_factory_defaults_to_double():
+    pedestal = Pedestal(2, 3)
+
+    assert isinstance(pedestal, Pedestal_d)
+    assert pedestal.n_samples == 1000
+
+
+def test_pedestal_factory_rejects_unbound_dtype():
+    with pytest.raises(ValueError, match="Unsupported dtype for Pedestal"):
+        Pedestal(2, 3, dtype=np.int32)
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("name", ["rows", "cols", "n_samples"])
+def test_pedestal_factory_rejects_negative_parameters(dtype, name):
+    parameters = dict(rows=1, cols=1, n_samples=10)
+    parameters[name] = -1
+    with pytest.raises(TypeError):
+        Pedestal(**parameters, dtype=dtype)
+
+
+@pytest.mark.parametrize("pedestal_type", [Pedestal_d, Pedestal_f, Pedestal_i16])
+@pytest.mark.parametrize(
+    "args", [(-1, 1), (1, -1), (-1, 1, 10), (1, -1, 10), (1, 1, -1)]
+)
+def test_pedestal_constructor_rejects_negative_parameters(pedestal_type, args):
+    with pytest.raises(TypeError):
+        pedestal_type(*args)
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+def test_pedestal_std_stays_finite_after_settling(dtype):
+    pedestal = Pedestal(1, 1, n_samples=10, dtype=dtype)
+    pedestal.push(np.array([[16382]], dtype=np.uint16))
+    frame = np.array([[16383]], dtype=np.uint16)
+    for _ in range(201):
+        pedestal.push(frame)
+
+    noise = pedestal.std().item()
+    assert np.isfinite(noise)
+    assert 0 <= noise < 1e-3
+
+
+@pytest.mark.parametrize(
+    ("pedestal_type", "expected_dtype"),
+    [(Pedestal_d, np.float64), (Pedestal_f, np.float32), (Pedestal_i16, np.int16)],
+)
+def test_double_precision_moments(pedestal_type, expected_dtype):
+    pedestal = pedestal_type(1, 1, 2)
+    for value in [30000, 30003, 30002]:
+        pedestal.push(np.array([[value]], dtype=np.uint16))
+
+    assert pedestal.mean().dtype == expected_dtype
+    assert pedestal.std().dtype == expected_dtype
+    np.testing.assert_array_equal(
+        pedestal.mean(), np.array([[30001.75]], dtype=expected_dtype)
+    )
+    np.testing.assert_allclose(
+        pedestal.std(),
+        np.array([[np.sqrt(1.1875)]], dtype=expected_dtype),
+        rtol=1e-6,
+    )
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("method", ["push", "push_with_threshold"])
+@pytest.mark.parametrize("shape", [(2, 2), (3, 2)])
+def test_pedestal_rejects_mismatched_frame_shapes(dtype, method, shape):
+    pedestal = Pedestal(2, 3, dtype=dtype)
+    initial = np.full((2, 3), 7, dtype=np.uint16)
+    pedestal.push(initial)
+    frame = np.zeros(shape, dtype=np.uint16)
+    args = (
+        (np.full((2, 3), 10, dtype=dtype),)
+        if method == "push_with_threshold"
+        else ()
+    )
+
+    with pytest.raises(
+        RuntimeError, match="Frame shape does not match pedestal shape"
+    ):
+        getattr(pedestal, method)(frame, *args)
+
+    np.testing.assert_array_equal(pedestal.mean(), initial)
 
 
 @pytest.mark.parametrize(
@@ -38,3 +137,68 @@ def test_pedestal_exposes_mean_as_read_only_buffer():
     np.testing.assert_array_equal(mean, pedestal.view())
     assert np.shares_memory(mean, pedestal.view())
     assert not mean.flags.writeable
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("method", ["push", "push_with_threshold"])
+def test_pedestal_push_accepts_contiguous_frames(dtype, method):
+    pedestal = Pedestal(2, 3, dtype=dtype)
+    frame = np.arange(6, dtype=np.uint16).reshape(2, 3)
+    args = (
+        (np.full((2, 3), 10, dtype=dtype),)
+        if method == "push_with_threshold"
+        else ()
+    )
+
+    getattr(pedestal, method)(frame, *args)
+
+    np.testing.assert_array_equal(pedestal.mean(), frame)
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("input_name", ["push_frame", "threshold_frame", "threshold"])
+@pytest.mark.parametrize("layout", ["transpose", "slice", "reverse"])
+def test_pedestal_rejects_noncontiguous_inputs(dtype, input_name, layout):
+    pedestal = Pedestal(2, 3, dtype=dtype)
+    frame = np.zeros((2, 3), dtype=np.uint16)
+    threshold = np.full((2, 3), 10, dtype=dtype)
+    input_dtype = dtype if input_name == "threshold" else np.uint16
+    if layout == "transpose":
+        invalid = np.ones((3, 2), dtype=input_dtype).T
+    elif layout == "slice":
+        invalid = np.ones((2, 6), dtype=input_dtype)[:, ::2]
+    else:
+        invalid = np.ones((4, 6), dtype=input_dtype)[:2, :3][:, ::-1]
+    assert not invalid.flags.c_contiguous
+
+    with pytest.raises(TypeError):
+        if input_name == "push_frame":
+            pedestal.push(invalid)
+        elif input_name == "threshold_frame":
+            pedestal.push_with_threshold(invalid, threshold)
+        else:
+            pedestal.push_with_threshold(frame, invalid)
+
+    np.testing.assert_array_equal(pedestal.mean(), np.zeros((2, 3)))
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("input_name", ["push_frame", "threshold_frame", "threshold"])
+@pytest.mark.parametrize("shape", [(), (6,), (2, 3, 2)])
+def test_pedestal_rejects_inputs_with_wrong_ndim(dtype, input_name, shape):
+    pedestal = Pedestal(2, 3, dtype=dtype)
+    frame = np.ones((2, 3), dtype=np.uint16)
+    threshold = np.full((2, 3), 10, dtype=dtype)
+    input_dtype = dtype if input_name == "threshold" else np.uint16
+    invalid = np.ones(shape, dtype=input_dtype)
+    name = "Threshold" if input_name == "threshold" else "Frame"
+
+    with pytest.raises(ValueError, match=f"{name} must be 2-dimensional"):
+        if input_name == "push_frame":
+            pedestal.push(invalid)
+        elif input_name == "threshold_frame":
+            pedestal.push_with_threshold(invalid, threshold)
+        else:
+            pedestal.push_with_threshold(frame, invalid)
+
+    np.testing.assert_array_equal(pedestal.mean(), np.zeros((2, 3)))

python/tests/test_ROI.py

@@ -0,0 +1,21 @@
+import pytest 
+
+import numpy as np
+
+from aare import ROI 
+
+
+def test_roi_slice():
+    """test the ROI slice method"""
+    roi = ROI(0, 4, 2, 4)
+
+    array = np.array([[1, 2, 3, 4],
+                      [5, 6, 7, 8],
+                      [9, 10, 11, 12],
+                      [13, 14, 15, 16]])
+    
+    sliced_array = array[roi.slice()]
+
+    assert sliced_array.shape == (2, 4)
+    assert np.array_equal(sliced_array, np.array([[9, 10, 11, 12],
+                                                   [13, 14, 15, 16]]))

python/tests/test_RawFile.py

@@ -1,8 +1,299 @@
 # SPDX-License-Identifier: MPL-2.0
 import pytest
-from aare import RawFile
+import json
+from aare import File, RawFile, RawSubFile, DetectorType, ROI, UDPPortPosition
 import numpy as np
 
+
+@pytest.fixture
+def small_raw_file(tmp_path):
+    master_path = tmp_path / "run_master_0.json"
+    master_path.write_text(json.dumps({
+        "Version": 7.2,
+        "Detector Type": "Jungfrau",
+        "Timing Mode": "auto",
+        "Geometry": {"x": 1, "y": 1},
+        "Image Size in bytes": 12,
+        "Pixels": {"x": 3, "y": 2},
+        "Max Frames Per File": 1,
+        "Total Frames": 2,
+        "Frames in File": 2,
+        "Frame Padding": 1,
+        "Frame Discard Policy": "nodiscard",
+    }))
+    for index in range(2):
+        data_path = tmp_path / f"run_d0_f{index}_0.raw"
+        # A detector header followed by six uint16 pixels.
+        data_path.write_bytes(bytes(112) + np.arange(6, dtype=np.uint16).tobytes())
+    return master_path
+
+
+@pytest.mark.parametrize("reader_type", [RawFile, File])
+@pytest.mark.parametrize("legacy_master", [False, True])
+@pytest.mark.parametrize("padding, policy, supported", [
+    (0, "nodiscard", False),
+    (0, "discard", False),
+    (0, "discardpartial", True),
+    (1, "nodiscard", True),
+    (1, "discard", True),
+    (1, "discardpartial", True),
+])
+def test_raw_frame_policy(small_raw_file, reader_type, legacy_master,
+                          padding, policy, supported):
+    master_path = small_raw_file
+    if legacy_master:
+        master_path = master_path.with_suffix(".raw")
+        master_path.write_text(
+            "Version : 6.4\n"
+            "Detector Type : Jungfrau\n"
+            "Timing Mode : auto\n"
+            "Geometry : [1, 1]\n"
+            "Image Size : 12\n"
+            "Pixels : [3, 2]\n"
+            "Dynamic Range : 16\n"
+            "Max Frames Per File : 1\n"
+            "Total Frames : 2\n"
+            "Frames in File : 2\n"
+            f"Frame Padding : {padding}\n"
+            f"Frame Discard Policy : {policy}\n"
+        )
+    else:
+        metadata = json.loads(master_path.read_text())
+        metadata["Frame Padding"] = padding
+        metadata["Frame Discard Policy"] = policy
+        master_path.write_text(json.dumps(metadata))
+
+    if supported:
+        reader = reader_type(master_path)
+        assert reader.total_frames == 2
+        frame = reader.read_frame()
+        if reader_type is RawFile:
+            _, frame = frame
+        np.testing.assert_array_equal(frame, np.arange(6).reshape(2, 3))
+    else:
+        message = "requires frame padding or discardpartial"
+        with pytest.raises(RuntimeError, match=message) as error:
+            reader_type(master_path)
+        assert str(master_path) in str(error.value)
+
+        for index in range(2):
+            (master_path.parent / f"run_d0_f{index}_0.raw").unlink()
+        with pytest.raises(RuntimeError, match=message) as error:
+            reader_type(master_path)
+        assert str(master_path) in str(error.value)
+
+
+@pytest.mark.parametrize("method, args, kwargs", [
+    ("read_frame", (), {}),
+    ("read_n", (1,), {}),
+    ("read_roi", (0,), {}),
+    ("read_rois", (), {}),
+    ("read_n_with_roi", (1,), {"roi_index": 0}),
+    ("frame_number", (2,), {}),
+])
+def test_raw_read_error_context(small_raw_file, method, args, kwargs):
+    reader = RawFile(small_raw_file)
+    reader.seek(2)
+    with pytest.raises(RuntimeError) as error:
+        getattr(reader, method)(*args, **kwargs)
+    assert "frame index 2" in str(error.value)
+    assert str(small_raw_file) in str(error.value)
+    assert ".raw" not in str(error.value)
+    if method == "read_frame":
+        assert "file contains 2 frames (indices are zero-based)" in str(error.value)
+
+
+@pytest.mark.parametrize("frame_index", [2, 17])
+def test_generic_raw_read_error_context(small_raw_file, frame_index):
+    reader = File(small_raw_file)
+    with pytest.raises(RuntimeError) as error:
+        reader.read_frame(frame_index)
+    assert f"frame index {frame_index}" in str(error.value)
+    assert str(small_raw_file) in str(error.value)
+
+
+@pytest.mark.parametrize("raw_subfile", [False, True])
+def test_truncated_raw_read_error_context(small_raw_file, raw_subfile):
+    data_path = small_raw_file.parent / "run_d0_f1_0.raw"
+    if raw_subfile:
+        reader = RawSubFile(
+            small_raw_file.parent / "run_d0_f0_0.raw",
+            DetectorType.Jungfrau, 2, 3, 16,
+        )
+    else:
+        reader = RawFile(small_raw_file)
+    data_path.write_bytes(bytes(112))
+    reader.seek(1)
+    with pytest.raises(RuntimeError) as error:
+        reader.read_frame()
+    assert "frame index 1" in str(error.value)
+    assert str(data_path) in str(error.value)
+    assert "End of file" in str(error.value)
+    assert "master" not in str(error.value)
+    assert "RawSubFile.cpp" not in str(error.value)
+
+
+@pytest.mark.parametrize("method, args", [
+    ("read_frame", ()),
+    ("read_n", (1,)),
+    ("read_roi", (0,)),
+    ("frame_number", (1,)),
+])
+def test_short_raw_subfile_sets_read_bounds(small_raw_file, method, args):
+    (small_raw_file.parent / "run_d0_f1_0.raw").unlink()
+    reader = RawFile(small_raw_file)
+    reader.seek(1)
+    with pytest.raises(RuntimeError) as error:
+        getattr(reader, method)(*args)
+    assert reader.total_frames == 1
+    assert "frame index 1" in str(error.value)
+    assert str(small_raw_file) in str(error.value)
+    assert ".raw" not in str(error.value)
+
+
+@pytest.mark.parametrize("counts, master_count, warn", [
+    ([2], 2, False),
+    ([2, 2], 2, False),
+    ([2, 2], 0, True),
+    ([2, 2], 1, True),
+    ([2, 2], 5, True),
+    ([1, 2], 2, True),
+    ([2, 1], 2, True),
+    ([1, 2], 1, True),
+    ([0, 2], 2, True),
+    ([0, 0], 0, False),
+])
+def test_raw_subfile_frame_counts(small_raw_file, capfd, counts, master_count, warn):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata["Geometry"]["y"] = len(counts)
+    metadata["Frames in File"] = master_count
+    metadata["Total Frames"] = 1000
+    small_raw_file.write_text(json.dumps(metadata))
+    data = (small_raw_file.parent / "run_d0_f0_0.raw").read_bytes()
+    for module, count in enumerate(counts):
+        for index in range(2):
+            path = small_raw_file.parent / f"run_d{module}_f{index}_0.raw"
+            if index < count:
+                path.write_bytes(data)
+            elif index == 0:
+                path.write_bytes(b"")
+            else:
+                path.unlink(missing_ok=True)
+
+    reader = RawFile(small_raw_file)
+    warning = capfd.readouterr().err
+    if warn:
+        assert warning.count("WARNING") == 1
+        assert str(small_raw_file) in warning
+    else:
+        assert "WARNING" not in warning
+
+    assert reader.total_frames == min(counts)
+    assert len(reader) == min(counts)
+    assert reader.master.frames_in_file == master_count
+    if min(counts):
+        _, frames = reader.read()
+        assert frames.shape == (min(counts), 2 * len(counts), 3)
+        reader.seek(0)
+        _, frames = reader.read_n(10)
+        assert len(frames) == min(counts)
+        assert reader.frame_number(min(counts) - 1) == 0
+    with pytest.raises(RuntimeError):
+        reader.frame_number(min(counts))
+    reader.seek(min(counts))
+    with pytest.raises(RuntimeError):
+        reader.read_frame()
+    assert "WARNING" not in capfd.readouterr().err
+
+
+@pytest.mark.parametrize("short_module", [0, 1])
+def test_raw_frame_count_across_rois(small_raw_file, short_module):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata.update({
+        "Version": 8.0,
+        "Geometry": {"x": 1, "y": 2},
+        "Image Size": 12,
+        "Receiver Rois": [
+            {"xmin": 0, "xmax": 2, "ymin": 0, "ymax": 1},
+            {"xmin": 0, "xmax": 2, "ymin": 2, "ymax": 3},
+        ],
+    })
+    small_raw_file.write_text(json.dumps(metadata))
+    for index in range(2):
+        data = (small_raw_file.parent / f"run_d0_f{index}_0.raw").read_bytes()
+        (small_raw_file.parent / f"run_d1_f{index}_0.raw").write_bytes(data)
+    (small_raw_file.parent / f"run_d{short_module}_f1_0.raw").unlink()
+
+    reader = RawFile(small_raw_file)
+    assert reader.total_frames == 1
+    _, rois = reader.read_rois()
+    assert len(rois) == 2
+    assert all(roi.shape == (2, 3) for roi in rois)
+    for roi_index in range(2):
+        reader.seek(0)
+        _, frames = reader.read_n_with_roi(10, roi_index=roi_index)
+        assert frames.shape == (1, 2, 3)
+        with pytest.raises(RuntimeError):
+            reader.read_roi(roi_index)
+
+
+@pytest.mark.parametrize("disabled_port", [0, 1])
+def test_raw_frame_count_with_disabled_udp_port(small_raw_file, disabled_port):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata.update({
+        "Geometry": {"x": 1, "y": 2},
+        "Frames in File": 0,
+        "Number of UDP Interfaces": 2,
+        "UDP Ports Type": ["bottom", "top"],
+        "UDP Ports Disabled": [disabled_port],
+    })
+    small_raw_file.write_text(json.dumps(metadata))
+    data = (small_raw_file.parent / "run_d0_f0_0.raw").read_bytes()
+    for index in range(2):
+        (small_raw_file.parent / f"run_d0_f{index}_0.raw").unlink()
+        path = small_raw_file.parent / f"run_d{1 - disabled_port}_f{index}_0.raw"
+        path.write_bytes(data)
+
+    reader = RawFile(small_raw_file)
+    assert reader.total_frames == 2
+    assert reader.master.frames_in_file == 0
+    assert reader.frame_number(1) == 0
+    _, frames = reader.read()
+    assert frames.shape == (2, 2, 3)
+
+
+def test_raw_frame_count_without_selected_subfiles(small_raw_file):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata.update({"Version": 8.0, "Image Size": 12, "Receiver Rois": []})
+    small_raw_file.write_text(json.dumps(metadata))
+    with pytest.raises(RuntimeError, match="No raw subfiles selected"):
+        RawFile(small_raw_file)
+
+
+@pytest.mark.parametrize("lagging_module", [0, 1])
+def test_raw_synchronization_error_identifies_subfile(
+    small_raw_file, lagging_module
+):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata["Geometry"]["y"] = 2
+    small_raw_file.write_text(json.dumps(metadata))
+    for index in range(2):
+        data = (small_raw_file.parent / f"run_d0_f{index}_0.raw").read_bytes()
+        (small_raw_file.parent / f"run_d1_f{index}_0.raw").write_bytes(data)
+
+    ahead_path = small_raw_file.parent / f"run_d{1 - lagging_module}_f1_0.raw"
+    with ahead_path.open("r+b") as output:
+        output.write(np.array([1], dtype=np.uint64).tobytes())
+
+    reader = RawFile(small_raw_file)
+    reader.seek(1)
+    with pytest.raises(RuntimeError) as error:
+        reader.read_frame()
+    assert "frame index 1" in str(error.value)
+    assert "synchroniz" in str(error.value)
+    lagging_path = small_raw_file.parent / f"run_d{lagging_module}_f1_0.raw"
+    assert str(lagging_path) in str(error.value)
+
 @pytest.mark.withdata
 def test_read_rawfile_with_roi_spanning_over_one_module(test_data_path):
 
@@ -110,3 +401,41 @@ def test_read_rawfile_eiger_and_compare_to_numpy(test_data_path):
         header, image1 = f.read_frame()
 
     assert (image == image1).all()
+
+
+@pytest.mark.withdata 
+def test_read_eiger_udp_port_disabled(test_data_path):
+    with RawFile(test_data_path / "raw/eiger/one_udp_port_disabled_master_0.json") as f:
+        _, frame = f.read_rois() 
+
+        assert(len(frame) == 2)
+        assert frame[0].shape == (256, 512)
+        assert frame[1].shape == (256, 1024)
+        assert f.master.udp_port_types == [UDPPortPosition.LEFT, UDPPortPosition.RIGHT]
+        rois = f.master.rois
+        assert len(rois) == 2
+        assert rois[0] == ROI(512, 1024, 0, 256)
+        assert rois[1] == ROI(0, 1024, 256, 512)
+
+    with RawFile(test_data_path / "raw/eiger/quad_eiger_disabled_bottom_port_master_0.json") as f:
+        _, frame = f.read_frame() 
+
+        assert frame.shape == (256, 512)
+        assert(f.master.disabled_udp_ports == [1])
+        assert f.master.udp_port_types == [UDPPortPosition.TOP, UDPPortPosition.BOTTOM]
+        rois = f.master.rois
+        assert len(rois) == 1
+        assert rois[0] == ROI(0, 512, 256, 512)
+
+    with RawFile(test_data_path / "raw/eiger/2_modules_eiger_disabled_udp_port_master_0.json") as f:
+        _, frame = f.read_rois() 
+
+        assert(len(frame) == 2)
+        assert frame[0].shape == (512, 512)
+        assert frame[1].shape == (512, 512)
+        assert (f.master.disabled_udp_ports == [1, 3, 5, 7])
+        assert f.master.udp_port_types == [UDPPortPosition.LEFT, UDPPortPosition.RIGHT]
+        rois = f.master.rois
+        assert len(rois) == 2
+        assert rois[0] == ROI(0, 512, 0, 512)
+        assert rois[1] == ROI(1024, 1536, 0, 512)

python/tests/test_RawMasterFile.py

@@ -1,7 +1,9 @@
 
 
+import json
+
 import pytest
-from aare import RawMasterFile, ReadoutMode, DetectorType
+from aare import RawMasterFile, ReadoutMode, DetectorType, FrameDiscardPolicy
 
 
 @pytest.mark.withdata
@@ -11,4 +13,37 @@ def test_read_rawfile_quad_eiger_and_compare_to_numpy(test_data_path):
    
     f = RawMasterFile(file_name)
     assert(f.reading_mode == ReadoutMode.UNKNOWN)
-    assert(f.detector_type == DetectorType.Jungfrau)
+    assert(f.detector_type == DetectorType.Jungfrau)
+
+
+@pytest.mark.parametrize(
+    "policy, expected_policy",
+    [
+        ("nodiscard", FrameDiscardPolicy.NoDiscard),
+        ("discard", FrameDiscardPolicy.Discard),
+        ("discardpartial", FrameDiscardPolicy.DiscardPartial),
+    ],
+)
+def test_raw_master_file_context_manager(tmp_path, policy, expected_policy):
+    file_name = tmp_path / "run_master_0.json"
+    file_name.write_text(json.dumps({
+        "Version": 7.2,
+        "Detector Type": "Jungfrau",
+        "Timing Mode": "auto",
+        "Geometry": {"x": 1, "y": 1},
+        "Image Size in bytes": 12,
+        "Pixels": {"x": 3, "y": 2},
+        "Max Frames Per File": 1,
+        "Total Frames": 2,
+        "Frames in File": 2,
+        "Frame Padding": 1,
+        "Frame Discard Policy": policy,
+    }))
+
+    
+    with RawMasterFile(file_name) as context_file:
+        assert context_file.reading_mode == ReadoutMode.UNKNOWN
+        assert context_file.detector_type == DetectorType.Jungfrau
+        frame_discard_policy = context_file.frame_discard_policy
+        assert isinstance(frame_discard_policy, FrameDiscardPolicy)
+        assert frame_discard_policy == expected_policy

python/tests/test_make_view.py

@@ -0,0 +1,74 @@
+# SPDX-License-Identifier: MPL-2.0
+"""Probe make_view in np_helper.hpp through the bindings in _aare.testing"""
+import numpy as np
+import pytest
+
+from aare._aare import testing
+
+
+
+
+def test_make_view_roundtrip():
+    """Send an array to C++ and check that we get the same array back"""
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    out = testing.view_roundtrip_2d(arr)
+    np.testing.assert_array_equal(out, arr)
+
+    #Sanity check to make sure we got a copy an not just the same thing back
+    arr[3,2] = 31655
+    assert arr[3,2] != out[3,2]
+
+def test_make_view_const_roundtrip():
+    """Const works if the array is not writeable"""
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr.flags.writeable = False
+    out = testing.const_view_roundtrip_2d(arr)
+    np.testing.assert_array_equal(out, arr)
+
+def test_make_view_const_roundtrip_writable():
+    """Const works also if the array is  writeable"""
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    out = testing.const_view_roundtrip_2d(arr)
+    np.testing.assert_array_equal(out, arr)
+
+def test_make_view_fill():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    testing.view_fill_2d(arr, 7)
+    assert (arr == 7).all()
+
+def test_make_view_fill_rejects_nonwriteable():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr.flags.writeable = False
+    with pytest.raises(ValueError, match="not writeable"):
+        testing.view_fill_2d(arr, 7)
+
+def test_make_view_2d_rejects_1d():
+    arr = np.arange(30, dtype=np.uint16)
+    with pytest.raises(ValueError, match="Expected 2D array, got 1D"):
+        testing.view_roundtrip_2d(arr)
+
+def test_make_view_2d_rejects_3d():
+    arr = np.arange(30, dtype=np.uint16).reshape(3, 5, 2)
+    with pytest.raises(ValueError, match="Expected 2D array, got 3D"):
+        testing.view_roundtrip_2d(arr)
+
+def test_make_view_rejects_noncontiguous():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr = arr[::2, ::3]  # non-contiguous
+    assert not arr.flags.c_contiguous
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        testing.view_roundtrip_2d(arr)
+
+def test_make_view_rejects_noncontiguous_const():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr = arr[::2, ::3]  # non-contiguous
+    assert not arr.flags.c_contiguous
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        testing.const_view_roundtrip_2d(arr)
+
+def test_make_view_rejects_fortran_order():
+    arr = np.asfortranarray(np.arange(30, dtype=np.float64).reshape(5, 6))
+    assert not arr.flags.c_contiguous
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        testing.view_roundtrip_2d(arr)
+