diff --git a/CMakeLists.txt b/CMakeLists.txt index 5f73d505..d84db06e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -420,6 +420,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 @@ -445,6 +446,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 @@ -514,6 +516,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 6d2e843a..d9fcaa23 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -10,6 +10,30 @@ - Added Windows wheels (``win_amd64``, Python 3.12 to 3.14) to the wheel build, so releases publish them to PyPI alongside the Linux and macOS wheels. +- added ``NDView::sub_view(index)`` operator returns a subview at index along dimension 0 +- Added ``GapPixels`` with ``insert_gap_pixels`` and ``gapped_shape`` to + insert gap pixels into Jungfrau and Eiger images in memory: 2 pixel gaps + at the chip boundaries inside each module and, only when ``ModuleGaps`` + is set, gaps 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. +- ``RawFile`` takes a ``GapPixels`` configuration, in Python the + ``gap_pixels`` argument, and inserts gap pixels while frames are + assembled. ``rows()``, ``cols()`` and every returned frame then include + the gaps, for the full frame and for each ROI. + +### 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. + +### Bugfixes: + +- Fixed ``Matterhorn10Transform`` with 4 bit dynamic range and + ``Mythen302Transform`` failing with ``Expected 1D array, got 2D`` on frames + read through ``CtbRawFile``. The transforms now flatten the frame before + calling the decoders. ### Build: diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index fab3fa19..84165139 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/benchmarks/rawfile_benchmark.cpp b/benchmarks/rawfile_benchmark.cpp new file mode 100644 index 00000000..cecd5183 --- /dev/null +++ b/benchmarks/rawfile_benchmark.cpp @@ -0,0 +1,112 @@ +// SPDX-License-Identifier: MPL-2.0 +#include "aare/GapPixels.hpp" +#include "aare/RawFile.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +/// One Jungfrau module, 512x1024 pixels at 16 bit, with every frame in a +/// single data file written to the temporary directory. +class RawFileFixture : public benchmark::Fixture { + protected: + static constexpr size_t rows = 512; + static constexpr size_t cols = 1024; + static constexpr size_t frame_count = 32; + std::filesystem::path m_directory; + + public: + void SetUp(const benchmark::State &) override { + const auto stamp = + std::chrono::steady_clock::now().time_since_epoch().count(); + m_directory = std::filesystem::temp_directory_path() / + ("aare-rawfile-benchmark-" + std::to_string(stamp)); + std::filesystem::create_directory(m_directory); + const nlohmann::json metadata = { + {"Version", 7.2}, + {"Detector Type", "Jungfrau"}, + {"Timing Mode", "auto"}, + {"Geometry", {{"x", 1}, {"y", 1}}}, + {"Image Size in bytes", rows * cols * sizeof(uint16_t)}, + {"Pixels", {{"x", cols}, {"y", rows}}}, + {"Max Frames Per File", frame_count}, + {"Total Frames", frame_count}, + {"Frames in File", frame_count}, + {"Frame Padding", 1}, + {"Frame Discard Policy", "nodiscard"}}; + std::ofstream(master()) << metadata; + + std::vector pixels(rows * cols); + for (size_t i = 0; i < pixels.size(); ++i) { + pixels[i] = static_cast(i % 1000 + 1); + } + std::ofstream data(m_directory / "run_d0_f0_0.raw", std::ios::binary); + for (size_t frame = 0; frame < frame_count; ++frame) { + aare::DetectorHeader header{}; + header.frameNumber = frame + 1; + data.write(reinterpret_cast(&header), sizeof(header)); + data.write( + reinterpret_cast(pixels.data()), + static_cast(pixels.size() * sizeof(uint16_t))); + } + } + + void TearDown(const benchmark::State &) override { + std::error_code error; + std::filesystem::remove_all(m_directory, error); + } + + std::filesystem::path master() const { + return m_directory / "run_master_0.json"; + } + + /// Read every frame of the file into one buffer per iteration + void read_all(benchmark::State &state, + std::optional gaps) const { + aare::RawFile file(master(), "r", gaps); + std::vector buffer(file.bytes_per_frame() * frame_count); + for (auto _ : state) { + file.seek(0); + file.read_into(buffer.data(), frame_count); + benchmark::DoNotOptimize(buffer.data()); + benchmark::ClobberMemory(); + } + state.SetItemsProcessed(static_cast(state.iterations()) * + static_cast(frame_count)); + state.SetBytesProcessed(static_cast(state.iterations()) * + static_cast(buffer.size())); + } +}; + +BENCHMARK_DEFINE_F(RawFileFixture, NoGaps)(benchmark::State &state) { + read_all(state, std::nullopt); +} + +BENCHMARK_DEFINE_F(RawFileFixture, GapsFillValue)(benchmark::State &state) { + read_all(state, aare::GapPixels{}); +} + +BENCHMARK_DEFINE_F(RawFileFixture, GapsSplitCounts)(benchmark::State &state) { + aare::GapPixels gaps; + gaps.split_counts = true; + gaps.seed = 1; + read_all(state, gaps); +} + +BENCHMARK_REGISTER_F(RawFileFixture, NoGaps)->Unit(benchmark::kMillisecond); +BENCHMARK_REGISTER_F(RawFileFixture, GapsFillValue) + ->Unit(benchmark::kMillisecond); +BENCHMARK_REGISTER_F(RawFileFixture, GapsSplitCounts) + ->Unit(benchmark::kMillisecond); + +} // namespace diff --git a/docs/src/GapPixels.rst b/docs/src/GapPixels.rst new file mode 100644 index 00000000..f5169e7a --- /dev/null +++ b/docs/src/GapPixels.rst @@ -0,0 +1,122 @@ +GapPixels +============= + +Jungfrau and Eiger modules are built from 256x256 pixel chips. The pixels +at a chip edge are physically double size, so an image that places the +chips next to each other is distorted. Gap pixels restore the geometry by +inserting pixels at every chip boundary inside a module and, optionally, +gaps between modules. + +Gap positions are determined in detector coordinates before any ROI. +A ROI only receives the gaps that fall strictly inside it, so a ROI that +starts at column 256 has no gap on its left edge. + +Configuration +------------- + +``GapPixels`` holds the configuration. The chip size (256 pixels) and the +module size (1024x512 pixels) are fixed. An Eiger quad (512x512 pixels) +needs no special handling, since its only chip boundaries are at pixel 256. +``RawFile`` rejects other detector types with ``std::invalid_argument``. +Gaps of up to 1024 pixels are supported, far beyond the 2 to 10 pixels of +a real detector, and detector coordinates are limited to 2^28. + +================ ================ =========================================== +Field Default Meaning +================ ================ =========================================== +``chip_gap`` 2 Pixels inserted at each chip boundary + inside a module. Source column 255 expands + to the right and column 256 to the left. +``module_gaps`` unset ``ModuleGaps{x, y}`` inserted at module + boundaries. Unset means no pixels are + inserted between modules. + ``ModuleGaps{8, 36}`` matches the + slsDetectorPackage GUI. +``fill_value`` 0 Value written to gap pixels, converted to + the pixel type. Finite values wrap for + integer pixels; NaN, infinity and + magnitudes of 2^63 or more raise for + integer pixels. Module gaps always + receive this value. +``split_counts`` false Split the counts of the double size edge + pixels between the pixel and its gap pixel + instead of writing the fill value. +``seed`` unset Seed for the random assignment of integer + remainders in split mode. +================ ================ =========================================== + +Image sizes with the default configuration and with ``ModuleGaps{8, 36}``: + +====================================== ========== ========= =================== +Detector range Ungapped Defaults ``ModuleGaps{8,36}`` +====================================== ========== ========= =================== +One module, columns 1024 1030 1030 +Two modules side by side, columns 2048 2060 2068 +Two modules stacked, rows 1024 1028 1064 +ROI columns 100 to 400 300 302 302 +ROI columns 256 to 512 256 256 256 +====================================== ========== ========= =================== + +Split counts +------------ + +With ``split_counts`` every double size pixel at a chip boundary is halved +between itself and its gap pixel. Integer pixels put an odd remainder on +one side at random; floating point pixels are halved exactly. Columns are +split first and rows second, so a pixel at a chip corner ends up in four +parts with its total preserved. Module gaps are not split and keep the +fill value. Splitting only happens where a gap pixel is inserted, so an +edge pixel on a ROI border keeps its full value. Split mode requires a +chip gap of 2. + +Reading files +------------- + +Pass a ``GapPixels`` to :cpp:class:`aare::RawFile`. Frame sizes, +``rows()``, ``cols()`` and ``bytes_per_frame()`` then include the gaps. +Module data is scattered into the gapped frame while it is assembled, so no +extra frame sized buffer is needed. Headers are unchanged. The generic +:cpp:class:`aare::File` does not insert gap pixels; use ``RawFile`` or +apply ``insert_gap_pixels`` to the frames it returns. + +.. code-block:: cpp + + #include "aare/GapPixels.hpp" + #include "aare/RawFile.hpp" + + aare::GapPixels gaps; + gaps.module_gaps = aare::ModuleGaps{8, 36}; + aare::RawFile file("run_master_0.json", "r", gaps); + auto frame = file.read_frame(); // 1064 x 1030 for two stacked modules + +Images in memory +---------------- + +``insert_gap_pixels`` inserts gap pixels into an image in memory, such as a +calibration map or a pedestal, and returns the gapped copy. The ROI gives +the detector coordinates of the image, so an image that covers part of the +detector receives exactly the gaps that fall inside it. + +.. code-block:: cpp + + #include "aare/GapPixels.hpp" + + aare::NDArray image = ...; // 512 x 1024 + auto gapped = aare::insert_gap_pixels( + image.view(), aare::ROI{0, 1024, 0, 512}, aare::GapPixels{}); + // gapped is 514 x 1030 + +``gapped_shape`` returns the shape of the gapped image of a ROI without +building it. + +.. doxygenstruct:: aare::ModuleGaps + :members: + :undoc-members: + +.. doxygenstruct:: aare::GapPixels + :members: + :undoc-members: + +.. doxygenfunction:: aare::gapped_shape + +.. doxygenfunction:: aare::insert_gap_pixels 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/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/pyGapPixels.rst b/docs/src/python/file/pyGapPixels.rst new file mode 100644 index 00000000..58e344a3 --- /dev/null +++ b/docs/src/python/file/pyGapPixels.rst @@ -0,0 +1,94 @@ +GapPixels +=================== + +.. py:currentmodule:: aare + +Jungfrau and Eiger modules are built from 256x256 pixel chips with double +size pixels at the chip edges. Gap pixels restore the geometry by inserting +pixels at every chip boundary inside a module and, optionally, gaps between +modules. +Gap positions are determined in detector coordinates before any ROI, so a +ROI only receives the gaps that fall strictly inside it. + +The chip size (256 pixels) and the module size (1024x512 pixels) are +fixed. An Eiger quad (512x512 pixels) needs no special handling, since its +only chip boundaries are at pixel 256. + +================ ================ =========================================== +Field Default Meaning +================ ================ =========================================== +``chip_gap`` 2 Pixels inserted at each chip boundary + inside a module. +``module_gaps`` None ``ModuleGaps(x, y)`` inserted at module + boundaries. None means no pixels are + inserted between modules. + ``ModuleGaps(8, 36)`` matches the + slsDetectorPackage GUI. +``fill_value`` 0 Value written to gap pixels, converted to + the pixel dtype. Module gaps always + receive this value. +``split_counts`` False Split the counts of the double size edge + pixels between the pixel and its gap + pixel. Integer remainders go to one side + at random. Requires a chip gap of 2. +``seed`` None Seed for the random remainders. +================ ================ =========================================== + +A single Jungfrau module of ``(512, 1024)`` pixels becomes ``(514, 1030)`` +with the defaults. Two stacked modules become ``(1028, 1030)`` with the +defaults and ``(1064, 1030)`` with ``ModuleGaps(8, 36)``. + +.. code-block:: python + + from aare import GapPixels, ModuleGaps, ROI, gapped_shape, insert_gap_pixels + + # defaults: chip gaps of 2 filled with 0, (512, 1024) -> (514, 1030) + gapped = insert_gap_pixels(image) + + # module gaps as in the slsDetectorPackage GUI, counts split at chip edges + gaps = GapPixels(module_gaps=ModuleGaps(8, 36), split_counts=True) + gapped = insert_gap_pixels(image, gaps) + + # an image that covers columns 100 to 400 of the first ten rows + roi = ROI(100, 400, 0, 10) + gapped_shape(roi) # (10, 302) + part = insert_gap_pixels(image[:10, 100:400], gaps, roi=roi) + +``insert_gap_pixels`` accepts 2D ``uint8``, ``uint16``, ``uint32``, +``int32``, ``float32`` and ``float64`` arrays and returns a new array of the +same dtype. Pass ``roi`` when the image covers part of the detector. +``gapped_shape`` returns the shape of the gapped image of a ROI. + +Reading files +------------- + +``RawFile`` takes ``gap_pixels``: ``True`` for the defaults or a +``GapPixels``. Frame shapes, ``rows`` and ``cols`` then include the gaps. +Only Jungfrau and Eiger raw files are supported; other detector types raise +``ValueError``. The generic ``File`` does not insert gap pixels; use +``RawFile`` or apply ``insert_gap_pixels`` to the frames it returns. + +.. code-block:: python + + from aare import GapPixels, ModuleGaps, RawFile + + # defaults: chip gaps of 2 filled with 0 + with RawFile("run_master_0.json", gap_pixels=True) as f: + header, image = f.read_frame() + + # module gaps as in the slsDetectorPackage GUI + gaps = GapPixels(module_gaps=ModuleGaps(8, 36), split_counts=False) + with RawFile("run_master_0.json", gap_pixels=gaps) as f: + header, image = f.read_frame() + +.. autoclass:: GapPixels + :members: + :undoc-members: + +.. autoclass:: ModuleGaps + :members: + :undoc-members: + +.. autofunction:: insert_gap_pixels + +.. autofunction:: gapped_shape 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/GapPixels.hpp b/include/aare/GapPixels.hpp new file mode 100644 index 00000000..84b452f7 --- /dev/null +++ b/include/aare/GapPixels.hpp @@ -0,0 +1,379 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" +#include "aare/ROI.hpp" +#include "aare/defs.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace aare { + +/** + * @brief Pixels inserted at each boundary between modules. + */ +struct ModuleGaps { + ssize_t x{}; + ssize_t y{}; + + bool operator==(const ModuleGaps &other) const { + return x == other.x && y == other.y; + } + bool operator!=(const ModuleGaps &other) const { return !(*this == other); } +}; + +/** + * @brief Configuration for inserting gap pixels into Jungfrau and Eiger + * images. Chip and module sizes are fixed: 256x256 pixel chips in 1024x512 + * pixel modules. Gaps are inserted at the chip boundaries inside each module + * and, only when module_gaps is set, between modules. An Eiger quad is + * 512x512 pixels, so its only chip boundaries are at pixel 256 and the + * module width never comes into play. + * + * Gap positions are determined in detector coordinates before any ROI, so + * a ROI only receives the gaps that fall strictly inside it. + */ +struct GapPixels { + /// @brief Pixels inserted at each chip boundary inside a module, at + /// most 1024. Source column 255 expands to the right and column 256 to + /// the left. + ssize_t chip_gap{2}; + + /// @brief Pixels inserted at module boundaries, at most 1024 each. + /// Unset means no pixels are inserted between modules. + std::optional module_gaps{}; + + /// @brief Value written to gap pixels, converted to the pixel type. + /// Finite values wrap for integer pixel types; NaN, infinity and + /// magnitudes of 2^63 or more are rejected. Module gaps always receive + /// this value. + double fill_value{0}; + + /// @brief Split the counts of the double-size pixels at chip boundaries + /// between the pixel and its gap pixel instead of writing the fill + /// value. Integer remainders go to one side at random. Requires a chip + /// gap of 2. + bool split_counts{false}; + + /// @brief Seed for the random assignment of remainders. Unset means the + /// generator is seeded from the system. + std::optional seed{}; + + bool operator==(const GapPixels &other) const { + return chip_gap == other.chip_gap && module_gaps == other.module_gaps && + fill_value == other.fill_value && + split_counts == other.split_counts && seed == other.seed; + } + bool operator!=(const GapPixels &other) const { return !(*this == other); } +}; + +namespace detail { + +/** + * @brief Gap layout along one axis in ungapped detector coordinates. + */ +struct AxisGaps { + ssize_t chip{}; + ssize_t module{}; + ssize_t chip_gap{}; + ssize_t module_gap{}; + + /// @brief Gapped coordinate of an ungapped detector coordinate + ssize_t gapped(ssize_t coordinate) const { + const ssize_t chips_before = coordinate / chip; + const ssize_t modules_before = coordinate / module; + return coordinate + chip_gap * (chips_before - modules_before) + + module_gap * modules_before; + } + + /// @brief Gapped size of the detector range [min, max) + ssize_t extent(ssize_t min, ssize_t max) const { + return max <= min ? 0 : gapped(max - 1) - gapped(min) + 1; + } + + bool is_module_boundary(ssize_t boundary) const { + return boundary % module == 0; + } + + /// @brief Chip boundaries strictly inside [min, max). A boundary is the + /// detector coordinate of the first pixel of the chip after it. + std::vector inner_boundaries(ssize_t min, ssize_t max) const { + std::vector boundaries; + for (ssize_t boundary = (min / chip + 1) * chip; boundary < max; + boundary += chip) { + boundaries.push_back(boundary); + } + return boundaries; + } + + /// @brief Gapped column range [first, last] of the gap at a boundary, + /// relative to the gapped coordinate of `origin` + std::pair gap_range(ssize_t boundary, + ssize_t origin) const { + const ssize_t start = gapped(origin); + return {gapped(boundary - 1) + 1 - start, gapped(boundary) - 1 - start}; + } +}; + +struct GapLayout { + AxisGaps x; + AxisGaps y; +}; + +/** + * @brief Resolve the gap layout of a configuration + * @throws std::invalid_argument for an inconsistent configuration + */ +GapLayout gap_layout(const GapPixels &gaps); + +/// @brief Gapped shape {rows, cols} of a detector coordinate ROI +std::array gapped_shape(const ROI &roi, const GapLayout &layout); + +void validate_roi(const ROI &roi); +void validate_destination(const std::array &shape, const ROI &roi, + const GapLayout &layout); + +std::mt19937_64 make_generator(std::optional seed); + +/** + * @brief Convert the fill value to the pixel type + * + * Finite values wrap into integer pixel types. NaN, infinity and magnitudes + * of 2^63 or more have no integer representation and are rejected, since + * converting them would be undefined behaviour. + * @throws std::invalid_argument for a value without an integer + * representation when T is an integer type + */ +template T fill_as(double value) { + if constexpr (std::is_floating_point_v) { + return static_cast(value); + } else { + if (!(std::fabs(value) < 0x1p63)) { + throw std::invalid_argument( + LOCATION + fmt::format("Fill value {} cannot be converted to " + "an integer pixel type", + value)); + } + return static_cast(static_cast(value)); + } +} + +/** + * @brief Write `value` to every gap row and gap column inside the ROI + * @param destination gapped image of `roi` + * @param roi detector coordinates of the destination image + * @param layout resolved gap layout + * @param value value written to the gap pixels + */ +template +void fill_gap_pixels(NDView destination, const ROI &roi, + const GapLayout &layout, T value) { + validate_destination(destination.shape(), roi, layout); + const ssize_t rows = destination.shape(0); + const ssize_t cols = destination.shape(1); + + for (const auto boundary : layout.x.inner_boundaries(roi.xmin, roi.xmax)) { + const auto [first, last] = layout.x.gap_range(boundary, roi.xmin); + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = first; col <= last; ++col) { + destination(row, col) = value; + } + } + } + for (const auto boundary : layout.y.inner_boundaries(roi.ymin, roi.ymax)) { + const auto [first, last] = layout.y.gap_range(boundary, roi.ymin); + for (ssize_t row = first; row <= last; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + destination(row, col) = value; + } + } + } +} + +/** + * @brief Copy a detector coordinate rectangle into its gapped position + * @param source ungapped pixels of `source_rect` + * @param source_rect detector coordinates of `source`, inside `roi` + * @param destination gapped image of `roi` + * @param roi detector coordinates of the destination image + * @param layout resolved gap layout + * @note Gap pixels are not written + */ +template +void copy_with_gaps(NDView source, const ROI &source_rect, + NDView destination, const ROI &roi, + const GapLayout &layout) { + validate_destination(destination.shape(), roi, layout); + validate_roi(source_rect); + if (source_rect.xmin < roi.xmin || source_rect.xmax > roi.xmax || + source_rect.ymin < roi.ymin || source_rect.ymax > roi.ymax) { + throw std::invalid_argument(LOCATION + + "Source rectangle extends outside the ROI"); + } + if (source.shape(0) != source_rect.height() || + source.shape(1) != source_rect.width()) { + throw std::invalid_argument( + LOCATION + "Source shape does not match the source rectangle"); + } + if (source.strides()[1] != 1 || destination.strides()[1] != 1) { + throw std::invalid_argument( + LOCATION + "Source and destination rows must be contiguous"); + } + if (source_rect.width() == 0 || source_rect.height() == 0) { + return; + } + + // contiguous column segments [begin, end) between chip boundaries + std::vector> segments; + ssize_t begin = source_rect.xmin; + for (const auto boundary : + layout.x.inner_boundaries(source_rect.xmin, source_rect.xmax)) { + segments.emplace_back(begin, boundary); + begin = boundary; + } + segments.emplace_back(begin, source_rect.xmax); + + const ssize_t x0 = layout.x.gapped(roi.xmin); + const ssize_t y0 = layout.y.gapped(roi.ymin); + for (ssize_t y = source_rect.ymin; y < source_rect.ymax; ++y) { + const T *source_row = &source(y - source_rect.ymin, 0); + T *destination_row = &destination(layout.y.gapped(y) - y0, 0); + for (const auto &[first, last] : segments) { + std::memcpy(destination_row + layout.x.gapped(first) - x0, + source_row + first - source_rect.xmin, + static_cast(last - first) * sizeof(T)); + } + } +} + +/** + * @brief Split the counts of the double-size pixels at every chip boundary + * inside the ROI between the pixel and its gap pixel. Columns are split + * first and rows second, so corner pixels end up in four parts with their + * total preserved. Module boundaries are left untouched. + * @param destination gapped image of `roi` after copy_with_gaps + * @param roi detector coordinates of the destination image + * @param layout resolved gap layout + * @param generator source of the random remainder assignment + */ +template +void split_edge_counts(NDView destination, const ROI &roi, + const GapLayout &layout, std::mt19937_64 &generator) { + validate_destination(destination.shape(), roi, layout); + if (layout.x.chip_gap != 2 || layout.y.chip_gap != 2) { + throw std::invalid_argument( + LOCATION + "Splitting counts requires a chip gap of 2"); + } + + auto halve = [&](T &pixel, T &gap) { + if constexpr (std::is_floating_point_v) { + const T half = pixel / 2; + pixel = half; + gap = half; + } else { + T low = pixel / 2; + T high = pixel - low; + if (low != high && (generator() & 1u)) { + std::swap(low, high); + } + pixel = low; + gap = high; + } + }; + + const ssize_t x0 = layout.x.gapped(roi.xmin); + const ssize_t y0 = layout.y.gapped(roi.ymin); + const ssize_t cols = destination.shape(1); + + // columns: every pixel row of every chip boundary that is not a module + // boundary + std::vector chip_columns; + for (const auto boundary : layout.x.inner_boundaries(roi.xmin, roi.xmax)) { + if (!layout.x.is_module_boundary(boundary)) { + chip_columns.push_back(boundary); + } + } + for (ssize_t y = roi.ymin; y < roi.ymax; ++y) { + const ssize_t row = layout.y.gapped(y) - y0; + for (const auto boundary : chip_columns) { + const ssize_t left = layout.x.gapped(boundary - 1) - x0; + const ssize_t right = layout.x.gapped(boundary) - x0; + halve(destination(row, left), destination(row, left + 1)); + halve(destination(row, right), destination(row, right - 1)); + } + } + + // rows: every gapped column except module gaps + std::vector module_gap_column(static_cast(cols), false); + for (const auto boundary : layout.x.inner_boundaries(roi.xmin, roi.xmax)) { + if (layout.x.is_module_boundary(boundary)) { + const auto [first, last] = layout.x.gap_range(boundary, roi.xmin); + for (ssize_t col = first; col <= last; ++col) { + module_gap_column[static_cast(col)] = true; + } + } + } + for (const auto boundary : layout.y.inner_boundaries(roi.ymin, roi.ymax)) { + if (layout.y.is_module_boundary(boundary)) { + continue; + } + const ssize_t top = layout.y.gapped(boundary - 1) - y0; + const ssize_t bottom = layout.y.gapped(boundary) - y0; + for (ssize_t col = 0; col < cols; ++col) { + if (module_gap_column[static_cast(col)]) { + continue; + } + halve(destination(top, col), destination(top + 1, col)); + halve(destination(bottom, col), destination(bottom - 1, col)); + } + } +} + +} // namespace detail + +/** + * @brief Shape {rows, cols} of the gapped image of a detector coordinate ROI + * @param roi detector coordinates of the image + * @param gaps gap configuration + * @throws std::invalid_argument for an invalid ROI or configuration + */ +std::array gapped_shape(const ROI &roi, const GapPixels &gaps = {}); + +/** + * @brief Insert gap pixels into an image + * @param source ungapped image of `roi` + * @param roi detector coordinates of `source` + * @param gaps gap configuration + * @return gapped image + * @throws std::invalid_argument for an invalid configuration, a fill value + * that cannot be converted to the pixel type or a source shape that does + * not match the ROI + */ +template +NDArray insert_gap_pixels(NDView source, const ROI &roi, + const GapPixels &gaps = {}) { + const auto layout = detail::gap_layout(gaps); + detail::validate_roi(roi); + const T fill = detail::fill_as(gaps.fill_value); + NDArray result(detail::gapped_shape(roi, layout)); + detail::fill_gap_pixels(result.view(), roi, layout, fill); + detail::copy_with_gaps(source, roi, result.view(), roi, layout); + if (gaps.split_counts) { + auto generator = detail::make_generator(gaps.seed); + detail::split_edge_counts(result.view(), roi, layout, generator); + } + return result; +} + +} // namespace aare diff --git a/include/aare/NDView.hpp b/include/aare/NDView.hpp index 715127c9..66ea7298 100644 --- a/include/aare/NDView.hpp +++ b/include/aare/NDView.hpp @@ -221,6 +221,21 @@ class NDView : public ArrayExpr, Ndim> { return NDView(buffer_ + first * strides_[0], new_shape); } + /** + * @brief Create a subview of a single index of the first dimension. + * This is useful for accessing a single frame in a batch of frames. + * @param index The index of the subview. + * @return A new NDView that is a subview of the current view. + * @throws std::runtime_error if the index is invalid. + */ + template 1) && D == Ndim, int> = 0> + NDView sub_view(const ssize_t index) const { + if (index < 0 || index >= shape_[0]) + throw std::runtime_error(LOCATION + "Invalid sub_view index"); + auto new_shape = drop_first_dim(shape_); + return NDView(buffer_ + index * strides_[0], new_shape); + } + private: template friend class NDView; diff --git a/include/aare/RawFile.hpp b/include/aare/RawFile.hpp index 4ef92962..3314c8cf 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 { @@ -26,6 +29,10 @@ class RawFile : public FileInterface { std::vector>> m_subfiles; // [ROI][modules_per_ROI] + /// @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; + RawMasterFile m_master; size_t m_current_frame{}; @@ -34,15 +41,32 @@ class RawFile : public FileInterface { /// @brief Minimum frame count across the selected raw subfile series. size_t m_frames_in_file{}; + /// @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 +126,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 +182,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/aare/transform.py b/python/aare/transform.py index 6d14231f..7739d69e 100644 --- a/python/aare/transform.py +++ b/python/aare/transform.py @@ -135,7 +135,7 @@ class Matterhorn10Transform: elif self.dynamic_range == 8: return np.take(data.view(np.uint8), self.pixel_map).reshape(self.num_counters, Matterhorn10.nRows, Matterhorn10.nCols) else: #dynamic range 4 - return np.take(_aare.expand4to8bit(data.view(np.uint8)), self.pixel_map).reshape(self.num_counters, Matterhorn10.nRows, Matterhorn10.nCols) + return np.take(_aare.expand4to8bit(data.view(np.uint8).reshape(-1)), self.pixel_map).reshape(self.num_counters, Matterhorn10.nRows, Matterhorn10.nCols) class Mythen302Transform: """ @@ -169,7 +169,7 @@ class Mythen302Transform: image : np.ndarray uint32 array of size 64, 3 """ - res = _aare.decode_my302(data, self.offset) + res = _aare.decode_my302(data.reshape(-1), self.offset) res = res.reshape( Mythen302Transform._n_channels, Mythen302Transform._n_counters ) diff --git a/python/src/bind_GapPixels.hpp b/python/src/bind_GapPixels.hpp new file mode 100644 index 00000000..41874ca5 --- /dev/null +++ b/python/src/bind_GapPixels.hpp @@ -0,0 +1,174 @@ +// SPDX-License-Identifier: MPL-2.0 +#pragma once +#include "aare/GapPixels.hpp" +#include "aare/ROI.hpp" +#include "aare/defs.hpp" +#include "np_helper.hpp" + +#include +#include +#include +#include +#include +#include + +namespace py = pybind11; +using namespace ::aare; + +namespace { + +template +py::array insert_gap_pixels_typed(const py::array &image, const ROI &roi, + const GapPixels &gaps) { + py::array_t typed(image); + auto source = make_const_view_2d(typed); + auto *result = new NDArray(insert_gap_pixels(source, roi, gaps)); + return return_image_data(result); +} + +} // namespace + +void define_gap_pixels_bindings(py::module &m) { + py::class_(m, "ModuleGaps", R"( + Pixels inserted at each boundary between modules. + )") + .def(py::init<>()) + .def(py::init(), py::arg("x"), py::arg("y")) + .def_readwrite("x", &ModuleGaps::x) + .def_readwrite("y", &ModuleGaps::y) + .def(py::self == py::self) + .def(py::self != py::self) + .def("__repr__", [](const ModuleGaps &self) { + return fmt::format("ModuleGaps(x={}, y={})", self.x, self.y); + }); + + py::class_(m, "GapPixels", R"( + Configuration for inserting gap pixels into Jungfrau and Eiger + images. Chip and module sizes are fixed: 256x256 pixel chips in + 1024x512 pixel modules. + + Parameters + ---------- + chip_gap : int + Pixels inserted at each chip boundary inside a module, at most + 1024. Source column 255 expands to the right and column 256 to + the left. + module_gaps : ModuleGaps, optional + Pixels inserted at module boundaries, at most 1024 each. None + means no pixels are inserted between modules. + fill_value : float + Value written to gap pixels, converted to the pixel type. Finite + values wrap for integer images; NaN, infinity and magnitudes of + 2**63 or more raise ValueError for integer images. Module gaps + always receive this value. + split_counts : bool + Split the counts of the double-size pixels at chip boundaries + between the pixel and its gap pixel. Integer remainders go to + one side at random. Requires a chip gap of 2. + seed : int, optional + Seed for the random assignment of remainders. + )") + .def( + py::init([](ssize_t chip_gap, std::optional module_gaps, + double fill_value, bool split_counts, + std::optional seed) { + GapPixels gaps; + gaps.chip_gap = chip_gap; + gaps.module_gaps = module_gaps; + gaps.fill_value = fill_value; + gaps.split_counts = split_counts; + gaps.seed = seed; + return gaps; + }), + py::arg("chip_gap") = 2, py::arg("module_gaps") = py::none(), + py::arg("fill_value") = 0.0, py::arg("split_counts") = false, + py::arg("seed") = py::none()) + .def_readwrite("chip_gap", &GapPixels::chip_gap) + .def_readwrite("module_gaps", &GapPixels::module_gaps) + .def_readwrite("fill_value", &GapPixels::fill_value) + .def_readwrite("split_counts", &GapPixels::split_counts) + .def_readwrite("seed", &GapPixels::seed) + .def(py::self == py::self) + .def(py::self != py::self) + .def("__repr__", [](const GapPixels &self) { + return fmt::format( + "GapPixels(chip_gap={}, module_gaps={}, fill_value={}, " + "split_counts={}, seed={})", + self.chip_gap, + self.module_gaps + ? fmt::format("ModuleGaps(x={}, y={})", self.module_gaps->x, + self.module_gaps->y) + : "None", + self.fill_value, self.split_counts ? "True" : "False", + self.seed ? fmt::format("{}", *self.seed) : "None"); + }); + + m.def( + "gapped_shape", + [](const ROI &roi, const GapPixels &gaps) { + const auto shape = gapped_shape(roi, gaps); + return py::make_tuple(shape[0], shape[1]); + }, + py::arg("roi"), py::arg("gaps") = GapPixels{}, R"( + Shape (rows, cols) of the gapped image of a detector coordinate ROI. + + Parameters + ---------- + roi : ROI + Detector coordinates of the image. + gaps : GapPixels + Gap configuration. + )"); + + m.def( + "insert_gap_pixels", + [](const py::array &image, const GapPixels &gaps, + std::optional roi) -> py::array { + if (image.ndim() != 2) { + throw py::value_error( + fmt::format("Expected a 2D image, got {}D", image.ndim())); + } + const ROI rect = + roi.value_or(ROI{0, image.shape(1), 0, image.shape(0)}); + if (py::isinstance>(image)) { + return insert_gap_pixels_typed(image, rect, gaps); + } + if (py::isinstance>(image)) { + return insert_gap_pixels_typed(image, rect, gaps); + } + if (py::isinstance>(image)) { + return insert_gap_pixels_typed(image, rect, gaps); + } + if (py::isinstance>(image)) { + return insert_gap_pixels_typed(image, rect, gaps); + } + if (py::isinstance>(image)) { + return insert_gap_pixels_typed(image, rect, gaps); + } + if (py::isinstance>(image)) { + return insert_gap_pixels_typed(image, rect, gaps); + } + throw py::type_error( + "Unsupported dtype, expected uint8, uint16, uint32, int32, " + "float32 or float64"); + }, + py::arg("image"), py::arg("gaps") = GapPixels{}, + py::arg("roi") = py::none(), R"( + Insert gap pixels into a 2D image. + + Parameters + ---------- + image : numpy.ndarray + Ungapped image of ``roi``. + gaps : GapPixels + Gap configuration. + roi : ROI, optional + Detector coordinates of the image. Defaults to the full image + starting at pixel (0, 0). + + Returns + ------- + numpy.ndarray + New array with the same dtype as ``image``. + )"); +} 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/module.cpp b/python/src/module.cpp index 2eb89d0e..66996151 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -51,6 +51,7 @@ PYBIND11_MODULE(_aare, m) { bind_eta_3x3(m); define_defs_bindings(m); + define_gap_pixels_bindings(m); bind_eta_types(m); } diff --git a/python/tests/test_GapPixels.py b/python/tests/test_GapPixels.py new file mode 100644 index 00000000..9fde15eb --- /dev/null +++ b/python/tests/test_GapPixels.py @@ -0,0 +1,130 @@ +# SPDX-License-Identifier: MPL-2.0 +import numpy as np +import pytest + +from aare import GapPixels, ModuleGaps, ROI, gapped_shape, insert_gap_pixels + + +def reference(image, boundaries_x, boundaries_y, fill=0): + """Insert gap columns and rows with numpy: boundaries map a source index to + the number of pixels inserted before it.""" + cols = [index for index, count in boundaries_x for _ in range(count)] + rows = [index for index, count in boundaries_y for _ in range(count)] + result = np.insert(image, cols, fill, axis=1) + return np.insert(result, rows, fill, axis=0) + + +def test_configuration_defaults_and_equality(): + gaps = GapPixels() + assert gaps.chip_gap == 2 + assert gaps.module_gaps is None + assert gaps.fill_value == 0 + assert gaps.split_counts is False + assert gaps.seed is None + assert gaps == GapPixels() + assert gaps != GapPixels(fill_value=1) + + gaps = GapPixels(chip_gap=4, module_gaps=ModuleGaps(8, 36), fill_value=-1, split_counts=False, seed=5) + assert gaps.module_gaps == ModuleGaps(x=8, y=36) + assert gaps.seed == 5 + assert "ModuleGaps(x=8, y=36)" in repr(gaps) + gaps.module_gaps = None + assert gaps.module_gaps is None + + +def test_gapped_shape(): + assert gapped_shape(ROI(0, 1024, 0, 512)) == (514, 1030) + assert gapped_shape(ROI(0, 1024, 0, 1024)) == (1028, 1030) + assert gapped_shape(ROI(0, 1024, 0, 1024), GapPixels(module_gaps=ModuleGaps(8, 36))) == (1064, 1030) + assert gapped_shape(ROI(0, 512, 0, 512)) == (514, 514) # Eiger quad + assert gapped_shape(ROI(256, 512, 0, 256)) == (256, 256) + assert gapped_shape(ROI(255, 257, 0, 1)) == (1, 4) + + +@pytest.mark.parametrize("dtype", [np.uint8, np.uint16, np.uint32, np.int32, np.float32, np.float64]) +def test_insert_gap_pixels_matches_numpy(dtype): + rng = np.random.default_rng(1) + image = rng.integers(0, 200, size=(512, 1024)).astype(dtype) + result = insert_gap_pixels(image, GapPixels(fill_value=3)) + expected = reference(image, [(256, 2), (512, 2), (768, 2)], [(256, 2)], fill=3) + assert result.dtype == image.dtype + assert result.shape == (514, 1030) + assert (result == expected).all() + + +def test_insert_gap_pixels_without_module_gaps(): + # two stacked modules: chip gaps at rows 256 and 768, none at the module + # boundary at row 512 + image = np.arange(1024 * 1024, dtype=np.uint32).reshape(1024, 1024) + result = insert_gap_pixels(image) + expected = reference(image, [(256, 2), (512, 2), (768, 2)], [(256, 2), (768, 2)]) + assert result.shape == (1028, 1030) + assert (result == expected).all() + + +def test_insert_gap_pixels_with_module_gaps(): + image = np.arange(1024 * 1024, dtype=np.uint32).reshape(1024, 1024) + gaps = GapPixels(module_gaps=ModuleGaps(8, 36), fill_value=9) + result = insert_gap_pixels(image, gaps) + expected = reference(image, [(256, 2), (512, 2), (768, 2)], [(256, 2), (512, 36), (768, 2)], fill=9) + assert result.shape == (1064, 1030) + assert (result == expected).all() + + +def test_insert_gap_pixels_in_roi(): + image = np.array([[11, 22]], dtype=np.uint16) + result = insert_gap_pixels(image, roi=ROI(255, 257, 10, 11)) + assert result.tolist() == [[11, 0, 0, 22]] + + # a gap on the ROI border is not inserted + image = np.arange(8, dtype=np.uint16).reshape(2, 4) + result = insert_gap_pixels(image, roi=ROI(256, 260, 0, 2)) + assert (result == image).all() + + +def test_insert_gap_pixels_split_counts(): + rng = np.random.default_rng(2) + image = rng.integers(0, 1000, size=(512, 1024)).astype(np.uint32) + gaps = GapPixels(split_counts=True, seed=3) + first = insert_gap_pixels(image, gaps) + second = insert_gap_pixels(image, gaps) + assert first.sum() == image.sum() + assert (first == second).all() + # edge pixel and its gap pixel hold floor and ceil of the half; row 255 + # is excluded because it is split again vertically + pair = np.stack([first[:255, 255], first[:255, 256]]) + assert (pair.sum(axis=0) == image[:255, 255]).all() + assert (pair.max(axis=0) == (image[:255, 255] + 1) // 2).all() + + halves = insert_gap_pixels(np.full((2, 2), 3.0), gaps, roi=ROI(255, 257, 255, 257)) + assert (halves == 0.75).all() + + +def test_insert_gap_pixels_rejects_bad_input(): + image = np.zeros((400, 400), dtype=np.uint16) + with pytest.raises(ValueError): + insert_gap_pixels(image, roi=ROI(0, 401, 0, 400)) + with pytest.raises(ValueError): + insert_gap_pixels(np.zeros((2, 2, 2), dtype=np.uint16)) + with pytest.raises(TypeError): + insert_gap_pixels(np.zeros((2, 2), dtype=np.int64)) + with pytest.raises(ValueError): + insert_gap_pixels(image, GapPixels(split_counts=True, chip_gap=4)) + with pytest.raises(ValueError): + insert_gap_pixels(image, GapPixels(fill_value=float("nan"))) + with pytest.raises(ValueError): + insert_gap_pixels(image, GapPixels(fill_value=1e20)) + with pytest.raises(ValueError): + insert_gap_pixels(image, GapPixels(chip_gap=6148914691236517205)) + with pytest.raises(ValueError): + insert_gap_pixels(image, GapPixels(module_gaps=ModuleGaps(8, 2000))) + with pytest.raises(ValueError): + insert_gap_pixels(image, roi=ROI(2**40, 2**40 + 400, 0, 400)) + + +def test_insert_gap_pixels_nan_fill_for_float_images(): + image = np.ones((2, 2), dtype=np.float32) + result = insert_gap_pixels(image, GapPixels(fill_value=float("nan")), roi=ROI(255, 257, 0, 2)) + assert result.shape == (2, 4) + assert np.isnan(result[:, 1:3]).all() + assert (result[:, [0, 3]] == 1).all() diff --git a/python/tests/test_RawFile.py b/python/tests/test_RawFile.py index 76988db9..d104a3c3 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,120 @@ 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 + # one chip boundary per module row, three per module column and no + # gap between modules + assert (reader.rows(), reader.cols()) == ((512 + 2) * layout_rows, (1024 + 6) * layout_cols) + header, image = reader.read_frame() + assert header.size == layout_rows * layout_cols + assert (image == insert_gap_pixels(expected, 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() + assert (reader.rows(), reader.cols()) == insert_gap_pixels(expected).shape + _, image = reader.read_frame() + assert (image == insert_gap_pixels(expected)).all() + + module_gaps = GapPixels(module_gaps=ModuleGaps(8, 36)) + with RawFile(master_path, gap_pixels=module_gaps) as reader: + assert (reader.rows(), reader.cols()) == insert_gap_pixels(expected, module_gaps).shape + _, image = reader.read_frame() + assert (image == insert_gap_pixels(expected, module_gaps)).all() + + +def test_raw_file_gap_pixels_split_counts(tmp_path): + # one Eiger module: four 256x512 ports in a 2x2 layout. Splitting is + # rejected on raw Jungfrau data. The plain read is the reference since + # Eiger ports are row flipped on reading. + master_path, _ = write_module_layout(tmp_path, 2, 2, module_rows=256, module_cols=512, + detector_type="Eiger") + gaps = GapPixels(split_counts=True, seed=11) + with RawFile(master_path) as reader: + _, plain = reader.read_frame() + with RawFile(master_path, gap_pixels=gaps) as reader: + _, image = reader.read_frame() + assert (image == insert_gap_pixels(plain, gaps)).all() + assert image.sum() == plain.sum() + + +def test_raw_file_gap_pixels_split_counts_rejected_for_jungfrau(tmp_path): + master_path, _ = write_module_layout(tmp_path, 1, 1, module_rows=512, module_cols=1024) + RawFile(master_path, gap_pixels=GapPixels()) + with pytest.raises(ValueError): + RawFile(master_path, gap_pixels=GapPixels(split_counts=True)) + + +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) + + +@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, 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, 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/python/tests/test_Transformation.py b/python/tests/test_Transformation.py index 0d91c10a..dad7d409 100644 --- a/python/tests/test_Transformation.py +++ b/python/tests/test_Transformation.py @@ -56,3 +56,34 @@ def test_matterhorn10_16bit_4counters(test_data_path): expected_data = np.tile(np.arange(255, -1, -1,dtype=np.uint16), (4, 256, 1)) # TODO: endianess issue ? assert np.all(frames == expected_data) + + +def test_matterhorn10_4bit_accepts_2d_frame(): + """CtbRawFile.read_frame returns a (1, n_bytes) array; the 4 bit path must accept it""" + n_bytes = (transform.Matterhorn10.nRows * transform.Matterhorn10.nCols * 4) // 8 + data = (np.arange(n_bytes) % 256).astype(np.uint8).reshape(1, n_bytes) + t = transform.Matterhorn10Transform(dynamic_range=4, num_counters=1) + + frames = t(data) + + assert frames.shape == (1, 256, 256) + assert frames.dtype == np.uint8 + + flat = data.reshape(-1) + nibbles = np.stack([flat & 0x0F, flat >> 4], axis=-1).reshape(-1) + expected = np.take(nibbles, t.pixel_map).reshape(1, 256, 256) + assert np.array_equal(frames, expected) + + +def test_mythen302_accepts_2d_frame(): + """CtbRawFile.read_frame returns a (1, n_bytes) array; Mythen302Transform must accept it""" + n_bytes = 64 * 3 * 3 + 2 # two extra bytes for the 4 bit offset + rng = np.random.default_rng(0) + data = rng.integers(0, 256, size=(1, n_bytes), dtype=np.uint8) + t = transform.Mythen302Transform(offset=4) + + res = t(data) + + assert res.shape == (64, 3) + assert res.dtype == np.uint32 + assert np.array_equal(res, t(data.reshape(-1))) diff --git a/src/GapPixels.cpp b/src/GapPixels.cpp new file mode 100644 index 00000000..109b6981 --- /dev/null +++ b/src/GapPixels.cpp @@ -0,0 +1,104 @@ +// SPDX-License-Identifier: MPL-2.0 +#include "aare/GapPixels.hpp" + +#include + +#include + +namespace aare { +namespace detail { + +namespace { + +constexpr ssize_t chip = 256; +constexpr ssize_t module_x = 1024; +constexpr ssize_t module_y = 512; + +// Gap sizes and detector coordinates are bounded so that gapped coordinates +// and gapped image sizes cannot overflow ssize_t. A gapped coordinate is at +// most (1 + 2 * max_gap / chip) times the detector coordinate, and the image +// size is the product of two such extents. +constexpr ssize_t max_gap = 1024; +constexpr ssize_t max_coordinate = ssize_t{1} << 28; +constexpr ssize_t max_extent = max_coordinate * (1 + 2 * max_gap / chip); +static_assert(max_extent <= std::numeric_limits::max() / max_extent, + "gapped image sizes must fit in ssize_t"); + +} // namespace + +GapLayout gap_layout(const GapPixels &gaps) { + if (gaps.chip_gap < 0) { + throw std::invalid_argument(LOCATION + "Chip gap must not be negative"); + } + const ModuleGaps module_gaps = gaps.module_gaps.value_or(ModuleGaps{}); + if (module_gaps.x < 0 || module_gaps.y < 0) { + throw std::invalid_argument(LOCATION + + "Module gaps must not be negative"); + } + if (gaps.chip_gap > max_gap) { + throw std::invalid_argument( + LOCATION + + fmt::format("Chip gap must not exceed {} pixels", max_gap)); + } + if (module_gaps.x > max_gap || module_gaps.y > max_gap) { + throw std::invalid_argument( + LOCATION + + fmt::format("Module gaps must not exceed {} pixels", max_gap)); + } + if (gaps.split_counts && gaps.chip_gap != 2) { + throw std::invalid_argument( + LOCATION + "Splitting counts requires a chip gap of 2"); + } + return {{chip, module_x, gaps.chip_gap, module_gaps.x}, + {chip, module_y, gaps.chip_gap, module_gaps.y}}; +} + +std::array gapped_shape(const ROI &roi, const GapLayout &layout) { + return {layout.y.extent(roi.ymin, roi.ymax), + layout.x.extent(roi.xmin, roi.xmax)}; +} + +void validate_roi(const ROI &roi) { + if (roi.xmin < 0 || roi.ymin < 0 || roi.xmax < roi.xmin || + roi.ymax < roi.ymin) { + throw std::invalid_argument( + LOCATION + fmt::format("Invalid ROI x [{}, {}) y [{}, {})", + roi.xmin, roi.xmax, roi.ymin, roi.ymax)); + } + if (roi.xmax > max_coordinate || roi.ymax > max_coordinate) { + throw std::invalid_argument( + LOCATION + + fmt::format("ROI coordinates must not exceed {}", max_coordinate)); + } +} + +void validate_destination(const std::array &shape, const ROI &roi, + const GapLayout &layout) { + validate_roi(roi); + const auto expected = gapped_shape(roi, layout); + if (shape != expected) { + throw std::invalid_argument( + LOCATION + + fmt::format("Destination shape ({}, {}) does not match the gapped " + "ROI shape ({}, {})", + shape[0], shape[1], expected[0], expected[1])); + } +} + +std::mt19937_64 make_generator(std::optional seed) { + if (seed) { + return std::mt19937_64(*seed); + } + std::random_device device; + return std::mt19937_64((static_cast(device()) << 32) ^ device()); +} + +} // namespace detail + +std::array gapped_shape(const ROI &roi, const GapPixels &gaps) { + const auto layout = detail::gap_layout(gaps); + detail::validate_roi(roi); + return detail::gapped_shape(roi, layout); +} + +} // namespace aare diff --git a/src/GapPixels.test.cpp b/src/GapPixels.test.cpp new file mode 100644 index 00000000..ae792319 --- /dev/null +++ b/src/GapPixels.test.cpp @@ -0,0 +1,588 @@ +// SPDX-License-Identifier: MPL-2.0 +#include "aare/GapPixels.hpp" + +#include +#include +#include + +#include +#include +#include + +using namespace aare; + +namespace { + +// chip 2, module 4, chip gap 1, module gap 3: boundaries at 2 (chip), 4 +// (module), 6 (chip) for an 8 pixel axis +const detail::AxisGaps tiny_axis{2, 4, 1, 3}; + +template NDArray pattern(ssize_t rows, ssize_t cols) { + NDArray image({rows, cols}); + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + image(row, col) = static_cast((row * cols + col) % 60000 + 1); + } + } + return image; +} + +template double total(NDView image) { + double sum = 0; + for (ssize_t row = 0; row < image.shape(0); ++row) { + for (ssize_t col = 0; col < image.shape(1); ++col) { + sum += static_cast(image(row, col)); + } + } + return sum; +} + +} // namespace + +TEST_CASE("Axis gap map inserts chip and module gaps", "[GapPixels]") { + const detail::AxisGaps &axis = tiny_axis; + CHECK(axis.gapped(0) == 0); + CHECK(axis.gapped(1) == 1); + CHECK(axis.gapped(2) == 3); + CHECK(axis.gapped(3) == 4); + CHECK(axis.gapped(4) == 8); + CHECK(axis.gapped(5) == 9); + CHECK(axis.gapped(6) == 11); + CHECK(axis.gapped(7) == 12); + + CHECK(axis.extent(0, 8) == 13); + CHECK(axis.extent(0, 2) == 2); + CHECK(axis.extent(2, 4) == 2); + CHECK(axis.extent(1, 3) == 3); + CHECK(axis.extent(3, 5) == 5); + CHECK(axis.extent(4, 4) == 0); + CHECK(axis.extent(5, 4) == 0); + + CHECK(axis.inner_boundaries(0, 8) == std::vector{2, 4, 6}); + CHECK(axis.inner_boundaries(1, 3) == std::vector{2}); + CHECK(axis.inner_boundaries(2, 4).empty()); + CHECK(axis.inner_boundaries(2, 5) == std::vector{4}); + CHECK(axis.inner_boundaries(0, 0).empty()); + + CHECK(axis.is_module_boundary(4)); + CHECK_FALSE(axis.is_module_boundary(2)); + CHECK(axis.gap_range(2, 0) == std::pair{2, 2}); + CHECK(axis.gap_range(4, 0) == std::pair{5, 7}); + CHECK(axis.gap_range(4, 3) == std::pair{1, 3}); +} + +TEST_CASE("Gap layout follows the configuration", "[GapPixels]") { + SECTION("Jungfrau defaults") { + const auto layout = detail::gap_layout(GapPixels{}); + CHECK(layout.x.chip == 256); + CHECK(layout.x.module == 1024); + CHECK(layout.x.chip_gap == 2); + CHECK(layout.x.module_gap == 0); + CHECK(layout.y.chip == 256); + CHECK(layout.y.module == 512); + CHECK(layout.y.chip_gap == 2); + CHECK(layout.y.module_gap == 0); + } + SECTION("explicit gaps") { + GapPixels gaps; + gaps.chip_gap = 4; + gaps.module_gaps = ModuleGaps{8, 36}; + const auto layout = detail::gap_layout(gaps); + CHECK(layout.x.chip_gap == 4); + CHECK(layout.y.chip_gap == 4); + CHECK(layout.x.module_gap == 8); + CHECK(layout.y.module_gap == 36); + } + SECTION("invalid configuration throws") { + GapPixels gaps; + gaps.chip_gap = -1; + CHECK_THROWS_AS(detail::gap_layout(gaps), std::invalid_argument); + gaps = GapPixels{}; + gaps.module_gaps = ModuleGaps{-8, 36}; + CHECK_THROWS_AS(detail::gap_layout(gaps), std::invalid_argument); + gaps = GapPixels{}; + gaps.split_counts = true; + gaps.chip_gap = 4; + CHECK_THROWS_AS(detail::gap_layout(gaps), std::invalid_argument); + } + SECTION("gaps are limited to 1024 pixels") { + GapPixels gaps; + gaps.chip_gap = 1024; + gaps.module_gaps = ModuleGaps{1024, 1024}; + CHECK_NOTHROW(detail::gap_layout(gaps)); + gaps.chip_gap = 1025; + CHECK_THROWS_AS(detail::gap_layout(gaps), std::invalid_argument); + gaps.chip_gap = 2; + gaps.module_gaps = ModuleGaps{8, 1025}; + CHECK_THROWS_AS(detail::gap_layout(gaps), std::invalid_argument); + } +} + +TEST_CASE("Gapped shape of detector ROIs", "[GapPixels]") { + GapPixels module_8_36; + module_8_36.module_gaps = ModuleGaps{8, 36}; + using Shape = std::array; + + CHECK(gapped_shape(ROI{0, 1024, 0, 512}) == Shape{514, 1030}); + CHECK(gapped_shape(ROI{0, 1024, 0, 512}, module_8_36) == Shape{514, 1030}); + CHECK(gapped_shape(ROI{0, 2048, 0, 512}) == Shape{514, 2060}); + CHECK(gapped_shape(ROI{0, 2048, 0, 512}, module_8_36) == Shape{514, 2068}); + CHECK(gapped_shape(ROI{0, 1024, 0, 1024}) == Shape{1028, 1030}); + CHECK(gapped_shape(ROI{0, 1024, 0, 1024}, module_8_36) == + Shape{1064, 1030}); + // an Eiger quad + CHECK(gapped_shape(ROI{0, 512, 0, 512}) == Shape{514, 514}); + + SECTION("ROI borders exclude gaps") { + CHECK(gapped_shape(ROI{100, 400, 0, 10}) == Shape{10, 302}); + CHECK(gapped_shape(ROI{256, 512, 0, 1}) == Shape{1, 256}); + CHECK(gapped_shape(ROI{255, 257, 0, 1}) == Shape{1, 4}); + CHECK(gapped_shape(ROI{200, 301, 300, 601}) == Shape{301, 103}); + CHECK(gapped_shape(ROI{200, 301, 300, 601}, module_8_36) == + Shape{337, 103}); + CHECK(gapped_shape(ROI{512, 768, 0, 256}) == Shape{256, 256}); + } + SECTION("invalid ROI throws") { + CHECK_THROWS_AS(gapped_shape(ROI{-1, 10, 0, 10}), + std::invalid_argument); + CHECK_THROWS_AS(gapped_shape(ROI{10, 5, 0, 10}), std::invalid_argument); + } + SECTION("coordinates are limited to 2^28") { + const ssize_t limit = ssize_t{1} << 28; + CHECK(gapped_shape(ROI{limit - 1, limit, 0, 1}) == Shape{1, 1}); + CHECK_THROWS_AS(gapped_shape(ROI{limit, limit + 1, 0, 1}), + std::invalid_argument); + CHECK_THROWS_AS(gapped_shape(ROI{0, 1, limit, limit + 1}), + std::invalid_argument); + } +} + +TEST_CASE("Copy and fill with a tiny layout", "[GapPixels]") { + const detail::GapLayout layout{tiny_axis, tiny_axis}; + // 4 rows (one chip boundary at 2) by 8 columns (chip, module, chip) + const ROI roi{0, 8, 0, 4}; + auto source = pattern(4, 8); + const auto shape = detail::gapped_shape(roi, layout); + REQUIRE(shape == std::array{5, 13}); + + NDArray destination(shape, uint16_t{7}); + detail::fill_gap_pixels(destination.view(), roi, layout, uint16_t{99}); + + SECTION("whole ROI") { + detail::copy_with_gaps(NDView(source.view()), roi, + destination.view(), roi, layout); + for (ssize_t row = 0; row < shape[0]; ++row) { + for (ssize_t col = 0; col < shape[1]; ++col) { + CAPTURE(row, col); + const bool gap_row = row == 2; + const bool gap_col = + col == 2 || (col >= 5 && col <= 7) || col == 10; + if (gap_row || gap_col) { + CHECK(destination(row, col) == 99); + } + } + } + for (ssize_t row = 0; row < 4; ++row) { + for (ssize_t col = 0; col < 8; ++col) { + CAPTURE(row, col); + CHECK(destination(layout.y.gapped(row), layout.x.gapped(col)) == + source(row, col)); + } + } + } + SECTION("two rectangles cover the ROI") { + const ROI left{0, 3, 0, 4}; + const ROI right{3, 8, 0, 4}; + NDArray left_pixels({4, 3}); + NDArray right_pixels({4, 5}); + for (ssize_t row = 0; row < 4; ++row) { + for (ssize_t col = 0; col < 8; ++col) { + if (col < 3) { + left_pixels(row, col) = source(row, col); + } else { + right_pixels(row, col - 3) = source(row, col); + } + } + } + detail::copy_with_gaps(NDView(left_pixels.view()), + left, destination.view(), roi, layout); + detail::copy_with_gaps(NDView(right_pixels.view()), + right, destination.view(), roi, layout); + for (ssize_t row = 0; row < 4; ++row) { + for (ssize_t col = 0; col < 8; ++col) { + CAPTURE(row, col); + CHECK(destination(layout.y.gapped(row), layout.x.gapped(col)) == + source(row, col)); + } + } + CHECK(destination(0, 2) == 99); + CHECK(destination(2, 0) == 99); + } + SECTION("mismatched shapes throw") { + NDArray wrong({5, 12}); + CHECK_THROWS_AS( + detail::fill_gap_pixels(wrong.view(), roi, layout, uint16_t{0}), + std::invalid_argument); + CHECK_THROWS_AS( + detail::copy_with_gaps(NDView(source.view()), + roi, wrong.view(), roi, layout), + std::invalid_argument); + CHECK_THROWS_AS(detail::copy_with_gaps( + NDView(source.view()), + ROI{0, 9, 0, 4}, destination.view(), roi, layout), + std::invalid_argument); + CHECK_THROWS_AS(detail::copy_with_gaps( + NDView(source.view()), + ROI{0, 7, 0, 4}, destination.view(), roi, layout), + std::invalid_argument); + } +} + +TEST_CASE("Gaps and coordinates that overflow the gapped shape are rejected", + "[GapPixels]") { + // chip_gap * 3 wrapped around in 64 bit arithmetic and reported a 1x1023 + // image for this source, so gap filling wrote past the allocation + auto source = pattern(1, 1024); + GapPixels gaps; + gaps.chip_gap = 6148914691236517205; + const ROI roi{0, 1024, 0, 1}; + CHECK_THROWS_AS(gapped_shape(roi, gaps), std::invalid_argument); + CHECK_THROWS_AS( + insert_gap_pixels(NDView(source.view()), roi, gaps), + std::invalid_argument); + + const ssize_t far = ssize_t{1} << 62; + const ROI far_roi{far, far + 1024, 0, 1}; + CHECK_THROWS_AS(gapped_shape(far_roi), std::invalid_argument); + CHECK_THROWS_AS( + insert_gap_pixels(NDView(source.view()), far_roi), + std::invalid_argument); +} + +TEST_CASE("insert_gap_pixels on a Jungfrau module", "[GapPixels]") { + const ROI roi{0, 1024, 0, 512}; + auto source = pattern(512, 1024); + GapPixels gaps; + gaps.fill_value = 5; + auto result = + insert_gap_pixels(NDView(source.view()), roi, gaps); + REQUIRE(result.shape() == std::array{514, 1030}); + + const auto layout = detail::gap_layout(gaps); + for (ssize_t row = 0; row < 512; ++row) { + for (ssize_t col = 0; col < 1024; ++col) { + if (result(layout.y.gapped(row), layout.x.gapped(col)) != + source(row, col)) { + CAPTURE(row, col); + REQUIRE(result(layout.y.gapped(row), layout.x.gapped(col)) == + source(row, col)); + } + } + } + for (ssize_t col : {256, 257, 514, 515, 772, 773}) { + for (ssize_t row = 0; row < 514; ++row) { + if (result(row, col) != 5) { + CAPTURE(row, col); + REQUIRE(result(row, col) == 5); + } + } + } + for (ssize_t row : {256, 257}) { + for (ssize_t col = 0; col < 1030; ++col) { + if (result(row, col) != 5) { + CAPTURE(row, col); + REQUIRE(result(row, col) == 5); + } + } + } + CHECK(total(NDView(result.view())) == + total(NDView(source.view())) + + 5.0 * (514 * 1030 - 512 * 1024)); +} + +TEST_CASE("insert_gap_pixels on ROIs", "[GapPixels]") { + SECTION("gap inside the ROI") { + const ROI roi{255, 257, 10, 11}; + NDArray source({1, 2}); + source(0, 0) = 11; + source(0, 1) = 22; + auto result = + insert_gap_pixels(NDView(source.view()), roi); + REQUIRE(result.shape() == std::array{1, 4}); + CHECK(result(0, 0) == 11); + CHECK(result(0, 1) == 0); + CHECK(result(0, 2) == 0); + CHECK(result(0, 3) == 22); + } + SECTION("gap on the ROI border is not inserted") { + const ROI roi{256, 260, 0, 2}; + auto source = pattern(2, 4); + auto result = + insert_gap_pixels(NDView(source.view()), roi); + REQUIRE(result.shape() == std::array{2, 4}); + for (ssize_t row = 0; row < 2; ++row) { + for (ssize_t col = 0; col < 4; ++col) { + CHECK(result(row, col) == source(row, col)); + } + } + } + SECTION("module boundary inside the ROI") { + GapPixels gaps; + gaps.module_gaps = ModuleGaps{8, 36}; + gaps.fill_value = -1; // wraps for unsigned types + const ROI roi{200, 301, 300, 601}; + auto source = pattern(301, 101); + auto result = insert_gap_pixels( + NDView(source.view()), roi, gaps); + REQUIRE(result.shape() == std::array{337, 103}); + // chip gap at column 256 -> gapped columns 56, 57 + CHECK(result(0, 55) == source(0, 55)); + CHECK(result(0, 56) == 65535); + CHECK(result(0, 57) == 65535); + CHECK(result(0, 58) == source(0, 56)); + // module gap at row 512 -> gapped rows 212 .. 247 + CHECK(result(211, 0) == source(211, 0)); + CHECK(result(212, 0) == 65535); + CHECK(result(247, 0) == 65535); + CHECK(result(248, 0) == source(212, 0)); + } + SECTION("source shape must match the ROI") { + auto source = pattern(2, 4); + CHECK_THROWS_AS( + insert_gap_pixels(NDView(source.view()), + ROI{0, 5, 0, 2}), + std::invalid_argument); + } + SECTION("floating point images") { + const ROI roi{255, 257, 0, 1}; + NDArray source({1, 2}); + source(0, 0) = 1.5; + source(0, 1) = 2.5; + GapPixels gaps; + gaps.fill_value = 0.25; + auto result = insert_gap_pixels(NDView(source.view()), + roi, gaps); + CHECK(result(0, 0) == 1.5); + CHECK(result(0, 1) == 0.25); + CHECK(result(0, 2) == 0.25); + CHECK(result(0, 3) == 2.5); + } + SECTION("fill values without an integer representation throw") { + const ROI roi{255, 257, 0, 1}; + auto source = pattern(1, 2); + GapPixels gaps; + for (const double value : + {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity(), + -std::numeric_limits::infinity(), 1e20, -1e20, 0x1p63}) { + gaps.fill_value = value; + CHECK_THROWS_AS( + insert_gap_pixels(NDView(source.view()), roi, + gaps), + std::invalid_argument); + } + + // the same values are valid for floating point images + gaps.fill_value = std::numeric_limits::quiet_NaN(); + auto floats = pattern(1, 2); + auto result = + insert_gap_pixels(NDView(floats.view()), roi, gaps); + CHECK(std::isnan(result(0, 1))); + CHECK(std::isnan(result(0, 2))); + } +} + +TEST_CASE("Module boundaries receive no gap unless module gaps are set", + "[GapPixels]") { + // rows 510..513 of two stacked modules hold the module boundary at 512 + const ROI roi{0, 2, 510, 514}; + auto source = pattern(4, 2); + GapPixels gaps; + gaps.fill_value = 9; + SECTION("fill value") {} + SECTION("split counts") { + gaps.split_counts = true; + gaps.seed = 1; + } + auto result = + insert_gap_pixels(NDView(source.view()), roi, gaps); + REQUIRE(result.shape() == std::array{4, 2}); + CHECK(result == source); + + gaps.module_gaps = ModuleGaps{0, 3}; + auto with_gap = + insert_gap_pixels(NDView(source.view()), roi, gaps); + REQUIRE(with_gap.shape() == std::array{7, 2}); + for (ssize_t row = 2; row < 5; ++row) { + CHECK(with_gap(row, 0) == 9); + CHECK(with_gap(row, 1) == 9); + } + CHECK(with_gap(1, 0) == source(1, 0)); + CHECK(with_gap(5, 0) == source(2, 0)); +} + +TEST_CASE("Split counts between edge pixels and their gaps", "[GapPixels]") { + GapPixels gaps; + gaps.split_counts = true; + gaps.seed = 42; + + SECTION("one boundary, integer counts") { + // columns 254..258 and rows 0..1: one chip boundary at 256 + const ROI roi{254, 258, 0, 2}; + NDArray source({2, 4}); + source(0, 0) = 10; + source(0, 1) = 7; // odd, splits 3 + 4 in random order + source(0, 2) = 9; + source(0, 3) = 12; + source(1, 0) = 1; + source(1, 1) = 8; + source(1, 2) = 2; + source(1, 3) = 3; + auto result = insert_gap_pixels( + NDView(source.view()), roi, gaps); + REQUIRE(result.shape() == std::array{2, 6}); + for (ssize_t row = 0; row < 2; ++row) { + CAPTURE(row); + CHECK(result(row, 0) == source(row, 0)); + CHECK(result(row, 5) == source(row, 3)); + // left edge pixel and its gap pixel to the right + CHECK(result(row, 1) + result(row, 2) == source(row, 1)); + CHECK(std::max(result(row, 1), result(row, 2)) == + (source(row, 1) + 1) / 2); + // right edge pixel and its gap pixel to the left + CHECK(result(row, 3) + result(row, 4) == source(row, 2)); + CHECK(std::max(result(row, 3), result(row, 4)) == + (source(row, 2) + 1) / 2); + } + CHECK(result(0, 1) + result(0, 2) == 7); + CHECK(result(1, 1) == 4); + CHECK(result(1, 2) == 4); + CHECK(result(1, 3) == 1); + CHECK(result(1, 4) == 1); + } + SECTION("corner pixels split four ways") { + const ROI roi{255, 257, 255, 257}; + NDArray source({2, 2}); + source(0, 0) = 9; + source(0, 1) = 8; + source(1, 0) = 4; + source(1, 1) = 1; + auto result = insert_gap_pixels( + NDView(source.view()), roi, gaps); + REQUIRE(result.shape() == std::array{4, 4}); + auto quad_sum = [&](ssize_t row, ssize_t col) { + return result(row, col) + result(row, col + 1) + + result(row + 1, col) + result(row + 1, col + 1); + }; + CHECK(quad_sum(0, 0) == 9); + CHECK(quad_sum(0, 2) == 8); + CHECK(quad_sum(2, 0) == 4); + CHECK(quad_sum(2, 2) == 1); + CHECK(result(0, 2) + result(0, 3) == 4); + CHECK(result(1, 2) + result(1, 3) == 4); + CHECK(result(2, 0) == 1); + CHECK(result(2, 1) == 1); + CHECK(result(3, 0) == 1); + CHECK(result(3, 1) == 1); + CHECK(total(NDView(result.view())) == 22); + } + SECTION("floating point halves exactly") { + const ROI roi{255, 257, 255, 257}; + NDArray source({2, 2}, 1.0); + source(0, 0) = 3.0; + auto result = insert_gap_pixels(NDView(source.view()), + roi, gaps); + for (ssize_t row = 0; row < 2; ++row) { + for (ssize_t col = 0; col < 2; ++col) { + CHECK(result(row, col) == 0.75); + } + } + CHECK(result(0, 2) == 0.25); + CHECK(result(3, 3) == 0.25); + } + SECTION("totals, remainders and reproducibility on a module") { + const ROI roi{0, 1024, 0, 512}; + NDArray source({512, 1024}); + for (ssize_t row = 0; row < 512; ++row) { + for (ssize_t col = 0; col < 1024; ++col) { + source(row, col) = + static_cast((row * 7 + col * 13) % 1000 + 1); + } + } + auto first = insert_gap_pixels(NDView(source.view()), + roi, gaps); + auto second = insert_gap_pixels( + NDView(source.view()), roi, gaps); + REQUIRE(first.shape() == std::array{514, 1030}); + CHECK(total(NDView(first.view())) == + total(NDView(source.view()))); + CHECK(first == second); + + // remainders of odd left edge pixels should land on either side + // about equally often + size_t odd = 0; + size_t kept_high = 0; + for (ssize_t row = 0; row < 512; ++row) { + if (row == 255 || row == 256) { + continue; // rows that are split again vertically + } + for (ssize_t boundary : {256, 512, 768}) { + const auto value = source(row, boundary - 1); + if (value % 2 == 0) { + continue; + } + ++odd; + const ssize_t gapped_row = row < 256 ? row : row + 2; + // gapped column of the left edge pixel, source column boundary + // - 1 + const ssize_t left = boundary - 1 + 2 * (boundary / 256 - 1); + if (first(gapped_row, left) > first(gapped_row, left + 1)) { + ++kept_high; + } + } + } + REQUIRE(odd > 500); + const double fraction = static_cast(kept_high) / odd; + CHECK(fraction > 0.4); + CHECK(fraction < 0.6); + + GapPixels other = gaps; + other.seed = 43; + auto third = insert_gap_pixels(NDView(source.view()), + roi, other); + CHECK(third != first); + } + SECTION("module gaps keep the fill value") { + GapPixels module_gaps = gaps; + module_gaps.module_gaps = ModuleGaps{8, 36}; + module_gaps.fill_value = 3; + const ROI roi{0, 1024, 0, 1024}; + NDArray source({1024, 1024}, uint16_t{7}); + auto result = insert_gap_pixels( + NDView(source.view()), roi, module_gaps); + REQUIRE(result.shape() == std::array{1064, 1030}); + for (ssize_t row = 514; row < 550; ++row) { + for (ssize_t col = 0; col < 1030; ++col) { + if (result(row, col) != 3) { + CAPTURE(row, col); + REQUIRE(result(row, col) == 3); + } + } + } + // the edge pixels next to the module gap are not split + CHECK(result(513, 0) == 7); + CHECK(result(550, 0) == 7); + // chip gaps are split: 7 -> 3 + 4 + CHECK(result(0, 255) + result(0, 256) == 7); + CHECK(total(NDView(result.view())) == + 7.0 * 1024 * 1024 + 3.0 * 36 * 1030); + } + SECTION("no split at a ROI border") { + const ROI roi{256, 260, 0, 1}; + NDArray source({1, 4}, uint16_t{9}); + auto result = insert_gap_pixels( + NDView(source.view()), roi, gaps); + REQUIRE(result.shape() == std::array{1, 4}); + CHECK(result(0, 0) == 9); + } +} diff --git a/src/NDView.test.cpp b/src/NDView.test.cpp index 59d231b3..808cb7a1 100644 --- a/src/NDView.test.cpp +++ b/src/NDView.test.cpp @@ -16,7 +16,7 @@ using aare::ssize_t; static_assert(std::is_convertible_v, NDView>); static_assert(!std::is_convertible_v, NDView>); -TEST_CASE("Calculate size from a shape") { +TEST_CASE("Calculate size from a shape", "[NDView]") { Shape<3> shape{2, 3, 4}; REQUIRE(num_elements(shape) == 24); @@ -34,7 +34,7 @@ TEST_CASE("Calculate size from a shape") { REQUIRE(num_elements(shape5) == 10485760000); } -TEST_CASE("Element reference 1D") { +TEST_CASE("Element reference 1D", "[NDView]") { std::vector vec; for (int i = 0; i != 10; ++i) { vec.push_back(i); @@ -47,7 +47,7 @@ TEST_CASE("Element reference 1D") { } } -TEST_CASE("Assign elements through () and []") { +TEST_CASE("Assign elements through () and []", "[NDView]") { std::vector vec; for (int i = 0; i != 10; ++i) { vec.push_back(i); @@ -80,7 +80,7 @@ TEST_CASE("Assign elements through () and []") { REQUIRE(data[9] == 9); } -TEST_CASE("Element reference 1D with a const NDView") { +TEST_CASE("Element reference 1D with a const NDView", "[NDView]") { std::vector vec; for (int i = 0; i != 10; ++i) { vec.push_back(i); @@ -93,7 +93,7 @@ TEST_CASE("Element reference 1D with a const NDView") { } } -TEST_CASE("Convert a mutable NDView to a const NDView") { +TEST_CASE("Convert a mutable NDView to a const NDView", "[NDView]") { std::vector vec{1, 2, 3, 4}; NDView mutable_view(vec.data(), Shape<2>{2, 2}); @@ -104,7 +104,7 @@ TEST_CASE("Convert a mutable NDView to a const NDView") { REQUIRE(const_view.strides() == mutable_view.strides()); } -TEST_CASE("Element reference 2D") { +TEST_CASE("Element reference 2D", "[NDView]") { std::vector vec(12); std::iota(vec.begin(), vec.end(), 0); @@ -120,7 +120,7 @@ TEST_CASE("Element reference 2D") { } } -TEST_CASE("Element reference 3D") { +TEST_CASE("Element reference 3D", "[NDView]") { std::vector vec; for (int i = 0; i != 24; ++i) { vec.push_back(i); @@ -139,7 +139,7 @@ TEST_CASE("Element reference 3D") { } } -TEST_CASE("Plus and minus with single value") { +TEST_CASE("Plus and minus with single value", "[NDView]") { std::vector vec(12); std::iota(vec.begin(), vec.end(), 0); NDView data(vec.data(), Shape<2>{3, 4}); @@ -161,7 +161,7 @@ TEST_CASE("Plus and minus with single value") { } } -TEST_CASE("Multiply and divide with single value") { +TEST_CASE("Multiply and divide with single value", "[NDView]") { std::vector vec; for (int i = 0; i != 12; ++i) { vec.push_back(i); @@ -185,7 +185,7 @@ TEST_CASE("Multiply and divide with single value") { } } -TEST_CASE("elementwise assign") { +TEST_CASE("elementwise assign", "[NDView]") { std::vector vec(25); NDView data(vec.data(), Shape<2>{5, 5}); @@ -195,7 +195,7 @@ TEST_CASE("elementwise assign") { } } -TEST_CASE("iterators") { +TEST_CASE("iterators", "[NDView]") { std::vector vec(12); std::iota(vec.begin(), vec.end(), 0); NDView data(vec.data(), Shape<1>{12}); @@ -220,7 +220,7 @@ TEST_CASE("iterators") { } } -TEST_CASE("divide with another NDView") { +TEST_CASE("divide with another NDView", "[NDView]") { std::vector vec0{9, 12, 3}; std::vector vec1{3, 2, 1}; std::vector result{3, 6, 3}; @@ -237,7 +237,7 @@ TEST_CASE("divide with another NDView") { } } -TEST_CASE("Retrieve shape") { +TEST_CASE("Retrieve shape", "[NDView]") { std::vector vec(12); std::iota(vec.begin(), vec.end(), 0); NDView data(vec.data(), Shape<2>{3, 4}); @@ -245,7 +245,7 @@ TEST_CASE("Retrieve shape") { REQUIRE(data.shape()[1] == 4); } -TEST_CASE("compare two views") { +TEST_CASE("compare two views", "[NDView]") { std::vector vec1(12); std::iota(vec1.begin(), vec1.end(), 0); NDView view1(vec1.data(), Shape<2>{3, 4}); @@ -257,7 +257,7 @@ TEST_CASE("compare two views") { REQUIRE((view1 == view2)); } -TEST_CASE("Compare two views with different size") { +TEST_CASE("Compare two views with different size", "[NDView]") { std::vector vec1(12); std::iota(vec1.begin(), vec1.end(), 0); NDView view1(vec1.data(), Shape<2>{3, 4}); @@ -269,7 +269,7 @@ TEST_CASE("Compare two views with different size") { REQUIRE_FALSE(view1 == view2); } -TEST_CASE("Compare two views with same size but different shape") { +TEST_CASE("Compare two views with same size but different shape", "[NDView]") { std::vector vec1(12); std::iota(vec1.begin(), vec1.end(), 0); NDView view1(vec1.data(), Shape<2>{3, 4}); @@ -281,7 +281,7 @@ TEST_CASE("Compare two views with same size but different shape") { REQUIRE_FALSE(view1 == view2); } -TEST_CASE("Create a view over a vector") { +TEST_CASE("Create a view over a vector", "[NDView]") { std::vector vec(12); std::iota(vec.begin(), vec.end(), 0); auto v = aare::make_view(vec); @@ -290,9 +290,26 @@ TEST_CASE("Create a view over a vector") { REQUIRE(v[11] == 11); } -TEST_CASE("NDView over byte") { +TEST_CASE("NDView over byte", "[NDView]") { std::vector buf(5); auto v = aare::make_view(buf); REQUIRE(v.shape()[0] == 5); REQUIRE(v[0] == std::byte{0}); } + +TEST_CASE("NDView subview", "[NDView]") { + std::vector vec(12); + std::iota(vec.begin(), vec.end(), 0); + NDView data(vec.data(), Shape<2>{3, 4}); + + auto subview = data.sub_view(1); + REQUIRE(subview.shape().size() == 1); + REQUIRE(subview.shape()[0] == 4); + CHECK(subview(0) == 4); + CHECK(subview(1) == 5); + CHECK(subview(2) == 6); + CHECK(subview(3) == 7); + + REQUIRE_THROWS_AS(data.sub_view(3), std::runtime_error); + REQUIRE_THROWS_AS(data.sub_view(-1), std::runtime_error); +} diff --git a/src/RawFile.cpp b/src/RawFile.cpp index e416a27e..34c336a2 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,69 @@ 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) { + const auto detector = m_master.detector_type(); + if (detector != DetectorType::Jungfrau && + detector != DetectorType::Eiger) { + throw std::invalid_argument( + LOCATION + "Gap pixels are only supported for Jungfrau " + "and Eiger detectors"); + } + if (gap_pixels->split_counts && + detector == DetectorType::Jungfrau) { + throw std::invalid_argument( + LOCATION + "split_counts cannot be used on raw Jungfrau " + "data: it would split the gain bits. Convert " + "to energy and use insert_gap_pixels"); + } + // bit depth and fill value are known at open, reject them here + // rather than on the first read + switch (m_master.bitdepth()) { + case 8: + detail::fill_as(gap_pixels->fill_value); + break; + case 16: + detail::fill_as(gap_pixels->fill_value); + break; + case 32: + detail::fill_as(gap_pixels->fill_value); + break; + default: + throw std::invalid_argument( + LOCATION + + fmt::format("Gap pixels are not supported for a bit " + "depth of {}", + m_master.bitdepth())); + } + m_gap_pixels = gap_pixels; + m_gap_layout = detail::gap_layout(*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 +257,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 +270,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 +301,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 +312,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 +385,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 +395,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 +447,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..a4d67ec9 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,354 @@ 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, + const GapPixels &gaps) { + auto expected = + insert_gap_pixels(NDView(plain.view()), roi, 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, + const GapPixels &gaps) { + switch (plain.bitdepth()) { + case 8: + return matches_gap_oracle(gapped, plain, roi, gaps); + case 16: + return matches_gap_oracle(gapped, plain, roi, gaps); + case 32: + return matches_gap_oracle(gapped, plain, roi, 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(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]") { + // one Eiger module: four 256x512 ports in a 2x2 layout. Splitting is + // rejected on raw Jungfrau data, so the synthetic file is Eiger. + TemporaryRawFiles files(xy{2, 2}, xy{256, 512}); + auto metadata = nlohmann::json::parse(std::ifstream(files.master_path())); + metadata["Detector Type"] = "Eiger"; + std::ofstream(files.master_path()) << metadata; + 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, 1024, 0, 512}; + REQUIRE(matches_gap_oracle(gapped_frame, plain_frame, detector, 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, + 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 split counts on raw Jungfrau data", + "[RawFile][GapPixels]") { + TemporaryRawFiles files(xy{1, 1}, xy{512, 1024}); + GapPixels gaps; + gaps.split_counts = true; + REQUIRE_NOTHROW(RawFile(files.master_path(), "r", GapPixels{})); + REQUIRE_THROWS_AS(RawFile(files.master_path(), "r", gaps), + std::invalid_argument); +} + +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", + {1028, 1030}, + {}}, + {"raw/virtual/920/jungfrau/two_modules_master_0.json", + {1064, 1030}, + ModuleGaps{8, 36}}, + {"raw/ROITestData/MultipleChipROI/run_master_2.json", {301, 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(), + gaps)); +} + +TEST_CASE("Split counts on a recorded Eiger file", + "[.with-data][RawFile][GapPixels]") { + const auto fpath = + test_data_path() / "raw/eiger/eiger_500k_16bit_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(), + gaps)); +} + +TEST_CASE("Split counts are rejected 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; + REQUIRE_THROWS_AS(RawFile(fpath, "r", gaps), std::invalid_argument); +} + +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], gaps)); + } +} + TEST_CASE("Read number of frames from a jungfrau raw file", "[.with-data][RawFile]") { @@ -535,7 +884,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 +906,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); }