Jungfraujoch/image_analysis/indexing/Indexer.cpp

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

#include "Indexer.h"

void Indexer::Setup(const DiffractionExperiment& experiment) {
    geom = experiment.GetDiffractionGeometry();
    indexing_tolerance = experiment.GetIndexingSettings().GetTolerance();
    dist_tolerance_vs_reference = experiment.GetIndexingSettings().GetUnitCellDistTolerance();
    step_down_to_reference = experiment.GetIndexingSettings().GetRotationIndexing();
    search_short_rows = experiment.IsRotationIndexing();
    viable_cell_min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots();
    index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings();
    wavelength_A = experiment.GetWavelength_A();
    spindle_axis = {};
    if (const auto goniometer = experiment.GetGoniometer())
        spindle_axis = goniometer->GetAxis();
    spindle_theta_max_deg = SpindleThetaMax_deg(wavelength_A, experiment.GetDetectorMaxResolution_A());
    centering = experiment.GetCentering();   // before SetupUnitCell, which reads it
    SetupUnitCell(experiment.GetUnitCell());
}

IndexerResult Indexer::Run(const std::vector<Coord> &coord, bool severity_only) {
    IndexerResult ret;
    auto start = std::chrono::steady_clock::now();
    std::optional<SpindleSeverity> severity;
    if (severity_only) {
        severity = RunSeverityOnly(coord);
    } else {
        ret.lattice = RunInternal(coord, coord.size());
        ret.executed = true;
        severity = GetSpindleSeverity();
    }

    auto end = std::chrono::steady_clock::now();
    std::chrono::duration<float> duration = end - start;
    ret.indexing_time_s = duration.count();
    if (severity) {
        ret.spindle_blind_fraction = severity->score;
        ret.spindle_row_length_A = severity->row_length_A;
        ret.spindle_miss_angle_deg = severity->miss_angle_deg;
    }
    return ret;
}

namespace {
    Eigen::Matrix3f Rows(const CrystalLattice &l) {
        Eigen::Matrix3f m;
        m << l.Vec0().x, l.Vec0().y, l.Vec0().z, l.Vec1().x, l.Vec1().y, l.Vec1().z,
                l.Vec2().x, l.Vec2().y, l.Vec2().z;
        return m;
    }

    CrystalLattice FromRows(const Eigen::Matrix3f &m) {
        return {Coord(m(0, 0), m(0, 1), m(0, 2)), Coord(m(1, 0), m(1, 1), m(1, 2)), Coord(m(2, 0), m(2, 1), m(2, 2))};
    }
}

CrystalLattice Indexer::IndexingCell(const UnitCell &cell) {
    const CrystalLattice conventional(cell);
    to_conventional = Eigen::Matrix3f::Identity();
    if (centering == 'P')
        return conventional;
    Eigen::Matrix3f reduced = Rows(conventional.ToPrimitive(centering).NiggliReduce());
    if (reduced.determinant() < 0.0f)
        reduced = -reduced;   // right-handed, as Refine returns its cells
    to_conventional = (Rows(conventional) * reduced.inverse()).array().round().matrix();
    return FromRows(reduced);
}

std::vector<CrystalLattice> Indexer::ToConventional(const std::vector<CrystalLattice> &lattices) const {
    std::vector<CrystalLattice> ret;
    ret.reserve(lattices.size());
    for (const auto &l: lattices)
        ret.push_back(FromRows(to_conventional * Rows(l)));
    return ret;
}