diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index 33bd4636..53ac4f73 100644 --- a/include/aare/ClusterFinder.hpp +++ b/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; } /** diff --git a/include/aare/ClusterFinderMT.hpp b/include/aare/ClusterFinderMT.hpp index a4995d0d..b4f4fd22 100644 --- a/include/aare/ClusterFinderMT.hpp +++ b/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>( - 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( diff --git a/python/aare/ClusterFinder.py b/python/aare/ClusterFinder.py index 93af877a..9c0151b8 100644 --- a/python/aare/ClusterFinder.py +++ b/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): diff --git a/python/src/bind_ClusterFinder.hpp b/python/src/bind_ClusterFinder.hpp index c13d3305..2e0bf8b3 100644 --- a/python/src/bind_ClusterFinder.hpp +++ b/python/src/bind_ClusterFinder.hpp @@ -32,8 +32,10 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) { py::class_>( m, class_name.c_str()) - .def(py::init, pd_type, size_t>(), py::arg("image_size"), - py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000) + .def(py::init, 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", diff --git a/python/src/bind_ClusterFinderMT.hpp b/python/src/bind_ClusterFinderMT.hpp index 252cc2f8..37ab8d47 100644 --- a/python/src/bind_ClusterFinderMT.hpp +++ b/python/src/bind_ClusterFinderMT.hpp @@ -32,10 +32,11 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) { py::class_>( m, class_name.c_str()) - .def(py::init, pd_type, size_t, size_t, size_t>(), + .def(py::init, 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 &self, py::array_t frame) { diff --git a/src/ClusterFinderMT.test.cpp b/src/ClusterFinderMT.test.cpp index fae8b247..ca630de4 100644 --- a/src/ClusterFinderMT.test.cpp +++ b/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( - 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; ClusterFinderMTWrapper cf( {static_cast(file.rows()), static_cast(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()); + + for (size_t i = 0; i < n_frames; ++i) { + frame = file.read_frame(); cf.find_clusters(frame.view()); } @@ -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 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 "