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 <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. Only where the panel is: its pixel range, axes,
// corner, pixel size, depth offset and camera length.
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
        return false;
    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;
    std::map<std::string, size_t> panel_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) && !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")) {
            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;
        }
        if (!SetPanel(filename, geom.panels[panel_index[name]], 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)");
    }
    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