Jungfraujoch/tests/XFELReaderTest.cpp

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <catch2/catch_all.hpp>
#include <algorithm>
#include <climits>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <hdf5.h>
#include <string>
#include <vector>

#include "../reader/JFJochXFELReader.h"
#include "../reader/JFJochXTCReader.h"

// Synthetic SACLA Cheetah and Cheetah CXI files, written with the HDF5 C API the way those pipelines
// write them: whole frames as single deflated (and byte-shuffled) chunks, so the direct-chunk path is
// what is read.
namespace {

void WriteFrames(hid_t loc, const std::string &name, hid_t type, const void *data, std::vector<hsize_t> dims,
                 bool shuffle = true) {
    const hid_t space = H5Screate_simple(static_cast<int>(dims.size()), dims.data(), nullptr);
    std::vector<hsize_t> chunk = dims;
    if (chunk.size() == 3)
        chunk[0] = 1;
    const hid_t dcpl = H5Pcreate(H5P_DATASET_CREATE);
    H5Pset_chunk(dcpl, static_cast<int>(chunk.size()), chunk.data());
    if (shuffle)
        H5Pset_shuffle(dcpl);
    H5Pset_deflate(dcpl, 4);
    const hid_t ds = H5Dcreate2(loc, name.c_str(), type, space, H5P_DEFAULT, dcpl, H5P_DEFAULT);
    H5Dwrite(ds, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
    H5Dclose(ds);
    H5Pclose(dcpl);
    H5Sclose(space);
}

void WriteVector(hid_t loc, const std::string &name, const std::vector<double> &v) {
    const hsize_t n = v.size();
    const hid_t space = H5Screate_simple(1, &n, nullptr);
    const hid_t ds = H5Dcreate2(loc, name.c_str(), H5T_NATIVE_DOUBLE, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(ds, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, v.data());
    H5Dclose(ds);
    H5Sclose(space);
}

void WriteScalar(hid_t loc, const std::string &name, double v) {
    const hid_t space = H5Screate(H5S_SCALAR);
    const hid_t ds = H5Dcreate2(loc, name.c_str(), H5T_NATIVE_DOUBLE, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    H5Dwrite(ds, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &v);
    H5Dclose(ds);
    H5Sclose(space);
}

// A pixel of the flat image, as the reader hands it out: a 16-bit signed frame without a gain in its own
// type, anything else as int32.
int32_t PixelAt(JFJochXFELReader &reader, int64_t image, size_t col, size_t row,
                CompressedImageMode mode = CompressedImageMode::Int32) {
    auto raw = reader.GetRawImage(image);
    REQUIRE(raw);
    REQUIRE(raw->image.GetMode() == mode);
    const size_t i = row * raw->image.GetWidth() + col;
    if (mode == CompressedImageMode::Int16)
        return reinterpret_cast<const int16_t *>(raw->image_buffer.data())[i];
    return reinterpret_cast<const int32_t *>(raw->image_buffer.data())[i];
}

}  // namespace

TEST_CASE("JFJochXFELReader_SACLA", "[portable]") {
    // Two sensors of 8 rows x 4 columns stacked in a 16 x 4 frame, 50 um pixels. The first lies square
    // with its first corner on the beam; the second is turned by 90 degrees with its corner 20 pixels
    // along x (posx = -1000 um), so its fast axis runs along y and its slow axis along -x.
    const std::string name = "xfel_reader_sacla.h5";
    {
        const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
        const hid_t meta = H5Gcreate2(f, "metadata", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        WriteVector(meta, "posx_in_um", {0, -1000});
        WriteVector(meta, "posy_in_um", {0, 0});
        WriteVector(meta, "posz_in_um", {0, 0});
        WriteVector(meta, "angle_in_rad", {0, 90});   // degrees, whatever the name says
        WriteVector(meta, "pixelsizex_in_um", {50, 50});
        WriteVector(meta, "pixelsizey_in_um", {50, 50});
        WriteScalar(meta, "distance_in_mm", 60.0);
        std::vector<uint16_t> pixelmask(16 * 4, 0);
        pixelmask[1] = 1;   // first sensor, row 0, column 1
        WriteFrames(meta, "pixelmask", H5T_NATIVE_UINT16, pixelmask.data(), {16, 4});
        H5Gclose(meta);
        for (int shot = 0; shot < 2; shot++) {
            const std::string group = "tag-" + std::to_string(100 + shot);
            const hid_t g = H5Gcreate2(f, group.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            std::vector<int16_t> frame(16 * 4, 10);
            frame[0] = static_cast<int16_t>(20 + shot);    // first sensor (fs 0, ss 0)
            frame[8 * 4] = -5;                             // second sensor (fs 0, ss 0): below zero
            frame[8 * 4 + 1] = 33;                         // second sensor (fs 1, ss 0)
            frame[9 * 4] = 44;                             // second sensor (fs 0, ss 1)
            frame[2 * 4] = INT16_MIN;                      // first sensor (fs 0, ss 2): no measurement
            WriteFrames(g, "data", H5T_NATIVE_INT16, frame.data(), {16, 4});
            WriteScalar(g, "photon_energy_ev", 10000.0 + 10.0 * shot);
            H5Gclose(g);
        }
        H5Fclose(f);
    }

    REQUIRE(JFJochXFELReader::CanReadPanelTable(name));
    CHECK_FALSE(JFJochXFELReader::CanReadCXI(name));
    JFJochXFELReader reader;
    REQUIRE_NOTHROW(reader.ReadPanelTable({name}));
    CHECK(reader.GetNumberOfImages() == 2);
    CHECK_FALSE(reader.CountsInPhotons());
    const auto x = reader.GetDataset()->experiment;
    CHECK(x.GetIncidentEnergy_keV() == Catch::Approx(10.005));
    CHECK(x.GetDetectorDistance_mm() == Catch::Approx(60.0));
    CHECK(x.GetPixelSize_mm() == Catch::Approx(0.05));
    // Two panels: the image is the stored frame, each sensor a module at its place on the detector plane.
    CHECK(x.GetXPixelsNum() == 4);
    CHECK(x.GetYPixelsNum() == 16);
    const auto &modules = x.GetModules();
    REQUIRE(modules.size() == 2);
    auto plane = [&](int64_t panel, double col, double row, double px, double py) {
        const auto [x_, y_] = modules.ImageToPlane(panel, col, row);
        CHECK(x_ == Catch::Approx(px).margin(1e-9));
        CHECK(y_ == Catch::Approx(py).margin(1e-9));
    };
    // The first sensor's first pixel centre is at (0.5, 0.5) from the beam, which the plane puts on a whole
    // column: the beam at 1.5, the pixel at column and row 2.
    CHECK(x.GetBeamX_pxl() == Catch::Approx(1.5));
    CHECK(x.GetBeamY_pxl() == Catch::Approx(1.5));
    plane(0, 0, 0, 2.0, 2.0);
    // Second sensor (rows 8-15): (fs 0, ss 0) centre at x = 20 - 0.5, y = 0.5; fs runs along +y, ss along -x.
    plane(1, 0, 8, 21.0, 2.0);
    plane(1, 1, 8, 21.0, 3.0);
    plane(1, 0, 9, 20.0, 2.0);
    constexpr auto int16 = CompressedImageMode::Int16;
    CHECK(PixelAt(reader, 0, 0, 0, int16) == 20);
    CHECK(PixelAt(reader, 1, 0, 0, int16) == 21);
    CHECK(PixelAt(reader, 0, 0, 8, int16) == -5);
    CHECK(PixelAt(reader, 0, 1, 8, int16) == 33);
    CHECK(PixelAt(reader, 0, 0, 9, int16) == 44);
    // The smallest int16 is no measurement, not a count (handed out as Int16, whose error value it is).
    CHECK(PixelAt(reader, 0, 0, 2, int16) == INT16_MIN);
    // The file's mask holds the first sensor's (fs 1, ss 0).
    const auto &mask = reader.GetDataset()->pixel_mask->GetMask();
    CHECK(mask[1] != 0);
    CHECK(mask[0] == 0);

    reader.ApplyGain(2.0f);
    CHECK(reader.CountsInPhotons());
    CHECK(PixelAt(reader, 0, 0, 0) == 10);
    CHECK(PixelAt(reader, 0, 0, 8) == -3);   // -2.5, rounded away from zero
    CHECK(PixelAt(reader, 0, 1, 8) == 17);   // 16.5, likewise
    reader.Close();
    std::filesystem::remove(name);
}

TEST_CASE("JFJochXFELReader_CXI", "[portable]") {
    // A single-panel Cheetah CXI file: three 6 x 6 frames of signed 32-bit counts, 75 um pixels,
    // 9 keV stated in joules, a mask_shared with one bad pixel.
    const std::string name = "xfel_reader_test.cxi";
    {
        const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
        for (const char *g: {"entry_1", "entry_1/data_1", "entry_1/instrument_1", "entry_1/instrument_1/detector_1",
                             "entry_1/instrument_1/source_1"})
            H5Gclose(H5Gcreate2(f, g, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT));
        std::vector<int32_t> frames(3 * 36, 5);
        frames[0] = -7;               // frame 0, pixel (0, 0): dark-subtracted below zero
        frames[36 + 7] = 70000;       // frame 1, pixel (1, 1)
        WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT32, frames.data(), {3, 6, 6});
        const std::string det = "entry_1/instrument_1/detector_1/";
        WriteVector(f, det + "x_pixel_size", {75e-6, 75e-6, 75e-6});
        WriteVector(f, det + "distance", {0.1, 0.1, 0.1});
        WriteVector(f, "entry_1/instrument_1/source_1/energy", std::vector<double>(3, 9000.0 * 1.602176634e-19));
        std::vector<uint16_t> mask(36, 0);
        mask[35] = 8;   // dead
        WriteFrames(f, det + "mask_shared", H5T_NATIVE_UINT16, mask.data(), {6, 6});
        H5Fclose(f);
    }

    REQUIRE(JFJochXFELReader::CanReadCXI(name));
    CHECK_FALSE(JFJochXFELReader::CanReadPanelTable(name));
    CHECK_FALSE(JFJochXFELReader::CXIRefusal(name).has_value());
    JFJochXFELReader reader;
    REQUIRE_NOTHROW(reader.ReadCXI({name}));
    CHECK(reader.GetNumberOfImages() == 3);
    const auto x = reader.GetDataset()->experiment;
    CHECK(x.GetIncidentEnergy_keV() == Catch::Approx(9.0));
    CHECK(x.GetPixelSize_mm() == Catch::Approx(0.075));
    CHECK(x.GetDetectorDistance_mm() == Catch::Approx(100.0));
    // Corner at (-3, -3): the pixel (fs, ss) lands at column and row 2 + fs, 2 + ss.
    CHECK(PixelAt(reader, 0, 2, 2) == -7);
    CHECK(PixelAt(reader, 1, 3, 3) == 70000);
    CHECK(PixelAt(reader, 2, 4, 4) == 5);
    const auto &m = reader.GetDataset()->pixel_mask->GetMask();
    CHECK(m[7 * x.GetXPixelsNum() + 7] != 0);
    CHECK(m[6 * x.GetXPixelsNum() + 6] == 0);
    reader.Close();
    std::filesystem::remove(name);
}

TEST_CASE("JFJochXFELReader_CXI_TiledRefused", "[portable]") {
    // A Cheetah CXI slab of a tiled detector is not square, and is not read as one flat panel.
    const std::string name = "xfel_reader_tiled.cxi";
    {
        const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
        for (const char *g: {"entry_1", "entry_1/data_1"})
            H5Gclose(H5Gcreate2(f, g, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT));
        std::vector<int16_t> frames(2 * 8 * 4, 1);
        WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT16, frames.data(), {2, 8, 4});
        H5Fclose(f);
    }
    CHECK_FALSE(JFJochXFELReader::CanReadCXI(name));
    const auto refusal = JFJochXFELReader::CXIRefusal(name);
    REQUIRE(refusal.has_value());
    CHECK(refusal->find("tiled detector") != std::string::npos);
    std::filesystem::remove(name);
}

TEST_CASE("JFJochXFELReader_CXI_Int16", "[portable]") {
    // A single-panel Cheetah CXI file as a Rayonix one is written: 16-bit signed frames, deflated without
    // a shuffle. The second frame's chunk is stored with its deflate skipped (filter mask bit 0), which the
    // direct-chunk decoding does not expect: HDF5 reads that one.
    const std::string name = "xfel_reader_int16.cxi";
    {
        const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
        for (const char *g: {"entry_1", "entry_1/data_1", "entry_1/instrument_1", "entry_1/instrument_1/detector_1",
                             "entry_1/instrument_1/source_1"})
            H5Gclose(H5Gcreate2(f, g, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT));
        std::vector<int16_t> frames(2 * 16, 3);
        frames[1] = 32767;          // frame 0, pixel (1, 0)
        frames[4 + 2] = -4;         // frame 0, pixel (2, 1)
        WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT16, frames.data(), {2, 4, 4}, false);
        std::vector<int16_t> raw(16, 9);
        raw[5] = 1234;              // frame 1, pixel (1, 1)
        const hid_t ds = H5Dopen2(f, "entry_1/data_1/data", H5P_DEFAULT);
        const hsize_t offset[3] = {1, 0, 0};
        REQUIRE(H5Dwrite_chunk(ds, H5P_DEFAULT, 1, offset, raw.size() * sizeof(int16_t), raw.data()) >= 0);
        H5Dclose(ds);
        const std::string det = "entry_1/instrument_1/detector_1/";
        WriteVector(f, det + "x_pixel_size", {100e-6, 100e-6});
        WriteVector(f, "entry_1/instrument_1/source_1/energy", std::vector<double>(2, 12000.0 * 1.602176634e-19));
        H5Fclose(f);
    }

    REQUIRE(JFJochXFELReader::CanReadCXI(name));
    JFJochXFELReader reader;
    REQUIRE_NOTHROW(reader.ReadCXI({name}));
    // Corner at (-2, -2): the pixel (fs, ss) lands at column and row 2 + fs, 2 + ss, in its stored type.
    constexpr auto int16 = CompressedImageMode::Int16;
    CHECK(PixelAt(reader, 0, 3, 2, int16) == 32767);
    CHECK(PixelAt(reader, 0, 4, 3, int16) == -4);
    CHECK(PixelAt(reader, 0, 5, 5, int16) == 3);
    CHECK(PixelAt(reader, 0, 1, 1, int16) == 0);   // the margin
    CHECK(PixelAt(reader, 1, 3, 3, int16) == 1234);
    CHECK(PixelAt(reader, 1, 2, 2, int16) == 9);
    CHECK(reader.GetDataset()->experiment.GetSaturationLimit() == 32767);
    reader.Close();
    std::filesystem::remove(name);
}

TEST_CASE("JFJochXFELReader_GeometryMasks", "[portable]") {
    // A Cheetah CXI file of two 3-column x 4-row panels side by side in a 4 x 6 frame, read through a
    // CrystFEL geometry file that names the masks: the per-frame mask with its own bad bits, a mask in a
    // file beside the geometry file, a box in the laboratory plane and a pixel range on panel p1. The
    // file's own mask_shared is not taken (the geometry names the masks); --mask-bad still adds its own.
    const auto dir = std::filesystem::path("xfel_reader_geom_masks");
    std::filesystem::create_directories(dir);
    const std::string name = (dir / "frames.cxi").string();
    {
        const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
        for (const char *g: {"entry_1", "entry_1/data_1", "entry_1/instrument_1", "entry_1/instrument_1/detector_1",
                             "entry_1/instrument_1/source_1"})
            H5Gclose(H5Gcreate2(f, g, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT));
        std::vector<int32_t> frames(2 * 24);
        for (size_t i = 0; i < frames.size(); i++)
            frames[i] = static_cast<int32_t>(i % 24) + 100;   // every pixel its own value, in both frames
        WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT32, frames.data(), {2, 4, 6});
        std::vector<uint16_t> mask(2 * 24, 0);
        mask[1 * 6 + 1] = 0x20;     // frame 0 only, (fs 1, ss 1): a bit Cheetah's own set does not include
        WriteFrames(f, "entry_1/data_1/mask", H5T_NATIVE_UINT16, mask.data(), {2, 4, 6});
        // Frame 1's mask chunk is stored with both its filters skipped: HDF5 reads that one.
        const hid_t ds = H5Dopen2(f, "entry_1/data_1/mask", H5P_DEFAULT);
        const hsize_t offset[3] = {1, 0, 0};
        REQUIRE(H5Dwrite_chunk(ds, H5P_DEFAULT, 0x3, offset, 24 * sizeof(uint16_t), mask.data() + 24) >= 0);
        H5Dclose(ds);
        const std::string det = "entry_1/instrument_1/detector_1/";
        std::vector<uint16_t> shared(24, 0);
        shared[3 * 6 + 5] = 8;      // (fs 5, ss 3) dead to Cheetah
        WriteFrames(f, det + "mask_shared", H5T_NATIVE_UINT16, shared.data(), {4, 6});
        WriteVector(f, "entry_1/instrument_1/source_1/energy", std::vector<double>(2, 9000.0 * 1.602176634e-19));
        H5Fclose(f);
    }
    {
        // Good where bit 0 is set; bit 1 is what --mask-bad below looks at.
        const hid_t f = H5Fcreate((dir / "pixels.h5").string().c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
        std::vector<uint16_t> pixels(24, 1);
        pixels[0 * 6 + 4] = 0;      // (fs 4, ss 0)
        pixels[2 * 6 + 3] = 3;      // (fs 3, ss 2)
        WriteFrames(f, "mask", H5T_NATIVE_UINT16, pixels.data(), {4, 6});
        H5Fclose(f);
    }
    const std::string geom = (dir / "frames.geom").string();
    std::ofstream(geom) << R"(
clen = 0.1
res = 10000
dim0 = %
dim1 = ss
dim2 = fs
data = /entry_1/data_1/data
mask = /entry_1/data_1/mask
mask_good = 0x0
mask_bad = 0xffff
mask1_data = /mask
mask1_file = pixels.h5
mask1_goodbits = 0x1
p0/min_fs = 0
p0/max_fs = 2
p0/min_ss = 0
p0/max_ss = 3
p0/fs = x
p0/ss = y
p0/corner_x = -3
p0/corner_y = -2
p1/min_fs = 3
p1/max_fs = 5
p1/min_ss = 0
p1/max_ss = 3
p1/fs = x
p1/ss = y
p1/corner_x = 0
p1/corner_y = -2
badbox/min_x = -3
badbox/max_x = -2
badbox/min_y = -2
badbox/max_y = -1
bad_p1/min_fs = 5
bad_p1/max_fs = 5
bad_p1/min_ss = 0
bad_p1/max_ss = 0
bad_p1/panel = p1
)";

    JFJochXFELReader reader;
    reader.ExtraMasks({{"/mask", (dir / "pixels.h5").string(), 0, 2}});
    REQUIRE_NOTHROW(reader.ReadFiles({name}, geom));
    REQUIRE(reader.GetNumberOfImages() == 2);

    // Where each stored pixel lands, found by its value in frame 1 (which the per-frame mask leaves alone).
    auto raw = reader.GetRawImage(1);
    REQUIRE(raw);
    const auto *frame1 = reinterpret_cast<const int32_t *>(raw->image_buffer.data());
    const size_t n = raw->image.GetWidth() * raw->image.GetHeight();
    auto where = [&](size_t fs, size_t ss) {
        const auto it = std::find(frame1, frame1 + n, static_cast<int32_t>(ss * 6 + fs) + 100);
        REQUIRE(it != frame1 + n);
        return static_cast<size_t>(it - frame1);
    };
    const auto &m = reader.GetDataset()->pixel_mask->GetMask();
    CHECK(m[where(2, 2)] == 0);
    CHECK(m[where(4, 0)] != 0);     // mask 1, beside the geometry file: lacks the good bit
    CHECK(m[where(3, 2)] != 0);     // --mask-bad: bit 1
    CHECK(m[where(5, 3)] == 0);     // the file's mask_shared: not taken
    CHECK(m[where(5, 1)] == 0);
    CHECK(m[where(0, 0)] != 0);     // the lab box: p0's first pixel (centre at x -2.5, y -1.5) and no other
    CHECK(m[where(1, 0)] == 0);
    CHECK(m[where(0, 1)] == 0);
    CHECK(m[where(5, 0)] != 0);     // bad_p1
    size_t usable = 0;
    for (size_t ss = 0; ss < 4; ss++)
        for (size_t fs = 0; fs < 6; fs++)
            usable += m[where(fs, ss)] == 0;
    CHECK(usable == 24 - 4);

    // The per-frame mask: (1, 1) is no measurement in frame 0 only.
    auto raw0 = reader.GetRawImage(0);
    REQUIRE(raw0);
    const auto *frame0 = reinterpret_cast<const int32_t *>(raw0->image_buffer.data());
    CHECK(frame0[where(1, 1)] == INT32_MIN);
    CHECK(frame0[where(2, 1)] == 1 * 6 + 2 + 100);
    reader.Close();
    std::filesystem::remove_all(dir);
}

TEST_CASE("JFJochXFELReader_DeflatedInPieces", "[portable]") {
    // Frames large enough to be inflated in several pieces (a 300 x 300 int32 frame and its uint16 per-frame
    // mask, deflated without the shuffle), read with a gain: every pixel is the stored value in photons,
    // rounded half away from zero, or no measurement where the per-frame mask marks it.
    const auto dir = std::filesystem::path("xfel_reader_in_pieces");
    std::filesystem::create_directories(dir);
    const std::string name = (dir / "frames.cxi").string();
    constexpr size_t side = 300, n = side * side;
    std::vector<int32_t> frames(2 * n);
    for (size_t i = 0; i < frames.size(); i++)
        frames[i] = static_cast<int32_t>(i % 1001) - 500;
    std::vector<uint16_t> mask(2 * n, 0);
    for (size_t i = 0; i < mask.size(); i += 997)
        mask[i] = 0x4;
    {
        const hid_t f = H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
        for (const char *g: {"entry_1", "entry_1/data_1", "entry_1/instrument_1", "entry_1/instrument_1/source_1"})
            H5Gclose(H5Gcreate2(f, g, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT));
        WriteFrames(f, "entry_1/data_1/data", H5T_NATIVE_INT32, frames.data(), {2, side, side}, false);
        WriteFrames(f, "entry_1/data_1/mask", H5T_NATIVE_UINT16, mask.data(), {2, side, side}, false);
        WriteVector(f, "entry_1/instrument_1/source_1/energy", std::vector<double>(2, 9000.0 * 1.602176634e-19));
        H5Fclose(f);
    }
    const std::string geom = (dir / "frames.geom").string();
    std::ofstream(geom) << R"(
clen = 0.1
res = 10000
dim0 = %
dim1 = ss
dim2 = fs
data = /entry_1/data_1/data
mask = /entry_1/data_1/mask
mask_good = 0x0
mask_bad = 0xffff
p0/min_fs = 0
p0/max_fs = 149
p0/min_ss = 0
p0/max_ss = 299
p0/fs = x
p0/ss = y
p0/corner_x = -150
p0/corner_y = -150
p1/min_fs = 150
p1/max_fs = 299
p1/min_ss = 0
p1/max_ss = 299
p1/fs = x
p1/ss = y
p1/corner_x = 0
p1/corner_y = -150
)";

    JFJochXFELReader reader;
    REQUIRE_NOTHROW(reader.ReadFiles({name}, geom, 2.0f));
    REQUIRE(reader.GetNumberOfImages() == 2);
    for (size_t image = 0; image < 2; image++) {
        auto raw = reader.GetRawImage(static_cast<int64_t>(image));
        REQUIRE(raw);
        REQUIRE(raw->image.GetWidth() * raw->image.GetHeight() == n);
        const auto *pixels = reinterpret_cast<const int32_t *>(raw->image_buffer.data());
        size_t wrong = 0;
        for (size_t i = 0; i < n; i++) {
            const size_t k = image * n + i;
            const int32_t expected = mask[k] ? INT32_MIN : static_cast<int32_t>(std::lround(frames[k] / 2.0f));
            wrong += pixels[i] != expected;
        }
        CHECK(wrong == 0);
    }
    reader.Close();
    std::filesystem::remove_all(dir);
}

namespace {

// An LCLS-I XTC stream of the shape JFJochXTCReader reads: a Configure datagram, then one event (an
// L1Accept) per image, each a container holding one 16-bit camera frame (Camera::FrameV1).
void WriteXTC(const std::string &name, uint32_t width, uint32_t height, uint32_t pixel_offset,
              const std::vector<std::vector<uint16_t>> &images) {
    std::ofstream f(name, std::ios::binary);
    auto words = [&](std::initializer_list<uint32_t> w) {
        for (uint32_t v: w)
            f.write(reinterpret_cast<const char *>(&v), 4);
    };
    constexpr uint32_t ID_XTC = 1, ID_FRAME = 2, XTC_HEADER = 20;
    // clock (2 words), timestamp (2 words, the transition in bits 24-27 of the first), env, Xtc header
    words({0, 0, 0u << 24, 0, 0, 0, 0, 0, ID_XTC, XTC_HEADER});
    const uint32_t frame_extent = XTC_HEADER + 16 + 2 * width * height;
    for (const auto &image: images) {
        words({0, 0, 12u << 24, 0, 0, 0, 0, 0, ID_XTC, XTC_HEADER + frame_extent});
        words({0, 0, 7, ID_FRAME | (1u << 16), frame_extent});
        words({width, height, 16, pixel_offset});
        f.write(reinterpret_cast<const char *>(image.data()), static_cast<std::streamsize>(2 * image.size()));
    }
}

}  // namespace

TEST_CASE("JFJochXTCReader_StoredPixels", "[portable]") {
    // Two streams of 3 x 2 images. Without an offset to take off, the main loop's read hands the pixels up
    // as stored (UINT16_WIDENED), the same values the decoded int32 image holds; with one, it is taken off.
    const std::string a = "xtc_reader_a.xtc", b = "xtc_reader_b.xtc", c = "xtc_reader_c.xtc";
    WriteXTC(a, 3, 2, 0, {{1, 2, 3, 4, 5, 65535}, {7, 7, 7, 7, 7, 7}});
    WriteXTC(b, 3, 2, 0, {{9, 8, 7, 6, 5, 4}});
    WriteXTC(c, 3, 2, 10, {{11, 12, 13, 14, 15, 9}});
    REQUIRE(JFJochXTCReader::CanRead(a));

    JFJochXTCReader reader;
    REQUIRE_NOTHROW(reader.ReadFiles({a, b}));
    REQUIRE(reader.GetNumberOfImages() == 3);
    auto raw = reader.GetRawImage(2);
    REQUIRE(raw);
    CHECK(raw->image.GetCompressionAlgorithm() == CompressionAlgorithm::UINT16_WIDENED);
    std::vector<uint8_t> buffer;
    auto pixels = reinterpret_cast<const int32_t *>(raw->image.GetUncompressedPtr(buffer));
    CHECK(pixels[0] == 9);
    CHECK(pixels[5] == 4);
    raw = reader.GetRawImage(0);
    REQUIRE(raw);
    pixels = reinterpret_cast<const int32_t *>(raw->image.GetUncompressedPtr(buffer));
    CHECK(pixels[5] == 65535);

    REQUIRE_NOTHROW(reader.ReadFiles({c}));
    raw = reader.GetRawImage(0);
    REQUIRE(raw);
    CHECK(raw->image.GetCompressionAlgorithm() == CompressionAlgorithm::NO_COMPRESSION);
    CHECK(raw->image.GetMode() == CompressedImageMode::Int32);
    pixels = reinterpret_cast<const int32_t *>(raw->image_buffer.data());
    CHECK(pixels[0] == 1);
    CHECK(pixels[5] == -1);
    reader.Close();
    for (const auto &n: {a, b, c})
        std::filesystem::remove(n);
}