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); + } +}