Jungfraujoch/reader/ImgCIF.h

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

#pragma once

#include <array>
#include <cstdint>
#include <map>
#include <optional>
#include <string>
#include <vector>

#include "../common/Coord.h"
#include "../common/DetectorOrientation.h"
#include "../common/GoniometerAxis.h"

// The parts of a full imgCIF header (the form CBFlib writes, and Bruker's APEX writes when it converts
// .sfrm frames) that say where the detector and the goniometer stand. Unlike a PILATUS miniCBF, which
// states its geometry in "# " comment lines, such a file states it only in CIF tables: an axis table
// (_axis), the setting of every axis at this frame (_diffrn_scan_frame_axis) and how it moves per
// frame (_diffrn_scan_axis), and which axes the image's two indices run along
// (_array_structure_list, _array_structure_list_axis).
//
// Every vector here is in the imgCIF laboratory frame (Z from the sample towards the source, Y
// opposite gravity, X completing a right-handed set) and every length in millimetres.
// Following Hammersley, Bernstein & Westbrook (2006) Int. Tables Cryst. G, 444-458
namespace imgcif {

// One "loop_" of a CIF block, as rows keyed by tag. Tags and values are both read as
// whitespace-separated tokens rather than by line, because a template packs several tags onto one
// line ("_axis.vector[1] _axis.vector[2] _axis.vector[3]") and the rows that follow are laid out to
// match the tags, not the lines.
std::vector<std::map<std::string, std::string>> ParseLoop(const std::string &text, const std::string &tag);

// A tag stated outside a loop ("_diffrn_radiation_wavelength.wavelength 0.71073"), or nothing.
std::optional<std::string> Item(const std::string &text, const std::string &tag);

struct Axis {
    std::string type;                   // "rotation" or "translation"
    std::string equipment;              // "goniometer", "detector", ...
    std::string depends_on;             // "." for an axis that hangs off nothing
    std::array<double, 3> vector{};
    std::array<double, 3> offset{};
    double setting = 0;                 // this frame's angle (deg) or displacement (mm)
    double increment = 0;               // its change over one frame
};

// One index of the stored image: the axis it runs along, and where pixel 1 and each next one sit on it.
struct ArrayIndex {
    std::string axis_id;
    int64_t dimension = 0;
    double displacement = 0;            // of index 1, mm
    double increment = 0;               // per index, mm
};

struct Description {
    std::map<std::string, Axis> axes;
    ArrayIndex fast;                    // _array_structure_list.index 1
    ArrayIndex slow;                    // _array_structure_list.index 2
    double wavelength_A = 0;
    double overload = 0;                // _array_intensities.overload; 0 where not stated
    double exposure_s = 0;
    std::string detector;
    // _diffrn_detector.detector names a CCD: its pixel values are ADU, not counted photons.
    bool ccd = false;
    // _diffrn_radiation_wavelength.polarizn_source_ratio, (Ip - In) / (Ip + In), and
    // polarizn_source_norm, the angle (deg) from the laboratory Y axis to the normal of the plane Ip
    // is polarized in. Nothing where the header states no ratio.
    std::optional<double> polarization_ratio;
    double polarization_norm_deg = 0;
};

// The polarization as DiffractionExperiment::PolarizationFactor takes it - the degree, positive for a
// beam polarized in the horizontal (laboratory X-Z) plane - or nothing where the header states none.
// A plane at an angle psi from the horizontal enters as its projection, the degree times cos(2 psi),
// as the correction knows only those two planes.
std::optional<float> PolarizationFactor(const Description &d);

// The description, or nothing when the header is not a full imgCIF one - no axis table, or no
// statement of where the axes stand. Throws on a table that is there but cannot be read.
std::optional<Description> Parse(const std::string &text);

struct Geometry {
    std::array<double, 3> origin{};     // centre of stored pixel (0, 0)
    std::array<double, 3> fast{};       // one stored column further, minus origin
    std::array<double, 3> slow{};       // one stored row further, minus origin
    std::string scan_axis;              // the goniometer axis that turns during the frame
    std::array<double, 3> scan_direction{};   // its laboratory direction at this frame's settings
    double scan_start_deg = 0;
    double scan_increment_deg = 0;
};

// Where the image and the rotation axis are at this frame, through the whole axis chain. Throws, naming
// what it found, on anything this cannot describe as one rotation about one axis: two axes moving at
// once, a detector or translation axis moving, or a chain that does not close.
Geometry Evaluate(const Description &d);

// An imgCIF laboratory direction in the internal frame. imgCIF puts Z from the sample towards the
// source and Y opposite gravity, while the internal frame has z along the beam and y along increasing
// row, so the two differ by a half turn about x - a rotation and not a mirror, so an axis carried
// through it turns the same way by the same angle.
Coord ToInternal(const std::array<double, 3> &v);

// The goniometer axis as GoniometerAxis takes it, for an imgCIF rotation axis. GoniometerAxis turns
// the crystal by minus the angle about its axis (IndexAndRefine, rot_to_image), the sense an NXmx
// master's vector carries, so a right-handed imgCIF rotation about v is the internal -v.
Coord GoniometerAxisOf(const std::array<double, 3> &v);

// The sweep as GoniometerAxis takes it, turning the scanned axis v from start_deg by increment_deg per
// frame. A sweep that turns backwards (Bruker frames state a negative increment) is the same rotation
// about -v by the opposite angles, and is stated that way: the analysis takes the oscillation width of a
// frame to be positive, as every other format states it.
GoniometerAxis SweepAxis(const std::string &name, double start_deg, double increment_deg,
                         const std::array<double, 3> &v);

// The detector plane in the terms the geometry is kept in: the point of normal incidence (PONI) in
// pixels, the distance along the normal, and the rotation that takes the detector frame to the
// laboratory, split into the discrete way the image is stored and the continuous rest. The table
// states the plane through three pixel centres, so this is only a change of description - a 2theta
// arm, a detector's small mounting tilts and a stored image turned or mirrored all come out of the same
// three points. Throws, naming path, on pixels that are not square or a grid that is not rectangular.
struct DetectorPlane {
    float pixel_mm = 0;
    float distance_mm = 0;
    float beam_x_pxl = 0;
    float beam_y_pxl = 0;
    DetectorOrientation orientation;
    float rot1 = 0, rot2 = 0, rot3 = 0;   // PONI angles, rad
};
DetectorPlane PlaneOf(const Geometry &g, const std::string &path);

}  // namespace imgcif