The score gated on 60 spots, but the seed escalation stops at the leanest seed that indexes - 30 spots on precisely the clean frames a grid scan produces - so the value was absent exactly where beamline automation most needs it, and absence maps to "engage": the protocol would have fired on every good frame, which degenerates the trigger into "always". The floor itself stays where it was calibrated; what changes is what it gates. When no escalation pass could answer, one severity-only pass runs over the full spot list - the row search alone, no reduction, no refinement - purely to produce the number. The same was true of the indexer choice: only the FFT family computes a row shortlist, so a deployment configured with the known-cell indexer - the ordinary online stills path - never produced the score at all. Where the severity-only pass has no row search to run, the severity is read off the rows of the winning lattice instead, which any indexer produces: the lattice's shortest few distinct directions, as many as the FFT shortlist resolves in practice, fed through the same window and scoring with equal magnitudes. The count parity is load-bearing - a worst case over every enumerable lattice direction fires on 100% of harmless mounts of a generic triclinic cell against 74% for this selection at theta_max = 15 deg, and an always-firing trigger decides nothing - while the diad-detection rate stays 1.00 on the monoclinic classes either way, a dropped axis row being recovered by the pair normals exactly as an invisible one is. A frame that neither indexed nor reached the spot floor still reports nothing, which is the honest answer and maps to the recoverable error. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
159 lines
8.4 KiB
C++
159 lines
8.4 KiB
C++
// 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());
|
|
}
|
|
}
|