Merge 'origin/main' into feature/cuda_clusterfinder

python/src/bind_Defs.hpp

@@ -22,4 +22,9 @@ void define_defs_bindings(py::module &m) {
     moench04.attr("nPixelsPerSuperColumn") = Moench04::nPixelsPerSuperColumn;
     moench04.attr("superColumnWidth") = Moench04::superColumnWidth;
     moench04.attr("adcNumbers") = Moench04::adcNumbers;
+
+    auto moench05 = py::class_<Moench05>(m, "Moench05");
+    moench05.attr("nRows") = Moench05::nRows;
+    moench05.attr("nCols") = Moench05::nCols;
+    moench05.attr("adcNumbers") = Moench05::adcNumbers;
 }

python/src/bind_Eta.hpp

@@ -13,12 +13,12 @@ void define_eta(py::module &m, const std::string &typestr) {
 
     py::class_<Eta2<T>>(m, class_name.c_str())
         .def(py::init<>())
-        .def_readonly("x", &Eta2<T>::x, "eta x value")
-        .def_readonly("y", &Eta2<T>::y, "eta y value")
-        .def_readonly("c", &Eta2<T>::c,
-                      "eta corner value cTopLeft, cTopRight, "
-                      "cBottomLeft, cBottomRight")
-        .def_readonly("sum", &Eta2<T>::sum, "photon energy of cluster");
+        .def_readwrite("x", &Eta2<T>::x, "eta x value")
+        .def_readwrite("y", &Eta2<T>::y, "eta y value")
+        .def_readwrite("c", &Eta2<T>::c,
+                       "eta corner value cTopLeft, cTopRight, "
+                       "cBottomLeft, cBottomRight")
+        .def_readwrite("sum", &Eta2<T>::sum, "photon energy of cluster");
 }
 
 void define_corner_enum(py::module &m) {

python/src/bind_Interpolator.hpp

@@ -11,17 +11,20 @@
 
 namespace py = pybind11;
 
+// clang-format off
 #define REGISTER_INTERPOLATOR_ETA2(T, N, M, U)                                 \
     register_interpolate<T, N, M, U, aare::calculate_full_eta2<T, N, M, U>>(   \
         interpolator, "_full_eta2", "full eta2");                              \
     register_interpolate<T, N, M, U, aare::calculate_eta2<T, N, M, U>>(        \
-        interpolator, "", "eta2");
+        interpolator, "", "eta2");                                             \
+    register_interpolate_custom_eta<T, N, M, U>(interpolator);
 
 #define REGISTER_INTERPOLATOR_ETA3(T, N, M, U)                                 \
     register_interpolate<T, N, M, U, aare::calculate_eta3<T, N, M, U>>(        \
         interpolator, "_eta3", "full eta3");                                   \
     register_interpolate<T, N, M, U, aare::calculate_cross_eta3<T, N, M, U>>(  \
         interpolator, "_cross_eta3", "cross eta3");
+// clang-format on 
 
 template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
           typename CoordType = uint16_t, auto EtaFunction>
@@ -48,6 +51,34 @@ void register_interpolate(py::class_<aare::Interpolator> &interpolator,
         docstring.c_str(), py::arg("cluster_vector"));
 }
 
+template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void register_interpolate_custom_eta(
+    py::class_<aare::Interpolator> &interpolator) {
+
+    using ClusterType = Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>;
+
+    interpolator.def(
+        "interpolate",
+        [](aare::Interpolator &self, const ClusterVector<ClusterType> &clusters,
+           const std::vector<Eta2<Type>> &etas) {
+            auto photons = self.interpolate<Type, ClusterSizeX, ClusterSizeY, CoordType>(clusters, etas);
+            auto *ptr = new std::vector<Photon>{std::move(photons)};
+            return return_vector(ptr);
+        },
+        R"(
+        Interpolation based on custom eta values provided by the user. 
+
+        Args:
+        cluster_vector: vector of clusters to interpolate
+        etas: vector of eta values for each cluster (must be in the same order as the clusters)
+
+        Returns:
+        interpolated photons
+        )",
+        py::arg("cluster_vector"), py::arg("etas"));
+}
+
 template <typename Type>
 void register_transform_eta_values(
     py::class_<aare::Interpolator> &interpolator) {
@@ -65,6 +96,8 @@ void define_interpolation_bindings(py::module &m) {
 
     PYBIND11_NUMPY_DTYPE(aare::Photon, x, y, energy);
 
+    PYBIND11_NUMPY_DTYPE(aare::Coordinate2D, x, y);
+
     auto interpolator =
         py::class_<aare::Interpolator>(m, "Interpolator")
             .def(py::init(

python/src/bind_PedestalTrackingPixelHistogram.hpp

@@ -0,0 +1,245 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/hist/PedestalTrackingPixelHistogram.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) {
+    py::class_<PedestalTrackingPixelHistogram>(
+        m, "PedestalTrackingPixelHistogram",
+        "A pixel-wise histogram of frame - pedestal residuals, with a "
+        "per-pixel running pedestal estimate sharded across worker threads")
+        .def(
+            py::init<int, int, int, double, double, int, std::size_t, double>(),
+            R"(
+             Initialize a PedestalTrackingPixelHistogram.
+
+             Args:
+                 rows: Number of rows in the detector
+                 cols: Number of columns in the detector
+                 n_bins: Number of histogram bins along the residual axis
+                 xmin: Minimum residual value (inclusive)
+                 xmax: Maximum residual value (exclusive)
+                 n_threads: Number of worker threads (default: 1). Each
+                     worker owns a disjoint row slice of both the
+                     pedestal and the histogram, so the partition
+                     determines per-thread memory usage.
+                 max_pending: Maximum number of frames that can be
+                     queued for asynchronous filling before
+                     fill_async() applies backpressure
+                     on the caller (default: 16).
+                 n_sigma: Sigma multiplier used as the gate for the
+                     pedestal-update side effect of
+                     fill_async(): a pixel sample is
+                     pushed back into the pedestal estimate iff
+                     ``abs(residual) < n_sigma * cached_std``. Set to
+                     ``0.0`` to disable the pedestal update and get
+                     histogram-only async behaviour (default: 1.0).
+                     Also exposed live via the ``n_sigma`` property.
+             )",
+            py::kw_only(), py::arg("rows"), py::arg("cols"), py::arg("n_bins"),
+            py::arg("xmin"), py::arg("xmax"), py::arg("n_threads") = 1,
+            py::arg("max_pending") = std::size_t{16}, py::arg("n_sigma") = 1.0)
+
+        .def(
+            "push_pedestal_no_update",
+            [](PedestalTrackingPixelHistogram &self,
+               py::array_t<PedestalTrackingPixelHistogram::FrameType, 0>
+                   frame) {
+                auto view = make_view_2d(frame);
+                self.push_pedestal_no_update(view);
+            },
+            R"(
+             Accumulate `frame` into the per-pixel running pedestal
+             estimate without refreshing the cached mean.
+
+             Use repeatedly while bootstrapping the pedestal, then call
+             update_mean() once before starting to fill the histogram.
+
+             Args:
+                 frame: A 2D numpy array of raw pixel values (dtype: uint16)
+             )",
+            py::arg("frame").noconvert())
+
+        .def("update_mean", &PedestalTrackingPixelHistogram::update_mean,
+             R"(
+             Refresh each partial pedestal's cached per-pixel mean from
+             its running sums. Drains pending async fills first, then
+             dispatches the update to the worker pool so the writes to
+             each shard happen on the same thread that reads them in
+             fill_async().
+             )",
+             py::call_guard<py::gil_scoped_release>())
+
+        .def(
+            "pedestal_mean",
+            [](const PedestalTrackingPixelHistogram &self) {
+                // pedestal_mean() flushes + locks + memcpys; do all of
+                // that without the GIL, only reacquire to wrap into a
+                // numpy array.
+                NDArray<PedestalTrackingPixelHistogram::AxisType, 2> *ptr =
+                    nullptr;
+                {
+                    py::gil_scoped_release release;
+                    ptr = new NDArray<PedestalTrackingPixelHistogram::AxisType,
+                                      2>(self.pedestal_mean());
+                }
+                return return_image_data(ptr);
+            },
+            R"(
+             Snapshot the per-pixel pedestal mean stitched together
+             from all shards.
+
+             Returns:
+                 A 2D numpy array (rows x cols, dtype: float64)
+                 containing the current cached pedestal mean.
+             )")
+
+        .def(
+            "fill_async",
+            [](PedestalTrackingPixelHistogram &self,
+               py::array_t<PedestalTrackingPixelHistogram::FrameType, 0>
+                   image) {
+                // Copy the numpy buffer into an owned NDArray while we
+                // still hold the GIL so we don't depend on the array's
+                // backing storage outliving this call.
+                auto view = make_view_2d(image);
+                NDArray<PedestalTrackingPixelHistogram::FrameType, 2> owned(
+                    view);
+                // Release the GIL while enqueueing -
+                // fill_async can block on backpressure
+                // when the queue is full.
+                py::gil_scoped_release release;
+                self.fill_async(std::move(owned));
+            },
+            R"(
+             Submit an image for asynchronous filling with sigma-clipped
+             pedestal tracking.
+
+             For each pixel the worker pool:
+               * histograms the pedestal-subtracted residual when it
+                 falls in ``[xmin, xmax)``, and
+               * additionally pushes the raw pixel value back into the
+                 per-thread pedestal estimate when
+                 ``abs(residual) < n_sigma * cached_std`` (the
+                 sigma-clipped pedestal-update gate).
+
+             The cached std is populated by ``update_mean()``, so
+             ``push_pedestal_no_update()`` + ``update_mean()`` must have
+             run at least once for the pedestal-update side effect to
+             fire. Setting ``n_sigma = 0`` disables the side effect and
+             recovers plain histogram-only async filling.
+
+             The image is copied into an internal buffer before this call
+             returns, so the caller may mutate or free the numpy array
+             immediately. If the internal queue is full this call blocks
+             (with the GIL released) until a slot becomes available.
+
+             Args:
+                 image: A 2D numpy array of raw pixel values (dtype: uint16)
+             )",
+            py::arg("image").noconvert())
+
+        .def("fill_from_file", &PedestalTrackingPixelHistogram::fill_from_file,
+             R"(
+             Fill the histogram from a file.
+
+             Args:
+                 file_path: Path to the file to fill from
+                 max_frames: Maximum number of frames to fill from the file (default: -1)
+             )",
+             py::call_guard<py::gil_scoped_release>(), py::arg("fname"),
+             py::arg("max_frames") = -1, py::arg("verbose") = false)
+        .def("process_pedestal_file",
+             &PedestalTrackingPixelHistogram::process_pedestal_file,
+             R"(
+             Process a pedestal file.
+
+             Args:
+                 file_path: Path to the file to process
+                 max_frames: Maximum number of frames to process from the file (default: -1)
+             )",
+             py::call_guard<py::gil_scoped_release>(), py::arg("fname"),
+             py::arg("max_frames") = -1, py::arg("verbose") = false)
+        .def_property("n_sigma", &PedestalTrackingPixelHistogram::n_sigma,
+                      &PedestalTrackingPixelHistogram::set_n_sigma,
+                      R"(
+             Sigma multiplier used as the pedestal-update gate in
+             fill_async(). Atomic; safe to read or write
+             from any thread. Setting it to 0.0 disables the pedestal
+             update entirely. The new value takes effect on subsequent
+             per-pixel evaluations inside the worker pool.
+             )")
+
+        .def("flush", &PedestalTrackingPixelHistogram::flush,
+             R"(
+             Block until all images submitted via
+             fill_async() have been merged into the
+             accumulators. Cheap when nothing is pending.
+             )",
+             py::call_guard<py::gil_scoped_release>())
+
+        .def(
+            "values",
+            [](const PedestalTrackingPixelHistogram &self) {
+                // values() implicitly flushes - release the GIL while it
+                // does so. Allocation/copy into the NDArray runs without
+                // the GIL too; only the numpy wrapping needs it.
+                NDArray<PedestalTrackingPixelHistogram::StorageType, 3> *ptr =
+                    nullptr;
+                {
+                    py::gil_scoped_release release;
+                    ptr =
+                        new NDArray<PedestalTrackingPixelHistogram::StorageType,
+                                    3>(self.values());
+                }
+                return return_image_data(ptr);
+            },
+            R"(
+             Get the histogram data as a numpy array.
+
+             Implicitly flushes any pending asynchronous fills before
+             returning, so the snapshot is consistent with everything
+             submitted up to this call.
+
+             Returns:
+                 A 3D numpy array (rows x cols x n_bins, dtype: uint16)
+                 containing the histogram bins for each pixel.
+             )")
+
+        .def(
+            "bin_centers",
+            [](const PedestalTrackingPixelHistogram &self) {
+                auto ptr =
+                    new NDArray<PedestalTrackingPixelHistogram::AxisType, 1>(
+                        self.bin_centers());
+                return return_image_data(ptr);
+            },
+            R"(
+             Get the bin centers along the residual axis.
+
+             Returns:
+                 A 1D numpy array (dtype: float32) of bin center values.
+             )")
+
+        .def(
+            "bin_edges",
+            [](const PedestalTrackingPixelHistogram &self) {
+                auto ptr =
+                    new NDArray<PedestalTrackingPixelHistogram::AxisType, 1>(
+                        self.bin_edges());
+                return return_image_data(ptr);
+            },
+            R"(
+             Get the bin edges along the residual axis.
+
+             Returns:
+                 A 1D numpy array (dtype: float32) of bin edge values.
+             )");
+}

python/src/bind_PixelHistogram.hpp

@@ -0,0 +1,147 @@
+// SPDX-License-Identifier: MPL-2.0
+#include "aare/hist/PixelHistogram.hpp"
+#include "np_helper.hpp"
+
+#include <cstddef>
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <string>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+namespace {
+
+template <typename StorageType>
+void define_pixel_histogram_binding(py::module &m, const char *class_name,
+                                    const char *storage_dtype) {
+    using Hist = PixelHistogram<StorageType, double>;
+
+    const std::string doc =
+        std::string("A histogram for pixel-wise statistics with float64 input "
+                    "axis and ") +
+        storage_dtype + " bin storage";
+
+    py::class_<Hist>(m, class_name, doc.c_str())
+        .def(py::init<int, int, int, double, double, int, std::size_t>(),
+             R"(
+             Initialize a PixelHistogram.
+
+             Args:
+                 rows: Number of rows in the detector
+                 cols: Number of columns in the detector
+                 n_bins: Number of histogram bins
+                 xmin: Minimum value for histogram range
+                 xmax: Maximum value for histogram range
+                 n_threads: Number of threads for parallel filling (default: 1)
+                 max_pending: Maximum number of images that can be queued for
+                     asynchronous filling before fill_async() applies
+                     backpressure on the caller (default: 16)
+             )",
+             py::kw_only(), py::arg("rows"), py::arg("cols"), py::arg("n_bins"),
+             py::arg("xmin"), py::arg("xmax"), py::arg("n_threads") = 1,
+             py::arg("max_pending") = std::size_t{16})
+
+        .def(
+            "fill_async",
+            [](Hist &self, py::array_t<double, 0> image) {
+                // Copy the numpy buffer into an owned NDArray while we
+                // still hold the GIL so we don't depend on the array's
+                // backing storage outliving this call.
+                auto view = make_view_2d(image);
+                NDArray<double, 2> owned(view);
+                // Release the GIL while enqueueing - fill_async can block
+                // on backpressure when the queue is full.
+                py::gil_scoped_release release;
+                self.fill_async(std::move(owned));
+            },
+            R"(
+             Submit an image for asynchronous filling.
+
+             The image is copied into an internal buffer before this call
+             returns, so the caller may mutate or free the numpy array
+             immediately. The actual histogram update happens on a
+             background thread. If the internal queue is full this call
+             blocks (with the GIL released) until a slot becomes available.
+
+             Args:
+                 image: A 2D numpy array of pixel values (dtype: float64)
+             )",
+            py::arg("image").noconvert())
+
+        .def("flush", &Hist::flush,
+             R"(
+             Block until all images submitted via fill_async() have been
+             merged into the accumulators. Cheap when nothing is pending.
+             )",
+             py::call_guard<py::gil_scoped_release>())
+
+        .def(
+            "values",
+            [](const Hist &self) {
+                // values() implicitly flushes - release the GIL while it
+                // does so. Allocation/copy into the NDArray runs without
+                // the GIL too; only the numpy wrapping needs it.
+                NDArray<StorageType, 3> *ptr = nullptr;
+                {
+                    py::gil_scoped_release release;
+                    ptr = new NDArray<StorageType, 3>(self.values());
+                }
+                return return_image_data(ptr);
+            },
+            R"(
+             Get the histogram data as a numpy array.
+
+             Implicitly flushes any pending asynchronous fills before
+             returning, so the snapshot is consistent with everything
+             submitted up to this call.
+
+             Returns:
+                 A 3D numpy array containing the histogram bins for each pixel
+             )")
+
+        .def(
+            "bin_centers",
+            [](const Hist &self) {
+                auto ptr = new NDArray<double, 1>(self.bin_centers());
+                return return_image_data(ptr);
+            },
+            R"(
+             Get the bin centers along the value axis.
+
+             Returns:
+                 A 1D numpy array containing the center values for each histogram bin
+             )")
+        .def(
+            "bin_edges",
+            [](const Hist &self) {
+                auto ptr = new NDArray<double, 1>(self.bin_edges());
+                return return_image_data(ptr);
+            },
+            R"(
+             Get the bin edges along the value axis.
+
+             Returns:
+                 A 1D numpy array containing the edge values for the histogram bins
+             )");
+}
+
+} // namespace
+
+void define_pixel_histogram_bindings(py::module &m) {
+    define_pixel_histogram_binding<double>(m, "PixelHistogram_d", "float64");
+    define_pixel_histogram_binding<float>(m, "PixelHistogram_f", "float32");
+    define_pixel_histogram_binding<std::uint64_t>(m, "PixelHistogram_u64",
+                                                  "uint64");
+    define_pixel_histogram_binding<std::uint32_t>(m, "PixelHistogram_u32",
+                                                  "uint32");
+    define_pixel_histogram_binding<std::uint16_t>(m, "PixelHistogram_u16",
+                                                  "uint16");
+    define_pixel_histogram_binding<std::uint8_t>(m, "PixelHistogram_u8",
+                                                 "uint8");
+
+    // Backwards-compatible alias for the generic Python class name.
+    m.attr("PixelHistogram") = m.attr("PixelHistogram_d");
+}

python/src/bind_RawFile.hpp

@@ -280,6 +280,7 @@ void define_raw_file_io_bindings(py::module &m) {
         .def("tell", &RawFile::tell, R"(
             Return the current frame number.)")
         .def_property_readonly("total_frames", &RawFile::total_frames)
+        .def("__len__", &RawFile::total_frames)
         .def("rows", static_cast<size_t (RawFile::*)() const>(&RawFile::rows))
         .def(
             "rows",

python/src/ctb_raw_file.hpp

@@ -226,5 +226,7 @@ void define_ctb_raw_file_io_bindings(py::module &m) {
         .def_property_readonly("image_size_in_bytes",
                                &CtbRawFile::image_size_in_bytes)
 
-        .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
+        .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file)
+        .def_property_readonly("total_frames", &CtbRawFile::total_frames)
+        .def("__len__", &CtbRawFile::total_frames);
 }

python/src/file.hpp

@@ -60,6 +60,7 @@ void define_file_io_bindings(py::module &m) {
         .def("seek", &File::seek)
         .def("tell", &File::tell)
         .def_property_readonly("total_frames", &File::total_frames)
+        .def("__len__", &File::total_frames)
         .def_property_readonly("rows", &File::rows)
         .def_property_readonly("cols", &File::cols)
         .def_property_readonly("bitdepth", &File::bitdepth)

python/src/fit.hpp

@@ -5,7 +5,6 @@
 #include <pybind11/stl.h>
 #include <pybind11/stl_bind.h>
 
-#include "aare/Chi2.hpp"
 #include "aare/Fit.hpp"
 #include "aare/FitModel.hpp"
 #include "aare/Models.hpp"
@@ -13,7 +12,7 @@
 namespace py = pybind11;
 using namespace pybind11::literals;
 
-template <typename Model, typename FCN>
+template <typename Model>
 py::object
 fit_dispatch(const aare::FitModel<Model> &model,
              py::array_t<double, py::array::c_style | py::array::forcecast> x,
@@ -22,7 +21,6 @@ fit_dispatch(const aare::FitModel<Model> &model,
 
 template <typename Model> void bind_fit_model(py::module &m, const char *name) {
     using FM = aare::FitModel<Model>;
-    using FCN = aare::func::Chi2Model1DGrad<Model>;
     py::class_<FM>(m, name)
         .def(py::init<unsigned int, unsigned int, double, bool>(),
              py::arg("strategy") = 0, py::arg("max_calls") = 100,
@@ -85,8 +83,7 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
                py::array_t<double, py::array::c_style | py::array::forcecast> x,
                py::array_t<double, py::array::c_style | py::array::forcecast> y,
                py::object y_err_obj, int n_threads) -> py::object {
-                return fit_dispatch<Model, FCN>(self, x, y, y_err_obj,
-                                                n_threads);
+                return fit_dispatch<Model>(self, x, y, y_err_obj, n_threads);
             },
             R"doc(
             Fit this model to 1D or 3D data using Minuit2.
@@ -145,7 +142,7 @@ py::dict pack_1d_result_dict(const aare::NDArray<double, 1> &result,
 }
 
 // Helper: typed dispatch for one Model, handles 1D/3D + y_err logic
-template <typename Model, typename FCN>
+template <typename Model>
 py::object
 fit_dispatch(const aare::FitModel<Model> &model,
              py::array_t<double, py::array::c_style | py::array::forcecast> x,
@@ -175,9 +172,9 @@ fit_dispatch(const aare::FitModel<Model> &model,
                 new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
             auto y_view_err = make_view_3d(y_err);
 
-            aare::fit_3d<Model, FCN>(model, x_view, y_view, y_view_err,
-                                     par_out->view(), err_out->view(),
-                                     chi2_out->view(), n_threads);
+            aare::fit_3d<Model>(model, x_view, y_view, y_view_err,
+                                par_out->view(), err_out->view(),
+                                chi2_out->view(), n_threads);
 
             if (model.compute_errors()) {
                 return py::dict("par"_a = return_image_data(par_out),
@@ -193,9 +190,9 @@ fit_dispatch(const aare::FitModel<Model> &model,
             NDView<double, 3> dummy_err{};
             NDView<double, 3> dummy_err_out{};
 
-            aare::fit_3d<Model, FCN>(model, x_view, y_view, dummy_err,
-                                     par_out->view(), dummy_err_out,
-                                     chi2_out->view(), n_threads);
+            aare::fit_3d<Model>(model, x_view, y_view, dummy_err,
+                                par_out->view(), dummy_err_out,
+                                chi2_out->view(), n_threads);
 
             return py::dict("par"_a = return_image_data(par_out),
                             "chi2"_a = return_image_data(chi2_out));
@@ -217,10 +214,9 @@ fit_dispatch(const aare::FitModel<Model> &model,
             }
 
             auto y_view_err = make_view_1d(y_err);
-            result =
-                aare::fit_pixel<Model, FCN>(model, x_view, y_view, y_view_err);
+            result = aare::fit_pixel<Model>(model, x_view, y_view, y_view_err);
         } else {
-            result = aare::fit_pixel<Model, FCN>(model, x_view, y_view);
+            result = aare::fit_pixel<Model>(model, x_view, y_view);
         }
 
         return pack_1d_result_dict<Model>(result, model.compute_errors());
@@ -679,6 +675,11 @@ void define_fit_bindings(py::module &m) {
 
     // ── Bind model classes ──────────────────────────────────────────
     bind_fit_model<aare::model::Gaussian>(m, "Gaussian");
+    bind_fit_model<aare::model::GaussianErfcPlateau>(m, "GaussianErfcPlateau");
+    bind_fit_model<aare::model::GaussianChargeSharing>(m,
+                                                       "GaussianChargeSharing");
+    bind_fit_model<aare::model::GaussianChargeSharingKb>(
+        m, "GaussianChargeSharingKb");
     bind_fit_model<aare::model::RisingScurve>(m, "RisingScurve");
     bind_fit_model<aare::model::FallingScurve>(m, "FallingScurve");
     bind_fit_model<aare::model::Pol1>(m, "Pol1");
@@ -691,28 +692,54 @@ void define_fit_bindings(py::module &m) {
            py::array_t<double, py::array::c_style | py::array::forcecast> y,
            py::object y_err_obj, int n_threads) -> py::object {
             using namespace aare::model;
-            using namespace aare::func;
 
             // ── Polynomial of degree 1 ───────
             if (py::isinstance<aare::FitModel<Pol1>>(model_obj)) {
                 const auto &mdl =
                     model_obj.cast<const aare::FitModel<Pol1> &>();
-                return fit_dispatch<Pol1, Chi2Pol1>(mdl, x, y, y_err_obj,
-                                                    n_threads);
+                return fit_dispatch<Pol1>(mdl, x, y, y_err_obj, n_threads);
             }
 
             // ── Polynomial of degree 2 ───────
             if (py::isinstance<aare::FitModel<Pol2>>(model_obj)) {
                 const auto &mdl =
                     model_obj.cast<const aare::FitModel<Pol2> &>();
-                return fit_dispatch<Pol2, Chi2Pol2>(mdl, x, y, y_err_obj,
-                                                    n_threads);
+                return fit_dispatch<Pol2>(mdl, x, y, y_err_obj, n_threads);
             }
+
             // ── Gaussian ───────
             if (py::isinstance<aare::FitModel<Gaussian>>(model_obj)) {
                 const auto &mdl =
                     model_obj.cast<const aare::FitModel<Gaussian> &>();
-                return fit_dispatch<Gaussian, Chi2Gaussian>(
+                return fit_dispatch<Gaussian>(mdl, x, y, y_err_obj, n_threads);
+            }
+
+            // ── GaussianErfcPlateau ───────
+            if (py::isinstance<aare::FitModel<GaussianErfcPlateau>>(
+                    model_obj)) {
+                const auto &mdl =
+                    model_obj
+                        .cast<const aare::FitModel<GaussianErfcPlateau> &>();
+                return fit_dispatch<GaussianErfcPlateau>(mdl, x, y, y_err_obj,
+                                                         n_threads);
+            }
+
+            // ── GaussianChargeSharing ───────
+            if (py::isinstance<aare::FitModel<GaussianChargeSharing>>(
+                    model_obj)) {
+                const auto &mdl =
+                    model_obj
+                        .cast<const aare::FitModel<GaussianChargeSharing> &>();
+                return fit_dispatch<GaussianChargeSharing>(mdl, x, y, y_err_obj,
+                                                           n_threads);
+            }
+
+            // ── GaussianChargeSharingKb ───────
+            if (py::isinstance<aare::FitModel<GaussianChargeSharingKb>>(
+                    model_obj)) {
+                const auto &mdl = model_obj.cast<
+                    const aare::FitModel<GaussianChargeSharingKb> &>();
+                return fit_dispatch<GaussianChargeSharingKb>(
                     mdl, x, y, y_err_obj, n_threads);
             }
 
@@ -720,16 +747,16 @@ void define_fit_bindings(py::module &m) {
             if (py::isinstance<aare::FitModel<RisingScurve>>(model_obj)) {
                 const auto &mdl =
                     model_obj.cast<const aare::FitModel<RisingScurve> &>();
-                return fit_dispatch<RisingScurve, Chi2RisingScurve>(
-                    mdl, x, y, y_err_obj, n_threads);
+                return fit_dispatch<RisingScurve>(mdl, x, y, y_err_obj,
+                                                  n_threads);
             }
 
             // ── Falling Scurve ───────
             if (py::isinstance<aare::FitModel<FallingScurve>>(model_obj)) {
                 const auto &mdl =
                     model_obj.cast<const aare::FitModel<FallingScurve> &>();
-                return fit_dispatch<FallingScurve, Chi2FallingScurve>(
-                    mdl, x, y, y_err_obj, n_threads);
+                return fit_dispatch<FallingScurve>(mdl, x, y, y_err_obj,
+                                                   n_threads);
             }
 
             throw std::runtime_error(

python/src/jungfrau_data_file.hpp

@@ -72,6 +72,7 @@ void define_jungfrau_data_file_io_bindings(py::module &m) {
         .def_property_readonly("bitdepth", &JungfrauDataFile::bitdepth)
         .def_property_readonly("current_file", &JungfrauDataFile::current_file)
         .def_property_readonly("total_frames", &JungfrauDataFile::total_frames)
+        .def("__len__", &JungfrauDataFile::total_frames)
         .def_property_readonly("n_files", &JungfrauDataFile::n_files)
         .def("read_frame", &read_dat_frame,
              R"(

python/src/module.cpp

@@ -12,6 +12,8 @@
 #include "bind_Defs.hpp"
 #include "bind_Eta.hpp"
 #include "bind_Interpolator.hpp"
+#include "bind_PedestalTrackingPixelHistogram.hpp"
+#include "bind_PixelHistogram.hpp"
 #include "bind_PixelMap.hpp"
 #include "bind_RawFile.hpp"
 #include "bind_calibration.hpp"
@@ -64,6 +66,8 @@ PYBIND11_MODULE(_aare, m) {
     define_raw_master_file_bindings(m);
     define_var_cluster_finder_bindings(m);
     define_pixel_map_bindings(m);
+    define_pixel_histogram_bindings(m);
+    define_pedestal_tracking_pixel_histogram_bindings(m);
     define_pedestal_bindings<double>(m, "Pedestal_d");
     define_pedestal_bindings<float>(m, "Pedestal_f");
     define_fit_bindings(m);

python/src/pedestal.hpp

@@ -5,6 +5,7 @@
 
 #include <cstdint>
 #include <filesystem>
+#include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 
@@ -12,7 +13,8 @@ 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::class_<Pedestal<SUM_TYPE>>(m, name.c_str(), py::buffer_protocol())
         .def(py::init<int, int, int>())
         .def(py::init<int, int>())
         .def("mean",
@@ -21,6 +23,23 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
                  *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>{};
@@ -33,6 +52,23 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
                  *std = self.std();
                  return return_image_data(std);
              })
+        .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)
@@ -49,6 +85,16 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
                  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,
@@ -57,5 +103,17 @@ void define_pedestal_bindings(py::module &m, const std::string &name) {
                 pedestal.push_no_update(v);
             },
             py::arg().noconvert())
-        .def("update_mean", &Pedestal<SUM_TYPE>::update_mean);
-}
+        .def("update_mean", &Pedestal<SUM_TYPE>::update_mean)
+        .def_buffer([](Pedestal<SUM_TYPE> &self) {
+            auto mean = self.view();
+            return py::buffer_info(
+                const_cast<SUM_TYPE *>(mean.data()), sizeof(SUM_TYPE),
+                py::format_descriptor<SUM_TYPE>::format(), 2,
+                {static_cast<py::ssize_t>(mean.shape(0)),
+                 static_cast<py::ssize_t>(mean.shape(1))},
+                {static_cast<py::ssize_t>(mean.strides()[0] * sizeof(SUM_TYPE)),
+                 static_cast<py::ssize_t>(mean.strides()[1] *
+                                          sizeof(SUM_TYPE))},
+                true);
+        });
+}

python/src/var_cluster.hpp

@@ -34,6 +34,10 @@ void define_var_cluster_finder_bindings(py::module &m) {
                  auto noise_map_span = make_view_2d(noise_map);
                  self.set_noiseMap(noise_map_span);
              })
+        .def("set_numberOfNeighbours",
+             &VarClusterFinder<double>::set_numberOfNeighbours)
+        .def("set_empty_surroundingPixels",
+             &VarClusterFinder<double>::set_empty_surroundingPixels)
         .def("set_peripheralThresholdFactor",
              &VarClusterFinder<double>::set_peripheralThresholdFactor)
         .def("find_clusters",