From 89962574ef5dddff7452ab8aa18bba2773b9952e Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 4 Sep 2026 09:58:43 +0200 Subject: [PATCH] spindle: the severity no longer rides on the indexing seed or on which indexer is configured 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 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- image_analysis/IndexAndRefine.cpp | 42 ++++++++++++++++-- image_analysis/indexing/FFTIndexer.cpp | 18 +++++--- image_analysis/indexing/FFTIndexer.h | 1 + image_analysis/indexing/Indexer.cpp | 14 ++++-- image_analysis/indexing/Indexer.h | 10 ++++- image_analysis/indexing/IndexerThreadPool.cpp | 11 ++--- image_analysis/indexing/IndexerThreadPool.h | 8 +++- .../indexing/SpindleBlindFraction.cpp | 43 +++++++++++++++++++ .../indexing/SpindleBlindFraction.h | 10 +++++ tests/IndexingUnitTest.cpp | 43 +++++++++++++++++++ tests/SpindleBlindFractionTest.cpp | 32 ++++++++++++++ 11 files changed, 211 insertions(+), 21 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 1c95af4d8..67c969705 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -175,10 +175,9 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data any_executed |= res.executed; // Kept from whichever run could measure it - the largest seed that answered - rather than from // the run whose lattice won: the severity describes the frame's lattice rows, not the cell that - // closed on them, so a frame that indexed nothing still has one. It stays absent while the - // escalation never leaves the lean seed, which is below the spot floor the score needs - // (SPINDLE_MIN_SPOTS in FFTIndexer.cpp) - so on frames that index cleanly at 30 spots there is - // no value at all. That is a gap, not a "no problem": see the note at SPINDLE_MIN_SPOTS. + // closed on them, so a frame that indexed nothing still has one. When the escalation never + // leaves a seed below the score's spot floor it stays absent here, and the severity-only pass + // after the loop supplies it instead. if (res.spindle_blind_fraction) spindle_blind_fraction = res.spindle_blind_fraction; if (!res.lattice.empty()) { @@ -211,6 +210,41 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data if (any_executed) msg.indexing_result = false; + // The severity's spot floor (SPINDLE_MIN_SPOTS) was calibrated on a frame's whole spot list, + // but the escalation stops at the leanest seed that indexes - 30 spots on precisely the clean + // frames a grid scan produces - so riding on the indexing seed left the score absent exactly + // where automation most needs it, and absence maps to "engage" (see SpindleBlindFraction.h). + // Decouple the two: when no escalation pass could answer, spend one severity-only row pass over + // the full spot list (no reduction, no refinement); where that path does not exist - an indexer + // without a row search - or refuses, read the severity off the rows of the winning lattice, + // which any indexer produces. + if (!spindle_blind_fraction && experiment.GetGoniometer().has_value()) { + std::vector recip; + recip.reserve(std::min(msg.spots.size(), FFT_MAX_SPOTS)); + for (const auto &i : msg.spots) { + if (index_ice_rings || !i.ice_ring) { + recip.push_back(i.ReciprocalCoord(geom_)); + if (recip.size() >= FFT_MAX_SPOTS) + break; + } + } + if (recip.size() >= SPINDLE_MIN_SPOTS) { + const auto algorithm = experiment.GetIndexingAlgorithm(); + if (algorithm == IndexingAlgorithmEnum::FFT || algorithm == IndexingAlgorithmEnum::FFTW) { + const auto res = indexer_->Run(experiment, recip, /*severity_only=*/true); + spindle_blind_fraction = res.spindle_blind_fraction; + } + if (!spindle_blind_fraction && !indexer_result.lattice.empty()) { + const float theta_max_deg = SpindleThetaMax_deg(experiment.GetWavelength_A(), + experiment.GetDetectorMaxResolution_A()); + if (const auto severity = SpindleBlindFractionFromLattice( + indexer_result.lattice.front(), experiment.GetGoniometer()->GetAxis(), + theta_max_deg)) + spindle_blind_fraction = severity->score; + } + } + } + msg.spindle_blind_fraction = spindle_blind_fraction; if (!indexer_result.lattice.empty()) { diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp index 773183545..7bec2f92d 100644 --- a/image_analysis/indexing/FFTIndexer.cpp +++ b/image_analysis/indexing/FFTIndexer.cpp @@ -466,6 +466,16 @@ std::optional FFTIndexer::SearchCap(const std::vector &coord, size return found; } +std::optional FFTIndexer::RunSeverityOnly(const std::vector &coord) { + if (!spindle_axis || spindle_theta_max_deg <= 0 || coord.size() < SPINDLE_MIN_SPOTS) + return {}; + const size_t nspots = std::min(coord.size(), static_cast(FFT_MAX_SPOTS)); + ExecuteFFT(coord, nspots); + std::vector magnitudes; + const auto shortlist = FilterFFTResults(30, &magnitudes); + return SpindleBlindFraction(shortlist, magnitudes, *spindle_axis, spindle_theta_max_deg); +} + std::vector FFTIndexer::RunInternal(const std::vector &coord, size_t nspots) { if (nspots > coord.size()) nspots = coord.size(); @@ -488,11 +498,9 @@ std::vector FFTIndexer::RunInternal(const std::vector &co // not the still's own spot resolution, which on a weak attenuated frame understates the sweep's // real loss. // The floor was calibrated on a pass over a frame's WHOLE spot list, but nspots here is whatever - // the caller fed this call. IndexAndRefine escalates 30 -> 80 -> all and usually stops at 30, so - // on a frame that indexes cleanly the score is simply never computed. Closing that gap means - // measuring the severity on the full list rather than riding on the indexing seed, which costs a - // pass this deliberately does not spend; until then the value is absent far more often than the - // floor alone implies. + // the caller fed this call - IndexAndRefine escalates 30 -> 80 -> all and usually stops at 30, + // so on a frame that indexes cleanly nothing is computed here; the caller then asks for the + // severity with a RunSeverityOnly pass over the full list instead of leaving it absent. spindle_severity = {}; if (spindle_axis && spindle_theta_max_deg > 0 && nspots >= SPINDLE_MIN_SPOTS) spindle_severity = SpindleBlindFraction(shortlist, magnitudes, *spindle_axis, diff --git a/image_analysis/indexing/FFTIndexer.h b/image_analysis/indexing/FFTIndexer.h index a483eb9ba..d29cce9bc 100644 --- a/image_analysis/indexing/FFTIndexer.h +++ b/image_analysis/indexing/FFTIndexer.h @@ -57,6 +57,7 @@ protected: // Filled by RunInternal from the shortlist it already computes; read out by Indexer::Run. std::optional spindle_severity; std::optional GetSpindleSeverity() const override { return spindle_severity; } + std::optional RunSeverityOnly(const std::vector &coord) override; virtual void ExecuteFFT(const std::vector &coord, size_t nspots) = 0; // Called after direction_vectors is rewritten, for implementations that keep a copy of it. diff --git a/image_analysis/indexing/Indexer.cpp b/image_analysis/indexing/Indexer.cpp index 82dc6dd90..e7725a8ce 100644 --- a/image_analysis/indexing/Indexer.cpp +++ b/image_analysis/indexing/Indexer.cpp @@ -17,16 +17,22 @@ void Indexer::Setup(const DiffractionExperiment& experiment) { SetupUnitCell(experiment.GetUnitCell()); } -IndexerResult Indexer::Run(const std::vector &coord) { +IndexerResult Indexer::Run(const std::vector &coord, bool severity_only) { IndexerResult ret; auto start = std::chrono::steady_clock::now(); - ret.lattice = RunInternal(coord, coord.size()); + std::optional 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 duration = end - start; ret.indexing_time_s = duration.count(); - ret.executed = true; - if (const auto severity = GetSpindleSeverity()) { + 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; diff --git a/image_analysis/indexing/Indexer.h b/image_analysis/indexing/Indexer.h index 4f5377c08..cd2190145 100644 --- a/image_analysis/indexing/Indexer.h +++ b/image_analysis/indexing/Indexer.h @@ -55,10 +55,18 @@ protected: virtual std::vector RunInternal(const std::vector &coord, size_t nspots) = 0; // Set by RunInternal when the implementation computes it; read out by Run. virtual std::optional GetSpindleSeverity() const { return {}; } + // The row pass alone, over the whole of `coord`, to produce the spindle severity without + // reducing or refining anything. Implemented by the FFT-family indexers, which own a row + // search; the others have nothing to answer with and return nothing. + virtual std::optional RunSeverityOnly(const std::vector &coord) { return {}; } public: virtual ~Indexer() = default; void Setup(const DiffractionExperiment& experiment); - IndexerResult Run(const std::vector &coord); + // severity_only = true runs RunSeverityOnly instead of indexing: no lattice comes back and + // `executed` stays false, since nothing was an indexing attempt. Used when the seed escalation + // indexed a frame from fewer spots than the severity's floor - the score must not be absent on + // precisely the frames that index cleanly. + IndexerResult Run(const std::vector &coord, bool severity_only = false); }; diff --git a/image_analysis/indexing/IndexerThreadPool.cpp b/image_analysis/indexing/IndexerThreadPool.cpp index 8684c3307..56b925c18 100644 --- a/image_analysis/indexing/IndexerThreadPool.cpp +++ b/image_analysis/indexing/IndexerThreadPool.cpp @@ -148,7 +148,7 @@ void IndexerThread::Worker(int threadid) { Indexer &indexer = **slot; indexer.Setup(input->experiment); - tmp_result = std::make_unique(indexer.Run(input->recip)); + tmp_result = std::make_unique(indexer.Run(input->recip, input->severity_only)); } catch (std::exception &e) { // Hand the failure back as a result carrying the reason. A nullptr here was // indistinguishable from a worker that was never dispatched, and both then read @@ -178,7 +178,7 @@ void IndexerThread::Finalize() { } std::unique_ptr IndexerThread::Run(const DiffractionExperiment &experiment, - const std::vector &recip) { + const std::vector &recip, bool severity_only) { std::unique_ptr tmp_result; { std::unique_lock lock(m); @@ -186,7 +186,7 @@ std::unique_ptr IndexerThread::Run(const DiffractionExperiment &e return nullptr; if (state != TaskState::IDLE) return nullptr; - task_input = std::make_unique(std::cref(experiment), std::cref(recip)); + task_input = std::make_unique(std::cref(experiment), std::cref(recip), severity_only); state = TaskState::READY; } c_start.notify_one(); @@ -231,7 +231,8 @@ int IndexerThreadPool::GetFreeWorker() { return -1; } -IndexerResult IndexerThreadPool::Run(const DiffractionExperiment &experiment, const std::vector &recip) { +IndexerResult IndexerThreadPool::Run(const DiffractionExperiment &experiment, const std::vector &recip, + bool severity_only) { const auto algorithm = experiment.GetIndexingAlgorithm(); if (algorithm == IndexingAlgorithmEnum::None) return IndexerResult{.lattice = {}, .indexing_time_s = 0, .executed = false}; @@ -285,7 +286,7 @@ IndexerResult IndexerThreadPool::Run(const DiffractionExperiment &experiment, co std::unique_ptr result; if (task >= 0) { try { - result = tasks[task]->Run(experiment, recip); + result = tasks[task]->Run(experiment, recip, severity_only); } catch (const std::exception &e) { spdlog::error("Indexer thread failed: {}", e.what()); result = std::make_unique(IndexerResult{ diff --git a/image_analysis/indexing/IndexerThreadPool.h b/image_analysis/indexing/IndexerThreadPool.h index 2fb319e9e..116444b44 100644 --- a/image_analysis/indexing/IndexerThreadPool.h +++ b/image_analysis/indexing/IndexerThreadPool.h @@ -44,6 +44,7 @@ class IndexerThread { struct TaskInput { const DiffractionExperiment &experiment; const std::vector &recip; + const bool severity_only; }; // Held by value: with IndexerConstruction::OnFirstUse the worker builds its indexer long after @@ -66,7 +67,8 @@ class IndexerThread { public: IndexerThread(const IndexingSettings& settings, int threadid, IndexerConstruction construction); ~IndexerThread(); - std::unique_ptr Run(const DiffractionExperiment &experiment, const std::vector &recip); + std::unique_ptr Run(const DiffractionExperiment &experiment, const std::vector &recip, + bool severity_only = false); void Finalize(); }; @@ -82,7 +84,9 @@ class IndexerThreadPool { public: IndexerThreadPool(const IndexingSettings& settings, IndexerConstruction construction = IndexerConstruction::Preconstruct); - IndexerResult Run(const DiffractionExperiment& experiment, const std::vector& recip); + // severity_only skips indexing and produces just the spindle severity - see Indexer::Run. + IndexerResult Run(const DiffractionExperiment& experiment, const std::vector& recip, + bool severity_only = false); }; diff --git a/image_analysis/indexing/SpindleBlindFraction.cpp b/image_analysis/indexing/SpindleBlindFraction.cpp index 9d61cbd78..d7e9bd623 100644 --- a/image_analysis/indexing/SpindleBlindFraction.cpp +++ b/image_analysis/indexing/SpindleBlindFraction.cpp @@ -134,3 +134,46 @@ std::optional SpindleBlindFraction(const std::vector &ro return ret; } + +std::optional SpindleBlindFractionFromLattice(const CrystalLattice &lattice, + const Coord &spindle, + float theta_max_deg) { + // Candidate rows: the direct lattice's shortest few distinct directions, drawn from the index + // box up to +/-2 - the range in which the symmetry axes of a reduced or conventional basis + // lie. The count matches what the FFT shortlist resolves in practice (four or five distinct + // rows - see FilterFFTResults), so the bound is taken over comparable evidence on either path. + // That parity is load-bearing: a worst case over every enumerable direction saturates towards + // "always engage" - measured on a generic triclinic cell it fires on 100% of harmless mounts, + // against 74% for this selection at theta_max = 15 deg - and an always-firing trigger decides + // nothing. The diad-detection rate stays 1.00 on the monoclinic classes either way, because a + // dropped axis row is recovered by the pair normals exactly as an invisible one is. + std::vector all; + all.reserve(62); + for (int u = 0; u <= 2; u++) + for (int v = (u == 0) ? 0 : -2; v <= 2; v++) + for (int w = (u == 0 && v == 0) ? 1 : -2; w <= 2; w++) + all.push_back(lattice.Vec0() * static_cast(u) + + lattice.Vec1() * static_cast(v) + + lattice.Vec2() * static_cast(w)); + std::sort(all.begin(), all.end(), + [](const Coord &a, const Coord &b) { return a.Length() < b.Length(); }); + + constexpr size_t MAX_LATTICE_ROWS = 6; + const float cos_5_deg = std::cos(5.0f * static_cast(M_PI) / 180.0f); + std::vector rows; + for (const auto &r : all) { + if (rows.size() >= MAX_LATTICE_ROWS + || r.Length() > MAX_ROW_LENGTH_RATIO * all.front().Length()) + break; + bool distinct = true; + for (const auto &k : rows) + if (std::fabs(r * k) / (r.Length() * k.Length()) > cos_5_deg) { + distinct = false; + break; + } + if (distinct) + rows.push_back(r); + } + const std::vector magnitudes(rows.size(), 1.0f); + return SpindleBlindFraction(rows, magnitudes, spindle, theta_max_deg); +} diff --git a/image_analysis/indexing/SpindleBlindFraction.h b/image_analysis/indexing/SpindleBlindFraction.h index 2f80366c5..9ac03a55d 100644 --- a/image_analysis/indexing/SpindleBlindFraction.h +++ b/image_analysis/indexing/SpindleBlindFraction.h @@ -8,6 +8,7 @@ #include #include "../../common/Coord.h" +#include "../../common/CrystalLattice.h" // How much of a single sweep's blind cone this orientation makes unrecoverable. // @@ -75,3 +76,12 @@ std::optional SpindleBlindFraction(const std::vector &ro const std::vector &magnitudes, const Coord &spindle, float theta_max_deg); + +// The same severity read off an indexed lattice instead of an FFT shortlist, for the paths that +// index without a row search (a known-cell indexer on the online stills path). The candidate rows +// are the lattice's shortest few distinct directions - as many as the FFT shortlist resolves in +// practice, so the bound is over comparable evidence on either path - with equal magnitudes: a +// lattice does not rank its rows, and all of them are equally real. +std::optional SpindleBlindFractionFromLattice(const CrystalLattice &lattice, + const Coord &spindle, + float theta_max_deg); diff --git a/tests/IndexingUnitTest.cpp b/tests/IndexingUnitTest.cpp index 545c68c19..49c0365d4 100644 --- a/tests/IndexingUnitTest.cpp +++ b/tests/IndexingUnitTest.cpp @@ -246,6 +246,49 @@ TEST_CASE("FFTIndexer","[Indexing]") { logger.Info("Time: {} ms", std::chrono::duration_cast(end - start).count()); } +TEST_CASE("FFTIndexer_SpindleSeverity", "[Indexing][Spindle]") { + // End to end through the real indexer: a crystal whose shortest row lies on the spindle must + // come back with a severity of 1 from a normal run, from a severity-only run - which must not + // index anything - and not at all when the frame is below the spot floor. + UnitCell uc{39, 45, 78, 90, 90, 90}; + CrystalLattice cl(uc); + + DiffractionExperiment experiment(DetJF4M()); + experiment.DetectorDistance_mm(75).BeamY_pxl(1136).BeamX_pxl(1090).IncidentEnergy_keV(12.4); + // The 39 A axis of the cell lies along x, and so does the spindle. + experiment.Goniometer(GoniometerAxis("omega", 0, 0.1f, Coord(1, 0, 0), {})); + + IndexingSettings settings; + settings.Algorithm(IndexingAlgorithmEnum::FFT) + .FFT_MaxUnitCell_A(250.0).FFT_HighResolution_A(2 * M_PI / 3.0); + experiment.ImportIndexingSettings(settings).SetUnitCell(uc); + + std::unique_ptr indexer = CreateIndexer(experiment); + REQUIRE(indexer); + indexer->Setup(experiment); + + std::vector vec; + for (int h = -2; h < 10; h++) + for (int k = -5; k < 10; k++) + for (int l = -3; l < 10; l++) + vec.push_back(h * cl.Astar() + k * cl.Bstar() + l * cl.Cstar()); + + auto full = indexer->Run(vec); + REQUIRE(full.spindle_blind_fraction.has_value()); + CHECK_THAT(*full.spindle_blind_fraction, Catch::Matchers::WithinAbs(1.0f, 1e-4)); + + auto severity_only = indexer->Run(vec, /*severity_only=*/true); + CHECK(severity_only.lattice.empty()); + CHECK_FALSE(severity_only.executed); + REQUIRE(severity_only.spindle_blind_fraction.has_value()); + CHECK_THAT(*severity_only.spindle_blind_fraction, Catch::Matchers::WithinAbs(1.0f, 1e-4)); + + // Below the spot floor there is no value - the CANNOT-SAY state - not a middling one. + const std::vector few(vec.begin(), vec.begin() + 40); + auto starved = indexer->Run(few, /*severity_only=*/true); + CHECK_FALSE(starved.spindle_blind_fraction.has_value()); +} + TEST_CASE("PostIndexingRefinement_MultiLattice_TwoCrystals_BraggPrediction","[Indexing]") { Logger logger("PostIndexingRefinement_MultiLattice_TwoCrystals_BraggPrediction"); diff --git a/tests/SpindleBlindFractionTest.cpp b/tests/SpindleBlindFractionTest.cpp index 67ee26063..4823ba27b 100644 --- a/tests/SpindleBlindFractionTest.cpp +++ b/tests/SpindleBlindFractionTest.cpp @@ -124,3 +124,35 @@ TEST_CASE("SpindleBlindFraction_Rows", "[Indexing][Spindle]") { 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()); + } +}