From f670ba77a28c3c5467b90932e81ddddfb4ba7f25 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Erik=20Fr=C3=B6jdh?= Date: Tue, 9 Jun 2026 09:08:48 +0200 Subject: [PATCH] PixelHistogram (#317) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Multi threaded filling of per pixel histograms for example for detector calibration 1. PixelHistogram - Generic variant expects already pedestal subtracted data 2. PedestalTrackingHistogram - Terrible name, useful class. Keeps it's own pedestal and does conversion and pedestal tracking in the worker threads. --------- Co-authored-by: Lars Erik Fröjd --- CMakeLists.txt | 21 +- RELEASE.md | 1 + docs/src/index.rst | 3 +- docs/src/python/histogram/index.rst | 9 + .../pyPedestalTrackingPixelHistogram.rst | 24 + .../src/python/histogram/pyPixelHistogram.rst | 41 ++ include/aare/Pedestal.hpp | 28 + .../hist/PedestalTrackingPixelHistogram.hpp | 133 ++++ include/aare/hist/PixelHistogram.hpp | 342 ++++++++++ include/aare/hist/PixelHistogramImpl.hpp | 143 +++++ python/aare/__init__.py | 12 +- python/aare/utils.py | 17 +- .../bind_PedestalTrackingPixelHistogram.hpp | 245 ++++++++ python/src/bind_PixelHistogram.hpp | 147 +++++ python/src/module.cpp | 4 + python/src/pedestal.hpp | 27 + src/hist/PedestalTrackingPixelHistogram.cpp | 584 ++++++++++++++++++ src/hist/PixelHistogram.test.cpp | 319 ++++++++++ src/hist/PixelHistogramImpl.test.cpp | 140 +++++ 19 files changed, 2231 insertions(+), 9 deletions(-) create mode 100644 docs/src/python/histogram/index.rst create mode 100644 docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst create mode 100644 docs/src/python/histogram/pyPixelHistogram.rst create mode 100644 include/aare/hist/PedestalTrackingPixelHistogram.hpp create mode 100644 include/aare/hist/PixelHistogram.hpp create mode 100644 include/aare/hist/PixelHistogramImpl.hpp create mode 100644 python/src/bind_PedestalTrackingPixelHistogram.hpp create mode 100644 python/src/bind_PixelHistogram.hpp create mode 100644 src/hist/PedestalTrackingPixelHistogram.cpp create mode 100644 src/hist/PixelHistogram.test.cpp create mode 100644 src/hist/PixelHistogramImpl.test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 36a17ae..7eac5ec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -379,6 +379,9 @@ set(PUBLICHEADERS include/aare/FileInterface.hpp include/aare/FilePtr.hpp include/aare/Frame.hpp + include/aare/hist/PixelHistogram.hpp + include/aare/hist/PixelHistogramImpl.hpp + include/aare/hist/PedestalTrackingPixelHistogram.hpp include/aare/GainMap.hpp include/aare/ROIGeometry.hpp include/aare/DetectorGeometry.hpp @@ -403,6 +406,7 @@ set(SourceFiles ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ROIGeometry.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PedestalTrackingPixelHistogram.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/FilePtr.cpp @@ -420,7 +424,7 @@ set(SourceFiles ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp) -add_library(aare_core STATIC ${SourceFiles}) +add_library(aare_core STATIC ${SourceFiles} ${PUBLICHEADERS}) target_include_directories( aare_core PUBLIC "$" "$") @@ -455,9 +459,8 @@ if(AARE_CUSTOM_ASSERT) target_compile_definitions(aare_core PUBLIC AARE_CUSTOM_ASSERT) endif() -set_target_properties( - aare_core PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR} - PUBLIC_HEADER "${PUBLICHEADERS}") +set_target_properties(aare_core PROPERTIES ARCHIVE_OUTPUT_DIRECTORY + ${CMAKE_BINARY_DIR}) if(AARE_TESTS) set(TestSources @@ -479,6 +482,8 @@ 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/hist/PixelHistogramImpl.test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PixelHistogram.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp @@ -494,8 +499,12 @@ if(AARE_MASTER_PROJECT) TARGETS aare_core aare_compiler_flags EXPORT "${TARGETS_EXPORT_NAME}" LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} - ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} - PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/aare) + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}) + install( + DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/include/aare/ + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/aare + FILES_MATCHING + PATTERN "*.hpp") endif() set(CMAKE_INSTALL_RPATH $ORIGIN) diff --git a/RELEASE.md b/RELEASE.md index cefe4de..708f4ba 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -9,6 +9,7 @@ - setter and getter for nSigma for ClusterFinder ``aare.ClusterFinder().nSigma = 2``, ``aare.ClusterFinderMT().set_nSigma(2)`` - mask opeartor for ClusterVector ``masked_clustervector = aare.ClusterVector()(mask)`` - passing pre computed eta values to ``aare.Interpolator.interpolate`` alongside clusters +- Added ``PixelHistogram`` and ``PedestalTrackingPixelHistogram`` ### Bugfixes: diff --git a/docs/src/index.rst b/docs/src/index.rst index 834c490..8b9af28 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -28,6 +28,7 @@ AARE pycalibration python/cluster/index python/file/index + python/histogram/index pyFit @@ -63,4 +64,4 @@ AARE Philosophy Workflow - Tests \ No newline at end of file + Tests diff --git a/docs/src/python/histogram/index.rst b/docs/src/python/histogram/index.rst new file mode 100644 index 0000000..5ca8b33 --- /dev/null +++ b/docs/src/python/histogram/index.rst @@ -0,0 +1,9 @@ +Histogram +========= + +.. toctree:: + :caption: Histogram + :maxdepth: 1 + + pyPixelHistogram + pyPedestalTrackingPixelHistogram diff --git a/docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst b/docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst new file mode 100644 index 0000000..18d438c --- /dev/null +++ b/docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst @@ -0,0 +1,24 @@ +PedestalTrackingPixelHistogram +============================== + +.. warning:: + + ``PedestalTrackingPixelHistogram`` is specifically designed for use in the Jungfrau calibration + pipeline. Make sure you understand the behaviour before using it in other contexts. + +``PedestalTrackingPixelHistogram`` accumulates a pixel-wise histogram of +``frame - pedestal`` residuals while maintaining a running per-pixel pedestal +estimate. + +Use ``push_pedestal_no_update()`` to seed the pedestal estimate, then +``update_mean()`` before submitting frames with ``fill_async()``. Pending +asynchronous fills are drained by ``flush()``, and snapshot methods such as +``values()`` and ``pedestal_mean()`` return numpy arrays. + +.. py:currentmodule:: aare + +.. autoclass:: PedestalTrackingPixelHistogram + :members: + :undoc-members: + :show-inheritance: + :inherited-members: diff --git a/docs/src/python/histogram/pyPixelHistogram.rst b/docs/src/python/histogram/pyPixelHistogram.rst new file mode 100644 index 0000000..2705732 --- /dev/null +++ b/docs/src/python/histogram/pyPixelHistogram.rst @@ -0,0 +1,41 @@ +PixelHistogram +============== + +``PixelHistogram`` accumulates one histogram per detector pixel from 2D +``float64`` images. The public ``PixelHistogram`` name is an alias for +``PixelHistogram_d``, which stores bin counts as ``float64``. + +.. note:: + + ``PixelHistogram`` is designed for fast filling from 2D images, utilizing multiple threads, + and contiguous storage of a specific type. No over/underflow bins are provided. For generic use cases consider + boost-histogram. (https://boost-histogram.readthedocs.io/en/latest/) + + +Storage-specific variants are also available when a smaller or integer count +type is preferred: + +* ``PixelHistogram_d``: ``float64`` storage +* ``PixelHistogram_f``: ``float32`` storage +* ``PixelHistogram_u64``: ``uint64`` storage +* ``PixelHistogram_u32``: ``uint32`` storage +* ``PixelHistogram_u16``: ``uint16`` storage +* ``PixelHistogram_u8``: ``uint8`` storage + +.. py:currentmodule:: aare + +.. autoclass:: PixelHistogram + :members: + :undoc-members: + :show-inheritance: + :inherited-members: + +Showing API for PixelHistogram_d, all variant share the same API. + +.. autoclass:: aare._aare.PixelHistogram_d + :members: + :undoc-members: + :show-inheritance: + :inherited-members: + :noindex: + diff --git a/include/aare/Pedestal.hpp b/include/aare/Pedestal.hpp index 5753176..6ea6a6f 100644 --- a/include/aare/Pedestal.hpp +++ b/include/aare/Pedestal.hpp @@ -45,6 +45,8 @@ template class Pedestal { NDArray mean() { return m_mean; } + const NDView view() const { return m_mean.view(); } + SUM_TYPE mean(const uint32_t row, const uint32_t col) const { return m_mean(row, col); } @@ -110,6 +112,32 @@ template class Pedestal { } } + template + void push_with_threshold(const NDView frame, + const NDView threshold) { + assert(frame.size() == m_rows * m_cols); + + // TODO! move away from m_rows, m_cols + if (frame.shape() != std::array{m_rows, m_cols}) { + throw std::runtime_error( + "Frame shape does not match pedestal shape"); + } + + if (threshold.shape() != std::array{m_rows, m_cols}) { + throw std::runtime_error( + "Threshold shape does not match pedestal shape"); + } + + for (size_t row = 0; row < m_rows; row++) { + for (size_t col = 0; col < m_cols; col++) { + if (fabs(frame(row, col) - mean(row, col)) < + threshold(row, col)) { + push(row, col, frame(row, col)); + } + } + } + } + /** * Push but don't update the cached mean. Speeds up the process * when initializing the pedestal. diff --git a/include/aare/hist/PedestalTrackingPixelHistogram.hpp b/include/aare/hist/PedestalTrackingPixelHistogram.hpp new file mode 100644 index 0000000..263ea45 --- /dev/null +++ b/include/aare/hist/PedestalTrackingPixelHistogram.hpp @@ -0,0 +1,133 @@ +#pragma once +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" +#include "aare/Pedestal.hpp" +#include "aare/ProducerConsumerQueue.hpp" +#include "aare/hist/PixelHistogramImpl.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace aare { + +class PedestalTrackingPixelHistogram { + public: + using StorageType = uint16_t; + using AxisType = float; // TODO: template on pedestal type if needed + using FrameType = uint16_t; + + private: + using Hist = PixelHistogramImpl; + using AsyncQueue = ProducerConsumerQueue>; + + // What kind of fan-out work the worker pool should currently do. + // Set under work_mutex_; read by worker_loop after wakeup. + enum class WorkKind { PushPedestal, UpdateMean, FillWithThreshold }; + + int rows_; + int cols_; + int n_threads_; + const AxisType xmin_; + const AxisType xmax_; + + std::vector row_offsets_; + std::vector partial_hists_; + // Per-thread pedestal sized [local_rows x cols]. Indexed by the + // worker using the LOCAL row index (i.e. 0..row_count(t)-1), NOT the + // global row index. Owned exclusively by worker `t` during a + // dispatched fan-out. + std::vector> partial_pedestals_; + std::vector> partial_std_; // cached for pedestal + // tracking + + // Thread pool members + std::vector workers_; + std::mutex work_mutex_; + std::condition_variable work_cv_; + std::condition_variable done_cv_; + WorkKind current_work_kind_; + const NDView *current_image_; + const std::vector> *current_images_; + int completed_threads_; + std::atomic stop_workers_; + int work_generation_; + + // Serialises all fan-outs: `push_pedestal_no_update`, `update_mean`, + // and the async coordinator's batch fill dispatches. + // Always the outermost lock; work_mutex_ is taken briefly inside it. + mutable std::mutex fill_mutex_; + + // Async producer/consumer pipeline. SPSC queue feeds the coordinator + // thread, which batches queued images before dispatching them. + std::unique_ptr async_queue_; + std::thread coordinator_; + std::atomic stop_coordinator_{false}; + std::atomic coordinator_busy_{false}; + std::atomic completed_async_fills_{0}; + std::chrono::microseconds async_wait_{100}; + std::size_t max_batch_size_; + std::atomic n_sigma_; + + // Private worker thread method + void worker_loop(int thread_id); + void coordinator_loop(); + int row_start(int thread_id) const; + int row_count(int thread_id) const; + + // Sets up the work generation, wakes workers, and waits for them. + // Caller MUST already hold `fill_mutex_`. `image` may be nullptr + // for work kinds that don't need it (e.g. UpdateMean). + void dispatch_(WorkKind kind, const NDView *image); + void dispatch_fill_batch_(const std::vector> &images); + + // Coordinator-facing entry point: takes fill_mutex_ and dispatches + // FillWithThreshold batches to the worker pool. Only ever called by the + // coordinator thread, on images already shape-checked by fill_async. + void fill_with_threshold_batch_(std::vector> &batch); + + public: + PedestalTrackingPixelHistogram(int rows, int cols, int n_bins, + AxisType xmin, AxisType xmax, + int n_threads = 1, + std::size_t max_pending = 16, + AxisType n_sigma = 1.0); + ~PedestalTrackingPixelHistogram(); + + void push_pedestal_no_update(const NDView &frame); + void update_mean(); + NDArray pedestal_mean() const; + + void fill_async(NDArray &&image); + + void fill_from_file(const std::filesystem::path &fname, + ssize_t max_frames = -1, bool verbose = false); + + void process_pedestal_file(const std::filesystem::path &fname, + ssize_t max_frames = -1, bool verbose = false); + + // Sigma multiplier for the pedestal-update gate in + // fill_async. Atomic; safe to read/write at any + // time (the new value takes effect on subsequent pixel evaluations). + PedestalTrackingPixelHistogram::AxisType n_sigma() const; + void set_n_sigma(PedestalTrackingPixelHistogram::AxisType n_sigma); + + // Wait for all queued async fills to complete. Cheap when the queue + // is already drained. + void flush() const; + + // Implicitly flushes pending async fills first so the snapshot is + // consistent with everything that was submitted up to the call. + NDArray values() const; + NDArray bin_centers() const; + NDArray bin_edges() const; +}; + +} // namespace aare diff --git a/include/aare/hist/PixelHistogram.hpp b/include/aare/hist/PixelHistogram.hpp new file mode 100644 index 0000000..05432d6 --- /dev/null +++ b/include/aare/hist/PixelHistogram.hpp @@ -0,0 +1,342 @@ +#pragma once +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" +#include "aare/ProducerConsumerQueue.hpp" +#include "aare/hist/PixelHistogramImpl.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace aare { +template +class PixelHistogram { + private: + using Hist = PixelHistogramImpl; + using AsyncQueue = ProducerConsumerQueue>; + + int rows_; + int cols_; + int n_threads_; + const AxisType xmin_; + const AxisType xmax_; + std::vector partial_hists_; + // Cumulative row offsets so that thread t owns rows + // [row_offsets_[t], row_offsets_[t + 1]) + // Length is n_threads_ + 1; partition is balanced (the first + // rows_ % n_threads_ threads get one extra row each). + std::vector row_offsets_; + + // Thread pool members + std::vector workers_; + std::mutex work_mutex_; + std::condition_variable work_cv_; + std::condition_variable done_cv_; + const NDView *current_image_; + int completed_threads_; + std::atomic stop_workers_; + int work_generation_; + + // Async producer/consumer pipeline. SPSC queue feeds the coordinator + // thread, which fans each image out to the worker pool one at a time. + // TODO: batch processing? + // TODO: FIFO to avoid allocations? + std::unique_ptr async_queue_; + std::thread coordinator_; + std::atomic stop_coordinator_{false}; + std::atomic coordinator_busy_{false}; + std::chrono::microseconds async_wait_{100}; + + // Private worker thread method + void worker_loop(int thread_id); + void coordinator_loop(); + // Fan a single image out to the worker pool and block until every + // worker has merged its row band. Only ever called by the coordinator + // thread, so no caller-serialisation lock is needed. + void dispatch(const NDView &image); + int row_start(int thread_id) const; + int row_count(int thread_id) const; + + public: + PixelHistogram(int rows, int cols, int n_bins, AxisType xmin, AxisType xmax, + int n_threads = 1, std::size_t max_pending = 16); + ~PixelHistogram(); + + // Asynchronous fill: takes ownership of `image`, enqueues it for the + // coordinator thread, and returns. Blocks the caller only if the queue + // is full (single-producer, single-consumer queue with a sleep-poll + // backpressure loop, matching the convention in ClusterFinderMT). + void fill_async(NDArray &&image); + + // Wait for all queued async fills to complete. Cheap when the queue + // is already drained. + void flush() const; + + // Implicitly flushes pending async fills first so the snapshot is + // consistent with everything that was submitted up to the call. + NDArray values() const; + NDArray bin_centers() const; + NDArray bin_edges() const; +}; + +template +PixelHistogram::PixelHistogram(int rows, int cols, + int n_bins, AxisType xmin, + AxisType xmax, + int n_threads, + std::size_t max_pending) + : rows_(rows), cols_(cols), n_threads_(n_threads), xmin_(xmin), xmax_(xmax), + current_image_(nullptr), completed_threads_(0), stop_workers_(false), + work_generation_(0) { + if (rows_ < 1 || cols_ < 1 || n_bins < 1) { + throw std::invalid_argument( + "PixelHistogram requires positive rows, cols and bins"); + } + if (n_threads < 1) { + throw std::invalid_argument( + "PixelHistogram requires at least one thread"); + } + if (max_pending < 1) { + throw std::invalid_argument("PixelHistogram requires max_pending >= 1"); + } + + n_threads_ = std::min(n_threads_, rows_); + + // Build a balanced row partition. With base = rows_ / n_threads_ and + // extra = rows_ % n_threads_, the first `extra` threads get base + 1 + // rows each and the rest get `base` rows. This avoids the + // ceil(rows/n_threads) scheme leaving trailing threads with zero or + // negative row counts (e.g. rows=17, n_threads=8). + row_offsets_.resize(n_threads_ + 1); + const int base = rows_ / n_threads_; + const int extra = rows_ % n_threads_; + int offset = 0; + for (int i = 0; i < n_threads_; ++i) { + row_offsets_[i] = offset; + offset += base + (i < extra ? 1 : 0); + } + row_offsets_[n_threads_] = offset; // == rows_ by construction + + // Initialize partial histograms for each thread + partial_hists_.reserve(n_threads_); + for (int i = 0; i < n_threads_; ++i) { + const auto local_rows = row_count(i); + partial_hists_.emplace_back(local_rows, cols, n_bins, xmin, xmax); + } + + // Spawn worker threads + for (int i = 0; i < n_threads_; ++i) { + workers_.emplace_back([this, i]() { this->worker_loop(i); }); + } + + // Async pipeline. The PCQ holds (size - 1) usable slots, so size up by + // one to honour the requested max_pending. + async_queue_ = std::make_unique( + static_cast(max_pending + 1)); + coordinator_ = std::thread([this]() { this->coordinator_loop(); }); +} + +template +PixelHistogram::~PixelHistogram() { + // Drain any pending async fills before tearing down the worker pool. + // The coordinator's loop keeps processing while stop_coordinator_ is + // true as long as the queue is non-empty (mirrors ClusterFinderMT). + if (coordinator_.joinable()) { + stop_coordinator_ = true; + coordinator_.join(); + } + + // Signal all workers to stop + { + std::unique_lock lock(work_mutex_); + stop_workers_ = true; + } + work_cv_.notify_all(); + + // Join all worker threads + for (auto &thread : workers_) { + if (thread.joinable()) { + thread.join(); + } + } +} + +template +int PixelHistogram::row_start(int thread_id) const { + return row_offsets_[thread_id]; +} + +template +int PixelHistogram::row_count(int thread_id) const { + return row_offsets_[thread_id + 1] - row_offsets_[thread_id]; +} + +template +void PixelHistogram::worker_loop(int thread_id) { + int last_generation = 0; + + while (true) { + std::unique_lock lock(work_mutex_); + work_cv_.wait(lock, [this, last_generation]() { + return work_generation_ != last_generation || stop_workers_; + }); + + if (stop_workers_) { + break; + } + + // Get work assignment + const NDView &image = *current_image_; + const int generation = work_generation_; + const int first_row = row_start(thread_id); + const int local_rows = row_count(thread_id); + + lock.unlock(); + + // Do the work: fill this thread's partial histogram. The + // [xmin, xmax) range gate lives inside PixelHistogramImpl::fill. + auto &my_hist = partial_hists_[thread_id]; + const auto cols = image.shape(1); + for (int local_row = 0; local_row < local_rows; ++local_row) { + const auto row = static_cast(first_row + local_row); + for (ssize_t col = 0; col < cols; ++col) { + const auto val = image(row, col); + my_hist.fill_unchecked(local_row, static_cast(col), val); + } + } + + // Signal completion + { + std::unique_lock done_lock(work_mutex_); + last_generation = generation; + completed_threads_++; + if (completed_threads_ == n_threads_) { + done_cv_.notify_one(); + } + } + } +} + +template +NDArray PixelHistogram::values() const { + // Make sure any pending async fills are merged in before we snapshot + // the partial histograms. Cheap when the queue is already drained. + flush(); + + const auto first_shard_view = partial_hists_.front().view(); + const auto cols = static_cast(first_shard_view.shape(1)); + const auto bins = static_cast(first_shard_view.shape(2)); + const auto rows = static_cast(rows_); + + NDArray data({rows, cols, bins}); + + // Each thread owns a disjoint, contiguous range of rows. The shard's + // dense storage layout [local_row][col][bin] is identical to the slice + // [first_row .. first_row + local_rows)[col][bin] of `data`, so the + // merge is just one bulk copy per thread; no per-element accumulation + // and no upfront zeroing of `data` is needed. + const size_t pixel_stride = static_cast(cols) * bins; + for (int t = 0; t < n_threads_; ++t) { + const auto first_row = static_cast(row_start(t)); + const auto local_rows = static_cast(row_count(t)); + if (local_rows == 0) + continue; + + const auto shard_view = partial_hists_[t].view(); + std::memcpy(data.data() + first_row * pixel_stride, shard_view.data(), + local_rows * pixel_stride * sizeof(StorageType)); + } + + return data; +} + +template +void PixelHistogram::dispatch( + const NDView &image) { + // Called only by the coordinator thread on images already shape-checked + // by fill_async, so there is no need to re-validate or to serialise + // against other callers. + + // Reset counters and set work + { + std::unique_lock lock(work_mutex_); + completed_threads_ = 0; + current_image_ = ℑ + ++work_generation_; + } + + // Signal all worker threads to start + work_cv_.notify_all(); + + // Wait for all workers to complete + { + std::unique_lock lock(work_mutex_); + done_cv_.wait(lock, + [this]() { return completed_threads_ == n_threads_; }); + current_image_ = nullptr; // Clear work assignment + } +} + +template +void PixelHistogram::fill_async( + NDArray &&image) { + if (image.shape(0) != rows_ || image.shape(1) != cols_) { + throw std::invalid_argument( + "PixelHistogram image shape does not match constructor shape"); + } + + // SPSC backpressure: spin with a short sleep until a slot frees up. + // The std::move only consumes `image` on the iteration that succeeds + // (placement-new inside write() runs only when the slot is free). + while (!async_queue_->write(std::move(image))) { + std::this_thread::sleep_for(async_wait_); + } +} + +template +void PixelHistogram::flush() const { + while (!async_queue_->isEmpty() || + coordinator_busy_.load(std::memory_order_acquire)) { + std::this_thread::sleep_for(async_wait_); + } +} + +template +void PixelHistogram::coordinator_loop() { + NDArray item; + while (!stop_coordinator_.load(std::memory_order_acquire) || + !async_queue_->isEmpty()) { + if (async_queue_->read(item)) { + coordinator_busy_.store(true, std::memory_order_release); + dispatch(item.view()); + coordinator_busy_.store(false, std::memory_order_release); + } else { + std::this_thread::sleep_for(async_wait_); + } + } +} + +template +NDArray +PixelHistogram::bin_centers() const { + // All shards share the same value-axis configuration, so any one will + // do; pick the first. + return partial_hists_.front().bin_centers(); +} + +template +NDArray PixelHistogram::bin_edges() const { + return partial_hists_.front().bin_edges(); +} + +} // namespace aare diff --git a/include/aare/hist/PixelHistogramImpl.hpp b/include/aare/hist/PixelHistogramImpl.hpp new file mode 100644 index 0000000..715537f --- /dev/null +++ b/include/aare/hist/PixelHistogramImpl.hpp @@ -0,0 +1,143 @@ +#pragma once +/* +Basic pixel histogram class with templated axis and storage type. +row, col are integers and the per pixel histogram axis is AxisType. + +Bin policy: regular binning on [xmin, xmax). Values outside this range are +silently dropped. +*/ + +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" + +#include +#include +#include +#include + +namespace aare { + +template class PixelHistogramImpl { + static_assert(std::is_floating_point_v, + "PixelHistogramImpl requires a floating-point axis type"); + NDArray m_values; + NDArray m_edges; + int m_rows; + int m_cols; + int m_n_bins; + T m_xmin; + T m_xmax; + // n_bins / (xmax - xmin), precomputed to keep the hot path + // division-free. + T m_scale; + + public: + PixelHistogramImpl(int rows, int cols, int n_bins, T xmin, T xmax); + + void fill(const NDView &frame); + void fill(int row, int col, T value); + void fill_unchecked(int row, int col, T value); + + NDArray values() const; + // Zero-copy view of the underlying [rows x cols x n_bins] storage. + // Lifetime is tied to *this. Use for low-level merge/stitching paths; + // prefer values() for the public API where you want an owned copy. + NDView view() const; + NDArray bin_centers() const; + NDArray bin_edges() const; +}; + +template +PixelHistogramImpl::PixelHistogramImpl(int rows, int cols, + int n_bins, T xmin, + T xmax) + : m_values(NDArray({static_cast(rows), + static_cast(cols), + static_cast(n_bins)}, + StorageType{0})), + m_edges(NDArray({static_cast(n_bins + 1)})), m_rows(rows), + m_cols(cols), m_n_bins(n_bins), m_xmin(xmin), m_xmax(xmax), + m_scale(static_cast(n_bins) / (xmax - xmin)) { + if (rows < 1 || cols < 1 || n_bins < 1) { + throw std::invalid_argument( + "PixelHistogramImpl requires positive rows, cols and bins"); + } + if (!(xmax > xmin)) { + throw std::invalid_argument("PixelHistogramImpl requires xmax > xmin"); + } + + const T range = xmax - xmin; + for (int i = 0; i <= n_bins; ++i) { + m_edges(i) = + xmin + (static_cast(i) * range) / static_cast(n_bins); + } +} + +template +void PixelHistogramImpl::fill(const NDView &frame) { + if (frame.shape(0) != m_rows || frame.shape(1) != m_cols) { + throw std::invalid_argument( + "PixelHistogramImpl::fill: frame shape does not match histogram"); + } + for (int row = 0; row < m_rows; ++row) { + for (int col = 0; col < m_cols; ++col) { + fill_unchecked(row, col, frame(row, col)); + } + } +} + +template +void PixelHistogramImpl::fill(int row, int col, T value) { + if (row < 0 || row >= m_rows || col < 0 || col >= m_cols) { + throw std::out_of_range( + "PixelHistogramImpl::fill: row or col out of range"); + } + fill_unchecked(row, col, value); +} + +template +void PixelHistogramImpl::fill_unchecked(int row, int col, + T value) { + if (value < m_xmin || value >= m_xmax) { + return; + } + 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; + } + if constexpr (std::is_integral_v) { + if (m_values(row, col, bin) >= + std::numeric_limits::max()) { + return; + } + } + ++m_values(row, col, bin); +} + +template +NDArray PixelHistogramImpl::values() const { + return m_values; +} + +template +NDView PixelHistogramImpl::view() const { + return m_values.view(); +} + +template +NDArray PixelHistogramImpl::bin_centers() const { + NDArray centers({static_cast(m_n_bins)}); + for (int i = 0; i < m_n_bins; ++i) { + centers(i) = (m_edges(i) + m_edges(i + 1)) / T(2); + } + return centers; +} + +template +NDArray PixelHistogramImpl::bin_edges() const { + return m_edges; +} + +} // namespace aare diff --git a/python/aare/__init__.py b/python/aare/__init__.py index 3e4e64d..10872f8 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -33,7 +33,7 @@ from .CtbRawFile import CtbRawFile from .RawFile import RawFile from .ScanParameters import ScanParameters -from .utils import random_pixels, random_pixel, flat_list, add_colorbar +from .utils import random_pixels, random_pixel, flat_list, add_colorbar, Timer #make functions available in the top level API @@ -44,3 +44,13 @@ from ._aare import apply_calibration, count_switching_pixels from ._aare import calculate_pedestal, calculate_pedestal_float, calculate_pedestal_g0, calculate_pedestal_g0_float from ._aare import VarClusterFinder +from ._aare import ( + PedestalTrackingPixelHistogram, + PixelHistogram, + PixelHistogram_d, + PixelHistogram_f, + PixelHistogram_u8, + PixelHistogram_u16, + PixelHistogram_u32, + PixelHistogram_u64, +) diff --git a/python/aare/utils.py b/python/aare/utils.py index 2e5d2d3..98d2349 100644 --- a/python/aare/utils.py +++ b/python/aare/utils.py @@ -2,6 +2,7 @@ import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import make_axes_locatable +import time def random_pixels(n_pixels, xmin=0, xmax=512, ymin=0, ymax=1024): """Return a list of random pixels. @@ -34,4 +35,18 @@ def add_colorbar(ax, im, size="5%", pad=0.05): divider = make_axes_locatable(ax) cax = divider.append_axes("right", size=size, pad=pad) plt.colorbar(im, cax=cax) - return ax, im, cax \ No newline at end of file + return ax, im, cax + + +class Timer: + def __init__(self, label="Elapsed time:", verbose=True): + self.label = label + self.verbose = verbose + def __enter__(self): + self.start = time.perf_counter() + return self + + def __exit__(self, exc_type, exc, tb): + self.elapsed = time.perf_counter() - self.start + if self.verbose: + print(f"{self.label} {self.elapsed:.3f}s") \ No newline at end of file diff --git a/python/src/bind_PedestalTrackingPixelHistogram.hpp b/python/src/bind_PedestalTrackingPixelHistogram.hpp new file mode 100644 index 0000000..7984273 --- /dev/null +++ b/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 +#include +#include +#include + +namespace py = pybind11; +using namespace ::aare; + +void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) { + py::class_( + m, "PedestalTrackingPixelHistogram", + "A pixel-wise histogram of frame - pedestal residuals, with a " + "per-pixel running pedestal estimate sharded across worker threads") + .def( + py::init(), + 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 + 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()) + + .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 *ptr = + nullptr; + { + py::gil_scoped_release release; + ptr = new NDArray(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 + 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 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::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::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()) + + .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 *ptr = + nullptr; + { + py::gil_scoped_release release; + ptr = + new NDArray(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( + 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( + 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. + )"); +} diff --git a/python/src/bind_PixelHistogram.hpp b/python/src/bind_PixelHistogram.hpp new file mode 100644 index 0000000..d205bce --- /dev/null +++ b/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 +#include +#include +#include +#include +#include + +namespace py = pybind11; +using namespace ::aare; + +namespace { + +template +void define_pixel_histogram_binding(py::module &m, const char *class_name, + const char *storage_dtype) { + using Hist = PixelHistogram; + + const std::string doc = + std::string("A histogram for pixel-wise statistics with float64 input " + "axis and ") + + storage_dtype + " bin storage"; + + py::class_(m, class_name, doc.c_str()) + .def(py::init(), + 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 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 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()) + + .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 *ptr = nullptr; + { + py::gil_scoped_release release; + ptr = new NDArray(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(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(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(m, "PixelHistogram_d", "float64"); + define_pixel_histogram_binding(m, "PixelHistogram_f", "float32"); + define_pixel_histogram_binding(m, "PixelHistogram_u64", + "uint64"); + define_pixel_histogram_binding(m, "PixelHistogram_u32", + "uint32"); + define_pixel_histogram_binding(m, "PixelHistogram_u16", + "uint16"); + define_pixel_histogram_binding(m, "PixelHistogram_u8", + "uint8"); + + // Backwards-compatible alias for the generic Python class name. + m.attr("PixelHistogram") = m.attr("PixelHistogram_d"); +} diff --git a/python/src/module.cpp b/python/src/module.cpp index 19bbe29..023a19e 100644 --- a/python/src/module.cpp +++ b/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(m, "Pedestal_d"); define_pedestal_bindings(m, "Pedestal_f"); define_fit_bindings(m); diff --git a/python/src/pedestal.hpp b/python/src/pedestal.hpp index 9bef948..6f8fc06 100644 --- a/python/src/pedestal.hpp +++ b/python/src/pedestal.hpp @@ -21,6 +21,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 &>(); + auto v = self.view(); + std::array shape{ + static_cast(v.shape(0)), + static_cast(v.shape(1))}; + std::array byte_strides{ + static_cast(v.strides()[0]) * + static_cast(sizeof(SUM_TYPE)), + static_cast(v.strides()[1]) * + static_cast(sizeof(SUM_TYPE))}; + auto arr = py::array_t(shape, byte_strides, v.data(), + self_py); + arr.attr("setflags")(py::arg("write") = false); + return arr; + }) .def("variance", [](Pedestal &self) { auto var = new NDArray{}; @@ -49,6 +66,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 &pedestal, + py::array_t &f, + py::array_t &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 &pedestal, diff --git a/src/hist/PedestalTrackingPixelHistogram.cpp b/src/hist/PedestalTrackingPixelHistogram.cpp new file mode 100644 index 0000000..b444cef --- /dev/null +++ b/src/hist/PedestalTrackingPixelHistogram.cpp @@ -0,0 +1,584 @@ +#include "aare/hist/PedestalTrackingPixelHistogram.hpp" +#include "aare/File.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include +namespace aare { + +PedestalTrackingPixelHistogram::PedestalTrackingPixelHistogram( + int rows, int cols, int n_bins, AxisType xmin, AxisType xmax, int n_threads, + std::size_t max_pending, AxisType n_sigma) + : rows_(rows), cols_(cols), n_threads_(n_threads), xmin_(xmin), xmax_(xmax), + current_work_kind_(WorkKind::FillWithThreshold), current_image_(nullptr), + current_images_(nullptr), completed_threads_(0), stop_workers_(false), + work_generation_(0), max_batch_size_(max_pending / 3), n_sigma_(n_sigma) { + if (rows_ < 1 || cols_ < 1 || n_bins < 1) { + throw std::invalid_argument("PedestalTrackingPixelHistogram requires " + "positive rows, cols and bins"); + } + if (n_threads < 1) { + throw std::invalid_argument( + "PedestalTrackingPixelHistogram requires at least one thread"); + } + if (max_pending < 1) { + throw std::invalid_argument( + "PedestalTrackingPixelHistogram requires max_pending >= 1"); + } + + n_threads_ = std::min(n_threads_, rows_); + + // Build a balanced row partition. With base = rows_ / n_threads_ and + // extra = rows_ % n_threads_, the first `extra` threads get base + 1 + // rows each and the rest get `base` rows. This avoids the + // ceil(rows/n_threads) scheme leaving trailing threads with zero or + // negative row counts (e.g. rows=17, n_threads=8). + row_offsets_.resize(n_threads_ + 1); + const int base = rows_ / n_threads_; + const int extra = rows_ % n_threads_; + int offset = 0; + for (int i = 0; i < n_threads_; ++i) { + row_offsets_[i] = offset; + offset += base + (i < extra ? 1 : 0); + } + row_offsets_[n_threads_] = offset; // == rows_ by construction + + // Initialize partial histograms, partial pedestals and the cached + // per-pixel std for each thread. All three are sized to the + // thread's row slice and indexed by local_row (0..local_rows-1), + // so the worker can address them with the same coordinates. + partial_hists_.reserve(n_threads_); + partial_pedestals_.reserve(n_threads_); + partial_std_.reserve(n_threads_); + for (int i = 0; i < n_threads_; ++i) { + const auto local_rows = row_count(i); + partial_hists_.emplace_back(local_rows, cols, n_bins, xmin, xmax); + partial_pedestals_.emplace_back(static_cast(local_rows), + static_cast(cols)); + partial_std_.emplace_back(NDArray( + {static_cast(local_rows), static_cast(cols)}, + 0.0)); + } + + // Spawn worker threads + for (int i = 0; i < n_threads_; ++i) { + workers_.emplace_back([this, i]() { this->worker_loop(i); }); + } + + // Async pipeline. The PCQ holds (size - 1) usable slots, so size up by + // one to honour the requested max_pending. + async_queue_ = std::make_unique( + static_cast(max_pending + 1)); + coordinator_ = std::thread([this]() { this->coordinator_loop(); }); +} + +PedestalTrackingPixelHistogram::~PedestalTrackingPixelHistogram() { + // Drain any pending async fills before tearing down the worker pool. + // The coordinator's loop keeps processing while stop_coordinator_ is + // true as long as the queue is non-empty (mirrors ClusterFinderMT). + if (coordinator_.joinable()) { + stop_coordinator_ = true; + coordinator_.join(); + } + + // Signal all workers to stop + { + std::unique_lock lock(work_mutex_); + stop_workers_ = true; + } + work_cv_.notify_all(); + + // Join all worker threads + for (auto &thread : workers_) { + if (thread.joinable()) { + thread.join(); + } + } +} + +int PedestalTrackingPixelHistogram::row_start(int thread_id) const { + return row_offsets_[thread_id]; +} + +int PedestalTrackingPixelHistogram::row_count(int thread_id) const { + return row_offsets_[thread_id + 1] - row_offsets_[thread_id]; +} + +void PedestalTrackingPixelHistogram::dispatch_( + WorkKind kind, const NDView *image) { + // Caller must already hold fill_mutex_. Reset counters, publish the + // new work item, wake the workers, wait for completion. + { + std::unique_lock lock(work_mutex_); + completed_threads_ = 0; + current_work_kind_ = kind; + current_image_ = image; + current_images_ = nullptr; + ++work_generation_; + } + work_cv_.notify_all(); + { + std::unique_lock lock(work_mutex_); + done_cv_.wait(lock, + [this]() { return completed_threads_ == n_threads_; }); + current_image_ = nullptr; + current_images_ = nullptr; + } +} + +void PedestalTrackingPixelHistogram::dispatch_fill_batch_( + const std::vector> &images) { + // Caller must already hold fill_mutex_. `images` is owned by + // fill_with_threshold_batch_ and stays alive until this blocking dispatch + // returns. + { + std::unique_lock lock(work_mutex_); + completed_threads_ = 0; + current_work_kind_ = WorkKind::FillWithThreshold; + current_image_ = nullptr; + current_images_ = &images; + ++work_generation_; + } + work_cv_.notify_all(); + { + std::unique_lock lock(work_mutex_); + done_cv_.wait(lock, + [this]() { return completed_threads_ == n_threads_; }); + current_images_ = nullptr; + } +} + +void PedestalTrackingPixelHistogram::push_pedestal_no_update( + const NDView &frame) { + if (frame.shape(0) != rows_ || frame.shape(1) != cols_) { + throw std::invalid_argument( + "PedestalTrackingPixelHistogram frame shape does not match " + "constructor shape"); + } + std::lock_guard fill_lock(fill_mutex_); + dispatch_(WorkKind::PushPedestal, &frame); +} + +void PedestalTrackingPixelHistogram::update_mean() { + // Drain any in-flight async fills first; the coordinator does NOT + // hold fill_mutex_ at that point, so we can grab it safely after. + flush(); + std::lock_guard fill_lock(fill_mutex_); + dispatch_(WorkKind::UpdateMean, nullptr); +} + +void PedestalTrackingPixelHistogram::worker_loop(int thread_id) { + int last_generation = 0; + + while (true) { + std::unique_lock lock(work_mutex_); + work_cv_.wait(lock, [this, last_generation]() { + return work_generation_ != last_generation || stop_workers_; + }); + + if (stop_workers_) { + break; + } + + // Snapshot the work assignment under the lock so we don't race + // against the next dispatch publishing new state. + const WorkKind kind = current_work_kind_; + const NDView *image = current_image_; + const std::vector> *images = current_images_; + const int generation = work_generation_; + const int first_row = row_start(thread_id); + const int local_rows = row_count(thread_id); + + lock.unlock(); + + auto &my_pedestal = partial_pedestals_[thread_id]; + auto &my_hist = partial_hists_[thread_id]; + + switch (kind) { + case WorkKind::PushPedestal: { + // Accumulate raw frame values into this thread's pedestal + // shard. Uses the pixel-level push_no_update which only + // touches m_sum/m_sum2/m_cur_samples (no m_mean writes). + for (int local_row = 0; local_row < local_rows; ++local_row) { + const auto row = static_cast(first_row + local_row); + for (ssize_t col = 0; col < image->shape(1); ++col) { + my_pedestal.template push_no_update( + static_cast(local_row), + static_cast(col), (*image)(row, col)); + } + } + break; + } + case WorkKind::UpdateMean: { + // Recompute m_mean from the running sums. Only touches this + // thread's shard. Also refresh the cached per-pixel std so + // FillWithThreshold can read it without recomputing on the + // hot path. + my_pedestal.update_mean(); + auto &my_std = partial_std_[thread_id]; + for (int local_row = 0; local_row < local_rows; ++local_row) { + for (int col = 0; col < cols_; ++col) { + my_std(local_row, col) = static_cast( + my_pedestal.std(static_cast(local_row), + static_cast(col))); + } + } + break; + } + case WorkKind::FillWithThreshold: { + // Histogram the pedestal-subtracted residual AND, for pixels + // whose residual is consistent with noise + // (|residual| < n_sigma * cached_std), feed the raw value + // back into the pedestal shard. With n_sigma <= 0, use a + // histogram-only hot path that skips the per-pixel pedestal + // tracking gate entirely. The [xmin, xmax) histogram gate + // lives inside PixelHistogramImpl::fill. + const auto n_sigma = n_sigma_.load(std::memory_order_relaxed); + if (n_sigma <= AxisType{0.0}) { + // Fill without pedestal tracking. + for (const auto &frame : *images) { + const auto cols = frame.shape(1); + for (int local_row = 0; local_row < local_rows; + ++local_row) { + const auto row = + static_cast(first_row + local_row); + for (ssize_t col = 0; col < cols; ++col) { + const FrameType raw = frame(row, col); + const AxisType val = + static_cast(raw) - + static_cast(my_pedestal.mean( + static_cast(local_row), + static_cast(col))); + my_hist.fill_unchecked(local_row, + static_cast(col), val); + } + } + } + break; + } else { + // Do pedestal tracking. Duplicated code for clean hot path. + + auto &my_std = partial_std_[thread_id]; + for (const auto &frame : *images) { + const auto cols = frame.shape(1); + for (int local_row = 0; local_row < local_rows; + ++local_row) { + const auto row = + static_cast(first_row + local_row); + for (ssize_t col = 0; col < cols; ++col) { + const FrameType raw = frame(row, col); + const AxisType val = + static_cast(raw) - + static_cast(my_pedestal.mean( + static_cast(local_row), + static_cast(col))); + my_hist.fill_unchecked(local_row, + static_cast(col), val); + const AxisType sigma = my_std(local_row, col); + if (sigma > AxisType{0.0} && + std::abs(static_cast(val)) < + n_sigma * sigma) { + my_pedestal.template push( + static_cast(local_row), + static_cast(col), raw); + } + } + } + } + break; + } + } + } + + // Signal completion + { + std::unique_lock done_lock(work_mutex_); + last_generation = generation; + completed_threads_++; + if (completed_threads_ == n_threads_) { + done_cv_.notify_one(); + } + } + } +} + +NDArray +PedestalTrackingPixelHistogram::values() const { + // Make sure any pending async fills are merged in before we snapshot + // the partial histograms. Cheap when the queue is already drained. + flush(); + + const auto first_shard_view = partial_hists_.front().view(); + const auto cols = static_cast(first_shard_view.shape(1)); + const auto bins = static_cast(first_shard_view.shape(2)); + const auto rows = static_cast(rows_); + + NDArray data({rows, cols, bins}); + + // Each thread owns a disjoint, contiguous range of rows. The shard's + // dense storage layout [local_row][col][bin] is identical to the slice + // [first_row .. first_row + local_rows)[col][bin] of `data`, so the + // merge is just one bulk copy per thread; no per-element accumulation + // and no upfront zeroing of `data` is needed. + const size_t pixel_stride = static_cast(cols) * bins; + for (int t = 0; t < n_threads_; ++t) { + const auto first_row = static_cast(row_start(t)); + const auto local_rows = static_cast(row_count(t)); + if (local_rows == 0) + continue; + + const auto shard_view = partial_hists_[t].view(); + std::memcpy(data.data() + first_row * pixel_stride, shard_view.data(), + local_rows * pixel_stride * sizeof(StorageType)); + } + + return data; +} + +NDArray +PedestalTrackingPixelHistogram::pedestal_mean() const { + // Drain in-flight async fills and serialise with all other fan-outs + // (Fill / PushPedestal / UpdateMean). m_mean is overwritten wholesale + // by Pedestal::update_mean, so without the lock we could read torn + // rows mid-update. + flush(); + std::lock_guard lock(fill_mutex_); + + NDArray data( + {static_cast(rows_), static_cast(cols_)}); + + // Each partial pedestal stores its slice of m_mean in C-order + // [local_rows x cols], identical in layout to the corresponding + // [first_row .. first_row + local_rows)[col] slice of `data`, so + // we can copy each shard with a single memcpy. + const size_t row_stride = static_cast(cols_); + for (int t = 0; t < n_threads_; ++t) { + const auto first_row = static_cast(row_start(t)); + const auto local_rows = static_cast(row_count(t)); + if (local_rows == 0) + continue; + + const auto view = partial_pedestals_[t].view(); + std::memcpy(data.data() + first_row * row_stride, view.data(), + local_rows * row_stride * sizeof(AxisType)); + } + + return data; +} + +void PedestalTrackingPixelHistogram::fill_with_threshold_batch_( + std::vector> &batch) { + // Called only by the coordinator thread on images already shape-checked + // by fill_async, so there is no need to re-validate. fill_mutex_ is + // still required: push_pedestal_no_update / update_mean / pedestal_mean + // can run concurrently and must not race this fan-out. + std::vector> views; + views.reserve(batch.size()); + for (auto &frame : batch) { + views.push_back(frame.view()); + } + + std::lock_guard fill_lock(fill_mutex_); + dispatch_fill_batch_(views); +} + +void PedestalTrackingPixelHistogram::fill_async(NDArray &&image) { + if (image.shape(0) != rows_ || image.shape(1) != cols_) { + throw std::invalid_argument( + "PedestalTrackingPixelHistogram image shape does not match " + "constructor shape"); + } + + // SPSC backpressure: spin with a short sleep until a slot frees up. + // The std::move only consumes `image` on the iteration that succeeds + // (placement-new inside write() runs only when the slot is free). + while (!async_queue_->write(std::move(image))) { + std::this_thread::sleep_for(async_wait_); + } +} + +PedestalTrackingPixelHistogram::AxisType +PedestalTrackingPixelHistogram::n_sigma() const { + return n_sigma_.load(std::memory_order_relaxed); +} + +void PedestalTrackingPixelHistogram::set_n_sigma( + PedestalTrackingPixelHistogram::AxisType n_sigma) { + n_sigma_.store(n_sigma, std::memory_order_relaxed); +} + +void PedestalTrackingPixelHistogram::flush() const { + while (!async_queue_->isEmpty() || + coordinator_busy_.load(std::memory_order_acquire)) { + std::this_thread::sleep_for(async_wait_); + } +} + +void PedestalTrackingPixelHistogram::coordinator_loop() { + std::vector> batch; + batch.reserve(max_batch_size_); + + while (!stop_coordinator_.load(std::memory_order_acquire) || + !async_queue_->isEmpty()) { + batch.clear(); + + auto *first_item = async_queue_->frontPtr(); + if (first_item == nullptr) { + std::this_thread::sleep_for(async_wait_); + continue; + } + + coordinator_busy_.store(true, std::memory_order_release); + batch.push_back(std::move(*first_item)); + async_queue_->popFront(); + + while (batch.size() < max_batch_size_) { + auto *item = async_queue_->frontPtr(); + if (item == nullptr) { + break; + } + batch.push_back(std::move(*item)); + async_queue_->popFront(); + } + + fill_with_threshold_batch_(batch); + completed_async_fills_.fetch_add(batch.size(), + std::memory_order_release); + coordinator_busy_.store(false, std::memory_order_release); + } +} + +NDArray +PedestalTrackingPixelHistogram::bin_centers() const { + // All shards share the same value-axis configuration, so any one will + // do; pick the first. + return partial_hists_.front().bin_centers(); +} + +NDArray +PedestalTrackingPixelHistogram::bin_edges() const { + return partial_hists_.front().bin_edges(); +} + +void PedestalTrackingPixelHistogram::fill_from_file( + const std::filesystem::path &fname, ssize_t max_frames, bool verbose) { + constexpr std::size_t progress_interval = 66; + auto last = std::chrono::steady_clock::now(); + const auto completed_start = + completed_async_fills_.load(std::memory_order_acquire); + std::size_t last_reported = 0; + const auto completed_for_this_file = [&]() { + return completed_async_fills_.load(std::memory_order_acquire) - + completed_start; + }; + + const auto wait_for_completed = [&](std::size_t target) { + while (completed_for_this_file() < target) { + std::this_thread::sleep_for(async_wait_); + } + }; + + File f(fname); + // check that row col matches constructor + if (f.rows() != static_cast(rows_) || + f.cols() != static_cast(cols_)) { + throw std::invalid_argument("PedestalTrackingPixelHistogram: Frame in " + "file {} has shape ({}, {}) does not match " + "constructor shape"); + } + + const ssize_t total_frames = f.total_frames(); + const ssize_t n_frames = + max_frames == -1 ? total_frames : std::min(max_frames, total_frames); + const auto print_progress = [&](std::size_t done) { + const auto now = std::chrono::steady_clock::now(); + const double dt = std::chrono::duration(now - last).count(); + const auto done_in_interval = done - last_reported; + const double fps = + dt > 0.0 ? static_cast(done_in_interval) / dt : 0.0; + + fmt::print( + "\rProgress: {}/{} ({:.1f}%) {:.1f} FPS ", done, n_frames, + 100.0 * static_cast(done) / static_cast(n_frames), + fps); + std::fflush(stdout); + last = now; + last_reported = done; + }; + + for (ssize_t i = 0; i < n_frames; ++i) { + aare::NDArray frame({rows_, cols_}); + f.read_into(reinterpret_cast(frame.data())); + fill_async(std::move(frame)); + + if (verbose && (i + 1) % progress_interval == 0) { + wait_for_completed(static_cast(i + 1)); + print_progress(static_cast(i + 1)); + } + } + flush(); + if (verbose) { + const auto done = completed_for_this_file(); + if (done > last_reported) { + print_progress(done); + } + fmt::print("\n\n"); + std::fflush(stdout); + } +} + +void PedestalTrackingPixelHistogram::process_pedestal_file( + const std::filesystem::path &fname, ssize_t max_frames, bool verbose) { + constexpr std::size_t progress_interval = 66; + auto last = std::chrono::steady_clock::now(); + + File f(fname); + // check that row col matches constructor + if (f.rows() != static_cast(rows_) || + f.cols() != static_cast(cols_)) { + throw std::invalid_argument("PedestalTrackingPixelHistogram: Frame in " + "file {} has shape ({}, {}) does not match " + "constructor shape"); + } + + const ssize_t total_frames = f.total_frames(); + const ssize_t n_frames = + max_frames == -1 ? total_frames : std::min(max_frames, total_frames); + + aare::NDArray frame({rows_, cols_}); + for (ssize_t i = 0; i < n_frames; ++i) { + f.read_into(reinterpret_cast(frame.data())); + push_pedestal_no_update(frame.view()); + if (verbose && + ((i + 1) % progress_interval == 0 || (i + 1 == n_frames))) { + const auto now = std::chrono::steady_clock::now(); + const double dt = std::chrono::duration(now - last).count(); + const std::size_t done_in_interval = + (i + 1) % progress_interval == 0 ? progress_interval + : (i + 1) % progress_interval; + const double fps = + dt > 0.0 ? static_cast(done_in_interval) / dt : 0.0; + fmt::print("\rProgress: {}/{} ({:.1f}%) {:.1f} FPS ", i + 1, + n_frames, + 100.0 * static_cast(i + 1) / + static_cast(n_frames), + fps); + std::fflush(stdout); + last = now; + } + } + update_mean(); + flush(); + if (verbose) { + fmt::print("\n\n"); + std::fflush(stdout); + } +} + +} // namespace aare diff --git a/src/hist/PixelHistogram.test.cpp b/src/hist/PixelHistogram.test.cpp new file mode 100644 index 0000000..db2f9ca --- /dev/null +++ b/src/hist/PixelHistogram.test.cpp @@ -0,0 +1,319 @@ +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "test_config.hpp" +#include "test_macros.hpp" + +#include "aare/hist/PixelHistogram.hpp" + +using aare::NDArray; +using aare::NDView; +using aare::PixelHistogram; + +namespace { +// The synchronous fill() has been removed; fill_async() is the only entry +// point. This helper submits one frame and blocks until it has been merged +// so the tests can keep their straightforward, ordered expectations. + +} // namespace + +TEST_CASE("Fill one pixel of a 5x10 histogram") { + PixelHistogram hist(5, 10, 20, 0.0, 10.0); + NDArray image( + {5, 10}, -1.0); // Need to fill with -1 to not generate counts + + image(2, 3) = 5.7; // This should go into bin 11 (since bins are [0-0.5), + // [0.5-1.0), ..., [9.5-10.0)) + + // fill_blocking(hist, image.view()); + hist.fill_async(NDArray(image)); + hist.flush(); // Wait for the async fill to complete before we check the + // results + + auto values = hist.values(); + REQUIRE(values.shape(0) == 5); + REQUIRE(values.shape(1) == 10); + REQUIRE(values.shape(2) == 20); + + // Check that the correct bin for pixel (2,3) has count 1 + CHECK(values(2, 3, 11) == 1); + + // Check that all other bins are zero + for (ssize_t row = 0; row < values.shape(0); ++row) { + for (ssize_t col = 0; col < values.shape(1); ++col) { + for (ssize_t bin = 0; bin < values.shape(2); ++bin) { + if (!(row == 2 && col == 3 && bin == 11)) { + CHECK(values(row, col, bin) == 0); + } + } + } + } +} + +TEST_CASE("Fill pixels with uneven partial histogram row slices") { + PixelHistogram hist(5, 4, 10, 0.0, 10.0, 3); + NDArray image({5, 4}, -1.0); + + image(0, 0) = 0.2; + image(1, 1) = 1.2; + image(2, 2) = 2.2; + image(3, 3) = 3.2; + image(4, 0) = 4.2; + + hist.fill_async(NDArray(image)); + hist.flush(); + + auto values = hist.values(); + REQUIRE(values.shape(0) == 5); + REQUIRE(values.shape(1) == 4); + REQUIRE(values.shape(2) == 10); + + CHECK(values(0, 0, 0) == 1); + CHECK(values(1, 1, 1) == 1); + CHECK(values(2, 2, 2) == 1); + CHECK(values(3, 3, 3) == 1); + CHECK(values(4, 0, 4) == 1); + + for (ssize_t row = 0; row < values.shape(0); ++row) { + for (ssize_t col = 0; col < values.shape(1); ++col) { + for (ssize_t bin = 0; bin < values.shape(2); ++bin) { + const bool expected = (row == 0 && col == 0 && bin == 0) || + (row == 1 && col == 1 && bin == 1) || + (row == 2 && col == 2 && bin == 2) || + (row == 3 && col == 3 && bin == 3) || + (row == 4 && col == 0 && bin == 4); + if (!expected) { + CHECK(values(row, col, bin) == 0); + } + } + } + } +} + +TEST_CASE("Row partitioning handles rows < n_threads * ceil(rows/n_threads)") { + // Regression test for the pre-existing row partitioning bug: with the + // old ceil(rows / n_threads) scheme, rows=17 and n_threads=8 left + // trailing threads with negative row counts and threw + // "stop >= start required" from boost::histogram during construction. + // The current scheme distributes the remainder so every thread gets + // at least floor(rows / n_threads) rows. We also exercise the case + // where n_threads > rows, which is clamped down to rows. + const auto p = GENERATE(table({ + {17, 8}, // previously broken: 8 threads * ceil(17/8) = 24 > 17 + {17, 5}, // previously broken: 5 threads * ceil(17/5) = 20 > 17 + {13, 4}, // previously broken: 4 * ceil(13/4) = 16 > 13 + {17, 32}, // n_threads > rows -> clamped to rows + {1, 8}, // n_threads > rows -> single thread, single row + {6, 6}, // balanced + })); + const int rows = std::get<0>(p); + const int n_threads = std::get<1>(p); + const int cols = 3; + const int n_bins = 4; + constexpr float xmin = 0.0f; + constexpr float xmax = 4.0f; + + PixelHistogram hist(rows, cols, n_bins, xmin, xmax, n_threads); + + // Put one in-range value per row so we can verify every row is covered. + NDArray img({rows, cols}, -1.0f); // -1 is out of range + for (ssize_t r = 0; r < rows; ++r) { + img(r, 0) = static_cast(r % n_bins) + 0.5f; + } + hist.fill_async(NDArray(img)); + hist.flush(); + + auto h = hist.values(); + REQUIRE(h.shape(0) == rows); + REQUIRE(h.shape(1) == cols); + REQUIRE(h.shape(2) == n_bins); + + for (ssize_t r = 0; r < rows; ++r) { + const int expected_bin = static_cast(r % n_bins); + for (ssize_t c = 0; c < cols; ++c) { + for (ssize_t b = 0; b < n_bins; ++b) { + const uint16_t want = + (c == 0 && b == expected_bin) ? uint16_t{1} : uint16_t{0}; + if (h(r, c, b) != want) { + INFO("rows=" << rows << " n_threads=" << n_threads + << " r=" << r << " c=" << c << " b=" << b + << " got=" << h(r, c, b) << " want=" << want); + CHECK(h(r, c, b) == want); + } + } + } + } +} + +TEST_CASE("Random fills match a reference implementation") { + // End-to-end correctness check: compare the implementation against a + // simple per-pixel reference for several thread counts. Values are + // sampled slightly wider than [xmin, xmax) so the out-of-range filter + // is also exercised. + constexpr int rows = 17; + constexpr int cols = 23; + constexpr int n_bins = 32; + constexpr float xmin = -1.5f; + constexpr float xmax = 4.5f; + const int n_threads = GENERATE(1, 2, 4, 8); + + PixelHistogram hist(rows, cols, n_bins, xmin, xmax, n_threads); + + std::mt19937 rng(0xC0FFEE); + std::uniform_real_distribution dist(xmin - 0.5f, xmax + 0.5f); + + constexpr int n_frames = 4; + std::vector> frames; + frames.reserve(n_frames); + for (int f = 0; f < n_frames; ++f) { + NDArray img({rows, cols}, 0.0f); + for (ssize_t r = 0; r < rows; ++r) { + for (ssize_t c = 0; c < cols; ++c) { + img(r, c) = dist(rng); + } + } + frames.push_back(std::move(img)); + } + + NDArray expected({rows, cols, n_bins}, 0u); + const float inv_range = static_cast(n_bins) / (xmax - xmin); + for (const auto &img : frames) { + for (ssize_t r = 0; r < rows; ++r) { + for (ssize_t c = 0; c < cols; ++c) { + const float v = img(r, c); + if (!(v >= xmin) || !(v < xmax)) + continue; + int bin = static_cast((v - xmin) * inv_range); + if (bin >= n_bins) + bin = n_bins - 1; + if (bin < 0) + bin = 0; + expected(r, c, bin) += 1; + } + } + } + + for (const auto &img : frames) { + hist.fill_async(NDArray(img)); + } + hist.flush(); + auto h = hist.values(); + + REQUIRE(h.shape(0) == rows); + REQUIRE(h.shape(1) == cols); + REQUIRE(h.shape(2) == n_bins); + + bool all_match = true; + for (ssize_t r = 0; r < rows && all_match; ++r) { + for (ssize_t c = 0; c < cols && all_match; ++c) { + for (ssize_t b = 0; b < n_bins && all_match; ++b) { + if (h(r, c, b) != expected(r, c, b)) { + all_match = false; + INFO("n_threads=" << n_threads << " r=" << r << " c=" << c + << " b=" << b << " got=" << h(r, c, b) + << " expected=" << expected(r, c, b)); + CHECK(h(r, c, b) == expected(r, c, b)); + } + } + } + } + CHECK(all_match); +} + +TEST_CASE("fill_async with mismatched shape throws") { + PixelHistogram hist(8, 8, 16, 0.0, 1.0, 2); + NDArray bad({4, 4}, 0.0f); + CHECK_THROWS_AS(hist.fill_async(std::move(bad)), std::invalid_argument); +} + +TEST_CASE("Destructor drains pending async fills") { + // Submit more frames than the queue can hold so backpressure kicks in, + // then immediately let the histogram go out of scope and verify that + // the merged values() matches the reference computed sequentially. + constexpr int rows = 11; + constexpr int cols = 7; + constexpr int n_bins = 8; + constexpr float xmin = 0.0f; + constexpr float xmax = 1.0f; + constexpr int n_frames = 12; + constexpr std::size_t max_pending = 2; + + std::vector> frames; + frames.reserve(n_frames); + std::mt19937 rng(0xDEADBEEF); + std::uniform_real_distribution dist(xmin, xmax); + for (int f = 0; f < n_frames; ++f) { + NDArray img({rows, cols}, 0.0f); + for (ssize_t r = 0; r < rows; ++r) + for (ssize_t c = 0; c < cols; ++c) + img(r, c) = dist(rng); + frames.push_back(std::move(img)); + } + + NDArray snapshot({rows, cols, n_bins}, uint16_t{0}); + { + PixelHistogram hist(rows, cols, n_bins, xmin, xmax, 2, max_pending); + for (auto &img : frames) { + // Move a copy so we can also build the reference below. + NDArray copy({rows, cols}, 0.0f); + std::memcpy(copy.data(), img.data(), copy.total_bytes()); + hist.fill_async(std::move(copy)); + } + // No explicit flush(); destructor must drain. + // Capture values() *after* the loop but inside the scope so it + // observes everything that was submitted (values flushes too). + snapshot = hist.values(); + } + + PixelHistogram reference(rows, cols, n_bins, xmin, xmax, 2); + for (const auto &img : frames) + reference.fill_async(NDArray(img.view())); + auto expected = reference.values(); + + bool all_match = true; + for (ssize_t r = 0; r < rows && all_match; ++r) { + for (ssize_t c = 0; c < cols && all_match; ++c) { + for (ssize_t b = 0; b < n_bins && all_match; ++b) { + if (snapshot(r, c, b) != expected(r, c, b)) { + all_match = false; + INFO("r=" << r << " c=" << c << " b=" << b + << " got=" << snapshot(r, c, b) + << " expected=" << expected(r, c, b)); + CHECK(snapshot(r, c, b) == expected(r, c, b)); + } + } + } + } + CHECK(all_match); +} + +TEST_CASE("PixelHistogram supports custom storage and axis types") { + PixelHistogram hist(2, 2, 4, 0.0, 4.0, 1); + NDArray image({2, 2}, -1.0); + + image(0, 0) = 0.25; + image(0, 1) = 1.25; + image(1, 0) = 2.25; + image(1, 1) = 3.25; + + hist.fill_async(NDArray(image)); + hist.flush(); + + auto values = hist.values(); + STATIC_REQUIRE(std::is_same_v); + CHECK(values(0, 0, 0) == uint32_t{1}); + CHECK(values(0, 1, 1) == uint32_t{1}); + CHECK(values(1, 0, 2) == uint32_t{1}); + CHECK(values(1, 1, 3) == uint32_t{1}); + + auto centers = hist.bin_centers(); + STATIC_REQUIRE(std::is_same_v); + CHECK(centers.shape(0) == 4); +} diff --git a/src/hist/PixelHistogramImpl.test.cpp b/src/hist/PixelHistogramImpl.test.cpp new file mode 100644 index 0000000..3cc2842 --- /dev/null +++ b/src/hist/PixelHistogramImpl.test.cpp @@ -0,0 +1,140 @@ +#include +#include + +#include "aare/NDArray.hpp" +#include "aare/hist/PixelHistogramImpl.hpp" + +TEST_CASE("PixelHistogramImpl construct a small histogram and fill with a few " + "values") { + int rows = 3; + int cols = 5; + int n_bins = 10; + double xmin = 0.0; + double xmax = 1.0; + aare::PixelHistogramImpl hist(rows, cols, n_bins, xmin, + xmax); + + // Check that the histogram is initialized correctly + REQUIRE(hist.view().shape(0) == rows); + REQUIRE(hist.view().shape(1) == cols); + REQUIRE(hist.view().shape(2) == n_bins); + + auto v = hist.view(); + for (int row = 0; row < rows; ++row) { + for (int col = 0; col < cols; ++col) { + for (int bin = 0; bin < n_bins; ++bin) { + REQUIRE(v(row, col, bin) == 0); + } + } + } + + auto edges = hist.bin_edges(); + for (int i = 0; i <= n_bins; ++i) { + REQUIRE_THAT(edges(i), Catch::Matchers::WithinRel(0.1 * i, 1e-12)); + } + + hist.fill(1, 1, 0.05); + REQUIRE(v(1, 1, 0) == 1); + + hist.fill(1, 1, 0.0999); + REQUIRE(v(1, 1, 0) == 2); + + hist.fill(1, 1, 0.1); + REQUIRE( + v(1, 1, 1) == + 1); // At the edge should go to the bin with the edge as the lower bound +} + +TEST_CASE("Fill a small histogram from an NDArray") { + int rows = 2; + int cols = 3; + int n_bins = 10; + double xmin = 0.0; + double xmax = 1.0; + aare::PixelHistogramImpl hist(rows, cols, n_bins, xmin, + xmax); + + aare::NDArray frame({rows, cols}); + frame(0, 0) = 0.05; // bin 0 + frame(0, 1) = 0.15; // bin 1 + frame(0, 2) = 0.25; // bin 2 + frame(1, 0) = 0.35; // bin 3 + frame(1, 1) = 0.45; // bin 4 + frame(1, 2) = 1.5; // out of range + + hist.fill(frame.view()); + + auto v = hist.view(); + REQUIRE(v(0, 0, 0) == 1); + REQUIRE(v(0, 1, 1) == 1); + REQUIRE(v(0, 2, 2) == 1); + REQUIRE(v(1, 0, 3) == 1); + REQUIRE(v(1, 1, 4) == 1); +} + +TEST_CASE("Check that pixel histogram does not overflow") { + int rows = 1; + int cols = 1; + int n_bins = 10; + double xmin = 0.0; + double xmax = 1.0; + aare::PixelHistogramImpl hist_u8(rows, cols, n_bins, xmin, + xmax); + + for (int i = 0; i < 255; ++i) { + hist_u8.fill(0, 0, 0.05); + } + + auto v0 = hist_u8.view(); + REQUIRE(v0(0, 0, 0) == 255); + + hist_u8.fill(0, 0, 0.05); + REQUIRE(v0(0, 0, 0) == 255); + + aare::PixelHistogramImpl hist_i8(rows, cols, n_bins, xmin, + xmax); + + for (int i = 0; i < 350; ++i) { + hist_i8.fill(0, 0, 0.05); + } + + auto v1 = hist_i8.view(); + REQUIRE(v1(0, 0, 0) == 127); +} + +TEST_CASE("Check that values outside the range do not affect the histogram") { + int rows = 1; + int cols = 1; + int n_bins = 10; + double xmin = 0.0; + double xmax = 1.0; + aare::PixelHistogramImpl hist(rows, cols, n_bins, xmin, + xmax); + hist.fill(0, 0, -0.1); // below range + hist.fill(0, 0, 1.0); // at upper edge (should be out of range) + hist.fill(0, 0, 1.1); // above range + + auto v = hist.view(); + auto total = std::accumulate(v.begin(), v.end(), 0); + REQUIRE(total == 0); +} + +TEST_CASE("Check that row and column bounds are checked") { + int rows = 1; + int cols = 1; + int n_bins = 10; + double xmin = 0.0; + double xmax = 1.0; + aare::PixelHistogramImpl hist(rows, cols, n_bins, xmin, + xmax); + REQUIRE_THROWS_AS(hist.fill(0, 1, -0.1), + std::out_of_range); // col out of range + REQUIRE_THROWS_AS(hist.fill(1, 0, 1.0), + std::out_of_range); // row out of range + REQUIRE_THROWS_AS(hist.fill(58, -1, 1.1), + std::out_of_range); // both out of range + + auto v = hist.view(); + auto total = std::accumulate(v.begin(), v.end(), 0); + REQUIRE(total == 0); +}