Merge branch 'main' into dev/cluster-vector

CMakeLists.txt

@@ -55,6 +55,10 @@ option(
   "Install the python extension in the install tree under CMAKE_INSTALL_PREFIX/aare/"
   OFF)
 option(AARE_ASAN "Enable AddressSanitizer" OFF)
+option(
+  AARE_TUNE_LOCAL
+  "Optimize for the building machine's CPU (-march=native -mtune=native). Not portable, the resulting binaries might not run on other machines."
+  OFF)
 
 # Configure which of the dependencies to use FetchContent for
 option(AARE_FETCH_FMT "Use FetchContent to download fmt" ON)
@@ -296,6 +300,22 @@ else()
     target_compile_options(aare_compiler_flags INTERFACE -Werror)
   endif()
 
+  if(AARE_TUNE_LOCAL)
+    include(CheckCXXCompilerFlag)
+    check_cxx_compiler_flag("-march=native" AARE_HAS_MARCH_NATIVE)
+    check_cxx_compiler_flag("-mtune=native" AARE_HAS_MTUNE_NATIVE)
+    if(AARE_HAS_MARCH_NATIVE AND AARE_HAS_MTUNE_NATIVE)
+      message(STATUS "Tuning for the local CPU: -march=native -mtune=native")
+      target_compile_options(aare_compiler_flags INTERFACE -march=native
+                                                           -mtune=native)
+    else()
+      message(
+        WARNING
+          "AARE_TUNE_LOCAL requested but the compiler does not support -march=native/-mtune=native. Ignoring."
+      )
+    endif()
+  endif()
+
 endif() # GCC/Clang specific
 
 if(AARE_PYTHON_BINDINGS)
@@ -339,6 +359,7 @@ set(PUBLICHEADERS
     include/aare/Models.hpp
     include/aare/FileInterface.hpp
     include/aare/FilePtr.hpp
+    include/aare/FastPedestal.hpp
     include/aare/Frame.hpp
     include/aare/MultiThreadedFileReader.hpp
     include/aare/hist/PixelHistogram.hpp
@@ -449,6 +470,7 @@ if(AARE_TESTS)
       ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PixelHistogramImpl.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PixelHistogram.test.cpp
+      ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PedestalTrackingPixelHistogram.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/MultiThreadedFileReader.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp

RELEASE.md

@@ -2,8 +2,29 @@
 
 ## Next
 
+### New Features:
+
+- Added ``FastPedestal`` in C++ and Python for per-pixel running mean,
+  population variance, and standard deviation. It supports exponentially
+  weighted updates, initialization from files, direct subtraction from NumPy
+  arrays, and ``float64``, ``float32``, and ``int16`` output types.
+- Added the ``Pedestal_i16`` Python binding alongside
+  ``FastPedestal_d``, ``FastPedestal_f``, and ``FastPedestal_i16``.
+- ``PedestalTrackingPixelHistogram.fill_from_file()`` now uses parallel,
+  double-buffered file reading and accepts ``reader_threads`` and
+  ``reader_chunk_size`` tuning parameters.
+- Added the ``AARE_TUNE_LOCAL`` CMake option to build with ``-march=native``
+  and ``-mtune=native`` when supported. Binaries built with this option are
+  specific to the local CPU and may not be portable.
+
 ### API Changes:
 
+- ``ClusterFinder`` now uses ``FastPedestal``. It must receive 1000 pedestal
+  frames before cluster finding; ``find_clusters()`` raises
+  an error until initialization is complete. Added ``update_threshold()`` to
+  recompute the per-pixel detection thresholds.
+- Added the ``queue_depth`` constructor argument to ``ClusterFinderMT`` to
+  configure the number of preallocated frame buffers per worker thread.
 - Exposed ``ClusterVector.empty()`` in the Python API.
 - Added ClusterVector.estimate_n_clusters
 - Removed the lmfit dependency and the legacy ``fit_gaus``, ``fit_pol1``,
@@ -169,3 +190,7 @@ alice.mazzoleni@psi.ch \
 dhanya.thattil@psi.ch
 
 
+
+
+
+

benchmarks/CMakeLists.txt

@@ -18,8 +18,8 @@ FetchContent_MakeAvailable(benchmark)
 add_executable(benchmarks)
 
 target_sources(
-  benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp
-                     reduce_benchmark.cpp)
+  benchmarks PRIVATE ndarray_benchmark.cpp ndview_benchmark.cpp
+                     calculateeta_benchmark.cpp reduce_benchmark.cpp)
 
 # Link Google Benchmark and other necessary libraries
 target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core

benchmarks/ndview_benchmark.cpp

@@ -0,0 +1,17 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include <benchmark/benchmark.h>
+
+using aare::NDArray;
+using aare::NDView;
+
+static void BM_CreateNDView(benchmark::State &st) {
+    NDArray<int, 2> arr{{1024, 1024}, 0};
+    for (auto _ : st) {
+        // This code gets timed
+        auto res = arr.view();
+        benchmark::DoNotOptimize(res);
+    }
+}
+BENCHMARK(BM_CreateNDView);

docs/src/index.rst

@@ -30,6 +30,7 @@ AARE
     python/file/index
     python/experimental/index
     python/histogram/index
+    python/pedestal/index
     pyFit
     
 

docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst

@@ -15,6 +15,14 @@ Use ``push_pedestal_no_update()`` to seed the pedestal estimate, then
 asynchronous fills are drained by ``flush()``, and snapshot methods such as
 ``values()`` and ``pedestal_mean()`` return numpy arrays.
 
+``fill_from_file()`` uses parallel file-reader workers and a double-buffered
+pipeline. After the initial batch has been read, histogram processing of one
+batch overlaps reading of the next. ``reader_threads`` and
+``reader_chunk_size`` tune the I/O stage independently of the histogram worker
+count. Two fixed-capacity buffers are allocated once and reused by alternating
+their read and histogram roles. Their approximate memory use is ``2 *
+reader_threads * reader_chunk_size * rows * cols * sizeof(uint16)``.
+
 .. py:currentmodule:: aare
 
 .. autoclass:: PedestalTrackingPixelHistogram

docs/src/python/pedestal/index.rst

@@ -0,0 +1,9 @@
+Pedestal
+========
+
+.. toctree::
+   :caption: Pedestal
+   :maxdepth: 1
+
+   pyFastPedestal
+   pyPedestal

docs/src/python/pedestal/pyFastPedestal.rst

@@ -0,0 +1,83 @@
+FastPedestal
+============
+
+``FastPedestal`` calculates a running mean, variance and standard deviation for each pixel in a
+series of frames. The python binding only exposes ``uint16`` input but the underlying
+C++ class is templated. Initialize it with ``n_samples`` frames using
+``add_init_frame()``. Once ``ready`` is true, use ``push_ema()`` to update the exponential 
+moving average initialized by the mean and with smoothing factor 1/n_samples.
+
+.. warning::
+  
+  FastPedestal is not usable until you have added ``n_samples`` initial frames with ``add_init_frame(raw)``.
+  You can check the state with ``ready``.
+
+The public factory selects the bound C++ specialization from ``dtype``:
+
+* ``numpy.float64`` creates ``FastPedestal_d``
+* ``numpy.float32`` creates ``FastPedestal_f``
+* ``numpy.int16`` creates ``FastPedestal_i16``
+
+The internal calculations are done with double, but the cached mean and on demand var and std are returned in the specified type.
+
+Factory
+-------
+
+.. py:currentmodule:: aare
+
+.. autofunction:: FastPedestal
+
+Loading from a file
+-------------------
+
+``FastPedestal.from_file()`` initializes the pedestal from ``n_samples``
+frames after ``skip_first``, then applies steady-state updates for any frames
+remaining in the file. The input frames must contain ``uint16`` data; ``dtype``
+selects the output type of the pedestal statistics.
+
+.. autofunction:: aare.FastPedestal.from_file
+
+.. code-block:: python
+
+   pedestal = FastPedestal.from_file(
+       "frames.npy", n_samples=100, skip_first=10, dtype=np.float32
+   )
+
+Example
+-------
+
+.. code-block:: python
+
+   import numpy as np
+   from aare import FastPedestal
+
+   pedestal = FastPedestal(512, 1024, n_samples=100, dtype=np.float32)
+
+   # Initialize with n_samples frames
+   for frame in initialization_frames:
+       pedestal.add_init_frame(frame)
+
+   # Now we can push a frame for pedestal update
+   if pedestal.ready:
+       pedestal.push_ema(next_frame)
+       
+   # Mean and std are  also ready
+   mean = pedestal.mean()
+   noise = pedestal.std()
+   
+   # Direct pedestal subtraction is also supported
+   for frame in raw_data:
+       image = frame - pedestal
+
+Complete API
+------------
+
+The API below is for the ``float64`` specialization. All dtype variants share
+the same API.
+
+.. autoclass:: aare._aare.FastPedestal_d
+   :special-members: __init__
+   :members:
+   :undoc-members:
+   :show-inheritance:
+   :inherited-members:

docs/src/python/pedestal/pyPedestal.rst

@@ -0,0 +1,42 @@
+Pedestal
+========
+
+``Pedestal`` calculates a running mean and variance for each pixel in a series
+of ``uint16`` frames. ``push()`` updates the cached mean immediately. For
+faster batch initialization, use ``push_no_update()`` for each frame and call
+``update_mean()`` after the batch.
+
+Three specializations are available from :mod:`aare`:
+
+* ``Pedestal_d`` uses ``float64`` storage
+* ``Pedestal_f`` uses ``float32`` storage
+* ``Pedestal_i16`` uses ``int16`` storage
+
+Example
+-------
+
+.. code-block:: python
+
+   from aare import Pedestal_d
+
+   pedestal = Pedestal_d(512, 1024, 100)
+
+   for frame in initialization_frames:
+       pedestal.push_no_update(frame)
+
+   pedestal.update_mean()
+   mean = pedestal.mean()
+   noise = pedestal.std()
+
+Complete API
+------------
+
+The API below is for the ``float64`` specialization. All dtype variants share
+the same API.
+
+.. autoclass:: aare._aare.Pedestal_d
+   :special-members: __init__
+   :members:
+   :undoc-members:
+   :show-inheritance:
+   :inherited-members:

include/aare/Backoff.hpp

@@ -0,0 +1,48 @@
+// SPDX-License-Identifier: MPL-2.0
+#pragma once
+
+#include <chrono>
+#include <thread>
+
+namespace aare {
+
+/**
+ * @brief Hint the CPU that this is a spin-wait. Prefer this over yield() when
+ * the wait is expected to be short.
+ */
+inline void cpu_relax() noexcept {
+#if defined(__x86_64__) || defined(__i386__)
+    __builtin_ia32_pause();
+#elif defined(__aarch64__)
+    asm volatile("yield" ::: "memory");
+#else
+    std::this_thread::yield();
+#endif
+}
+
+/**
+ * @brief Escalating wait for producer/consumer idle and backpressure loops.
+ *
+ * Starts with a pause instruction (sub-microsecond, no syscall), then yield,
+ * then a short sleep. Call reset() whenever work arrives so a busy pipeline
+ * never leaves the pause tier.
+ */
+class Backoff {
+    int m_count{0};
+
+  public:
+    void reset() noexcept { m_count = 0; }
+
+    void pause() noexcept {
+        if (m_count < 64) {
+            cpu_relax();
+        } else if (m_count < 256) {
+            std::this_thread::yield();
+        } else {
+            std::this_thread::sleep_for(std::chrono::microseconds(50));
+        }
+        ++m_count;
+    }
+};
+
+} // namespace aare

include/aare/CircularFifo.hpp

@@ -1,11 +1,13 @@
 // SPDX-License-Identifier: MPL-2.0
 #pragma once
 
+#include <cassert>
 #include <chrono>
 #include <fmt/color.h>
 #include <fmt/format.h>
 #include <memory>
 #include <thread>
+#include <utility>
 
 #include "aare/ProducerConsumerQueue.hpp"
 
@@ -29,6 +31,23 @@ template <class ItemType> class CircularFifo {
         }
     }
 
+    /**
+     * @brief Construct a fifo and seed the free list using a factory.
+     * @param size number of items circulating in the fifo
+     * @param make_item callable invoked as make_item(i) for each slot index i
+     *
+     * Use this instead of CircularFifo(size) when the items need to be
+     * initialized, for example to hold a preallocated buffer or to carry
+     * their own slot index.
+     */
+    template <class F>
+    CircularFifo(uint32_t size, F make_item)
+        : fifo_size(size), free_slots(size + 1), filled_slots(size + 1) {
+        for (size_t i = 0; i < fifo_size; ++i) {
+            free_slots.write(make_item(i));
+        }
+    }
+
     bool next() {
         // TODO! avoid default constructing ItemType
         ItemType it;
@@ -47,6 +66,13 @@ template <class ItemType> class CircularFifo {
     auto numFreeSlots() const noexcept { return free_slots.sizeGuess(); }
     auto isFull() const noexcept { return filled_slots.isFull(); }
 
+    /**
+     * @brief True if there are no filled slots waiting to be consumed.
+     * @note Prefer this over numFilledSlots() == 0 since sizeGuess() may
+     * under-report when called from the producing thread.
+     */
+    auto isEmpty() const noexcept { return filled_slots.isEmpty(); }
+
     ItemType pop_free() {
         ItemType v;
         while (!free_slots.read(v))
@@ -75,8 +101,22 @@ template <class ItemType> class CircularFifo {
 
     ItemType *frontPtr() { return filled_slots.frontPtr(); }
 
-    // TODO! Add function to move item from filled to free to be used
-    // with the frontPtr function
+    /**
+     * @brief Return the front filled item to the free list. To be used
+     * together with frontPtr() once the item has been consumed in place.
+     * @warning The fifo must not be empty when calling this.
+     *
+     * The item is written to the free list before it is popped from the
+     * filled list, so it can never be dropped. The write cannot fail: both
+     * queues hold size + 1 slots while only size items circulate.
+     */
+    void recycle_front() {
+        ItemType *it = filled_slots.frontPtr();
+        assert(it != nullptr);
+        [[maybe_unused]] const bool ok = free_slots.write(std::move(*it));
+        assert(ok);
+        filled_slots.popFront();
+    }
 
     template <class... Args> void push_value(Args &&...recordArgs) {
         while (!filled_slots.write(std::forward<Args>(recordArgs)...))

include/aare/ClusterCollector.hpp

@@ -2,7 +2,11 @@
 #pragma once
 #include <atomic>
 #include <thread>
+#include <type_traits>
+#include <utility>
+#include <vector>
 
+#include "aare/Backoff.hpp"
 #include "aare/ClusterFinderMT.hpp"
 #include "aare/ClusterVector.hpp"
 #include "aare/ProducerConsumerQueue.hpp"
@@ -13,48 +17,64 @@ namespace aare {
 template <typename ClusterType,
           typename = std::enable_if_t<is_cluster_v<ClusterType>>>
 class ClusterCollector {
-    ProducerConsumerQueue<ClusterVector<ClusterType>> *m_source;
+    using SourceQueue = ProducerConsumerQueue<ClusterVector<ClusterType>>;
+
+    SourceQueue *m_source;
     std::atomic<bool> m_stop_requested{false};
     std::atomic<bool> m_stopped{true};
-    std::chrono::milliseconds m_default_wait{1};
     std::thread m_thread;
     std::vector<ClusterVector<ClusterType>> m_clusters;
 
     void process() {
-        m_stopped = false;
+        m_stopped.store(false, std::memory_order_release);
         fmt::print("ClusterCollector started\n");
-        while (!m_stop_requested || !m_source->isEmpty()) {
+        Backoff backoff;
+        while (true) {
             if (ClusterVector<ClusterType> *clusters = m_source->frontPtr();
                 clusters != nullptr) {
+                backoff.reset();
                 m_clusters.push_back(std::move(*clusters));
                 m_source->popFront();
-            } else {
-                std::this_thread::sleep_for(m_default_wait);
+                continue;
             }
+
+            if (m_stop_requested.load(std::memory_order_acquire))
+                break;
+
+            backoff.pause();
         }
         fmt::print("ClusterCollector stopped\n");
-        m_stopped = true;
+        m_stopped.store(true, std::memory_order_release);
     }
 
   public:
-    ClusterCollector(ClusterFinderMT<ClusterType, uint16_t, double> *source) {
-        m_source = source->sink();
+    explicit ClusterCollector(SourceQueue *source) : m_source(source) {
         m_thread =
             std::thread(&ClusterCollector::process,
                         this); // only one process does that so why isnt it
                                // automatically written to m_cluster in collect
                                // - instead of writing first to m_sink?
     }
+
+    template <typename Finder,
+              std::enable_if_t<
+                  std::is_same_v<decltype(std::declval<Finder &>().sink()),
+                                 SourceQueue *>,
+                  int> = 0>
+    explicit ClusterCollector(Finder *source)
+        : ClusterCollector(source->sink()) {}
+
     void stop() {
-        m_stop_requested = true;
-        m_thread.join();
+        m_stop_requested.store(true, std::memory_order_release);
+        if (m_thread.joinable())
+            m_thread.join();
     }
     std::vector<ClusterVector<ClusterType>> steal_clusters() {
-        if (!m_stopped) {
+        if (!m_stopped.load(std::memory_order_acquire)) {
             throw std::runtime_error("ClusterCollector is still running");
         }
         return std::move(m_clusters);
     }
 };
 
-} // namespace aare
+} // namespace aare

include/aare/ClusterFinder.hpp

@@ -3,6 +3,7 @@
 #include "aare/ClusterFile.hpp"
 #include "aare/ClusterVector.hpp"
 #include "aare/Dtype.hpp"
+#include "aare/FastPedestal.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
 #include "aare/Pedestal.hpp"
@@ -18,6 +19,13 @@ struct no_2x2_cluster {
         ClusterType::cluster_size_x > 2 && ClusterType::cluster_size_y > 2;
 };
 
+/**
+ * @brief Find fixed-size photon clusters using a per-pixel pedestal and noise
+ * threshold.
+ * @tparam ClusterType Output cluster type; both dimensions must exceed 2.
+ * @tparam FRAME_TYPE Input pixel type.
+ * @tparam PEDESTAL_TYPE Type used for pedestal and threshold calculations.
+ */
 template <typename ClusterType = Cluster<int32_t, 3, 3>,
           typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
           typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
@@ -26,51 +34,86 @@ class ClusterFinder {
     PEDESTAL_TYPE m_nSigma;
     const PEDESTAL_TYPE c2;
     const PEDESTAL_TYPE c3;
-    Pedestal<PEDESTAL_TYPE> m_pedestal;
+    FastPedestal<PEDESTAL_TYPE> m_pedestal;
     ClusterVector<ClusterType> m_clusters;
 
     static const uint8_t ClusterSizeX = ClusterType::cluster_size_x;
     static const uint8_t ClusterSizeY = ClusterType::cluster_size_y;
     using CT = typename ClusterType::value_type;
 
+    NDArray<PEDESTAL_TYPE, 2> m_threshold;
+    NDArray<PEDESTAL_TYPE, 2> m_pd_corrected_frame;
+
   public:
     /**
-     * @brief Construct a new ClusterFinder object
-     * @param image_size size of the image
-     * @param cluster_size size of the cluster (x, y)
-     * @param nSigma number of sigma above the pedestal to consider a photon
-     * @param capacity initial capacity of the cluster vector
-     *
+     * @brief Construct a cluster finder with an empty pedestal.
+     * @param image_size Image shape as (rows, columns).
+     * @param nSigma Per-pixel noise threshold multiplier.
+     * @param capacity Initial cluster-vector capacity.
      */
     ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
                   size_t capacity = 1000000)
         : m_image_size(image_size), m_nSigma(nSigma),
           c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)),
           c3(sqrt(ClusterSizeX * ClusterSizeY)),
-          m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {
+          m_pedestal(image_size[0], image_size[1]), m_clusters(capacity),
+          m_threshold({image_size[0], image_size[1]}, 0),
+          m_pd_corrected_frame({image_size[0], image_size[1]}, 0) {
         LOG(logDEBUG) << "ClusterFinder: "
                       << "image_size: " << image_size[0] << "x" << image_size[1]
                       << ", nSigma: " << nSigma << ", capacity: " << capacity;
     }
 
-    void set_nSigma(PEDESTAL_TYPE nSigma) { m_nSigma = nSigma; }
-
-    PEDESTAL_TYPE get_nSigma() const { return m_nSigma; }
-
-    void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
-        m_pedestal.push(frame);
+    /**
+     * @brief Set the noise multiplier used for threshold calculation and
+     * recompute the threshold.
+     * @param nSigma New per-pixel noise multiplier.
+     */
+    void set_nSigma(PEDESTAL_TYPE nSigma) {
+        m_nSigma = nSigma;
+        update_threshold();
     }
 
-    NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
-    NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
-    void clear_pedestal() { m_pedestal.clear(); }
+    /** @brief Return the current noise multiplier used for threshold
+     * calculation. */
+    PEDESTAL_TYPE get_nSigma() const { return m_nSigma; }
 
     /**
-     * @brief Move the clusters from the ClusterVector in the ClusterFinder to a
-     * new ClusterVector and return it.
-     * @param realloc_same_capacity if true the new ClusterVector will have the
-     * same capacity as the old one
+     * @brief Add a dark frame to the pedestal estimator.
      *
+     * The threshold is initialized automatically when the pedestal first
+     * becomes ready. Later frames update the ready pedestal.
+     * @param frame Dark frame matching the configured image shape.
+     * @throws std::runtime_error if the frame shape does not match.
+     */
+    void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
+        if (!m_pedestal.ready()) {
+            m_pedestal.add_init_frame(frame);
+            // Initialize the threshold when the pedestal becomes ready.
+            if (m_pedestal.ready()) {
+                update_threshold();
+            }
+        } else {
+            m_pedestal.push_ema(frame);
+        }
+    }
+
+    /** @brief Return a copy of the per-pixel pedestal mean. */
+    NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
+
+    /** @brief Return the per-pixel pedestal standard deviation (noise). */
+    NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
+
+    /** @brief Clear the pedestal and mark it as not ready. */
+    void clear_pedestal() { m_pedestal.clear(); }
+
+    /** @brief Recompute the threshold as noise multiplied by nSigma. */
+    void update_threshold() { m_threshold = m_pedestal.std() * m_nSigma; }
+
+    /**
+     * @brief Move out all accumulated clusters and reset the internal vector.
+     * @param realloc_same_capacity Preserve the previous capacity when true.
+     * @return The accumulated clusters and their frame metadata.
      */
     ClusterVector<ClusterType>
     steal_clusters(bool realloc_same_capacity = false) {
@@ -81,109 +124,179 @@ class ClusterFinder {
             m_clusters = ClusterVector<ClusterType>{};
         return tmp;
     }
+
+  private:
+    /**
+     * @brief Process a single pixel: scan its cluster window, decide whether it
+     * is a photon or a pedestal value, and store the cluster if needed.
+     * @tparam CheckBounds Skip out-of-image neighbours when true; assume the
+     * complete window is in bounds when false. Skipped cluster values remain 0.
+     */
+    template <bool CheckBounds>
+    void process_pixel(const NDView<FRAME_TYPE, 2> &frame, const int iy,
+                       const int ix) {
+
+        constexpr int dy = ClusterSizeY / 2;
+        constexpr int dx = ClusterSizeX / 2;
+        constexpr int has_center_pixel_x = ClusterSizeX % 2;
+        constexpr int has_center_pixel_y = ClusterSizeY % 2;
+
+        PEDESTAL_TYPE max = std::numeric_limits<PEDESTAL_TYPE>::lowest();
+        PEDESTAL_TYPE total = 0;
+
+        const int cols = static_cast<int>(frame.shape(1));
+        const int rows = static_cast<int>(frame.shape(0));
+        const auto center =
+            (static_cast<std::size_t>(iy) * static_cast<std::size_t>(cols)) +
+            static_cast<std::size_t>(ix);
+        const auto *corrected = m_pd_corrected_frame.data();
+        const PEDESTAL_TYPE threshold = m_threshold.data()[center];
+        const PEDESTAL_TYPE value = corrected[center];
+
+        if (value < -threshold)
+            return; // NEGATIVE_PEDESTAL, nothing to do for this pixel
+                    // TODO! No pedestal update???
+
+        for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
+            const int y = iy + ir;
+            if constexpr (CheckBounds) {
+                if (y < 0 || y >= rows)
+                    continue;
+            }
+
+            const auto *row = corrected + static_cast<std::size_t>(y) * cols;
+            for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
+                const int x = ix + ic;
+                if constexpr (CheckBounds) {
+                    if (x < 0 || x >= cols)
+                        continue;
+                }
+
+                const PEDESTAL_TYPE val = row[x];
+                total += val;
+                max = std::max(max, val);
+            }
+        }
+
+        if ((max > threshold)) {
+            if (value < max)
+                return; // Not max go to the next pixel, no pedestal update
+        } else if (total > c3 * threshold) {
+            // pass, store the cluster below
+        } else {
+            m_pedestal.push_ema_unchecked(center, frame.data()[center]);
+            return; // It was a pedestal value nothing to store
+        }
+
+        // Store cluster
+        if (value == max) {
+            ClusterType cluster{};
+            cluster.x = ix;
+            cluster.y = iy;
+
+            int i = 0;
+            for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
+                const int y = iy + ir;
+                for (int ic = -dx; ic < dx + has_center_pixel_x; ic++, i++) {
+                    const int x = ix + ic;
+                    if constexpr (CheckBounds) {
+                        if (x < 0 || x >= cols || y < 0 || y >= rows)
+                            continue;
+                    }
+
+                    const PEDESTAL_TYPE corrected_value =
+                        corrected[(static_cast<std::size_t>(y) * cols) + x];
+                    // If the cluster type is an integral type, and the
+                    // pedestal is a floating point type then we need to
+                    // round the value before storing it
+                    if constexpr (std::is_integral_v<CT> &&
+                                  std::is_floating_point_v<PEDESTAL_TYPE>) {
+                        cluster.data[i] =
+                            static_cast<CT>(std::lround(corrected_value));
+                    }
+                    // On the other hand if both are floating point or both
+                    // are integral then we can just static cast directly
+                    else {
+                        cluster.data[i] = static_cast<CT>(corrected_value);
+                    }
+                }
+            }
+
+            // Add the cluster to the output ClusterVector
+            m_clusters.push_back(cluster);
+        }
+    }
+
+  public:
+    /**
+     * @brief Find clusters in one frame and update eligible pedestal pixels.
+     * @param frame Input frame matching the configured image shape.
+     * @param frame_number Metadata assigned to the accumulated clusters.
+     * @pre frame has the same shape as image_size passed to the constructor.
+     * @throws std::runtime_error if the pedestal is not ready.
+     * @note Clusters accumulate until steal_clusters() is called.
+     */
     void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
+
         // // TODO! deal with even size clusters
         // // currently 3,3 -> +/- 1
         // //  4,4 -> +/- 2
-        int dy = ClusterSizeY / 2;
-        int dx = ClusterSizeX / 2;
-        int has_center_pixel_x =
-            ClusterSizeX %
-            2; // for even sized clusters there is no proper cluster center and
-               // even amount of pixels around the center
-        int has_center_pixel_y = ClusterSizeY % 2;
+        if (!m_pedestal.ready()) {
+            throw std::runtime_error(
+                "Pedestal is not ready, cannot find clusters");
+        }
+
+        constexpr int dy = ClusterSizeY / 2;
+        constexpr int dx = ClusterSizeX / 2;
+        constexpr int has_center_pixel_x = ClusterSizeX % 2;
+        constexpr int has_center_pixel_y = ClusterSizeY % 2;
+
+        // Largest neighbour offset below/right of the current pixel. Pixels
+        // further than this from an edge have their whole window in bounds.
+        constexpr int down = dy + has_center_pixel_y - 1;
+        constexpr int right = dx + has_center_pixel_x - 1;
 
         m_clusters.set_frame_number(frame_number);
-        for (int iy = 0; iy < frame.shape(0); iy++) {
-            for (int ix = 0; ix < frame.shape(1); ix++) {
 
-                PEDESTAL_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
-                PEDESTAL_TYPE total = 0;
+        const int rows = static_cast<int>(frame.shape(0));
+        const int cols = static_cast<int>(frame.shape(1));
 
-                // What can we short circuit here?
-                PEDESTAL_TYPE rms = m_pedestal.std(iy, ix);
-                PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix));
+        // TODO! See if we can get the same performace using the operator-
+        // m_pd_corrected_frame = frame - m_pedestal.view();
 
-                if (value < -m_nSigma * rms)
-                    continue; // NEGATIVE_PEDESTAL go to next pixel
-                              // TODO! No pedestal update???
+        // here we should be able to safely assume that the frame and corrected
+        // frame have the same size
+        auto n_pixels = frame.size();
+        auto pd = m_pedestal.view().data();
+        auto corrected = m_pd_corrected_frame.data();
+        auto frame_data = frame.data();
+        for (ssize_t i = 0; i < n_pixels; i++) {
+            corrected[i] = static_cast<PEDESTAL_TYPE>(frame_data[i]) - pd[i];
+        }
 
-                for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
-                    for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
-                        if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
-                            iy + ir >= 0 && iy + ir < frame.shape(0)) {
-                            PEDESTAL_TYPE val =
-                                frame(iy + ir, ix + ic) -
-                                m_pedestal.mean(iy + ir, ix + ic);
+        // Interior pixels can skip the per-neighbour bounds checks; pixels
+        // within dx/dy of an edge take the bounds-checked path. Iteration order
+        // (row-major, increasing ix) is preserved so results are identical.
+        const int ix_begin = dx;
+        const int ix_end = cols - right; // exclusive
 
-                            total += val;
-                            max = std::max(max, val);
-                        }
-                    }
-                }
+        for (int iy = 0; iy < rows; iy++) {
+            const bool interior_row = iy >= dy && iy < rows - down;
 
-                if ((max > m_nSigma * rms)) {
-                    if (value < max)
-                        continue; // Not max go to the next pixel
-                                  // but also no pedestal update
-                } else if (total > c3 * m_nSigma * rms) {
-                    // pass
-                } else {
-                    // m_pedestal.push(iy, ix, frame(iy, ix));   // Safe option
-                    m_pedestal.push_fast(
-                        iy, ix,
-                        frame(iy,
-                              ix)); // Assume we have reached n_samples in the
-                                    // pedestal, slight performance improvement
-                    continue;       // It was a pedestal value nothing to store
-                }
-
-                // Store cluster
-                if (value == max) {
-                    ClusterType cluster{};
-                    cluster.x = ix;
-                    cluster.y = iy;
-
-                    // Fill the cluster data since we have a photon to store
-                    // It's worth redoing the look since most of the time we
-                    // don't have a photon
-                    int i = 0;
-                    for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
-                        for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
-                            if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
-                                iy + ir >= 0 && iy + ir < frame.shape(0)) {
-
-                                // If the cluster type is an integral type, and
-                                // the pedestal is a floating point type then we
-                                // need to round the value before storing it
-                                if constexpr (std::is_integral_v<CT> &&
-                                              std::is_floating_point_v<
-                                                  PEDESTAL_TYPE>) {
-                                    auto tmp = std::lround(
-                                        frame(iy + ir, ix + ic) -
-                                        m_pedestal.mean(iy + ir, ix + ic));
-                                    cluster.data[i] = static_cast<CT>(tmp);
-                                }
-                                // On the other hand if both are floating point
-                                // or both are integral then we can just static
-                                // cast directly
-                                else {
-                                    auto tmp =
-                                        frame(iy + ir, ix + ic) -
-                                        m_pedestal.mean(iy + ir, ix + ic);
-                                    cluster.data[i] = static_cast<CT>(tmp);
-                                }
-                            }
-                            i++;
-                        }
-                    }
-
-                    // Add the cluster to the output ClusterVector
-                    m_clusters.push_back(cluster);
-                }
+            if (!interior_row || ix_begin >= ix_end) {
+                for (int ix = 0; ix < cols; ix++)
+                    process_pixel<true>(frame, iy, ix);
+                continue;
             }
+
+            for (int ix = 0; ix < ix_begin; ix++)
+                process_pixel<true>(frame, iy, ix);
+            for (int ix = ix_begin; ix < ix_end; ix++)
+                process_pixel<false>(frame, iy, ix);
+            for (int ix = ix_end; ix < cols; ix++)
+                process_pixel<true>(frame, iy, ix);
         }
     }
 };
 
-} // namespace aare
+} // namespace aare

include/aare/ClusterFinderMT.hpp

@@ -6,11 +6,13 @@
 #include <thread>
 #include <vector>
 
+#include "aare/Backoff.hpp"
+#include "aare/CircularFifo.hpp"
 #include "aare/ClusterFinder.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/ProducerConsumerQueue.hpp"
 #include "aare/logger.hpp"
-
+#include <ctime>
 namespace aare {
 
 enum class FrameType {
@@ -18,10 +20,40 @@ enum class FrameType {
     PEDESTAL,
 };
 
-struct FrameWrapper {
-    FrameType type;
-    uint64_t frame_number;
-    NDArray<uint16_t, 2> data;
+/**
+ * @brief Ticket identifying a frame buffer in a FramePool. Trivially
+ * copyable, the buffer itself never travels through the queues.
+ */
+struct FrameRef {
+    FrameType type{};
+    uint32_t slot{};
+    uint64_t frame_number{};
+};
+
+/**
+ * @brief Fixed set of frame buffers allocated once at construction.
+ *
+ * Buffers are addressed by slot index and are never moved or reallocated, so
+ * the queues only need to pass a FrameRef around. This keeps the per frame
+ * allocation out of the hot path entirely.
+ */
+class FramePool {
+    std::vector<NDArray<uint16_t, 2>> m_buffers;
+
+  public:
+    FramePool(size_t depth, Shape<2> shape) {
+        m_buffers.reserve(depth); // no reallocation, slots stay stable
+        for (size_t i = 0; i < depth; ++i) {
+            m_buffers.emplace_back(shape);
+        }
+    }
+
+    NDArray<uint16_t, 2> &operator[](uint32_t slot) { return m_buffers[slot]; }
+    const NDArray<uint16_t, 2> &operator[](uint32_t slot) const {
+        return m_buffers[slot];
+    }
+
+    size_t size() const { return m_buffers.size(); }
 };
 
 /**
@@ -42,17 +74,18 @@ class ClusterFinderMT {
     size_t m_current_thread{0};
     size_t m_n_threads{0};
     using Finder = ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>;
-    using InputQueue = ProducerConsumerQueue<FrameWrapper>;
+    using InputQueue = CircularFifo<FrameRef>;
     using OutputQueue = ProducerConsumerQueue<ClusterVector<ClusterType>>;
     std::vector<std::unique_ptr<InputQueue>> m_input_queues;
     std::vector<std::unique_ptr<OutputQueue>> m_output_queues;
+    std::vector<std::unique_ptr<FramePool>> m_frame_pools;
 
     OutputQueue m_sink{1000}; // All clusters go into this queue
 
     std::vector<std::unique_ptr<Finder>> m_cluster_finders;
     std::vector<std::thread> m_threads;
     std::thread m_collect_thread;
-    std::chrono::milliseconds m_default_wait{1};
+    std::chrono::microseconds m_default_wait{50};
 
   private:
     std::atomic<bool> m_stop_requested{false};
@@ -65,31 +98,35 @@ class ClusterFinderMT {
     void process(int thread_id) {
         auto cf = m_cluster_finders[thread_id].get();
         auto q = m_input_queues[thread_id].get();
-        bool realloc_same_capacity = true;
+        auto *pool = m_frame_pools[thread_id].get();
+        Backoff backoff;
 
         while (!m_stop_requested || !q->isEmpty()) {
-            if (FrameWrapper *frame = q->frontPtr(); frame != nullptr) {
+            if (FrameRef *ref = q->frontPtr(); ref != nullptr) {
+                backoff.reset();
+                auto view = (*pool)[ref->slot].view();
 
-                switch (frame->type) {
-                case FrameType::DATA:
-                    cf->find_clusters(frame->data.view(), frame->frame_number);
+                switch (ref->type) {
+                case FrameType::DATA: {
+                    cf->find_clusters(view, ref->frame_number);
+                    // Steal before the write so a failed write cannot drop the
+                    // clusters by re-stealing an empty vector on retry.
+                    auto clusters = cf->steal_clusters(true);
                     while (!m_output_queues[thread_id]->write(
-                        cf->steal_clusters(realloc_same_capacity))) {
-                        std::this_thread::sleep_for(m_default_wait);
+                        std::move(clusters))) {
+                        backoff.pause();
                     }
-
                     break;
-
+                }
                 case FrameType::PEDESTAL:
-                    m_cluster_finders[thread_id]->push_pedestal_frame(
-                        frame->data.view());
+                    m_cluster_finders[thread_id]->push_pedestal_frame(view);
                     break;
                 }
 
-                // frame is processed now discard it
-                m_input_queues[thread_id]->popFront();
+                // frame is processed, hand the buffer back to the free list
+                q->recycle_front();
             } else {
-                std::this_thread::sleep_for(m_default_wait);
+                backoff.pause();
             }
         }
     }
@@ -100,18 +137,24 @@ class ClusterFinderMT {
      */
     void collect() {
         bool empty = true;
+        Backoff backoff;
         while (!m_stop_requested || !empty || !m_processing_threads_stopped) {
-            empty = true;
+            bool moved_any = false;
             for (auto &queue : m_output_queues) {
-                if (!queue->isEmpty()) {
-
-                    while (!m_sink.write(std::move(*queue->frontPtr()))) {
-                        std::this_thread::sleep_for(m_default_wait);
+                while (auto *front = queue->frontPtr()) {
+                    while (!m_sink.write(std::move(*front))) {
+                        backoff.pause();
                     }
                     queue->popFront();
-                    empty = false;
+                    moved_any = true;
                 }
             }
+            empty = !moved_any;
+            if (moved_any) {
+                backoff.reset();
+            } else {
+                backoff.pause();
+            }
         }
     }
 
@@ -124,16 +167,22 @@ class ClusterFinderMT {
      * @param capacity initial capacity of the cluster vector. Should match
      * expected number of clusters in a frame per frame.
      * @param n_threads number of threads to use
+     * @param queue_depth number of frame buffers per thread. These are
+     * allocated once and recycled, so the total resident frame memory is
+     * n_threads * queue_depth * frame size. Keeping the in flight data below
+     * the L3 size keeps the per frame copy cheap.
      */
     ClusterFinderMT(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
-                    size_t capacity = 2000, size_t n_threads = 3)
+                    size_t capacity = 2000, size_t n_threads = 3,
+                    size_t queue_depth = 16)
         : m_n_threads(n_threads) {
 
         LOG(logDEBUG1) << "ClusterFinderMT: "
                        << "image_size: " << image_size[0] << "x"
                        << image_size[1] << ", nSigma: " << nSigma
                        << ", capacity: " << capacity
-                       << ", n_threads: " << n_threads;
+                       << ", n_threads: " << n_threads
+                       << ", queue_depth: " << queue_depth;
 
         for (size_t i = 0; i < n_threads; i++) {
             m_cluster_finders.push_back(
@@ -142,7 +191,13 @@ class ClusterFinderMT {
                     image_size, nSigma, capacity));
         }
         for (size_t i = 0; i < n_threads; i++) {
-            m_input_queues.emplace_back(std::make_unique<InputQueue>(200));
+            m_frame_pools.emplace_back(
+                std::make_unique<FramePool>(queue_depth, image_size));
+            m_input_queues.emplace_back(std::make_unique<InputQueue>(
+                static_cast<uint32_t>(queue_depth), [](size_t slot) {
+                    return FrameRef{FrameType::DATA,
+                                    static_cast<uint32_t>(slot), 0};
+                }));
             m_output_queues.emplace_back(std::make_unique<OutputQueue>(200));
         }
         // TODO! Should we start automatically?
@@ -212,13 +267,22 @@ class ClusterFinderMT {
      * expected to be dark. No photon finding is done. Just pedestal update.
      */
     void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
-        FrameWrapper fw{FrameType::PEDESTAL, 0,
-                        NDArray(frame)}; // TODO! copies the data!
+        for (size_t i = 0; i < m_n_threads; ++i) {
+            auto *q = m_input_queues[i].get();
+            Backoff backoff;
 
-        for (auto &queue : m_input_queues) {
-            while (!queue->write(fw)) {
-                std::this_thread::sleep_for(m_default_wait);
+            FrameRef ref;
+            while (!q->try_pop_free(ref)) {
+                backoff.pause();
             }
+
+            ref.type = FrameType::PEDESTAL;
+            ref.frame_number = 0;
+            (*m_frame_pools[i])[ref.slot].copy_from(frame);
+
+            // Cannot fail, the free list is what limits how many frames are
+            // in flight so there is always room in the filled list.
+            q->try_push_value(ref);
         }
     }
 
@@ -228,11 +292,29 @@ class ClusterFinderMT {
      * @note Spin locks with a default wait if the queue is full.
      */
     void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
-        FrameWrapper fw{FrameType::DATA, frame_number,
-                        NDArray(frame)}; // TODO! copies the data!
-        while (!m_input_queues[m_current_thread % m_n_threads]->write(fw)) {
-            std::this_thread::sleep_for(m_default_wait);
+        const size_t tid = m_current_thread % m_n_threads;
+        auto *q = m_input_queues[tid].get();
+        Backoff backoff;
+
+        FrameRef ref;
+        while (!q->try_pop_free(ref)) {
+            backoff.pause();
         }
+
+        ref.type = FrameType::DATA;
+        ref.frame_number = frame_number;
+
+        // DualTimer dt;
+        (*m_frame_pools[tid])[ref.slot].copy_from(frame);
+        // auto [wall_ns, cpu_ns] = dt.elapsed();
+        // std::cerr << "ClusterFinderMT: find_clusters: copied frame "
+        //           << frame_number << " took " << wall_ns/1000.0 << " wall_us
+        //           and "
+        //           << cpu_ns/1000.0 << " cpu_us" << std::endl;
+
+        // Cannot fail, the free list is what limits how many frames are in
+        // flight so there is always room in the filled list.
+        q->try_push_value(ref);
         m_current_thread++;
     }
 
@@ -245,6 +327,19 @@ class ClusterFinderMT {
         }
     }
 
+    /**
+     * @brief Recompute the threshold (nSigma * pedestal std) on all cluster
+     * finders. Requires the processing threads to be stopped.
+     */
+    void update_threshold() {
+        if (!m_processing_threads_stopped) {
+            throw std::runtime_error("ClusterFinderMT is still running");
+        }
+        for (auto &cf : m_cluster_finders) {
+            cf->update_threshold();
+        }
+    }
+
     /**
      * @brief Return the pedestal currently used by the cluster finder
      * @param thread_index index of the thread

include/aare/ClusterVector.hpp

@@ -59,6 +59,7 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
     }
 
     ClusterVector(ClusterVector &&other) noexcept = default;
+    ClusterVector &operator=(ClusterVector &&other) noexcept = default;
 
     /**
      * @brief Return a filtered copy selected by a one-dimensional Boolean mask.

include/aare/FastPedestal.hpp

@@ -0,0 +1,424 @@
+// SPDX-License-Identifier: MPL-2.0
+#pragma once
+#include "aare/File.hpp"
+#include "aare/Frame.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include <cstddef>
+#include <sys/types.h>
+
+namespace aare {
+
+/**
+ * @brief Maintain per-pixel mean, population variance, and standard deviation.
+ *
+ * Initialization accumulates exactly n_samples frames. Subsequent push_ema()
+ * update the exponential moving average initialized with the mean using
+ * a smoothing factor of (1/ n_samples). Internal moments are stored
+ * in double precision.
+ * @tparam PEDESTAL_TYPE Type returned for the mean, variance, and standard
+ * deviation.
+ */
+template <typename PEDESTAL_TYPE> class FastPedestal {
+
+    // Did we accumulate enough samples and updated the mean?
+    bool m_ready = false;
+
+    uint32_t m_rows;
+    uint32_t m_cols;
+
+    uint32_t m_samples;
+    double m_inv_samples;       // precompute 1/m_samples for faster division
+    uint32_t m_cur_samples = 0; // number of samples accumulated so far
+
+    // For cache locality we want to keep sum and sum2 close. Should
+    // improve performance for random access.
+    struct Entry {
+        double sum;
+        double sum2;
+    };
+    NDArray<Entry, 2> m_sum;
+
+    // Cache mean since it is used over and over in the ClusterFinder
+    // This optimization is related to the access pattern of the ClusterFinder
+    // Relies on having more reads than pushes to the pedestal
+    // But also makes sense when subtracting the pedestal from the frame
+    NDArray<PEDESTAL_TYPE, 2> m_mean;
+
+    // Helper function to convert row and column indices to a flat index
+    // used to provide both row column and flat index access to the pedestal
+    size_t rc_to_index(uint32_t row, uint32_t col) const {
+        return (static_cast<std::size_t>(row) * m_cols) + col;
+    }
+
+  public:
+    /**
+     * @brief Construct an empty pedestal that becomes ready after n_samples
+     * initialization frames.
+     * @param rows Number of image rows.
+     * @param cols Number of image columns.
+     * @param n_samples Number of initialization frames and reciprocal of the
+     * weight assigned to each subsequent value.
+     * @throws std::runtime_error if rows, cols, or n_samples is zero.
+     */
+    FastPedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
+        : m_rows(rows), m_cols(cols), m_samples(n_samples),
+          m_inv_samples(1.0 / n_samples), m_sum({rows, cols}, Entry{0, 0}),
+          m_mean({rows, cols}, PEDESTAL_TYPE(0)) {
+        if (!(rows > 0 && cols > 0 && n_samples > 0)) {
+            throw std::runtime_error(
+                fmt::format("Invalid parameters for FastPedestal: rows={}, "
+                            "cols={}, n_samples={} need to be positive",
+                            rows, cols, n_samples));
+        }
+    }
+
+    ~FastPedestal() = default;
+
+    /**
+     * @brief Return a non-owning view of the cached mean.
+     * @throws std::runtime_error if ready() is false.
+     * @note The caller must treat the data as read-only and must not retain the
+     * view after this object is destroyed, moved, or assigned.
+     */
+    NDView<const PEDESTAL_TYPE, 2> view() const {
+        if (!ready()) {
+            throw std::runtime_error(
+                "Pedestal is not ready, cannot return view");
+        }
+        return m_mean.view();
+    }
+
+    /**
+     * @brief Return a copy of the cached mean.
+     * @throws std::runtime_error if ready() is false.
+     */
+    NDArray<PEDESTAL_TYPE, 2> mean() {
+        if (!ready()) {
+            throw std::runtime_error(
+                "Pedestal is not ready, cannot return mean");
+        }
+        return m_mean;
+    }
+
+    /**
+     * @brief Return the cached mean at (row, col).
+     * @throws std::runtime_error if ready() is false or either index is out of
+     * range.
+     */
+    PEDESTAL_TYPE mean(uint32_t row, uint32_t col) const {
+        if (!ready()) {
+            throw std::runtime_error(
+                "Pedestal is not ready, cannot return mean");
+        }
+        if (row >= m_rows || col >= m_cols) {
+            throw std::runtime_error(
+                fmt::format("Invalid indices for FastPedestal mean: row={}, "
+                            "col={} must be in [0, {}), [0, {})",
+                            row, col, m_rows, m_cols));
+        }
+        return m_mean(row, col);
+    }
+
+    /**
+     * @brief Return the cached mean at a flat row-major index.
+     * @pre ready() is true, the cache is current, and index is valid; the index
+     * is not checked.
+     */
+    PEDESTAL_TYPE mean_unchecked(ssize_t index) const { return m_mean[index]; }
+
+    /**
+     * @brief Calculate and return the population variance of every pixel.
+     * @throws std::runtime_error if ready() is false.
+     * @note The result is normalized by n_samples.
+     */
+    NDArray<PEDESTAL_TYPE, 2> variance() {
+        if (!ready()) {
+            throw std::runtime_error(
+                "Pedestal is not ready, cannot return variance");
+        }
+        NDArray<PEDESTAL_TYPE, 2> res({m_rows, m_cols});
+        for (ssize_t i = 0; i < m_sum.size(); ++i) {
+            res[i] = variance_unchecked(i);
+        }
+        return res;
+    }
+
+    /**
+     * @brief Calculate the population variance at (row, col).
+     * @throws std::runtime_error if ready() is false or either index is out of
+     * range.
+     */
+    PEDESTAL_TYPE variance(const uint32_t row, const uint32_t col) const {
+        if (!ready()) {
+            throw std::runtime_error(
+                "Pedestal is not ready, cannot return variance");
+        }
+        if (row >= m_rows || col >= m_cols) {
+            throw std::runtime_error(fmt::format(
+                "Invalid indices for FastPedestal variance: row={}, "
+                "col={} must be in [0, {}), [0, {})",
+                row, col, m_rows, m_cols));
+        }
+        return variance_unchecked(rc_to_index(row, col));
+    }
+
+    /**
+     * @brief Calculate the population variance at a flat row-major index.
+     * @pre ready() is true and index is valid; the index is not checked.
+     */
+    PEDESTAL_TYPE variance_unchecked(ssize_t index) const {
+        const auto &entry = m_sum[index];
+        const auto m = entry.sum * m_inv_samples;
+        return std::fma(-m, m, entry.sum2 * m_inv_samples);
+    }
+
+    /**
+     * @brief Calculate and return the population standard deviation of every
+     * pixel.
+     * @throws std::runtime_error if ready() is false.
+     */
+    NDArray<PEDESTAL_TYPE, 2> std() {
+        if (!ready()) {
+            throw std::runtime_error(
+                "Pedestal is not ready, cannot return std");
+        }
+        NDArray<PEDESTAL_TYPE, 2> res({m_rows, m_cols});
+        for (ssize_t i = 0; i < m_sum.size(); ++i) {
+            res[i] = std_unchecked(i);
+        }
+        return res;
+    }
+
+    /**
+     * @brief Calculate the population standard deviation at (row, col).
+     * @throws std::runtime_error if ready() is false or either index is out of
+     * range.
+     */
+    PEDESTAL_TYPE std(const uint32_t row, const uint32_t col) const {
+        if (!ready()) {
+            throw std::runtime_error(
+                "Pedestal is not ready, cannot return std");
+        }
+        if (row >= m_rows || col >= m_cols) {
+            throw std::runtime_error(
+                fmt::format("Invalid indices for FastPedestal std: row={}, "
+                            "col={} must be in [0, {}), [0, {})",
+                            row, col, m_rows, m_cols));
+        }
+        return std::sqrt(variance(row, col));
+    }
+
+    /**
+     * @brief Calculate the population standard deviation at a flat row-major
+     * index.
+     * @pre ready() is true and index is valid; the index is not checked.
+     */
+    PEDESTAL_TYPE std_unchecked(ssize_t index) const {
+        return std::sqrt(variance_unchecked(index));
+    }
+
+    /**
+     * @brief Return whether initialization is complete (cur_samples() equals
+     * n_samples()).
+     */
+    bool ready() const { return m_ready; }
+
+    /**
+     * @brief Return the stored number of accumulated initialization frames.
+     * @note The value is in [0, n_samples] and does not change during
+     * steady-state pushes.
+     */
+    uint32_t cur_samples() const { return m_cur_samples; }
+
+    /**
+     * @brief Zero the moments and cached mean, and mark the pedestal not ready.
+     */
+    void clear() {
+        m_sum = Entry{0., 0.};
+        m_mean = PEDESTAL_TYPE(0.);
+        m_cur_samples = 0;
+        m_ready = false;
+    }
+
+    /**
+     * @brief Update every pixel using the steady-state exponential estimator,
+     * giving the new value weight 1 / n_samples.
+     * @param frame Frame whose shape must exactly match the pedestal.
+     * @throws std::runtime_error if the shape differs or ready() is false.
+     */
+    template <typename T> void push_ema(NDView<T, 2> frame) {
+        if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
+            throw std::runtime_error(
+                "Frame shape does not match pedestal shape");
+        }
+
+        if (!ready()) {
+            throw std::runtime_error("Pedestal is not ready, cannot push");
+        }
+
+        const auto size = static_cast<std::size_t>(m_rows) * m_cols;
+        const auto *data = frame.data();
+        for (std::size_t index = 0; index < size; ++index) {
+            push_ema_unchecked(index, data[index]);
+        }
+    }
+
+    /**
+     * @brief Update every pixel from a Frame using the steady-state estimator.
+     * @tparam T Actual pixel type stored in frame; this is not runtime-checked.
+     * @param frame Frame whose shape must exactly match the pedestal.
+     * @throws std::runtime_error if the shape differs or ready() is false.
+     */
+    template <typename T> void push_ema(Frame &frame) {
+        push_ema<T>(frame.view<T>());
+    }
+
+    /**
+     * @brief Update the exponential moving average with smoothing factor
+     * 1/n_samples
+     * @param row Pixel row.
+     * @param col Pixel column.
+     * @param val New pixel value.
+     * @pre row and col are valid; indices are not checked.
+     * @throws std::runtime_error if ready() is false.
+     */
+    template <typename T>
+    void push_ema(const uint32_t row, const uint32_t col, const T val) {
+        if (!ready()) {
+            throw std::runtime_error("Pedestal is not ready, cannot push");
+        }
+
+        push_ema_unchecked(rc_to_index(row, col), val);
+    }
+
+    /**
+     * @brief Update one pixel and its cached mean without runtime checks in
+     * release builds.
+     * @param index Flat row-major pixel index.
+     * @param value New pixel value, with weight 1 / n_samples.
+     * @pre ready() is true and index is a valid flat row-major index. These
+     * preconditions are asserted only in debug builds.
+     */
+    template <typename T>
+    void push_ema_unchecked(const std::size_t index, const T value) noexcept {
+        assert(m_ready);
+        assert(index < static_cast<std::size_t>(m_sum.size()));
+
+        const auto val = static_cast<double>(value);
+        auto &entry = m_sum[index];
+        entry.sum += val - entry.sum * m_inv_samples;
+        entry.sum2 += val * val - entry.sum2 * m_inv_samples;
+        m_mean[index] = static_cast<PEDESTAL_TYPE>(entry.sum * m_inv_samples);
+    }
+
+    /**
+     * @brief Accumulate one initialization frame.
+     * @param frame Frame whose shape must exactly match the pedestal.
+     * @throws std::runtime_error if the shape differs or n_samples frames have
+     * already been accumulated.
+     * @note The statistics can be accessed and ready() becomes true only after
+     * the n_samples frames have been added.
+     */
+    template <typename T> void add_init_frame(NDView<T, 2> frame) {
+        if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
+            throw std::runtime_error(
+                "Frame shape does not match pedestal shape");
+        }
+
+        // if the pedestal is already initialized we cannot add more frames
+        if (ready()) {
+            throw std::runtime_error("Pedestal initialization is already done");
+        }
+
+        for (ssize_t i = 0; i < m_sum.size(); ++i) {
+            const auto val = static_cast<double>(frame[i]);
+            auto &entry = m_sum[i];
+            entry.sum += val;
+            entry.sum2 += val * val;
+        }
+        m_cur_samples += 1;
+
+        if (m_cur_samples == m_samples) {
+            update_mean();
+            m_ready = true;
+        }
+    }
+
+    /**
+     * @brief Initialize from n_samples frames after skip_first, then apply all
+     * remaining file frames as steady-state updates.
+     * @tparam T Pixel representation stored in the file; this is not checked.
+     * @param filename Input image file. Its dimensions define the pedestal
+     * shape.
+     * @param n_samples Number of initialization frames.
+     * @param skip_first Number of leading frames to ignore.
+     * @throws std::runtime_error if fewer than n_samples frames remain after
+     * skip_first, or if any constructor argument is invalid.
+     */
+    template <typename T>
+    static FastPedestal from_file(const std::filesystem::path &filename,
+                                  uint32_t n_samples = 1000,
+                                  uint32_t skip_first = 0) {
+        File f(filename);
+
+        const auto total_frames = f.total_frames();
+        const auto first_frame = static_cast<size_t>(skip_first);
+        const auto initialization_frames = static_cast<size_t>(n_samples);
+        if (first_frame > total_frames ||
+            initialization_frames > total_frames - first_frame) {
+            throw std::runtime_error(
+                "File has less frames than the number of samples needed to "
+                "initialize the pedestal");
+        }
+
+        if (skip_first > 0) {
+            f.seek(static_cast<size_t>(skip_first));
+        }
+        const auto rows = static_cast<uint32_t>(f.rows());
+        const auto cols = static_cast<uint32_t>(f.cols());
+        FastPedestal pedestal(rows, cols, n_samples);
+        NDArray<T, 2> frame({rows, cols});
+
+        auto frame_index = first_frame;
+        const auto initialization_end = first_frame + initialization_frames;
+        while (frame_index < initialization_end) {
+            f.read_into(frame.buffer());
+            pedestal.template add_init_frame<T>(frame.view());
+            frame_index++;
+        }
+
+        // read the rest of the file
+        while (frame_index < total_frames) {
+            f.read_into(frame.buffer());
+            pedestal.template push_ema<T>(frame.view());
+            frame_index++;
+        }
+        return pedestal;
+    }
+
+    /** @brief Return the number of image rows. */
+    uint32_t rows() const { return m_rows; }
+
+    /** @brief Return the number of image columns. */
+    uint32_t cols() const { return m_cols; }
+
+    /**
+     * @brief Return the initialization frame count and steady-state
+     * update-weight denominator.
+     */
+    uint32_t n_samples() const { return m_samples; }
+
+  private:
+    /**
+     * @brief Write the cached mean after the final add_init_frame. All other
+     * (non initialization) pushes update the cached mean immediately.
+     */
+    void update_mean() {
+        for (ssize_t i = 0; i < m_sum.size(); i++) {
+            auto &entry = m_sum[i];
+            m_mean[i] = static_cast<PEDESTAL_TYPE>(entry.sum * m_inv_samples);
+        }
+    }
+};
+} // namespace aare

include/aare/NDArray.hpp

@@ -146,6 +146,23 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      */
     ~NDArray() { delete[] data_; }
 
+    /**
+     * @brief Copy data from a view of matching shape into this array without
+     * reallocating.
+     * @param v view to copy from, must have the same shape as this array
+     * @throws std::runtime_error if the shapes differ
+     *
+     * Use this instead of assigning a new NDArray when the buffer needs to be
+     * kept, for example when the array is part of a preallocated pool.
+     */
+    void copy_from(const NDView<T, Ndim> v) {
+        if (v.shape() != shape_) {
+            throw std::runtime_error(LOCATION +
+                                     "Shape mismatch in NDArray::copy_from");
+        }
+        std::copy(v.begin(), v.end(), begin());
+    }
+
     ///////////////////////////////////////////////////////////////////////////////
     // Iterators and indexing
     //

include/aare/Pedestal.hpp

@@ -29,17 +29,22 @@ template <typename SUM_TYPE = double> class Pedestal {
     // Relies on having more reads than pushes to the pedestal
     NDArray<SUM_TYPE, 2> m_mean;
 
+    // Cache std. Only refreshed via update_std() to keep push() cheap.
+    NDArray<SUM_TYPE, 2> m_std;
+
   public:
     Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
         : m_rows(rows), m_cols(cols), m_samples(n_samples),
           m_cur_samples(NDArray<uint32_t, 2>({rows, cols}, 0)),
           m_sum(NDArray<SUM_TYPE, 2>({rows, cols})),
           m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})),
-          m_mean(NDArray<SUM_TYPE, 2>({rows, cols})) {
+          m_mean(NDArray<SUM_TYPE, 2>({rows, cols})),
+          m_std(NDArray<SUM_TYPE, 2>({rows, cols})) {
         assert(rows > 0 && cols > 0 && n_samples > 0);
         m_sum = 0;
         m_sum2 = 0;
         m_mean = 0;
+        m_std = 0;
     }
     ~Pedestal() = default;
 
@@ -51,6 +56,12 @@ template <typename SUM_TYPE = double> class Pedestal {
         return m_mean(row, col);
     }
 
+    NDArray<SUM_TYPE, 2> cached_std() { return m_std; }
+
+    SUM_TYPE cached_std(const uint32_t row, const uint32_t col) const {
+        return m_std(row, col);
+    }
+
     SUM_TYPE std(const uint32_t row, const uint32_t col) const {
         return std::sqrt(variance(row, col));
     }
@@ -87,6 +98,7 @@ template <typename SUM_TYPE = double> class Pedestal {
         m_sum2 = 0;
         m_cur_samples = 0;
         m_mean = 0;
+        m_std = 0;
     }
 
     void clear(const uint32_t row, const uint32_t col) {
@@ -94,6 +106,7 @@ template <typename SUM_TYPE = double> class Pedestal {
         m_sum2(row, col) = 0;
         m_cur_samples(row, col) = 0;
         m_mean(row, col) = 0;
+        m_std(row, col) = 0;
     }
 
     template <typename T> void push(NDView<T, 2> frame) {
@@ -211,6 +224,17 @@ template <typename SUM_TYPE = double> class Pedestal {
      */
     void update_mean() { m_mean = m_sum / m_cur_samples; }
 
+    /**
+     * @brief Refresh the cached std for all pixels from the current sums.
+     * Kept separate from push() so pushes stay cheap; call before reading
+     * cached_std() (analogous to update_mean()).
+     */
+    void update_std() {
+        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+            m_std(i / m_cols, i % m_cols) = std(i / m_cols, i % m_cols);
+        }
+    }
+
     template <typename T>
     void push_fast(const uint32_t row, const uint32_t col, const T val_) {
         // Assume we reached the steady state where all pixels have

include/aare/hist/PedestalTrackingPixelHistogram.hpp

@@ -1,4 +1,5 @@
 #pragma once
+#include "aare/FastPedestal.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
 #include "aare/Pedestal.hpp"
@@ -44,7 +45,7 @@ class PedestalTrackingPixelHistogram {
     // worker using the LOCAL row index (i.e. 0..row_count(t)-1), NOT the
     // global row index. Owned exclusively by worker `t` during a
     // dispatched fan-out.
-    std::vector<Pedestal<AxisType>> partial_pedestals_;
+    std::vector<FastPedestal<AxisType>> partial_pedestals_;
     std::vector<NDArray<AxisType, 2>> partial_std_; // cached for pedestal
                                                     // tracking
 
@@ -65,6 +66,11 @@ class PedestalTrackingPixelHistogram {
     // Always the outermost lock; work_mutex_ is taken briefly inside it.
     mutable std::mutex fill_mutex_;
 
+    // Serialises producers of async_queue_ (which is SPSC) and gives
+    // fill_from_file exclusive ownership of the ingestion path while its
+    // direct, double-buffered dispatch is active.
+    std::mutex ingestion_mutex_;
+
     // Async producer/consumer pipeline. SPSC queue feeds the coordinator
     // thread, which batches queued images before dispatching them.
     std::unique_ptr<AsyncQueue> async_queue_;
@@ -107,8 +113,23 @@ class PedestalTrackingPixelHistogram {
 
     void fill_async(NDArray<FrameType, 2> &&image);
 
+    /**
+     * @brief Fill from ordered, parallel-read batches.
+     *
+     * After the first batch, reading batch N+1 overlaps histogram processing
+     * for batch N. Two fixed-capacity read buffers are allocated once and
+     * reused for the duration of the call.
+     *
+     * @param fname input file accepted by File
+     * @param max_frames maximum frames to process, or all frames for -1
+     * @param verbose print periodic progress information
+     * @param reader_threads number of MultiThreadedFileReader workers
+     * @param reader_chunk_size frames claimed per reader worker and batch
+     */
     void fill_from_file(const std::filesystem::path &fname,
-                        ssize_t max_frames = -1, bool verbose = false);
+                        ssize_t max_frames = -1, bool verbose = false,
+                        std::size_t reader_threads = 2,
+                        std::size_t reader_chunk_size = 4);
 
     void process_pedestal_file(const std::filesystem::path &fname,
                                ssize_t max_frames = -1, bool verbose = false);

include/aare/hist/PixelHistogramImpl.hpp

@@ -10,6 +10,7 @@ silently dropped.
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
 
+#include <algorithm>
 #include <cstddef>
 #include <limits>
 #include <stdexcept>
@@ -37,6 +38,9 @@ template <typename T, typename StorageType> class PixelHistogramImpl {
     void fill(const NDView<T, 2> &frame);
     void fill(int row, int col, T value);
     void fill_unchecked(int row, int col, T value);
+    // Fill using a row-major pixel index, avoiding repeated 3D stride
+    // calculation in tiled worker loops. The index is intentionally unchecked.
+    void fill_flat_unchecked(std::size_t pixel, T value);
 
     NDArray<StorageType, 3> values() const;
     // Zero-copy view of the underlying [rows x cols x n_bins] storage.
@@ -98,22 +102,29 @@ void PixelHistogramImpl<T, StorageType>::fill(int row, int col, T value) {
 template <typename T, typename StorageType>
 void PixelHistogramImpl<T, StorageType>::fill_unchecked(int row, int col,
                                                         T value) {
+    const auto pixel =
+        static_cast<std::size_t>(row) * static_cast<std::size_t>(m_cols) +
+        static_cast<std::size_t>(col);
+    fill_flat_unchecked(pixel, value);
+}
+
+template <typename T, typename StorageType>
+void PixelHistogramImpl<T, StorageType>::fill_flat_unchecked(std::size_t pixel,
+                                                             T value) {
     if (value < m_xmin || value >= m_xmax) {
         return;
     }
     int bin = static_cast<int>((value - m_xmin) * m_scale);
     // Guard against floating-point rounding pushing val just below
-    // xmax to bin == n_bins.
-    if (bin >= m_n_bins) {
-        bin = m_n_bins - 1;
-    }
+    bin = std::clamp(bin, 0, m_n_bins - 1);
+    auto &cell = m_values.data()[pixel * static_cast<std::size_t>(m_n_bins) +
+                                 static_cast<std::size_t>(bin)];
     if constexpr (std::is_integral_v<StorageType>) {
-        if (m_values(row, col, bin) >=
-            std::numeric_limits<StorageType>::max()) {
+        if (cell >= std::numeric_limits<StorageType>::max()) {
             return;
         }
     }
-    ++m_values(row, col, bin);
+    ++cell;
 }
 
 template <typename T, typename StorageType>

python/CMakeLists.txt

@@ -39,7 +39,9 @@ set(PYTHON_FILES
     aare/ClusterFinder.py
     aare/ClusterVector.py
     aare/Cluster.py
+    aare/FastPedestal.py
     aare/calibration.py
+    aare/factory.py
     aare/experimental.py
     aare/RawFile.py
     aare/transform.py

python/aare/Cluster.py

@@ -1,6 +1,6 @@
 from . import _aare 
 import numpy as np
-from .ClusterFinder import _type_to_char
+from .factory import _type_to_char
 
 
 def Cluster(x : int, y : int, data, cluster_size=(3,3), dtype = np.int32):
@@ -21,4 +21,4 @@ def Cluster(x : int, y : int, data, cluster_size=(3,3), dtype = np.int32):
     except AttributeError:
         raise ValueError(f"Unsupported combination of type and cluster size: {dtype}/{cluster_size} when requesting {class_name}")
    
-    return cls(x, y, data) 
+    return cls(x, y, data)

python/aare/ClusterFinder.py

@@ -1,21 +1,10 @@
 # SPDX-License-Identifier: MPL-2.0
 from . import _aare
 import numpy as np
+from .factory import _type_to_char
 
 _supported_cluster_sizes = [(2,2), (3,3), (5,5), (7,7), (9,9),]
 
-def _type_to_char(dtype):
-    if dtype == np.int32:
-        return 'i'
-    elif dtype == np.float32:
-        return 'f'
-    elif dtype == np.float64:
-        return 'd'
-    elif dtype == np.int16:
-        return 'i16'
-    else:
-        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32, np.float32, and np.float64 are supported.")
-
 def _get_class(name, cluster_size, dtype):
     """
     Helper function to get the class based on the name, cluster size, and dtype.

python/aare/FastPedestal.py

@@ -0,0 +1,56 @@
+# SPDX-License-Identifier: MPL-2.0
+
+import numpy as np
+
+from .factory import _get_typed_class
+
+
+def _get_fast_pedestal_class(dtype):
+    return _get_typed_class("FastPedestal", dtype)
+
+
+def FastPedestal(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 ``add_init_frame()`` exactly ``n_samples`` times before using the
+    statistics or calling ``push_ema()``. Subsequent frames have weight
+    ``1 / n_samples`` in the running mean and population variance.
+
+    Args:
+        rows: Number of image rows.
+        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.
+            Supported values are ``np.float64``, ``np.float32``, and
+            ``np.int16``.
+    """
+    cls = _get_fast_pedestal_class(dtype)
+    return cls(rows, cols, n_samples)
+
+
+def from_file(filename, n_samples=1000, skip_first=0, dtype=np.float64):
+    """Create a FastPedestal from frames in a file.
+
+    After ignoring ``skip_first`` frames, the next ``n_samples`` frames
+    initialize the pedestal. Every remaining frame is then applied as a
+    steady-state update. Input frames are read as uint16 data.
+
+    Args:
+        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.
+
+    Raises:
+        RuntimeError: If fewer than ``n_samples`` frames remain after
+            ``skip_first`` or ``n_samples`` is zero.
+    """
+    cls = _get_fast_pedestal_class(dtype)
+    return cls.from_file(
+        filename, n_samples=n_samples, skip_first=skip_first
+    )
+
+
+FastPedestal.from_file = from_file

python/aare/__init__.py

@@ -5,6 +5,16 @@ from . import _aare
 from . import transform
 from . import experimental
 
+from ._aare import (
+    FastPedestal_d,
+    FastPedestal_f,
+    FastPedestal_i16,
+    Pedestal_d,
+    Pedestal_f,
+    Pedestal_i16,
+    ClusterFinder_Cluster3x3i,
+    VarClusterFinder,
+)
 from ._aare import (
     File,
     JungfrauDataFile,
@@ -20,6 +30,7 @@ from ._aare import corner
 # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
 
 from ._version import __version__
+from .FastPedestal import FastPedestal
 from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile
 from .ClusterVector import ClusterVector
 from .Cluster import Cluster

python/aare/factory.py

@@ -0,0 +1,36 @@
+# SPDX-License-Identifier: MPL-2.0
+
+import numpy as np
+
+from . import _aare
+
+
+_TYPE_TO_CHAR = {
+    np.dtype(np.int32): "i",
+    np.dtype(np.float32): "f",
+    np.dtype(np.float64): "d",
+    np.dtype(np.int16): "i16",
+}
+
+
+def _type_to_char(dtype):
+    """Return the suffix used by bindings instantiated for ``dtype``."""
+    try:
+        return _TYPE_TO_CHAR[np.dtype(dtype)]
+    except (KeyError, TypeError):
+        supported = ", ".join(str(dtype) for dtype in _TYPE_TO_CHAR)
+        raise ValueError(
+            f"Unsupported dtype: {dtype}. Supported dtypes are {supported}."
+        ) from None
+
+
+def _get_typed_class(name, dtype):
+    """Return a bound class named ``<name>_<dtype suffix>``."""
+    class_name = f"{name}_{_type_to_char(dtype)}"
+    try:
+        return getattr(_aare, class_name)
+    except AttributeError:
+        raise ValueError(
+            f"Unsupported dtype for {name}: {dtype} "
+            f"(binding {class_name} is not available)."
+        ) from None

python/src/bind_Cluster.hpp

@@ -10,7 +10,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 

python/src/bind_ClusterCollector.hpp

@@ -6,6 +6,8 @@
 #include "aare/ClusterVector.hpp"
 #include "aare/NDView.hpp"
 #include "aare/Pedestal.hpp"
+
+#include "module_config.hpp"
 #include "np_helper.hpp"
 
 #include <cstdint>
@@ -15,7 +17,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 
@@ -30,7 +31,7 @@ void define_ClusterCollector(py::module &m, const std::string &typestr) {
     using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
 
     py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
-        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *>())
+        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, pd_type> *>())
         .def("stop", &ClusterCollector<ClusterType>::stop)
         .def(
             "steal_clusters",

python/src/bind_ClusterFileSink.hpp

@@ -15,8 +15,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
-
 using namespace aare;
 
 #pragma GCC diagnostic push
@@ -29,6 +27,8 @@ void define_ClusterFileSink(py::module &m, const std::string &typestr) {
 
     using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
 
+    // TODO! adapt to set pedestal type (needs templating of ClusterFileSink)
+    // or maybe access through base class?
     py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
         .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *,
                       const std::filesystem::path &>())

python/src/bind_ClusterFinder.hpp

@@ -6,6 +6,8 @@
 #include "aare/ClusterVector.hpp"
 #include "aare/NDView.hpp"
 #include "aare/Pedestal.hpp"
+
+#include "module_config.hpp"
 #include "np_helper.hpp"
 
 #include <cstdint>
@@ -15,7 +17,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 
@@ -48,6 +49,8 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
              })
         .def("clear_pedestal",
              &ClusterFinder<ClusterType, uint16_t, pd_type>::clear_pedestal)
+        .def("update_threshold",
+             &ClusterFinder<ClusterType, uint16_t, pd_type>::update_threshold)
         .def_property_readonly(
             "pedestal",
             [](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {

python/src/bind_ClusterFinderMT.hpp

@@ -6,6 +6,8 @@
 #include "aare/ClusterVector.hpp"
 #include "aare/NDView.hpp"
 #include "aare/Pedestal.hpp"
+
+#include "module_config.hpp"
 #include "np_helper.hpp"
 
 #include <cstdint>
@@ -15,7 +17,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 
@@ -31,9 +32,10 @@ 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>(),
+        .def(py::init<Shape<2>, pd_type, 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("capacity") = 2048, py::arg("n_threads") = 3,
+             py::arg("queue_depth") = 16)
         .def("push_pedestal_frame",
              [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
                 py::array_t<uint16_t> frame) {
@@ -46,7 +48,6 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
                py::array_t<uint16_t> frame, uint64_t frame_number) {
                 auto view = make_view_2d(frame);
                 self.find_clusters(view, frame_number);
-                return;
             },
             py::arg(), py::arg("frame_number") = 0)
         .def_property_readonly(
@@ -56,6 +57,8 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
             })
         .def("clear_pedestal",
              &ClusterFinderMT<ClusterType, uint16_t, pd_type>::clear_pedestal)
+        .def("update_threshold",
+             &ClusterFinderMT<ClusterType, uint16_t, pd_type>::update_threshold)
         .def("sync", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::sync)
         .def("stop", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::stop)
         .def("start", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::start)

python/src/bind_ClusterVector.hpp

@@ -15,7 +15,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 

python/src/bind_FastPedestal.hpp

@@ -0,0 +1,147 @@
+// SPDX-License-Identifier: MPL-2.0
+
+#include "aare/FastPedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+template <typename SUM_TYPE>
+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.",
+        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. It becomes ready after n_samples "
+             "calls to add_init_frame().")
+        .def(py::init<uint32_t, uint32_t>(), py::arg("rows"), py::arg("cols"),
+             "Construct an empty pedestal with n_samples=1000.")
+
+        .def(
+            "mean",
+            [](FastPedestal<SUM_TYPE> &self) {
+                auto mean = new NDArray<SUM_TYPE, 2>{};
+                *mean = self.mean();
+                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) {
+                auto standard_deviation = new NDArray<SUM_TYPE, 2>{};
+                *standard_deviation = self.std();
+                return return_image_data(standard_deviation);
+            },
+            "Return the population standard deviation as a NumPy array. The "
+            "pedestal must be ready.")
+        .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. "
+            "The pedestal must be ready.")
+
+        // We need to buffer protocol to allow for numpy operations using the
+        // pedestal mean
+        .def_buffer([](FastPedestal<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);
+        })
+        // Subtracting a FastPedestal from a NumPy array
+        .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 FastPedestal from a NumPy array.")
+        .def("clear", py::overload_cast<>(&FastPedestal<SUM_TYPE>::clear),
+             "Reset all statistics and initialization state to zero.")
+        .def_property_readonly("rows", &FastPedestal<SUM_TYPE>::rows,
+                               "Number of image rows.")
+        .def_property_readonly("cols", &FastPedestal<SUM_TYPE>::cols,
+                               "Number of image columns.")
+        .def_property_readonly("cur_samples",
+                               &FastPedestal<SUM_TYPE>::cur_samples,
+                               "Number of initialization frames accumulated. "
+                               "Steady-state pushes do not change it.")
+        .def_property_readonly(
+            "ready", &FastPedestal<SUM_TYPE>::ready,
+            "Whether n_samples initialization frames have been accumulated.")
+        .def_property_readonly(
+            "n_samples", &FastPedestal<SUM_TYPE>::n_samples,
+            "Initialization frame count and steady-state update-weight "
+            "denominator.")
+        .def(
+            "clone",
+            [](FastPedestal<SUM_TYPE> &pedestal) {
+                return FastPedestal<SUM_TYPE>(pedestal);
+            },
+            "Return an independent copy of the pedestal and its state.")
+        .def(
+            "push_ema",
+            [](FastPedestal<SUM_TYPE> &pedestal,
+               py::array_t<uint16_t, py::array::c_style> &frame) {
+                pedestal.push_ema(make_view_2d(frame));
+            },
+            py::arg("frame").noconvert(),
+            "Update exponential moving average. The pedstal must  "
+            "already be ready for this update.")
+        .def(
+            "add_init_frame",
+            [](FastPedestal<SUM_TYPE> &pedestal,
+               py::array_t<uint16_t, py::array::c_style> &frame) {
+                pedestal.add_init_frame(make_view_2d(frame));
+            },
+            py::arg("frame").noconvert(),
+            "Accumulate one uint16 initialization frame. Call exactly "
+            "n_samples times to make the pedestal ready.")
+        .def_static(
+            "from_file",
+            [](const std::filesystem::path &filename, uint32_t n_samples,
+               uint32_t skip_first) {
+                return FastPedestal<SUM_TYPE>::template from_file<uint16_t>(
+                    filename, n_samples, skip_first);
+            },
+            py::arg("filename"), py::arg("n_samples") = 1000,
+            py::arg("skip_first") = 0,
+            "Create a pedestal from a uint16 file. Skip skip_first frames, use "
+            "the next n_samples for initialization, then apply every remaining "
+            "frame as a steady-state update.");
+}

python/src/bind_PedestalTrackingPixelHistogram.hpp

@@ -153,9 +153,13 @@ void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) {
              Args:
                  file_path: Path to the file to fill from
                  max_frames: Maximum number of frames to fill from the file (default: -1)
+                 reader_threads: Number of parallel file reader workers (default: 2)
+                 reader_chunk_size: Frames claimed by each reader worker per batch (default: 4)
              )",
              py::call_guard<py::gil_scoped_release>(), py::arg("fname"),
-             py::arg("max_frames") = -1, py::arg("verbose") = false)
+             py::arg("max_frames") = -1, py::arg("verbose") = false,
+             py::arg("reader_threads") = std::size_t{2},
+             py::arg("reader_chunk_size") = std::size_t{4})
         .def("process_pedestal_file",
              &PedestalTrackingPixelHistogram::process_pedestal_file,
              R"(

python/src/module.cpp

@@ -1,6 +1,8 @@
 // SPDX-License-Identifier: MPL-2.0
 // Files with bindings to the different classes
 
+#include "module_config.hpp"
+
 // New style file naming
 #include "bind_Cluster.hpp"
 #include "bind_ClusterCollector.hpp"
@@ -11,6 +13,7 @@
 #include "bind_ClusterVector.hpp"
 #include "bind_Defs.hpp"
 #include "bind_Eta.hpp"
+#include "bind_FastPedestal.hpp"
 #include "bind_Interpolator.hpp"
 #include "bind_MultiThreadedFileReader.hpp"
 #include "bind_PedestalTrackingPixelHistogram.hpp"
@@ -30,6 +33,7 @@
 #include "var_cluster.hpp"
 
 // Pybind stuff
+#include <cstdint>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 
@@ -75,6 +79,10 @@ PYBIND11_MODULE(_aare, m) {
     define_pedestal_tracking_pixel_histogram_bindings(m);
     define_pedestal_bindings<double>(m, "Pedestal_d");
     define_pedestal_bindings<float>(m, "Pedestal_f");
+    define_pedestal_bindings<int16_t>(m, "Pedestal_i16");
+    define_fast_pedestal_bindings<double>(m, "FastPedestal_d");
+    define_fast_pedestal_bindings<float>(m, "FastPedestal_f");
+    define_fast_pedestal_bindings<int16_t>(m, "FastPedestal_i16");
     define_fit_bindings(m);
     define_interpolation_bindings(m);
     define_jungfrau_data_file_io_bindings(m);
@@ -106,6 +114,7 @@ PYBIND11_MODULE(_aare, m) {
     DEFINE_BINDINGS_CLUSTERFINDER(int, 3, 3, uint16_t, i);
     DEFINE_BINDINGS_CLUSTERFINDER(double, 3, 3, uint16_t, d);
     DEFINE_BINDINGS_CLUSTERFINDER(float, 3, 3, uint16_t, f);
+    DEFINE_BINDINGS_CLUSTERFINDER(int16_t, 3, 3, uint16_t, i16);
 
     DEFINE_BINDINGS_CLUSTERFINDER(int, 5, 5, uint16_t, i);
     DEFINE_BINDINGS_CLUSTERFINDER(double, 5, 5, uint16_t, d);

python/src/module_config.hpp

@@ -0,0 +1,5 @@
+#pragma once
+#include <cstdint>
+
+// Configure module wide pedestal type for cluster finding
+using pd_type = double;

python/src/pedestal.hpp

@@ -52,6 +52,12 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
                  *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,

python/tests/test_Cluster.py

@@ -103,15 +103,18 @@ def test_max_sum():
 
 def test_cluster_finder(): 
     """Test ClusterFinder""" 
+    shape = [100,100]
+    cf = _aare.ClusterFinder_Cluster3x3i(shape)
 
-    clusterfinder = _aare.ClusterFinder_Cluster3x3i([100,100])
+    #Push 1000 frames to the pedestal
+    for i in range(1000):
+        frame = np.random.normal(loc = 100, scale = 5, size = shape).astype(np.uint16)
+        cf.push_pedestal_frame(frame)
+    cf.update_threshold()
+    frame = np.zeros(shape=shape, dtype=np.uint16)
+    cf.find_clusters(frame)
 
-    #frame = np.random.rand(100,100)
-    frame = np.zeros(shape=[100,100])
-
-    clusterfinder.find_clusters(frame)
-
-    clusters = clusterfinder.steal_clusters(False) #conversion does not work
+    clusters = cf.steal_clusters(False) #conversion does not work
 
     assert clusters.size == 0
 

python/tests/test_FastPedestal.py

@@ -0,0 +1,145 @@
+import numpy as np
+import pytest
+
+from aare import (
+    FastPedestal,
+    FastPedestal_d,
+    FastPedestal_f,
+    FastPedestal_i16,
+)
+
+
+@pytest.mark.parametrize(
+    ("dtype", "pedestal_type"),
+    [
+        (np.float64, FastPedestal_d),
+        (np.float32, FastPedestal_f),
+        (np.int16, FastPedestal_i16),
+    ],
+)
+def test_fast_pedestal_factory(dtype, pedestal_type):
+    pedestal = FastPedestal(2, 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_fast_pedestal_factory_defaults_to_double():
+    assert isinstance(FastPedestal(2, 3), FastPedestal_d)
+
+
+def test_fast_pedestal_factory_rejects_unbound_dtype():
+    with pytest.raises(ValueError, match="Unsupported dtype for FastPedestal"):
+        FastPedestal(2, 3, dtype=np.int32)
+
+
+@pytest.mark.parametrize(
+    ("kwargs", "expected_n_samples"),
+    [
+        ({"rows": 2, "cols": 3}, 1000),
+        ({"rows": 2, "cols": 3, "n_samples": 4}, 4),
+    ],
+)
+def test_fast_pedestal_binding_accepts_constructor_keywords(
+    kwargs, expected_n_samples
+):
+    pedestal = FastPedestal_d(**kwargs)
+
+    assert pedestal.rows == 2
+    assert pedestal.cols == 3
+    assert pedestal.n_samples == expected_n_samples
+
+
+@pytest.mark.parametrize(
+    ("dtype", "pedestal_type", "expected_dtype"),
+    [
+        (np.float64, FastPedestal_d, np.float64),
+        (np.float32, FastPedestal_f, np.float32),
+        (np.int16, FastPedestal_i16, np.int16),
+    ],
+)
+def test_fast_pedestal_factory_from_file(
+    tmp_path, dtype, pedestal_type, expected_dtype
+):
+    frames = np.array(
+        [[[100, 100]], [[2, 4]], [[4, 6]], [[5, 7]]], dtype=np.uint16
+    )
+    filename = tmp_path / "frames.npy"
+    np.save(filename, frames)
+
+    pedestal = FastPedestal.from_file(
+        filename, n_samples=2, skip_first=1, dtype=dtype
+    )
+
+    assert isinstance(pedestal, pedestal_type)
+    assert pedestal.ready
+    assert pedestal.cur_samples == 2
+    assert pedestal.mean().dtype == expected_dtype
+    np.testing.assert_array_equal(pedestal.mean(), [[4, 6]])
+
+
+def test_fast_pedestal_factory_from_file_rejects_unbound_dtype():
+    with pytest.raises(ValueError, match="Unsupported dtype for FastPedestal"):
+        FastPedestal.from_file("unused.npy", dtype=np.int32)
+
+
+def test_fast_pedestal_from_file_rejects_skip_beyond_end(tmp_path):
+    filename = tmp_path / "frames.npy"
+    np.save(filename, np.zeros((1, 1, 1), dtype=np.uint16))
+
+    with pytest.raises(RuntimeError, match="less frames"):
+        FastPedestal.from_file(filename, n_samples=1, skip_first=2)
+
+
+@pytest.mark.parametrize(
+    ("pedestal_type", "expected_dtype"),
+    [(FastPedestal_d, np.float64), (FastPedestal_f, np.float32)],
+)
+def test_fast_pedestal_initialization(pedestal_type, expected_dtype):
+    pedestal = pedestal_type(2, 3, 2)
+    first = np.array([[2, 4, 6], [8, 10, 12]], dtype=np.uint16)
+    second = np.array([[4, 6, 8], [10, 12, 14]], dtype=np.uint16)
+
+    pedestal.add_init_frame(first)
+    pedestal.add_init_frame(second)
+
+
+    expected_mean = np.array(
+        [[3, 5, 7], [9, 11, 13]], dtype=expected_dtype
+    )
+    np.testing.assert_array_equal(pedestal.mean(), expected_mean)
+    np.testing.assert_array_equal(pedestal.std(), np.ones((2, 3)))
+
+
+def test_fast_pedestal_steady_state_push_ema():
+    pedestal = FastPedestal_d(1, 2, 2)
+    pedestal.add_init_frame(np.array([[2, 4]], dtype=np.uint16))
+    pedestal.add_init_frame(np.array([[4, 6]], dtype=np.uint16))
+
+
+    pedestal.push_ema(np.array([[6, 8]], dtype=np.uint16))
+
+    np.testing.assert_array_equal(pedestal.mean(), [[4.5, 6.5]])
+
+
+def test_fast_pedestal_exposes_read_only_buffer_and_subtraction():
+    pedestal = FastPedestal_d(1, 2, 1)
+    pedestal.add_init_frame(np.array([[2, 4]], dtype=np.uint16))
+
+
+    view = np.asarray(pedestal)
+    result = np.array([[12, 14]], dtype=np.uint16) - pedestal
+
+    np.testing.assert_array_equal(view, [[2, 4]])
+    np.testing.assert_array_equal(result, [[10, 10]])
+    assert np.shares_memory(view, pedestal.view())
+    assert not view.flags.writeable
+
+
+def test_fast_pedestal_rejects_wrong_shape():
+    pedestal = FastPedestal_d(2, 3)
+
+    with pytest.raises(RuntimeError, match="shape"):
+        pedestal.add_init_frame(np.zeros((2, 2), dtype=np.uint16))

src/ClusterFinderMT.test.cpp

@@ -21,14 +21,21 @@ class ClusterFinderMTWrapper
 
   public:
     ClusterFinderMTWrapper(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
-                           size_t capacity = 2000, size_t n_threads = 3)
+                           size_t capacity = 2000, size_t n_threads = 3,
+                           size_t queue_depth = 16)
         : ClusterFinderMT<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>(
-              image_size, nSigma, capacity, n_threads) {}
+              image_size, nSigma, capacity, n_threads, queue_depth) {}
 
     size_t get_m_input_queues_size() const {
         return this->m_input_queues.size();
     }
 
+    size_t get_m_frame_pools_size() const { return this->m_frame_pools.size(); }
+
+    size_t get_frame_pool_depth(size_t thread_index) const {
+        return this->m_frame_pools[thread_index]->size();
+    }
+
     size_t get_m_output_queues_size() const {
         return this->m_output_queues.size();
     }
@@ -100,3 +107,73 @@ TEST_CASE("multithreaded cluster finder", "[.with-data]") {
     auto clustervec = clustercollector.steal_clusters();
     // CHECK(clustervec.size() == ) //dont know how many clusters to expect
 }
+
+TEST_CASE("frame buffers are recycled when pushing more frames than the pool "
+          "holds",
+          "[.files]") {
+    using ClusterType = Cluster<int32_t, 3, 3>;
+
+    const size_t n_threads = 2;
+    const size_t queue_depth = 4;
+    const size_t n_frames = 5 * queue_depth;
+    const Shape<2> image_size{10, 10};
+
+    ClusterFinderMTWrapper<ClusterType> cf(image_size, 5, 200, n_threads,
+                                           queue_depth);
+
+    CHECK(cf.get_m_frame_pools_size() == n_threads);
+    for (size_t i = 0; i < n_threads; ++i) {
+        CHECK(cf.get_frame_pool_depth(i) == queue_depth);
+    }
+
+    NDArray<uint16_t, 2> frame(image_size, 0);
+
+    // More frames than the pool holds, so this only completes if the workers
+    // return the buffers to the free list.
+    for (size_t i = 0; i < n_frames; ++i) {
+        cf.find_clusters(frame.view(), i);
+    }
+
+    cf.stop();
+
+    CHECK(cf.m_input_queues_are_empty() == true);
+}
+
+TEST_CASE("cluster collector accepts finders with a matching cluster type") {
+    using ClusterType = Cluster<int32_t, 3, 3>;
+    using Finder = ClusterFinderMTWrapper<ClusterType, uint16_t, float>;
+    using OtherFinder =
+        ClusterFinderMTWrapper<Cluster<int32_t, 5, 5>, uint16_t, float>;
+
+    static_assert(
+        std::is_constructible_v<ClusterCollector<ClusterType>, Finder *>);
+    static_assert(
+        !std::is_constructible_v<ClusterCollector<ClusterType>, OtherFinder *>);
+
+    Finder cf({10, 10});
+    cf.stop();
+
+    ClusterCollector<ClusterType> collector(&cf);
+    collector.stop();
+
+    CHECK(collector.steal_clusters().empty());
+}
+
+TEST_CASE("cluster collector drains queued clusters when stopped") {
+    using ClusterType = Cluster<int32_t, 3, 3>;
+    ProducerConsumerQueue<ClusterVector<ClusterType>> source(4);
+
+    for (uint64_t frame_number = 1; frame_number <= 3; ++frame_number) {
+        REQUIRE(source.write(ClusterVector<ClusterType>(4, frame_number)));
+    }
+
+    ClusterCollector<ClusterType> collector(&source);
+    collector.stop();
+
+    auto clusters = collector.steal_clusters();
+    REQUIRE(clusters.size() == 3);
+    CHECK(source.isEmpty());
+    for (size_t i = 0; i < clusters.size(); ++i) {
+        CHECK(clusters[i].frame_number() == static_cast<int32_t>(i + 1));
+    }
+}

src/hist/PedestalTrackingPixelHistogram.cpp

@@ -1,11 +1,16 @@
 #include "aare/hist/PedestalTrackingPixelHistogram.hpp"
 #include "aare/File.hpp"
+#include "aare/MultiThreadedFileReader.hpp"
 
 #include <algorithm>
+#include <array>
 #include <chrono>
 #include <cmath>
 #include <cstring>
+#include <functional>
+#include <future>
 #include <stdexcept>
+#include <type_traits>
 #include <utility>
 #include <vector>
 
@@ -202,17 +207,8 @@ void PedestalTrackingPixelHistogram::worker_loop(int thread_id) {
 
         switch (kind) {
         case WorkKind::PushPedestal: {
-            // Accumulate raw frame values into this thread's pedestal
-            // shard. Uses the pixel-level push_no_update which only
-            // touches m_sum/m_sum2/m_cur_samples (no m_mean writes).
-            for (int local_row = 0; local_row < local_rows; ++local_row) {
-                const auto row = static_cast<ssize_t>(first_row + local_row);
-                for (ssize_t col = 0; col < image->shape(1); ++col) {
-                    my_pedestal.template push_no_update<FrameType>(
-                        static_cast<uint32_t>(local_row),
-                        static_cast<uint32_t>(col), (*image)(row, col));
-                }
-            }
+            auto frame = image->sub_view(first_row, first_row + local_rows);
+            my_pedestal.add_init_frame(frame);
             break;
         }
         case WorkKind::UpdateMean: {
@@ -220,7 +216,6 @@ void PedestalTrackingPixelHistogram::worker_loop(int thread_id) {
             // thread's shard. Also refresh the cached per-pixel std so
             // FillWithThreshold can read it without recomputing on the
             // hot path.
-            my_pedestal.update_mean();
             auto &my_std = partial_std_[thread_id];
             for (int local_row = 0; local_row < local_rows; ++local_row) {
                 for (int col = 0; col < cols_; ++col) {
@@ -240,59 +235,63 @@ void PedestalTrackingPixelHistogram::worker_loop(int thread_id) {
             // tracking gate entirely. The [xmin, xmax) histogram gate
             // lives inside PixelHistogramImpl::fill.
             const auto n_sigma = n_sigma_.load(std::memory_order_relaxed);
-            if (n_sigma <= AxisType{0.0}) {
-                // Fill without pedestal tracking.
-                for (const auto &frame : *images) {
-                    const auto cols = frame.shape(1);
-                    for (int local_row = 0; local_row < local_rows;
-                         ++local_row) {
-                        const auto row =
-                            static_cast<ssize_t>(first_row + local_row);
-                        for (ssize_t col = 0; col < cols; ++col) {
-                            const FrameType raw = frame(row, col);
-                            const AxisType val =
-                                static_cast<AxisType>(raw) -
-                                static_cast<AxisType>(my_pedestal.mean(
-                                    static_cast<uint32_t>(local_row),
-                                    static_cast<uint32_t>(col)));
-                            my_hist.fill_unchecked(local_row,
-                                                   static_cast<int>(col), val);
-                        }
-                    }
-                }
-                break;
-            } else {
-                // Do pedestal tracking. Duplicated code for clean hot path.
+            constexpr std::size_t pixel_tile_size = 256;
+            const auto local_pixels = static_cast<std::size_t>(local_rows) *
+                                      static_cast<std::size_t>(cols_);
+            const auto global_pixel_begin =
+                static_cast<std::size_t>(first_row) *
+                static_cast<std::size_t>(cols_);
+            const auto &frames = *images;
 
-                auto &my_std = partial_std_[thread_id];
-                for (const auto &frame : *images) {
-                    const auto cols = frame.shape(1);
-                    for (int local_row = 0; local_row < local_rows;
-                         ++local_row) {
-                        const auto row =
-                            static_cast<ssize_t>(first_row + local_row);
-                        for (ssize_t col = 0; col < cols; ++col) {
-                            const FrameType raw = frame(row, col);
+            // Compile the shared traversal once for each mode. if constexpr
+            // removes the tracking gate entirely from the histogram-only hot
+            // path, while keeping the runtime mode decision outside the
+            // per-pixel loops.
+            const auto fill_tiles = [&](auto tracking,
+                                        const AxisType *std_data) {
+                constexpr bool track_pedestal = decltype(tracking)::value;
+
+                // Each worker retains exclusive ownership of its row shard.
+                // Pixel tiling keeps a bounded part of that shard's
+                // [pixel x bin] storage hot while input remains contiguous.
+                for (std::size_t tile_begin = 0; tile_begin < local_pixels;
+                     tile_begin += pixel_tile_size) {
+                    const auto tile_end =
+                        std::min(tile_begin + pixel_tile_size, local_pixels);
+                    for (const auto &frame : frames) {
+                        const auto *input =
+                            frame.data() + global_pixel_begin + tile_begin;
+                        for (std::size_t local_pixel = tile_begin;
+                             local_pixel < tile_end; ++local_pixel, ++input) {
+                            const FrameType raw = *input;
                             const AxisType val =
                                 static_cast<AxisType>(raw) -
-                                static_cast<AxisType>(my_pedestal.mean(
-                                    static_cast<uint32_t>(local_row),
-                                    static_cast<uint32_t>(col)));
-                            my_hist.fill_unchecked(local_row,
-                                                   static_cast<int>(col), val);
-                            const AxisType sigma = my_std(local_row, col);
-                            if (sigma > AxisType{0.0} &&
-                                std::abs(static_cast<AxisType>(val)) <
-                                    n_sigma * sigma) {
-                                my_pedestal.template push<FrameType>(
-                                    static_cast<uint32_t>(local_row),
-                                    static_cast<uint32_t>(col), raw);
+                                my_pedestal.mean_unchecked(
+                                    static_cast<ssize_t>(local_pixel));
+                            my_hist.fill_flat_unchecked(local_pixel, val);
+
+                            if constexpr (track_pedestal) {
+                                const AxisType sigma = std_data[local_pixel];
+                                if (sigma > AxisType{0.0} &&
+                                    std::abs(val) < n_sigma * sigma) {
+                                    my_pedestal.push_ema_unchecked<FrameType>(
+                                        local_pixel, raw);
+                                }
                             }
                         }
                     }
                 }
-                break;
+            };
+
+            if (n_sigma <= AxisType{0.0}) {
+                fill_tiles(std::false_type{}, nullptr);
+            } else {
+                // Frames stay chronological for every pixel. Accepted raw
+                // values update that pixel's sums and cached mean immediately;
+                // pixels are otherwise independent.
+                fill_tiles(std::true_type{}, partial_std_[thread_id].data());
             }
+            break;
         }
         }
 
@@ -395,6 +394,11 @@ void PedestalTrackingPixelHistogram::fill_async(NDArray<FrameType, 2> &&image) {
             "constructor shape");
     }
 
+    // ProducerConsumerQueue is SPSC. Serialising this short producer-side
+    // operation also prevents fill_async from racing fill_from_file's direct
+    // batch dispatch.
+    std::lock_guard<std::mutex> ingestion_lock(ingestion_mutex_);
+
     // SPSC backpressure: spin with a short sleep until a slot frees up.
     // The std::move only consumes `image` on the iteration that succeeds
     // (placement-new inside write() runs only when the slot is free).
@@ -467,35 +471,42 @@ PedestalTrackingPixelHistogram::bin_edges() const {
 }
 
 void PedestalTrackingPixelHistogram::fill_from_file(
-    const std::filesystem::path &fname, ssize_t max_frames, bool verbose) {
+    const std::filesystem::path &fname, ssize_t max_frames, bool verbose,
+    std::size_t reader_threads, std::size_t reader_chunk_size) {
     constexpr std::size_t progress_interval = 66;
     auto last = std::chrono::steady_clock::now();
-    const auto completed_start =
-        completed_async_fills_.load(std::memory_order_acquire);
     std::size_t last_reported = 0;
-    const auto completed_for_this_file = [&]() {
-        return completed_async_fills_.load(std::memory_order_acquire) -
-               completed_start;
-    };
 
-    const auto wait_for_completed = [&](std::size_t target) {
-        while (completed_for_this_file() < target) {
-            std::this_thread::sleep_for(async_wait_);
-        }
-    };
+    if (max_frames < -1) {
+        throw std::invalid_argument(
+            "PedestalTrackingPixelHistogram max_frames must be -1 or "
+            "non-negative");
+    }
 
-    File f(fname);
-    // check that row col matches constructor
-    if (f.rows() != static_cast<size_t>(rows_) ||
-        f.cols() != static_cast<size_t>(cols_)) {
+    // Preserve the old max_frames behaviour: values beyond EOF are clamped
+    // rather than rejected by MultiThreadedFileReader.
+    const auto source_total_frames = File(fname).total_frames();
+    const auto n_frames = max_frames == -1
+                              ? source_total_frames
+                              : std::min(static_cast<std::size_t>(max_frames),
+                                         source_total_frames);
+
+    experimental::MultiThreadedFileReader reader(fname, reader_threads,
+                                                 reader_chunk_size, n_frames);
+    if (reader.rows() != static_cast<size_t>(rows_) ||
+        reader.cols() != static_cast<size_t>(cols_)) {
         throw std::invalid_argument("PedestalTrackingPixelHistogram: Frame in "
                                     "file {} has shape ({}, {}) does not match "
                                     "constructor shape");
     }
+    if (reader.dtype() != Dtype::UINT16 ||
+        reader.bytes_per_frame() != static_cast<std::size_t>(rows_) *
+                                        static_cast<std::size_t>(cols_) *
+                                        sizeof(FrameType)) {
+        throw std::invalid_argument(
+            "PedestalTrackingPixelHistogram requires uint16 file frames");
+    }
 
-    const ssize_t total_frames = f.total_frames();
-    const ssize_t n_frames =
-        max_frames == -1 ? total_frames : std::min(max_frames, total_frames);
     const auto print_progress = [&](std::size_t done) {
         const auto now = std::chrono::steady_clock::now();
         const double dt = std::chrono::duration<double>(now - last).count();
@@ -503,30 +514,89 @@ void PedestalTrackingPixelHistogram::fill_from_file(
         const double fps =
             dt > 0.0 ? static_cast<double>(done_in_interval) / dt : 0.0;
 
-        fmt::print(
-            "\rProgress: {}/{} ({:.1f}%)  {:.1f} FPS    ", done, n_frames,
-            100.0 * static_cast<double>(done) / static_cast<double>(n_frames),
-            fps);
+        fmt::print("\rProgress: {}/{} ({:.1f}%)  {:.1f} FPS    ", done,
+                   n_frames,
+                   n_frames == 0 ? 100.0
+                                 : 100.0 * static_cast<double>(done) /
+                                       static_cast<double>(n_frames),
+                   fps);
         std::fflush(stdout);
         last = now;
         last_reported = done;
     };
 
-    for (ssize_t i = 0; i < n_frames; ++i) {
-        aare::NDArray<uint16_t> frame({rows_, cols_});
-        f.read_into(reinterpret_cast<std::byte *>(frame.data()));
-        fill_async(std::move(frame));
+    using ReadBuffer = NDArray<FrameType, 3>;
+    const auto read_into_buffer = [&reader](ReadBuffer &buffer) {
+        const auto frame_count = reader.next_read_frames();
+        const auto frames_read = reader.read_into(buffer.buffer());
+        if (frames_read != frame_count) {
+            throw std::runtime_error(
+                "MultiThreadedFileReader returned an incomplete batch");
+        }
+        return frames_read;
+    };
 
-        if (verbose && (i + 1) % progress_interval == 0) {
-            wait_for_completed(static_cast<std::size_t>(i + 1));
-            print_progress(static_cast<std::size_t>(i + 1));
+    const auto fill_batch = [this](ReadBuffer &buffer,
+                                   std::size_t frame_count) {
+        std::vector<NDView<FrameType, 2>> views;
+        views.reserve(frame_count);
+        auto batch_view = buffer.view();
+        for (std::size_t i = 0; i < frame_count; ++i) {
+            views.push_back(batch_view(i));
+        }
+        std::lock_guard<std::mutex> fill_lock(fill_mutex_);
+        dispatch_fill_batch_(views);
+    };
+
+    // Exclude other queue producers for the duration. Drain anything already
+    // submitted before bypassing the coordinator with direct batch dispatch.
+    std::lock_guard<std::mutex> ingestion_lock(ingestion_mutex_);
+    flush();
+
+    std::size_t completed = 0;
+    if (n_frames != 0) {
+        // Allocate the maximum wave size once per buffer. The last read may
+        // contain fewer frames; its returned count limits the views passed to
+        // the histogram workers, leaving the unused tail untouched.
+        const auto buffer_capacity = reader.next_read_frames();
+        std::array<ReadBuffer, 2> buffers{
+            ReadBuffer({static_cast<ssize_t>(buffer_capacity), rows_, cols_}),
+            ReadBuffer({static_cast<ssize_t>(buffer_capacity), rows_, cols_})};
+
+        std::size_t current_index = 0;
+        auto current_frames = read_into_buffer(buffers[current_index]);
+        while (true) {
+            const bool has_next = completed + current_frames < n_frames;
+            const std::size_t next_index = current_index ^ std::size_t{1};
+
+            // MultiThreadedFileReader::read_into is blocking. Run the next
+            // read on a dedicated prefetch task while the current batch is
+            // processed by the histogram worker pool.
+            std::future<std::size_t> next;
+            if (has_next) {
+                next = std::async(std::launch::async, read_into_buffer,
+                                  std::ref(buffers[next_index]));
+            }
+
+            fill_batch(buffers[current_index], current_frames);
+            completed += current_frames;
+
+            if (verbose && (completed - last_reported >= progress_interval ||
+                            completed == n_frames)) {
+                print_progress(completed);
+            }
+
+            if (!has_next) {
+                break;
+            }
+            current_frames = next.get();
+            current_index = next_index;
         }
     }
-    flush();
+
     if (verbose) {
-        const auto done = completed_for_this_file();
-        if (done > last_reported) {
-            print_progress(done);
+        if (completed > last_reported || n_frames == 0) {
+            print_progress(completed);
         }
         fmt::print("\n\n");
         std::fflush(stdout);

src/hist/PedestalTrackingPixelHistogram.test.cpp

@@ -0,0 +1,289 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/hist/PedestalTrackingPixelHistogram.hpp"
+
+#include "aare/File.hpp"
+#include "aare/Frame.hpp"
+#include "aare/NumpyFile.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <algorithm>
+#include <chrono>
+#include <cstdint>
+#include <filesystem>
+#include <functional>
+#include <stdexcept>
+#include <string>
+#include <utility>
+
+using aare::FileConfig;
+using aare::Frame;
+using aare::NumpyFile;
+using aare::PedestalTrackingPixelHistogram;
+
+namespace {
+
+class TemporaryHistogramFile {
+  public:
+    using Generator =
+        std::function<std::uint16_t(std::size_t, std::size_t, std::size_t)>;
+
+    TemporaryHistogramFile(std::size_t rows = 2, std::size_t cols = 3,
+                           std::size_t frames = 10, Generator generator = {}) {
+        const auto unique =
+            std::chrono::steady_clock::now().time_since_epoch().count();
+        path_ = std::filesystem::temp_directory_path() /
+                ("aare-pedestal-hist-" + std::to_string(unique) + ".npy");
+
+        FileConfig config;
+        config.dtype = aare::Dtype::UINT16;
+        config.rows = rows;
+        config.cols = cols;
+        NumpyFile file(path_, "w", config);
+        for (std::size_t frame_index = 0; frame_index < frames; ++frame_index) {
+            Frame frame(rows, cols, config.dtype);
+            auto image = frame.view<std::uint16_t>();
+            for (ssize_t row = 0; row < image.shape(0); ++row) {
+                for (ssize_t col = 0; col < image.shape(1); ++col) {
+                    image(row, col) =
+                        generator ? generator(frame_index,
+                                              static_cast<std::size_t>(row),
+                                              static_cast<std::size_t>(col))
+                                  : static_cast<std::uint16_t>(frame_index);
+                }
+            }
+            file.write(frame);
+        }
+    }
+
+    ~TemporaryHistogramFile() { std::filesystem::remove(path_); }
+
+    TemporaryHistogramFile(const TemporaryHistogramFile &) = delete;
+    TemporaryHistogramFile &operator=(const TemporaryHistogramFile &) = delete;
+
+    const std::filesystem::path &path() const { return path_; }
+
+  private:
+    std::filesystem::path path_;
+};
+
+} // namespace
+
+TEST_CASE("Pedestal tracking histogram fills ordered multi-reader batches",
+          "[PedestalTrackingPixelHistogram]") {
+    TemporaryHistogramFile file;
+    PedestalTrackingPixelHistogram histogram(2, 3, 10, 0.0f, 10.0f, 2, 4, 0.0f);
+
+    // Two reader workers claiming two frames each produces a full batch of
+    // four followed by a partial batch of three.
+    histogram.fill_from_file(file.path(), 7, false, 2, 2);
+    const auto values = histogram.values();
+
+    for (ssize_t row = 0; row < 2; ++row) {
+        for (ssize_t col = 0; col < 3; ++col) {
+            for (ssize_t bin = 0; bin < 10; ++bin) {
+                CHECK(values(row, col, bin) == (bin < 7 ? 1 : 0));
+            }
+        }
+    }
+}
+
+TEST_CASE("Pedestal tracking file fill handles limits and reader options",
+          "[PedestalTrackingPixelHistogram]") {
+    TemporaryHistogramFile file;
+    PedestalTrackingPixelHistogram histogram(2, 3, 10, 0.0f, 10.0f, 2, 4, 0.0f);
+
+    CHECK_NOTHROW(histogram.fill_from_file(file.path(), 0, false, 2, 2));
+    CHECK_THROWS_AS(histogram.fill_from_file(file.path(), -2, false, 2, 2),
+                    std::invalid_argument);
+    CHECK_THROWS_AS(histogram.fill_from_file(file.path(), -1, false, 0, 2),
+                    std::invalid_argument);
+    CHECK_THROWS_AS(histogram.fill_from_file(file.path(), -1, false, 2, 0),
+                    std::invalid_argument);
+
+    // Preserve the previous API's clamp-at-EOF behaviour.
+    histogram.fill_from_file(file.path(), 100, false, 2, 3);
+    const auto values = histogram.values();
+    for (ssize_t row = 0; row < 2; ++row) {
+        for (ssize_t col = 0; col < 3; ++col) {
+            for (ssize_t bin = 0; bin < 10; ++bin) {
+                CHECK(values(row, col, bin) == 1);
+            }
+        }
+    }
+}
+
+TEST_CASE("Pedestal tracking file fill validates frame metadata",
+          "[PedestalTrackingPixelHistogram]") {
+    TemporaryHistogramFile wrong_shape(3, 3, 1);
+    PedestalTrackingPixelHistogram histogram(2, 3, 10, 0.0f, 10.0f, 1, 4, 0.0f);
+
+    CHECK_THROWS_AS(
+        histogram.fill_from_file(wrong_shape.path(), -1, false, 2, 1),
+        std::invalid_argument);
+}
+
+TEST_CASE("Histogram-only fill crosses a pixel tile without changing pedestal",
+          "[PedestalTrackingPixelHistogram]") {
+    constexpr int rows = 1;
+    constexpr int cols = 513;
+    constexpr std::size_t pedestal_samples = 1000;
+
+    const auto baseline = [](std::size_t col) {
+        return static_cast<std::uint16_t>(100 + col % 17);
+    };
+    TemporaryHistogramFile file(
+        rows, cols, 2,
+        [baseline](std::size_t frame, std::size_t, std::size_t col) {
+            return static_cast<std::uint16_t>(baseline(col) + frame + 1);
+        });
+
+    // Negative n_sigma exercises the complete histogram-only condition; zero
+    // is covered by the ordered batch tests above.
+    PedestalTrackingPixelHistogram histogram(rows, cols, 4, 0.0f, 4.0f, 1, 6,
+                                             -1.0f);
+
+    for (std::size_t seed = 0; seed < pedestal_samples; ++seed) {
+        aare::NDArray<std::uint16_t, 2> frame({rows, cols});
+        const int offset = seed % 2 == 0 ? -1 : 1;
+        for (ssize_t col = 0; col < cols; ++col) {
+            frame(0, col) = static_cast<std::uint16_t>(
+                static_cast<int>(baseline(static_cast<std::size_t>(col))) +
+                offset);
+        }
+        histogram.push_pedestal_no_update(frame.view());
+    }
+    histogram.update_mean();
+    const auto mean_before = histogram.pedestal_mean();
+
+    // One two-frame batch traverses the 512-pixel tile and its one-pixel tail.
+    histogram.fill_from_file(file.path(), -1, false, 1, 2);
+
+    const auto values = histogram.values();
+    for (ssize_t col = 0; col < cols; ++col) {
+        for (ssize_t bin = 0; bin < 4; ++bin) {
+            CHECK(values(0, col, bin) == ((bin == 1 || bin == 2) ? 1 : 0));
+        }
+    }
+
+    const auto mean_after = histogram.pedestal_mean();
+    CHECK(
+        std::equal(mean_before.begin(), mean_before.end(), mean_after.begin()));
+}
+
+TEST_CASE("Pedestal tracking threshold is strict and rejects zero sigma",
+          "[PedestalTrackingPixelHistogram]") {
+    constexpr std::size_t pedestal_samples = 1000;
+    PedestalTrackingPixelHistogram histogram(1, 3, 8, -4.0f, 4.0f, 1, 4, 2.0f);
+
+    // Pixels 0 and 1 have mean 100 and population sigma 1. Pixel 2 has the
+    // same mean and zero sigma.
+    for (std::size_t seed = 0; seed < pedestal_samples; ++seed) {
+        aare::NDArray<std::uint16_t, 2> frame({1, 3});
+        const auto noisy = static_cast<std::uint16_t>(seed % 2 == 0 ? 99 : 101);
+        frame(0, 0) = noisy;
+        frame(0, 1) = noisy;
+        frame(0, 2) = 100;
+        histogram.push_pedestal_no_update(frame.view());
+    }
+    histogram.update_mean();
+
+    const auto mean_before = histogram.pedestal_mean();
+    REQUIRE(mean_before(0, 0) == 100.0f);
+    REQUIRE(mean_before(0, 1) == 100.0f);
+    REQUIRE(mean_before(0, 2) == 100.0f);
+
+    aare::NDArray<std::uint16_t, 2> frame({1, 3});
+    frame(0, 0) = 102; // residual == 2 * sigma: excluded
+    frame(0, 1) = 101; // residual < 2 * sigma: included
+    frame(0, 2) = 101; // sigma == 0: excluded
+    histogram.fill_async(std::move(frame));
+    histogram.flush();
+
+    const auto mean_after = histogram.pedestal_mean();
+    CHECK(mean_after(0, 0) == 100.0f);
+    CHECK(mean_after(0, 1) == static_cast<float>(100001.0 * (1.0 / 1000.0)));
+    CHECK(mean_after(0, 2) == 100.0f);
+
+    // Histogramming uses each residual before a possible EMA update.
+    const auto values = histogram.values();
+    for (ssize_t bin = 0; bin < 8; ++bin) {
+        CHECK(values(0, 0, bin) == (bin == 6 ? 1 : 0));
+        CHECK(values(0, 1, bin) == (bin == 5 ? 1 : 0));
+        CHECK(values(0, 2, bin) == (bin == 5 ? 1 : 0));
+    }
+}
+
+TEST_CASE("Tiled pedestal tracking matches chronological single-frame fills",
+          "[PedestalTrackingPixelHistogram]") {
+    constexpr int rows = 5;
+    constexpr int cols = 257;
+    constexpr std::size_t frames = 40;
+    constexpr int bins = 32;
+    constexpr float xmin = -16.0f;
+    constexpr float xmax = 16.0f;
+
+    const auto baseline = [](std::size_t row, std::size_t col) {
+        return static_cast<std::uint16_t>(100 + (row + col) % 5);
+    };
+    const auto sample = [baseline](std::size_t frame, std::size_t row,
+                                   std::size_t col) {
+        const int offsets[] = {1, -1, 8, 0};
+        return static_cast<std::uint16_t>(static_cast<int>(baseline(row, col)) +
+                                          offsets[frame % 4]);
+    };
+
+    TemporaryHistogramFile file(rows, cols, frames, sample);
+    PedestalTrackingPixelHistogram tiled(rows, cols, bins, xmin, xmax, 2, 16,
+                                         2.0f);
+    PedestalTrackingPixelHistogram chronological(rows, cols, bins, xmin, xmax,
+                                                 2, 16, 2.0f);
+
+    // Seed a non-zero cached standard deviation. Five rows split over two
+    // workers make both row shards cross the 512-pixel tile boundary. Forty
+    // data frames exercise chronological processing across a sizable batch.
+    for (std::size_t seed = 0; seed < 1000; ++seed) {
+        aare::NDArray<std::uint16_t, 2> frame({rows, cols});
+        const int offset = seed % 2 == 0 ? -2 : 2;
+        for (ssize_t row = 0; row < rows; ++row) {
+            for (ssize_t col = 0; col < cols; ++col) {
+                frame(row, col) = static_cast<std::uint16_t>(
+                    static_cast<int>(baseline(row, col)) + offset);
+            }
+        }
+        tiled.push_pedestal_no_update(frame.view());
+        chronological.push_pedestal_no_update(frame.view());
+    }
+    tiled.update_mean();
+    chronological.update_mean();
+
+    // One forty-frame reader wave exercises tiled batch processing.
+    tiled.fill_from_file(file.path(), -1, false, 2, 20);
+
+    // The reference path establishes the same result one chronological frame
+    // at a time, without any cross-frame traversal reordering.
+    for (std::size_t frame_index = 0; frame_index < frames; ++frame_index) {
+        aare::NDArray<std::uint16_t, 2> frame({rows, cols});
+        for (ssize_t row = 0; row < rows; ++row) {
+            for (ssize_t col = 0; col < cols; ++col) {
+                frame(row, col) =
+                    sample(frame_index, static_cast<std::size_t>(row),
+                           static_cast<std::size_t>(col));
+            }
+        }
+        chronological.fill_async(std::move(frame));
+        chronological.flush();
+    }
+
+    const auto tiled_values = tiled.values();
+    const auto chronological_values = chronological.values();
+    REQUIRE(tiled_values.shape() == chronological_values.shape());
+    CHECK(std::equal(tiled_values.begin(), tiled_values.end(),
+                     chronological_values.begin()));
+
+    const auto tiled_mean = tiled.pedestal_mean();
+    const auto chronological_mean = chronological.pedestal_mean();
+    REQUIRE(tiled_mean.shape() == chronological_mean.shape());
+    CHECK(std::equal(tiled_mean.begin(), tiled_mean.end(),
+                     chronological_mean.begin()));
+}

src/hist/PixelHistogramImpl.test.cpp

@@ -72,6 +72,17 @@ TEST_CASE("Fill a small histogram from an NDArray") {
     REQUIRE(v(1, 1, 4) == 1);
 }
 
+TEST_CASE("Flat pixel filling uses pixel-major histogram storage") {
+    aare::PixelHistogramImpl<float, uint16_t> hist(2, 3, 4, 0.0f, 4.0f);
+
+    hist.fill_flat_unchecked(0, 0.5f);
+    hist.fill_flat_unchecked(4, 2.5f);
+
+    const auto values = hist.view();
+    CHECK(values(0, 0, 0) == 1);
+    CHECK(values(1, 1, 2) == 1);
+}
+
 TEST_CASE("Check that pixel histogram does not overflow") {
     int rows = 1;
     int cols = 1;