diff --git a/CMakeLists.txt b/CMakeLists.txt index e044b3f8..5cb45ea8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,6 +55,10 @@ option( "Install the python extension in the install tree under CMAKE_INSTALL_PREFIX/aare/" OFF) option(AARE_ASAN "Enable AddressSanitizer" OFF) +option( + AARE_TUNE_LOCAL + "Optimize for the building machine's CPU (-march=native -mtune=native). Not portable, the resulting binaries might not run on other machines." + OFF) # Configure which of the dependencies to use FetchContent for option(AARE_FETCH_FMT "Use FetchContent to download fmt" ON) @@ -296,6 +300,22 @@ else() target_compile_options(aare_compiler_flags INTERFACE -Werror) endif() + if(AARE_TUNE_LOCAL) + include(CheckCXXCompilerFlag) + check_cxx_compiler_flag("-march=native" AARE_HAS_MARCH_NATIVE) + check_cxx_compiler_flag("-mtune=native" AARE_HAS_MTUNE_NATIVE) + if(AARE_HAS_MARCH_NATIVE AND AARE_HAS_MTUNE_NATIVE) + message(STATUS "Tuning for the local CPU: -march=native -mtune=native") + target_compile_options(aare_compiler_flags INTERFACE -march=native + -mtune=native) + else() + message( + WARNING + "AARE_TUNE_LOCAL requested but the compiler does not support -march=native/-mtune=native. Ignoring." + ) + endif() + endif() + endif() # GCC/Clang specific if(AARE_PYTHON_BINDINGS) @@ -339,6 +359,7 @@ set(PUBLICHEADERS include/aare/Models.hpp include/aare/FileInterface.hpp include/aare/FilePtr.hpp + include/aare/FastPedestal.hpp include/aare/Frame.hpp include/aare/MultiThreadedFileReader.hpp include/aare/hist/PixelHistogram.hpp @@ -449,6 +470,7 @@ if(AARE_TESTS) ${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/hist/PedestalTrackingPixelHistogram.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/MultiThreadedFileReader.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp diff --git a/RELEASE.md b/RELEASE.md index 638b8b3b..f4ed7745 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -2,8 +2,29 @@ ## Next +### New Features: + +- Added ``FastPedestal`` in C++ and Python for per-pixel running mean, + population variance, and standard deviation. It supports exponentially + weighted updates, initialization from files, direct subtraction from NumPy + arrays, and ``float64``, ``float32``, and ``int16`` output types. +- Added the ``Pedestal_i16`` Python binding alongside + ``FastPedestal_d``, ``FastPedestal_f``, and ``FastPedestal_i16``. +- ``PedestalTrackingPixelHistogram.fill_from_file()`` now uses parallel, + double-buffered file reading and accepts ``reader_threads`` and + ``reader_chunk_size`` tuning parameters. +- Added the ``AARE_TUNE_LOCAL`` CMake option to build with ``-march=native`` + and ``-mtune=native`` when supported. Binaries built with this option are + specific to the local CPU and may not be portable. + ### API Changes: +- ``ClusterFinder`` now uses ``FastPedestal``. It must receive 1000 pedestal + frames before cluster finding; ``find_clusters()`` raises + an error until initialization is complete. Added ``update_threshold()`` to + recompute the per-pixel detection thresholds. +- Added the ``queue_depth`` constructor argument to ``ClusterFinderMT`` to + configure the number of preallocated frame buffers per worker thread. - Exposed ``ClusterVector.empty()`` in the Python API. - Added ClusterVector.estimate_n_clusters - Removed the lmfit dependency and the legacy ``fit_gaus``, ``fit_pol1``, @@ -167,5 +188,3 @@ dhanya.thattil@psi.ch - - diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index ef378e4d..78031277 100644 --- a/benchmarks/CMakeLists.txt +++ b/benchmarks/CMakeLists.txt @@ -18,8 +18,8 @@ FetchContent_MakeAvailable(benchmark) add_executable(benchmarks) target_sources( - benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp - reduce_benchmark.cpp) + benchmarks PRIVATE ndarray_benchmark.cpp ndview_benchmark.cpp + calculateeta_benchmark.cpp reduce_benchmark.cpp) # Link Google Benchmark and other necessary libraries target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core diff --git a/benchmarks/ndview_benchmark.cpp b/benchmarks/ndview_benchmark.cpp new file mode 100644 index 00000000..b0653193 --- /dev/null +++ b/benchmarks/ndview_benchmark.cpp @@ -0,0 +1,17 @@ +// SPDX-License-Identifier: MPL-2.0 +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" +#include + +using aare::NDArray; +using aare::NDView; + +static void BM_CreateNDView(benchmark::State &st) { + NDArray arr{{1024, 1024}, 0}; + for (auto _ : st) { + // This code gets timed + auto res = arr.view(); + benchmark::DoNotOptimize(res); + } +} +BENCHMARK(BM_CreateNDView); \ No newline at end of file diff --git a/docs/src/index.rst b/docs/src/index.rst index d60fc78e..ab22e57d 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -30,6 +30,7 @@ AARE python/file/index python/experimental/index python/histogram/index + python/pedestal/index pyFit diff --git a/docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst b/docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst index 18d438cd..66fe092d 100644 --- a/docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst +++ b/docs/src/python/histogram/pyPedestalTrackingPixelHistogram.rst @@ -15,6 +15,14 @@ Use ``push_pedestal_no_update()`` to seed the pedestal estimate, then asynchronous fills are drained by ``flush()``, and snapshot methods such as ``values()`` and ``pedestal_mean()`` return numpy arrays. +``fill_from_file()`` uses parallel file-reader workers and a double-buffered +pipeline. After the initial batch has been read, histogram processing of one +batch overlaps reading of the next. ``reader_threads`` and +``reader_chunk_size`` tune the I/O stage independently of the histogram worker +count. Two fixed-capacity buffers are allocated once and reused by alternating +their read and histogram roles. Their approximate memory use is ``2 * +reader_threads * reader_chunk_size * rows * cols * sizeof(uint16)``. + .. py:currentmodule:: aare .. autoclass:: PedestalTrackingPixelHistogram diff --git a/docs/src/python/pedestal/index.rst b/docs/src/python/pedestal/index.rst new file mode 100644 index 00000000..4e0a05e4 --- /dev/null +++ b/docs/src/python/pedestal/index.rst @@ -0,0 +1,9 @@ +Pedestal +======== + +.. toctree:: + :caption: Pedestal + :maxdepth: 1 + + pyFastPedestal + pyPedestal diff --git a/docs/src/python/pedestal/pyFastPedestal.rst b/docs/src/python/pedestal/pyFastPedestal.rst new file mode 100644 index 00000000..d442d298 --- /dev/null +++ b/docs/src/python/pedestal/pyFastPedestal.rst @@ -0,0 +1,83 @@ +FastPedestal +============ + +``FastPedestal`` calculates a running mean, variance and standard deviation for each pixel in a +series of frames. The python binding only exposes ``uint16`` input but the underlying +C++ class is templated. Initialize it with ``n_samples`` frames using +``add_init_frame()``. Once ``ready`` is true, use ``push_ema()`` to update the exponential +moving average initialized by the mean and with smoothing factor 1/n_samples. + +.. warning:: + + FastPedestal is not usable until you have added ``n_samples`` initial frames with ``add_init_frame(raw)``. + You can check the state with ``ready``. + +The public factory selects the bound C++ specialization from ``dtype``: + +* ``numpy.float64`` creates ``FastPedestal_d`` +* ``numpy.float32`` creates ``FastPedestal_f`` +* ``numpy.int16`` creates ``FastPedestal_i16`` + +The internal calculations are done with double, but the cached mean and on demand var and std are returned in the specified type. + +Factory +------- + +.. py:currentmodule:: aare + +.. autofunction:: FastPedestal + +Loading from a file +------------------- + +``FastPedestal.from_file()`` initializes the pedestal from ``n_samples`` +frames after ``skip_first``, then applies steady-state updates for any frames +remaining in the file. The input frames must contain ``uint16`` data; ``dtype`` +selects the output type of the pedestal statistics. + +.. autofunction:: aare.FastPedestal.from_file + +.. code-block:: python + + pedestal = FastPedestal.from_file( + "frames.npy", n_samples=100, skip_first=10, dtype=np.float32 + ) + +Example +------- + +.. code-block:: python + + import numpy as np + from aare import FastPedestal + + pedestal = FastPedestal(512, 1024, n_samples=100, dtype=np.float32) + + # Initialize with n_samples frames + for frame in initialization_frames: + pedestal.add_init_frame(frame) + + # Now we can push a frame for pedestal update + if pedestal.ready: + pedestal.push_ema(next_frame) + + # Mean and std are also ready + mean = pedestal.mean() + noise = pedestal.std() + + # Direct pedestal subtraction is also supported + for frame in raw_data: + image = frame - pedestal + +Complete API +------------ + +The API below is for the ``float64`` specialization. All dtype variants share +the same API. + +.. autoclass:: aare._aare.FastPedestal_d + :special-members: __init__ + :members: + :undoc-members: + :show-inheritance: + :inherited-members: diff --git a/docs/src/python/pedestal/pyPedestal.rst b/docs/src/python/pedestal/pyPedestal.rst new file mode 100644 index 00000000..7d6aecb2 --- /dev/null +++ b/docs/src/python/pedestal/pyPedestal.rst @@ -0,0 +1,42 @@ +Pedestal +======== + +``Pedestal`` calculates a running mean and variance for each pixel in a series +of ``uint16`` frames. ``push()`` updates the cached mean immediately. For +faster batch initialization, use ``push_no_update()`` for each frame and call +``update_mean()`` after the batch. + +Three specializations are available from :mod:`aare`: + +* ``Pedestal_d`` uses ``float64`` storage +* ``Pedestal_f`` uses ``float32`` storage +* ``Pedestal_i16`` uses ``int16`` storage + +Example +------- + +.. code-block:: python + + from aare import Pedestal_d + + pedestal = Pedestal_d(512, 1024, 100) + + for frame in initialization_frames: + pedestal.push_no_update(frame) + + pedestal.update_mean() + mean = pedestal.mean() + noise = pedestal.std() + +Complete API +------------ + +The API below is for the ``float64`` specialization. All dtype variants share +the same API. + +.. autoclass:: aare._aare.Pedestal_d + :special-members: __init__ + :members: + :undoc-members: + :show-inheritance: + :inherited-members: diff --git a/include/aare/Backoff.hpp b/include/aare/Backoff.hpp new file mode 100644 index 00000000..617c8aae --- /dev/null +++ b/include/aare/Backoff.hpp @@ -0,0 +1,48 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once + +#include +#include + +namespace aare { + +/** + * @brief Hint the CPU that this is a spin-wait. Prefer this over yield() when + * the wait is expected to be short. + */ +inline void cpu_relax() noexcept { +#if defined(__x86_64__) || defined(__i386__) + __builtin_ia32_pause(); +#elif defined(__aarch64__) + asm volatile("yield" ::: "memory"); +#else + std::this_thread::yield(); +#endif +} + +/** + * @brief Escalating wait for producer/consumer idle and backpressure loops. + * + * Starts with a pause instruction (sub-microsecond, no syscall), then yield, + * then a short sleep. Call reset() whenever work arrives so a busy pipeline + * never leaves the pause tier. + */ +class Backoff { + int m_count{0}; + + public: + void reset() noexcept { m_count = 0; } + + void pause() noexcept { + if (m_count < 64) { + cpu_relax(); + } else if (m_count < 256) { + std::this_thread::yield(); + } else { + std::this_thread::sleep_for(std::chrono::microseconds(50)); + } + ++m_count; + } +}; + +} // namespace aare diff --git a/include/aare/CircularFifo.hpp b/include/aare/CircularFifo.hpp index e0f4842b..fb275b72 100644 --- a/include/aare/CircularFifo.hpp +++ b/include/aare/CircularFifo.hpp @@ -1,11 +1,13 @@ // SPDX-License-Identifier: MPL-2.0 #pragma once +#include #include #include #include #include #include +#include #include "aare/ProducerConsumerQueue.hpp" @@ -29,6 +31,23 @@ template class CircularFifo { } } + /** + * @brief Construct a fifo and seed the free list using a factory. + * @param size number of items circulating in the fifo + * @param make_item callable invoked as make_item(i) for each slot index i + * + * Use this instead of CircularFifo(size) when the items need to be + * initialized, for example to hold a preallocated buffer or to carry + * their own slot index. + */ + template + CircularFifo(uint32_t size, F make_item) + : fifo_size(size), free_slots(size + 1), filled_slots(size + 1) { + for (size_t i = 0; i < fifo_size; ++i) { + free_slots.write(make_item(i)); + } + } + bool next() { // TODO! avoid default constructing ItemType ItemType it; @@ -47,6 +66,13 @@ template class CircularFifo { auto numFreeSlots() const noexcept { return free_slots.sizeGuess(); } auto isFull() const noexcept { return filled_slots.isFull(); } + /** + * @brief True if there are no filled slots waiting to be consumed. + * @note Prefer this over numFilledSlots() == 0 since sizeGuess() may + * under-report when called from the producing thread. + */ + auto isEmpty() const noexcept { return filled_slots.isEmpty(); } + ItemType pop_free() { ItemType v; while (!free_slots.read(v)) @@ -75,8 +101,22 @@ template class CircularFifo { ItemType *frontPtr() { return filled_slots.frontPtr(); } - // TODO! Add function to move item from filled to free to be used - // with the frontPtr function + /** + * @brief Return the front filled item to the free list. To be used + * together with frontPtr() once the item has been consumed in place. + * @warning The fifo must not be empty when calling this. + * + * The item is written to the free list before it is popped from the + * filled list, so it can never be dropped. The write cannot fail: both + * queues hold size + 1 slots while only size items circulate. + */ + void recycle_front() { + ItemType *it = filled_slots.frontPtr(); + assert(it != nullptr); + [[maybe_unused]] const bool ok = free_slots.write(std::move(*it)); + assert(ok); + filled_slots.popFront(); + } template void push_value(Args &&...recordArgs) { while (!filled_slots.write(std::forward(recordArgs)...)) diff --git a/include/aare/ClusterCollector.hpp b/include/aare/ClusterCollector.hpp index 23666b7d..aa4e334c 100644 --- a/include/aare/ClusterCollector.hpp +++ b/include/aare/ClusterCollector.hpp @@ -2,7 +2,11 @@ #pragma once #include #include +#include +#include +#include +#include "aare/Backoff.hpp" #include "aare/ClusterFinderMT.hpp" #include "aare/ClusterVector.hpp" #include "aare/ProducerConsumerQueue.hpp" @@ -13,48 +17,64 @@ namespace aare { template >> class ClusterCollector { - ProducerConsumerQueue> *m_source; + using SourceQueue = ProducerConsumerQueue>; + + SourceQueue *m_source; std::atomic m_stop_requested{false}; std::atomic m_stopped{true}; - std::chrono::milliseconds m_default_wait{1}; std::thread m_thread; std::vector> m_clusters; void process() { - m_stopped = false; + m_stopped.store(false, std::memory_order_release); fmt::print("ClusterCollector started\n"); - while (!m_stop_requested || !m_source->isEmpty()) { + Backoff backoff; + while (true) { if (ClusterVector *clusters = m_source->frontPtr(); clusters != nullptr) { + backoff.reset(); m_clusters.push_back(std::move(*clusters)); m_source->popFront(); - } else { - std::this_thread::sleep_for(m_default_wait); + continue; } + + if (m_stop_requested.load(std::memory_order_acquire)) + break; + + backoff.pause(); } fmt::print("ClusterCollector stopped\n"); - m_stopped = true; + m_stopped.store(true, std::memory_order_release); } public: - ClusterCollector(ClusterFinderMT *source) { - m_source = source->sink(); + explicit ClusterCollector(SourceQueue *source) : m_source(source) { m_thread = std::thread(&ClusterCollector::process, this); // only one process does that so why isnt it // automatically written to m_cluster in collect // - instead of writing first to m_sink? } + + template ().sink()), + SourceQueue *>, + int> = 0> + explicit ClusterCollector(Finder *source) + : ClusterCollector(source->sink()) {} + void stop() { - m_stop_requested = true; - m_thread.join(); + m_stop_requested.store(true, std::memory_order_release); + if (m_thread.joinable()) + m_thread.join(); } std::vector> steal_clusters() { - if (!m_stopped) { + if (!m_stopped.load(std::memory_order_acquire)) { throw std::runtime_error("ClusterCollector is still running"); } return std::move(m_clusters); } }; -} // namespace aare \ No newline at end of file +} // namespace aare diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index 83b5c780..33bd4636 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -3,6 +3,7 @@ #include "aare/ClusterFile.hpp" #include "aare/ClusterVector.hpp" #include "aare/Dtype.hpp" +#include "aare/FastPedestal.hpp" #include "aare/NDArray.hpp" #include "aare/NDView.hpp" #include "aare/Pedestal.hpp" @@ -18,6 +19,13 @@ struct no_2x2_cluster { ClusterType::cluster_size_x > 2 && ClusterType::cluster_size_y > 2; }; +/** + * @brief Find fixed-size photon clusters using a per-pixel pedestal and noise + * threshold. + * @tparam ClusterType Output cluster type; both dimensions must exceed 2. + * @tparam FRAME_TYPE Input pixel type. + * @tparam PEDESTAL_TYPE Type used for pedestal and threshold calculations. + */ template , typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double, typename = std::enable_if_t::value>> @@ -26,51 +34,86 @@ 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 - * @param image_size size of the image - * @param cluster_size size of the cluster (x, y) - * @param nSigma number of sigma above the pedestal to consider a photon - * @param capacity initial capacity of the cluster vector - * + * @brief Construct a cluster finder with an empty pedestal. + * @param image_size Image shape as (rows, columns). + * @param nSigma Per-pixel noise threshold multiplier. + * @param capacity Initial cluster-vector capacity. */ ClusterFinder(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0, size_t capacity = 1000000) : 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_threshold({image_size[0], image_size[1]}, 0), + m_pd_corrected_frame({image_size[0], image_size[1]}, 0) { LOG(logDEBUG) << "ClusterFinder: " << "image_size: " << image_size[0] << "x" << image_size[1] << ", nSigma: " << nSigma << ", capacity: " << capacity; } - void set_nSigma(PEDESTAL_TYPE nSigma) { m_nSigma = nSigma; } - - PEDESTAL_TYPE get_nSigma() const { return m_nSigma; } - - void push_pedestal_frame(NDView frame) { - m_pedestal.push(frame); + /** + * @brief Set the noise multiplier used for threshold calculation and + * recompute the threshold. + * @param nSigma New per-pixel noise multiplier. + */ + void set_nSigma(PEDESTAL_TYPE nSigma) { + m_nSigma = nSigma; + update_threshold(); } - NDArray pedestal() { return m_pedestal.mean(); } - NDArray noise() { return m_pedestal.std(); } - void clear_pedestal() { m_pedestal.clear(); } + /** @brief Return the current noise multiplier used for threshold + * calculation. */ + PEDESTAL_TYPE get_nSigma() const { return m_nSigma; } /** - * @brief Move the clusters from the ClusterVector in the ClusterFinder to a - * new ClusterVector and return it. - * @param realloc_same_capacity if true the new ClusterVector will have the - * same capacity as the old one + * @brief Add a dark frame to the pedestal estimator. * + * The threshold is initialized automatically when the pedestal first + * becomes ready. Later frames update the ready pedestal. + * @param frame Dark frame matching the configured image shape. + * @throws std::runtime_error if the frame shape does not match. + */ + void push_pedestal_frame(NDView frame) { + if (!m_pedestal.ready()) { + m_pedestal.add_init_frame(frame); + // Initialize the threshold when the pedestal becomes ready. + if (m_pedestal.ready()) { + update_threshold(); + } + } else { + m_pedestal.push_ema(frame); + } + } + + /** @brief Return a copy of the per-pixel pedestal mean. */ + NDArray pedestal() { return m_pedestal.mean(); } + + /** @brief Return the per-pixel pedestal standard deviation (noise). */ + NDArray noise() { return m_pedestal.std(); } + + /** @brief Clear the pedestal and mark it as not ready. */ + void clear_pedestal() { m_pedestal.clear(); } + + /** @brief Recompute the threshold as noise multiplied by nSigma. */ + void update_threshold() { m_threshold = m_pedestal.std() * m_nSigma; } + + /** + * @brief Move out all accumulated clusters and reset the internal vector. + * @param realloc_same_capacity Preserve the previous capacity when true. + * @return The accumulated clusters and their frame metadata. */ ClusterVector steal_clusters(bool realloc_same_capacity = false) { @@ -81,109 +124,179 @@ class ClusterFinder { m_clusters = ClusterVector{}; return tmp; } + + private: + /** + * @brief Process a single pixel: scan its cluster window, decide whether it + * is a photon or a pedestal value, and store the cluster if needed. + * @tparam CheckBounds Skip out-of-image neighbours when true; assume the + * complete window is in bounds when false. Skipped cluster values remain 0. + */ + template + void process_pixel(const NDView &frame, const int iy, + const int ix) { + + constexpr int dy = ClusterSizeY / 2; + constexpr int dx = ClusterSizeX / 2; + constexpr int has_center_pixel_x = ClusterSizeX % 2; + constexpr int has_center_pixel_y = ClusterSizeY % 2; + + PEDESTAL_TYPE max = std::numeric_limits::lowest(); + PEDESTAL_TYPE total = 0; + + const int cols = static_cast(frame.shape(1)); + const int rows = static_cast(frame.shape(0)); + const auto center = + (static_cast(iy) * static_cast(cols)) + + static_cast(ix); + const auto *corrected = m_pd_corrected_frame.data(); + const PEDESTAL_TYPE threshold = m_threshold.data()[center]; + const PEDESTAL_TYPE value = corrected[center]; + + if (value < -threshold) + return; // NEGATIVE_PEDESTAL, nothing to do for this pixel + // TODO! No pedestal update??? + + for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { + const int y = iy + ir; + if constexpr (CheckBounds) { + if (y < 0 || y >= rows) + continue; + } + + const auto *row = corrected + static_cast(y) * cols; + for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) { + const int x = ix + ic; + if constexpr (CheckBounds) { + if (x < 0 || x >= cols) + continue; + } + + const PEDESTAL_TYPE val = row[x]; + total += val; + max = std::max(max, val); + } + } + + if ((max > threshold)) { + if (value < max) + return; // Not max go to the next pixel, no pedestal update + } else if (total > c3 * threshold) { + // pass, store the cluster below + } else { + m_pedestal.push_ema_unchecked(center, frame.data()[center]); + return; // It was a pedestal value nothing to store + } + + // Store cluster + if (value == max) { + ClusterType cluster{}; + cluster.x = ix; + cluster.y = iy; + + int i = 0; + for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { + const int y = iy + ir; + for (int ic = -dx; ic < dx + has_center_pixel_x; ic++, i++) { + const int x = ix + ic; + if constexpr (CheckBounds) { + if (x < 0 || x >= cols || y < 0 || y >= rows) + continue; + } + + const PEDESTAL_TYPE corrected_value = + corrected[(static_cast(y) * cols) + x]; + // If the cluster type is an integral type, and the + // pedestal is a floating point type then we need to + // round the value before storing it + if constexpr (std::is_integral_v && + std::is_floating_point_v) { + cluster.data[i] = + static_cast(std::lround(corrected_value)); + } + // On the other hand if both are floating point or both + // are integral then we can just static cast directly + else { + cluster.data[i] = static_cast(corrected_value); + } + } + } + + // Add the cluster to the output ClusterVector + m_clusters.push_back(cluster); + } + } + + public: + /** + * @brief Find clusters in one frame and update eligible pedestal pixels. + * @param frame Input frame matching the configured image shape. + * @param frame_number Metadata assigned to the accumulated clusters. + * @pre frame has the same shape as image_size passed to the constructor. + * @throws std::runtime_error if the pedestal is not ready. + * @note Clusters accumulate until steal_clusters() is called. + */ void find_clusters(NDView frame, uint64_t frame_number = 0) { + // // 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 = - 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; + if (!m_pedestal.ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot find clusters"); + } + + constexpr int dy = ClusterSizeY / 2; + constexpr int dx = ClusterSizeX / 2; + constexpr int has_center_pixel_x = ClusterSizeX % 2; + constexpr int has_center_pixel_y = ClusterSizeY % 2; + + // Largest neighbour offset below/right of the current pixel. Pixels + // further than this from an edge have their whole window in bounds. + constexpr int down = dy + has_center_pixel_y - 1; + constexpr int right = dx + has_center_pixel_x - 1; m_clusters.set_frame_number(frame_number); - for (int iy = 0; iy < frame.shape(0); iy++) { - for (int ix = 0; ix < frame.shape(1); ix++) { - PEDESTAL_TYPE max = std::numeric_limits::min(); - PEDESTAL_TYPE total = 0; + const int rows = static_cast(frame.shape(0)); + const int cols = static_cast(frame.shape(1)); - // 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)); + // TODO! See if we can get the same performace using the operator- + // m_pd_corrected_frame = frame - m_pedestal.view(); - if (value < -m_nSigma * rms) - continue; // NEGATIVE_PEDESTAL go to next pixel - // TODO! No pedestal update??? + // here we should be able to safely assume that the frame and corrected + // frame have the same size + auto n_pixels = frame.size(); + auto pd = m_pedestal.view().data(); + auto corrected = m_pd_corrected_frame.data(); + auto frame_data = frame.data(); + for (ssize_t i = 0; i < n_pixels; i++) { + corrected[i] = static_cast(frame_data[i]) - pd[i]; + } - for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { - for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) { - 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); + // Interior pixels can skip the per-neighbour bounds checks; pixels + // within dx/dy of an edge take the bounds-checked path. Iteration order + // (row-major, increasing ix) is preserved so results are identical. + const int ix_begin = dx; + const int ix_end = cols - right; // exclusive - total += val; - max = std::max(max, val); - } - } - } + for (int iy = 0; iy < rows; iy++) { + const bool interior_row = iy >= dy && iy < rows - down; - if ((max > m_nSigma * rms)) { - if (value < max) - continue; // Not max go to the next pixel - // but also no pedestal update - } else if (total > c3 * m_nSigma * rms) { - // pass - } else { - // m_pedestal.push(iy, ix, frame(iy, ix)); // Safe option - m_pedestal.push_fast( - iy, ix, - frame(iy, - ix)); // Assume we have reached n_samples in the - // pedestal, slight performance improvement - continue; // It was a pedestal value nothing to store - } - - // Store cluster - if (value == max) { - ClusterType cluster{}; - cluster.x = ix; - cluster.y = iy; - - // Fill the cluster data since we have a photon to store - // It's worth redoing the look since most of the time we - // don't have a photon - int i = 0; - for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { - for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) { - if (ix + ic >= 0 && ix + ic < frame.shape(1) && - iy + ir >= 0 && iy + ir < frame.shape(0)) { - - // If the cluster type is an integral type, and - // the pedestal is a floating point type then we - // need to round the value before storing it - if constexpr (std::is_integral_v && - std::is_floating_point_v< - PEDESTAL_TYPE>) { - auto tmp = std::lround( - frame(iy + ir, ix + ic) - - m_pedestal.mean(iy + ir, ix + ic)); - cluster.data[i] = static_cast(tmp); - } - // On the other hand if both are floating point - // or both are integral then we can just static - // cast directly - else { - auto tmp = - frame(iy + ir, ix + ic) - - m_pedestal.mean(iy + ir, ix + ic); - cluster.data[i] = static_cast(tmp); - } - } - i++; - } - } - - // Add the cluster to the output ClusterVector - m_clusters.push_back(cluster); - } + if (!interior_row || ix_begin >= ix_end) { + for (int ix = 0; ix < cols; ix++) + process_pixel(frame, iy, ix); + continue; } + + for (int ix = 0; ix < ix_begin; ix++) + process_pixel(frame, iy, ix); + for (int ix = ix_begin; ix < ix_end; ix++) + process_pixel(frame, iy, ix); + for (int ix = ix_end; ix < cols; ix++) + process_pixel(frame, iy, ix); } } }; -} // namespace aare \ No newline at end of file +} // namespace aare diff --git a/include/aare/ClusterFinderMT.hpp b/include/aare/ClusterFinderMT.hpp index d4ff0018..a4995d0d 100644 --- a/include/aare/ClusterFinderMT.hpp +++ b/include/aare/ClusterFinderMT.hpp @@ -6,11 +6,13 @@ #include #include +#include "aare/Backoff.hpp" +#include "aare/CircularFifo.hpp" #include "aare/ClusterFinder.hpp" #include "aare/NDArray.hpp" #include "aare/ProducerConsumerQueue.hpp" #include "aare/logger.hpp" - +#include namespace aare { enum class FrameType { @@ -18,10 +20,40 @@ enum class FrameType { PEDESTAL, }; -struct FrameWrapper { - FrameType type; - uint64_t frame_number; - NDArray data; +/** + * @brief Ticket identifying a frame buffer in a FramePool. Trivially + * copyable, the buffer itself never travels through the queues. + */ +struct FrameRef { + FrameType type{}; + uint32_t slot{}; + uint64_t frame_number{}; +}; + +/** + * @brief Fixed set of frame buffers allocated once at construction. + * + * Buffers are addressed by slot index and are never moved or reallocated, so + * the queues only need to pass a FrameRef around. This keeps the per frame + * allocation out of the hot path entirely. + */ +class FramePool { + std::vector> m_buffers; + + public: + FramePool(size_t depth, Shape<2> shape) { + m_buffers.reserve(depth); // no reallocation, slots stay stable + for (size_t i = 0; i < depth; ++i) { + m_buffers.emplace_back(shape); + } + } + + NDArray &operator[](uint32_t slot) { return m_buffers[slot]; } + const NDArray &operator[](uint32_t slot) const { + return m_buffers[slot]; + } + + size_t size() const { return m_buffers.size(); } }; /** @@ -42,17 +74,18 @@ class ClusterFinderMT { size_t m_current_thread{0}; size_t m_n_threads{0}; using Finder = ClusterFinder; - using InputQueue = ProducerConsumerQueue; + using InputQueue = CircularFifo; using OutputQueue = ProducerConsumerQueue>; std::vector> m_input_queues; std::vector> m_output_queues; + std::vector> m_frame_pools; OutputQueue m_sink{1000}; // All clusters go into this queue std::vector> m_cluster_finders; std::vector m_threads; std::thread m_collect_thread; - std::chrono::milliseconds m_default_wait{1}; + std::chrono::microseconds m_default_wait{50}; private: std::atomic m_stop_requested{false}; @@ -65,31 +98,35 @@ class ClusterFinderMT { void process(int thread_id) { auto cf = m_cluster_finders[thread_id].get(); auto q = m_input_queues[thread_id].get(); - bool realloc_same_capacity = true; + auto *pool = m_frame_pools[thread_id].get(); + Backoff backoff; while (!m_stop_requested || !q->isEmpty()) { - if (FrameWrapper *frame = q->frontPtr(); frame != nullptr) { + if (FrameRef *ref = q->frontPtr(); ref != nullptr) { + backoff.reset(); + auto view = (*pool)[ref->slot].view(); - switch (frame->type) { - case FrameType::DATA: - cf->find_clusters(frame->data.view(), frame->frame_number); + switch (ref->type) { + case FrameType::DATA: { + cf->find_clusters(view, ref->frame_number); + // Steal before the write so a failed write cannot drop the + // clusters by re-stealing an empty vector on retry. + auto clusters = cf->steal_clusters(true); while (!m_output_queues[thread_id]->write( - cf->steal_clusters(realloc_same_capacity))) { - std::this_thread::sleep_for(m_default_wait); + std::move(clusters))) { + backoff.pause(); } - break; - + } case FrameType::PEDESTAL: - m_cluster_finders[thread_id]->push_pedestal_frame( - frame->data.view()); + m_cluster_finders[thread_id]->push_pedestal_frame(view); break; } - // frame is processed now discard it - m_input_queues[thread_id]->popFront(); + // frame is processed, hand the buffer back to the free list + q->recycle_front(); } else { - std::this_thread::sleep_for(m_default_wait); + backoff.pause(); } } } @@ -100,18 +137,24 @@ class ClusterFinderMT { */ void collect() { bool empty = true; + Backoff backoff; while (!m_stop_requested || !empty || !m_processing_threads_stopped) { - empty = true; + bool moved_any = false; for (auto &queue : m_output_queues) { - if (!queue->isEmpty()) { - - while (!m_sink.write(std::move(*queue->frontPtr()))) { - std::this_thread::sleep_for(m_default_wait); + while (auto *front = queue->frontPtr()) { + while (!m_sink.write(std::move(*front))) { + backoff.pause(); } queue->popFront(); - empty = false; + moved_any = true; } } + empty = !moved_any; + if (moved_any) { + backoff.reset(); + } else { + backoff.pause(); + } } } @@ -124,16 +167,22 @@ class ClusterFinderMT { * @param capacity initial capacity of the cluster vector. Should match * expected number of clusters in a frame per frame. * @param n_threads number of threads to use + * @param queue_depth number of frame buffers per thread. These are + * allocated once and recycled, so the total resident frame memory is + * n_threads * queue_depth * frame size. Keeping the in flight data below + * the L3 size keeps the per frame copy cheap. */ ClusterFinderMT(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0, - size_t capacity = 2000, size_t n_threads = 3) + size_t capacity = 2000, size_t n_threads = 3, + size_t queue_depth = 16) : m_n_threads(n_threads) { LOG(logDEBUG1) << "ClusterFinderMT: " << "image_size: " << image_size[0] << "x" << image_size[1] << ", nSigma: " << nSigma << ", capacity: " << capacity - << ", n_threads: " << n_threads; + << ", n_threads: " << n_threads + << ", queue_depth: " << queue_depth; for (size_t i = 0; i < n_threads; i++) { m_cluster_finders.push_back( @@ -142,7 +191,13 @@ class ClusterFinderMT { image_size, nSigma, capacity)); } for (size_t i = 0; i < n_threads; i++) { - m_input_queues.emplace_back(std::make_unique(200)); + m_frame_pools.emplace_back( + std::make_unique(queue_depth, image_size)); + m_input_queues.emplace_back(std::make_unique( + static_cast(queue_depth), [](size_t slot) { + return FrameRef{FrameType::DATA, + static_cast(slot), 0}; + })); m_output_queues.emplace_back(std::make_unique(200)); } // TODO! Should we start automatically? @@ -212,13 +267,22 @@ class ClusterFinderMT { * expected to be dark. No photon finding is done. Just pedestal update. */ void push_pedestal_frame(NDView frame) { - FrameWrapper fw{FrameType::PEDESTAL, 0, - NDArray(frame)}; // TODO! copies the data! + for (size_t i = 0; i < m_n_threads; ++i) { + auto *q = m_input_queues[i].get(); + Backoff backoff; - for (auto &queue : m_input_queues) { - while (!queue->write(fw)) { - std::this_thread::sleep_for(m_default_wait); + FrameRef ref; + while (!q->try_pop_free(ref)) { + backoff.pause(); } + + ref.type = FrameType::PEDESTAL; + ref.frame_number = 0; + (*m_frame_pools[i])[ref.slot].copy_from(frame); + + // Cannot fail, the free list is what limits how many frames are + // in flight so there is always room in the filled list. + q->try_push_value(ref); } } @@ -228,11 +292,29 @@ class ClusterFinderMT { * @note Spin locks with a default wait if the queue is full. */ void find_clusters(NDView frame, uint64_t frame_number = 0) { - FrameWrapper fw{FrameType::DATA, frame_number, - NDArray(frame)}; // TODO! copies the data! - while (!m_input_queues[m_current_thread % m_n_threads]->write(fw)) { - std::this_thread::sleep_for(m_default_wait); + const size_t tid = m_current_thread % m_n_threads; + auto *q = m_input_queues[tid].get(); + Backoff backoff; + + FrameRef ref; + while (!q->try_pop_free(ref)) { + backoff.pause(); } + + ref.type = FrameType::DATA; + ref.frame_number = frame_number; + + // DualTimer dt; + (*m_frame_pools[tid])[ref.slot].copy_from(frame); + // auto [wall_ns, cpu_ns] = dt.elapsed(); + // std::cerr << "ClusterFinderMT: find_clusters: copied frame " + // << frame_number << " took " << wall_ns/1000.0 << " wall_us + // and " + // << cpu_ns/1000.0 << " cpu_us" << std::endl; + + // Cannot fail, the free list is what limits how many frames are in + // flight so there is always room in the filled list. + q->try_push_value(ref); m_current_thread++; } @@ -245,6 +327,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/ClusterVector.hpp b/include/aare/ClusterVector.hpp index 0adcb5e3..a4141857 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -53,11 +53,8 @@ class ClusterVector> { m_data.reserve(capacity); } - // Move constructor - ClusterVector(ClusterVector &&other) noexcept - : m_data(other.m_data), m_frame_number(other.m_frame_number) { - other.m_data.clear(); - } + ClusterVector(ClusterVector &&other) noexcept = default; + ClusterVector &operator=(ClusterVector &&other) noexcept = default; /** * @brief Create a copy of the clustervector by filtering clusters in the @@ -80,16 +77,16 @@ class ClusterVector> { return result; } - // Move assignment operator - ClusterVector &operator=(ClusterVector &&other) noexcept { - if (this != &other) { - m_data = other.m_data; - m_frame_number = other.m_frame_number; - other.m_data.clear(); - other.m_frame_number = 0; - } - return *this; - } + // // Move assignment operator + // ClusterVector &operator=(ClusterVector &&other) noexcept { + // if (this != &other) { + // m_data = other.m_data; + // m_frame_number = other.m_frame_number; + // other.m_data.clear(); + // other.m_frame_number = 0; + // } + // return *this; + // } /** * @brief Sum the pixels in each cluster diff --git a/include/aare/FastPedestal.hpp b/include/aare/FastPedestal.hpp new file mode 100644 index 00000000..48aec4bf --- /dev/null +++ b/include/aare/FastPedestal.hpp @@ -0,0 +1,424 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once +#include "aare/File.hpp" +#include "aare/Frame.hpp" +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" +#include +#include + +namespace aare { + +/** + * @brief Maintain per-pixel mean, population variance, and standard deviation. + * + * Initialization accumulates exactly n_samples frames. Subsequent push_ema() + * update the exponential moving average initialized with the mean using + * a smoothing factor of (1/ n_samples). Internal moments are stored + * in double precision. + * @tparam PEDESTAL_TYPE Type returned for the mean, variance, and standard + * deviation. + */ +template class FastPedestal { + + // Did we accumulate enough samples and updated the mean? + bool m_ready = false; + + uint32_t m_rows; + uint32_t m_cols; + + uint32_t m_samples; + double m_inv_samples; // precompute 1/m_samples for faster division + uint32_t m_cur_samples = 0; // number of samples accumulated so far + + // For cache locality we want to keep sum and sum2 close. Should + // improve performance for random access. + struct Entry { + double sum; + double sum2; + }; + NDArray m_sum; + + // Cache mean since it is used over and over in the ClusterFinder + // This optimization is related to the access pattern of the ClusterFinder + // Relies on having more reads than pushes to the pedestal + // But also makes sense when subtracting the pedestal from the frame + NDArray m_mean; + + // Helper function to convert row and column indices to a flat index + // used to provide both row column and flat index access to the pedestal + size_t rc_to_index(uint32_t row, uint32_t col) const { + return (static_cast(row) * m_cols) + col; + } + + public: + /** + * @brief Construct an empty pedestal that becomes ready after n_samples + * initialization frames. + * @param rows Number of image rows. + * @param cols Number of image columns. + * @param n_samples Number of initialization frames and reciprocal of the + * weight assigned to each subsequent value. + * @throws std::runtime_error if rows, cols, or n_samples is zero. + */ + 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({rows, cols}, Entry{0, 0}), + m_mean({rows, cols}, PEDESTAL_TYPE(0)) { + if (!(rows > 0 && cols > 0 && n_samples > 0)) { + throw std::runtime_error( + fmt::format("Invalid parameters for FastPedestal: rows={}, " + "cols={}, n_samples={} need to be positive", + rows, cols, n_samples)); + } + } + + ~FastPedestal() = default; + + /** + * @brief Return a non-owning view of the cached mean. + * @throws std::runtime_error if ready() is false. + * @note The caller must treat the data as read-only and must not retain the + * view after this object is destroyed, moved, or assigned. + */ + NDView view() const { + if (!ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot return view"); + } + return m_mean.view(); + } + + /** + * @brief Return a copy of the cached mean. + * @throws std::runtime_error if ready() is false. + */ + NDArray mean() { + if (!ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot return mean"); + } + return m_mean; + } + + /** + * @brief Return the cached mean at (row, col). + * @throws std::runtime_error if ready() is false or either index is out of + * range. + */ + PEDESTAL_TYPE mean(uint32_t row, uint32_t col) const { + if (!ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot return mean"); + } + if (row >= m_rows || col >= m_cols) { + throw std::runtime_error( + fmt::format("Invalid indices for FastPedestal mean: row={}, " + "col={} must be in [0, {}), [0, {})", + row, col, m_rows, m_cols)); + } + return m_mean(row, col); + } + + /** + * @brief Return the cached mean at a flat row-major index. + * @pre ready() is true, the cache is current, and index is valid; the index + * is not checked. + */ + PEDESTAL_TYPE mean_unchecked(ssize_t index) const { return m_mean[index]; } + + /** + * @brief Calculate and return the population variance of every pixel. + * @throws std::runtime_error if ready() is false. + * @note The result is normalized by n_samples. + */ + NDArray variance() { + if (!ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot return variance"); + } + NDArray res({m_rows, m_cols}); + for (ssize_t i = 0; i < m_sum.size(); ++i) { + res[i] = variance_unchecked(i); + } + return res; + } + + /** + * @brief Calculate the population variance at (row, col). + * @throws std::runtime_error if ready() is false or either index is out of + * range. + */ + PEDESTAL_TYPE variance(const uint32_t row, const uint32_t col) const { + if (!ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot return variance"); + } + if (row >= m_rows || col >= m_cols) { + throw std::runtime_error(fmt::format( + "Invalid indices for FastPedestal variance: row={}, " + "col={} must be in [0, {}), [0, {})", + row, col, m_rows, m_cols)); + } + return variance_unchecked(rc_to_index(row, col)); + } + + /** + * @brief Calculate the population variance at a flat row-major index. + * @pre ready() is true and index is valid; the index is not checked. + */ + PEDESTAL_TYPE variance_unchecked(ssize_t index) const { + const auto &entry = m_sum[index]; + const auto m = entry.sum * m_inv_samples; + return std::fma(-m, m, entry.sum2 * m_inv_samples); + } + + /** + * @brief Calculate and return the population standard deviation of every + * pixel. + * @throws std::runtime_error if ready() is false. + */ + NDArray std() { + if (!ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot return std"); + } + NDArray res({m_rows, m_cols}); + for (ssize_t i = 0; i < m_sum.size(); ++i) { + res[i] = std_unchecked(i); + } + return res; + } + + /** + * @brief Calculate the population standard deviation at (row, col). + * @throws std::runtime_error if ready() is false or either index is out of + * range. + */ + PEDESTAL_TYPE std(const uint32_t row, const uint32_t col) const { + if (!ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot return std"); + } + if (row >= m_rows || col >= m_cols) { + throw std::runtime_error( + fmt::format("Invalid indices for FastPedestal std: row={}, " + "col={} must be in [0, {}), [0, {})", + row, col, m_rows, m_cols)); + } + return std::sqrt(variance(row, col)); + } + + /** + * @brief Calculate the population standard deviation at a flat row-major + * index. + * @pre ready() is true and index is valid; the index is not checked. + */ + PEDESTAL_TYPE std_unchecked(ssize_t index) const { + return std::sqrt(variance_unchecked(index)); + } + + /** + * @brief Return whether initialization is complete (cur_samples() equals + * n_samples()). + */ + bool ready() const { return m_ready; } + + /** + * @brief Return the stored number of accumulated initialization frames. + * @note The value is in [0, n_samples] and does not change during + * steady-state pushes. + */ + uint32_t cur_samples() const { return m_cur_samples; } + + /** + * @brief Zero the moments and cached mean, and mark the pedestal not ready. + */ + void clear() { + m_sum = Entry{0., 0.}; + m_mean = PEDESTAL_TYPE(0.); + m_cur_samples = 0; + m_ready = false; + } + + /** + * @brief Update every pixel using the steady-state exponential estimator, + * giving the new value weight 1 / n_samples. + * @param frame Frame whose shape must exactly match the pedestal. + * @throws std::runtime_error if the shape differs or ready() is false. + */ + template void push_ema(NDView frame) { + if (frame.shape() != std::array{m_rows, m_cols}) { + throw std::runtime_error( + "Frame shape does not match pedestal shape"); + } + + if (!ready()) { + throw std::runtime_error("Pedestal is not ready, cannot push"); + } + + const auto size = static_cast(m_rows) * m_cols; + const auto *data = frame.data(); + for (std::size_t index = 0; index < size; ++index) { + push_ema_unchecked(index, data[index]); + } + } + + /** + * @brief Update every pixel from a Frame using the steady-state estimator. + * @tparam T Actual pixel type stored in frame; this is not runtime-checked. + * @param frame Frame whose shape must exactly match the pedestal. + * @throws std::runtime_error if the shape differs or ready() is false. + */ + template void push_ema(Frame &frame) { + push_ema(frame.view()); + } + + /** + * @brief Update the exponential moving average with smoothing factor + * 1/n_samples + * @param row Pixel row. + * @param col Pixel column. + * @param val New pixel value. + * @pre row and col are valid; indices are not checked. + * @throws std::runtime_error if ready() is false. + */ + template + void push_ema(const uint32_t row, const uint32_t col, const T val) { + if (!ready()) { + throw std::runtime_error("Pedestal is not ready, cannot push"); + } + + push_ema_unchecked(rc_to_index(row, col), val); + } + + /** + * @brief Update one pixel and its cached mean without runtime checks in + * release builds. + * @param index Flat row-major pixel index. + * @param value New pixel value, with weight 1 / n_samples. + * @pre ready() is true and index is a valid flat row-major index. These + * preconditions are asserted only in debug builds. + */ + template + void push_ema_unchecked(const std::size_t index, const T value) noexcept { + assert(m_ready); + assert(index < static_cast(m_sum.size())); + + const auto val = static_cast(value); + auto &entry = m_sum[index]; + entry.sum += val - entry.sum * m_inv_samples; + entry.sum2 += val * val - entry.sum2 * m_inv_samples; + m_mean[index] = static_cast(entry.sum * m_inv_samples); + } + + /** + * @brief Accumulate one initialization frame. + * @param frame Frame whose shape must exactly match the pedestal. + * @throws std::runtime_error if the shape differs or n_samples frames have + * already been accumulated. + * @note The statistics can be accessed and ready() becomes true only after + * the n_samples frames have been added. + */ + template void add_init_frame(NDView frame) { + if (frame.shape() != std::array{m_rows, m_cols}) { + throw std::runtime_error( + "Frame shape does not match pedestal shape"); + } + + // if the pedestal is already initialized we cannot add more frames + if (ready()) { + throw std::runtime_error("Pedestal initialization is already done"); + } + + for (ssize_t i = 0; i < m_sum.size(); ++i) { + const auto val = static_cast(frame[i]); + auto &entry = m_sum[i]; + entry.sum += val; + entry.sum2 += val * val; + } + m_cur_samples += 1; + + if (m_cur_samples == m_samples) { + update_mean(); + m_ready = true; + } + } + + /** + * @brief Initialize from n_samples frames after skip_first, then apply all + * remaining file frames as steady-state updates. + * @tparam T Pixel representation stored in the file; this is not checked. + * @param filename Input image file. Its dimensions define the pedestal + * shape. + * @param n_samples Number of initialization frames. + * @param skip_first Number of leading frames to ignore. + * @throws std::runtime_error if fewer than n_samples frames remain after + * skip_first, or if any constructor argument is invalid. + */ + template + static FastPedestal from_file(const std::filesystem::path &filename, + uint32_t n_samples = 1000, + uint32_t skip_first = 0) { + File f(filename); + + const auto total_frames = f.total_frames(); + const auto first_frame = static_cast(skip_first); + const auto initialization_frames = static_cast(n_samples); + if (first_frame > total_frames || + initialization_frames > total_frames - first_frame) { + throw std::runtime_error( + "File has less frames than the number of samples needed to " + "initialize the pedestal"); + } + + if (skip_first > 0) { + f.seek(static_cast(skip_first)); + } + const auto rows = static_cast(f.rows()); + const auto cols = static_cast(f.cols()); + FastPedestal pedestal(rows, cols, n_samples); + NDArray frame({rows, cols}); + + auto frame_index = first_frame; + const auto initialization_end = first_frame + initialization_frames; + while (frame_index < initialization_end) { + f.read_into(frame.buffer()); + pedestal.template add_init_frame(frame.view()); + frame_index++; + } + + // read the rest of the file + while (frame_index < total_frames) { + f.read_into(frame.buffer()); + pedestal.template push_ema(frame.view()); + frame_index++; + } + return pedestal; + } + + /** @brief Return the number of image rows. */ + uint32_t rows() const { return m_rows; } + + /** @brief Return the number of image columns. */ + uint32_t cols() const { return m_cols; } + + /** + * @brief Return the initialization frame count and steady-state + * update-weight denominator. + */ + uint32_t n_samples() const { return m_samples; } + + private: + /** + * @brief Write the cached mean after the final add_init_frame. All other + * (non initialization) pushes update the cached mean immediately. + */ + void update_mean() { + for (ssize_t i = 0; i < m_sum.size(); i++) { + auto &entry = m_sum[i]; + m_mean[i] = static_cast(entry.sum * m_inv_samples); + } + } +}; +} // namespace aare diff --git a/include/aare/NDArray.hpp b/include/aare/NDArray.hpp index 43964028..fb6c3bb5 100644 --- a/include/aare/NDArray.hpp +++ b/include/aare/NDArray.hpp @@ -146,6 +146,23 @@ class NDArray : public ArrayExpr, Ndim> { */ ~NDArray() { delete[] data_; } + /** + * @brief Copy data from a view of matching shape into this array without + * reallocating. + * @param v view to copy from, must have the same shape as this array + * @throws std::runtime_error if the shapes differ + * + * Use this instead of assigning a new NDArray when the buffer needs to be + * kept, for example when the array is part of a preallocated pool. + */ + void copy_from(const NDView v) { + if (v.shape() != shape_) { + throw std::runtime_error(LOCATION + + "Shape mismatch in NDArray::copy_from"); + } + std::copy(v.begin(), v.end(), begin()); + } + /////////////////////////////////////////////////////////////////////////////// // Iterators and indexing // diff --git a/include/aare/Pedestal.hpp b/include/aare/Pedestal.hpp index 48bb6568..4e68384f 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/PedestalTrackingPixelHistogram.hpp b/include/aare/hist/PedestalTrackingPixelHistogram.hpp index 263ea450..88ca7f2f 100644 --- a/include/aare/hist/PedestalTrackingPixelHistogram.hpp +++ b/include/aare/hist/PedestalTrackingPixelHistogram.hpp @@ -1,4 +1,5 @@ #pragma once +#include "aare/FastPedestal.hpp" #include "aare/NDArray.hpp" #include "aare/NDView.hpp" #include "aare/Pedestal.hpp" @@ -44,7 +45,7 @@ class PedestalTrackingPixelHistogram { // 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_pedestals_; std::vector> partial_std_; // cached for pedestal // tracking @@ -65,6 +66,11 @@ class PedestalTrackingPixelHistogram { // Always the outermost lock; work_mutex_ is taken briefly inside it. mutable std::mutex fill_mutex_; + // Serialises producers of async_queue_ (which is SPSC) and gives + // fill_from_file exclusive ownership of the ingestion path while its + // direct, double-buffered dispatch is active. + std::mutex ingestion_mutex_; + // Async producer/consumer pipeline. SPSC queue feeds the coordinator // thread, which batches queued images before dispatching them. std::unique_ptr async_queue_; @@ -107,8 +113,23 @@ class PedestalTrackingPixelHistogram { void fill_async(NDArray &&image); + /** + * @brief Fill from ordered, parallel-read batches. + * + * After the first batch, reading batch N+1 overlaps histogram processing + * for batch N. Two fixed-capacity read buffers are allocated once and + * reused for the duration of the call. + * + * @param fname input file accepted by File + * @param max_frames maximum frames to process, or all frames for -1 + * @param verbose print periodic progress information + * @param reader_threads number of MultiThreadedFileReader workers + * @param reader_chunk_size frames claimed per reader worker and batch + */ void fill_from_file(const std::filesystem::path &fname, - ssize_t max_frames = -1, bool verbose = false); + ssize_t max_frames = -1, bool verbose = false, + std::size_t reader_threads = 2, + std::size_t reader_chunk_size = 4); void process_pedestal_file(const std::filesystem::path &fname, ssize_t max_frames = -1, bool verbose = false); diff --git a/include/aare/hist/PixelHistogramImpl.hpp b/include/aare/hist/PixelHistogramImpl.hpp index 30b468b9..d68795c8 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 @@ -37,6 +38,9 @@ template class PixelHistogramImpl { void fill(const NDView &frame); void fill(int row, int col, T value); void fill_unchecked(int row, int col, T value); + // Fill using a row-major pixel index, avoiding repeated 3D stride + // calculation in tiled worker loops. The index is intentionally unchecked. + void fill_flat_unchecked(std::size_t pixel, T value); NDArray values() const; // Zero-copy view of the underlying [rows x cols x n_bins] storage. @@ -98,22 +102,29 @@ void PixelHistogramImpl::fill(int row, int col, T value) { template void PixelHistogramImpl::fill_unchecked(int row, int col, T value) { + const auto pixel = + static_cast(row) * static_cast(m_cols) + + static_cast(col); + fill_flat_unchecked(pixel, value); +} + +template +void PixelHistogramImpl::fill_flat_unchecked(std::size_t pixel, + 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; - } + bin = std::clamp(bin, 0, m_n_bins - 1); + auto &cell = m_values.data()[pixel * static_cast(m_n_bins) + + static_cast(bin)]; if constexpr (std::is_integral_v) { - if (m_values(row, col, bin) >= - std::numeric_limits::max()) { + if (cell >= std::numeric_limits::max()) { return; } } - ++m_values(row, col, bin); + ++cell; } template diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 401bf26b..74f35ce3 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -39,7 +39,9 @@ set(PYTHON_FILES aare/ClusterFinder.py aare/ClusterVector.py aare/Cluster.py + aare/FastPedestal.py aare/calibration.py + aare/factory.py aare/experimental.py aare/RawFile.py aare/transform.py diff --git a/python/aare/Cluster.py b/python/aare/Cluster.py index f2ea07eb..7946567d 100644 --- a/python/aare/Cluster.py +++ b/python/aare/Cluster.py @@ -1,6 +1,6 @@ from . import _aare import numpy as np -from .ClusterFinder import _type_to_char +from .factory import _type_to_char def Cluster(x : int, y : int, data, cluster_size=(3,3), dtype = np.int32): @@ -21,4 +21,4 @@ def Cluster(x : int, y : int, data, cluster_size=(3,3), dtype = np.int32): except AttributeError: raise ValueError(f"Unsupported combination of type and cluster size: {dtype}/{cluster_size} when requesting {class_name}") - return cls(x, y, data) \ No newline at end of file + return cls(x, y, data) diff --git a/python/aare/ClusterFinder.py b/python/aare/ClusterFinder.py index eae7ad7d..93af877a 100644 --- a/python/aare/ClusterFinder.py +++ b/python/aare/ClusterFinder.py @@ -1,21 +1,10 @@ # SPDX-License-Identifier: MPL-2.0 from . import _aare import numpy as np +from .factory import _type_to_char _supported_cluster_sizes = [(2,2), (3,3), (5,5), (7,7), (9,9),] -def _type_to_char(dtype): - if dtype == np.int32: - return 'i' - elif dtype == np.float32: - return 'f' - elif dtype == np.float64: - return 'd' - elif dtype == np.int16: - return 'i16' - else: - raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32, np.float32, and np.float64 are supported.") - def _get_class(name, cluster_size, dtype): """ Helper function to get the class based on the name, cluster size, and dtype. diff --git a/python/aare/FastPedestal.py b/python/aare/FastPedestal.py new file mode 100644 index 00000000..38bbeb5b --- /dev/null +++ b/python/aare/FastPedestal.py @@ -0,0 +1,56 @@ +# SPDX-License-Identifier: MPL-2.0 + +import numpy as np + +from .factory import _get_typed_class + + +def _get_fast_pedestal_class(dtype): + return _get_typed_class("FastPedestal", dtype) + + +def FastPedestal(rows, cols, n_samples=1000, dtype=np.float64): + """Create an empty per-pixel running pedestal. + + This factory hides the dtype suffix used by the templated C++ bindings. + Call ``add_init_frame()`` exactly ``n_samples`` times before using the + statistics or calling ``push_ema()``. Subsequent frames have weight + ``1 / n_samples`` in the running mean and population variance. + + Args: + rows: Number of image rows. + cols: Number of image columns. + n_samples: Initialization frame count and steady-state update-weight + denominator. + dtype: Output dtype for the mean, variance, and standard deviation. + Supported values are ``np.float64``, ``np.float32``, and + ``np.int16``. + """ + cls = _get_fast_pedestal_class(dtype) + return cls(rows, cols, n_samples) + + +def from_file(filename, n_samples=1000, skip_first=0, dtype=np.float64): + """Create a FastPedestal from frames in a file. + + After ignoring ``skip_first`` frames, the next ``n_samples`` frames + initialize the pedestal. Every remaining frame is then applied as a + steady-state update. Input frames are read as uint16 data. + + Args: + filename: Input image file. + n_samples: Number of frames used for initialization. + skip_first: Number of leading frames to ignore. + dtype: Output dtype for the mean, variance, and standard deviation. + + Raises: + RuntimeError: If fewer than ``n_samples`` frames remain after + ``skip_first`` or ``n_samples`` is zero. + """ + cls = _get_fast_pedestal_class(dtype) + return cls.from_file( + filename, n_samples=n_samples, skip_first=skip_first + ) + + +FastPedestal.from_file = from_file diff --git a/python/aare/__init__.py b/python/aare/__init__.py index 84224613..72c2573b 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -5,6 +5,16 @@ from . import _aare from . import transform from . import experimental +from ._aare import ( + FastPedestal_d, + FastPedestal_f, + FastPedestal_i16, + Pedestal_d, + Pedestal_f, + Pedestal_i16, + ClusterFinder_Cluster3x3i, + VarClusterFinder, +) from ._aare import ( File, JungfrauDataFile, @@ -20,6 +30,7 @@ from ._aare import corner # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i from ._version import __version__ +from .FastPedestal import FastPedestal from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile from .ClusterVector import ClusterVector from .Cluster import Cluster diff --git a/python/aare/factory.py b/python/aare/factory.py new file mode 100644 index 00000000..4da3fdbe --- /dev/null +++ b/python/aare/factory.py @@ -0,0 +1,36 @@ +# SPDX-License-Identifier: MPL-2.0 + +import numpy as np + +from . import _aare + + +_TYPE_TO_CHAR = { + np.dtype(np.int32): "i", + np.dtype(np.float32): "f", + np.dtype(np.float64): "d", + np.dtype(np.int16): "i16", +} + + +def _type_to_char(dtype): + """Return the suffix used by bindings instantiated for ``dtype``.""" + try: + return _TYPE_TO_CHAR[np.dtype(dtype)] + except (KeyError, TypeError): + supported = ", ".join(str(dtype) for dtype in _TYPE_TO_CHAR) + raise ValueError( + f"Unsupported dtype: {dtype}. Supported dtypes are {supported}." + ) from None + + +def _get_typed_class(name, dtype): + """Return a bound class named ``_``.""" + class_name = f"{name}_{_type_to_char(dtype)}" + try: + return getattr(_aare, class_name) + except AttributeError: + raise ValueError( + f"Unsupported dtype for {name}: {dtype} " + f"(binding {class_name} is not available)." + ) from None diff --git a/python/src/bind_Cluster.hpp b/python/src/bind_Cluster.hpp index eca0881b..5416d04f 100644 --- a/python/src/bind_Cluster.hpp +++ b/python/src/bind_Cluster.hpp @@ -10,7 +10,6 @@ #include namespace py = pybind11; -using pd_type = double; using namespace aare; diff --git a/python/src/bind_ClusterCollector.hpp b/python/src/bind_ClusterCollector.hpp index 2e621eda..01dd03c2 100644 --- a/python/src/bind_ClusterCollector.hpp +++ b/python/src/bind_ClusterCollector.hpp @@ -6,6 +6,8 @@ #include "aare/ClusterVector.hpp" #include "aare/NDView.hpp" #include "aare/Pedestal.hpp" + +#include "module_config.hpp" #include "np_helper.hpp" #include @@ -15,7 +17,6 @@ #include namespace py = pybind11; -using pd_type = double; using namespace aare; @@ -30,7 +31,7 @@ void define_ClusterCollector(py::module &m, const std::string &typestr) { using ClusterType = Cluster; py::class_>(m, class_name.c_str()) - .def(py::init *>()) + .def(py::init *>()) .def("stop", &ClusterCollector::stop) .def( "steal_clusters", diff --git a/python/src/bind_ClusterFileSink.hpp b/python/src/bind_ClusterFileSink.hpp index ade574d4..b366886a 100644 --- a/python/src/bind_ClusterFileSink.hpp +++ b/python/src/bind_ClusterFileSink.hpp @@ -15,8 +15,6 @@ #include namespace py = pybind11; -using pd_type = double; - using namespace aare; #pragma GCC diagnostic push @@ -29,6 +27,8 @@ void define_ClusterFileSink(py::module &m, const std::string &typestr) { using ClusterType = Cluster; + // TODO! adapt to set pedestal type (needs templating of ClusterFileSink) + // or maybe access through base class? py::class_>(m, class_name.c_str()) .def(py::init *, const std::filesystem::path &>()) diff --git a/python/src/bind_ClusterFinder.hpp b/python/src/bind_ClusterFinder.hpp index da479dc5..c13d3305 100644 --- a/python/src/bind_ClusterFinder.hpp +++ b/python/src/bind_ClusterFinder.hpp @@ -6,6 +6,8 @@ #include "aare/ClusterVector.hpp" #include "aare/NDView.hpp" #include "aare/Pedestal.hpp" + +#include "module_config.hpp" #include "np_helper.hpp" #include @@ -15,7 +17,6 @@ #include namespace py = pybind11; -using pd_type = double; using namespace aare; @@ -48,6 +49,8 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) { }) .def("clear_pedestal", &ClusterFinder::clear_pedestal) + .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 4bb81012..252cc2f8 100644 --- a/python/src/bind_ClusterFinderMT.hpp +++ b/python/src/bind_ClusterFinderMT.hpp @@ -6,6 +6,8 @@ #include "aare/ClusterVector.hpp" #include "aare/NDView.hpp" #include "aare/Pedestal.hpp" + +#include "module_config.hpp" #include "np_helper.hpp" #include @@ -15,7 +17,6 @@ #include namespace py = pybind11; -using pd_type = double; using namespace aare; @@ -31,9 +32,10 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) { py::class_>( m, class_name.c_str()) - .def(py::init, pd_type, size_t, size_t>(), + .def(py::init, pd_type, size_t, size_t, size_t>(), py::arg("image_size"), py::arg("n_sigma") = 5.0, - py::arg("capacity") = 2048, py::arg("n_threads") = 3) + py::arg("capacity") = 2048, py::arg("n_threads") = 3, + py::arg("queue_depth") = 16) .def("push_pedestal_frame", [](ClusterFinderMT &self, py::array_t frame) { @@ -46,7 +48,6 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) { py::array_t frame, uint64_t frame_number) { auto view = make_view_2d(frame); self.find_clusters(view, frame_number); - return; }, py::arg(), py::arg("frame_number") = 0) .def_property_readonly( @@ -56,6 +57,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/bind_ClusterVector.hpp b/python/src/bind_ClusterVector.hpp index a3d59f57..8c863183 100644 --- a/python/src/bind_ClusterVector.hpp +++ b/python/src/bind_ClusterVector.hpp @@ -15,7 +15,6 @@ #include namespace py = pybind11; -using pd_type = double; using namespace aare; diff --git a/python/src/bind_FastPedestal.hpp b/python/src/bind_FastPedestal.hpp new file mode 100644 index 00000000..c1911e52 --- /dev/null +++ b/python/src/bind_FastPedestal.hpp @@ -0,0 +1,147 @@ +// SPDX-License-Identifier: MPL-2.0 + +#include "aare/FastPedestal.hpp" +#include "np_helper.hpp" + +#include +#include +#include +#include + +namespace py = pybind11; + +template +void define_fast_pedestal_bindings(py::module &m, const std::string &name) { + + py::class_>( + m, name.c_str(), + "Maintain a per-pixel running mean, population variance, and " + "standard deviation.", + py::buffer_protocol()) + .def(py::init(), py::arg("rows"), + py::arg("cols"), py::arg("n_samples"), + "Construct an empty pedestal. It becomes ready after n_samples " + "calls to add_init_frame().") + .def(py::init(), py::arg("rows"), py::arg("cols"), + "Construct an empty pedestal with n_samples=1000.") + + .def( + "mean", + [](FastPedestal &self) { + auto mean = new NDArray{}; + *mean = self.mean(); + return return_image_data(mean); + }, + "Return a copy of the cached mean. The pedestal must be ready.") + .def( + "var", + [](FastPedestal &self) { + auto variance = new NDArray{}; + *variance = self.variance(); + return return_image_data(variance); + }, + "Return the population variance, normalized by n_samples, as a " + "NumPy array. The pedestal must be ready.") + .def( + "std", + [](FastPedestal &self) { + auto standard_deviation = new NDArray{}; + *standard_deviation = self.std(); + return return_image_data(standard_deviation); + }, + "Return the population standard deviation as a NumPy array. The " + "pedestal must be ready.") + .def( + "view", + [](py::object self_py) { + return py::module_::import("numpy").attr("asarray")(self_py); + }, + "Return a non-owning, non-writable NumPy view of the cached mean. " + "The pedestal must be ready.") + + // We need to buffer protocol to allow for numpy operations using the + // pedestal mean + .def_buffer([](FastPedestal &self) { + auto mean = self.view(); + return py::buffer_info( + const_cast(mean.data()), sizeof(SUM_TYPE), + py::format_descriptor::format(), 2, + {static_cast(mean.shape(0)), + static_cast(mean.shape(1))}, + {static_cast(mean.strides()[0] * sizeof(SUM_TYPE)), + static_cast(mean.strides()[1] * + sizeof(SUM_TYPE))}, + true); + }) + // Subtracting a FastPedestal from a NumPy array + .def( + "__array_ufunc__", + [](py::object self, py::object ufunc, const std::string &method, + py::args inputs, py::kwargs kwargs) -> py::object { + if (method != "__call__" || inputs.size() != 2 || + inputs[1].ptr() != self.ptr() || + py::cast(ufunc.attr("__name__")) != + "subtract") { + return py::reinterpret_borrow( + Py_NotImplemented); + } + + auto mean = + py::module_::import("builtins").attr("memoryview")(self); + return ufunc(inputs[0], mean, **kwargs); + }, + "Support subtracting a FastPedestal from a NumPy array.") + .def("clear", py::overload_cast<>(&FastPedestal::clear), + "Reset all statistics and initialization state to zero.") + .def_property_readonly("rows", &FastPedestal::rows, + "Number of image rows.") + .def_property_readonly("cols", &FastPedestal::cols, + "Number of image columns.") + .def_property_readonly("cur_samples", + &FastPedestal::cur_samples, + "Number of initialization frames accumulated. " + "Steady-state pushes do not change it.") + .def_property_readonly( + "ready", &FastPedestal::ready, + "Whether n_samples initialization frames have been accumulated.") + .def_property_readonly( + "n_samples", &FastPedestal::n_samples, + "Initialization frame count and steady-state update-weight " + "denominator.") + .def( + "clone", + [](FastPedestal &pedestal) { + return FastPedestal(pedestal); + }, + "Return an independent copy of the pedestal and its state.") + .def( + "push_ema", + [](FastPedestal &pedestal, + py::array_t &frame) { + pedestal.push_ema(make_view_2d(frame)); + }, + py::arg("frame").noconvert(), + "Update exponential moving average. The pedstal must " + "already be ready for this update.") + .def( + "add_init_frame", + [](FastPedestal &pedestal, + py::array_t &frame) { + pedestal.add_init_frame(make_view_2d(frame)); + }, + py::arg("frame").noconvert(), + "Accumulate one uint16 initialization frame. Call exactly " + "n_samples times to make the pedestal ready.") + .def_static( + "from_file", + [](const std::filesystem::path &filename, uint32_t n_samples, + uint32_t skip_first) { + return FastPedestal::template from_file( + filename, n_samples, skip_first); + }, + py::arg("filename"), py::arg("n_samples") = 1000, + py::arg("skip_first") = 0, + "Create a pedestal from a uint16 file. Skip skip_first frames, use " + "the next n_samples for initialization, then apply every remaining " + "frame as a steady-state update."); +} diff --git a/python/src/bind_PedestalTrackingPixelHistogram.hpp b/python/src/bind_PedestalTrackingPixelHistogram.hpp index 79842732..35b8b24a 100644 --- a/python/src/bind_PedestalTrackingPixelHistogram.hpp +++ b/python/src/bind_PedestalTrackingPixelHistogram.hpp @@ -153,9 +153,13 @@ void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) { Args: file_path: Path to the file to fill from max_frames: Maximum number of frames to fill from the file (default: -1) + reader_threads: Number of parallel file reader workers (default: 2) + reader_chunk_size: Frames claimed by each reader worker per batch (default: 4) )", py::call_guard(), py::arg("fname"), - py::arg("max_frames") = -1, py::arg("verbose") = false) + py::arg("max_frames") = -1, py::arg("verbose") = false, + py::arg("reader_threads") = std::size_t{2}, + py::arg("reader_chunk_size") = std::size_t{4}) .def("process_pedestal_file", &PedestalTrackingPixelHistogram::process_pedestal_file, R"( diff --git a/python/src/module.cpp b/python/src/module.cpp index 52c830b9..d53baf4f 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -1,6 +1,8 @@ // SPDX-License-Identifier: MPL-2.0 // Files with bindings to the different classes +#include "module_config.hpp" + // New style file naming #include "bind_Cluster.hpp" #include "bind_ClusterCollector.hpp" @@ -11,6 +13,7 @@ #include "bind_ClusterVector.hpp" #include "bind_Defs.hpp" #include "bind_Eta.hpp" +#include "bind_FastPedestal.hpp" #include "bind_Interpolator.hpp" #include "bind_MultiThreadedFileReader.hpp" #include "bind_PedestalTrackingPixelHistogram.hpp" @@ -30,6 +33,7 @@ #include "var_cluster.hpp" // Pybind stuff +#include #include #include @@ -75,6 +79,10 @@ PYBIND11_MODULE(_aare, m) { define_pedestal_tracking_pixel_histogram_bindings(m); define_pedestal_bindings(m, "Pedestal_d"); define_pedestal_bindings(m, "Pedestal_f"); + define_pedestal_bindings(m, "Pedestal_i16"); + define_fast_pedestal_bindings(m, "FastPedestal_d"); + define_fast_pedestal_bindings(m, "FastPedestal_f"); + define_fast_pedestal_bindings(m, "FastPedestal_i16"); define_fit_bindings(m); define_interpolation_bindings(m); define_jungfrau_data_file_io_bindings(m); @@ -106,6 +114,7 @@ PYBIND11_MODULE(_aare, m) { DEFINE_BINDINGS_CLUSTERFINDER(int, 3, 3, uint16_t, i); DEFINE_BINDINGS_CLUSTERFINDER(double, 3, 3, uint16_t, d); DEFINE_BINDINGS_CLUSTERFINDER(float, 3, 3, uint16_t, f); + DEFINE_BINDINGS_CLUSTERFINDER(int16_t, 3, 3, uint16_t, i16); DEFINE_BINDINGS_CLUSTERFINDER(int, 5, 5, uint16_t, i); DEFINE_BINDINGS_CLUSTERFINDER(double, 5, 5, uint16_t, d); diff --git a/python/src/module_config.hpp b/python/src/module_config.hpp new file mode 100644 index 00000000..e2f4e8a7 --- /dev/null +++ b/python/src/module_config.hpp @@ -0,0 +1,5 @@ +#pragma once +#include + +// Configure module wide pedestal type for cluster finding +using pd_type = double; \ No newline at end of file diff --git a/python/src/pedestal.hpp b/python/src/pedestal.hpp index 238848b7..229f3adf 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, diff --git a/python/tests/test_Cluster.py b/python/tests/test_Cluster.py index cafb118e..fece7273 100644 --- a/python/tests/test_Cluster.py +++ b/python/tests/test_Cluster.py @@ -103,15 +103,18 @@ def test_max_sum(): def test_cluster_finder(): """Test ClusterFinder""" + shape = [100,100] + cf = _aare.ClusterFinder_Cluster3x3i(shape) - clusterfinder = _aare.ClusterFinder_Cluster3x3i([100,100]) + #Push 1000 frames to the pedestal + for i in range(1000): + frame = np.random.normal(loc = 100, scale = 5, size = shape).astype(np.uint16) + cf.push_pedestal_frame(frame) + cf.update_threshold() + frame = np.zeros(shape=shape, dtype=np.uint16) + cf.find_clusters(frame) - #frame = np.random.rand(100,100) - frame = np.zeros(shape=[100,100]) - - clusterfinder.find_clusters(frame) - - clusters = clusterfinder.steal_clusters(False) #conversion does not work + clusters = cf.steal_clusters(False) #conversion does not work assert clusters.size == 0 diff --git a/python/tests/test_FastPedestal.py b/python/tests/test_FastPedestal.py new file mode 100644 index 00000000..24c115fc --- /dev/null +++ b/python/tests/test_FastPedestal.py @@ -0,0 +1,145 @@ +import numpy as np +import pytest + +from aare import ( + FastPedestal, + FastPedestal_d, + FastPedestal_f, + FastPedestal_i16, +) + + +@pytest.mark.parametrize( + ("dtype", "pedestal_type"), + [ + (np.float64, FastPedestal_d), + (np.float32, FastPedestal_f), + (np.int16, FastPedestal_i16), + ], +) +def test_fast_pedestal_factory(dtype, pedestal_type): + pedestal = FastPedestal(2, 3, n_samples=4, dtype=dtype) + + assert isinstance(pedestal, pedestal_type) + assert pedestal.rows == 2 + assert pedestal.cols == 3 + assert pedestal.n_samples == 4 + + +def test_fast_pedestal_factory_defaults_to_double(): + assert isinstance(FastPedestal(2, 3), FastPedestal_d) + + +def test_fast_pedestal_factory_rejects_unbound_dtype(): + with pytest.raises(ValueError, match="Unsupported dtype for FastPedestal"): + FastPedestal(2, 3, dtype=np.int32) + + +@pytest.mark.parametrize( + ("kwargs", "expected_n_samples"), + [ + ({"rows": 2, "cols": 3}, 1000), + ({"rows": 2, "cols": 3, "n_samples": 4}, 4), + ], +) +def test_fast_pedestal_binding_accepts_constructor_keywords( + kwargs, expected_n_samples +): + pedestal = FastPedestal_d(**kwargs) + + assert pedestal.rows == 2 + assert pedestal.cols == 3 + assert pedestal.n_samples == expected_n_samples + + +@pytest.mark.parametrize( + ("dtype", "pedestal_type", "expected_dtype"), + [ + (np.float64, FastPedestal_d, np.float64), + (np.float32, FastPedestal_f, np.float32), + (np.int16, FastPedestal_i16, np.int16), + ], +) +def test_fast_pedestal_factory_from_file( + tmp_path, dtype, pedestal_type, expected_dtype +): + frames = np.array( + [[[100, 100]], [[2, 4]], [[4, 6]], [[5, 7]]], dtype=np.uint16 + ) + filename = tmp_path / "frames.npy" + np.save(filename, frames) + + pedestal = FastPedestal.from_file( + filename, n_samples=2, skip_first=1, dtype=dtype + ) + + assert isinstance(pedestal, pedestal_type) + assert pedestal.ready + assert pedestal.cur_samples == 2 + assert pedestal.mean().dtype == expected_dtype + np.testing.assert_array_equal(pedestal.mean(), [[4, 6]]) + + +def test_fast_pedestal_factory_from_file_rejects_unbound_dtype(): + with pytest.raises(ValueError, match="Unsupported dtype for FastPedestal"): + FastPedestal.from_file("unused.npy", dtype=np.int32) + + +def test_fast_pedestal_from_file_rejects_skip_beyond_end(tmp_path): + filename = tmp_path / "frames.npy" + np.save(filename, np.zeros((1, 1, 1), dtype=np.uint16)) + + with pytest.raises(RuntimeError, match="less frames"): + FastPedestal.from_file(filename, n_samples=1, skip_first=2) + + +@pytest.mark.parametrize( + ("pedestal_type", "expected_dtype"), + [(FastPedestal_d, np.float64), (FastPedestal_f, np.float32)], +) +def test_fast_pedestal_initialization(pedestal_type, expected_dtype): + pedestal = pedestal_type(2, 3, 2) + first = np.array([[2, 4, 6], [8, 10, 12]], dtype=np.uint16) + second = np.array([[4, 6, 8], [10, 12, 14]], dtype=np.uint16) + + pedestal.add_init_frame(first) + pedestal.add_init_frame(second) + + + expected_mean = np.array( + [[3, 5, 7], [9, 11, 13]], dtype=expected_dtype + ) + np.testing.assert_array_equal(pedestal.mean(), expected_mean) + np.testing.assert_array_equal(pedestal.std(), np.ones((2, 3))) + + +def test_fast_pedestal_steady_state_push_ema(): + pedestal = FastPedestal_d(1, 2, 2) + pedestal.add_init_frame(np.array([[2, 4]], dtype=np.uint16)) + pedestal.add_init_frame(np.array([[4, 6]], dtype=np.uint16)) + + + pedestal.push_ema(np.array([[6, 8]], dtype=np.uint16)) + + np.testing.assert_array_equal(pedestal.mean(), [[4.5, 6.5]]) + + +def test_fast_pedestal_exposes_read_only_buffer_and_subtraction(): + pedestal = FastPedestal_d(1, 2, 1) + pedestal.add_init_frame(np.array([[2, 4]], dtype=np.uint16)) + + + view = np.asarray(pedestal) + result = np.array([[12, 14]], dtype=np.uint16) - pedestal + + np.testing.assert_array_equal(view, [[2, 4]]) + np.testing.assert_array_equal(result, [[10, 10]]) + assert np.shares_memory(view, pedestal.view()) + assert not view.flags.writeable + + +def test_fast_pedestal_rejects_wrong_shape(): + pedestal = FastPedestal_d(2, 3) + + with pytest.raises(RuntimeError, match="shape"): + pedestal.add_init_frame(np.zeros((2, 2), dtype=np.uint16)) diff --git a/src/ClusterFinderMT.test.cpp b/src/ClusterFinderMT.test.cpp index e8904814..fae8b247 100644 --- a/src/ClusterFinderMT.test.cpp +++ b/src/ClusterFinderMT.test.cpp @@ -21,14 +21,21 @@ class ClusterFinderMTWrapper public: ClusterFinderMTWrapper(Shape<2> image_size, PEDESTAL_TYPE nSigma = 5.0, - size_t capacity = 2000, size_t n_threads = 3) + size_t capacity = 2000, size_t n_threads = 3, + size_t queue_depth = 16) : ClusterFinderMT( - image_size, nSigma, capacity, n_threads) {} + image_size, nSigma, capacity, n_threads, queue_depth) {} size_t get_m_input_queues_size() const { return this->m_input_queues.size(); } + size_t get_m_frame_pools_size() const { return this->m_frame_pools.size(); } + + size_t get_frame_pool_depth(size_t thread_index) const { + return this->m_frame_pools[thread_index]->size(); + } + size_t get_m_output_queues_size() const { return this->m_output_queues.size(); } @@ -100,3 +107,73 @@ TEST_CASE("multithreaded cluster finder", "[.with-data]") { auto clustervec = clustercollector.steal_clusters(); // CHECK(clustervec.size() == ) //dont know how many clusters to expect } + +TEST_CASE("frame buffers are recycled when pushing more frames than the pool " + "holds", + "[.files]") { + using ClusterType = Cluster; + + const size_t n_threads = 2; + const size_t queue_depth = 4; + const size_t n_frames = 5 * queue_depth; + const Shape<2> image_size{10, 10}; + + ClusterFinderMTWrapper cf(image_size, 5, 200, n_threads, + queue_depth); + + CHECK(cf.get_m_frame_pools_size() == n_threads); + for (size_t i = 0; i < n_threads; ++i) { + CHECK(cf.get_frame_pool_depth(i) == queue_depth); + } + + NDArray frame(image_size, 0); + + // More frames than the pool holds, so this only completes if the workers + // return the buffers to the free list. + for (size_t i = 0; i < n_frames; ++i) { + cf.find_clusters(frame.view(), i); + } + + cf.stop(); + + CHECK(cf.m_input_queues_are_empty() == true); +} + +TEST_CASE("cluster collector accepts finders with a matching cluster type") { + using ClusterType = Cluster; + using Finder = ClusterFinderMTWrapper; + using OtherFinder = + ClusterFinderMTWrapper, uint16_t, float>; + + static_assert( + std::is_constructible_v, Finder *>); + static_assert( + !std::is_constructible_v, OtherFinder *>); + + Finder cf({10, 10}); + cf.stop(); + + ClusterCollector collector(&cf); + collector.stop(); + + CHECK(collector.steal_clusters().empty()); +} + +TEST_CASE("cluster collector drains queued clusters when stopped") { + using ClusterType = Cluster; + ProducerConsumerQueue> source(4); + + for (uint64_t frame_number = 1; frame_number <= 3; ++frame_number) { + REQUIRE(source.write(ClusterVector(4, frame_number))); + } + + ClusterCollector collector(&source); + collector.stop(); + + auto clusters = collector.steal_clusters(); + REQUIRE(clusters.size() == 3); + CHECK(source.isEmpty()); + for (size_t i = 0; i < clusters.size(); ++i) { + CHECK(clusters[i].frame_number() == static_cast(i + 1)); + } +} diff --git a/src/ClusterVector.test.cpp b/src/ClusterVector.test.cpp index 33b0b98c..2c098e14 100644 --- a/src/ClusterVector.test.cpp +++ b/src/ClusterVector.test.cpp @@ -1,6 +1,7 @@ // SPDX-License-Identifier: MPL-2.0 #include "aare/ClusterVector.hpp" #include +#include #include #include @@ -21,6 +22,29 @@ TEST_CASE("After pushing back one element the ClusterVector is not empty") { REQUIRE(!cv.empty()); } +TEST_CASE("ClusterVector move constructor preserves contents") { + ClusterVector source(4, 42); + source.push_back(C1{1, 2, {3, 4, 5, 6}}); + + ClusterVector moved(std::move(source)); + + REQUIRE(moved.size() == 1); + CHECK(moved.frame_number() == 42); + CHECK(moved[0].data[3] == 6); +} + +TEST_CASE("ClusterVector move assignment preserves contents") { + ClusterVector source(4, 42); + source.push_back(C1{1, 2, {3, 4, 5, 6}}); + ClusterVector moved(1); + + moved = std::move(source); + + REQUIRE(moved.size() == 1); + CHECK(moved.frame_number() == 42); + CHECK(moved[0].data[3] == 6); +} + TEST_CASE("item_size return the size of the cluster stored") { ClusterVector cv(4); CHECK(cv.item_size() == sizeof(C1)); @@ -275,4 +299,4 @@ TEST_CASE("Gain Map Calculation Index Map") { CHECK(index_map_x == clustertestdata.index_map_x); CHECK(index_map_y == clustertestdata.index_map_y); -} \ No newline at end of file +} diff --git a/src/hist/PedestalTrackingPixelHistogram.cpp b/src/hist/PedestalTrackingPixelHistogram.cpp index 0695c2f3..f02004ab 100644 --- a/src/hist/PedestalTrackingPixelHistogram.cpp +++ b/src/hist/PedestalTrackingPixelHistogram.cpp @@ -1,11 +1,16 @@ #include "aare/hist/PedestalTrackingPixelHistogram.hpp" #include "aare/File.hpp" +#include "aare/MultiThreadedFileReader.hpp" #include +#include #include #include #include +#include +#include #include +#include #include #include @@ -202,17 +207,8 @@ void PedestalTrackingPixelHistogram::worker_loop(int 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)); - } - } + auto frame = image->sub_view(first_row, first_row + local_rows); + my_pedestal.add_init_frame(frame); break; } case WorkKind::UpdateMean: { @@ -220,7 +216,6 @@ void PedestalTrackingPixelHistogram::worker_loop(int thread_id) { // 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) { @@ -240,59 +235,63 @@ void PedestalTrackingPixelHistogram::worker_loop(int thread_id) { // 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. + constexpr std::size_t pixel_tile_size = 256; + const auto local_pixels = static_cast(local_rows) * + static_cast(cols_); + const auto global_pixel_begin = + static_cast(first_row) * + static_cast(cols_); + const auto &frames = *images; - 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); + // Compile the shared traversal once for each mode. if constexpr + // removes the tracking gate entirely from the histogram-only hot + // path, while keeping the runtime mode decision outside the + // per-pixel loops. + const auto fill_tiles = [&](auto tracking, + const AxisType *std_data) { + constexpr bool track_pedestal = decltype(tracking)::value; + + // Each worker retains exclusive ownership of its row shard. + // Pixel tiling keeps a bounded part of that shard's + // [pixel x bin] storage hot while input remains contiguous. + for (std::size_t tile_begin = 0; tile_begin < local_pixels; + tile_begin += pixel_tile_size) { + const auto tile_end = + std::min(tile_begin + pixel_tile_size, local_pixels); + for (const auto &frame : frames) { + const auto *input = + frame.data() + global_pixel_begin + tile_begin; + for (std::size_t local_pixel = tile_begin; + local_pixel < tile_end; ++local_pixel, ++input) { + const FrameType raw = *input; 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); + my_pedestal.mean_unchecked( + static_cast(local_pixel)); + my_hist.fill_flat_unchecked(local_pixel, val); + + if constexpr (track_pedestal) { + const AxisType sigma = std_data[local_pixel]; + if (sigma > AxisType{0.0} && + std::abs(val) < n_sigma * sigma) { + my_pedestal.push_ema_unchecked( + local_pixel, raw); + } } } } } - break; + }; + + if (n_sigma <= AxisType{0.0}) { + fill_tiles(std::false_type{}, nullptr); + } else { + // Frames stay chronological for every pixel. Accepted raw + // values update that pixel's sums and cached mean immediately; + // pixels are otherwise independent. + fill_tiles(std::true_type{}, partial_std_[thread_id].data()); } + break; } } @@ -395,6 +394,11 @@ void PedestalTrackingPixelHistogram::fill_async(NDArray &&image) { "constructor shape"); } + // ProducerConsumerQueue is SPSC. Serialising this short producer-side + // operation also prevents fill_async from racing fill_from_file's direct + // batch dispatch. + std::lock_guard ingestion_lock(ingestion_mutex_); + // 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). @@ -467,35 +471,42 @@ PedestalTrackingPixelHistogram::bin_edges() const { } void PedestalTrackingPixelHistogram::fill_from_file( - const std::filesystem::path &fname, ssize_t max_frames, bool verbose) { + const std::filesystem::path &fname, ssize_t max_frames, bool verbose, + std::size_t reader_threads, std::size_t reader_chunk_size) { 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_); - } - }; + if (max_frames < -1) { + throw std::invalid_argument( + "PedestalTrackingPixelHistogram max_frames must be -1 or " + "non-negative"); + } - File f(fname); - // check that row col matches constructor - if (f.rows() != static_cast(rows_) || - f.cols() != static_cast(cols_)) { + // Preserve the old max_frames behaviour: values beyond EOF are clamped + // rather than rejected by MultiThreadedFileReader. + const auto source_total_frames = File(fname).total_frames(); + const auto n_frames = max_frames == -1 + ? source_total_frames + : std::min(static_cast(max_frames), + source_total_frames); + + experimental::MultiThreadedFileReader reader(fname, reader_threads, + reader_chunk_size, n_frames); + if (reader.rows() != static_cast(rows_) || + reader.cols() != static_cast(cols_)) { throw std::invalid_argument("PedestalTrackingPixelHistogram: Frame in " "file {} has shape ({}, {}) does not match " "constructor shape"); } + if (reader.dtype() != Dtype::UINT16 || + reader.bytes_per_frame() != static_cast(rows_) * + static_cast(cols_) * + sizeof(FrameType)) { + throw std::invalid_argument( + "PedestalTrackingPixelHistogram requires uint16 file frames"); + } - 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(); @@ -503,30 +514,89 @@ void PedestalTrackingPixelHistogram::fill_from_file( 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); + fmt::print("\rProgress: {}/{} ({:.1f}%) {:.1f} FPS ", done, + n_frames, + n_frames == 0 ? 100.0 + : 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)); + using ReadBuffer = NDArray; + const auto read_into_buffer = [&reader](ReadBuffer &buffer) { + const auto frame_count = reader.next_read_frames(); + const auto frames_read = reader.read_into(buffer.buffer()); + if (frames_read != frame_count) { + throw std::runtime_error( + "MultiThreadedFileReader returned an incomplete batch"); + } + return frames_read; + }; - if (verbose && (i + 1) % progress_interval == 0) { - wait_for_completed(static_cast(i + 1)); - print_progress(static_cast(i + 1)); + const auto fill_batch = [this](ReadBuffer &buffer, + std::size_t frame_count) { + std::vector> views; + views.reserve(frame_count); + auto batch_view = buffer.view(); + for (std::size_t i = 0; i < frame_count; ++i) { + views.push_back(batch_view(i)); + } + std::lock_guard fill_lock(fill_mutex_); + dispatch_fill_batch_(views); + }; + + // Exclude other queue producers for the duration. Drain anything already + // submitted before bypassing the coordinator with direct batch dispatch. + std::lock_guard ingestion_lock(ingestion_mutex_); + flush(); + + std::size_t completed = 0; + if (n_frames != 0) { + // Allocate the maximum wave size once per buffer. The last read may + // contain fewer frames; its returned count limits the views passed to + // the histogram workers, leaving the unused tail untouched. + const auto buffer_capacity = reader.next_read_frames(); + std::array buffers{ + ReadBuffer({static_cast(buffer_capacity), rows_, cols_}), + ReadBuffer({static_cast(buffer_capacity), rows_, cols_})}; + + std::size_t current_index = 0; + auto current_frames = read_into_buffer(buffers[current_index]); + while (true) { + const bool has_next = completed + current_frames < n_frames; + const std::size_t next_index = current_index ^ std::size_t{1}; + + // MultiThreadedFileReader::read_into is blocking. Run the next + // read on a dedicated prefetch task while the current batch is + // processed by the histogram worker pool. + std::future next; + if (has_next) { + next = std::async(std::launch::async, read_into_buffer, + std::ref(buffers[next_index])); + } + + fill_batch(buffers[current_index], current_frames); + completed += current_frames; + + if (verbose && (completed - last_reported >= progress_interval || + completed == n_frames)) { + print_progress(completed); + } + + if (!has_next) { + break; + } + current_frames = next.get(); + current_index = next_index; } } - flush(); + if (verbose) { - const auto done = completed_for_this_file(); - if (done > last_reported) { - print_progress(done); + if (completed > last_reported || n_frames == 0) { + print_progress(completed); } fmt::print("\n\n"); std::fflush(stdout); diff --git a/src/hist/PedestalTrackingPixelHistogram.test.cpp b/src/hist/PedestalTrackingPixelHistogram.test.cpp new file mode 100644 index 00000000..2c9f1590 --- /dev/null +++ b/src/hist/PedestalTrackingPixelHistogram.test.cpp @@ -0,0 +1,289 @@ +// SPDX-License-Identifier: MPL-2.0 +#include "aare/hist/PedestalTrackingPixelHistogram.hpp" + +#include "aare/File.hpp" +#include "aare/Frame.hpp" +#include "aare/NumpyFile.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using aare::FileConfig; +using aare::Frame; +using aare::NumpyFile; +using aare::PedestalTrackingPixelHistogram; + +namespace { + +class TemporaryHistogramFile { + public: + using Generator = + std::function; + + TemporaryHistogramFile(std::size_t rows = 2, std::size_t cols = 3, + std::size_t frames = 10, Generator generator = {}) { + const auto unique = + std::chrono::steady_clock::now().time_since_epoch().count(); + path_ = std::filesystem::temp_directory_path() / + ("aare-pedestal-hist-" + std::to_string(unique) + ".npy"); + + FileConfig config; + config.dtype = aare::Dtype::UINT16; + config.rows = rows; + config.cols = cols; + NumpyFile file(path_, "w", config); + for (std::size_t frame_index = 0; frame_index < frames; ++frame_index) { + Frame frame(rows, cols, config.dtype); + auto image = frame.view(); + for (ssize_t row = 0; row < image.shape(0); ++row) { + for (ssize_t col = 0; col < image.shape(1); ++col) { + image(row, col) = + generator ? generator(frame_index, + static_cast(row), + static_cast(col)) + : static_cast(frame_index); + } + } + file.write(frame); + } + } + + ~TemporaryHistogramFile() { std::filesystem::remove(path_); } + + TemporaryHistogramFile(const TemporaryHistogramFile &) = delete; + TemporaryHistogramFile &operator=(const TemporaryHistogramFile &) = delete; + + const std::filesystem::path &path() const { return path_; } + + private: + std::filesystem::path path_; +}; + +} // namespace + +TEST_CASE("Pedestal tracking histogram fills ordered multi-reader batches", + "[PedestalTrackingPixelHistogram]") { + TemporaryHistogramFile file; + PedestalTrackingPixelHistogram histogram(2, 3, 10, 0.0f, 10.0f, 2, 4, 0.0f); + + // Two reader workers claiming two frames each produces a full batch of + // four followed by a partial batch of three. + histogram.fill_from_file(file.path(), 7, false, 2, 2); + const auto values = histogram.values(); + + for (ssize_t row = 0; row < 2; ++row) { + for (ssize_t col = 0; col < 3; ++col) { + for (ssize_t bin = 0; bin < 10; ++bin) { + CHECK(values(row, col, bin) == (bin < 7 ? 1 : 0)); + } + } + } +} + +TEST_CASE("Pedestal tracking file fill handles limits and reader options", + "[PedestalTrackingPixelHistogram]") { + TemporaryHistogramFile file; + PedestalTrackingPixelHistogram histogram(2, 3, 10, 0.0f, 10.0f, 2, 4, 0.0f); + + CHECK_NOTHROW(histogram.fill_from_file(file.path(), 0, false, 2, 2)); + CHECK_THROWS_AS(histogram.fill_from_file(file.path(), -2, false, 2, 2), + std::invalid_argument); + CHECK_THROWS_AS(histogram.fill_from_file(file.path(), -1, false, 0, 2), + std::invalid_argument); + CHECK_THROWS_AS(histogram.fill_from_file(file.path(), -1, false, 2, 0), + std::invalid_argument); + + // Preserve the previous API's clamp-at-EOF behaviour. + histogram.fill_from_file(file.path(), 100, false, 2, 3); + const auto values = histogram.values(); + for (ssize_t row = 0; row < 2; ++row) { + for (ssize_t col = 0; col < 3; ++col) { + for (ssize_t bin = 0; bin < 10; ++bin) { + CHECK(values(row, col, bin) == 1); + } + } + } +} + +TEST_CASE("Pedestal tracking file fill validates frame metadata", + "[PedestalTrackingPixelHistogram]") { + TemporaryHistogramFile wrong_shape(3, 3, 1); + PedestalTrackingPixelHistogram histogram(2, 3, 10, 0.0f, 10.0f, 1, 4, 0.0f); + + CHECK_THROWS_AS( + histogram.fill_from_file(wrong_shape.path(), -1, false, 2, 1), + std::invalid_argument); +} + +TEST_CASE("Histogram-only fill crosses a pixel tile without changing pedestal", + "[PedestalTrackingPixelHistogram]") { + constexpr int rows = 1; + constexpr int cols = 513; + constexpr std::size_t pedestal_samples = 1000; + + const auto baseline = [](std::size_t col) { + return static_cast(100 + col % 17); + }; + TemporaryHistogramFile file( + rows, cols, 2, + [baseline](std::size_t frame, std::size_t, std::size_t col) { + return static_cast(baseline(col) + frame + 1); + }); + + // Negative n_sigma exercises the complete histogram-only condition; zero + // is covered by the ordered batch tests above. + PedestalTrackingPixelHistogram histogram(rows, cols, 4, 0.0f, 4.0f, 1, 6, + -1.0f); + + for (std::size_t seed = 0; seed < pedestal_samples; ++seed) { + aare::NDArray frame({rows, cols}); + const int offset = seed % 2 == 0 ? -1 : 1; + for (ssize_t col = 0; col < cols; ++col) { + frame(0, col) = static_cast( + static_cast(baseline(static_cast(col))) + + offset); + } + histogram.push_pedestal_no_update(frame.view()); + } + histogram.update_mean(); + const auto mean_before = histogram.pedestal_mean(); + + // One two-frame batch traverses the 512-pixel tile and its one-pixel tail. + histogram.fill_from_file(file.path(), -1, false, 1, 2); + + const auto values = histogram.values(); + for (ssize_t col = 0; col < cols; ++col) { + for (ssize_t bin = 0; bin < 4; ++bin) { + CHECK(values(0, col, bin) == ((bin == 1 || bin == 2) ? 1 : 0)); + } + } + + const auto mean_after = histogram.pedestal_mean(); + CHECK( + std::equal(mean_before.begin(), mean_before.end(), mean_after.begin())); +} + +TEST_CASE("Pedestal tracking threshold is strict and rejects zero sigma", + "[PedestalTrackingPixelHistogram]") { + constexpr std::size_t pedestal_samples = 1000; + PedestalTrackingPixelHistogram histogram(1, 3, 8, -4.0f, 4.0f, 1, 4, 2.0f); + + // Pixels 0 and 1 have mean 100 and population sigma 1. Pixel 2 has the + // same mean and zero sigma. + for (std::size_t seed = 0; seed < pedestal_samples; ++seed) { + aare::NDArray frame({1, 3}); + const auto noisy = static_cast(seed % 2 == 0 ? 99 : 101); + frame(0, 0) = noisy; + frame(0, 1) = noisy; + frame(0, 2) = 100; + histogram.push_pedestal_no_update(frame.view()); + } + histogram.update_mean(); + + const auto mean_before = histogram.pedestal_mean(); + REQUIRE(mean_before(0, 0) == 100.0f); + REQUIRE(mean_before(0, 1) == 100.0f); + REQUIRE(mean_before(0, 2) == 100.0f); + + aare::NDArray frame({1, 3}); + frame(0, 0) = 102; // residual == 2 * sigma: excluded + frame(0, 1) = 101; // residual < 2 * sigma: included + frame(0, 2) = 101; // sigma == 0: excluded + histogram.fill_async(std::move(frame)); + histogram.flush(); + + const auto mean_after = histogram.pedestal_mean(); + CHECK(mean_after(0, 0) == 100.0f); + CHECK(mean_after(0, 1) == static_cast(100001.0 * (1.0 / 1000.0))); + CHECK(mean_after(0, 2) == 100.0f); + + // Histogramming uses each residual before a possible EMA update. + const auto values = histogram.values(); + for (ssize_t bin = 0; bin < 8; ++bin) { + CHECK(values(0, 0, bin) == (bin == 6 ? 1 : 0)); + CHECK(values(0, 1, bin) == (bin == 5 ? 1 : 0)); + CHECK(values(0, 2, bin) == (bin == 5 ? 1 : 0)); + } +} + +TEST_CASE("Tiled pedestal tracking matches chronological single-frame fills", + "[PedestalTrackingPixelHistogram]") { + constexpr int rows = 5; + constexpr int cols = 257; + constexpr std::size_t frames = 40; + constexpr int bins = 32; + constexpr float xmin = -16.0f; + constexpr float xmax = 16.0f; + + const auto baseline = [](std::size_t row, std::size_t col) { + return static_cast(100 + (row + col) % 5); + }; + const auto sample = [baseline](std::size_t frame, std::size_t row, + std::size_t col) { + const int offsets[] = {1, -1, 8, 0}; + return static_cast(static_cast(baseline(row, col)) + + offsets[frame % 4]); + }; + + TemporaryHistogramFile file(rows, cols, frames, sample); + PedestalTrackingPixelHistogram tiled(rows, cols, bins, xmin, xmax, 2, 16, + 2.0f); + PedestalTrackingPixelHistogram chronological(rows, cols, bins, xmin, xmax, + 2, 16, 2.0f); + + // Seed a non-zero cached standard deviation. Five rows split over two + // workers make both row shards cross the 512-pixel tile boundary. Forty + // data frames exercise chronological processing across a sizable batch. + for (std::size_t seed = 0; seed < 1000; ++seed) { + aare::NDArray frame({rows, cols}); + const int offset = seed % 2 == 0 ? -2 : 2; + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + frame(row, col) = static_cast( + static_cast(baseline(row, col)) + offset); + } + } + tiled.push_pedestal_no_update(frame.view()); + chronological.push_pedestal_no_update(frame.view()); + } + tiled.update_mean(); + chronological.update_mean(); + + // One forty-frame reader wave exercises tiled batch processing. + tiled.fill_from_file(file.path(), -1, false, 2, 20); + + // The reference path establishes the same result one chronological frame + // at a time, without any cross-frame traversal reordering. + for (std::size_t frame_index = 0; frame_index < frames; ++frame_index) { + aare::NDArray frame({rows, cols}); + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + frame(row, col) = + sample(frame_index, static_cast(row), + static_cast(col)); + } + } + chronological.fill_async(std::move(frame)); + chronological.flush(); + } + + const auto tiled_values = tiled.values(); + const auto chronological_values = chronological.values(); + REQUIRE(tiled_values.shape() == chronological_values.shape()); + CHECK(std::equal(tiled_values.begin(), tiled_values.end(), + chronological_values.begin())); + + const auto tiled_mean = tiled.pedestal_mean(); + const auto chronological_mean = chronological.pedestal_mean(); + REQUIRE(tiled_mean.shape() == chronological_mean.shape()); + CHECK(std::equal(tiled_mean.begin(), tiled_mean.end(), + chronological_mean.begin())); +} diff --git a/src/hist/PixelHistogramImpl.test.cpp b/src/hist/PixelHistogramImpl.test.cpp index 3cc28424..f3526c17 100644 --- a/src/hist/PixelHistogramImpl.test.cpp +++ b/src/hist/PixelHistogramImpl.test.cpp @@ -72,6 +72,17 @@ TEST_CASE("Fill a small histogram from an NDArray") { REQUIRE(v(1, 1, 4) == 1); } +TEST_CASE("Flat pixel filling uses pixel-major histogram storage") { + aare::PixelHistogramImpl hist(2, 3, 4, 0.0f, 4.0f); + + hist.fill_flat_unchecked(0, 0.5f); + hist.fill_flat_unchecked(4, 2.5f); + + const auto values = hist.view(); + CHECK(values(0, 0, 0) == 1); + CHECK(values(1, 1, 2) == 1); +} + TEST_CASE("Check that pixel histogram does not overflow") { int rows = 1; int cols = 1;