Fix/clusterfinder test (#359)

- fixed test in ClusterFinderMT (use minimum one pedestal to accumulate
dummy pedestal)
- added minimum_pedestal_samples as additional constructor argument in
ClusterFinder and ClusterFinderMT
- bugfix: race conditions for empty queues in collect() - check directly
when needed

include/aare/ClusterFinder.hpp

@@ -50,18 +50,21 @@ class ClusterFinder {
      * @param image_size Image shape as (rows, columns).
      * @param nSigma Per-pixel noise threshold multiplier.
      * @param capacity Initial cluster-vector capacity.
+     * @param min_pedestal_samples Minimum number of pedestal frames to
+     * accumulate to get reasonable statistics
      */
     ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
-                  size_t capacity = 1000000)
+                  size_t capacity = 1000000, size_t min_pedestal_samples = 1000)
         : m_image_size(image_size), m_nSigma(nSigma),
           c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)),
           c3(sqrt(ClusterSizeX * ClusterSizeY)),
-          m_pedestal(image_size[0], image_size[1]), m_clusters(capacity),
-          m_threshold({image_size[0], image_size[1]}, 0),
+          m_pedestal(image_size[0], image_size[1], min_pedestal_samples),
+          m_clusters(capacity), m_threshold({image_size[0], image_size[1]}, 0),
           m_pd_corrected_frame({image_size[0], image_size[1]}, 0) {
         LOG(logDEBUG) << "ClusterFinder: "
                       << "image_size: " << image_size[0] << "x" << image_size[1]
-                      << ", nSigma: " << nSigma << ", capacity: " << capacity;
+                      << ", nSigma: " << nSigma << ", capacity: " << capacity
+                      << ", min_pedestal_samples: " << min_pedestal_samples;
     }
 
     /**

include/aare/ClusterFinderMT.hpp

@@ -131,17 +131,26 @@ class ClusterFinderMT {
         }
     }
 
+    bool output_queues_are_empty() const {
+        for (auto &q : m_output_queues) {
+            if (!q->isEmpty()) {
+                return false;
+            }
+        }
+        return true;
+    }
+
     /**
      * @brief Collect all the clusters from the output queues and write them to
      * the sink
      */
     void collect() {
-        bool empty = true;
         Backoff backoff;
-        while (!m_stop_requested || !empty || !m_processing_threads_stopped) {
+        while (!m_stop_requested || !output_queues_are_empty() ||
+               !m_processing_threads_stopped) {
             bool moved_any = false;
             for (auto &queue : m_output_queues) {
-                while (auto *front = queue->frontPtr()) {
+                if (auto *front = queue->frontPtr(); front != nullptr) {
                     while (!m_sink.write(std::move(*front))) {
                         backoff.pause();
                     }
@@ -149,7 +158,6 @@ class ClusterFinderMT {
                     moved_any = true;
                 }
             }
-            empty = !moved_any;
             if (moved_any) {
                 backoff.reset();
             } else {
@@ -171,16 +179,19 @@ class ClusterFinderMT {
      * allocated once and recycled, so the total resident frame memory is
      * n_threads * queue_depth * frame size. Keeping the in flight data below
      * the L3 size keeps the per frame copy cheap.
+     * @param min_pedestal_samples minimum number of pedestal samples to
+     * accumulate before using the pedestal
      */
     ClusterFinderMT(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
                     size_t capacity = 2000, size_t n_threads = 3,
-                    size_t queue_depth = 16)
+                    size_t queue_depth = 16, size_t min_pedestal_samples = 1000)
         : m_n_threads(n_threads) {
 
         LOG(logDEBUG1) << "ClusterFinderMT: "
                        << "image_size: " << image_size[0] << "x"
                        << image_size[1] << ", nSigma: " << nSigma
                        << ", capacity: " << capacity
+                       << ", min_pedestal_samples: " << min_pedestal_samples
                        << ", n_threads: " << n_threads
                        << ", queue_depth: " << queue_depth;
 
@@ -188,7 +199,7 @@ class ClusterFinderMT {
             m_cluster_finders.push_back(
                 std::make_unique<
                     ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>>(
-                    image_size, nSigma, capacity));
+                    image_size, nSigma, capacity, min_pedestal_samples));
         }
         for (size_t i = 0; i < n_threads; i++) {
             m_frame_pools.emplace_back(

python/aare/ClusterFinder.py

@@ -18,24 +18,55 @@ def _get_class(name, cluster_size, dtype):
 
 
 
-def ClusterFinder(image_size, cluster_size=(3,3), n_sigma=5, dtype = np.int32, capacity = 1024):
+def ClusterFinder(image_size, *, cluster_size=(3,3), n_sigma=5, dtype = np.int32, capacity = 1024, min_pedestal_samples = 1000):
     """
     Factory function to create a ClusterFinder object. Provides a cleaner syntax for 
     the templated ClusterFinder in C++.
+
+    Parameters
+    ----------
+    image_size : tuple
+        The size of the image as a tuple (height, width).
+    cluster_size : tuple, optional
+        The size of the cluster to find as a tuple (height, width). Default is (3,3).
+    n_sigma : int, optional
+        Multiplier of the standard deviation used as a threshold to identify potential photon pixels. Default is 5.
+    dtype : data-type, optional
+        The data type of the image. Default is np.int32.
+    capacity : int, optional
+        The maximum number of clusters than can be stored before reallocating. Default is 1024.
+    min_pedestal_samples : int, optional
+        The minimum number of pedestal samples to accumulate before using the pedestal. Default is 1000.
     """
     cls = _get_class("ClusterFinder", cluster_size, dtype)
-    return cls(image_size, n_sigma=n_sigma, capacity=capacity)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity, min_pedestal_samples=min_pedestal_samples)
 
 
-
-def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3): 
+def ClusterFinderMT(image_size, *, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3, min_pedestal_samples = 1000): 
     """ 
     Factory function to create a ClusterFinderMT object. Provides a cleaner syntax for 
     the templated ClusterFinderMT in C++.
+
+    Parameters
+    ----------
+    image_size : tuple
+        The size of the image as a tuple (height, width).
+    cluster_size : tuple, optional
+        The size of the cluster to find as a tuple (height, width). Default is (3,3).
+    n_sigma : int, optional
+        Multiplier of the standard deviation used as a threshold to identify potential photon pixels. Default is 5.
+    dtype : data-type, optional
+        The data type of the image. Default is np.int32.
+    capacity : int, optional
+        The maximum number of clusters than can be stored before reallocating. Default is 1024.
+    n_threads : int, optional
+        The number of threads to use for processing. Default is 3.
+    min_pedestal_samples : int, optional
+        The minimum number of pedestal samples to accumulate before using the pedestal. Default is 1000.
     """
 
     cls = _get_class("ClusterFinderMT", cluster_size, dtype)
-    return cls(image_size, n_sigma=n_sigma, capacity=capacity, n_threads=n_threads)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity, min_pedestal_samples=min_pedestal_samples, n_threads=n_threads)
 
 
 def ClusterCollector(clusterfindermt, dtype=np.int32): 

python/src/bind_ClusterFinder.hpp

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

python/src/bind_ClusterFinderMT.hpp

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

src/ClusterFinderMT.test.cpp

@@ -22,9 +22,11 @@ class ClusterFinderMTWrapper
   public:
     ClusterFinderMTWrapper(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0,
                            size_t capacity = 2000, size_t n_threads = 3,
-                           size_t queue_depth = 16)
+                           size_t queue_depth = 16,
+                           size_t minimum_pedestal_samples = 1000)
         : ClusterFinderMT<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>(
-              image_size, nSigma, capacity, n_threads, queue_depth) {}
+              image_size, nSigma, capacity, n_threads, queue_depth,
+              minimum_pedestal_samples) {}
 
     size_t get_m_input_queues_size() const {
         return this->m_input_queues.size();
@@ -73,13 +75,14 @@ TEST_CASE("multithreaded cluster finder", "[.with-data]") {
     File file(fpath);
 
     size_t n_threads = 2;
-    size_t n_frames_pd = 10;
+    size_t minimum_pedestal_samples = 1;
+    size_t n_frames = 10;
 
     using ClusterType = Cluster<int32_t, 3, 3>;
 
     ClusterFinderMTWrapper<ClusterType> cf(
         {static_cast<int64_t>(file.rows()), static_cast<int64_t>(file.cols())},
-        5, 2000, n_threads); // no idea what frame type is!!! default uint16_t
+        5, 2000, n_threads, 16, minimum_pedestal_samples);
 
     CHECK(cf.get_m_input_queues_size() == n_threads);
     CHECK(cf.get_m_output_queues_size() == n_threads);
@@ -87,8 +90,11 @@ TEST_CASE("multithreaded cluster finder", "[.with-data]") {
     CHECK(cf.m_output_queues_are_empty() == true);
     CHECK(cf.m_input_queues_are_empty() == true);
 
-    for (size_t i = 0; i < n_frames_pd; ++i) {
-        auto frame = file.read_frame();
+    auto frame = file.read_frame();
+    cf.push_pedestal_frame(frame.view<uint16_t>());
+
+    for (size_t i = 0; i < n_frames; ++i) {
+        frame = file.read_frame();
         cf.find_clusters(frame.view<uint16_t>());
     }
 
@@ -97,7 +103,7 @@ TEST_CASE("multithreaded cluster finder", "[.with-data]") {
     CHECK(cf.m_output_queues_are_empty() == true);
     CHECK(cf.m_input_queues_are_empty() == true);
 
-    CHECK(cf.m_sink_size() == n_frames_pd);
+    CHECK(cf.m_sink_size() == n_frames);
     ClusterCollector<ClusterType> clustercollector(&cf);
 
     clustercollector.stop();
@@ -105,7 +111,6 @@ TEST_CASE("multithreaded cluster finder", "[.with-data]") {
     CHECK(cf.m_sink_size() == 0);
 
     auto clustervec = clustercollector.steal_clusters();
-    // CHECK(clustervec.size() == ) //dont know how many clusters to expect
 }
 
 TEST_CASE("frame buffers are recycled when pushing more frames than the pool "