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()); + } +}