From 406aade69625fdb95884ff19f660ba94f71e1baa Mon Sep 17 00:00:00 2001 From: Erik Frojdh Date: Thu, 8 Oct 2026 15:36:46 +0200 Subject: [PATCH] first attempt at gappixels --- CMakeLists.txt | 3 + RELEASE.md | 15 + benchmarks/CMakeLists.txt | 2 +- docs/src/File.rst | 5 + docs/src/MultiThreadedFileReader.rst | 5 + docs/src/RawFile.rst | 8 + docs/src/index.rst | 1 + .../experimental/MultiThreadedFileReader.rst | 5 + docs/src/python/file/index.rst | 1 + docs/src/python/file/pyFile.rst | 6 + docs/src/python/file/pyRawFile.rst | 6 + include/aare/FileInterface.hpp | 10 +- include/aare/MultiThreadedFileReader.hpp | 11 +- include/aare/RawFile.hpp | 52 ++- python/aare/RawFile.py | 14 +- python/aare/__init__.py | 1 + python/src/bind_MultiThreadedFileReader.hpp | 20 +- python/src/bind_RawFile.hpp | 43 ++- python/src/file.hpp | 39 +- python/src/module.cpp | 2 + python/tests/test_MultiThreadedFileReader.py | 37 ++ python/tests/test_RawFile.py | 147 ++++++++ src/File.cpp | 8 +- src/MultiThreadedFileReader.cpp | 18 +- src/RawFile.cpp | 308 ++++++++-------- src/RawFile.test.cpp | 337 +++++++++++++++++- tests/raw_file_helpers.hpp | 77 +++- 27 files changed, 995 insertions(+), 186 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9a6ba241..62a2a368 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -388,6 +388,7 @@ set(PUBLICHEADERS include/aare/ROI.hpp include/aare/ROIGeometry.hpp include/aare/DetectorGeometry.hpp + include/aare/GapPixels.hpp include/aare/JungfrauDataFile.hpp include/aare/logger.hpp include/aare/NDArray.hpp @@ -413,6 +414,7 @@ set(SourceFiles ${CMAKE_CURRENT_SOURCE_DIR}/src/ROI.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ROIGeometry.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/GapPixels.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PedestalTrackingPixelHistogram.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp @@ -482,6 +484,7 @@ if(AARE_TESTS) ${CMAKE_CURRENT_SOURCE_DIR}/src/fit/Fit.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/GapPixels.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolation.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp diff --git a/RELEASE.md b/RELEASE.md index 86d00ed2..3bf40b65 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -4,6 +4,15 @@ ### New Features: +- Added gap pixel insertion for Jungfrau and Eiger raw files. ``RawFile``, + ``File`` (through ``FileConfig::gap_pixels``) and + ``MultiThreadedFileReader`` take a ``GapPixels`` configuration, in Python + the ``gap_pixels`` argument, and assemble frames with 2 pixel chip gaps + and optional ``ModuleGaps`` between modules. Gap positions follow + detector coordinates before any ROI. Gap pixels receive a fill value or, + with ``split_counts``, half of the neighbouring double size pixel. The + new ``insert_gap_pixels`` and ``gapped_shape`` apply the same transform + to images in memory. - Added the Python ``xy`` type with ``row`` and ``col`` fields, enabling access to ``RawFile.geometry``, ``RawMasterFile.detector_layout``, and ``RawMasterFile.udp_interfaces_per_module``. It unpacks as ``row, col``. @@ -11,6 +20,12 @@ build, so releases publish them to PyPI alongside the Linux and macOS wheels. +### Other: + +- ``RawFile`` assembles multi-module frames through one placement path with + a scratch buffer allocated once at open, instead of allocating a module + buffer for every frame. + ## 2026.9.30 ### New Features: diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index e874a003..f89b57d3 100644 --- a/benchmarks/CMakeLists.txt +++ b/benchmarks/CMakeLists.txt @@ -20,7 +20,7 @@ add_executable(benchmarks) target_sources( benchmarks PRIVATE ndarray_benchmark.cpp ndview_benchmark.cpp calculateeta_benchmark.cpp - reduce_benchmark.cpp clusterfile_benchmark.cpp) + reduce_benchmark.cpp clusterfile_benchmark.cpp rawfile_benchmark.cpp) # Link Google Benchmark and other necessary libraries target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core diff --git a/docs/src/File.rst b/docs/src/File.rst index 784a2722..0fa1ddc4 100644 --- a/docs/src/File.rst +++ b/docs/src/File.rst @@ -1,6 +1,11 @@ File ============= +``File`` opens ``.json`` and ``.raw`` master files, ``.npy`` files and +Jungfrau ``.dat`` files behind one interface. ``FileConfig::gap_pixels`` +inserts :doc:`gap pixels ` when reading raw files; other file +types reject a set value with ``std::invalid_argument``. + .. doxygenclass:: aare::File :members: diff --git a/docs/src/MultiThreadedFileReader.rst b/docs/src/MultiThreadedFileReader.rst index a2bdebf8..63765a2a 100644 --- a/docs/src/MultiThreadedFileReader.rst +++ b/docs/src/MultiThreadedFileReader.rst @@ -25,6 +25,11 @@ an allocation when the caller already owns a buffer of at least remaining from the current position, and ``seek()`` to reposition the reader. Call ``close()`` to release all worker file handles early. +The optional fifth argument inserts :doc:`gap pixels ` when +reading raw files. ``rows()``, ``cols()`` and ``bytes_per_frame()`` then +include the gaps. Each worker reader uses the configured seed plus its +worker index for split counts. + .. note:: Multiple workers do not guarantee faster reads. Performance depends on the diff --git a/docs/src/RawFile.rst b/docs/src/RawFile.rst index 8ecd30e9..db10ae78 100644 --- a/docs/src/RawFile.rst +++ b/docs/src/RawFile.rst @@ -42,6 +42,14 @@ The top-level frame bounds error states the number of frames in the file and that indices are zero-based. +Pass a :doc:`GapPixels ` configuration as the third constructor +argument to insert gap pixels while frames are assembled. ``rows()``, +``cols()``, ``pixels_per_frame()``, ``bytes_per_frame()`` and every returned +frame then include the gaps, for the full frame and for each ROI. The option +is only supported for Jungfrau and Eiger and throws +``std::invalid_argument`` for other detector types. Headers and frame +numbers are unchanged. + .. doxygenclass:: aare::RawFile :members: :undoc-members: diff --git a/docs/src/index.rst b/docs/src/index.rst index 4f0151b7..c869d999 100644 --- a/docs/src/index.rst +++ b/docs/src/index.rst @@ -44,6 +44,7 @@ AARE NDView Frame File + GapPixels MultiThreadedFileReader Dtype Cluster diff --git a/docs/src/python/experimental/MultiThreadedFileReader.rst b/docs/src/python/experimental/MultiThreadedFileReader.rst index d8553e4f..827e217c 100644 --- a/docs/src/python/experimental/MultiThreadedFileReader.rst +++ b/docs/src/python/experimental/MultiThreadedFileReader.rst @@ -24,6 +24,11 @@ iteration position. The context manager closes all worker files on exit. ``close()`` is also available for explicit cleanup and may be called repeatedly. +``gap_pixels`` inserts :doc:`gap pixels <../file/pyGapPixels>` when reading +raw files. Pass ``True`` for the defaults or a ``GapPixels``. ``rows`` and +``cols`` then include the gaps. Each worker uses the configured seed plus its +worker index for split counts. + .. autoclass:: MultiThreadedFileReader :members: :undoc-members: diff --git a/docs/src/python/file/index.rst b/docs/src/python/file/index.rst index 0bcc4f9c..54e3fed2 100644 --- a/docs/src/python/file/index.rst +++ b/docs/src/python/file/index.rst @@ -8,6 +8,7 @@ File I/O pyClusterFile pyCtbRawFile pyFile + pyGapPixels pyJungfrauDataFile pyRawFile pyRawMasterFile diff --git a/docs/src/python/file/pyFile.rst b/docs/src/python/file/pyFile.rst index 5c180d06..28d5e783 100644 --- a/docs/src/python/file/pyFile.rst +++ b/docs/src/python/file/pyFile.rst @@ -4,6 +4,12 @@ File .. py:currentmodule:: aare +``File`` opens raw master files, ``.npy`` files and Jungfrau ``.dat`` files +behind one interface. ``gap_pixels`` inserts +:doc:`gap pixels ` when reading raw files; pass ``True`` for +the defaults or a ``GapPixels``. Other file types raise ``ValueError`` when +it is set. + .. autoclass:: File :members: :undoc-members: diff --git a/docs/src/python/file/pyRawFile.rst b/docs/src/python/file/pyRawFile.rst index a128d2fe..480560e2 100644 --- a/docs/src/python/file/pyRawFile.rst +++ b/docs/src/python/file/pyRawFile.rst @@ -40,6 +40,12 @@ For a file containing five frames, the valid indices are 0 through 4. Attempting to read index 5 raises an error stating: ``Frame index 5 is out of range: file contains 5 frames (indices are zero-based)``. +``gap_pixels`` inserts :doc:`gap pixels ` while frames are +assembled. Pass ``True`` for the defaults or a ``GapPixels``. ``rows()``, +``cols()`` and every returned frame then include the gaps, for the full frame +and for each ROI. Only Jungfrau and Eiger are supported; other detector +types raise ``ValueError``. Headers are unchanged. + .. autoclass:: RawFile :members: :undoc-members: diff --git a/include/aare/FileInterface.hpp b/include/aare/FileInterface.hpp index 87ca3dbf..3219ca4f 100644 --- a/include/aare/FileInterface.hpp +++ b/include/aare/FileInterface.hpp @@ -2,9 +2,11 @@ #pragma once #include "aare/Dtype.hpp" #include "aare/Frame.hpp" +#include "aare/GapPixels.hpp" #include "aare/defs.hpp" #include +#include #include namespace aare { @@ -15,6 +17,7 @@ namespace aare { * rows: number of rows in the file * cols: number of columns in the file * geometry: geometry of the file + * gap_pixels: insert gap pixels when reading raw files */ struct FileConfig { aare::Dtype dtype{typeid(uint16_t)}; @@ -24,7 +27,8 @@ struct FileConfig { return dtype == other.dtype && rows == other.rows && cols == other.cols && geometry == other.geometry && detector_type == other.detector_type && - max_frames_per_file == other.max_frames_per_file; + max_frames_per_file == other.max_frames_per_file && + gap_pixels == other.gap_pixels; } bool operator!=(const FileConfig &other) const { return !(*this == other); } @@ -34,6 +38,10 @@ struct FileConfig { DetectorType detector_type{DetectorType::Unknown}; int max_frames_per_file{}; size_t total_frames{}; + + /// @brief Insert gap pixels into every frame. Only RawFile supports + /// this; other readers reject a set value. + std::optional gap_pixels{}; // std::string to_string() const { // return "{ dtype: " + dtype.to_string() + // ", rows: " + std::to_string(rows) + diff --git a/include/aare/MultiThreadedFileReader.hpp b/include/aare/MultiThreadedFileReader.hpp index d7a930b8..c15040ef 100644 --- a/include/aare/MultiThreadedFileReader.hpp +++ b/include/aare/MultiThreadedFileReader.hpp @@ -3,6 +3,7 @@ #include "aare/Dtype.hpp" #include "aare/File.hpp" +#include "aare/GapPixels.hpp" #include #include @@ -25,10 +26,13 @@ class MultiThreadedFileReader { * @param n_threads maximum number of worker threads * @param chunk_size number of frames claimed by a worker at a time * @param total_frames number of frames to read, or all frames when omitted + * @param gap_pixels insert gap pixels, raw files only. Each worker + * reader uses the configured seed plus its worker index. */ MultiThreadedFileReader(std::filesystem::path fname, size_t n_threads, size_t chunk_size, - std::optional total_frames = std::nullopt); + std::optional total_frames = std::nullopt, + std::optional gap_pixels = std::nullopt); MultiThreadedFileReader(const MultiThreadedFileReader &) = delete; MultiThreadedFileReader & @@ -66,6 +70,9 @@ class MultiThreadedFileReader { bool is_open() const noexcept { return !m_files.empty(); } size_t n_threads() const noexcept { return m_n_threads; } + const std::optional &gap_pixels() const noexcept { + return m_gap_pixels; + } size_t chunk_size() const noexcept { return m_chunk_size; } size_t total_frames() const noexcept { return m_total_frames; } size_t source_total_frames() const noexcept { @@ -96,9 +103,11 @@ class MultiThreadedFileReader { size_t m_bytes_per_frame; size_t m_total_bytes; size_t m_current_frame; + std::optional m_gap_pixels; std::vector m_files; void ensure_open() const; + FileConfig worker_config(size_t worker_index) const; }; } // namespace aare::experimental diff --git a/include/aare/RawFile.hpp b/include/aare/RawFile.hpp index 4ef92962..bda86c7e 100644 --- a/include/aare/RawFile.hpp +++ b/include/aare/RawFile.hpp @@ -2,6 +2,7 @@ #pragma once #include "aare/FileInterface.hpp" #include "aare/Frame.hpp" +#include "aare/GapPixels.hpp" #include "aare/NDArray.hpp" //for pixel map #include "aare/ROIGeometry.hpp" #include "aare/RawMasterFile.hpp" @@ -11,7 +12,9 @@ #include "../tests/friend_test.hpp" #endif +#include #include +#include namespace aare { @@ -29,20 +32,39 @@ class RawFile : public FileInterface { RawMasterFile m_master; size_t m_current_frame{}; - std::vector m_ROI_geometries; - /// @brief Minimum frame count across the selected raw subfile series. size_t m_frames_in_file{}; + /// @brief Scratch space for one module part that is not contiguous in + /// the assembled frame. Sized at open to the largest part in any ROI. + std::vector m_part_buffer; + + /// @brief Gap pixel configuration, empty when reading without gaps + std::optional m_gap_pixels; + detail::GapLayout m_gap_layout{}; + std::mt19937_64 m_generator; + + /// @brief Detector coordinates of each ROI, in roi_geometries() order + std::vector m_roi_rects; + /// @brief {rows, cols} of the assembled frame of each ROI, including + /// gap pixels when enabled + std::vector> m_roi_shapes; + public: /** * @brief RawFile constructor * @param fname path to the master file (.json) * @param mode file mode (only "r" is supported at the moment) + * @param gap_pixels insert gap pixels into every assembled frame. Frame + * sizes, rows() and cols() then include the gaps. Only supported for + * Jungfrau and Eiger. * @throws std::runtime_error if frame padding is disabled and the frame * discard policy is not DiscardPartial. + * @throws std::invalid_argument if gap pixels are requested for another + * detector type or the gap configuration is invalid */ - RawFile(const std::filesystem::path &fname, const std::string &mode = "r"); + RawFile(const std::filesystem::path &fname, const std::string &mode = "r", + std::optional gap_pixels = std::nullopt); virtual ~RawFile() override = default; Frame read_frame() override; @@ -102,6 +124,10 @@ class RawFile : public FileInterface { */ size_t pixels_per_frame(const size_t roi_index); size_t bytes_per_pixel() const; + + /// @brief Gap pixel configuration, empty when reading without gaps + const std::optional &gap_pixels() const { return m_gap_pixels; } + void seek(size_t frame_index) override; size_t tell() override; /// @brief Minimum actual frame count across all subfiles and ROIs. @@ -154,6 +180,26 @@ class RawFile : public FileInterface { std::runtime_error frame_error(size_t frame_index, const std::string &message) const; + /** + * @brief Frame index to read from each module part of a ROI so that all + * parts hold the same detector frame number. Parts whose frame number + * lags are advanced; the index never moves backwards. + * @throws std::runtime_error if a part runs out of frames while + * synchronizing + */ + std::vector synchronized_frame_indices(size_t frame_index, + size_t roi_index); + + /** + * @brief Assemble a frame with gap pixels: fill the gap lines, scatter + * every module part through the scratch buffer and, if configured, + * split the counts of the chip edge pixels. + */ + template + void assemble_with_gaps(std::byte *frame_buffer, size_t roi_index, + const std::vector &frame_indices, + DetectorHeader *header); + /** * @brief read the frame at the given frame index into the image buffer * @param frame_index frame number to read diff --git a/python/aare/RawFile.py b/python/aare/RawFile.py index 253fcc53..56eccf46 100644 --- a/python/aare/RawFile.py +++ b/python/aare/RawFile.py @@ -4,8 +4,18 @@ import numpy as np from .ScanParameters import ScanParameters class RawFile(_aare.RawFile): - def __init__(self, fname, chunk_size = 1): - super().__init__(fname) + def __init__(self, fname, chunk_size = 1, gap_pixels = False): + """Open a raw file for reading. + + Args: + fname: Master file (.json or .raw). + chunk_size: Frames returned per iteration step. + gap_pixels: Insert gap pixels into every frame. True uses the + default GapPixels configuration, or pass a GapPixels. Frame + shapes then include the gaps. Only Jungfrau and Eiger are + supported. + """ + super().__init__(fname, gap_pixels) self._chunk_size = chunk_size diff --git a/python/aare/__init__.py b/python/aare/__init__.py index ea3f4c39..6394489c 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -25,6 +25,7 @@ from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarCluster from ._aare import DetectorType, FrameDiscardPolicy, ReadoutMode from ._aare import hitmap from ._aare import ROI +from ._aare import GapPixels, ModuleGaps, insert_gap_pixels, gapped_shape from ._aare import xy from ._aare import corner diff --git a/python/src/bind_MultiThreadedFileReader.hpp b/python/src/bind_MultiThreadedFileReader.hpp index 0f39cb6c..425fad77 100644 --- a/python/src/bind_MultiThreadedFileReader.hpp +++ b/python/src/bind_MultiThreadedFileReader.hpp @@ -1,6 +1,7 @@ // SPDX-License-Identifier: MPL-2.0 #pragma once +#include "aare/GapPixels.hpp" #include "aare/MultiThreadedFileReader.hpp" #include @@ -79,9 +80,21 @@ inline void define_multi_threaded_file_reader_bindings(py::module_ &m) { reader.attr("__module__") = "aare.experimental"; reader .def(py::init>(), + std::optional, std::optional>(), py::arg("fname"), py::arg("n_threads"), py::arg("chunk_size"), py::arg("total_frames") = py::none(), + py::arg("gap_pixels") = py::none()) + .def(py::init([](std::filesystem::path fname, size_t n_threads, + size_t chunk_size, std::optional total_frames, + bool gap_pixels) { + return MultiThreadedFileReader( + std::move(fname), n_threads, chunk_size, total_frames, + gap_pixels ? std::optional(GapPixels{}) + : std::nullopt); + }), + py::arg("fname"), py::arg("n_threads"), py::arg("chunk_size"), + py::arg("total_frames") = py::none(), + py::arg("gap_pixels") = false, R"doc( Read chunks of detector frames concurrently. @@ -93,7 +106,12 @@ inline void define_multi_threaded_file_reader_bindings(py::module_ &m) { n_threads: Maximum number of worker threads. chunk_size: Number of frames read per claimed chunk. total_frames: Optional frame limit. None reads all frames. + gap_pixels: Insert gap pixels, raw files only. True uses + the default GapPixels configuration. Each worker + uses the configured seed plus its worker index. )doc") + .def_property_readonly("gap_pixels", + &MultiThreadedFileReader::gap_pixels) .def( "read", [](MultiThreadedFileReader &self) { diff --git a/python/src/bind_RawFile.hpp b/python/src/bind_RawFile.hpp index 27e6511c..f3b2d672 100644 --- a/python/src/bind_RawFile.hpp +++ b/python/src/bind_RawFile.hpp @@ -3,6 +3,7 @@ #include "aare/CtbRawFile.hpp" #include "aare/File.hpp" #include "aare/Frame.hpp" +#include "aare/GapPixels.hpp" #include "aare/RawFile.hpp" #include "aare/RawMasterFile.hpp" #include "aare/RawSubFile.hpp" @@ -13,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -34,7 +36,33 @@ std::runtime_error raw_file_read_error(RawFile &file, void define_raw_file_io_bindings(py::module &m) { py::class_(m, "RawFile") - .def(py::init()) + .def(py::init([](const std::filesystem::path &fname, + const GapPixels &gap_pixels) { + return std::make_unique(fname, "r", gap_pixels); + }), + py::arg("fname"), py::arg("gap_pixels")) + .def(py::init([](const std::filesystem::path &fname, bool gap_pixels) { + return std::make_unique( + fname, "r", + gap_pixels ? std::optional(GapPixels{}) + : std::nullopt); + }), + py::arg("fname"), py::arg("gap_pixels") = false, R"( + Open a raw file for reading. + + Parameters + ---------- + fname : path + Master file (.json or .raw). + gap_pixels : bool or GapPixels + Insert gap pixels into every frame. True uses the default + GapPixels configuration. Frame shapes, rows() and cols() + then include the gaps. Only Jungfrau and Eiger are + supported. + )") + .def_property_readonly("gap_pixels", &RawFile::gap_pixels, R"( + Gap pixel configuration, or None when reading without gaps. + )") .def("read_frame", [](RawFile &self) { if (self.n_modules_in_roi().size() > 1) { @@ -118,8 +146,8 @@ void define_raw_file_io_bindings(py::module &m) { std::vector shape; shape.reserve(2); - shape.push_back(self.roi_geometries(roi_index).pixels_y()); - shape.push_back(self.roi_geometries(roi_index).pixels_x()); + shape.push_back(self.rows(roi_index)); + shape.push_back(self.cols(roi_index)); py::array image = allocate_image_data(self.bytes_per_pixel(), shape); @@ -177,8 +205,8 @@ void define_raw_file_io_bindings(py::module &m) { for (size_t r = 0; r < number_of_ROIs; r++) { std::vector shape; shape.reserve(2); - shape.push_back(self.roi_geometries(r).pixels_y()); - shape.push_back(self.roi_geometries(r).pixels_x()); + shape.push_back(self.rows(r)); + shape.push_back(self.cols(r)); images[r] = allocate_image_data(self.bytes_per_pixel(), shape); @@ -229,9 +257,8 @@ void define_raw_file_io_bindings(py::module &m) { if (n_frames == 0) { throw raw_file_read_error(self, "No frames left in file"); } - std::vector shape{ - n_frames, self.roi_geometries(roi_index).pixels_y(), - self.roi_geometries(roi_index).pixels_x()}; + std::vector shape{n_frames, self.rows(roi_index), + self.cols(roi_index)}; // return headers from all subfiles auto n_mod = self.n_modules_in_roi()[roi_index]; diff --git a/python/src/file.hpp b/python/src/file.hpp index 96960d91..3814832b 100644 --- a/python/src/file.hpp +++ b/python/src/file.hpp @@ -3,6 +3,7 @@ #include "aare/CtbRawFile.hpp" #include "aare/File.hpp" #include "aare/Frame.hpp" +#include "aare/GapPixels.hpp" #include "aare/ROI.hpp" #include "aare/RawFile.hpp" #include "aare/RawMasterFile.hpp" @@ -29,13 +30,37 @@ void define_file_io_bindings(py::module &m) { debug, roundRNumber, detType, version, packetMask); py::class_(m, "File") - .def(py::init([](const std::filesystem::path &fname) { - return File(fname, "r", {}); - })) - .def(py::init( - [](const std::filesystem::path &fname, const std::string &mode) { - return File(fname, mode, {}); - })) + .def(py::init([](const std::filesystem::path &fname, + const std::string &mode, + std::optional gap_pixels) { + FileConfig cfg; + cfg.gap_pixels = std::move(gap_pixels); + return File(fname, mode, cfg); + }), + py::arg("fname"), py::arg("mode") = "r", + py::arg("gap_pixels") = py::none()) + .def(py::init([](const std::filesystem::path &fname, + const std::string &mode, bool gap_pixels) { + FileConfig cfg; + if (gap_pixels) { + cfg.gap_pixels = GapPixels{}; + } + return File(fname, mode, cfg); + }), + py::arg("fname"), py::arg("mode") = "r", + py::arg("gap_pixels") = false, R"( + Open a file for reading. + + Parameters + ---------- + fname : path + Master file (.json, .raw), .npy or .dat file. + mode : str + Only "r" is supported. + gap_pixels : bool or GapPixels + Insert gap pixels when reading raw files. True uses the + default GapPixels configuration. + )") .def(py::init()) diff --git a/python/src/module.cpp b/python/src/module.cpp index b83deaca..0df7fe8d 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -14,6 +14,7 @@ #include "bind_Defs.hpp" #include "bind_Eta.hpp" #include "bind_FastPedestal.hpp" +#include "bind_GapPixels.hpp" #include "bind_Interpolator.hpp" #include "bind_MultiThreadedFileReader.hpp" #include "bind_Pedestal.hpp" @@ -156,6 +157,7 @@ PYBIND11_MODULE(_aare, m) { register_calculate_3x3eta(m); define_defs_bindings(m); + define_gap_pixels_bindings(m); using Sum_index_pair_d = Sum_index_pair; PYBIND11_NUMPY_DTYPE(Sum_index_pair_d, sum, index); diff --git a/python/tests/test_MultiThreadedFileReader.py b/python/tests/test_MultiThreadedFileReader.py index 0dc9c654..0b0215c2 100644 --- a/python/tests/test_MultiThreadedFileReader.py +++ b/python/tests/test_MultiThreadedFileReader.py @@ -3,7 +3,9 @@ import numpy as np import pytest import aare +from aare import GapPixels, RawFile, insert_gap_pixels from aare.experimental import MultiThreadedFileReader +from test_RawFile import write_module_layout @pytest.fixture @@ -145,3 +147,38 @@ def test_invalid_configuration( chunk_size=chunk_size, total_frames=total_frames, ) + + +def test_gap_pixels_on_raw_files(tmp_path): + master_path, expected = write_module_layout(tmp_path, 1, 1, frames=4, + module_rows=512, module_cols=1024) + gaps = GapPixels(fill_value=2) + reader = MultiThreadedFileReader(master_path, n_threads=2, chunk_size=1, gap_pixels=gaps) + assert reader.gap_pixels == gaps + assert (reader.rows, reader.cols) == (514, 1030) + frames = reader.read_all() + assert frames.shape == (4, 514, 1030) + assert (frames == insert_gap_pixels(expected, aare.DetectorType.Jungfrau, gaps)).all() + + reader = MultiThreadedFileReader(master_path, n_threads=2, chunk_size=1, gap_pixels=True) + assert reader.gap_pixels == GapPixels() + assert reader.read_all().shape == (4, 514, 1030) + + reader = MultiThreadedFileReader(master_path, n_threads=2, chunk_size=1) + assert reader.gap_pixels is None + assert reader.read_all().shape == (4, 512, 1024) + + +def test_gap_pixels_split_counts_preserve_totals(tmp_path): + master_path, expected = write_module_layout(tmp_path, 1, 1, frames=4, + module_rows=512, module_cols=1024) + gaps = GapPixels(split_counts=True, seed=1) + frames = MultiThreadedFileReader(master_path, n_threads=2, chunk_size=1, gap_pixels=gaps).read_all() + assert frames.shape == (4, 514, 1030) + assert (frames.sum(axis=(1, 2)) == expected.sum()).all() + + +def test_gap_pixels_rejected_for_numpy_files(frame_file): + path, _ = frame_file + with pytest.raises(ValueError): + MultiThreadedFileReader(path, n_threads=1, chunk_size=1, gap_pixels=True) diff --git a/python/tests/test_RawFile.py b/python/tests/test_RawFile.py index 76988db9..85353393 100644 --- a/python/tests/test_RawFile.py +++ b/python/tests/test_RawFile.py @@ -2,6 +2,7 @@ import pytest import json from aare import File, RawFile, RawSubFile, DetectorType, ROI, UDPPortPosition, xy +from aare import GapPixels, ModuleGaps, insert_gap_pixels import numpy as np @@ -28,6 +29,46 @@ def small_raw_file(tmp_path): return master_path +MODULE_ROWS = 2 +MODULE_COLS = 3 + + +def write_module_layout(directory, layout_rows, layout_cols, frames=2, + module_rows=MODULE_ROWS, module_cols=MODULE_COLS, + detector_type="Jungfrau"): + """Write a raw file series with modules in the given layout. + + Modules are numbered as DetectorGeometry orders them, down each column + first, and every module holds distinct values so misplacement is visible. + Returns the master path and the expected image of frame 0. + """ + master_path = directory / "run_master_0.json" + master_path.write_text(json.dumps({ + "Version": 7.2, + "Detector Type": detector_type, + "Timing Mode": "auto", + "Geometry": {"x": layout_cols, "y": layout_rows}, + "Image Size in bytes": module_rows * module_cols * 2, + "Pixels": {"x": module_cols, "y": module_rows}, + "Max Frames Per File": 1, + "Total Frames": frames, + "Frames in File": frames, + "Frame Padding": 1, + "Frame Discard Policy": "nodiscard", + })) + expected = np.zeros((module_rows * layout_rows, module_cols * layout_cols), dtype=np.uint16) + for module in range(layout_rows * layout_cols): + module_row, module_col = module % layout_rows, module // layout_rows + pixels = ((np.arange(module_rows * module_cols, dtype=np.uint32) + module * 10 + 1) % 65536).astype(np.uint16) + expected[module_row * module_rows:(module_row + 1) * module_rows, + module_col * module_cols:(module_col + 1) * module_cols] = pixels.reshape(module_rows, module_cols) + for index in range(frames): + data_path = directory / f"run_d{module}_f{index}_0.raw" + # 112 byte detector header followed by the module pixels + data_path.write_bytes(bytes(112) + pixels.tobytes()) + return master_path, expected + + @pytest.mark.parametrize("reader_type", [RawFile, File]) @pytest.mark.parametrize("legacy_master", [False, True]) @pytest.mark.parametrize("padding, policy, supported", [ @@ -453,3 +494,109 @@ def test_read_eiger_udp_port_disabled(test_data_path): assert len(rois) == 2 assert rois[0] == ROI(0, 512, 0, 512) assert rois[1] == ROI(1024, 1536, 0, 512) + + +@pytest.mark.parametrize("layout", [(1, 1), (2, 1), (1, 2), (2, 2), (3, 2)]) +@pytest.mark.parametrize("reader_type", [RawFile, File]) +def test_raw_file_assembles_module_layout(tmp_path, layout, reader_type): + layout_rows, layout_cols = layout + master_path, expected = write_module_layout(tmp_path, layout_rows, layout_cols) + + with reader_type(master_path) as reader: + # RawFile binds rows/cols as methods (ROI overloads), File as properties + rows = reader.rows() if reader_type is RawFile else reader.rows + cols = reader.cols() if reader_type is RawFile else reader.cols + assert (rows, cols) == expected.shape + result = reader.read_frame() + image = result[1] if reader_type is RawFile else result + assert image.shape == expected.shape + assert (image == expected).all() + + reader.seek(0) + result = reader.read_n(2) + images = result[1] if reader_type is RawFile else result + assert images.shape == (2,) + expected.shape + assert (images == expected).all() + + +@pytest.mark.parametrize("layout", [(1, 1), (2, 1), (1, 2)]) +def test_raw_file_gap_pixels(tmp_path, layout): + layout_rows, layout_cols = layout + master_path, expected = write_module_layout(tmp_path, layout_rows, layout_cols, + module_rows=512, module_cols=1024) + gaps = GapPixels(fill_value=5) + with RawFile(master_path, gap_pixels=gaps) as reader: + assert reader.gap_pixels == gaps + assert (reader.rows(), reader.cols()) == (512 * layout_rows + 2 * (layout_rows * 2 - 1), + 1024 * layout_cols + 2 * (layout_cols * 4 - 1)) + header, image = reader.read_frame() + assert header.size == layout_rows * layout_cols + assert (image == insert_gap_pixels(expected, DetectorType.Jungfrau, gaps)).all() + reader.seek(0) + _, images = reader.read_n(2) + assert images.shape == (2, reader.rows(), reader.cols()) + assert (images[1] == image).all() + + with RawFile(master_path) as reader: + assert reader.gap_pixels is None + assert (reader.rows(), reader.cols()) == expected.shape + + with RawFile(master_path, gap_pixels=True) as reader: + assert reader.gap_pixels == GapPixels() + _, image = reader.read_frame() + assert (image == insert_gap_pixels(expected, DetectorType.Jungfrau)).all() + + with File(master_path, gap_pixels=True) as reader: + assert (reader.rows, reader.cols) == insert_gap_pixels(expected, DetectorType.Jungfrau).shape + assert (reader.read_frame() == insert_gap_pixels(expected, DetectorType.Jungfrau)).all() + + module_gaps = GapPixels(module_gaps=ModuleGaps(8, 36)) + with File(master_path, gap_pixels=module_gaps) as reader: + assert (reader.read_frame() == insert_gap_pixels(expected, DetectorType.Jungfrau, module_gaps)).all() + + +def test_raw_file_gap_pixels_split_counts(tmp_path): + master_path, expected = write_module_layout(tmp_path, 1, 1, module_rows=512, module_cols=1024) + gaps = GapPixels(split_counts=True, seed=11) + with RawFile(master_path, gap_pixels=gaps) as reader: + _, image = reader.read_frame() + assert (image == insert_gap_pixels(expected, DetectorType.Jungfrau, gaps)).all() + assert image.sum() == expected.sum() + + +def test_raw_file_gap_pixels_unsupported_detector(tmp_path): + master_path, _ = write_module_layout(tmp_path, 1, 1, detector_type="Moench") + RawFile(master_path) + with pytest.raises(ValueError): + RawFile(master_path, gap_pixels=True) + with pytest.raises(ValueError): + File(master_path, gap_pixels=True) + + +@pytest.mark.withdata +def test_raw_file_gap_pixels_recorded_jungfrau(test_data_path): + path = test_data_path / "raw/jungfrau/jungfrau_single_master_0.json" + gaps = GapPixels(fill_value=3) + with RawFile(path, gap_pixels=gaps) as gapped, RawFile(path) as plain: + assert (gapped.rows(), gapped.cols()) == (514, 1030) + _, image = gapped.read_frame() + _, reference = plain.read_frame() + assert image.shape == (514, 1030) + assert (image == insert_gap_pixels(reference, DetectorType.Jungfrau, gaps)).all() + + +@pytest.mark.withdata +def test_raw_file_gap_pixels_recorded_rois(test_data_path): + path = test_data_path / "raw/ROITestData/MultipleROIs/run_master_0.json" + with RawFile(path, gap_pixels=True) as gapped, RawFile(path) as plain: + assert gapped.num_rois == 2 + _, images = gapped.read_rois() + _, references = plain.read_rois() + rois = plain.master.rois + assert images[0].shape == (303, 103) + assert images[1].shape == (101, 103) + for image, reference, roi in zip(images, references, rois): + assert (image == insert_gap_pixels(reference, DetectorType.Jungfrau, roi=roi)).all() + gapped.seek(0) + _, frames = gapped.read_n_with_roi(2, roi_index=1) + assert frames.shape == (2, 101, 103) diff --git a/src/File.cpp b/src/File.cpp index 53cc158c..de71660b 100644 --- a/src/File.cpp +++ b/src/File.cpp @@ -23,10 +23,12 @@ File::File(const std::filesystem::path &fname, const std::string &mode, // Assuming we are pointing at a master file? // TODO! How do we read raw files directly? if (fname.extension() == ".raw" || fname.extension() == ".json") { - // file_impl = new RawFile(fname, mode, cfg); - file_impl = std::make_unique(fname, mode); + file_impl = std::make_unique(fname, mode, cfg.gap_pixels); + } else if (cfg.gap_pixels) { + throw std::invalid_argument(fmt::format( + "Gap pixels are only supported for raw files, cannot open '{}'", + fname.string())); } else if (fname.extension() == ".npy") { - // file_impl = new NumpyFile(fname, mode, cfg); file_impl = std::make_unique(fname, mode, cfg); } else if (fname.extension() == ".dat") { file_impl = std::make_unique(fname); diff --git a/src/MultiThreadedFileReader.cpp b/src/MultiThreadedFileReader.cpp index 6810e9bf..ca210922 100644 --- a/src/MultiThreadedFileReader.cpp +++ b/src/MultiThreadedFileReader.cpp @@ -24,11 +24,12 @@ size_t checked_product(size_t lhs, size_t rhs) { MultiThreadedFileReader::MultiThreadedFileReader( std::filesystem::path fname, size_t n_threads, size_t chunk_size, - std::optional total_frames) + std::optional total_frames, std::optional gap_pixels) : m_fname(std::move(fname)), m_n_threads(n_threads), m_chunk_size(chunk_size), m_total_frames(0), m_source_total_frames(0), m_rows(0), m_cols(0), m_bitdepth(0), m_dtype(Dtype::NONE), - m_bytes_per_frame(0), m_total_bytes(0), m_current_frame(0) { + m_bytes_per_frame(0), m_total_bytes(0), m_current_frame(0), + m_gap_pixels(std::move(gap_pixels)) { if (m_n_threads == 0) { throw std::invalid_argument( "MultiThreadedFileReader requires at least one thread"); @@ -38,7 +39,7 @@ MultiThreadedFileReader::MultiThreadedFileReader( "MultiThreadedFileReader chunk size must be greater than zero"); } - File file(m_fname); + File file(m_fname, "r", worker_config(0)); m_source_total_frames = file.total_frames(); m_total_frames = total_frames.value_or(m_source_total_frames); if (m_total_frames > m_source_total_frames) { @@ -56,10 +57,19 @@ MultiThreadedFileReader::MultiThreadedFileReader( m_files.reserve(m_n_threads); m_files.push_back(std::move(file)); for (size_t i = 1; i < m_n_threads; ++i) { - m_files.emplace_back(m_fname); + m_files.emplace_back(m_fname, "r", worker_config(i)); } } +FileConfig MultiThreadedFileReader::worker_config(size_t worker_index) const { + FileConfig config; + config.gap_pixels = m_gap_pixels; + if (config.gap_pixels && config.gap_pixels->seed) { + *config.gap_pixels->seed += worker_index; + } + return config; +} + size_t MultiThreadedFileReader::remaining_frames() const noexcept { return m_total_frames - m_current_frame; } diff --git a/src/RawFile.cpp b/src/RawFile.cpp index e416a27e..f12a2a55 100644 --- a/src/RawFile.cpp +++ b/src/RawFile.cpp @@ -9,6 +9,7 @@ #include "aare/logger.hpp" #include +#include #include #include @@ -23,7 +24,8 @@ std::runtime_error RawFile::frame_error(size_t frame_index, frame_index, m_master.master_fname().string(), message)); } -RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode) +RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode, + std::optional gap_pixels) : m_master(fname) { m_mode = mode; @@ -49,11 +51,37 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode) std::optional min_frames; size_t max_frames = 0; + size_t max_part_bytes = 0; for (const auto &subfiles : m_subfiles) { for (const auto &subfile : subfiles) { const auto count = subfile->frames_in_file(); min_frames = min_frames ? std::min(*min_frames, count) : count; max_frames = std::max(max_frames, count); + max_part_bytes = + std::max(max_part_bytes, subfile->bytes_per_frame()); + } + } + m_part_buffer.resize(max_part_bytes); + + m_roi_rects = m_master.rois(); + if (m_roi_rects.size() != num_rois) { + throw std::runtime_error( + LOCATION + "ROI count differs from ROI geometry count"); + } + m_roi_shapes.reserve(num_rois); + if (gap_pixels) { + m_gap_pixels = gap_pixels; + m_gap_layout = detail::gap_layout( + m_master.detector_type(), m_master.quad() == 1, *gap_pixels); + m_generator = detail::make_generator(gap_pixels->seed); + for (const auto &rect : m_roi_rects) { + m_roi_shapes.push_back( + detail::gapped_shape(rect, m_gap_layout)); + } + } else { + for (const auto &roi : m_master.roi_geometries()) { + m_roi_shapes.push_back({static_cast(roi.pixels_y()), + static_cast(roi.pixels_x())}); } } if (!min_frames) { @@ -197,9 +225,7 @@ size_t RawFile::bytes_per_frame() { } size_t RawFile::bytes_per_frame(const size_t roi_index) { - return m_master.roi_geometries().at(roi_index).pixels_x() * - m_master.roi_geometries().at(roi_index).pixels_y() * - m_master.bitdepth() / bits_per_byte; + return pixels_per_frame(roi_index) * m_master.bitdepth() / bits_per_byte; } size_t RawFile::pixels_per_frame() { @@ -212,8 +238,7 @@ size_t RawFile::pixels_per_frame() { } size_t RawFile::pixels_per_frame(const size_t roi_index) { - return m_master.roi_geometries().at(roi_index).pixels_x() * - m_master.roi_geometries().at(roi_index).pixels_y(); + return rows(roi_index) * cols(roi_index); } DetectorType RawFile::detector_type() const { return m_master.detector_type(); } @@ -244,7 +269,7 @@ size_t RawFile::rows() const { return rows(0); } size_t RawFile::rows(const size_t roi_index) const { - return m_master.roi_geometries().at(roi_index).pixels_y(); + return static_cast(m_roi_shapes.at(roi_index)[0]); } size_t RawFile::cols() const { if (m_master.roi_geometries().size() > 1) { @@ -255,7 +280,7 @@ size_t RawFile::cols() const { return cols(0); } size_t RawFile::cols(const size_t roi_index) const { - return m_master.roi_geometries().at(roi_index).pixels_x(); + return static_cast(m_roi_shapes.at(roi_index)[1]); } size_t RawFile::bitdepth() const { return m_master.bitdepth(); } @@ -328,8 +353,8 @@ DetectorHeader RawFile::read_header(const std::filesystem::path &fname) { RawMasterFile RawFile::master() const { return m_master; } Frame RawFile::get_frame(size_t frame_index, const size_t roi_index) { - auto f = Frame(m_master.roi_geometries().at(roi_index).pixels_y(), - m_master.roi_geometries().at(roi_index).pixels_x(), + auto f = Frame(static_cast(rows(roi_index)), + static_cast(cols(roi_index)), Dtype::from_bitdepth(m_master.bitdepth())); std::byte *frame_buffer = f.data(); get_frame_into(frame_index, frame_buffer, roi_index); @@ -338,6 +363,47 @@ Frame RawFile::get_frame(size_t frame_index, const size_t roi_index) { size_t RawFile::bytes_per_pixel() const { return m_master.bitdepth() / 8; } +std::vector RawFile::synchronized_frame_indices(size_t frame_index, + size_t roi_index) { + const auto &roi = m_master.roi_geometries().at(roi_index); + const size_t n_parts = roi.num_modules_in_roi(); + std::vector frame_indices(n_parts, frame_index); + if (n_parts == 1) { + return frame_indices; + } + + std::vector frame_numbers(n_parts); + for (size_t part_idx = 0; part_idx != n_parts; ++part_idx) { + frame_numbers[part_idx] = + m_subfiles[roi_index][part_idx]->frame_number(frame_index); + } + + while (!all_equal(frame_numbers)) { + // advance the part with the lowest frame number until all agree + const auto min_frame_idx = static_cast(std::distance( + frame_numbers.begin(), + std::min_element(frame_numbers.begin(), frame_numbers.end()))); + + frame_indices[min_frame_idx]++; + if (frame_indices[min_frame_idx] >= total_frames()) { + throw frame_error( + frame_index, + LOCATION + + fmt::format( + "Frame index {} out of range while synchronizing " + "module {} for ROI {}; last data file '{}'", + frame_indices[min_frame_idx], min_frame_idx, roi_index, + m_subfiles[roi_index][min_frame_idx] + ->current_path() + .string())); + } + frame_numbers[min_frame_idx] = + m_subfiles[roi_index][min_frame_idx]->frame_number( + frame_indices[min_frame_idx]); + } + return frame_indices; +} + void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, const size_t roi_index, DetectorHeader *header) { LOG(logDEBUG) << "RawFile::get_frame_into(" << frame_index << ")"; @@ -349,148 +415,106 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, "frames (indices are zero-based)", frame_index, total_frames())); } - std::vector frame_numbers( - m_master.roi_geometries().at(roi_index).num_modules_in_roi()); - std::vector frame_indices( - m_master.roi_geometries().at(roi_index).num_modules_in_roi(), - frame_index); - // sync the frame numbers + const auto &roi = m_master.roi_geometries().at(roi_index); + const auto frame_indices = + synchronized_frame_indices(frame_index, roi_index); - if (m_master.roi_geometries().at(roi_index).num_modules_in_roi() != - 1) { // if we have more than one module - for (size_t part_idx = 0; - part_idx != - m_master.roi_geometries().at(roi_index).num_modules_in_roi(); - ++part_idx) { - frame_numbers[part_idx] = - m_subfiles[roi_index][part_idx]->frame_number(frame_index); - } - - // 1. if frame number vector is the same break - while (!all_equal(frame_numbers)) { - - // 2. find the index of the minimum frame number, - auto min_frame_idx = std::distance( - frame_numbers.begin(), - std::min_element(frame_numbers.begin(), frame_numbers.end())); - - // 3. increase its index and update its respective frame - // number - frame_indices[min_frame_idx]++; - - // 4. if we can't increase its index => throw error - if (frame_indices[min_frame_idx] >= total_frames()) { - throw frame_error( - frame_index, - LOCATION + - fmt::format( - "Frame index {} out of range while synchronizing " - "module {} for ROI {}; last data file '{}'", - frame_indices[min_frame_idx], min_frame_idx, - roi_index, - m_subfiles[roi_index][min_frame_idx] - ->current_path() - .string())); - } - - frame_numbers[min_frame_idx] = - m_subfiles[roi_index][min_frame_idx]->frame_number( - frame_indices[min_frame_idx]); + if (m_gap_pixels) { + switch (m_master.bitdepth()) { + case 8: + assemble_with_gaps(frame_buffer, roi_index, frame_indices, + header); + break; + case 16: + assemble_with_gaps(frame_buffer, roi_index, frame_indices, + header); + break; + case 32: + assemble_with_gaps(frame_buffer, roi_index, frame_indices, + header); + break; + default: + throw frame_error(frame_index, + LOCATION + + fmt::format("Gap pixels are not supported " + "for a bit depth of {}", + m_master.bitdepth())); } + return; } - if (m_master.detector_layout().col == 1) { - // get the part from each subfile and copy it to the frame - for (size_t part_idx = 0; - part_idx != - m_master.roi_geometries().at(roi_index).num_modules_in_roi(); - ++part_idx) { - auto corrected_idx = frame_indices[part_idx]; + const size_t bytes_per_pixel = m_master.bitdepth() / bits_per_byte; + const size_t roi_width = roi.pixels_x(); - // This is where we start writing - auto offset = - (m_master.geometry() - .get_module_geometries( - m_master.roi_geometries() - .at(roi_index) - .module_indices_in_roi(part_idx)) - .origin_y * - m_master.roi_geometries().at(roi_index).pixels_x() + - m_master.geometry() - .get_module_geometries( - m_master.roi_geometries() - .at(roi_index) - .module_indices_in_roi(part_idx)) - .origin_x) * - m_master.bitdepth() / 8; + for (size_t part_idx = 0; part_idx != roi.num_modules_in_roi(); + ++part_idx) { + // origin and size of the part relative to the ROI + const auto &pos = m_master.geometry().get_module_geometries( + roi.module_indices_in_roi(part_idx)); + const auto width = static_cast(pos.width); + const auto height = static_cast(pos.height); + auto &subfile = *m_subfiles[roi_index][part_idx]; + subfile.seek(frame_indices[part_idx]); - if (m_master.geometry() - .get_module_geometries(m_master.roi_geometries() - .at(roi_index) - .module_indices_in_roi(part_idx)) - .origin_x != 0) - throw frame_error( - frame_index, - LOCATION + - " Implementation error. x pos not 0."); // TODO: - // origin - // can still - // change if - // roi - // changes - // TODO! What if the files don't match? - m_subfiles[roi_index][part_idx]->seek(corrected_idx); - m_subfiles[roi_index][part_idx]->read_into(frame_buffer + offset, - header); - if (header) - ++header; - } + std::byte *dest = + frame_buffer + (static_cast(pos.origin_y) * roi_width + + static_cast(pos.origin_x)) * + bytes_per_pixel; - } else { - // TODO! should we read row by row? - - // create a buffer large enough to hold a full module - auto bytes_per_part = - m_master.pixels_y() * m_master.pixels_x() * m_master.bitdepth() / - 8; // TODO! replace with image_size_in_bytes // TODO - // shouldnt it only be the module size? - check - - auto *part_buffer = new std::byte[bytes_per_part]; - - // TODO! if we have many submodules we should reorder them on - // the module level - - for (size_t part_idx = 0; - part_idx != - m_master.roi_geometries().at(roi_index).num_modules_in_roi(); - ++part_idx) { - auto pos = m_master.geometry().get_module_geometries( - m_master.roi_geometries().at(roi_index).module_indices_in_roi( - part_idx)); - auto corrected_idx = frame_indices[part_idx]; - - m_subfiles[roi_index][part_idx]->seek(corrected_idx); - m_subfiles[roi_index][part_idx]->read_into(part_buffer, header); - if (header) - ++header; - - for (size_t cur_row = 0; cur_row < static_cast(pos.height); - cur_row++) { - - auto irow = (pos.origin_y + cur_row); - auto icol = pos.origin_x; - auto dest = - (irow * m_master.roi_geometries().at(roi_index).pixels_x() + - icol); - dest = dest * m_master.bitdepth() / 8; - memcpy(frame_buffer + dest, - part_buffer + - cur_row * pos.width * m_master.bitdepth() / 8, - pos.width * m_master.bitdepth() / 8); + if (width == roi_width) { + // the part spans full rows of the ROI and is contiguous in the + // assembled frame + subfile.read_into(dest, header); + } else { + subfile.read_into(m_part_buffer.data(), header); + const size_t row_bytes = width * bytes_per_pixel; + const size_t roi_row_bytes = roi_width * bytes_per_pixel; + for (size_t row = 0; row < height; ++row) { + std::memcpy(dest + row * roi_row_bytes, + m_part_buffer.data() + row * row_bytes, row_bytes); } } - delete[] part_buffer; + if (header) { + ++header; + } + } +} + +template +void RawFile::assemble_with_gaps(std::byte *frame_buffer, size_t roi_index, + const std::vector &frame_indices, + DetectorHeader *header) { + const auto &roi = m_master.roi_geometries().at(roi_index); + const ROI &rect = m_roi_rects.at(roi_index); + NDView destination(reinterpret_cast(frame_buffer), + m_roi_shapes.at(roi_index)); + detail::fill_gap_pixels(destination, rect, m_gap_layout, + detail::fill_as(m_gap_pixels->fill_value)); + + for (size_t part_idx = 0; part_idx != roi.num_modules_in_roi(); + ++part_idx) { + const auto &pos = m_master.geometry().get_module_geometries( + roi.module_indices_in_roi(part_idx)); + auto &subfile = *m_subfiles[roi_index][part_idx]; + subfile.seek(frame_indices[part_idx]); + subfile.read_into(m_part_buffer.data(), header); + if (header) { + ++header; + } + // detector coordinates of the part + const ROI part_rect{ + rect.xmin + pos.origin_x, rect.xmin + pos.origin_x + pos.width, + rect.ymin + pos.origin_y, rect.ymin + pos.origin_y + pos.height}; + NDView source( + reinterpret_cast(m_part_buffer.data()), + {pos.height, pos.width}); + detail::copy_with_gaps(source, part_rect, destination, rect, + m_gap_layout); + } + + if (m_gap_pixels->split_counts) { + detail::split_edge_counts(destination, rect, m_gap_layout, m_generator); } } diff --git a/src/RawFile.test.cpp b/src/RawFile.test.cpp index ee0b9e24..f4dc29e0 100644 --- a/src/RawFile.test.cpp +++ b/src/RawFile.test.cpp @@ -1,6 +1,7 @@ // SPDX-License-Identifier: MPL-2.0 #include "aare/RawFile.hpp" #include "aare/File.hpp" +#include "aare/GapPixels.hpp" #include "aare/RawMasterFile.hpp" //needed for ROI #include @@ -222,6 +223,338 @@ TEST_CASE("RawFile frame count comes from the complete subfile series", REQUIRE_NOTHROW(reader.read_frame(1)); } +TEST_CASE("RawFile assembles every module layout", "[RawFile][geometry]") { + const auto [layout_rows, layout_cols] = GENERATE(table({ + {1, 1}, + {2, 1}, + {1, 2}, + {2, 2}, + {3, 2}, + })); + CAPTURE(layout_rows, layout_cols); + TemporaryRawFiles files(xy{layout_rows, layout_cols}); + const size_t rows = TemporaryRawFiles::module_rows * layout_rows; + const size_t cols = TemporaryRawFiles::module_cols * layout_cols; + + auto check = [&](NDView image) { + REQUIRE(image.shape(0) == static_cast(rows)); + REQUIRE(image.shape(1) == static_cast(cols)); + for (size_t row = 0; row < rows; ++row) { + for (size_t col = 0; col < cols; ++col) { + CAPTURE(row, col); + REQUIRE(image(row, col) == files.expected(row, col)); + } + } + }; + + RawFile file(files.master_path()); + REQUIRE(file.rows() == rows); + REQUIRE(file.cols() == cols); + REQUIRE(file.n_modules() == files.n_modules()); + REQUIRE(file.bytes_per_frame() == rows * cols * sizeof(uint16_t)); + + SECTION("read_frame") { + auto frame = file.read_frame(); + check(frame.view()); + } + SECTION("read_into with headers") { + NDArray image( + {static_cast(rows), static_cast(cols)}); + std::vector headers(files.n_modules()); + file.seek(1); + file.read_into(reinterpret_cast(image.data()), + headers.data()); + check(image.view()); + for (size_t module = 0; module < files.n_modules(); ++module) { + CAPTURE(module); + REQUIRE(headers[module].frameNumber == 101); + REQUIRE(headers[module].row == module % layout_rows); + REQUIRE(headers[module].column == module / layout_rows); + } + } + SECTION("read_n") { + auto frames = file.read_n(2); + REQUIRE(frames.size() == 2); + for (auto &frame : frames) { + check(frame.view()); + } + } + SECTION("generic File") { + File reader(files.master_path()); + auto frame = reader.read_frame(1); + check(frame.view()); + } +} + +namespace { + +/// Compare a frame read with gap pixels against the free function applied +/// to the same frame read without gaps. +template +bool matches_gap_oracle(Frame &gapped, Frame &plain, const ROI &roi, + DetectorType detector, bool quad, + const GapPixels &gaps) { + auto expected = insert_gap_pixels(NDView(plain.view()), roi, + detector, quad, gaps); + auto actual = gapped.view(); + if (expected.shape() != actual.shape()) { + UNSCOPED_INFO("shape " << actual.shape(0) << "x" << actual.shape(1) + << " expected " << expected.shape(0) << "x" + << expected.shape(1)); + return false; + } + for (ssize_t row = 0; row < expected.shape(0); ++row) { + for (ssize_t col = 0; col < expected.shape(1); ++col) { + if (expected(row, col) != actual(row, col)) { + UNSCOPED_INFO("pixel (" << row << ", " << col << ") " + << actual(row, col) << " expected " + << expected(row, col)); + return false; + } + } + } + return true; +} + +bool matches_gap_oracle(Frame &gapped, Frame &plain, const ROI &roi, + DetectorType detector, bool quad, + const GapPixels &gaps) { + switch (plain.bitdepth()) { + case 8: + return matches_gap_oracle(gapped, plain, roi, detector, quad, + gaps); + case 16: + return matches_gap_oracle(gapped, plain, roi, detector, quad, + gaps); + case 32: + return matches_gap_oracle(gapped, plain, roi, detector, quad, + gaps); + default: + return false; + } +} + +} // namespace + +TEST_CASE("RawFile inserts gap pixels into synthetic Jungfrau modules", + "[RawFile][GapPixels]") { + const auto [layout_rows, layout_cols] = GENERATE(table({ + {1, 1}, + {2, 1}, + {1, 2}, + })); + CAPTURE(layout_rows, layout_cols); + TemporaryRawFiles files(xy{layout_rows, layout_cols}, xy{512, 1024}); + const ROI detector{0, static_cast(files.cols()), 0, + static_cast(files.rows())}; + + GapPixels gaps; + SECTION("defaults") {} + SECTION("module gaps and fill value") { + gaps.module_gaps = ModuleGaps{8, 36}; + gaps.fill_value = 7; + } + const auto layout = detail::gap_layout(DetectorType::Jungfrau, false, gaps); + const auto shape = detail::gapped_shape(detector, layout); + const auto fill = static_cast(gaps.fill_value); + + RawFile file(files.master_path(), "r", gaps); + REQUIRE(file.gap_pixels() == gaps); + REQUIRE(file.rows() == static_cast(shape[0])); + REQUIRE(file.cols() == static_cast(shape[1])); + REQUIRE(file.pixels_per_frame() == + static_cast(shape[0] * shape[1])); + REQUIRE(file.bytes_per_frame() == file.pixels_per_frame() * 2); + + std::vector headers(files.n_modules()); + NDArray image(shape); + file.read_into(reinterpret_cast(image.data()), headers.data()); + for (const auto &header : headers) { + REQUIRE(header.frameNumber == 100); + } + + for (size_t row = 0; row < files.rows(); ++row) { + for (size_t col = 0; col < files.cols(); ++col) { + const auto value = + image(layout.y.gapped(static_cast(row)), + layout.x.gapped(static_cast(col))); + if (value != files.expected(row, col)) { + CAPTURE(row, col); + REQUIRE(value == files.expected(row, col)); + } + } + } + for (const auto boundary : + layout.x.inner_boundaries(detector.xmin, detector.xmax)) { + const auto [first, last] = layout.x.gap_range(boundary, 0); + for (ssize_t col = first; col <= last; ++col) { + for (ssize_t row = 0; row < shape[0]; ++row) { + if (image(row, col) != fill) { + CAPTURE(row, col); + REQUIRE(image(row, col) == fill); + } + } + } + } + for (const auto boundary : + layout.y.inner_boundaries(detector.ymin, detector.ymax)) { + const auto [first, last] = layout.y.gap_range(boundary, 0); + for (ssize_t row = first; row <= last; ++row) { + for (ssize_t col = 0; col < shape[1]; ++col) { + if (image(row, col) != fill) { + CAPTURE(row, col); + REQUIRE(image(row, col) == fill); + } + } + } + } + + auto frame = file.read_frame(); + REQUIRE(frame.rows() == static_cast(shape[0])); + REQUIRE(frame.cols() == static_cast(shape[1])); + REQUIRE(frame.view() == image.view()); +} + +TEST_CASE("RawFile splits counts like the free function", + "[RawFile][GapPixels]") { + TemporaryRawFiles files(xy{1, 2}, xy{512, 1024}); + GapPixels gaps; + gaps.split_counts = true; + gaps.seed = 7; + RawFile gapped(files.master_path(), "r", gaps); + RawFile plain(files.master_path()); + auto gapped_frame = gapped.read_frame(); + auto plain_frame = plain.read_frame(); + const ROI detector{0, 2048, 0, 512}; + REQUIRE(matches_gap_oracle(gapped_frame, plain_frame, detector, + DetectorType::Jungfrau, false, gaps)); + + // the next frame continues the generator, so the oracle for frame 1 is + // the free function applied after the same number of draws + auto frames = gapped.read_n(1); + auto second_plain = plain.read_frame(); + auto expected = insert_gap_pixels( + NDView(second_plain.view()), detector, + DetectorType::Jungfrau, false, gaps); + double gapped_total = 0; + double expected_total = 0; + auto view = frames[0].view(); + for (ssize_t row = 0; row < expected.shape(0); ++row) { + for (ssize_t col = 0; col < expected.shape(1); ++col) { + gapped_total += view(row, col); + expected_total += expected(row, col); + } + } + REQUIRE(gapped_total == expected_total); +} + +TEST_CASE("RawFile rejects gap pixels for unsupported detectors", + "[RawFile][GapPixels]") { + TemporaryRawFiles files; + auto metadata = nlohmann::json::parse(std::ifstream(files.master_path())); + metadata["Detector Type"] = "Moench"; + std::ofstream(files.master_path()) << metadata; + REQUIRE_NOTHROW(RawFile(files.master_path())); + REQUIRE_THROWS_AS(RawFile(files.master_path(), "r", GapPixels{}), + std::invalid_argument); +} + +TEST_CASE("Gap pixels on recorded Jungfrau and Eiger files", + "[.with-data][RawFile][GapPixels]") { + struct Case { + const char *path; + std::array shape; + std::optional module_gaps; + }; + const auto test = GENERATE(values({ + {"raw/jungfrau/jungfrau_single_master_0.json", {514, 1030}, {}}, + {"raw/jungfrau/jungfrau_double_master_0.json", {514, 1030}, {}}, + {"raw/eiger/Lab6_20500eV_2deg_20240629_master_7.json", {514, 1030}, {}}, + {"raw/eiger/eiger_500k_16bit_master_0.json", {514, 1030}, {}}, + {"raw/eiger_quad_data/" + "W13_vrpreampscan_m21C_300V_800eV_vthre2000_master_0.json", + {514, 514}, + {}}, + {"raw/virtual/920/jungfrau/two_modules_master_0.json", + {1030, 1030}, + {}}, + {"raw/virtual/920/jungfrau/two_modules_master_0.json", + {1064, 1030}, + ModuleGaps{8, 36}}, + {"raw/ROITestData/MultipleChipROI/run_master_2.json", {303, 103}, {}}, + {"raw/ROITestData/MultipleChipROI/run_master_2.json", + {337, 103}, + ModuleGaps{8, 36}}, + {"raw/virtual/920/jungfrau/" + "single_module_roi_512_767_0_255_master_0.json", + {256, 256}, + {}}, + })); + CAPTURE(test.path); + const auto fpath = test_data_path() / test.path; + REQUIRE(std::filesystem::exists(fpath)); + + GapPixels gaps; + gaps.module_gaps = test.module_gaps; + gaps.fill_value = 3; + RawFile gapped(fpath, "r", gaps); + RawFile plain(fpath); + REQUIRE(gapped.rows() == static_cast(test.shape[0])); + REQUIRE(gapped.cols() == static_cast(test.shape[1])); + REQUIRE(gapped.bytes_per_frame() == + gapped.pixels_per_frame() * gapped.bytes_per_pixel()); + + auto gapped_frame = gapped.read_frame(); + auto plain_frame = plain.read_frame(); + REQUIRE(gapped_frame.rows() == static_cast(test.shape[0])); + REQUIRE(gapped_frame.cols() == static_cast(test.shape[1])); + REQUIRE(matches_gap_oracle(gapped_frame, plain_frame, plain.master().roi(), + plain.detector_type(), + plain.master().quad() == 1, gaps)); +} + +TEST_CASE("Split counts on a recorded Jungfrau file", + "[.with-data][RawFile][GapPixels]") { + const auto fpath = + test_data_path() / "raw/jungfrau/jungfrau_single_master_0.json"; + GapPixels gaps; + gaps.split_counts = true; + gaps.seed = 2024; + RawFile gapped(fpath, "r", gaps); + RawFile plain(fpath); + auto gapped_frame = gapped.read_frame(); + auto plain_frame = plain.read_frame(); + REQUIRE(matches_gap_oracle(gapped_frame, plain_frame, plain.master().roi(), + DetectorType::Jungfrau, false, gaps)); +} + +TEST_CASE("Gap pixels on a file with two ROIs", + "[.with-data][RawFile][GapPixels]") { + const auto fpath = + test_data_path() / "raw/ROITestData/MultipleROIs/run_master_0.json"; + GapPixels gaps; + RawFile gapped(fpath, "r", gaps); + RawFile plain(fpath); + REQUIRE(gapped.num_rois() == 2); + // x 200..300 holds the chip boundary at 256; y 100..400 holds the one at + // 256 while y 600..700 holds none + REQUIRE(gapped.rows(0) == 303); + REQUIRE(gapped.cols(0) == 103); + REQUIRE(gapped.rows(1) == 101); + REQUIRE(gapped.cols(1) == 103); + + auto gapped_frames = gapped.read_rois(); + auto plain_frames = plain.read_rois(); + const auto rois = plain.master().rois(); + REQUIRE(gapped_frames.size() == 2); + for (size_t index = 0; index < 2; ++index) { + CAPTURE(index); + REQUIRE(matches_gap_oracle(gapped_frames[index], plain_frames[index], + rois[index], DetectorType::Jungfrau, false, + gaps)); + } +} + TEST_CASE("Read number of frames from a jungfrau raw file", "[.with-data][RawFile]") { @@ -535,7 +868,7 @@ TEST_CASE("Open multi module file with ROI", } TEST_CASE("Open multi module file with ROI spanning over multiple modules", - "[.width-data][read_rois][RawFile]") { + "[.with-data][read_rois][RawFile]") { auto fpath = test_data_path() / "raw/ROITestData/MultipleChipROI/run_master_2.json"; @@ -557,7 +890,7 @@ TEST_CASE("Open multi module file with ROI spanning over multiple modules", } TEST_CASE("Open multi module file with two ROIs", - "[.width-data][read_rois][RawFile]") { + "[.with-data][read_rois][RawFile]") { auto fpath = test_data_path() / "raw/ROITestData/MultipleROIs/run_master_0.json"; diff --git a/tests/raw_file_helpers.hpp b/tests/raw_file_helpers.hpp index b1a9b17c..4b8733e0 100644 --- a/tests/raw_file_helpers.hpp +++ b/tests/raw_file_helpers.hpp @@ -8,44 +8,99 @@ #include #include #include +#include +/** + * @brief Throwaway raw file series with a 2x3 pixel Jungfrau "module" per + * data file. Modules are laid out as DetectorGeometry orders them (module + * index runs down each column first) and every module holds distinct pixel + * values so that misplaced modules are detected. + */ class TemporaryRawFiles { std::filesystem::path m_directory; + aare::xy m_layout; + aare::xy m_module_size; public: - explicit TemporaryRawFiles(size_t modules = 1) { + static constexpr size_t module_rows = 2; + static constexpr size_t module_cols = 3; + static constexpr size_t frames = 2; + + /// Value written to pixel index `pixel` (row major) of module `module` + static uint16_t pixel_value(size_t module, size_t pixel) { + return static_cast((module * 10 + pixel + 1) % 65536); + } + + /// Modules stacked in one column + explicit TemporaryRawFiles(size_t modules = 1) + : TemporaryRawFiles(aare::xy{static_cast(modules), 1}) {} + + /// Modules of `module_size` pixels arranged in `layout.row` rows and + /// `layout.col` columns. The default module is 2x3 pixels; pass + /// {512, 1024} for a full Jungfrau module. + explicit TemporaryRawFiles(aare::xy layout, + aare::xy module_size = {module_rows, + module_cols}) + : m_layout(layout), m_module_size(module_size) { const auto unique = std::chrono::steady_clock::now().time_since_epoch().count(); m_directory = std::filesystem::temp_directory_path() / ("aare-raw-" + std::to_string(unique)); std::filesystem::create_directory(m_directory); + const size_t pixels_per_module = module_size.row * module_size.col; const nlohmann::json metadata = { {"Version", 7.2}, {"Detector Type", "Jungfrau"}, {"Timing Mode", "auto"}, - {"Geometry", {{"x", 1}, {"y", modules}}}, - {"Image Size in bytes", 12}, - {"Pixels", {{"x", 3}, {"y", 2}}}, + {"Geometry", {{"x", layout.col}, {"y", layout.row}}}, + {"Image Size in bytes", pixels_per_module * 2}, + {"Pixels", {{"x", module_size.col}, {"y", module_size.row}}}, {"Max Frames Per File", 1}, - {"Total Frames", 2}, - {"Frames in File", 2}, + {"Total Frames", frames}, + {"Frames in File", frames}, {"Frame Padding", 1}, {"Frame Discard Policy", "nodiscard"}}; std::ofstream(master_path()) << metadata; - for (size_t module = 0; module < modules; ++module) { - for (size_t file = 0; file < 2; ++file) { + std::vector pixels(pixels_per_module); + for (size_t module = 0; module < n_modules(); ++module) { + for (size_t pixel = 0; pixel < pixels_per_module; ++pixel) { + pixels[pixel] = pixel_value(module, pixel); + } + for (size_t file = 0; file < frames; ++file) { std::ofstream output(data_path(module, file), std::ios::binary); aare::DetectorHeader header{}; header.frameNumber = 100 + file; - const uint16_t pixels[6] = {1, 2, 3, 4, 5, 6}; + header.row = static_cast(module % layout.row); + header.column = static_cast(module / layout.row); output.write(reinterpret_cast(&header), sizeof(header)); - output.write(reinterpret_cast(pixels), - sizeof(pixels)); + output.write(reinterpret_cast(pixels.data()), + static_cast(pixels.size() * + sizeof(uint16_t))); } } } + aare::xy layout() const { return m_layout; } + aare::xy module_size() const { return m_module_size; } + size_t n_modules() const { return m_layout.row * m_layout.col; } + size_t rows() const { return m_module_size.row * m_layout.row; } + size_t cols() const { return m_module_size.col * m_layout.col; } + + /// Module index (data file number) of the module at a layout position + size_t module_index(size_t module_row, size_t module_col) const { + return module_col * m_layout.row + module_row; + } + + /// Expected value of pixel (row, col) in the assembled image + uint16_t expected(size_t row, size_t col) const { + const size_t module = + module_index(row / m_module_size.row, col / m_module_size.col); + const size_t pixel = (row % m_module_size.row) * m_module_size.col + + col % m_module_size.col; + return pixel_value(module, pixel); + } + TemporaryRawFiles(const TemporaryRawFiles &) = delete; TemporaryRawFiles &operator=(const TemporaryRawFiles &) = delete; ~TemporaryRawFiles() { std::filesystem::remove_all(m_directory); }