first attempt at gappixels

This commit is contained in:
Erik Frojdh committed 2026-10-08 15:36:46 +02:00
1 parent bf0aa4a3bf
commit 406aade696
27 files changed
+995 -186

No files matched your search

+3
View File
@@ -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
+15
View File
@@ -4,6 +4,15 @@
### New Features:
- Added gap pixel insertion for Jungfrau and Eiger raw files. ``RawFile``,
``File`` (through ``FileConfig::gap_pixels``) and
``MultiThreadedFileReader`` take a ``GapPixels`` configuration, in Python
the ``gap_pixels`` argument, and assemble frames with 2 pixel chip gaps
and optional ``ModuleGaps`` between modules. Gap positions follow
detector coordinates before any ROI. Gap pixels receive a fill value or,
with ``split_counts``, half of the neighbouring double size pixel. The
new ``insert_gap_pixels`` and ``gapped_shape`` apply the same transform
to images in memory.
- Added the Python ``xy`` type with ``row`` and ``col`` fields, enabling access
to ``RawFile.geometry``, ``RawMasterFile.detector_layout``, and
``RawMasterFile.udp_interfaces_per_module``. It unpacks as ``row, col``.
@@ -11,6 +20,12 @@
build, so releases publish them to PyPI alongside the Linux and macOS
wheels.
### Other:
- ``RawFile`` assembles multi-module frames through one placement path with
a scratch buffer allocated once at open, instead of allocating a module
buffer for every frame.
## 2026.9.30
### New Features:
+1 -1
View File
@@ -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
+5
View File
@@ -1,6 +1,11 @@
File
=============
``File`` opens ``.json`` and ``.raw`` master files, ``.npy`` files and
Jungfrau ``.dat`` files behind one interface. ``FileConfig::gap_pixels``
inserts :doc:`gap pixels <GapPixels>` when reading raw files; other file
types reject a set value with ``std::invalid_argument``.
.. doxygenclass:: aare::File
:members:
+5
View File
@@ -25,6 +25,11 @@ an allocation when the caller already owns a buffer of at least
remaining from the current position, and ``seek()`` to reposition the reader.
Call ``close()`` to release all worker file handles early.
The optional fifth argument inserts :doc:`gap pixels <GapPixels>` when
reading raw files. ``rows()``, ``cols()`` and ``bytes_per_frame()`` then
include the gaps. Each worker reader uses the configured seed plus its
worker index for split counts.
.. note::
Multiple workers do not guarantee faster reads. Performance depends on the
+8
View File
@@ -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 <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:
+1
View File
@@ -44,6 +44,7 @@ AARE
NDView
Frame
File
GapPixels
MultiThreadedFileReader
Dtype
Cluster
@@ -24,6 +24,11 @@ iteration position. The context manager closes all worker files on exit.
``close()`` is also available for explicit cleanup and may be called
repeatedly.
``gap_pixels`` inserts :doc:`gap pixels <../file/pyGapPixels>` when reading
raw files. Pass ``True`` for the defaults or a ``GapPixels``. ``rows`` and
``cols`` then include the gaps. Each worker uses the configured seed plus its
worker index for split counts.
.. autoclass:: MultiThreadedFileReader
:members:
:undoc-members:
+1
View File
@@ -8,6 +8,7 @@ File I/O
pyClusterFile
pyCtbRawFile
pyFile
pyGapPixels
pyJungfrauDataFile
pyRawFile
pyRawMasterFile
+6
View File
@@ -4,6 +4,12 @@ File
.. py:currentmodule:: aare
``File`` opens raw master files, ``.npy`` files and Jungfrau ``.dat`` files
behind one interface. ``gap_pixels`` inserts
:doc:`gap pixels <pyGapPixels>` when reading raw files; pass ``True`` for
the defaults or a ``GapPixels``. Other file types raise ``ValueError`` when
it is set.
.. autoclass:: File
:members:
:undoc-members:
+6
View File
@@ -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 <pyGapPixels>` 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:
+9 -1
View File
@@ -2,9 +2,11 @@
#pragma once
#include "aare/Dtype.hpp"
#include "aare/Frame.hpp"
#include "aare/GapPixels.hpp"
#include "aare/defs.hpp"
#include <filesystem>
#include <optional>
#include <vector>
namespace aare {
@@ -15,6 +17,7 @@ namespace aare {
* rows: number of rows in the file
* cols: number of columns in the file
* geometry: geometry of the file
* gap_pixels: insert gap pixels when reading raw files
*/
struct FileConfig {
aare::Dtype dtype{typeid(uint16_t)};
@@ -24,7 +27,8 @@ struct FileConfig {
return dtype == other.dtype && rows == other.rows &&
cols == other.cols && geometry == other.geometry &&
detector_type == other.detector_type &&
max_frames_per_file == other.max_frames_per_file;
max_frames_per_file == other.max_frames_per_file &&
gap_pixels == other.gap_pixels;
}
bool operator!=(const FileConfig &other) const { return !(*this == other); }
@@ -34,6 +38,10 @@ struct FileConfig {
DetectorType detector_type{DetectorType::Unknown};
int max_frames_per_file{};
size_t total_frames{};
/// @brief Insert gap pixels into every frame. Only RawFile supports
/// this; other readers reject a set value.
std::optional<GapPixels> gap_pixels{};
// std::string to_string() const {
// return "{ dtype: " + dtype.to_string() +
// ", rows: " + std::to_string(rows) +
+10 -1
View File
@@ -3,6 +3,7 @@
#include "aare/Dtype.hpp"
#include "aare/File.hpp"
#include "aare/GapPixels.hpp"
#include <cstddef>
#include <filesystem>
@@ -25,10 +26,13 @@ class MultiThreadedFileReader {
* @param n_threads maximum number of worker threads
* @param chunk_size number of frames claimed by a worker at a time
* @param total_frames number of frames to read, or all frames when omitted
* @param gap_pixels insert gap pixels, raw files only. Each worker
* reader uses the configured seed plus its worker index.
*/
MultiThreadedFileReader(std::filesystem::path fname, size_t n_threads,
size_t chunk_size,
std::optional<size_t> total_frames = std::nullopt);
std::optional<size_t> total_frames = std::nullopt,
std::optional<GapPixels> gap_pixels = std::nullopt);
MultiThreadedFileReader(const MultiThreadedFileReader &) = delete;
MultiThreadedFileReader &
@@ -66,6 +70,9 @@ class MultiThreadedFileReader {
bool is_open() const noexcept { return !m_files.empty(); }
size_t n_threads() const noexcept { return m_n_threads; }
const std::optional<GapPixels> &gap_pixels() const noexcept {
return m_gap_pixels;
}
size_t chunk_size() const noexcept { return m_chunk_size; }
size_t total_frames() const noexcept { return m_total_frames; }
size_t source_total_frames() const noexcept {
@@ -96,9 +103,11 @@ class MultiThreadedFileReader {
size_t m_bytes_per_frame;
size_t m_total_bytes;
size_t m_current_frame;
std::optional<GapPixels> m_gap_pixels;
std::vector<File> m_files;
void ensure_open() const;
FileConfig worker_config(size_t worker_index) const;
};
} // namespace aare::experimental
+49 -3
View File
@@ -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 <array>
#include <optional>
#include <random>
namespace aare {
@@ -29,20 +32,39 @@ class RawFile : public FileInterface {
RawMasterFile m_master;
size_t m_current_frame{};
std::vector<ROIGeometry> m_ROI_geometries;
/// @brief Minimum frame count across the selected raw subfile series.
size_t m_frames_in_file{};
/// @brief Scratch space for one module part that is not contiguous in
/// the assembled frame. Sized at open to the largest part in any ROI.
std::vector<std::byte> m_part_buffer;
/// @brief Gap pixel configuration, empty when reading without gaps
std::optional<GapPixels> 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<ROI> m_roi_rects;
/// @brief {rows, cols} of the assembled frame of each ROI, including
/// gap pixels when enabled
std::vector<std::array<ssize_t, 2>> 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<GapPixels> gap_pixels = std::nullopt);
virtual ~RawFile() override = default;
Frame read_frame() override;
@@ -102,6 +124,10 @@ class RawFile : public FileInterface {
*/
size_t pixels_per_frame(const size_t roi_index);
size_t bytes_per_pixel() const;
/// @brief Gap pixel configuration, empty when reading without gaps
const std::optional<GapPixels> &gap_pixels() const { return m_gap_pixels; }
void seek(size_t frame_index) override;
size_t tell() override;
/// @brief Minimum actual frame count across all subfiles and ROIs.
@@ -154,6 +180,26 @@ class RawFile : public FileInterface {
std::runtime_error frame_error(size_t frame_index,
const std::string &message) const;
/**
* @brief Frame index to read from each module part of a ROI so that all
* parts hold the same detector frame number. Parts whose frame number
* lags are advanced; the index never moves backwards.
* @throws std::runtime_error if a part runs out of frames while
* synchronizing
*/
std::vector<size_t> 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 <typename T>
void assemble_with_gaps(std::byte *frame_buffer, size_t roi_index,
const std::vector<size_t> &frame_indices,
DetectorHeader *header);
/**
* @brief read the frame at the given frame index into the image buffer
* @param frame_index frame number to read
+12 -2
View File
@@ -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
+1
View File
@@ -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
+19 -1
View File
@@ -1,6 +1,7 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/GapPixels.hpp"
#include "aare/MultiThreadedFileReader.hpp"
#include <cstddef>
@@ -79,9 +80,21 @@ inline void define_multi_threaded_file_reader_bindings(py::module_ &m) {
reader.attr("__module__") = "aare.experimental";
reader
.def(py::init<std::filesystem::path, size_t, size_t,
std::optional<size_t>>(),
std::optional<size_t>, std::optional<GapPixels>>(),
py::arg("fname"), py::arg("n_threads"), py::arg("chunk_size"),
py::arg("total_frames") = py::none(),
py::arg("gap_pixels") = py::none())
.def(py::init([](std::filesystem::path fname, size_t n_threads,
size_t chunk_size, std::optional<size_t> total_frames,
bool gap_pixels) {
return MultiThreadedFileReader(
std::move(fname), n_threads, chunk_size, total_frames,
gap_pixels ? std::optional<GapPixels>(GapPixels{})
: std::nullopt);
}),
py::arg("fname"), py::arg("n_threads"), py::arg("chunk_size"),
py::arg("total_frames") = py::none(),
py::arg("gap_pixels") = false,
R"doc(
Read chunks of detector frames concurrently.
@@ -93,7 +106,12 @@ inline void define_multi_threaded_file_reader_bindings(py::module_ &m) {
n_threads: Maximum number of worker threads.
chunk_size: Number of frames read per claimed chunk.
total_frames: Optional frame limit. None reads all frames.
gap_pixels: Insert gap pixels, raw files only. True uses
the default GapPixels configuration. Each worker
uses the configured seed plus its worker index.
)doc")
.def_property_readonly("gap_pixels",
&MultiThreadedFileReader::gap_pixels)
.def(
"read",
[](MultiThreadedFileReader &self) {
+35 -8
View File
@@ -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 <cstdint>
#include <filesystem>
#include <fmt/format.h>
#include <memory>
#include <pybind11/iostream.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
@@ -34,7 +36,33 @@ std::runtime_error raw_file_read_error(RawFile &file,
void define_raw_file_io_bindings(py::module &m) {
py::class_<RawFile>(m, "RawFile")
.def(py::init<const std::filesystem::path &>())
.def(py::init([](const std::filesystem::path &fname,
const GapPixels &gap_pixels) {
return std::make_unique<RawFile>(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<RawFile>(
fname, "r",
gap_pixels ? std::optional<GapPixels>(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<size_t> 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<size_t> 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<size_t> shape{
n_frames, self.roi_geometries(roi_index).pixels_y(),
self.roi_geometries(roi_index).pixels_x()};
std::vector<size_t> 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];
+32 -7
View File
@@ -3,6 +3,7 @@
#include "aare/CtbRawFile.hpp"
#include "aare/File.hpp"
#include "aare/Frame.hpp"
#include "aare/GapPixels.hpp"
#include "aare/ROI.hpp"
#include "aare/RawFile.hpp"
#include "aare/RawMasterFile.hpp"
@@ -29,13 +30,37 @@ void define_file_io_bindings(py::module &m) {
debug, roundRNumber, detType, version, packetMask);
py::class_<File>(m, "File")
.def(py::init([](const std::filesystem::path &fname) {
return File(fname, "r", {});
}))
.def(py::init(
[](const std::filesystem::path &fname, const std::string &mode) {
return File(fname, mode, {});
}))
.def(py::init([](const std::filesystem::path &fname,
const std::string &mode,
std::optional<GapPixels> gap_pixels) {
FileConfig cfg;
cfg.gap_pixels = std::move(gap_pixels);
return File(fname, mode, cfg);
}),
py::arg("fname"), py::arg("mode") = "r",
py::arg("gap_pixels") = py::none())
.def(py::init([](const std::filesystem::path &fname,
const std::string &mode, bool gap_pixels) {
FileConfig cfg;
if (gap_pixels) {
cfg.gap_pixels = GapPixels{};
}
return File(fname, mode, cfg);
}),
py::arg("fname"), py::arg("mode") = "r",
py::arg("gap_pixels") = false, R"(
Open a file for reading.
Parameters
----------
fname : path
Master file (.json, .raw), .npy or .dat file.
mode : str
Only "r" is supported.
gap_pixels : bool or GapPixels
Insert gap pixels when reading raw files. True uses the
default GapPixels configuration.
)")
.def(py::init<const std::filesystem::path &, const std::string &,
const FileConfig &>())
+2
View File
@@ -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<int16_t, 3, 3, uint16_t>(m);
define_defs_bindings(m);
define_gap_pixels_bindings(m);
using Sum_index_pair_d = Sum_index_pair<double, corner>;
PYBIND11_NUMPY_DTYPE(Sum_index_pair_d, sum, index);
@@ -3,7 +3,9 @@ import numpy as np
import pytest
import aare
from aare import GapPixels, RawFile, insert_gap_pixels
from aare.experimental import MultiThreadedFileReader
from test_RawFile import write_module_layout
@pytest.fixture
@@ -145,3 +147,38 @@ def test_invalid_configuration(
chunk_size=chunk_size,
total_frames=total_frames,
)
def test_gap_pixels_on_raw_files(tmp_path):
master_path, expected = write_module_layout(tmp_path, 1, 1, frames=4,
module_rows=512, module_cols=1024)
gaps = GapPixels(fill_value=2)
reader = MultiThreadedFileReader(master_path, n_threads=2, chunk_size=1, gap_pixels=gaps)
assert reader.gap_pixels == gaps
assert (reader.rows, reader.cols) == (514, 1030)
frames = reader.read_all()
assert frames.shape == (4, 514, 1030)
assert (frames == insert_gap_pixels(expected, aare.DetectorType.Jungfrau, gaps)).all()
reader = MultiThreadedFileReader(master_path, n_threads=2, chunk_size=1, gap_pixels=True)
assert reader.gap_pixels == GapPixels()
assert reader.read_all().shape == (4, 514, 1030)
reader = MultiThreadedFileReader(master_path, n_threads=2, chunk_size=1)
assert reader.gap_pixels is None
assert reader.read_all().shape == (4, 512, 1024)
def test_gap_pixels_split_counts_preserve_totals(tmp_path):
master_path, expected = write_module_layout(tmp_path, 1, 1, frames=4,
module_rows=512, module_cols=1024)
gaps = GapPixels(split_counts=True, seed=1)
frames = MultiThreadedFileReader(master_path, n_threads=2, chunk_size=1, gap_pixels=gaps).read_all()
assert frames.shape == (4, 514, 1030)
assert (frames.sum(axis=(1, 2)) == expected.sum()).all()
def test_gap_pixels_rejected_for_numpy_files(frame_file):
path, _ = frame_file
with pytest.raises(ValueError):
MultiThreadedFileReader(path, n_threads=1, chunk_size=1, gap_pixels=True)
+147
View File
@@ -2,6 +2,7 @@
import pytest
import json
from aare import File, RawFile, RawSubFile, DetectorType, ROI, UDPPortPosition, xy
from aare import GapPixels, ModuleGaps, insert_gap_pixels
import numpy as np
@@ -28,6 +29,46 @@ def small_raw_file(tmp_path):
return master_path
MODULE_ROWS = 2
MODULE_COLS = 3
def write_module_layout(directory, layout_rows, layout_cols, frames=2,
module_rows=MODULE_ROWS, module_cols=MODULE_COLS,
detector_type="Jungfrau"):
"""Write a raw file series with modules in the given layout.
Modules are numbered as DetectorGeometry orders them, down each column
first, and every module holds distinct values so misplacement is visible.
Returns the master path and the expected image of frame 0.
"""
master_path = directory / "run_master_0.json"
master_path.write_text(json.dumps({
"Version": 7.2,
"Detector Type": detector_type,
"Timing Mode": "auto",
"Geometry": {"x": layout_cols, "y": layout_rows},
"Image Size in bytes": module_rows * module_cols * 2,
"Pixels": {"x": module_cols, "y": module_rows},
"Max Frames Per File": 1,
"Total Frames": frames,
"Frames in File": frames,
"Frame Padding": 1,
"Frame Discard Policy": "nodiscard",
}))
expected = np.zeros((module_rows * layout_rows, module_cols * layout_cols), dtype=np.uint16)
for module in range(layout_rows * layout_cols):
module_row, module_col = module % layout_rows, module // layout_rows
pixels = ((np.arange(module_rows * module_cols, dtype=np.uint32) + module * 10 + 1) % 65536).astype(np.uint16)
expected[module_row * module_rows:(module_row + 1) * module_rows,
module_col * module_cols:(module_col + 1) * module_cols] = pixels.reshape(module_rows, module_cols)
for index in range(frames):
data_path = directory / f"run_d{module}_f{index}_0.raw"
# 112 byte detector header followed by the module pixels
data_path.write_bytes(bytes(112) + pixels.tobytes())
return master_path, expected
@pytest.mark.parametrize("reader_type", [RawFile, File])
@pytest.mark.parametrize("legacy_master", [False, True])
@pytest.mark.parametrize("padding, policy, supported", [
@@ -453,3 +494,109 @@ def test_read_eiger_udp_port_disabled(test_data_path):
assert len(rois) == 2
assert rois[0] == ROI(0, 512, 0, 512)
assert rois[1] == ROI(1024, 1536, 0, 512)
@pytest.mark.parametrize("layout", [(1, 1), (2, 1), (1, 2), (2, 2), (3, 2)])
@pytest.mark.parametrize("reader_type", [RawFile, File])
def test_raw_file_assembles_module_layout(tmp_path, layout, reader_type):
layout_rows, layout_cols = layout
master_path, expected = write_module_layout(tmp_path, layout_rows, layout_cols)
with reader_type(master_path) as reader:
# RawFile binds rows/cols as methods (ROI overloads), File as properties
rows = reader.rows() if reader_type is RawFile else reader.rows
cols = reader.cols() if reader_type is RawFile else reader.cols
assert (rows, cols) == expected.shape
result = reader.read_frame()
image = result[1] if reader_type is RawFile else result
assert image.shape == expected.shape
assert (image == expected).all()
reader.seek(0)
result = reader.read_n(2)
images = result[1] if reader_type is RawFile else result
assert images.shape == (2,) + expected.shape
assert (images == expected).all()
@pytest.mark.parametrize("layout", [(1, 1), (2, 1), (1, 2)])
def test_raw_file_gap_pixels(tmp_path, layout):
layout_rows, layout_cols = layout
master_path, expected = write_module_layout(tmp_path, layout_rows, layout_cols,
module_rows=512, module_cols=1024)
gaps = GapPixels(fill_value=5)
with RawFile(master_path, gap_pixels=gaps) as reader:
assert reader.gap_pixels == gaps
assert (reader.rows(), reader.cols()) == (512 * layout_rows + 2 * (layout_rows * 2 - 1),
1024 * layout_cols + 2 * (layout_cols * 4 - 1))
header, image = reader.read_frame()
assert header.size == layout_rows * layout_cols
assert (image == insert_gap_pixels(expected, DetectorType.Jungfrau, gaps)).all()
reader.seek(0)
_, images = reader.read_n(2)
assert images.shape == (2, reader.rows(), reader.cols())
assert (images[1] == image).all()
with RawFile(master_path) as reader:
assert reader.gap_pixels is None
assert (reader.rows(), reader.cols()) == expected.shape
with RawFile(master_path, gap_pixels=True) as reader:
assert reader.gap_pixels == GapPixels()
_, image = reader.read_frame()
assert (image == insert_gap_pixels(expected, DetectorType.Jungfrau)).all()
with File(master_path, gap_pixels=True) as reader:
assert (reader.rows, reader.cols) == insert_gap_pixels(expected, DetectorType.Jungfrau).shape
assert (reader.read_frame() == insert_gap_pixels(expected, DetectorType.Jungfrau)).all()
module_gaps = GapPixels(module_gaps=ModuleGaps(8, 36))
with File(master_path, gap_pixels=module_gaps) as reader:
assert (reader.read_frame() == insert_gap_pixels(expected, DetectorType.Jungfrau, module_gaps)).all()
def test_raw_file_gap_pixels_split_counts(tmp_path):
master_path, expected = write_module_layout(tmp_path, 1, 1, module_rows=512, module_cols=1024)
gaps = GapPixels(split_counts=True, seed=11)
with RawFile(master_path, gap_pixels=gaps) as reader:
_, image = reader.read_frame()
assert (image == insert_gap_pixels(expected, DetectorType.Jungfrau, gaps)).all()
assert image.sum() == expected.sum()
def test_raw_file_gap_pixels_unsupported_detector(tmp_path):
master_path, _ = write_module_layout(tmp_path, 1, 1, detector_type="Moench")
RawFile(master_path)
with pytest.raises(ValueError):
RawFile(master_path, gap_pixels=True)
with pytest.raises(ValueError):
File(master_path, gap_pixels=True)
@pytest.mark.withdata
def test_raw_file_gap_pixels_recorded_jungfrau(test_data_path):
path = test_data_path / "raw/jungfrau/jungfrau_single_master_0.json"
gaps = GapPixels(fill_value=3)
with RawFile(path, gap_pixels=gaps) as gapped, RawFile(path) as plain:
assert (gapped.rows(), gapped.cols()) == (514, 1030)
_, image = gapped.read_frame()
_, reference = plain.read_frame()
assert image.shape == (514, 1030)
assert (image == insert_gap_pixels(reference, DetectorType.Jungfrau, gaps)).all()
@pytest.mark.withdata
def test_raw_file_gap_pixels_recorded_rois(test_data_path):
path = test_data_path / "raw/ROITestData/MultipleROIs/run_master_0.json"
with RawFile(path, gap_pixels=True) as gapped, RawFile(path) as plain:
assert gapped.num_rois == 2
_, images = gapped.read_rois()
_, references = plain.read_rois()
rois = plain.master.rois
assert images[0].shape == (303, 103)
assert images[1].shape == (101, 103)
for image, reference, roi in zip(images, references, rois):
assert (image == insert_gap_pixels(reference, DetectorType.Jungfrau, roi=roi)).all()
gapped.seek(0)
_, frames = gapped.read_n_with_roi(2, roi_index=1)
assert frames.shape == (2, 101, 103)
+5 -3
View File
@@ -23,10 +23,12 @@ File::File(const std::filesystem::path &fname, const std::string &mode,
// Assuming we are pointing at a master file?
// TODO! How do we read raw files directly?
if (fname.extension() == ".raw" || fname.extension() == ".json") {
// file_impl = new RawFile(fname, mode, cfg);
file_impl = std::make_unique<RawFile>(fname, mode);
file_impl = std::make_unique<RawFile>(fname, mode, cfg.gap_pixels);
} else if (cfg.gap_pixels) {
throw std::invalid_argument(fmt::format(
"Gap pixels are only supported for raw files, cannot open '{}'",
fname.string()));
} else if (fname.extension() == ".npy") {
// file_impl = new NumpyFile(fname, mode, cfg);
file_impl = std::make_unique<NumpyFile>(fname, mode, cfg);
} else if (fname.extension() == ".dat") {
file_impl = std::make_unique<JungfrauDataFile>(fname);
+14 -4
View File
@@ -24,11 +24,12 @@ size_t checked_product(size_t lhs, size_t rhs) {
MultiThreadedFileReader::MultiThreadedFileReader(
std::filesystem::path fname, size_t n_threads, size_t chunk_size,
std::optional<size_t> total_frames)
std::optional<size_t> total_frames, std::optional<GapPixels> gap_pixels)
: m_fname(std::move(fname)), m_n_threads(n_threads),
m_chunk_size(chunk_size), m_total_frames(0), m_source_total_frames(0),
m_rows(0), m_cols(0), m_bitdepth(0), m_dtype(Dtype::NONE),
m_bytes_per_frame(0), m_total_bytes(0), m_current_frame(0) {
m_bytes_per_frame(0), m_total_bytes(0), m_current_frame(0),
m_gap_pixels(std::move(gap_pixels)) {
if (m_n_threads == 0) {
throw std::invalid_argument(
"MultiThreadedFileReader requires at least one thread");
@@ -38,7 +39,7 @@ MultiThreadedFileReader::MultiThreadedFileReader(
"MultiThreadedFileReader chunk size must be greater than zero");
}
File file(m_fname);
File file(m_fname, "r", worker_config(0));
m_source_total_frames = file.total_frames();
m_total_frames = total_frames.value_or(m_source_total_frames);
if (m_total_frames > m_source_total_frames) {
@@ -56,10 +57,19 @@ MultiThreadedFileReader::MultiThreadedFileReader(
m_files.reserve(m_n_threads);
m_files.push_back(std::move(file));
for (size_t i = 1; i < m_n_threads; ++i) {
m_files.emplace_back(m_fname);
m_files.emplace_back(m_fname, "r", worker_config(i));
}
}
FileConfig MultiThreadedFileReader::worker_config(size_t worker_index) const {
FileConfig config;
config.gap_pixels = m_gap_pixels;
if (config.gap_pixels && config.gap_pixels->seed) {
*config.gap_pixels->seed += worker_index;
}
return config;
}
size_t MultiThreadedFileReader::remaining_frames() const noexcept {
return m_total_frames - m_current_frame;
}
+166 -142
View File
@@ -9,6 +9,7 @@
#include "aare/logger.hpp"
#include <algorithm>
#include <cstring>
#include <fmt/format.h>
#include <nlohmann/json.hpp>
@@ -23,7 +24,8 @@ std::runtime_error RawFile::frame_error(size_t frame_index,
frame_index, m_master.master_fname().string(), message));
}
RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode)
RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode,
std::optional<GapPixels> gap_pixels)
: m_master(fname) {
m_mode = mode;
@@ -49,11 +51,37 @@ RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode)
std::optional<size_t> min_frames;
size_t max_frames = 0;
size_t max_part_bytes = 0;
for (const auto &subfiles : m_subfiles) {
for (const auto &subfile : subfiles) {
const auto count = subfile->frames_in_file();
min_frames = min_frames ? std::min(*min_frames, count) : count;
max_frames = std::max(max_frames, count);
max_part_bytes =
std::max(max_part_bytes, subfile->bytes_per_frame());
}
}
m_part_buffer.resize(max_part_bytes);
m_roi_rects = m_master.rois();
if (m_roi_rects.size() != num_rois) {
throw std::runtime_error(
LOCATION + "ROI count differs from ROI geometry count");
}
m_roi_shapes.reserve(num_rois);
if (gap_pixels) {
m_gap_pixels = gap_pixels;
m_gap_layout = detail::gap_layout(
m_master.detector_type(), m_master.quad() == 1, *gap_pixels);
m_generator = detail::make_generator(gap_pixels->seed);
for (const auto &rect : m_roi_rects) {
m_roi_shapes.push_back(
detail::gapped_shape(rect, m_gap_layout));
}
} else {
for (const auto &roi : m_master.roi_geometries()) {
m_roi_shapes.push_back({static_cast<ssize_t>(roi.pixels_y()),
static_cast<ssize_t>(roi.pixels_x())});
}
}
if (!min_frames) {
@@ -197,9 +225,7 @@ size_t RawFile::bytes_per_frame() {
}
size_t RawFile::bytes_per_frame(const size_t roi_index) {
return m_master.roi_geometries().at(roi_index).pixels_x() *
m_master.roi_geometries().at(roi_index).pixels_y() *
m_master.bitdepth() / bits_per_byte;
return pixels_per_frame(roi_index) * m_master.bitdepth() / bits_per_byte;
}
size_t RawFile::pixels_per_frame() {
@@ -212,8 +238,7 @@ size_t RawFile::pixels_per_frame() {
}
size_t RawFile::pixels_per_frame(const size_t roi_index) {
return m_master.roi_geometries().at(roi_index).pixels_x() *
m_master.roi_geometries().at(roi_index).pixels_y();
return rows(roi_index) * cols(roi_index);
}
DetectorType RawFile::detector_type() const { return m_master.detector_type(); }
@@ -244,7 +269,7 @@ size_t RawFile::rows() const {
return rows(0);
}
size_t RawFile::rows(const size_t roi_index) const {
return m_master.roi_geometries().at(roi_index).pixels_y();
return static_cast<size_t>(m_roi_shapes.at(roi_index)[0]);
}
size_t RawFile::cols() const {
if (m_master.roi_geometries().size() > 1) {
@@ -255,7 +280,7 @@ size_t RawFile::cols() const {
return cols(0);
}
size_t RawFile::cols(const size_t roi_index) const {
return m_master.roi_geometries().at(roi_index).pixels_x();
return static_cast<size_t>(m_roi_shapes.at(roi_index)[1]);
}
size_t RawFile::bitdepth() const { return m_master.bitdepth(); }
@@ -328,8 +353,8 @@ DetectorHeader RawFile::read_header(const std::filesystem::path &fname) {
RawMasterFile RawFile::master() const { return m_master; }
Frame RawFile::get_frame(size_t frame_index, const size_t roi_index) {
auto f = Frame(m_master.roi_geometries().at(roi_index).pixels_y(),
m_master.roi_geometries().at(roi_index).pixels_x(),
auto f = Frame(static_cast<uint32_t>(rows(roi_index)),
static_cast<uint32_t>(cols(roi_index)),
Dtype::from_bitdepth(m_master.bitdepth()));
std::byte *frame_buffer = f.data();
get_frame_into(frame_index, frame_buffer, roi_index);
@@ -338,6 +363,47 @@ Frame RawFile::get_frame(size_t frame_index, const size_t roi_index) {
size_t RawFile::bytes_per_pixel() const { return m_master.bitdepth() / 8; }
std::vector<size_t> 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<size_t> frame_indices(n_parts, frame_index);
if (n_parts == 1) {
return frame_indices;
}
std::vector<size_t> 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<size_t>(std::distance(
frame_numbers.begin(),
std::min_element(frame_numbers.begin(), frame_numbers.end())));
frame_indices[min_frame_idx]++;
if (frame_indices[min_frame_idx] >= total_frames()) {
throw frame_error(
frame_index,
LOCATION +
fmt::format(
"Frame index {} out of range while synchronizing "
"module {} for ROI {}; last data file '{}'",
frame_indices[min_frame_idx], min_frame_idx, roi_index,
m_subfiles[roi_index][min_frame_idx]
->current_path()
.string()));
}
frame_numbers[min_frame_idx] =
m_subfiles[roi_index][min_frame_idx]->frame_number(
frame_indices[min_frame_idx]);
}
return frame_indices;
}
void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer,
const size_t roi_index, DetectorHeader *header) {
LOG(logDEBUG) << "RawFile::get_frame_into(" << frame_index << ")";
@@ -349,148 +415,106 @@ void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer,
"frames (indices are zero-based)",
frame_index, total_frames()));
}
std::vector<size_t> frame_numbers(
m_master.roi_geometries().at(roi_index).num_modules_in_roi());
std::vector<size_t> frame_indices(
m_master.roi_geometries().at(roi_index).num_modules_in_roi(),
frame_index);
// sync the frame numbers
const auto &roi = m_master.roi_geometries().at(roi_index);
const auto frame_indices =
synchronized_frame_indices(frame_index, roi_index);
if (m_master.roi_geometries().at(roi_index).num_modules_in_roi() !=
1) { // if we have more than one module
for (size_t part_idx = 0;
part_idx !=
m_master.roi_geometries().at(roi_index).num_modules_in_roi();
++part_idx) {
frame_numbers[part_idx] =
m_subfiles[roi_index][part_idx]->frame_number(frame_index);
}
// 1. if frame number vector is the same break
while (!all_equal(frame_numbers)) {
// 2. find the index of the minimum frame number,
auto min_frame_idx = std::distance(
frame_numbers.begin(),
std::min_element(frame_numbers.begin(), frame_numbers.end()));
// 3. increase its index and update its respective frame
// number
frame_indices[min_frame_idx]++;
// 4. if we can't increase its index => throw error
if (frame_indices[min_frame_idx] >= total_frames()) {
throw frame_error(
frame_index,
LOCATION +
fmt::format(
"Frame index {} out of range while synchronizing "
"module {} for ROI {}; last data file '{}'",
frame_indices[min_frame_idx], min_frame_idx,
roi_index,
m_subfiles[roi_index][min_frame_idx]
->current_path()
.string()));
}
frame_numbers[min_frame_idx] =
m_subfiles[roi_index][min_frame_idx]->frame_number(
frame_indices[min_frame_idx]);
if (m_gap_pixels) {
switch (m_master.bitdepth()) {
case 8:
assemble_with_gaps<uint8_t>(frame_buffer, roi_index, frame_indices,
header);
break;
case 16:
assemble_with_gaps<uint16_t>(frame_buffer, roi_index, frame_indices,
header);
break;
case 32:
assemble_with_gaps<uint32_t>(frame_buffer, roi_index, frame_indices,
header);
break;
default:
throw frame_error(frame_index,
LOCATION +
fmt::format("Gap pixels are not supported "
"for a bit depth of {}",
m_master.bitdepth()));
}
return;
}
if (m_master.detector_layout().col == 1) {
// get the part from each subfile and copy it to the frame
for (size_t part_idx = 0;
part_idx !=
m_master.roi_geometries().at(roi_index).num_modules_in_roi();
++part_idx) {
auto corrected_idx = frame_indices[part_idx];
const size_t bytes_per_pixel = m_master.bitdepth() / bits_per_byte;
const size_t roi_width = roi.pixels_x();
// This is where we start writing
auto offset =
(m_master.geometry()
.get_module_geometries(
m_master.roi_geometries()
.at(roi_index)
.module_indices_in_roi(part_idx))
.origin_y *
m_master.roi_geometries().at(roi_index).pixels_x() +
m_master.geometry()
.get_module_geometries(
m_master.roi_geometries()
.at(roi_index)
.module_indices_in_roi(part_idx))
.origin_x) *
m_master.bitdepth() / 8;
for (size_t part_idx = 0; part_idx != roi.num_modules_in_roi();
++part_idx) {
// origin and size of the part relative to the ROI
const auto &pos = m_master.geometry().get_module_geometries(
roi.module_indices_in_roi(part_idx));
const auto width = static_cast<size_t>(pos.width);
const auto height = static_cast<size_t>(pos.height);
auto &subfile = *m_subfiles[roi_index][part_idx];
subfile.seek(frame_indices[part_idx]);
if (m_master.geometry()
.get_module_geometries(m_master.roi_geometries()
.at(roi_index)
.module_indices_in_roi(part_idx))
.origin_x != 0)
throw frame_error(
frame_index,
LOCATION +
" Implementation error. x pos not 0."); // TODO:
// origin
// can still
// change if
// roi
// changes
// TODO! What if the files don't match?
m_subfiles[roi_index][part_idx]->seek(corrected_idx);
m_subfiles[roi_index][part_idx]->read_into(frame_buffer + offset,
header);
if (header)
++header;
}
std::byte *dest =
frame_buffer + (static_cast<size_t>(pos.origin_y) * roi_width +
static_cast<size_t>(pos.origin_x)) *
bytes_per_pixel;
} else {
// TODO! should we read row by row?
// create a buffer large enough to hold a full module
auto bytes_per_part =
m_master.pixels_y() * m_master.pixels_x() * m_master.bitdepth() /
8; // TODO! replace with image_size_in_bytes // TODO
// shouldnt it only be the module size? - check
auto *part_buffer = new std::byte[bytes_per_part];
// TODO! if we have many submodules we should reorder them on
// the module level
for (size_t part_idx = 0;
part_idx !=
m_master.roi_geometries().at(roi_index).num_modules_in_roi();
++part_idx) {
auto pos = m_master.geometry().get_module_geometries(
m_master.roi_geometries().at(roi_index).module_indices_in_roi(
part_idx));
auto corrected_idx = frame_indices[part_idx];
m_subfiles[roi_index][part_idx]->seek(corrected_idx);
m_subfiles[roi_index][part_idx]->read_into(part_buffer, header);
if (header)
++header;
for (size_t cur_row = 0; cur_row < static_cast<size_t>(pos.height);
cur_row++) {
auto irow = (pos.origin_y + cur_row);
auto icol = pos.origin_x;
auto dest =
(irow * m_master.roi_geometries().at(roi_index).pixels_x() +
icol);
dest = dest * m_master.bitdepth() / 8;
memcpy(frame_buffer + dest,
part_buffer +
cur_row * pos.width * m_master.bitdepth() / 8,
pos.width * m_master.bitdepth() / 8);
if (width == roi_width) {
// the part spans full rows of the ROI and is contiguous in the
// assembled frame
subfile.read_into(dest, header);
} else {
subfile.read_into(m_part_buffer.data(), header);
const size_t row_bytes = width * bytes_per_pixel;
const size_t roi_row_bytes = roi_width * bytes_per_pixel;
for (size_t row = 0; row < height; ++row) {
std::memcpy(dest + row * roi_row_bytes,
m_part_buffer.data() + row * row_bytes, row_bytes);
}
}
delete[] part_buffer;
if (header) {
++header;
}
}
}
template <typename T>
void RawFile::assemble_with_gaps(std::byte *frame_buffer, size_t roi_index,
const std::vector<size_t> &frame_indices,
DetectorHeader *header) {
const auto &roi = m_master.roi_geometries().at(roi_index);
const ROI &rect = m_roi_rects.at(roi_index);
NDView<T, 2> destination(reinterpret_cast<T *>(frame_buffer),
m_roi_shapes.at(roi_index));
detail::fill_gap_pixels(destination, rect, m_gap_layout,
detail::fill_as<T>(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<const T, 2> source(
reinterpret_cast<const T *>(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);
}
}
+335 -2
View File
@@ -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 <catch2/catch_test_macros.hpp>
@@ -222,6 +223,338 @@ TEST_CASE("RawFile frame count comes from the complete subfile series",
REQUIRE_NOTHROW(reader.read_frame(1));
}
TEST_CASE("RawFile assembles every module layout", "[RawFile][geometry]") {
const auto [layout_rows, layout_cols] = GENERATE(table<uint32_t, uint32_t>({
{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<uint16_t, 2> image) {
REQUIRE(image.shape(0) == static_cast<ssize_t>(rows));
REQUIRE(image.shape(1) == static_cast<ssize_t>(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<uint16_t>());
}
SECTION("read_into with headers") {
NDArray<uint16_t, 2> image(
{static_cast<ssize_t>(rows), static_cast<ssize_t>(cols)});
std::vector<DetectorHeader> headers(files.n_modules());
file.seek(1);
file.read_into(reinterpret_cast<std::byte *>(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<uint16_t>());
}
}
SECTION("generic File") {
File reader(files.master_path());
auto frame = reader.read_frame(1);
check(frame.view<uint16_t>());
}
}
namespace {
/// Compare a frame read with gap pixels against the free function applied
/// to the same frame read without gaps.
template <typename T>
bool matches_gap_oracle(Frame &gapped, Frame &plain, const ROI &roi,
DetectorType detector, bool quad,
const GapPixels &gaps) {
auto expected = insert_gap_pixels(NDView<const T, 2>(plain.view<T>()), roi,
detector, quad, gaps);
auto actual = gapped.view<T>();
if (expected.shape() != actual.shape()) {
UNSCOPED_INFO("shape " << actual.shape(0) << "x" << actual.shape(1)
<< " expected " << expected.shape(0) << "x"
<< expected.shape(1));
return false;
}
for (ssize_t row = 0; row < expected.shape(0); ++row) {
for (ssize_t col = 0; col < expected.shape(1); ++col) {
if (expected(row, col) != actual(row, col)) {
UNSCOPED_INFO("pixel (" << row << ", " << col << ") "
<< actual(row, col) << " expected "
<< expected(row, col));
return false;
}
}
}
return true;
}
bool matches_gap_oracle(Frame &gapped, Frame &plain, const ROI &roi,
DetectorType detector, bool quad,
const GapPixels &gaps) {
switch (plain.bitdepth()) {
case 8:
return matches_gap_oracle<uint8_t>(gapped, plain, roi, detector, quad,
gaps);
case 16:
return matches_gap_oracle<uint16_t>(gapped, plain, roi, detector, quad,
gaps);
case 32:
return matches_gap_oracle<uint32_t>(gapped, plain, roi, detector, quad,
gaps);
default:
return false;
}
}
} // namespace
TEST_CASE("RawFile inserts gap pixels into synthetic Jungfrau modules",
"[RawFile][GapPixels]") {
const auto [layout_rows, layout_cols] = GENERATE(table<uint32_t, uint32_t>({
{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<ssize_t>(files.cols()), 0,
static_cast<ssize_t>(files.rows())};
GapPixels gaps;
SECTION("defaults") {}
SECTION("module gaps and fill value") {
gaps.module_gaps = ModuleGaps{8, 36};
gaps.fill_value = 7;
}
const auto layout = detail::gap_layout(DetectorType::Jungfrau, false, gaps);
const auto shape = detail::gapped_shape(detector, layout);
const auto fill = static_cast<uint16_t>(gaps.fill_value);
RawFile file(files.master_path(), "r", gaps);
REQUIRE(file.gap_pixels() == gaps);
REQUIRE(file.rows() == static_cast<size_t>(shape[0]));
REQUIRE(file.cols() == static_cast<size_t>(shape[1]));
REQUIRE(file.pixels_per_frame() ==
static_cast<size_t>(shape[0] * shape[1]));
REQUIRE(file.bytes_per_frame() == file.pixels_per_frame() * 2);
std::vector<DetectorHeader> headers(files.n_modules());
NDArray<uint16_t, 2> image(shape);
file.read_into(reinterpret_cast<std::byte *>(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<ssize_t>(row)),
layout.x.gapped(static_cast<ssize_t>(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<size_t>(shape[0]));
REQUIRE(frame.cols() == static_cast<size_t>(shape[1]));
REQUIRE(frame.view<uint16_t>() == image.view());
}
TEST_CASE("RawFile splits counts like the free function",
"[RawFile][GapPixels]") {
TemporaryRawFiles files(xy{1, 2}, xy{512, 1024});
GapPixels gaps;
gaps.split_counts = true;
gaps.seed = 7;
RawFile gapped(files.master_path(), "r", gaps);
RawFile plain(files.master_path());
auto gapped_frame = gapped.read_frame();
auto plain_frame = plain.read_frame();
const ROI detector{0, 2048, 0, 512};
REQUIRE(matches_gap_oracle(gapped_frame, plain_frame, detector,
DetectorType::Jungfrau, false, gaps));
// the next frame continues the generator, so the oracle for frame 1 is
// the free function applied after the same number of draws
auto frames = gapped.read_n(1);
auto second_plain = plain.read_frame();
auto expected = insert_gap_pixels(
NDView<const uint16_t, 2>(second_plain.view<uint16_t>()), detector,
DetectorType::Jungfrau, false, gaps);
double gapped_total = 0;
double expected_total = 0;
auto view = frames[0].view<uint16_t>();
for (ssize_t row = 0; row < expected.shape(0); ++row) {
for (ssize_t col = 0; col < expected.shape(1); ++col) {
gapped_total += view(row, col);
expected_total += expected(row, col);
}
}
REQUIRE(gapped_total == expected_total);
}
TEST_CASE("RawFile rejects gap pixels for unsupported detectors",
"[RawFile][GapPixels]") {
TemporaryRawFiles files;
auto metadata = nlohmann::json::parse(std::ifstream(files.master_path()));
metadata["Detector Type"] = "Moench";
std::ofstream(files.master_path()) << metadata;
REQUIRE_NOTHROW(RawFile(files.master_path()));
REQUIRE_THROWS_AS(RawFile(files.master_path(), "r", GapPixels{}),
std::invalid_argument);
}
TEST_CASE("Gap pixels on recorded Jungfrau and Eiger files",
"[.with-data][RawFile][GapPixels]") {
struct Case {
const char *path;
std::array<ssize_t, 2> shape;
std::optional<ModuleGaps> module_gaps;
};
const auto test = GENERATE(values<Case>({
{"raw/jungfrau/jungfrau_single_master_0.json", {514, 1030}, {}},
{"raw/jungfrau/jungfrau_double_master_0.json", {514, 1030}, {}},
{"raw/eiger/Lab6_20500eV_2deg_20240629_master_7.json", {514, 1030}, {}},
{"raw/eiger/eiger_500k_16bit_master_0.json", {514, 1030}, {}},
{"raw/eiger_quad_data/"
"W13_vrpreampscan_m21C_300V_800eV_vthre2000_master_0.json",
{514, 514},
{}},
{"raw/virtual/920/jungfrau/two_modules_master_0.json",
{1030, 1030},
{}},
{"raw/virtual/920/jungfrau/two_modules_master_0.json",
{1064, 1030},
ModuleGaps{8, 36}},
{"raw/ROITestData/MultipleChipROI/run_master_2.json", {303, 103}, {}},
{"raw/ROITestData/MultipleChipROI/run_master_2.json",
{337, 103},
ModuleGaps{8, 36}},
{"raw/virtual/920/jungfrau/"
"single_module_roi_512_767_0_255_master_0.json",
{256, 256},
{}},
}));
CAPTURE(test.path);
const auto fpath = test_data_path() / test.path;
REQUIRE(std::filesystem::exists(fpath));
GapPixels gaps;
gaps.module_gaps = test.module_gaps;
gaps.fill_value = 3;
RawFile gapped(fpath, "r", gaps);
RawFile plain(fpath);
REQUIRE(gapped.rows() == static_cast<size_t>(test.shape[0]));
REQUIRE(gapped.cols() == static_cast<size_t>(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<size_t>(test.shape[0]));
REQUIRE(gapped_frame.cols() == static_cast<size_t>(test.shape[1]));
REQUIRE(matches_gap_oracle(gapped_frame, plain_frame, plain.master().roi(),
plain.detector_type(),
plain.master().quad() == 1, gaps));
}
TEST_CASE("Split counts on a recorded Jungfrau file",
"[.with-data][RawFile][GapPixels]") {
const auto fpath =
test_data_path() / "raw/jungfrau/jungfrau_single_master_0.json";
GapPixels gaps;
gaps.split_counts = true;
gaps.seed = 2024;
RawFile gapped(fpath, "r", gaps);
RawFile plain(fpath);
auto gapped_frame = gapped.read_frame();
auto plain_frame = plain.read_frame();
REQUIRE(matches_gap_oracle(gapped_frame, plain_frame, plain.master().roi(),
DetectorType::Jungfrau, false, gaps));
}
TEST_CASE("Gap pixels on a file with two ROIs",
"[.with-data][RawFile][GapPixels]") {
const auto fpath =
test_data_path() / "raw/ROITestData/MultipleROIs/run_master_0.json";
GapPixels gaps;
RawFile gapped(fpath, "r", gaps);
RawFile plain(fpath);
REQUIRE(gapped.num_rois() == 2);
// x 200..300 holds the chip boundary at 256; y 100..400 holds the one at
// 256 while y 600..700 holds none
REQUIRE(gapped.rows(0) == 303);
REQUIRE(gapped.cols(0) == 103);
REQUIRE(gapped.rows(1) == 101);
REQUIRE(gapped.cols(1) == 103);
auto gapped_frames = gapped.read_rois();
auto plain_frames = plain.read_rois();
const auto rois = plain.master().rois();
REQUIRE(gapped_frames.size() == 2);
for (size_t index = 0; index < 2; ++index) {
CAPTURE(index);
REQUIRE(matches_gap_oracle(gapped_frames[index], plain_frames[index],
rois[index], DetectorType::Jungfrau, false,
gaps));
}
}
TEST_CASE("Read number of frames from a jungfrau raw file",
"[.with-data][RawFile]") {
@@ -535,7 +868,7 @@ TEST_CASE("Open multi module file with ROI",
}
TEST_CASE("Open multi module file with ROI spanning over multiple modules",
"[.width-data][read_rois][RawFile]") {
"[.with-data][read_rois][RawFile]") {
auto fpath =
test_data_path() / "raw/ROITestData/MultipleChipROI/run_master_2.json";
@@ -557,7 +890,7 @@ TEST_CASE("Open multi module file with ROI spanning over multiple modules",
}
TEST_CASE("Open multi module file with two ROIs",
"[.width-data][read_rois][RawFile]") {
"[.with-data][read_rois][RawFile]") {
auto fpath =
test_data_path() / "raw/ROITestData/MultipleROIs/run_master_0.json";
+66 -11
View File
@@ -8,44 +8,99 @@
#include <fstream>
#include <nlohmann/json.hpp>
#include <string>
#include <vector>
/**
* @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<uint16_t>((module * 10 + pixel + 1) % 65536);
}
/// Modules stacked in one column
explicit TemporaryRawFiles(size_t modules = 1)
: TemporaryRawFiles(aare::xy{static_cast<uint32_t>(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<uint16_t> 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<uint16_t>(module % layout.row);
header.column = static_cast<uint16_t>(module / layout.row);
output.write(reinterpret_cast<const char *>(&header),
sizeof(header));
output.write(reinterpret_cast<const char *>(pixels),
sizeof(pixels));
output.write(reinterpret_cast<const char *>(pixels.data()),
static_cast<std::streamsize>(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); }