From 4fe587c76a83404669848705476063e7dcec914b Mon Sep 17 00:00:00 2001 From: Erik Frojdh Date: Thu, 30 Jul 2026 18:18:13 +0200 Subject: [PATCH] optimizations --- include/aare/ClusterFinder.hpp | 55 ++++++++++++++++-------- include/aare/ClusterFinderMT.hpp | 13 ++++++ include/aare/FastPedestal.hpp | 55 +++++++++++++++++++++--- include/aare/Pedestal.hpp | 26 ++++++++++- include/aare/hist/PixelHistogramImpl.hpp | 12 +++--- python/src/bind_ClusterFinder.hpp | 4 ++ python/src/bind_ClusterFinderMT.hpp | 2 + python/src/fast_pedestal.hpp | 7 +++ python/src/pedestal.hpp | 7 +++ 9 files changed, 148 insertions(+), 33 deletions(-) diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index 83b5c78..4b5c98d 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -6,6 +6,7 @@ #include "aare/NDArray.hpp" #include "aare/NDView.hpp" #include "aare/Pedestal.hpp" +#include "aare/FastPedestal.hpp" #include "aare/defs.hpp" #include @@ -26,13 +27,16 @@ class ClusterFinder { PEDESTAL_TYPE m_nSigma; const PEDESTAL_TYPE c2; const PEDESTAL_TYPE c3; - Pedestal m_pedestal; + FastPedestal m_pedestal; ClusterVector m_clusters; static const uint8_t ClusterSizeX = ClusterType::cluster_size_x; static const uint8_t ClusterSizeY = ClusterType::cluster_size_y; using CT = typename ClusterType::value_type; + NDArray m_threshold; + NDArray m_pd_corrected_frame; + public: /** * @brief Construct a new ClusterFinder object @@ -47,7 +51,8 @@ class ClusterFinder { : m_image_size(image_size), m_nSigma(nSigma), c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)), c3(sqrt(ClusterSizeX * ClusterSizeY)), - m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) { + m_pedestal(image_size[0], image_size[1]), m_clusters(capacity), + m_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; @@ -58,13 +63,25 @@ class ClusterFinder { PEDESTAL_TYPE get_nSigma() const { return m_nSigma; } void push_pedestal_frame(NDView frame) { - m_pedestal.push(frame); + if (!m_pedestal.ready()) { + m_pedestal.push_init(frame); + } else { + m_pedestal.push(frame); + } } NDArray pedestal() { return m_pedestal.mean(); } NDArray noise() { return m_pedestal.std(); } void clear_pedestal() { m_pedestal.clear(); } + /** + * @brief Refresh the cached std of the underlying pedestal. Call before + * reading the pedestal's cached std. + */ + void update_std() { m_pedestal.update_std(); } + + void update_threshold() { m_threshold = m_pedestal.std() * m_nSigma; } + /** * @brief Move the clusters from the ClusterVector in the ClusterFinder to a * new ClusterVector and return it. @@ -85,15 +102,17 @@ class ClusterFinder { // // TODO! deal with even size clusters // // currently 3,3 -> +/- 1 // // 4,4 -> +/- 2 - int dy = ClusterSizeY / 2; - int dx = ClusterSizeX / 2; - int has_center_pixel_x = + constexpr int dy = ClusterSizeY / 2; + constexpr int dx = ClusterSizeX / 2; + constexpr int has_center_pixel_x = ClusterSizeX % 2; // for even sized clusters there is no proper cluster center and // even amount of pixels around the center - int has_center_pixel_y = ClusterSizeY % 2; + constexpr int has_center_pixel_y = ClusterSizeY % 2; m_clusters.set_frame_number(frame_number); + + m_pd_corrected_frame = frame - m_pedestal.view(); for (int iy = 0; iy < frame.shape(0); iy++) { for (int ix = 0; ix < frame.shape(1); ix++) { @@ -101,10 +120,11 @@ class ClusterFinder { PEDESTAL_TYPE total = 0; // What can we short circuit here? - PEDESTAL_TYPE rms = m_pedestal.std(iy, ix); - PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix)); + // PEDESTAL_TYPE rms = m_pedestal.cached_std(iy, ix); + PEDESTAL_TYPE threshold = m_threshold(iy, ix); + PEDESTAL_TYPE value = m_pd_corrected_frame(iy, ix); - if (value < -m_nSigma * rms) + if (value < -threshold) continue; // NEGATIVE_PEDESTAL go to next pixel // TODO! No pedestal update??? @@ -113,8 +133,7 @@ class ClusterFinder { if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) { PEDESTAL_TYPE val = - frame(iy + ir, ix + ic) - - m_pedestal.mean(iy + ir, ix + ic); + m_pd_corrected_frame(iy + ir, ix + ic); total += val; max = std::max(max, val); @@ -122,15 +141,15 @@ class ClusterFinder { } } - if ((max > m_nSigma * rms)) { + if ((max > threshold)) { if (value < max) continue; // Not max go to the next pixel // but also no pedestal update - } else if (total > c3 * m_nSigma * rms) { + } else if (total > c3 * threshold) { // pass } else { // m_pedestal.push(iy, ix, frame(iy, ix)); // Safe option - m_pedestal.push_fast( + m_pedestal.push( iy, ix, frame(iy, ix)); // Assume we have reached n_samples in the @@ -160,8 +179,7 @@ class ClusterFinder { std::is_floating_point_v< PEDESTAL_TYPE>) { auto tmp = std::lround( - frame(iy + ir, ix + ic) - - m_pedestal.mean(iy + ir, ix + ic)); + m_pd_corrected_frame(iy + ir, ix + ic)); cluster.data[i] = static_cast(tmp); } // On the other hand if both are floating point @@ -169,8 +187,7 @@ class ClusterFinder { // cast directly else { auto tmp = - frame(iy + ir, ix + ic) - - m_pedestal.mean(iy + ir, ix + ic); + m_pd_corrected_frame(iy + ir, ix + ic); cluster.data[i] = static_cast(tmp); } } diff --git a/include/aare/ClusterFinderMT.hpp b/include/aare/ClusterFinderMT.hpp index d4ff001..2c54bf6 100644 --- a/include/aare/ClusterFinderMT.hpp +++ b/include/aare/ClusterFinderMT.hpp @@ -245,6 +245,19 @@ class ClusterFinderMT { } } + /** + * @brief Recompute the threshold (nSigma * pedestal std) on all cluster + * finders. Requires the processing threads to be stopped. + */ + void update_threshold() { + if (!m_processing_threads_stopped) { + throw std::runtime_error("ClusterFinderMT is still running"); + } + for (auto &cf : m_cluster_finders) { + cf->update_threshold(); + } + } + /** * @brief Return the pedestal currently used by the cluster finder * @param thread_index index of the thread diff --git a/include/aare/FastPedestal.hpp b/include/aare/FastPedestal.hpp index ab3aa96..0d245e8 100644 --- a/include/aare/FastPedestal.hpp +++ b/include/aare/FastPedestal.hpp @@ -17,6 +17,8 @@ template class FastPedestal { // TODO! Force floating point sum type? // how does the internal calculation work with integers? + bool m_ready = false; + uint32_t m_rows; uint32_t m_cols; @@ -37,15 +39,20 @@ template class FastPedestal { // Relies on having more reads than pushes to the pedestal NDArray m_mean; + // Cache std. Only refreshed via update_std() to keep push() cheap. + NDArray m_std; + public: FastPedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000) : m_rows(rows), m_cols(cols), m_samples(n_samples), m_inv_samples(1.0 / n_samples), m_sum(NDArray({rows, cols})), - m_mean(NDArray({rows, cols})) { + m_mean(NDArray({rows, cols})), + m_std(NDArray({rows, cols})) { assert(rows > 0 && cols > 0 && n_samples > 0); m_sum = Entry{SUM_TYPE(0), SUM_TYPE(0)}; m_mean = SUM_TYPE(0); + m_std = SUM_TYPE(0); } ~FastPedestal() = default; @@ -57,9 +64,15 @@ template class FastPedestal { return m_mean(row, col); } + NDArray cached_std() { return m_std; } + + SUM_TYPE cached_std(const uint32_t row, const uint32_t col) const { + return m_std(row, col); + } + SUM_TYPE variance(const uint32_t row, const uint32_t col) const { auto &entry = m_sum(row, col); - auto m2 = (entry.sum * m_inv_samples) * ((entry.sum * m_inv_samples)); + auto m2 = entry.sum * m_inv_samples * entry.sum * m_inv_samples; return entry.sum2 * m_inv_samples - m2; } @@ -87,19 +100,17 @@ template class FastPedestal { return res; } - bool ready() { return m_cur_samples == m_samples; } + bool ready() { return m_ready; } uint32_t cur_samples() { return m_cur_samples; } void clear() { m_sum = Entry{SUM_TYPE(0), SUM_TYPE(0)}; m_mean = SUM_TYPE(0); + m_std = SUM_TYPE(0); + m_ready = false; } - void clear(const uint32_t row, const uint32_t col) { - m_sum(row, col) = Entry{SUM_TYPE(0), SUM_TYPE(0)}; - m_mean(row, col) = 0; - } template void push(NDView frame) { assert(frame.size() == m_rows * m_cols); @@ -125,6 +136,11 @@ template class FastPedestal { "Frame shape does not match pedestal shape"); } + // if already full, throw an error + if (m_cur_samples == m_samples) { + throw std::runtime_error("Pedestal is full"); + } + for (size_t row = 0; row < m_rows; row++) { for (size_t col = 0; col < m_cols; col++) { const auto val = static_cast(frame(row, col)); @@ -134,6 +150,12 @@ template class FastPedestal { } } m_cur_samples += 1; + + if (m_cur_samples == m_samples) { + update_std(); + update_mean(); + m_ready = true; + } } template void push(Frame &frame) { @@ -151,6 +173,9 @@ template class FastPedestal { // their own pixel level operations) template void push(const uint32_t row, const uint32_t col, const T val_) { + if (!ready()) { + throw std::runtime_error("Pedestal is not ready, cannot push"); + } SUM_TYPE val = static_cast(val_); auto &entry = m_sum(row, col); entry.sum += val - entry.sum * m_inv_samples; @@ -160,6 +185,9 @@ template class FastPedestal { template void push_no_update(const uint32_t row, const uint32_t col, const T val_) { + if (!ready()) { + throw std::runtime_error("Pedestal is not ready, cannot push"); + } SUM_TYPE val = static_cast(val_); auto &entry = m_sum(row, col); entry.sum += val - entry.sum * m_inv_samples; @@ -178,5 +206,18 @@ template class FastPedestal { } } } + + /** + * @brief Refresh the cached std for all pixels from the current sums. + * Kept separate from push() so that pushes stay cheap; call this before + * reading cached_std(). + */ + void update_std() { + for (size_t row = 0; row < m_rows; row++) { + for (size_t col = 0; col < m_cols; col++) { + m_std(row, col) = std(row, col); + } + } + } }; } // namespace aare diff --git a/include/aare/Pedestal.hpp b/include/aare/Pedestal.hpp index 6ea6a6f..1688888 100644 --- a/include/aare/Pedestal.hpp +++ b/include/aare/Pedestal.hpp @@ -29,17 +29,22 @@ template class Pedestal { // Relies on having more reads than pushes to the pedestal NDArray m_mean; + // Cache std. Only refreshed via update_std() to keep push() cheap. + NDArray m_std; + public: Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000) : m_rows(rows), m_cols(cols), m_samples(n_samples), m_cur_samples(NDArray({rows, cols}, 0)), m_sum(NDArray({rows, cols})), m_sum2(NDArray({rows, cols})), - m_mean(NDArray({rows, cols})) { + m_mean(NDArray({rows, cols})), + m_std(NDArray({rows, cols})) { assert(rows > 0 && cols > 0 && n_samples > 0); m_sum = 0; m_sum2 = 0; m_mean = 0; + m_std = 0; } ~Pedestal() = default; @@ -51,6 +56,12 @@ template class Pedestal { return m_mean(row, col); } + NDArray cached_std() { return m_std; } + + SUM_TYPE cached_std(const uint32_t row, const uint32_t col) const { + return m_std(row, col); + } + SUM_TYPE std(const uint32_t row, const uint32_t col) const { return std::sqrt(variance(row, col)); } @@ -87,6 +98,7 @@ template class Pedestal { m_sum2 = 0; m_cur_samples = 0; m_mean = 0; + m_std = 0; } void clear(const uint32_t row, const uint32_t col) { @@ -94,6 +106,7 @@ template class Pedestal { m_sum2(row, col) = 0; m_cur_samples(row, col) = 0; m_mean(row, col) = 0; + m_std(row, col) = 0; } template void push(NDView frame) { @@ -211,6 +224,17 @@ template class Pedestal { */ void update_mean() { m_mean = m_sum / m_cur_samples; } + /** + * @brief Refresh the cached std for all pixels from the current sums. + * Kept separate from push() so pushes stay cheap; call before reading + * cached_std() (analogous to update_mean()). + */ + void update_std() { + for (uint32_t i = 0; i < m_rows * m_cols; i++) { + m_std(i / m_cols, i % m_cols) = std(i / m_cols, i % m_cols); + } + } + template void push_fast(const uint32_t row, const uint32_t col, const T val_) { // Assume we reached the steady state where all pixels have diff --git a/include/aare/hist/PixelHistogramImpl.hpp b/include/aare/hist/PixelHistogramImpl.hpp index 715537f..8008ec0 100644 --- a/include/aare/hist/PixelHistogramImpl.hpp +++ b/include/aare/hist/PixelHistogramImpl.hpp @@ -10,6 +10,7 @@ silently dropped. #include "aare/NDArray.hpp" #include "aare/NDView.hpp" +#include #include #include #include @@ -103,17 +104,16 @@ void PixelHistogramImpl::fill_unchecked(int row, int col, } int bin = static_cast((value - m_xmin) * m_scale); // Guard against floating-point rounding pushing val just below - // xmax to bin == n_bins. - if (bin >= m_n_bins) { - bin = m_n_bins - 1; - } + bin = std::clamp(bin, 0, m_n_bins - 1); + auto& cell = m_values(row, col, bin); if constexpr (std::is_integral_v) { - if (m_values(row, col, bin) >= + if (cell >= std::numeric_limits::max()) { return; } } - ++m_values(row, col, bin); + ++cell; + } template diff --git a/python/src/bind_ClusterFinder.hpp b/python/src/bind_ClusterFinder.hpp index da479dc..4a72db7 100644 --- a/python/src/bind_ClusterFinder.hpp +++ b/python/src/bind_ClusterFinder.hpp @@ -48,6 +48,10 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) { }) .def("clear_pedestal", &ClusterFinder::clear_pedestal) + .def("update_std", + &ClusterFinder::update_std) + .def("update_threshold", + &ClusterFinder::update_threshold) .def_property_readonly( "pedestal", [](ClusterFinder &self) { diff --git a/python/src/bind_ClusterFinderMT.hpp b/python/src/bind_ClusterFinderMT.hpp index 4bb8101..1a303ba 100644 --- a/python/src/bind_ClusterFinderMT.hpp +++ b/python/src/bind_ClusterFinderMT.hpp @@ -56,6 +56,8 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) { }) .def("clear_pedestal", &ClusterFinderMT::clear_pedestal) + .def("update_threshold", + &ClusterFinderMT::update_threshold) .def("sync", &ClusterFinderMT::sync) .def("stop", &ClusterFinderMT::stop) .def("start", &ClusterFinderMT::start) diff --git a/python/src/fast_pedestal.hpp b/python/src/fast_pedestal.hpp index b8815bf..73651b3 100644 --- a/python/src/fast_pedestal.hpp +++ b/python/src/fast_pedestal.hpp @@ -51,6 +51,12 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) { *standard_deviation = self.std(); return return_image_data(standard_deviation); }) + .def("cached_std", + [](FastPedestal &self) { + auto standard_deviation = new NDArray{}; + *standard_deviation = self.cached_std(); + return return_image_data(standard_deviation); + }) .def( "__array_ufunc__", [](py::object self, py::object ufunc, const std::string &method, @@ -102,6 +108,7 @@ void define_fast_pedestal_bindings(py::module &m, const std::string &name) { // }, // py::arg("frame").noconvert()) .def("update_mean", &FastPedestal::update_mean) + .def("update_std", &FastPedestal::update_std) .def_buffer([](FastPedestal &self) { auto mean = self.view(); return py::buffer_info( diff --git a/python/src/pedestal.hpp b/python/src/pedestal.hpp index 238848b..e38ba62 100644 --- a/python/src/pedestal.hpp +++ b/python/src/pedestal.hpp @@ -52,6 +52,12 @@ void define_pedestal_bindings(py::module &m, const std::string &name) { *std = self.std(); return return_image_data(std); }) + .def("cached_std", + [](Pedestal &self) { + auto standard_deviation = new NDArray{}; + *standard_deviation = self.cached_std(); + return return_image_data(standard_deviation); + }) .def( "__array_ufunc__", [](py::object self, py::object ufunc, const std::string &method, @@ -104,6 +110,7 @@ void define_pedestal_bindings(py::module &m, const std::string &name) { }, py::arg().noconvert()) .def("update_mean", &Pedestal::update_mean) + .def("update_std", &Pedestal::update_std) .def_buffer([](Pedestal &self) { auto mean = self.view(); return py::buffer_info(