Jungfraujoch/reader/CrystFELGeometry.cpp

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

#include "CrystFELGeometry.h"

#include <algorithm>
#include <array>
#include <filesystem>
#include <fstream>
#include <map>
#include <regex>
#include <set>
#include <sstream>

#include "../common/JFJochException.h"

namespace crystfel {

using xfel::Value;
using xfel::Panel;

namespace {

std::string Trim(const std::string &s) {
    const auto b = s.find_first_not_of(" \t\r");
    if (b == std::string::npos)
        return "";
    return s.substr(b, s.find_last_not_of(" \t\r") - b + 1);
}

[[noreturn]] void Refuse(const std::string &filename, const std::string &key, const std::string &why) {
    throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                          filename + ": '" + key + "' " + why + "; this geometry is not read");
}

// "0.1 m", "70.5 mm", "/LCLS/detector_1/EncoderValue", "9.3 keV": a number or a path, with an optional
// unit after it. units maps each accepted unit to its factor into the unit held, and "" to the default.
Value ParseValue(const std::string &filename, const std::string &key, const std::string &text,
                 const std::map<std::string, double> &units, double path_default) {
    std::istringstream ss(text);
    std::string token, unit, rest;
    ss >> token >> unit >> rest;
    if (token.empty() || !rest.empty())
        Refuse(filename, key, "has a value that is not a number or a path with one unit (" + text + ")");
    const bool is_path = token[0] == '/';
    double scale = is_path ? path_default : units.at("");
    if (!unit.empty()) {
        const auto it = units.find(unit);
        if (it == units.end() || unit.empty())
            Refuse(filename, key, "is stated in a unit not understood (" + unit + ")");
        scale = it->second;
    }
    Value v;
    v.scale = scale;
    if (is_path)
        v.path = token;
    else {
        size_t used = 0;
        double number = 0;
        try {
            number = std::stod(token, &used);
        } catch (const std::exception &) {
            used = 0;
        }
        if (used != token.size())
            Refuse(filename, key, "has a value that is not a number or a path (" + text + ")");
        v.value = number * scale;
        v.scale = 1;
    }
    return v;
}

const std::map<std::string, double> LENGTH_UNITS = {{"", 1.0}, {"m", 1.0}, {"mm", 1e-3}};

// "+0.000390x -1.000000y", "-x", "0.5x+0.5y": the in-plane direction of a panel axis. A panel tilted
// out of the detector plane (a z component) is not read: it would be placed flat.
void ParseDirection(const std::string &filename, const std::string &key, const std::string &txt,
                    double &x, double &y) {
    x = 0;
    y = 0;
    static const std::regex term(R"(([+-]?)\s*([0-9.]*(?:[eE][+-]?[0-9]+)?)\s*([xyz]))");
    for (auto it = std::sregex_iterator(txt.begin(), txt.end(), term); it != std::sregex_iterator(); ++it) {
        double v = (*it)[2].str().empty() ? 1.0 : std::stod((*it)[2].str());
        if ((*it)[1] == "-")
            v = -v;
        if ((*it)[3] == "x")
            x = v;
        else if ((*it)[3] == "y")
            y = v;
        else if (v != 0)
            Refuse(filename, key, "tilts the panel out of the detector plane (" + txt + ")");
    }
}

// Keywords that only matter to refinement or peak search, read past without a word.
bool Bookkeeping(const std::string &key) {
    return key.starts_with("rigid_group") || key.starts_with("group_") || key.starts_with("peak_list")
           || key == "badrow_direction" || key == "peak_sep";
}

// Keywords that change where a pixel lands or what it is worth, which this reader does not do.
const std::set<std::string> REFUSED = {"detector_shift_x", "detector_shift_y"};

// Sets one panel property; false if the key is not one. Where the panel is - its pixel range, axes,
// corner, pixel size, depth offset and camera length - and which of its pixels are left out.
bool SetPanel(const std::string &filename, Panel &p, const std::string &key, const std::string &value) {
    if (key == "min_fs") p.min_fs = std::stoll(value);
    else if (key == "max_fs") p.max_fs = std::stoll(value);
    else if (key == "min_ss") p.min_ss = std::stoll(value);
    else if (key == "max_ss") p.max_ss = std::stoll(value);
    else if (key == "fs") ParseDirection(filename, key, value, p.fs_x, p.fs_y);
    else if (key == "ss") ParseDirection(filename, key, value, p.ss_x, p.ss_y);
    else if (key == "corner_x") p.corner_x = std::stod(value);
    else if (key == "corner_y") p.corner_y = std::stod(value);
    else if (key == "res") p.res = std::stod(value);
    else if (key == "clen") p.clen = ParseValue(filename, key, value, LENGTH_UNITS, 1e-3);
    else if (key == "coffset") p.coffset = std::stod(value);
    else if (key == "mask_edge_pixels") p.mask_edge_pixels = std::stoll(value);
    else if (key == "no_index") p.no_index = value == "1" || value == "true";
    else
        return false;
    return true;
}

// The up to eight masks a panel has (maskN_data, maskN_file, maskN_goodbits, maskN_badbits, N = 0..7;
// mask, mask_file, mask_good and mask_bad are mask 0's). Sets the one the key names; false if the key
// is not a mask key.
using Masks = std::array<xfel::Mask, 8>;
bool SetMask(Masks &masks, const std::string &key, const std::string &value) {
    static const std::map<std::string, std::string> OLD = {
        {"mask", "data"}, {"mask_file", "file"}, {"mask_good", "goodbits"}, {"mask_bad", "badbits"}};
    static const std::regex CURRENT(R"(mask([0-7])_(data|file|goodbits|badbits))");
    size_t n = 0;
    std::string field;
    std::smatch match;
    if (OLD.contains(key))
        field = OLD.at(key);
    else if (std::regex_match(key, match, CURRENT)) {
        n = std::stoul(match[1].str());
        field = match[2].str();
    } else
        return false;
    if (field == "data") masks[n].data = value;
    else if (field == "file") masks[n].file = value;
    else if (field == "goodbits") masks[n].good = static_cast<uint32_t>(std::stoul(value, nullptr, 0));
    else masks[n].bad = static_cast<uint32_t>(std::stoul(value, nullptr, 0));
    return true;
}

// Sets one property of a bad region: a box in the laboratory plane (min_x .. max_y, in pixels) or a
// pixel range in the stored frame (min_fs .. max_ss, inclusive) on one panel - or, in the old syntax,
// with no panel, on the whole frame. False if the key is not one.
bool SetBad(xfel::BadRegion &b, const std::string &key, const std::string &value) {
    if (key == "min_x") b.min_x = std::stod(value);
    else if (key == "max_x") b.max_x = std::stod(value);
    else if (key == "min_y") b.min_y = std::stod(value);
    else if (key == "max_y") b.max_y = std::stod(value);
    else if (key == "min_fs") b.min_fs = std::stoll(value);
    else if (key == "max_fs") b.max_fs = std::stoll(value);
    else if (key == "min_ss") b.min_ss = std::stoll(value);
    else if (key == "max_ss") b.max_ss = std::stoll(value);
    else if (key == "panel") b.panel = value;
    else
        return false;
    if (key.ends_with("_x") || key.ends_with("_y"))
        b.in_lab = true;
    return true;
}

}  // namespace

xfel::Layout Read(const std::string &filename, std::vector<std::string> &ignored) {
    std::ifstream in(filename);
    if (!in)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Cannot open " + filename);

    xfel::Layout geom;
    // Global panel keys hold for every panel that appears after them, as in the old syntax; a file in
    // the current syntax states them first, so this reads both the same.
    Panel defaults;
    Masks default_masks;
    std::vector<Masks> panel_masks;    // each panel's, beside geom.panels
    std::map<std::string, size_t> panel_index, bad_index;
    std::map<std::string, std::string> dims;

    std::string line;
    while (std::getline(in, line)) {
        line = line.substr(0, line.find(';'));
        const auto eq = line.find('=');
        if (eq == std::string::npos)
            continue;
        const std::string key = Trim(line.substr(0, eq));
        const std::string value = Trim(line.substr(eq + 1));
        const auto slash = key.find('/');

        if (REFUSED.contains(key.substr(slash == std::string::npos ? 0 : slash + 1)))
            Refuse(filename, key, "is not supported");

        if (slash == std::string::npos) {
            if (key == "data") geom.data = value;
            else if (key.starts_with("dim")) dims[key] = value;
            else if (!SetPanel(filename, defaults, key, value) && !SetMask(default_masks, key, value)
                     && !Bookkeeping(key))
                ignored.push_back(key);
            continue;
        }

        const std::string name = key.substr(0, slash);
        const std::string prop = key.substr(slash + 1);

        if (name.starts_with("bad")) {
            if (!bad_index.contains(name)) {
                bad_index[name] = geom.bad.size();
                geom.bad.emplace_back();
            }
            if (!SetBad(geom.bad[bad_index[name]], prop, value))
                ignored.push_back(key);
            continue;
        }

        if (prop == "data") {
            if (value != geom.data)
                Refuse(filename, key, "puts a panel in its own dataset");
            continue;
        }
        if (prop.starts_with("dim")) {
            dims[prop] = value;
            continue;
        }
        if (!panel_index.contains(name)) {
            panel_index[name] = geom.panels.size();
            geom.panels.push_back(defaults);
            geom.panels.back().name = name;
            panel_masks.push_back(default_masks);
        }
        const size_t n = panel_index[name];
        if (!SetPanel(filename, geom.panels[n], prop, value) && !SetMask(panel_masks[n], prop, value)
            && !Bookkeeping(prop))
            ignored.push_back(key);
    }

    // The frame layouts read: a stack of whole frames (%, ss, fs) or one frame (ss, fs). A dimension
    // naming a module number lays the panels out in a way the stored frame here cannot express.
    for (const auto &[key, value]: dims) {
        const bool ok = (key == "dim0" && (value == "%" || (dims.size() == 2 && value == "ss")))
                        || (key == "dim1" && (value == "ss" || (dims.size() == 2 && value == "fs")))
                        || (key == "dim2" && value == "fs" && dims.size() == 3);
        if (!ok)
            Refuse(filename, key + " = " + value, "lays the frame out in a way not supported");
    }

    if (geom.panels.empty())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, filename + " defines no panels");
    for (const auto &p: geom.panels) {
        if (!(p.res > 0))
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                  filename + ": panel " + p.name + " has no pixel size (res)");
    }

    // Each mask once for the whole frame where every panel has it alike (the usual case: stated once,
    // before the panels), otherwise once for each panel that has it, on that panel's pixels only.
    for (size_t n = 0; n < default_masks.size(); n++) {
        const auto &first = panel_masks[0][n];
        const bool alike = std::all_of(panel_masks.begin(), panel_masks.end(),
                                       [&](const Masks &m) { return m[n] == first; });
        if (alike) {
            if (!first.data.empty())
                geom.masks.push_back(first);
            continue;
        }
        for (size_t i = 0; i < geom.panels.size(); i++) {
            if (panel_masks[i][n].data.empty())
                continue;
            geom.masks.push_back(panel_masks[i][n]);
            geom.masks.back().panel = geom.panels[i].name;
        }
    }
    // A mask file is named as it is, or relative to the working directory, as CrystFEL reads it; a
    // relative one not found there is looked for beside the geometry file.
    for (auto &m: geom.masks) {
        const std::filesystem::path file(m.file);
        const auto beside = std::filesystem::path(filename).parent_path() / file;
        if (!m.file.empty() && file.is_relative() && !std::filesystem::exists(file) && std::filesystem::exists(beside))
            m.file = beside.string();
    }
    return geom;
}

DetectorModules ModulesOnImage(const xfel::Layout &geom, double beam_x, double beam_y) {
    if (geom.panels.size() < 2)
        return {};
    std::vector<DetectorModule> out;
    for (const auto &p: geom.panels) {
        DetectorModule m;
        m.x0 = p.min_fs;
        m.y0 = p.min_ss;
        m.width = p.max_fs - p.min_fs + 1;
        m.height = p.max_ss - p.min_ss + 1;
        // The panel's first pixel has its centre half a pixel along both axes from the corner.
        m.ox = beam_x + p.corner_x + 0.5 * (p.fs_x + p.ss_x);
        m.oy = beam_y + p.corner_y + 0.5 * (p.fs_y + p.ss_y);
        m.fx = p.fs_x;
        m.fy = p.fs_y;
        m.sx = p.ss_x;
        m.sy = p.ss_y;
        out.push_back(m);
    }
    return DetectorModules(out);
}

}  // namespace crystfel