Jungfraujoch/image_analysis/RotationIndexer.cpp

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

#include "RotationIndexer.h"

#include "geom_refinement/XtalOptimizer.h"
#include "indexing/FFTIndexer.h"
#include "lattice_search/LatticeSearch.h"
#include <iostream>

RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer)
    : experiment(x),
      index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()),
      v_(experiment.GetImageNum()),
      axis_(x.GetGoniometer()),
      geom_(x.GetDiffractionGeometry()),
      updated_geom_(geom_),
      indexer_(indexer) {
    if (axis_) {
        angle_norm_deg = std::fabs(axis_->GetIncrement_deg());
        if (angle_norm_deg < 1e-6) {
            // Guard against rotation close to zero
            axis_ = std::nullopt;
        } else {
            if (x.GetImageNum() < min_images_for_indexing) {
                // For short measurements - only indexing at the end
                first_image_to_try_indexing = std::max<int64_t>(0, x.GetImageNum() - 1);
                next_image_to_try_indexing = first_image_to_try_indexing;
                image_stride = 1;
            } else {
                first_image_to_try_indexing = std::max<int64_t>(min_images_for_indexing,
                    x.GetIndexingSettings().GetRotationIndexingMinAngularRange_deg() / angle_norm_deg);
                next_image_to_try_indexing = first_image_to_try_indexing;
                image_stride = std::ceil(x.GetIndexingSettings().GetRotationIndexingAngularStride_deg() / angle_norm_deg);
                if (image_stride == 0)
                    image_stride = 1;
            }
        }
    }
}

IndexerResult RotationIndexer::RunIndexing() const {
    std::vector<Coord> coords;
    coords.reserve(max_spots);

    for (int i = 0; i < v_.size(); i++) {
        const float angle_deg = axis_->GetAngle_deg(i) + axis_->GetWedge_deg() / 2.0f;
        const auto rot = axis_->GetTransformationAngle(angle_deg);

        for (const auto &s: v_[i])
            coords.emplace_back(rot * s.ReciprocalCoord(geom_));
    }

    return indexer_.Run(experiment, coords);
}

void RotationIndexer::TryIndex() {
    auto indexer_result = RunIndexing();

    if (!indexer_result.lattice.empty() && indexer_result.lattice[0].CalcVolume() > 1.0) {
        auto sg = experiment.GetGemmiSpaceGroup();
        if (sg) {
            search_result_ = LatticeSearchResult{
                .niggli_class = 0, // Since Niggli class was not searched for, we don't know which one
                .conventional = indexer_result.lattice[0], // If lattice provided, it is for now primitive == conventional
                .system = sg->crystal_system(),
                .centering = sg->centring_type(),
            };
        } else {
            // Find lattice type based on cell
            search_result_ = LatticeSearch(indexer_result.lattice[0]);
        }

        // Run refinement
        DiffractionExperiment experiment_copy(experiment);
        XtalOptimizerData data{
            .geom = experiment_copy.GetDiffractionGeometry(),
            .latt = search_result_.conventional,
            .crystal_system = search_result_.system,
            .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
            .refine_beam_center = true,
            .refine_distance_mm =  false,
            .refine_detector_angles = true,
            .refine_rotation_axis = true,
            .index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(),
            .axis = axis_
        };

        if (data.crystal_system == gemmi::CrystalSystem::Trigonal)
            data.crystal_system = gemmi::CrystalSystem::Hexagonal;

        if (data.crystal_system == gemmi::CrystalSystem::Monoclinic)
            data.latt.ReorderMonoclinic();

        if (XtalOptimizer(data, v_)) {
            indexed_lattice = data.latt;
            updated_geom_ = data.geom;
            axis_ = data.axis;
            return;
        }
    }

    if (indexing_range_multiplier < max_indexing_range_multiplier) {
        indexing_range_multiplier++;
        next_image_to_try_indexing = first_image_to_try_indexing * indexing_range_multiplier;
    } else {
        next_image_to_try_indexing = INT64_MAX;
    }
}

void RotationIndexer::ProcessImage(int64_t image, const std::vector<SpotToSave> &spots) {
    std::unique_lock ul(m);

    // For non-rotation just ignore the whole procedure
    if (!axis_)
        return;

    if (!indexed_lattice && image >= last_accumulated_image + image_stride) {
        v_[image].reserve(spots.size());

        for (const auto &s: spots) {
            if (index_ice_rings || !s.ice_ring)
                v_[image].emplace_back(s);
        }

        accumulated_images++;
        last_accumulated_image = image;

        const bool short_scan_last_image =
            (experiment.GetImageNum() < min_images_for_indexing) &&
            (image >= experiment.GetImageNum() - 1);

        if ((accumulated_images >= min_images_for_indexing || short_scan_last_image) &&
            image >= next_image_to_try_indexing)
            TryIndex();
    }
}

std::optional<RotationIndexerResult> RotationIndexer::GetLattice() const {
    std::unique_lock ul(m);

    if (!indexed_lattice)
        return {};
    return RotationIndexerResult{
        .lattice = indexed_lattice.value(),
        .search_result = search_result_,
        .geom = updated_geom_,
        .axis = axis_
    };
}