Jungfraujoch/tests/SpindleBlindFractionTest.cpp

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

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>

#include "../image_analysis/indexing/SpindleBlindFraction.h"

using Catch::Matchers::WithinAbs;

TEST_CASE("SpindleBlindFraction_Overlap", "[Indexing][Spindle]") {
    // x is the FOLDED miss-angle over theta_max: 0 both on the spindle and perpendicular to it.
    // A row at either end leaves the whole cone unrecoverable; one on the cone edge leaves none.
    CHECK_THAT(BlindConeSelfOverlap(0.0f), WithinAbs(1.0f, 1e-6));
    CHECK_THAT(BlindConeSelfOverlap(1.0f), WithinAbs(0.0f, 1e-6));
    CHECK_THAT(BlindConeSelfOverlap(2.0f), WithinAbs(0.0f, 1e-6));
    // Half way into the cone the two caps still share 39% of their area.
    CHECK_THAT(BlindConeSelfOverlap(0.5f), WithinAbs(0.3910f, 1e-3));
    // Monotone decreasing
    for (int i = 0; i < 20; i++)
        CHECK(BlindConeSelfOverlap(i / 20.0f) >= BlindConeSelfOverlap((i + 1) / 20.0f));
}

TEST_CASE("SpindleBlindFraction_Rows", "[Indexing][Spindle]") {
    const Coord spindle(0, 0, 1);
    const float theta_max = 20.0f;

    SECTION("a short row on the spindle is the worst case") {
        const std::vector<Coord> rows = {Coord(0, 0, 50), Coord(60, 0, 0), Coord(0, 70, 0)};
        const std::vector<float> mag = {100, 90, 80};
        const auto s = SpindleBlindFraction(rows, mag, spindle, theta_max);
        REQUIRE(s.has_value());
        CHECK_THAT(s->score, WithinAbs(1.0f, 1e-5));
        CHECK_THAT(s->miss_angle_deg, WithinAbs(0.0f, 1e-3));
        CHECK_THAT(s->row_length_A, WithinAbs(50.0f, 1e-3));
    }

    SECTION("a row perpendicular to the spindle is the worst case too") {
        // The 2-fold about it carries the blind cone onto the cone's opposite lobe, which the sweep
        // leaves equally unmeasured. Reporting this as harmless was the bug the fold fixes.
        const std::vector<Coord> rows = {Coord(50, 0, 0), Coord(0, 60, 0), Coord(0, 0, 70)};
        const std::vector<float> mag = {100, 90, 80};
        const auto s = SpindleBlindFraction(rows, mag, spindle, theta_max);
        REQUIRE(s.has_value());
        CHECK_THAT(s->score, WithinAbs(1.0f, 1e-5));
        CHECK_THAT(s->miss_angle_deg, WithinAbs(90.0f, 1e-3));
    }

    SECTION("no short row near either end of the range scores zero") {
        // Rows well away from both the spindle and its perpendicular plane: any 2-fold about them
        // swings the cone clear of itself, and one sweep loses nothing symmetry could have returned.
        const std::vector<Coord> rows = {Coord(50, 0, 50), Coord(0, 60, 60), Coord(40, 40, 56)};
        const std::vector<float> mag = {100, 90, 80};
        const auto s = SpindleBlindFraction(rows, mag, spindle, theta_max);
        REQUIRE(s.has_value());
        CHECK_THAT(s->score, WithinAbs(0.0f, 1e-6));
    }

    SECTION("a row too long to be a symmetry axis is ignored") {
        // 300 A along the spindle, in a crystal whose own rows are 70-85 A: 4x the shortest row is
        // not a plausible symmetry axis, and the cone it sits in is not the crystal's problem.
        const std::vector<Coord> rows = {Coord(50, 0, 50), Coord(0, 60, 60), Coord(0, 0, 300)};
        const std::vector<float> mag = {100, 90, 80};
        const auto s = SpindleBlindFraction(rows, mag, spindle, theta_max);
        REQUIRE(s.has_value());
        CHECK_THAT(s->score, WithinAbs(0.0f, 1e-6));
    }

    SECTION("the same row in a crystal that IS that big is not ignored") {
        const std::vector<Coord> rows = {Coord(180, 0, 180), Coord(0, 200, 200), Coord(0, 0, 300)};
        const std::vector<float> mag = {100, 90, 80};
        const auto s = SpindleBlindFraction(rows, mag, spindle, theta_max);
        REQUIRE(s.has_value());
        CHECK_THAT(s->score, WithinAbs(1.0f, 1e-5));
        CHECK_THAT(s->miss_angle_deg, WithinAbs(0.0f, 1e-3));
    }

    SECTION("a weak spurious short row does not shrink the length window") {
        // What a long-cell still produces: the real rows near 300 A plus a weaker short peak. Taking
        // the window off that peak would hide the aligned row and report the orientation harmless.
        const std::vector<Coord> rows = {Coord(70, 70, 30), Coord(180, 0, 180), Coord(0, 0, 300)};
        const std::vector<float> mag = {40, 100, 95};
        const auto s = SpindleBlindFraction(rows, mag, spindle, theta_max);
        REQUIRE(s.has_value());
        CHECK_THAT(s->score, WithinAbs(1.0f, 1e-5));
    }

    SECTION("a lone 2-fold on an axis too long to see is recovered from the visible rows' normal") {
        // A monoclinic-like crystal: the unique axis is far beyond the length window, so no
        // shortlist row points along it, but every visible row is perpendicular to it, and the
        // normal of any two of them is its direction - a symmetry axis is parallel in the direct
        // and reciprocal bases. Here that direction is perpendicular to the spindle, the case the
        // fold exists for; before the pair-normal search this scored 0, a silent "safe".
        const std::vector<Coord> rows = {Coord(0, 45, 45), Coord(0, 60, 25)};
        const std::vector<float> mag = {100, 90};
        const auto s = SpindleBlindFraction(rows, mag, spindle, theta_max);
        REQUIRE(s.has_value());
        CHECK_THAT(s->score, WithinAbs(1.0f, 1e-5));
        CHECK_THAT(s->miss_angle_deg, WithinAbs(90.0f, 1e-3));
        // 0 marks a direction inferred from a pair of rows rather than a measured row.
        CHECK_THAT(s->row_length_A, WithinAbs(0.0f, 1e-6));
    }

    SECTION("a shortlist the pass could not resolve gives no answer at all") {
        // Strong rows ten times longer than the shortest entry: the grid has lost the crystal's real
        // rows and is returning spurious short ones. Reporting zero here would be a silent "safe".
        const std::vector<Coord> rows = {Coord(12, 5, 0), Coord(250, 0, 0), Coord(0, 0, 300)};
        const std::vector<float> mag = {40, 100, 95};
        CHECK_FALSE(SpindleBlindFraction(rows, mag, spindle, theta_max).has_value());
    }

    SECTION("no rows, no answer") {
        CHECK_FALSE(SpindleBlindFraction({}, {}, spindle, theta_max).has_value());
        CHECK_FALSE(SpindleBlindFraction({Coord(0, 0, 50)}, {1.0f}, Coord(0, 0, 0), theta_max).has_value());
    }

    SECTION("a wider cone at long wavelength makes the same miss-angle worse") {
        const std::vector<Coord> rows = {Coord(0, 20, 50), Coord(50, 0, 50)};
        const std::vector<float> mag = {100, 90};
        const auto narrow = SpindleBlindFraction(rows, mag, spindle, 10.0f);
        const auto wide = SpindleBlindFraction(rows, mag, spindle, 35.0f);
        REQUIRE(narrow.has_value());
        REQUIRE(wide.has_value());
        CHECK(wide->score > narrow->score);
    }
}

TEST_CASE("SpindleBlindFraction_FromLattice", "[Indexing][Spindle]") {
    const Coord spindle(0, 0, 1);
    const float theta_max = 20.0f;

    SECTION("a known-cell frame answers from its lattice rows") {
        // Monoclinic-like basis with the unique axis along the spindle. The severity needs no FFT
        // shortlist: the lattice's own short rows carry the answer, here a worst case twice over
        // (a row on the spindle and rows perpendicular to it).
        const CrystalLattice latt(Coord(50, 0, 0), Coord(0, 0, 60), Coord(20, 70, 0));
        const auto s = SpindleBlindFractionFromLattice(latt, spindle, theta_max);
        REQUIRE(s.has_value());
        CHECK_THAT(s->score, WithinAbs(1.0f, 1e-5));
    }

    SECTION("a long unique axis outside the window is recovered from the pair normals") {
        // The 300 A axis is excluded by the length window as a row, but every kept row is
        // perpendicular to it, so the normal of any pair recovers its direction - perpendicular
        // to the spindle, the lone-diad worst case.
        const CrystalLattice latt(Coord(0, 45, 45), Coord(300, 0, 0), Coord(0, -60, 25));
        const auto s = SpindleBlindFractionFromLattice(latt, spindle, theta_max);
        REQUIRE(s.has_value());
        CHECK_THAT(s->score, WithinAbs(1.0f, 1e-5));
        CHECK_THAT(s->miss_angle_deg, WithinAbs(90.0f, 1e-3));
        CHECK_THAT(s->row_length_A, WithinAbs(0.0f, 1e-6));
    }

    SECTION("no cone, no answer") {
        const CrystalLattice latt(Coord(50, 0, 0), Coord(0, 0, 60), Coord(20, 70, 0));
        CHECK_FALSE(SpindleBlindFractionFromLattice(latt, spindle, 0.0f).has_value());
    }
}

TEST_CASE("SpindleTrigger_States", "[Indexing][Spindle]") {
    // The 0.5 threshold is a fold-angle gate: the score is monotone in the folded miss-angle and
    // crosses 0.5 at fold = 0.4040 * theta_max (verified root of the overlap closed form).
    CHECK_THAT(BlindConeSelfOverlap(0.4040f), WithinAbs(0.5f, 1e-3));

    CHECK(SpindleTriggerState(0.5f) == SpindleTrigger::Engage);
    CHECK(SpindleTriggerState(1.0f) == SpindleTrigger::Engage);
    CHECK(SpindleTriggerState(0.49f) == SpindleTrigger::DontEngage);
    CHECK(SpindleTriggerState(0.0f) == SpindleTrigger::DontEngage);
    // No value is a third state, and the one automation must treat as Engage: a measured 0 says
    // "one sweep loses nothing", absence says "nobody could look" - taking the second for the
    // first is the unrecoverable error.
    CHECK(SpindleTriggerState(std::nullopt) == SpindleTrigger::CannotSay);
}