Cleaned up Pedestal (#370)

- renamed pedestal.hpp to bind_Pedestal.hpp for python bindings
- variance is now always calculated in double and made private
- cost correctness in a few places
- removed push_fast and explicit updates of mean from Pedestal. If
performance is needed use FastPedestal
- Pedestal no longer caches std

CMakeLists.txt

@@ -473,6 +473,7 @@ if(AARE_TESTS)
       ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinderMT.test.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
+      ${CMAKE_CURRENT_SOURCE_DIR}/src/FastPedestal.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

RELEASE.md

@@ -4,6 +4,10 @@
 
 ### New Features:
 
+- Added the Python ``Pedestal`` factory with ``dtype`` selection, matching
+  ``FastPedestal`` and defaulting to ``float64`` output.
+- Added ``FastPedestal`` in C++ and Python for per-pixel running mean
+  and population standard deviation. It supports exponentially
 - Added ``ClusterFile.frames()`` and ``ClusterFile.chunks()`` in C++ and
   Python for iteration from the current file position. Frames preserve empty
   frames and frame numbers; chunks support an optional positive size override.
@@ -26,6 +30,18 @@
 
 ### API Changes:
 
+- ``FastPedestal`` variance is now a private ``double`` intermediate. Removed
+  the C++ ``variance()``/``variance_unchecked()`` APIs and Python ``var()``.
+  Standard deviation is calculated before conversion to the output type,
+  avoiding overflow of intermediate variance for ``int16`` output. Negative
+  variance from floating-point roundoff is clamped to zero.
+- ``Pedestal`` now always accumulates sums and sums of squares in ``double``,
+  like ``FastPedestal``. Mean and standard deviation output types
+  are unchanged. Removed the C++ ``get_sum()``/``get_sum2()`` getters and
+  Python ``sum``/``sum2`` properties; internal sums are no longer exposed.
+- Removed the public ``Pedestal.variance()`` and ``cached_std()`` APIs and
+  C++ ``update_std()``. Use ``std()`` to calculate the current population
+  standard deviation; variance is now an internal implementation detail.
 - Added ``ClusterFile::read_frame(ClusterVector&)`` for allocation-reusing C++
   reads. It returns ``false`` at a clean end of file. The value-returning C++
   overload now returns ``std::optional<ClusterVector<ClusterType>>`` and
@@ -58,6 +74,16 @@
 - ``TimingMode::Auto`` changed to ``TimingMode::AUTO_TIMING``, ``TimingMode::Trigger`` changed to ``TimingMode::TRIGGER_EXPOSURE``
 
 ### Bugfixes:
+- ``Pedestal`` reports mismatched frame shapes with exceptions in all push
+  overloads, including Debug builds, instead of aborting on assertions.
+- Python ``Pedestal`` constructors reject negative dimensions and sample
+  counts instead of converting them to large unsigned values.
+- ``Pedestal`` clamps negative variance from floating-point roundoff to zero,
+  preventing NaN standard deviations for nearly constant inputs.
+- Python ``Pedestal.push()`` now requires C-contiguous ``uint16`` frames
+  without implicit conversion. Both ``push()`` and ``push_with_threshold()``
+  validate that frames and thresholds are two-dimensional before constructing
+  views, preventing incorrect results from unsupported array layouts or ranks.
 
 - Fixed a leaked empty ``ClusterVector`` at the end of Python ``ClusterFile``
   iteration. Chunk iteration now rejects a zero chunk size.

docs/src/python/pedestal/pyFastPedestal.rst

@@ -1,7 +1,7 @@
 FastPedestal
 ============
 
-``FastPedestal`` calculates a running mean, variance and standard deviation for each pixel in a
+``FastPedestal`` calculates a running mean 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 
@@ -18,7 +18,10 @@ The public factory selects the bound C++ specialization from ``dtype``:
 * ``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.
+Internal moments and variance are calculated in double precision. Variance is
+private and stays in double precision through the square root; the cached mean
+and on-demand standard deviation are returned in the specified type. Negative
+variance caused by floating-point roundoff is clamped to zero.
 
 Factory
 -------

docs/src/python/pedestal/pyPedestal.rst

@@ -1,30 +1,57 @@
 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.
+``Pedestal`` calculates a running mean and population standard deviation for
+each pixel in a series of ``uint16`` frames. ``push()`` updates the cached mean
+immediately; ``std()`` calculates the noise from the current statistics.
 
-Three specializations are available from :mod:`aare`:
+``push()`` and ``push_with_threshold()`` require C-contiguous, two-dimensional
+NumPy frames with dtype ``uint16`` and shape matching the pedestal.
+``push_with_threshold()`` also requires a C-contiguous, two-dimensional
+threshold array with the pedestal's output dtype and shape. Noncontiguous
+inputs raise ``TypeError``; inputs with the wrong number of dimensions raise
+``ValueError``. Mismatched frame or threshold shapes raise ``RuntimeError``.
+Use ``numpy.ascontiguousarray()`` to copy a sliced or transposed array into the
+required layout when needed.
 
-* ``Pedestal_d`` uses ``float64`` storage
-* ``Pedestal_f`` uses ``float32`` storage
-* ``Pedestal_i16`` uses ``int16`` storage
+Internal sums and sums of squares always use ``float64``. Three
+specializations are available from :mod:`aare` for the mean and standard
+deviation output types:
+
+* ``Pedestal_d`` returns ``float64``
+* ``Pedestal_f`` returns ``float32``
+* ``Pedestal_i16`` returns ``int16``
+
+The public ``Pedestal`` factory selects the specialization from ``dtype``,
+defaulting to ``numpy.float64``.
+
+Constructor dimensions and ``n_samples`` must be positive integers. Negative
+values raise ``TypeError`` and zero values raise ``RuntimeError``.
+Internally, negative variance caused by floating-point roundoff is clamped to
+zero before taking its square root, keeping the standard deviation finite for
+nearly constant inputs. Only the final standard deviation is converted to the
+output dtype.
+
+Factory
+-------
+
+.. py:currentmodule:: aare
+
+.. autofunction:: Pedestal
 
 Example
 -------
 
 .. code-block:: python
 
-   from aare import Pedestal_d
+   import numpy as np
+   from aare import Pedestal
 
-   pedestal = Pedestal_d(512, 1024, 100)
+   pedestal = Pedestal(512, 1024, n_samples=100, dtype=np.float32)
 
    for frame in initialization_frames:
-       pedestal.push_no_update(frame)
+       pedestal.push(frame)
 
-   pedestal.update_mean()
    mean = pedestal.mean()
    noise = pedestal.std()
 

include/aare/FastPedestal.hpp

@@ -4,20 +4,20 @@
 #include "aare/Frame.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
+#include <algorithm>
 #include <cstddef>
 #include <sys/types.h>
 
 namespace aare {
 
 /**
- * @brief Maintain per-pixel mean, population variance, and standard deviation.
+ * @brief Maintain per-pixel mean and population 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.
+ * in double precision. Variance is a private double-precision intermediate.
+ * @tparam PEDESTAL_TYPE Type returned for the mean and standard deviation.
  */
 template <typename PEDESTAL_TYPE> class FastPedestal {
 
@@ -93,7 +93,7 @@ template <typename PEDESTAL_TYPE> class FastPedestal {
      * @brief Return a copy of the cached mean.
      * @throws std::runtime_error if ready() is false.
      */
-    NDArray<PEDESTAL_TYPE, 2> mean() {
+    NDArray<PEDESTAL_TYPE, 2> mean() const {
         if (!ready()) {
             throw std::runtime_error(
                 "Pedestal is not ready, cannot return mean");
@@ -127,58 +127,12 @@ template <typename PEDESTAL_TYPE> class FastPedestal {
      */
     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() {
+    NDArray<PEDESTAL_TYPE, 2> std() const {
         if (!ready()) {
             throw std::runtime_error(
                 "Pedestal is not ready, cannot return std");
@@ -206,7 +160,7 @@ template <typename PEDESTAL_TYPE> class FastPedestal {
                             "col={} must be in [0, {}), [0, {})",
                             row, col, m_rows, m_cols));
         }
-        return std::sqrt(variance(row, col));
+        return std_unchecked(rc_to_index(row, col));
     }
 
     /**
@@ -410,6 +364,14 @@ template <typename PEDESTAL_TYPE> class FastPedestal {
     uint32_t n_samples() const { return m_samples; }
 
   private:
+    double variance_unchecked(ssize_t index) const {
+        const auto &entry = m_sum[index];
+        const auto mean = entry.sum * m_inv_samples;
+        const auto variance = std::fma(-mean, mean, entry.sum2 * m_inv_samples);
+        // Roundoff in the moments can make a near-zero variance negative.
+        return std::max(variance, 0.0);
+    }
+
     /**
      * @brief Write the cached mean after the final add_init_frame. All other
      * (non initialization) pushes update the cached mean immediately.

include/aare/Pedestal.hpp

@@ -3,15 +3,19 @@
 #include "aare/Frame.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
+#include <algorithm>
 #include <cstddef>
 
 namespace aare {
 
 /**
- * @brief Calculate the pedestal of a series of frames. Can be used as
- * standalone but mostly used in the ClusterFinder.
+ * @brief Maintain per-pixel mean and population standard deviation.
  *
- * @tparam SUM_TYPE type of the sum
+ * Each pixel accumulates its first n_samples values. Subsequent pushes update
+ * the exponential moving average with a smoothing factor of 1 / n_samples.
+ * Statistics are available during initialization and are zero for empty pixels.
+ * Internal moments and the private variance intermediate use double precision.
+ * @tparam SUM_TYPE Type returned for the mean and standard deviation.
  */
 template <typename SUM_TYPE = double> class Pedestal {
     uint32_t m_rows;
@@ -20,98 +24,110 @@ template <typename SUM_TYPE = double> class Pedestal {
     uint32_t m_samples;
     NDArray<uint32_t, 2> m_cur_samples;
 
-    // TODO! in case of int needs to be changed to uint64_t
-    NDArray<SUM_TYPE, 2> m_sum;
-    NDArray<SUM_TYPE, 2> m_sum2;
+    NDArray<double, 2> m_sum;
+    NDArray<double, 2> m_sum2;
 
     // 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
     NDArray<SUM_TYPE, 2> m_mean;
 
-    // Cache std. Only refreshed via update_std() to keep push() cheap.
-    NDArray<SUM_TYPE, 2> m_std;
-
   public:
+    /**
+     * @brief Construct an empty pedestal with zero mean and standard deviation.
+     * @param rows Number of image rows.
+     * @param cols Number of image columns.
+     * @param n_samples Number of initialization samples per pixel and
+     * reciprocal of the weight assigned to each subsequent value.
+     * @throws std::runtime_error if rows, cols, or n_samples is zero.
+     */
     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_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;
+          m_sum(NDArray<double, 2>({rows, cols}, 0.0)),
+          m_sum2(NDArray<double, 2>({rows, cols}, 0.0)),
+          m_mean(NDArray<SUM_TYPE, 2>({rows, cols}, SUM_TYPE(0))) {
+        if (!(rows > 0 && cols > 0 && n_samples > 0)) {
+            throw std::runtime_error(
+                fmt::format("Invalid parameters for Pedestal: rows={}, "
+                            "cols={}, n_samples={} need to be positive",
+                            rows, cols, n_samples));
+        }
     }
+
     ~Pedestal() = default;
 
-    NDArray<SUM_TYPE, 2> mean() { return m_mean; }
-
+    /**
+     * @brief Return a non-owning view of the cached mean.
+     * @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 SUM_TYPE, 2> view() const { return m_mean.view(); }
 
+    /**
+     * @brief Return the cached mean at (row, col), or zero for an empty pixel.
+     * @pre row and col are valid pixel indices.
+     */
     SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
         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);
-    }
+    /** @brief Return a copy of the cached mean. */
+    NDArray<SUM_TYPE, 2> mean() const { return m_mean; }
 
+    /**
+     * @brief Calculate the population standard deviation at (row, col).
+     * @pre row and col are valid pixel indices.
+     * @note Variance is normalized by the pixel's current sample count. Empty
+     * pixels return zero, and negative variance from roundoff is clamped to
+     * zero.
+     */
     SUM_TYPE std(const uint32_t row, const uint32_t col) const {
         return std::sqrt(variance(row, col));
     }
 
-    SUM_TYPE variance(const uint32_t row, const uint32_t col) const {
-        if (m_cur_samples(row, col) == 0) {
-            return 0.0;
+    /**
+     * @brief Calculate and return the population standard deviation of every
+     * pixel. Empty pixels return zero.
+     */
+    NDArray<SUM_TYPE, 2> std() const {
+        NDArray<SUM_TYPE, 2> res({m_rows, m_cols});
+        for (uint32_t row = 0; row < m_rows; ++row) {
+            for (uint32_t col = 0; col < m_cols; ++col) {
+                res(row, col) = std(row, col);
+            }
         }
-        return m_sum2(row, col) / m_cur_samples(row, col) -
-               mean(row, col) * mean(row, col);
-    }
-
-    NDArray<SUM_TYPE, 2> variance() {
-        NDArray<SUM_TYPE, 2> variance_array({m_rows, m_cols});
-        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
-            variance_array(i / m_cols, i % m_cols) =
-                variance(i / m_cols, i % m_cols);
-        }
-        return variance_array;
-    }
-
-    NDArray<SUM_TYPE, 2> std() {
-        NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
-        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
-            standard_deviation_array(i / m_cols, i % m_cols) =
-                std(i / m_cols, i % m_cols);
-        }
-
-        return standard_deviation_array;
+        return res;
     }
 
+    /**
+     * @brief Zero the moments, cached mean, and sample counts of every pixel.
+     */
     void clear() {
         m_sum = 0;
         m_sum2 = 0;
         m_cur_samples = 0;
         m_mean = 0;
-        m_std = 0;
     }
 
+    /**
+     * @brief Zero the moments, cached mean, and sample count at (row, col).
+     * @pre row and col are valid pixel indices.
+     */
     void clear(const uint32_t row, const uint32_t col) {
         m_sum(row, col) = 0;
         m_sum2(row, col) = 0;
         m_cur_samples(row, col) = 0;
         m_mean(row, col) = 0;
-        m_std(row, col) = 0;
     }
 
+    /**
+     * @brief Accumulate or exponentially update every pixel and its cached
+     * mean.
+     * @param frame Frame whose shape must exactly match the pedestal.
+     * @throws std::runtime_error if the shape differs.
+     */
     template <typename T> void push(NDView<T, 2> frame) {
-        assert(frame.size() == m_rows * m_cols);
-
         // TODO! move away from m_rows, m_cols
         if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
             throw std::runtime_error(
@@ -125,11 +141,18 @@ template <typename SUM_TYPE = double> class Pedestal {
         }
     }
 
+    /**
+     * @brief Push only pixels whose absolute difference from the cached mean is
+     * strictly less than their threshold.
+     * @param frame Frame whose shape must exactly match the pedestal.
+     * @param threshold Per-pixel thresholds with the same shape as the
+     * pedestal.
+     * @throws std::runtime_error if either shape differs.
+     * @note Rejected pixels keep their statistics and sample counts unchanged.
+     */
     template <typename T>
     void push_with_threshold(const NDView<T, 2> frame,
                              const NDView<SUM_TYPE, 2> threshold) {
-        assert(frame.size() == m_rows * m_cols);
-
         // TODO! move away from m_rows, m_cols
         if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
             throw std::runtime_error(
@@ -152,45 +175,43 @@ template <typename SUM_TYPE = double> class Pedestal {
     }
 
     /**
-     * Push but don't update the cached mean. Speeds up the process
-     * when initializing the pedestal.
-     *
+     * @brief Accumulate or exponentially update every pixel from a Frame.
+     * @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.
      */
-    template <typename T> void push_no_update(NDView<T, 2> frame) {
-        assert(frame.size() == m_rows * m_cols);
+    template <typename T> void push(Frame &frame) { push<T>(frame.view<T>()); }
 
-        // TODO! move away from m_rows, m_cols
-        if (frame.shape() != std::array<ssize_t, 2>{m_rows, m_cols}) {
-            throw std::runtime_error(
-                "Frame shape does not match pedestal shape");
-        }
-
-        for (size_t row = 0; row < m_rows; row++) {
-            for (size_t col = 0; col < m_cols; col++) {
-                push_no_update<T>(row, col, frame(row, col));
-            }
-        }
-    }
-
-    template <typename T> void push(Frame &frame) {
-        assert(frame.rows() == static_cast<size_t>(m_rows) &&
-               frame.cols() == static_cast<size_t>(m_cols));
-        push<T>(frame.view<T>());
-    }
-
-    // getter functions
+    /** @brief Return the number of image rows. */
     uint32_t rows() const { return m_rows; }
-    uint32_t cols() const { return m_cols; }
-    uint32_t n_samples() const { return m_samples; }
-    NDArray<uint32_t, 2> cur_samples() const { return m_cur_samples; }
-    NDArray<SUM_TYPE, 2> get_sum() const { return m_sum; }
-    NDArray<SUM_TYPE, 2> get_sum2() const { return m_sum2; }
 
-    // pixel level operations (should be refactored to allow users to implement
-    // their own pixel level operations)
+    /** @brief Return the number of image columns. */
+    uint32_t cols() const { return m_cols; }
+
+    /**
+     * @brief Return the initialization sample count per pixel and steady-state
+     * update-weight denominator.
+     */
+    uint32_t n_samples() const { return m_samples; }
+
+    /**
+     * @brief Return a copy of the per-pixel initialization sample counts.
+     * @note Each count is in [0, n_samples] and does not change during
+     * steady-state pushes. Thresholded pushes can leave counts unequal.
+     */
+    NDArray<uint32_t, 2> cur_samples() const { return m_cur_samples; }
+
+    /**
+     * @brief Accumulate or exponentially update one pixel and its cached mean.
+     * @param row Pixel row.
+     * @param col Pixel column.
+     * @param val_ New pixel value, with weight 1 / n_samples after
+     * initialization.
+     * @pre row and col are valid pixel indices.
+     */
     template <typename T>
     void push(const uint32_t row, const uint32_t col, const T val_) {
-        SUM_TYPE val = static_cast<SUM_TYPE>(val_);
+        const auto val = static_cast<double>(val_);
         if (m_cur_samples(row, col) < m_samples) {
             m_sum(row, col) += val;
             m_sum2(row, col) += val * val;
@@ -204,45 +225,16 @@ template <typename SUM_TYPE = double> class Pedestal {
         m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col);
     }
 
-    template <typename T>
-    void push_no_update(const uint32_t row, const uint32_t col, const T val_) {
-        SUM_TYPE val = static_cast<SUM_TYPE>(val_);
-        if (m_cur_samples(row, col) < m_samples) {
-            m_sum(row, col) += val;
-            m_sum2(row, col) += val * val;
-            m_cur_samples(row, col)++;
-        } else {
-            m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col);
-            m_sum2(row, col) +=
-                val * val - m_sum2(row, col) / m_cur_samples(row, col);
+  private:
+    double variance(const uint32_t row, const uint32_t col) const {
+        if (m_cur_samples(row, col) == 0) {
+            return 0.0;
         }
-    }
-
-    /**
-     * @brief Update the mean of the pedestal. This is used after having done
-     * push_no_update. It is not necessary to call this function after push.
-     */
-    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
-        // m_samples samples
-        SUM_TYPE val = static_cast<SUM_TYPE>(val_);
-        m_sum(row, col) += val - m_sum(row, col) / m_samples;
-        m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples;
-        m_mean(row, col) = m_sum(row, col) / m_samples;
+        const auto mean = m_sum(row, col) / m_cur_samples(row, col);
+        const auto var =
+            m_sum2(row, col) / m_cur_samples(row, col) - mean * mean;
+        // Roundoff in the moments can make a near-zero variance negative.
+        return std::max(var, 0.0);
     }
 };
 } // namespace aare

python/CMakeLists.txt

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

python/aare/FastPedestal.py

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

python/aare/Pedestal.py

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

python/aare/__init__.py

@@ -33,6 +33,7 @@ from ._aare import UDPPortPosition
 
 from ._version import __version__
 from .FastPedestal import FastPedestal
+from .Pedestal import Pedestal
 from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile
 from .ClusterVector import ClusterVector
 from .Cluster import Cluster

python/src/bind_FastPedestal.hpp

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

python/src/bind_Pedestal.hpp

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

python/src/module.cpp

@@ -16,6 +16,7 @@
 #include "bind_FastPedestal.hpp"
 #include "bind_Interpolator.hpp"
 #include "bind_MultiThreadedFileReader.hpp"
+#include "bind_Pedestal.hpp"
 #include "bind_PedestalTrackingPixelHistogram.hpp"
 #include "bind_PixelHistogram.hpp"
 #include "bind_PixelMap.hpp"
@@ -27,7 +28,6 @@
 #include "file.hpp"
 #include "fit.hpp"
 #include "jungfrau_data_file.hpp"
-#include "pedestal.hpp"
 #include "raw_master_file.hpp"
 #include "raw_sub_file.hpp"
 #include "var_cluster.hpp"

python/src/pedestal.hpp

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

python/tests/test_FastPedestal.py

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

python/tests/test_Pedestal.py

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

src/FastPedestal.test.cpp

@@ -0,0 +1,60 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/FastPedestal.hpp"
+
+#include <catch2/catch_template_test_macros.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+using namespace aare;
+
+TEMPLATE_TEST_CASE("fast pedestal keeps variance in double precision",
+                   "[fast-pedestal]", double, float, int16_t) {
+    FastPedestal<TestType> pedestal(1, 1, 2);
+    NDArray<uint16_t, 2> frame({1, 1}, 0);
+    pedestal.add_init_frame(frame.view());
+    frame[0] = 1000;
+    pedestal.add_init_frame(frame.view());
+
+    const auto &ready = pedestal;
+    REQUIRE(ready.std(0, 0) == static_cast<TestType>(500));
+    REQUIRE(ready.std_unchecked(0) == ready.std(0, 0));
+    REQUIRE(ready.std()[0] == ready.std(0, 0));
+
+    pedestal.push_ema(frame.view());
+    REQUIRE(ready.std(0, 0) == static_cast<TestType>(std::sqrt(187500.0)));
+    REQUIRE(ready.std_unchecked(0) == ready.std(0, 0));
+    REQUIRE(ready.std()[0] == ready.std(0, 0));
+}
+
+TEST_CASE("fast pedestal validates standard deviation access",
+          "[fast-pedestal]") {
+    FastPedestal<double> pedestal(1, 1, 1);
+    REQUIRE_THROWS(pedestal.std());
+    REQUIRE_THROWS(pedestal.std(0, 0));
+
+    NDArray<uint16_t, 2> frame({1, 1}, 42);
+    pedestal.add_init_frame(frame.view());
+    REQUIRE(pedestal.std(0, 0) == 0.0);
+    REQUIRE_THROWS(pedestal.std(1, 0));
+    REQUIRE_THROWS(pedestal.std(0, 1));
+}
+
+TEMPLATE_TEST_CASE("fast pedestal std stays finite after settling",
+                   "[fast-pedestal]", double, float, int16_t) {
+    FastPedestal<TestType> pedestal(1, 1, 10);
+    NDArray<uint16_t, 2> frame({1, 1}, 16382);
+    pedestal.add_init_frame(frame.view());
+    frame[0] = 16383;
+    for (int i = 0; i < 9; ++i) {
+        pedestal.add_init_frame(frame.view());
+    }
+    for (int i = 0; i < 300; ++i) {
+        pedestal.push_ema(frame.view());
+    }
+
+    const auto noise = pedestal.std(0, 0);
+    REQUIRE(std::isfinite(noise));
+    REQUIRE(noise >= 0);
+    REQUIRE(static_cast<double>(noise) < 1e-3);
+    REQUIRE(pedestal.std_unchecked(0) == noise);
+    REQUIRE(pedestal.std()[0] == noise);
+}

src/Pedestal.test.cpp

@@ -1,12 +1,35 @@
 // SPDX-License-Identifier: MPL-2.0
 #include "aare/Pedestal.hpp"
 
+#include <catch2/catch_template_test_macros.hpp>
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
 #include <chrono>
 #include <random>
 
 using namespace aare;
+
+TEMPLATE_TEST_CASE("pedestal uses double precision moments", "[pedestal]",
+                   double, float, int16_t) {
+    Pedestal<TestType> pedestal(1, 1, 2);
+    pedestal.push(0, 0, uint16_t{30000});
+    pedestal.push(0, 0, uint16_t{30003});
+
+    REQUIRE(pedestal.mean(0, 0) == static_cast<TestType>(30001.5));
+    REQUIRE(pedestal.std(0, 0) == static_cast<TestType>(1.5));
+
+    pedestal.push(0, 0, uint16_t{30002});
+
+    REQUIRE(pedestal.mean(0, 0) == static_cast<TestType>(30001.75));
+    // Static cast to double is needed to avoid a Catch2 warning about
+    // float to double conversion on mac.
+    REQUIRE_THAT(
+        static_cast<double>(pedestal.std(0, 0)),
+        Catch::Matchers::WithinAbs(static_cast<double>(static_cast<TestType>(
+                                       std::sqrt(TestType(1.1875)))),
+                                   static_cast<double>(TestType(1e-6))));
+}
+
 TEST_CASE("test pedestal constructor") {
     aare::Pedestal pedestal(10, 10, 5);
     REQUIRE(pedestal.rows() == 10);
@@ -14,13 +37,47 @@ TEST_CASE("test pedestal constructor") {
     REQUIRE(pedestal.n_samples() == 5);
     for (int i = 0; i < 10; i++) {
         for (int j = 0; j < 10; j++) {
-            REQUIRE(pedestal.get_sum()(i, j) == 0);
-            REQUIRE(pedestal.get_sum2()(i, j) == 0);
+            REQUIRE(pedestal.mean(i, j) == 0);
             REQUIRE(pedestal.cur_samples()(i, j) == 0);
         }
     }
 }
 
+TEST_CASE("pedestal rejects mismatched frame shapes", "[pedestal]") {
+    Pedestal<> pedestal(2, 3);
+    pedestal.push(0, 0, uint16_t{7});
+    NDArray<double, 2> threshold({2, 3}, 10.0);
+
+    for (const auto shape :
+         {std::array<uint32_t, 2>{2, 2}, std::array<uint32_t, 2>{3, 2}}) {
+        Frame frame(shape[0], shape[1], Dtype::UINT16);
+        REQUIRE_THROWS_WITH(pedestal.push(frame.view<uint16_t>()),
+                            "Frame shape does not match pedestal shape");
+        REQUIRE_THROWS_WITH(pedestal.push_with_threshold(frame.view<uint16_t>(),
+                                                         threshold.view()),
+                            "Frame shape does not match pedestal shape");
+        REQUIRE_THROWS_WITH(pedestal.push<uint16_t>(frame),
+                            "Frame shape does not match pedestal shape");
+        REQUIRE(pedestal.mean(0, 0) == 7);
+        REQUIRE(pedestal.cur_samples()(0, 0) == 1);
+    }
+}
+
+TEMPLATE_TEST_CASE("pedestal std stays finite after settling", "[pedestal]",
+                   double, float, int16_t) {
+    Pedestal<TestType> pedestal(1, 1, 10);
+    pedestal.push(0, 0, uint16_t{16382});
+    for (int i = 0; i < 201; ++i) {
+        pedestal.push(0, 0, uint16_t{16383});
+    }
+
+    const auto noise = pedestal.std(0, 0);
+    REQUIRE(std::isfinite(noise));
+    REQUIRE(noise >= 0);
+    REQUIRE(static_cast<double>(noise) < 1e-3);
+    REQUIRE(pedestal.std()(0, 0) == noise);
+}
+
 TEST_CASE("test pedestal push") {
     aare::Pedestal pedestal(10, 10, 5);
     aare::Frame frame(10, 10, Dtype::UINT16);
@@ -34,8 +91,7 @@ TEST_CASE("test pedestal push") {
     pedestal.push<uint16_t>(frame);
     for (int i = 0; i < 10; i++) {
         for (int j = 0; j < 10; j++) {
-            REQUIRE(pedestal.get_sum()(i, j) == i + j);
-            REQUIRE(pedestal.get_sum2()(i, j) == (i + j) * (i + j));
+            REQUIRE(pedestal.mean(i, j) == i + j);
             REQUIRE(pedestal.cur_samples()(i, j) == 1);
         }
     }
@@ -44,8 +100,7 @@ TEST_CASE("test pedestal push") {
     pedestal.clear();
     for (int i = 0; i < 10; i++) {
         for (int j = 0; j < 10; j++) {
-            REQUIRE(pedestal.get_sum()(i, j) == 0);
-            REQUIRE(pedestal.get_sum2()(i, j) == 0);
+            REQUIRE(pedestal.mean(i, j) == 0);
             REQUIRE(pedestal.cur_samples()(i, j) == 0);
         }
     }
@@ -57,16 +112,10 @@ TEST_CASE("test pedestal push") {
             for (uint32_t j = 0; j < 10; j++) {
                 if (k < 5) {
                     REQUIRE(pedestal.cur_samples()(i, j) == k + 1);
-                    REQUIRE(pedestal.get_sum()(i, j) == (k + 1) * (i + j));
-                    REQUIRE(pedestal.get_sum2()(i, j) ==
-                            (k + 1) * (i + j) * (i + j));
                 } else {
                     REQUIRE(pedestal.cur_samples()(i, j) == 5);
-                    REQUIRE(pedestal.get_sum()(i, j) == 5 * (i + j));
-                    REQUIRE(pedestal.get_sum2()(i, j) == 5 * (i + j) * (i + j));
                 }
                 REQUIRE(pedestal.mean(i, j) == (i + j));
-                REQUIRE(pedestal.variance(i, j) == 0);
                 REQUIRE(pedestal.std(i, j) == 0);
             }
         }
@@ -74,7 +123,7 @@ TEST_CASE("test pedestal push") {
 }
 
 TEST_CASE("test pedestal with normal distribution") {
-    const double MEAN = 5.0, STD = 2.0, VAR = STD * STD, TOLERANCE = 0.1;
+    const double MEAN = 5.0, STD = 2.0, TOLERANCE = 0.1;
 
     unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
     std::default_random_engine generator(seed);
@@ -91,17 +140,14 @@ TEST_CASE("test pedestal with normal distribution") {
         pedestal.push<double>(frame);
     }
     auto mean = pedestal.mean();
-    auto variance = pedestal.variance();
     auto standard_deviation = pedestal.std();
 
     for (int i = 0; i < 3; i++) {
         for (int j = 0; j < 5; j++) {
             REQUIRE_THAT(mean(i, j),
                          Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE));
-            REQUIRE_THAT(variance(i, j),
-                         Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE));
             REQUIRE_THAT(standard_deviation(i, j),
                          Catch::Matchers::WithinAbs(STD, STD * TOLERANCE));
         }
     }
-}
+}