From f43e48cb901049c9adfb44ea5f9628f938aa69d7 Mon Sep 17 00:00:00 2001 From: AliceMazzoleni99 Date: Mon, 12 Oct 2026 17:21:25 +0200 Subject: [PATCH 1/5] added subview (#419) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - adds subview(ssize_t index) opeartor - returns subview with dropped dimension --------- Co-authored-by: Erik Fröjdh --- RELEASE.md | 1 + include/aare/NDView.hpp | 15 ++++++++++++ src/NDView.test.cpp | 53 +++++++++++++++++++++++++++-------------- 3 files changed, 51 insertions(+), 18 deletions(-) diff --git a/RELEASE.md b/RELEASE.md index 86d00ed2..adb8405a 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -10,6 +10,7 @@ - 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 ## 2026.9.30 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/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); +} From 4408fe95b084ef87a808413d2d00d510e4d04d8a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Erik=20Fr=C3=B6jdh?= Date: Tue, 13 Oct 2026 15:00:45 +0200 Subject: [PATCH 2/5] Fixed failing tests due to updated bindings (#420) - New python bindings return a (1,n_bytes) array instead of a 1d array - Flatten on python side before transform --- RELEASE.md | 7 +++++++ python/aare/transform.py | 4 ++-- python/tests/test_Transformation.py | 31 +++++++++++++++++++++++++++++ 3 files changed, 40 insertions(+), 2 deletions(-) diff --git a/RELEASE.md b/RELEASE.md index adb8405a..907c2668 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -12,6 +12,13 @@ wheels. - added ``NDView::sub_view(index)`` operator returns a subview at index along dimension 0 +### 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. + ## 2026.9.30 ### New Features: 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/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))) From 5ccda7955d5fff9e0d2b3a4bc6d63383e998abf8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Erik=20Fr=C3=B6jdh?= Date: Tue, 13 Oct 2026 17:24:05 +0200 Subject: [PATCH 3/5] Assemble raw frames through one placement path (#421) First part of splitting up dev/gap. No gap pixels in here. - Extract the frame number sync loop into `RawFile::synchronized_frame_indices`. - Replace the two placement branches keyed on the detector layout with one loop. A part as wide as the ROI is read straight into the frame; a narrower part goes through a scratch buffer sized once at open, instead of a buffer allocated per frame. - The "x pos not 0" implementation error is gone, since the offset is computed for both paths. --- RELEASE.md | 5 + include/aare/RawFile.hpp | 14 +++ src/RawFile.cpp | 213 ++++++++++++++------------------------- 3 files changed, 97 insertions(+), 135 deletions(-) diff --git a/RELEASE.md b/RELEASE.md index 907c2668..cd35a69a 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -12,6 +12,11 @@ wheels. - added ``NDView::sub_view(index)`` operator returns a subview at index along dimension 0 +### 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 diff --git a/include/aare/RawFile.hpp b/include/aare/RawFile.hpp index 4ef92962..e2666577 100644 --- a/include/aare/RawFile.hpp +++ b/include/aare/RawFile.hpp @@ -26,6 +26,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{}; @@ -154,6 +158,16 @@ 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 read the frame at the given frame index into the image buffer * @param frame_index frame number to read diff --git a/src/RawFile.cpp b/src/RawFile.cpp index e416a27e..54adada6 100644 --- a/src/RawFile.cpp +++ b/src/RawFile.cpp @@ -9,6 +9,7 @@ #include "aare/logger.hpp" #include +#include #include #include @@ -49,13 +50,17 @@ 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); if (!min_frames) { throw std::runtime_error(fmt::format( "No raw subfiles selected by '{}'", fname.string())); @@ -338,6 +343,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 +395,45 @@ 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); - } + const size_t bytes_per_pixel = m_master.bitdepth() / bits_per_byte; + const size_t roi_width = roi.pixels_x(); - // 1. if frame number vector is the same break - while (!all_equal(frame_numbers)) { + 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]); - // 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())); + std::byte *dest = + frame_buffer + (static_cast(pos.origin_y) * roi_width + + static_cast(pos.origin_x)) * + bytes_per_pixel; - // 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_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]; - - // 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; - - 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; - } - - } 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; + } } } From d1e432f8cabccf3693c7bb7b19fe53afac1d41d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Erik=20Fr=C3=B6jdh?= Date: Thu, 15 Oct 2026 10:31:45 +0200 Subject: [PATCH 4/5] Add gap pixel insertion for images in memory (#423) Second part of splitting up dev/gap. Only the in-memory transform; no reader changes. - `GapPixels` and `ModuleGaps` configure 2 pixel gaps at the chip boundaries inside each Jungfrau or Eiger module and, only when module gaps are set, gaps between modules. Gap pixels get a fill value or, with `split_counts`, half of the neighbouring double size pixel. - `insert_gap_pixels` and `gapped_shape` apply the transform to an image in memory. Gap positions follow detector coordinates, so a ROI receives exactly the gaps that fall inside it. - Python bindings, tests and docs for the above. Chip and module sizes are fixed at 256 and 1024x512, so there is no detector type or quad argument: an Eiger quad is 512x512 and its only chip boundaries are at pixel 256. Reading raw files with gaps comes in the next PR. --- CMakeLists.txt | 3 + RELEASE.md | 7 + docs/src/GapPixels.rst | 101 +++++ docs/src/index.rst | 1 + docs/src/python/file/index.rst | 1 + docs/src/python/file/pyGapPixels.rst | 72 ++++ include/aare/GapPixels.hpp | 379 +++++++++++++++++ python/aare/__init__.py | 1 + python/src/bind_GapPixels.hpp | 174 ++++++++ python/src/module.cpp | 2 + python/tests/test_GapPixels.py | 130 ++++++ src/GapPixels.cpp | 104 +++++ src/GapPixels.test.cpp | 588 +++++++++++++++++++++++++++ 13 files changed, 1563 insertions(+) create mode 100644 docs/src/GapPixels.rst create mode 100644 docs/src/python/file/pyGapPixels.rst create mode 100644 include/aare/GapPixels.hpp create mode 100644 python/src/bind_GapPixels.hpp create mode 100644 python/tests/test_GapPixels.py create mode 100644 src/GapPixels.cpp create mode 100644 src/GapPixels.test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 9a6ba241..62a2a368 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -388,6 +388,7 @@ set(PUBLICHEADERS include/aare/ROI.hpp include/aare/ROIGeometry.hpp include/aare/DetectorGeometry.hpp + include/aare/GapPixels.hpp include/aare/JungfrauDataFile.hpp include/aare/logger.hpp include/aare/NDArray.hpp @@ -413,6 +414,7 @@ set(SourceFiles ${CMAKE_CURRENT_SOURCE_DIR}/src/ROI.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ROIGeometry.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/GapPixels.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/hist/PedestalTrackingPixelHistogram.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp @@ -482,6 +484,7 @@ if(AARE_TESTS) ${CMAKE_CURRENT_SOURCE_DIR}/src/fit/Fit.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/GapPixels.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolation.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp diff --git a/RELEASE.md b/RELEASE.md index cd35a69a..0cd63bae 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -11,12 +11,19 @@ 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. ### 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 diff --git a/docs/src/GapPixels.rst b/docs/src/GapPixels.rst new file mode 100644 index 00000000..0470878d --- /dev/null +++ b/docs/src/GapPixels.rst @@ -0,0 +1,101 @@ +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. +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. + +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/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..020144c3 --- /dev/null +++ b/docs/src/python/file/pyGapPixels.rst @@ -0,0 +1,72 @@ +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, seed=1) + 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. + +.. autoclass:: GapPixels + :members: + :undoc-members: + +.. autoclass:: ModuleGaps + :members: + :undoc-members: + +.. autofunction:: insert_gap_pixels + +.. autofunction:: gapped_shape 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/python/aare/__init__.py b/python/aare/__init__.py index ea3f4c39..6394489c 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -25,6 +25,7 @@ from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarCluster from ._aare import DetectorType, FrameDiscardPolicy, ReadoutMode from ._aare import hitmap from ._aare import ROI +from ._aare import GapPixels, ModuleGaps, insert_gap_pixels, gapped_shape from ._aare import xy from ._aare import corner diff --git a/python/src/bind_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/module.cpp b/python/src/module.cpp index b83deaca..0df7fe8d 100644 --- a/python/src/module.cpp +++ b/python/src/module.cpp @@ -14,6 +14,7 @@ #include "bind_Defs.hpp" #include "bind_Eta.hpp" #include "bind_FastPedestal.hpp" +#include "bind_GapPixels.hpp" #include "bind_Interpolator.hpp" #include "bind_MultiThreadedFileReader.hpp" #include "bind_Pedestal.hpp" @@ -156,6 +157,7 @@ PYBIND11_MODULE(_aare, m) { register_calculate_3x3eta(m); define_defs_bindings(m); + define_gap_pixels_bindings(m); using Sum_index_pair_d = Sum_index_pair; PYBIND11_NUMPY_DTYPE(Sum_index_pair_d, sum, index); diff --git a/python/tests/test_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/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); + } +} From 7abbf3cbb5a945374d11adf80089050ae55b10e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Erik=20Fr=C3=B6jdh?= Date: Mon, 26 Oct 2026 14:39:13 +0100 Subject: [PATCH 5/5] Read raw files with gap pixels (#426) Third part of splitting up dev/gap. The multithreaded reader comes in the next PR. - `RawFile` takes an optional `GapPixels`. Frames are assembled with the gaps inserted: fill the gap lines, scatter each module part through the scratch buffer into its gapped position, and split the chip edge counts if configured. `rows()`, `cols()`, `pixels_per_frame()` and `bytes_per_frame()` include the gaps, for the full frame and for each ROI. Only Jungfrau and Eiger are supported. - Python `RawFile` accepts `gap_pixels` as a bool or a `GapPixels`, exposes it as a property, and takes frame shapes from `rows()` and `cols()`. - The generic `File` is unchanged: it does not know the detector geometry of `.npy` and `.dat` files, so gap pixels stay a `RawFile` option. Use `RawFile`, or apply `insert_gap_pixels` to frames from `File`. - Tests: generalized raw file helper with arbitrary module layouts, module layout assembly tests (these cover #421, which went in without tests), synthetic gap tests, a split-counts oracle against `insert_gap_pixels`, unsupported detector rejection, and data-backed recorded-file and two-ROI cases. - A raw file benchmark with no-gaps, fill-value and split-counts cases. - Docs and release note. --- RELEASE.md | 4 + benchmarks/CMakeLists.txt | 2 +- benchmarks/rawfile_benchmark.cpp | 112 +++++++++ docs/src/GapPixels.rst | 21 ++ docs/src/RawFile.rst | 8 + docs/src/python/file/pyGapPixels.rst | 24 +- docs/src/python/file/pyRawFile.rst | 6 + include/aare/RawFile.hpp | 36 ++- python/aare/RawFile.py | 14 +- python/src/bind_RawFile.hpp | 43 +++- python/tests/test_RawFile.py | 158 ++++++++++++ src/RawFile.cpp | 133 +++++++++- src/RawFile.test.cpp | 353 ++++++++++++++++++++++++++- tests/raw_file_helpers.hpp | 77 +++++- 14 files changed, 955 insertions(+), 36 deletions(-) create mode 100644 benchmarks/rawfile_benchmark.cpp diff --git a/RELEASE.md b/RELEASE.md index 0cd63bae..e8b90b5f 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -17,6 +17,10 @@ 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: diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index e874a003..f89b57d3 100644 --- a/benchmarks/CMakeLists.txt +++ b/benchmarks/CMakeLists.txt @@ -20,7 +20,7 @@ add_executable(benchmarks) target_sources( benchmarks PRIVATE ndarray_benchmark.cpp ndview_benchmark.cpp calculateeta_benchmark.cpp - reduce_benchmark.cpp clusterfile_benchmark.cpp) + reduce_benchmark.cpp clusterfile_benchmark.cpp rawfile_benchmark.cpp) # Link Google Benchmark and other necessary libraries target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core diff --git a/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 index 0470878d..f5169e7a 100644 --- a/docs/src/GapPixels.rst +++ b/docs/src/GapPixels.rst @@ -17,6 +17,7 @@ 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. @@ -68,6 +69,26 @@ 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 ---------------- 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/python/file/pyGapPixels.rst b/docs/src/python/file/pyGapPixels.rst index 020144c3..58e344a3 100644 --- a/docs/src/python/file/pyGapPixels.rst +++ b/docs/src/python/file/pyGapPixels.rst @@ -46,7 +46,7 @@ defaults and ``(1064, 1030)`` with ``ModuleGaps(8, 36)``. 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, seed=1) + 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 @@ -59,6 +59,28 @@ defaults and ``(1064, 1030)`` with ``ModuleGaps(8, 36)``. 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: 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/RawFile.hpp b/include/aare/RawFile.hpp index e2666577..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 { @@ -38,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; @@ -106,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. @@ -168,6 +192,16 @@ class RawFile : public FileInterface { 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/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/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/src/RawFile.cpp b/src/RawFile.cpp index 54adada6..34c336a2 100644 --- a/src/RawFile.cpp +++ b/src/RawFile.cpp @@ -24,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; @@ -61,6 +62,60 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode) } } 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) { throw std::runtime_error(fmt::format( "No raw subfiles selected by '{}'", fname.string())); @@ -202,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() { @@ -217,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(); } @@ -249,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) { @@ -260,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(); } @@ -333,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); @@ -400,6 +452,30 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, const auto frame_indices = synchronized_frame_indices(frame_index, roi_index); + 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; + } + const size_t bytes_per_pixel = m_master.bitdepth() / bits_per_byte; const size_t roi_width = roi.pixels_x(); @@ -437,6 +513,43 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer, } } +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); + } +} + std::vector RawFile::read_n(size_t n_frames) { // TODO: implement this in a more efficient way if (num_rois() > 1) { 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); }