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

#include "JFJochXFELReader.h"
#include "CrystFELGeometry.h"

#include <algorithm>
#include <climits>
#include <cmath>
#include <type_traits>
#include <limits>
#include <cstring>
#include <filesystem>
#include <numeric>

#include <zlib.h>

#include "../common/JFJochException.h"
#include "../common/Logger.h"
#include "../writer/HDF5Objects.h"
#include "../compression/JFJochDecompress.h"

namespace {

constexpr double JOULE_PER_EV = 1.602176634e-19;
constexpr double PI_D = 3.14159265358979323846;
const std::string CXI_DATA = "/entry_1/data_1/data";
const std::string CXI_DETECTOR = "/entry_1/instrument_1/detector_1/";

// The bits of a Cheetah pixel mask that mark a pixel unusable in every frame: invalid, hot, dead,
// shadowed, bad, missing and noisy (Cheetah's PIXEL_IS_* flags). The peak mask, saturation and the
// resolution limits are left alone.
constexpr uint32_t CHEETAH_BAD_BITS = 1u | 4u | 8u | 16u | 128u | 512u | 1024u;

std::string ReplaceEvent(std::string path, const std::string &group) {
    const auto pos = path.find('%');
    if (pos != std::string::npos)
        path.replace(pos, 1, group);
    return path;
}

// A running mean, and the spread, of a value the geometry takes from the image files.
struct Mean {
    double sum = 0, min = INFINITY, max = -INFINITY;
    size_t n = 0;
    void Add(double v) { sum += v; n++; min = std::min(min, v); max = std::max(max, v); }
    double Get() const { return sum / static_cast<double>(n); }
};

// Adds what one file states at path: one value per shot group (a path with %), a per-frame vector, or
// one scalar for the file.
void Accumulate(HDF5ReadOnlyFile &file, const std::string &path, const std::vector<std::string> &groups,
                Mean &mean) {
    if (path.find('%') != std::string::npos) {
        for (const auto &g: groups)
            mean.Add(HDF5DataSet(file, ReplaceEvent(path, g)).ReadScalar<double>());
    } else if (file.GetDimension(path).empty())
        mean.Add(HDF5DataSet(file, path).ReadScalar<double>());
    else
        for (const double v: file.ReadVector<double>(path))
            mean.Add(v);
}

// b where use_b, otherwise a - as arithmetic rather than a branch, which the compiler would move the work
// computing a into, and then not vectorise the loop over the pixels.
int32_t Pick(bool use_b, int32_t a, int32_t b) {
    const int32_t all = -static_cast<int32_t>(use_b);
    return (a & ~all) | (b & all);
}

// Pixels are decoded and converted a piece at a time where they can be (ReadStored): a piece fits the core's
// own cache with the window inflate keeps, so it is converted before it leaves.
constexpr size_t PIECE_BYTES = 128 * 1024;

// The pixels [first, first + count) of a byte-shuffled frame of n pixels of size bytes each (byte b of pixel
// i is shuffled[b * n + i]), into out as stored: 2 or 4 bytes (plain_type), little-endian, as the frames are read.
void Unshuffle(const uint8_t *shuffled, size_t n, size_t size, size_t first, size_t count, uint8_t *out) {
    const uint8_t *b0 = shuffled + first, *b1 = b0 + n;
    if (size == 2) {
        auto *o = reinterpret_cast<uint16_t *>(out);
        for (size_t i = 0; i < count; i++)
            o[i] = static_cast<uint16_t>(b0[i] | b1[i] << 8);
    } else {
        const uint8_t *b2 = b1 + n, *b3 = b2 + n;
        auto *o = reinterpret_cast<uint32_t *>(out);
        for (size_t i = 0; i < count; i++)
            o[i] = b0[i] | b1[i] << 8 | b2[i] << 16 | static_cast<uint32_t>(b3[i]) << 24;
    }
}

}  // namespace

JFJochXFELReader::JFJochXFELReader() = default;

JFJochXFELReader::~JFJochXFELReader() {
    Close();
}

// Whether the file is HDF5 at all, asked without HDF5 printing its error stack for one that is not (an
// XTC stream offered to every reader in turn). Under hdf5_mutex.
static bool IsHDF5(const std::string &path) {
    htri_t r = -1;
    H5E_BEGIN_TRY {
        r = H5Fis_accessible(path.c_str(), H5P_DEFAULT);
    } H5E_END_TRY
    return r > 0;
}

// The HDF5 filter IDs registered for bitshuffle and for plain LZ4.
constexpr H5Z_filter_t BITSHUFFLE_FILTER_ID = 32008;
constexpr H5Z_filter_t LZ4_FILTER_ID = 32004;

std::shared_ptr<const HDF5ImageSource::RawFile> JFJochXFELReader::RawFileAt(size_t file) const {
    if (raw_files_.size() != files_.size()) {
        raw_files_.assign(files_.size(), nullptr);
        user_block_.assign(files_.size(), 0);
        raw_order_.clear();
    }
    if (!raw_files_[file]) {
        const hid_t fcpl = H5Fget_create_plist(FileAt(file).GetID());
        hsize_t user_block = 0;
        if (fcpl >= 0) {
            H5Pget_userblock(fcpl, &user_block);
            H5Pclose(fcpl);
        }
        auto raw = std::make_shared<const HDF5ImageSource::RawFile>(files_[file]);
        if (!raw->IsOpen())
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot open " + files_[file]);
        raw_files_[file] = raw;
        user_block_[file] = user_block;
        raw_order_.push_back(file);
        if (raw_order_.size() > MAX_OPEN_FILES) {
            raw_files_[raw_order_.front()].reset();
            raw_order_.erase(raw_order_.begin());
        }
    }
    return raw_files_[file];
}

HDF5ReadOnlyFile &JFJochXFELReader::FileAt(size_t file) const {
    auto it = std::find_if(open_files_.begin(), open_files_.end(), [&](const auto &e) { return e.first == file; });
    if (it == open_files_.end()) {
        if (open_files_.size() >= MAX_OPEN_FILES)
            open_files_.erase(open_files_.begin());
        open_files_.emplace_back(file, std::make_unique<HDF5ReadOnlyFile>(files_[file]));
    } else if (it != open_files_.end() - 1)
        std::rotate(it, it + 1, open_files_.end());
    return *open_files_.back().second;
}

// The files given, in the order given; a directory is not read (the caller expands a glob).
static void CheckFiles(const std::vector<std::string> &files) {
    if (files.empty())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No file to read");
}

void JFJochXFELReader::ReadFiles(const std::vector<std::string> &files, const std::string &geometry_file,
                                     std::optional<float> gain) {
    Close();
    CheckFiles(files);
    std::vector<std::string> ignored;
    geom_ = crystfel::Read(geometry_file, ignored);
    source_ = geometry_file;
    if (!ignored.empty()) {
        std::string list;
        for (const auto &key: ignored)
            list += (list.empty() ? "" : ", ") + key;
        Logger("XFELReader").Info("{}: the module geometry, camera length and masks are taken; ignored: {}",
                                  geometry_file, list);
    }
    Logger("XFELReader").Info("{}: {} mask array(s), {} bad region(s)", geometry_file, geom_.masks.size(),
                              geom_.bad.size());

    // What the geometry file does not state comes from the data files themselves, where they are one of
    // the layouts this reader knows: a Cheetah CXI file (photon energy, camera length, the shared mask
    // and the per-shot one) or a SACLA Cheetah file (photon energy per shot, camera length, pixel mask).
    // The file's own masks are taken only where the geometry names none: a geometry that does names the
    // mask arrays with the bits that matter for this detector. Otherwise the photon energy is the
    // command line's to give, and the camera length too where the geometry has none.
    {
        std::unique_lock ul(hdf5_mutex);
        HDF5ReadOnlyFile file(files[0]);
        const bool file_masks = geom_.masks.empty();
        if (file.Exists(CXI_DATA)) {
            if (file.Exists("/entry_1/instrument_1/source_1/energy"))
                geom_.photon_energy = xfel::Value{"/entry_1/instrument_1/source_1/energy", 0, 1.0 / JOULE_PER_EV};
            if (file_masks && file.Exists(CXI_DETECTOR + "mask_shared"))
                geom_.masks.push_back({CXI_DETECTOR + "mask_shared", "", 0, CHEETAH_BAD_BITS});
            if (file_masks && file.Exists("/entry_1/data_1/mask"))
                geom_.masks.push_back({"/entry_1/data_1/mask", "", 0, CHEETAH_BAD_BITS});
            if (!geom_.panels[0].clen && file.Exists(CXI_DETECTOR + "distance"))
                for (auto &p: geom_.panels)
                    p.clen = xfel::Value{CXI_DETECTOR + "distance", 0, 1.0};
        } else if (file.Exists("/metadata/distance_in_mm")) {
            geom_.photon_energy = xfel::Value{"/%/photon_energy_ev", 0, 1.0};
            if (file_masks)
                geom_.masks.push_back({"/metadata/pixelmask", "", 0, 0xffffffffu});
            if (!geom_.panels[0].clen)
                for (auto &p: geom_.panels)
                    p.clen = xfel::Value{"/metadata/distance_in_mm", 0, 1e-3};
        }
    }
    for (auto &p: geom_.panels) {
        if (!p.clen)
            p.clen = xfel::Value{"", 0.1, 1.0};   // a placeholder: the distance is given on the command line
        if (gain)
            p.adu_per_photon = *gain;
    }
    Open(files);
    Place();
}

// A Cheetah CXI file carries what one flat panel needs - the frames at /entry_1/data_1/data, the pixel
// size, the photon energy and a camera length, all in SI units - but not where the beam is. The panel is
// put with its middle on the beam, which is as good as no centre at all: the run measures it.
//
// Only a square frame is taken this way. A Cheetah file from a tiled detector (CSPAD, AGIPD, ePix,
// JUNGFRAU, MPCCD) holds its panels stacked in a slab, which read as one flat panel would place every
// pixel wrongly and say nothing; those slabs are not square, the single-panel CCDs are.
bool JFJochXFELReader::CanReadCXI(const std::string &path) {
    try {
        if (!std::filesystem::is_regular_file(path))
            return false;
        std::unique_lock ul(hdf5_mutex);
        if (!IsHDF5(path))
            return false;
        HDF5ReadOnlyFile file(path);
        if (!file.Exists(CXI_DATA))
            return false;
        const auto dim = file.GetDimension(CXI_DATA);
        return dim.size() == 3 && dim[1] == dim[2]
               && file.Exists(CXI_DETECTOR + "x_pixel_size")
               && file.Exists("/entry_1/instrument_1/source_1/energy");
    } catch (...) {
        return false;
    }
}

std::optional<std::string> JFJochXFELReader::CXIRefusal(const std::string &path) {
    try {
        if (!std::filesystem::is_regular_file(path) || CanReadCXI(path))
            return std::nullopt;
        std::unique_lock ul(hdf5_mutex);
        if (!IsHDF5(path))
            return std::nullopt;
        HDF5ReadOnlyFile file(path);
        if (!file.Exists(CXI_DATA))
            return std::nullopt;
        const auto dim = file.GetDimension(CXI_DATA);
        if (dim.size() == 3 && dim[1] != dim[2])
            return fmt::format("{} is a Cheetah CXI file whose {} x {} frame is the panels of a tiled detector "
                               "stacked in a slab; their layout is not in the file - give it with "
                               "--crystfel-geometry", path, dim[2], dim[1]);
        return fmt::format("{} is a CXI file without the pixel size and photon energy of a single-panel "
                           "detector; give its geometry with --crystfel-geometry", path);
    } catch (...) {
        return std::nullopt;
    }
}

void JFJochXFELReader::ReadCXI(const std::vector<std::string> &files, std::optional<float> gain) {
    Close();
    CheckFiles(files);
    std::vector<hsize_t> dim;
    bool has_mask = false;
    {
        std::unique_lock ul(hdf5_mutex);
        HDF5ReadOnlyFile file(files[0]);
        dim = file.GetDimension(CXI_DATA);
        has_mask = file.Exists(CXI_DETECTOR + "mask_shared");
        const auto pixel = file.ReadVector<double>(CXI_DETECTOR + "x_pixel_size");
        xfel::Panel p;
        p.name = "cxi";
        p.max_ss = static_cast<int64_t>(dim[1]) - 1;
        p.max_fs = static_cast<int64_t>(dim[2]) - 1;
        p.corner_x = -0.5 * static_cast<double>(dim[2]);
        p.corner_y = -0.5 * static_cast<double>(dim[1]);
        p.res = 1.0 / pixel.at(0);
        if (file.Exists(CXI_DETECTOR + "distance"))
            p.clen = xfel::Value{CXI_DETECTOR + "distance", 0, 1.0};
        else
            p.clen = xfel::Value{"", 0.1, 1.0};   // a placeholder: the distance is given or found
        if (gain)
            p.adu_per_photon = *gain;
        geom_ = {};
        geom_.panels.push_back(p);
        geom_.data = CXI_DATA;
        geom_.photon_energy = xfel::Value{"/entry_1/instrument_1/source_1/energy", 0, 1.0 / JOULE_PER_EV};
        if (has_mask)
            geom_.masks.push_back({CXI_DETECTOR + "mask_shared", "", 0, CHEETAH_BAD_BITS});
    }
    source_ = "the CXI metadata";
    Open(files);
    Place();
}

// A SACLA Cheetah file: one group per shot (tag-N) with the frame and the shot's photon energy, and
// /metadata stating where every sensor of the MPCCD lies (posx/posy in um, the in-plane angle - in
// degrees, whatever its name says - and the pixel size), the camera length and a pixel mask that
// already holds the sensors' edges and the rows behind each amplifier's shield. The sensors are stacked
// one after the other along the frame's slow axis.
bool JFJochXFELReader::CanReadPanelTable(const std::string &path) {
    try {
        if (!std::filesystem::is_regular_file(path))
            return false;
        std::unique_lock ul(hdf5_mutex);
        if (!IsHDF5(path))
            return false;
        HDF5ReadOnlyFile file(path);
        return file.Exists("/metadata/posx_in_um") && file.Exists("/metadata/posy_in_um")
               && file.Exists("/metadata/angle_in_rad") && file.Exists("/metadata/pixelsizex_in_um")
               && file.Exists("/metadata/distance_in_mm");
    } catch (...) {
        return false;
    }
}

void JFJochXFELReader::ReadPanelTable(const std::vector<std::string> &files, std::optional<float> gain) {
    Close();
    CheckFiles(files);
    struct Metrology {
        std::vector<double> posx, posy, angle, pixel;
        double distance_mm = 0;
        bool operator==(const Metrology &) const = default;
    };
    auto read_metrology = [](HDF5ReadOnlyFile &file) {
        Metrology m;
        m.posx = file.ReadVector<double>("/metadata/posx_in_um");
        m.posy = file.ReadVector<double>("/metadata/posy_in_um");
        m.angle = file.ReadVector<double>("/metadata/angle_in_rad");
        m.pixel = file.ReadVector<double>("/metadata/pixelsizex_in_um");
        m.distance_mm = HDF5DataSet(file, "/metadata/distance_in_mm").ReadScalar<double>();
        return m;
    };

    Metrology m;
    std::string first_group;
    {
        std::unique_lock ul(hdf5_mutex);
        HDF5ReadOnlyFile file(files[0]);
        m = read_metrology(file);
        // The sensors of one run are placed once: every file has to state the same metrology.
        for (size_t f = 1; f < files.size(); f++) {
            HDF5ReadOnlyFile other(files[f]);
            if (!(read_metrology(other) == m))
                throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                      files[f] + " places the sensors or the detector differently from "
                                      + files[0] + "; read the runs separately");
        }
        for (const auto &g: file.FindLeafs("/"))
            if (file.Exists("/" + g + "/data")) {
                first_group = g;
                break;
            }
        if (first_group.empty())
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No shot groups in " + files[0]);
    }
    const size_t nsensor = m.posx.size();
    if (nsensor == 0 || m.posy.size() != nsensor || m.angle.size() != nsensor || m.pixel.size() != nsensor)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              files[0] + ": /metadata does not state the same number of positions, angles and "
                                         "pixel sizes");
    if (!(m.distance_mm > 0))
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, files[0] + " states no camera length");

    std::vector<hsize_t> dim;
    {
        std::unique_lock ul(hdf5_mutex);
        HDF5ReadOnlyFile file(files[0]);
        dim = file.GetDimension("/" + first_group + "/data");
    }
    if (dim.size() != 2 || dim[0] % nsensor != 0)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              files[0] + ": the frame does not hold " + std::to_string(nsensor) + " sensors stacked");
    const auto rows = static_cast<int64_t>(dim[0] / nsensor);

    geom_ = {};
    geom_.data = "/%/data";
    geom_.photon_energy = xfel::Value{"/%/photon_energy_ev", 0, 1.0};
    geom_.masks.push_back({"/metadata/pixelmask", "", 0, 0xffffffffu});
    // The placement the SACLA pipeline itself writes into its CrystFEL geometry: the corner of the
    // sensor's first pixel at (-posx, posy) and its fast axis at the stated angle.
    for (size_t i = 0; i < nsensor; i++) {
        const double a = m.angle[i] * PI_D / 180.0;
        xfel::Panel p;
        p.name = "q" + std::to_string(i + 1);
        p.min_ss = static_cast<int64_t>(i) * rows;
        p.max_ss = p.min_ss + rows - 1;
        p.max_fs = static_cast<int64_t>(dim[1]) - 1;
        p.corner_x = -m.posx[i] / m.pixel[i];
        p.corner_y = m.posy[i] / m.pixel[i];
        p.fs_x = std::cos(a);
        p.fs_y = std::sin(a);
        p.ss_x = -std::sin(a);
        p.ss_y = std::cos(a);
        p.res = 1e6 / m.pixel[i];
        p.clen = xfel::Value{"", m.distance_mm * 1e-3, 1.0};
        if (gain)
            p.adu_per_photon = *gain;
        geom_.panels.push_back(p);
    }
    source_ = "the /metadata panel table";
    Open(files);
    Place();
}

JFJochXFELReader::Layout JFJochXFELReader::ReadLayout(HDF5ReadOnlyFile &file, const Frame &frame) {
    Layout l;
    HDF5DataSet dataset(file, frame.dataset);
    const auto dim = file.GetDimension(frame.dataset);
    if (dim.size() < 2 || dim.size() > 3)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              frame.dataset + " is not a stack of frames or one frame");
    l.height = dim[dim.size() - 2];
    l.width = dim[dim.size() - 1];

    const hid_t dcpl = H5Dget_create_plist(dataset.GetID());
    const int rank = static_cast<int>(dim.size());
    std::vector<hsize_t> chunk(rank, 0);
    bool ok = H5Pget_layout(dcpl) == H5D_CHUNKED && H5Pget_chunk(dcpl, rank, chunk.data()) == rank
              && chunk[rank - 1] == l.width && chunk[rank - 2] == l.height && (rank == 2 || chunk[0] == 1);
    bool deflate = false, bitshuffle = false;
    for (int i = 0; ok && i < H5Pget_nfilters(dcpl); i++) {
        unsigned int flags;
        size_t nelements = 0;
        const H5Z_filter_t filter = H5Pget_filter2(dcpl, i, &flags, &nelements, nullptr, 0, nullptr, nullptr);
        if (filter == H5Z_FILTER_DEFLATE)
            deflate = true;
        else if (filter == H5Z_FILTER_SHUFFLE)
            l.shuffled = true;
        else if (filter == BITSHUFFLE_FILTER_ID || filter == LZ4_FILTER_ID)
            bitshuffle = true;
        else
            ok = false;
    }
    H5Pclose(dcpl);
    // Bitshuffle/LZ4 (or zstd, or plain LZ4) is the one filter of its chunk: CompressedImage decodes it.
    if (ok && bitshuffle && !deflate && !l.shuffled)
        l.algorithm = HDF5Dcpl(dataset).GetCompression();

    const hid_t type = H5Dget_type(dataset.GetID());
    l.stored_size = H5Tget_size(type);
    l.stored_float = H5Tget_class(type) == H5T_FLOAT;
    l.stored_signed = H5Tget_sign(type) == H5T_SGN_2;
    l.plain_type = ((H5Tget_class(type) == H5T_INTEGER && (l.stored_size == 2 || l.stored_size == 4))
                    || (l.stored_float && l.stored_size == 4))
                   && H5Tget_order(type) == H5T_ORDER_LE;
    ok &= l.plain_type;
    H5Tclose(type);
    l.deflated = deflate;
    l.direct_chunks = ok && (deflate || (bitshuffle && l.algorithm != CompressionAlgorithm::NO_COMPRESSION));
    return l;
}

// A per-frame mask is read as the frames are where it is stored in plain unsigned integers - read as
// stored, they are what HDF5 would convert them to - of the frame's size, in every file the way the
// first file stores it. Any other is read by HDF5.
JFJochXFELReader::Layout JFJochXFELReader::FrameMaskLayout(const xfel::Mask &m) const {
    Layout layout;
    for (size_t i = 0; i < frames_.size(); i++) {
        if (i > 0 && frames_[i].file == frames_[i - 1].file)
            continue;   // a file's first frame stands for all of its frames, as for the frames themselves
        const Frame mask{frames_[i].file, ReplaceEvent(m.data, frames_[i].group), frames_[i].index, frames_[i].group};
        if (!FileAt(mask.file).Exists(mask.dataset))
            return {};
        const Layout l = ReadLayout(FileAt(mask.file), mask);
        if (i == 0)
            layout = l;
        else if (!(l == layout))
            return {};
    }
    if (layout.stored_float || layout.stored_signed || layout.height != frame_height_ || layout.width != frame_width_)
        return {};
    return layout;
}

void JFJochXFELReader::Open(const std::vector<std::string> &files) {
    files_ = files;
    geom_.masks.insert(geom_.masks.end(), extra_masks_.begin(), extra_masks_.end());

    // Every frame of every file, in file order, with the mean of every value taken from the files.
    Mean energy, clen;
    const auto &p0 = geom_.panels.at(0);
    const auto pct = geom_.data.find('%');
    for (const auto &p: geom_.panels)
        if (p.res != p0.res || p.clen->path != p0.clen->path || p.clen->value != p0.clen->value)
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                  source_ + ": panels with different camera lengths or pixel sizes cannot be put "
                                            "on one flat detector");
    const xfel::Value *energy_value = geom_.photon_energy ? &*geom_.photon_energy
                                                              : (geom_.wavelength ? &*geom_.wavelength : nullptr);

    std::unique_lock ul(hdf5_mutex);
    for (size_t f = 0; f < files_.size(); f++) {
        HDF5ReadOnlyFile file(files_[f]);
        const size_t first = frames_.size();
        std::vector<std::string> groups;
        if (pct != std::string::npos) {
            const std::string parent = geom_.data.substr(0, pct);
            for (const auto &group: file.FindLeafs(parent.empty() ? "/" : parent)) {
                const auto dataset = ReplaceEvent(geom_.data, group);
                if (file.Exists(dataset)) {
                    frames_.push_back({f, dataset, -1, group});
                    groups.push_back(group);
                }
            }
        } else if (file.Exists(geom_.data)) {
            const auto dim = file.GetDimension(geom_.data);
            if (dim.size() == 3) {
                for (hsize_t i = 0; i < dim[0]; i++)
                    frames_.push_back({f, geom_.data, static_cast<int64_t>(i), ""});
            } else if (dim.size() == 2)
                frames_.push_back({f, geom_.data, -1, ""});
        }
        if (frames_.size() == first)
            continue;

        // Every file has to store its frames the way the first does: the decoding is set up once.
        const Layout l = ReadLayout(file, frames_[first]);
        if (first == 0)
            layout_ = l;
        else if (!(l == layout_))
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                  files_[f] + " stores its frames differently (size, type or compression) from "
                                  + files_[0] + "; read them separately");

        if (energy_value != nullptr && !energy_value->path.empty())
            Accumulate(file, energy_value->path, groups, energy);
        if (!p0.clen->path.empty())
            Accumulate(file, p0.clen->path, groups, clen);
    }
    if (frames_.empty())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "No frames at " + geom_.data + " in the files given");
    frame_height_ = layout_.height;
    frame_width_ = layout_.width;

    clen_m_ = p0.clen->path.empty() ? p0.clen->value : clen.Get() * p0.clen->scale;
    double coffset_sum = 0;
    for (const auto &p: geom_.panels)
        coffset_sum += p.coffset;
    distance_mm_ = (clen_m_ + coffset_sum / static_cast<double>(geom_.panels.size())) * 1000.0;
    // Nothing states the photon energy: it is the command line's to give (--wavelength).
    photon_energy_eV_ = 0;
    if (energy_value != nullptr) {
        const double v = energy_value->path.empty() ? energy_value->value : energy.Get() * energy_value->scale;
        photon_energy_eV_ = geom_.photon_energy ? v : WVL_1A_IN_KEV * 1000.0 / v;
        if (!(photon_energy_eV_ > 0) || !std::isfinite(photon_energy_eV_))
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                  source_ + ": no usable photon energy in the files given");
    }
    if (energy.n > 1)
        Logger("XFELReader").Info("Photon energy {:.1f} eV, the mean of {} values from {:.1f} to {:.1f} (the run "
                                      "is read at one wavelength)", photon_energy_eV_, energy.n,
                                      geom_.photon_energy ? energy.min * energy_value->scale
                                                          : WVL_1A_IN_KEV * 1000.0 / (energy.max * energy_value->scale),
                                      geom_.photon_energy ? energy.max * energy_value->scale
                                                          : WVL_1A_IN_KEV * 1000.0 / (energy.min * energy_value->scale));
}

void JFJochXFELReader::Place() {
    std::vector<uint32_t> file_mask(frame_width_ * frame_height_, 0);
    std::vector<uint8_t> masked_by_file(frame_width_ * frame_height_, 0);
    {
        std::unique_lock ul(hdf5_mutex);
        frame_masks_.clear();
        frame_mask_layouts_.clear();
        for (const auto &m: geom_.masks) {
            std::vector<uint32_t> v;
            if (m.file.empty() && (m.data.find('%') != std::string::npos
                                   || FileAt(0).GetDimension(m.data).size() == 3)) {
                frame_masks_.push_back(m);    // one per shot or one plane per frame: applied as each is read
                frame_mask_layouts_.push_back(FrameMaskLayout(m));
                continue;
            }
            if (!m.file.empty()) {
                HDF5ReadOnlyFile mask_file(m.file);
                v = mask_file.ReadVector<uint32_t>(m.data, {0, 0}, {frame_height_, frame_width_});
            } else
                v = FileAt(0).ReadVector<uint32_t>(m.data, {0, 0}, {frame_height_, frame_width_});
            MarkBad(m, v.data(), 0, v.size(), masked_by_file);
        }
    }

    // CrystFEL's laboratory frame and the detector plane's (column, row, beam) are both right-handed with
    // the beam along z, so the point (X, Y) lies at column beam_x + X, row beam_y + Y of the detector plane
    // - and CrystFEL's horizontal x stays the image's horizontal, where the polarisation correction expects
    // it. A pixel centre is at its index. A panel set nearer or further than the plane (its own coffset)
    // is projected onto it along the ray from the sample, so each pixel keeps its scattering direction.
    auto position = [this](const xfel::Panel &p, double fs, double ss) {
        const double scale = distance_mm_ / ((clen_m_ + p.coffset) * 1000.0);
        return std::pair{scale * (p.corner_x + fs * p.fs_x + ss * p.ss_x),
                         scale * (p.corner_y + fs * p.fs_y + ss * p.ss_y)};
    };
    // A detector of several panels is analysed as stored: the image is the stored frame, and each panel a
    // module that maps its box of the frame onto the detector plane exactly where the layout puts it.
    // A single panel goes onto an image of its own whole, as a block: its axes turned to the nearest image
    // axis and its first pixel put at a whole pixel. Where it really lies - the fraction of a pixel, the small
    // turn - is then its module's placement, so no pixel is moved on its own.
    const bool stored = geom_.panels.size() >= 2;
    struct Block {
        int64_t fx = 1, fy = 0, sx = 0, sy = 1;   // image step per +1 fs / +1 ss (single panel)
        double x0 = 0, y0 = 0;                    // plane position of the first pixel's centre, beam at 0
    };
    auto snap = [](double x, double y, int64_t &ix, int64_t &iy) {
        if (std::fabs(x) >= std::fabs(y)) { ix = x > 0 ? 1 : -1; iy = 0; }
        else { ix = 0; iy = y > 0 ? 1 : -1; }
    };
    std::vector<Block> blocks;
    double min_x = 1e9, min_y = 1e9, max_x = -1e9, max_y = -1e9;
    for (const auto &p: geom_.panels) {
        Block bl;
        if (!stored) {
            snap(p.fs_x, p.fs_y, bl.fx, bl.fy);
            snap(p.ss_x, p.ss_y, bl.sx, bl.sy);
            if (bl.fx * bl.sy - bl.fy * bl.sx == 0)
                throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                      "Panel " + p.name + " has its fast and slow axes along the same line");
        }
        std::tie(bl.x0, bl.y0) = position(p, 0.5, 0.5);
        const double w = static_cast<double>(p.max_fs - p.min_fs), h = static_cast<double>(p.max_ss - p.min_ss);
        for (double fs: {0.0, w})
            for (double ss: {0.0, h}) {
                const auto [x, y] = stored
                    ? position(p, fs + 0.5, ss + 0.5)
                    : std::pair{bl.x0 + fs * static_cast<double>(bl.fx) + ss * static_cast<double>(bl.sx),
                                bl.y0 + fs * static_cast<double>(bl.fy) + ss * static_cast<double>(bl.sy)};
                min_x = std::min(min_x, x); max_x = std::max(max_x, x);
                min_y = std::min(min_y, y); max_y = std::max(max_y, y);
            }
        blocks.push_back(bl);
    }
    // The first panel's first pixel centre on a whole column and row, two pixels of margin all round.
    const double beam_x = std::ceil(blocks[0].x0 - min_x) + 2 - blocks[0].x0;
    const double beam_y = std::ceil(blocks[0].y0 - min_y) + 2 - blocks[0].y0;
    width_ = stored ? frame_width_ : static_cast<size_t>(std::ceil(beam_x + max_x)) + 3;
    height_ = stored ? frame_height_ : static_cast<size_t>(std::ceil(beam_y + max_y)) + 3;

    // Everything starts as a gap; the panels' pixels are put on top, less the bad regions.
    std::vector<uint32_t> mask(width_ * height_, 1u << PixelMask::ModuleGapPixelBit);
    std::vector<uint8_t> taken(width_ * height_, 0);
    blocks_.clear();
    // One gain for every panel: what the run is given (ReadFiles, ApplyGain) is set on every panel alike.
    const auto &p0 = geom_.panels[0];
    counts_in_photons_ = p0.adu_per_photon.has_value() || p0.adu_per_eV.has_value();
    const double gain = p0.adu_per_photon.has_value() ? p0.adu_per_photon.value()
                        : p0.adu_per_eV.has_value() ? p0.adu_per_eV.value() * photon_energy_eV_ : 1.0;
    adu_per_photon_ = static_cast<float>(gain);
    std::optional<double> saturation;
    std::vector<DetectorModule> modules;

    for (size_t n = 0; n < geom_.panels.size(); n++) {
        const auto &p = geom_.panels[n];
        const Block &bl = blocks[n];
        if (p.max_fs >= static_cast<int64_t>(frame_width_) || p.max_ss >= static_cast<int64_t>(frame_height_))
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                  "Panel " + p.name + " lies outside the stored frame");
        if (p.max_adu.has_value())
            saturation = std::min(saturation.value_or(INT64_MAX), p.max_adu.value() / gain);
        if (p.no_index)
            continue;

        const int64_t w = p.max_fs - p.min_fs, h = p.max_ss - p.min_ss;
        const double scale = distance_mm_ / ((clen_m_ + p.coffset) * 1000.0);
        // The image pixel the panel's first pixel lands on (single panel).
        const int64_t c0 = std::lround(beam_x + bl.x0), r0 = std::lround(beam_y + bl.y0);
        DetectorModule m;
        if (stored) {
            // Its box in the stored frame, and the plane position and steps of its pixels: the panel's own axes,
            // scaled by its projection.
            m.x0 = p.min_fs;
            m.y0 = p.min_ss;
            m.width = w + 1;
            m.height = h + 1;
            m.ox = beam_x + bl.x0;
            m.oy = beam_y + bl.y0;
            m.fx = scale * p.fs_x;
            m.fy = scale * p.fs_y;
            m.sx = scale * p.ss_x;
            m.sy = scale * p.ss_y;
        } else {
            // The block it fills, and its placement: the plane position of image pixel (bx0, by0), and the
            // plane steps of +1 column and +1 row - the panel's own axes, scaled by its projection, taken
            // through the block's turn.
            const int64_t bx0 = c0 + std::min<int64_t>(0, bl.fx * w) + std::min<int64_t>(0, bl.sx * h);
            const int64_t by0 = r0 + std::min<int64_t>(0, bl.fy * w) + std::min<int64_t>(0, bl.sy * h);
            m.x0 = bx0;
            m.y0 = by0;
            m.width = std::abs(bl.fx * w + bl.sx * h) + 1;
            m.height = std::abs(bl.fy * w + bl.sy * h) + 1;
            // image (dc, dr) = fs (fx, fy) + ss (sx, sy); the inverse of a signed permutation is its transpose.
            auto plane_of = [&](double dc, double dr) {
                const double fs = dc * static_cast<double>(bl.fx) + dr * static_cast<double>(bl.fy);
                const double ss = dc * static_cast<double>(bl.sx) + dr * static_cast<double>(bl.sy);
                const auto [x, y] = position(p, fs + 0.5, ss + 0.5);
                return std::pair{beam_x + x, beam_y + y};
            };
            std::tie(m.ox, m.oy) = plane_of(static_cast<double>(bx0 - c0), static_cast<double>(by0 - r0));
            m.fx = scale * (static_cast<double>(bl.fx) * p.fs_x + static_cast<double>(bl.sx) * p.ss_x);
            m.fy = scale * (static_cast<double>(bl.fx) * p.fs_y + static_cast<double>(bl.sx) * p.ss_y);
            m.sx = scale * (static_cast<double>(bl.fy) * p.fs_x + static_cast<double>(bl.sy) * p.ss_x);
            m.sy = scale * (static_cast<double>(bl.fy) * p.fs_y + static_cast<double>(bl.sy) * p.ss_y);
            blocks_.push_back({static_cast<size_t>(p.min_ss) * frame_width_ + static_cast<size_t>(p.min_fs), w + 1,
                               h + 1, r0 * static_cast<int64_t>(width_) + c0,
                               bl.fy * static_cast<int64_t>(width_) + bl.fx,
                               bl.sy * static_cast<int64_t>(width_) + bl.sx});
        }
        modules.push_back(m);

        const int64_t e = p.mask_edge_pixels;
        for (int64_t ss = p.min_ss; ss <= p.max_ss; ss++) {
            for (int64_t fs = p.min_fs; fs <= p.max_fs; fs++) {
                const int64_t col = stored ? fs : c0 + (fs - p.min_fs) * bl.fx + (ss - p.min_ss) * bl.sx;
                const int64_t row = stored ? ss : r0 + (fs - p.min_fs) * bl.fy + (ss - p.min_ss) * bl.sy;
                const size_t dst = static_cast<size_t>(row) * width_ + static_cast<size_t>(col);
                if (taken[dst])
                    throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                          "Panel " + p.name + " overlaps another on the image");
                taken[dst] = 1;
                if (fs < p.min_fs + e || fs > p.max_fs - e || ss < p.min_ss + e || ss > p.max_ss - e)
                    continue;
                const auto [x, y] = position(p, fs - p.min_fs + 0.5, ss - p.min_ss + 0.5);
                bool bad = false;
                for (const auto &b: geom_.bad) {
                    if (b.in_lab)
                        bad |= x >= b.min_x && x <= b.max_x && y >= b.min_y && y <= b.max_y;
                    else if (b.panel.empty() || b.panel == p.name)
                        bad |= fs >= b.min_fs && fs <= b.max_fs && ss >= b.min_ss && ss <= b.max_ss;
                }
                if (bad)
                    continue;

                mask[dst] = masked_by_file[ss * frame_width_ + fs] ? (1u << PixelMask::ErrorPixelBit) : 0;
            }
        }
    }

    dataset_ = std::make_shared<JFJochReaderDataset>();
    dataset_->experiment = default_experiment;
    DetectorSetup detector = DetDECTRIS(static_cast<int64_t>(width_), static_cast<int64_t>(height_),
                                        "placed by " + source_, {});
    detector.PixelSize_um(static_cast<float>(1e6 / geom_.panels[0].res));
    // Nothing a parallax or sensor-efficiency correction needs is stated in these files.
    detector.SensorThickness_um(0);
    // Where nothing states a limit, the stored integer type still has one - a Cheetah file keeps its
    // frames as int16, and a value clipped at 32767 is no measurement either.
    if (!saturation.has_value() && !layout_.stored_float && layout_.stored_size > 0) {
        const double type_max = std::exp2(8.0 * static_cast<double>(layout_.stored_size)
                                          - (layout_.stored_signed ? 1 : 0)) - 1;
        saturation = type_max / adu_per_photon_;
    }
    if (saturation.has_value())
        detector.SaturationLimit(static_cast<int64_t>(std::floor(saturation.value())));
    detector.BitDepthImage(32);
    detector.MinFrameTime(std::chrono::microseconds(0));
    detector.MinCountTime(std::chrono::microseconds(0));
    detector.ReadOutTime(std::chrono::nanoseconds(0));
    dataset_->experiment.Detector(detector);
    // Dark-subtracted counts go below zero; the images are handed out as signed int32.
    dataset_->experiment.BitDepthImage(32);
    dataset_->experiment.PixelSigned(true);
    dataset_->experiment.BeamX_pxl(static_cast<float>(beam_x));
    dataset_->experiment.BeamY_pxl(static_cast<float>(beam_y));
    dataset_->experiment.DetectorDistance_mm(static_cast<float>(distance_mm_));
    if (photon_energy_eV_ > 0)
        dataset_->experiment.IncidentEnergy_keV(static_cast<float>(photon_energy_eV_ / 1000.0));
    dataset_->experiment.ImagesPerTrigger(static_cast<int64_t>(frames_.size()));
    // Each panel is a module, placed where the layout puts it; one panel square with the image needs none.
    const DetectorModules placed(modules);
    if (stored || placed.Moved())
        dataset_->experiment.Modules(placed);
    dataset_->pixel_mask = std::make_shared<const PixelMask>(mask);

    Logger("XFELReader").Info("{} frames in {} files, {} panels placed by {}; {} x {} image, beam at ({:.2f}, "
                                  "{:.2f}) px, {:.3f} mm, {:.1f} eV, {}{}", frames_.size(), files_.size(),
                                  geom_.panels.size(), source_, width_, height_, beam_x, beam_y, distance_mm_,
                                  photon_energy_eV_,
                                  counts_in_photons_ ? "converted to photons" : "in detector units (no gain known)",
                                  layout_.direct_chunks ? ", chunks decoded in parallel" : "");
    SetStartMessage(dataset_);
}

template <class T>
void JFJochXFELReader::MarkBad(const xfel::Mask &m, const T *v, size_t first, size_t n, std::vector<uint8_t> &bad) const {
    int64_t min_fs = 0, max_fs = static_cast<int64_t>(frame_width_) - 1;
    int64_t min_ss = 0, max_ss = static_cast<int64_t>(frame_height_) - 1;
    for (const auto &p: geom_.panels)
        if (!m.panel.empty() && p.name == m.panel) {
            min_fs = p.min_fs; max_fs = p.max_fs;
            min_ss = p.min_ss; max_ss = p.max_ss;
        }
    // The mask's bits in locals and plain pointers, a row at a time: a byte stored through bad could otherwise
    // be m's or bad's own, every pixel would read them again, and the loop would not be vectorised.
    const uint32_t good = m.good, bad_bits = m.bad;
    const bool nonzero_good = m.nonzero_good;
    const size_t last = first + n;   // the pixels [first, last) of the frame are those v holds
    for (int64_t ss = std::max<int64_t>(min_ss, first / frame_width_);
         ss <= max_ss && static_cast<size_t>(ss) * frame_width_ < last; ss++) {
        const size_t begin = std::max<size_t>(static_cast<size_t>(ss) * frame_width_ + min_fs, first);
        const size_t end = std::min<size_t>(static_cast<size_t>(ss) * frame_width_ + max_fs + 1, last);
        if (end <= begin)
            continue;
        const T *in = v + (begin - first);
        uint8_t *out = bad.data() + begin;
        for (size_t i = 0; i < end - begin; i++)
            out[i] |= xfel::Mask::Bad(in[i], good, bad_bits, nonzero_good);
    }
}

void JFJochXFELReader::ApplyGain(float adu_per_photon) {
    for (auto &p: geom_.panels)
        p.adu_per_photon = adu_per_photon;
    Place();
}

void JFJochXFELReader::ClearGain() {
    for (auto &p: geom_.panels)
        p.adu_per_photon.reset();
    Place();
}

bool JFJochXFELReader::StoredFloat() const {
    return layout_.stored_float || !layout_.plain_type;
}

bool JFJochXFELReader::HasPhotonEnergy() const {
    return photon_energy_eV_ > 0;
}

bool JFJochXFELReader::CountsInPhotons() const {
    return counts_in_photons_;
}

uint64_t JFJochXFELReader::GetNumberOfImages() const {
    return frames_.size();
}

void JFJochXFELReader::Close() {
    std::unique_lock ul(hdf5_mutex);
    open_files_.clear();
    raw_files_.clear();
    user_block_.clear();
    raw_order_.clear();
    files_.clear();
    frames_.clear();
    dataset_.reset();
}

template <class Buffer, class Piece>
bool JFJochXFELReader::ReadStored(const Frame &frame, const Layout &layout, Buffer &stored,
                                  bool keep_compressed, Piece piece) const {
    constexpr bool in_pieces = !std::is_same_v<Piece, std::nullptr_t>;
    const size_t n = frame_width_ * frame_height_;
    const size_t size = layout.stored_size;
    const std::vector<hsize_t> start = frame.index >= 0
                                       ? std::vector<hsize_t>{static_cast<hsize_t>(frame.index), 0, 0}
                                       : std::vector<hsize_t>{0, 0};
    const std::vector<hsize_t> count = frame.index >= 0 ? std::vector<hsize_t>{1, frame_height_, frame_width_}
                                                        : std::vector<hsize_t>{frame_height_, frame_width_};
    // Under the lock only where the chunk is (a lookup in its index); the read and the decoding not.
    std::shared_ptr<const HDF5ImageSource::RawFile> raw_file;
    uint64_t address = 0, chunk_size = 0;
    bool read_by_hdf5 = false;
    {
        std::unique_lock ul(hdf5_mutex);
        HDF5DataSet dataset(FileAt(frame.file), frame.dataset);
        unsigned filter_mask = 0;
        if (layout.direct_chunks) {
            haddr_t addr = HADDR_UNDEF;
            hsize_t size = 0;
            if (H5Dget_chunk_info_by_coord(dataset.GetID(), start.data(), &filter_mask, &addr, &size) < 0
                || addr == HADDR_UNDEF || size == 0)
                throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                      "No stored chunk for " + frame.dataset + " of " + files_[frame.file]);
            address = addr;
            chunk_size = size;
        }
        // A chunk stored with one of its filters skipped (its filter mask) is not the encoding decoded
        // below: HDF5 reads it, as it reads a frame that is not stored as one chunk.
        if (!layout.direct_chunks || filter_mask != 0) {
            if (layout.plain_type)
                dataset.ReadVectorToU8(stored, start, count);
            else {
                std::vector<float> v(frame_width_ * frame_height_);
                dataset.ReadVector(v, start, count);
                stored.resize(v.size() * sizeof(float));
                std::memcpy(stored.data(), v.data(), stored.size());
            }
            read_by_hdf5 = true;
        } else {
            raw_file = RawFileAt(frame.file);
            address += user_block_[frame.file];   // HDF5's chunk addresses count from the end of the user block
        }
    }
    if (read_by_hdf5) {
        if constexpr (in_pieces)
            piece(0, stored.data(), n);
        return false;
    }
    if (keep_compressed && !layout.deflated) {
        stored.resize(chunk_size);
        raw_file->ReadAt(stored.data(), stored.size(), address);
        return true;
    }
    thread_local std::vector<uint8_t> compressed, shuffled;
    compressed.resize(chunk_size);
    raw_file->ReadAt(compressed.data(), compressed.size(), address);

    if (layout.deflated && !layout.shuffled && in_pieces) {
        if constexpr (in_pieces)
            InflateInPieces(compressed, n * size, size, frame, piece);
        return false;
    }
    if (layout.deflated) {
        if (!layout.shuffled)
            stored.resize(n * size);
        else
            shuffled.resize(n * size);
        uint8_t *inflated = layout.shuffled ? shuffled.data() : stored.data();
        uLongf inflated_size = n * size;
        if (uncompress(inflated, &inflated_size, compressed.data(), compressed.size()) != Z_OK
            || inflated_size != n * size)
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                  "Cannot inflate " + frame.dataset + " of " + files_[frame.file]);
        if (layout.shuffled) {
            if constexpr (in_pieces) {
                // A piece at a time, as for an unshuffled chunk.
                thread_local std::vector<uint8_t> out(PIECE_BYTES);
                const size_t piece_pixels = PIECE_BYTES / size;
                for (size_t first = 0; first < n; first += piece_pixels) {
                    const size_t count = std::min(piece_pixels, n - first);
                    Unshuffle(shuffled.data(), n, size, first, count, out.data());
                    piece(first, out.data(), count);
                }
                return false;
            }
            stored.resize(n * size);
            Unshuffle(shuffled.data(), n, size, 0, n, stored.data());
        }
    } else {
        stored.resize(n * size);
        JFJochDecompressPtr(stored.data(), layout.algorithm, compressed.data(), compressed.size(), n, size);
    }
    if constexpr (in_pieces)
        piece(0, stored.data(), n);
    return false;
}

template <class Piece>
void JFJochXFELReader::InflateInPieces(const std::vector<uint8_t> &compressed, size_t bytes, size_t size,
                                       const Frame &frame, Piece piece) const {
    thread_local std::vector<uint8_t> out(PIECE_BYTES);
    z_stream z{};
    if (inflateInit(&z) != Z_OK)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot start to inflate");
    z.next_in = const_cast<uint8_t *>(compressed.data());
    z.avail_in = compressed.size();
    size_t done = 0;
    int ret = Z_OK;
    while (ret == Z_OK) {
        z.next_out = out.data();
        z.avail_out = PIECE_BYTES;
        ret = inflate(&z, Z_NO_FLUSH);
        const size_t got = PIECE_BYTES - z.avail_out;
        // A piece is whole pixels: inflate fills it before it returns, except at the end of the frame.
        if ((ret != Z_OK && ret != Z_STREAM_END) || done + got > bytes || got % size != 0)
            break;
        if (got > 0)
            piece(done / size, out.data(), got / size);
        done += got;
    }
    inflateEnd(&z);
    if (ret != Z_STREAM_END || done != bytes)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Cannot inflate " + frame.dataset + " of " + files_[frame.file]);
}

template <class Out, class In, class Buffer, class Convert>
void JFJochXFELReader::PlaceInto(const Frame &frame, const std::vector<uint8_t> &masked, Out err, Buffer &buffer,
                                 Convert convert) const {
    buffer.resize(width_ * height_ * sizeof(Out));
    auto out = reinterpret_cast<Out *>(buffer.data());
    thread_local std::vector<uint8_t> stored;
    if (blocks_.empty()) {   // the stored frame, as stored: converted a piece at a time, as it is read
        ReadStored(frame, layout_, stored, false, [&](size_t first, const uint8_t *bytes, size_t n) {
            const auto *in = reinterpret_cast<const In *>(bytes);
            Out *o = out + first;
            if (masked.empty())
                for (size_t i = 0; i < n; i++)
                    o[i] = static_cast<Out>(convert(in[i]));
            else {
                const uint8_t *m = masked.data() + first;
                for (size_t i = 0; i < n; i++)
                    o[i] = static_cast<Out>(Pick(m[i] != 0, convert(in[i]), err));
            }
        });
        return;
    }
    ReadStored(frame, layout_, stored, false);
    const auto *in = reinterpret_cast<const In *>(stored.data());
    std::fill(out, out + width_ * height_, Out{0});   // between the panels
    for (const auto &b: blocks_)
        for (int64_t ss = 0; ss < b.height; ss++) {
            const size_t src = b.src + static_cast<size_t>(ss) * frame_width_;
            Out *dst = out + b.dst + ss * b.step_ss;
            if (masked.empty() && b.step_fs == 1)
                std::transform(in + src, in + src + b.width, dst, [&](In v) { return static_cast<Out>(convert(v)); });
            else
                for (int64_t fs = 0; fs < b.width; fs++)
                    dst[fs * b.step_fs] = !masked.empty() && masked[src + fs] ? err : static_cast<Out>(convert(in[src + fs]));
        }
}

template <class Buffer>
CompressedImage JFJochXFELReader::DecodeInto(int64_t image_number, Buffer &buffer) const {
    const auto &frame = frames_.at(image_number);
    const CompressedImageMode stored_mode = CalcImageMode(layout_.stored_size, layout_.stored_float,
                                                          layout_.stored_signed);
    // A frame of several panels in signed integers, read in its own units, goes up as stored, with no copy:
    // a bitshuffle/LZ4 chunk still compressed, to be decoded where it is analysed (on the GPU, as
    // Jungfraujoch's own images are), anything else decoded once. Its smallest value is no measurement
    // (below), which the analysis takes it as anyway.
    const bool as_stored = blocks_.empty() && layout_.plain_type && !layout_.stored_float && layout_.stored_signed
                           && adu_per_photon_ == 1.0f;
    if (as_stored && frame_masks_.empty()) {
        const bool compressed = ReadStored(frame, layout_, buffer, true);
        return CompressedImage(buffer.data(), buffer.size(), width_, height_, stored_mode,
                               compressed ? layout_.algorithm : CompressionAlgorithm::NO_COMPRESSION);
    }
    // A pixel a per-frame mask marks in this frame is handed out as no measurement.
    std::vector<uint8_t> masked;
    if (!frame_masks_.empty()) {
        masked.assign(frame_width_ * frame_height_, 0);
        for (size_t i = 0; i < frame_masks_.size(); i++) {
            const auto &m = frame_masks_[i];
            const auto &mask_layout = frame_mask_layouts_[i];
            const Frame mask{frame.file, ReplaceEvent(m.data, frame.group), frame.index, frame.group};
            if (mask_layout.plain_type) {
                thread_local std::vector<uint8_t> mask_stored;
                ReadStored(mask, mask_layout, mask_stored, false, [&](size_t first, const uint8_t *bytes, size_t n) {
                    if (mask_layout.stored_size == 2)
                        MarkBad(m, reinterpret_cast<const uint16_t *>(bytes), first, n, masked);
                    else
                        MarkBad(m, reinterpret_cast<const uint32_t *>(bytes), first, n, masked);
                });
            } else {
                std::vector<uint32_t> v;
                {
                    std::unique_lock ul(hdf5_mutex);
                    v = frame.index >= 0
                            ? FileAt(frame.file).ReadVector<uint32_t>(mask.dataset, {static_cast<hsize_t>(frame.index), 0, 0},
                                                                      {1, frame_height_, frame_width_})
                            : FileAt(frame.file).ReadVector<uint32_t>(mask.dataset, {0, 0}, {frame_height_, frame_width_});
                }
                MarkBad(m, v.data(), 0, v.size(), masked);
            }
        }
    }

    const float gain = adu_per_photon_;
    // A float that is no number (NaN, inf) is no measurement; anything past int32 is clamped to it.
    // Written without a branch or a library call, so that the loop over the pixels is vectorised. The clamp
    // keeps the sign of an infinity; std::lround is spelled out: t is c truncated, c - t is exact, and a
    // fraction of a half or more rounds away from zero.
    const auto photons = [gain](auto v) {
        const float x = static_cast<float>(v) / gain;
        constexpr float limit = 2147483520.0f;
        const float c = std::fabs(x) < limit ? x : std::copysign(limit, x);
        const auto t = static_cast<int32_t>(c);
        const float fraction = c - static_cast<float>(t);
        const int32_t rounded = t + (fraction >= 0.5f) - (fraction <= -0.5f);
        return Pick(!std::isfinite(x), rounded, INT32_MIN);
    };
    const auto as_is = [](auto v) { return v; };
    const auto image = [&](CompressedImageMode mode) {
        return CompressedImage(buffer.data(), buffer.size(), width_, height_, mode, CompressionAlgorithm::NO_COMPRESSION);
    };
    // A signed integer frame's smallest value is no measurement: what a converted file
    // (tools/battery/still/agipd_to_int16.py) writes for a masked pixel or a gap.
    const auto integers = [&](auto to_int32) {
        const auto convert = [to_int32](auto v) -> int32_t {
            if constexpr (std::is_signed_v<decltype(v)>)
                return Pick(v == std::numeric_limits<decltype(v)>::min(), to_int32(v), INT32_MIN);
            return to_int32(v);
        };
        if (layout_.stored_size == 2 && layout_.stored_signed)
            PlaceInto<int32_t, int16_t>(frame, masked, INT32_MIN, buffer, convert);
        else if (layout_.stored_size == 2)
            PlaceInto<int32_t, uint16_t>(frame, masked, INT32_MIN, buffer, convert);
        else if (layout_.stored_signed)
            PlaceInto<int32_t, int32_t>(frame, masked, INT32_MIN, buffer, convert);
        else
            PlaceInto<int32_t, uint32_t>(frame, masked, INT32_MIN, buffer, convert);
    };

    // Floats become photons, rounded. Integers stay integers: divided only by a gain other than one, and a
    // 16-bit signed frame - as Cheetah writes - is handed out in its own type.
    if (layout_.stored_float || !layout_.plain_type)
        PlaceInto<int32_t, float>(frame, masked, INT32_MIN, buffer, photons);
    else if (gain != 1.0f)
        integers(photons);
    else if (layout_.stored_size == 2 && layout_.stored_signed) {
        PlaceInto<int16_t, int16_t>(frame, masked, INT16_MIN, buffer, as_is);
        return image(CompressedImageMode::Int16);
    } else
        integers([](auto v) { return static_cast<int32_t>(v); });
    return image(CompressedImageMode::Int32);
}

bool JFJochXFELReader::LoadImage_i(std::shared_ptr<JFJochReaderDataset> &dataset,
                                       DataMessage &message,
                                       std::vector<uint8_t> &buffer,
                                       int64_t image_number,
                                       bool update_dataset) {
    (void) update_dataset;
    if (!dataset || image_number < 0 || static_cast<size_t>(image_number) >= frames_.size())
        return false;
    message.image = DecodeInto(image_number, buffer);
    message.number = image_number;
    return true;
}

bool JFJochXFELReader::ReadRawImage(int64_t image_number, JFJochReaderRawImage &image) {
    if (image_number < 0 || static_cast<size_t>(image_number) >= frames_.size())
        return false;
    image.image = DecodeInto(image_number, image.image_buffer);
    return true;
}

std::vector<SpotToSave> JFJochXFELReader::ReadSpots(int64_t) const {
    return {};
}