diff --git a/docs/CPU_DATA_ANALYSIS_INDEXING.md b/docs/CPU_DATA_ANALYSIS_INDEXING.md index 4fa756e0d..c4a86b77c 100644 --- a/docs/CPU_DATA_ANALYSIS_INDEXING.md +++ b/docs/CPU_DATA_ANALYSIS_INDEXING.md @@ -39,7 +39,7 @@ A spot is indexed if $\delta^2 < \tau^2$, where $\tau$ is the configured toleran For indexed spots, the reciprocal lattice point $\mathbf{p} = h\mathbf{a}^*+k\mathbf{b}^*+l\mathbf{c}^*$ is used to compute $\Delta_\mathrm{Ewald}(\mathbf{p})$ (stored as a diagnostic and later used in profile-radius estimation). -A frame is taken to be this crystal's when at least a fraction $g = 0.20$ of its in-resolution, non-ice spots index. On rotation data that decision is what admits the frame to integration, so its denominator matters: every spot handed to it that is not a reflection of this crystal argues against the frame. Where it cannot do that job — a lattice whose pooled spots fall below $g$ and which fewer than half the validation frames clear — every frame is integrated instead: the frames that clear $g$ are then only the upper tail of the same sparse population, not the frames the crystal was in, and admitting them alone dropped most of a small-molecule sweep and its completeness with it. +A frame is taken to be this crystal's when at least a fraction $g = 0.20$ of its in-resolution, non-ice spots index. On rotation data that decision is what admits the frame to integration, so its denominator matters: every spot handed to it that is not a reflection of this crystal argues against the frame. Where it cannot do that job — a lattice whose pooled spots fall below $g$ and which fewer than half the validation frames clear — every frame is integrated instead: the frames that clear $g$ are then only the upper tail of the same sparse population, not the frames the crystal was in, and admitting them alone dropped most of a small-molecule sweep and its completeness with it. A rotation frame is refined on the sweep's one cell, and a crystal whose cell grows with the dose leaves that cell behind, so a frame that fails $g$ is tried once more on the sweep's cell scaled isotropically (steps of 0.1 %, out to ±2 %) to the scale that puts most of its non-ice spots on the lattice, and is integrated on that cell where it then clears $g$. The scale has to gain more spots than the count's own noise, $\sqrt{n}$, so forty trial scales cannot turn a chance spot into a frame. A frame the sweep's cell already admits is integrated exactly as before. Measured on a $P2_1$ crystal losing 46 Ų over the sweep, the best scale walks from −0.1 % on the first frame to beyond +0.6 % at 115°, and without it the frames past that point — and the axial row that decides its screw axis — were refused. That test decides a **frame**. Whether a rotation run has a lattice **at all** is decided on the spots instead. A frame count comes from serial crystallography, where each image is its own experiment; a rotation sweep is one crystal and one orientation matrix, its frames are not independent of each other, and what such a count mostly measures is how many spots happen to land on a frame — a sweep carrying four spots an image cannot reach a six-spot bar on three frames in four however right the lattice is. The refusal therefore compares the fraction of *all* spots in the sampled frames that the lattice explains against what the same lattice explains when each frame's spots are put at **another frame's angle**: same lattice, same spots, same detector, same refinement, with only the claim that these spots were seen at *these* angles removed. That difference is the evidence, and it carries no spots-per-frame number anywhere, so nothing has to be chosen for a crystal that diffracts weakly. Measured over a hundred datasets the permuted level never exceeds 2.6 % and the smallest true margin is seventeen points. `--min-indexed-spots` (default 6, floor 4 — four is where a lattice stops being fitted by any three spots) still sets the reported indexing rate and the count the first pass *ranks* candidate lattices by; every rescue and every arbiter still counts frames. diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 1e717dde5..61d9ee36c 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -698,7 +698,7 @@ void IndexAndRefine::ProbeSupercellFrame(const DataMessage &msg, BraggPrediction std::optional IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettings &spot_finding_settings, - int64_t *spots_on_lattice) { + int64_t *spots_on_lattice, bool follow_cell_drift) { if (!indexer_ || !spot_finding_settings.indexing) return std::nullopt; @@ -712,7 +712,12 @@ IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettin if (!outcome.lattice_candidate) return std::nullopt; - if (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() != GeomRefinementAlgorithmEnum::None) + // Kept for the cell-drift retry below, on the only path that makes it. + std::optional unrefined; + if (rotation_indexer && follow_cell_drift) + unrefined = outcome; + const bool refine = experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() != GeomRefinementAlgorithmEnum::None; + if (refine) RefineGeometryIfNeeded(msg, outcome); if (!outcome.lattice_candidate.has_value()) @@ -721,8 +726,59 @@ IndexAndRefine::DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettin // AnalyzeIndexing answers "is this frame worth integrating"; msg.indexing_result carries the // stricter "does this frame index on its own", which on rotation is not the same question. if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate, outcome.extra_lattice_candidates, - spots_on_lattice, integrate_every_frame_)) - return std::nullopt; + spots_on_lattice, integrate_every_frame_)) { + // A rotation frame is refined on the sweep's cell, and a crystal whose cell grows with the dose + // leaves that cell behind: the frames it misses are refused as if the crystal had gone. So the + // other hypothesis - this frame's cell is the sweep's, scaled - is tried on a frame that is + // refused, and the frame is integrated on it where it then indexes. Measured on a P2_1 crystal + // losing B by 46 A^2 over the sweep: the scale that puts most spots on the lattice walks from + // -0.1 % on the first frame to +0.4 % at 50 deg and past +0.6 % at 115 deg, where the sweep's + // cell puts 60-80 of ~400 spots on the lattice and that one 110-150, so the frames past it were + // refused - and with them the axial row that decides the screw axis. A frame the sweep's cell + // already integrates is left exactly as it was. + if (!unrefined) + return std::nullopt; + const float tol = experiment.GetIndexingSettings().GetTolerance(); + // The spots the gate counts: an ice ring sits at one resolution, and a cell scaled to put its + // nodes on the ring would otherwise be scored on it. + const bool index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(); + std::vector counted; + for (const auto &sp : msg.spots) + if (index_ice_rings || !sp.ice_ring) + counted.push_back(sp); + const CrystalLattice &refined = *outcome.lattice_candidate; + const auto geom = outcome.experiment.GetDiffractionGeometry(); + auto scaled = [](const CrystalLattice &l, float s) { + return CrystalLattice(l.Vec0() * s, l.Vec1() * s, l.Vec2() * s); + }; + // Steps of 0.1 %, out to +-2 %, nearest first, so a tie goes to the smaller change. The scale + // has to put more spots on the lattice than the count's own noise, sqrt(count), or forty trial + // scales would turn one chance spot on a sparse frame into a frame that passes. + const int sweep_count = CountIndexedSpots(geom, refined, counted, tol * tol); + float best_scale = 1.0f; + int best_count = sweep_count; + for (int i = 1; i <= 20; i++) + for (const int sign : {1, -1}) { + const float s = 1.0f + 0.001f * static_cast(sign * i); + const int n = CountIndexedSpots(geom, scaled(refined, s), counted, tol * tol); + if (n > best_count) { + best_count = n; + best_scale = s; + } + } + if (best_count <= sweep_count + std::sqrt(static_cast(best_count))) + return std::nullopt; + outcome = *unrefined; + outcome.lattice_candidate = scaled(*unrefined->lattice_candidate, best_scale); + for (auto &el : outcome.extra_lattice_candidates) + el = scaled(el, best_scale); + if (refine) + RefineGeometryIfNeeded(msg, outcome); + if (!outcome.lattice_candidate.has_value() + || !AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate, + outcome.extra_lattice_candidates, spots_on_lattice, integrate_every_frame_)) + return std::nullopt; + } { std::unique_lock ul(reflections_mutex); @@ -737,7 +793,7 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, BraggPrediction &prediction, const BraggIntegrateFn &integrate, const BraggIntegrateFn &probe_integrate) { - auto outcome = DetermineRefineAnalyze(msg, spot_finding_settings); + auto outcome = DetermineRefineAnalyze(msg, spot_finding_settings, nullptr, /*follow_cell_drift=*/true); if (outcome && spot_finding_settings.quick_integration) QuickPredictAndIntegrate(msg, spot_finding_settings, prediction, integrate, *outcome, probe_integrate); } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index a5eb15a28..fa69e17fd 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -139,10 +139,12 @@ class IndexAndRefine { // Shared indexing path: determine the lattice/symmetry, refine geometry, and run AnalyzeIndexing. // Returns the outcome (ready for integration) when the frame indexes, nullopt otherwise. Both the // real per-image ProcessImage and the first-pass scheme validation go through this, so they cannot - // diverge. + // diverge. follow_cell_drift (ProcessImage only): a rotation frame the sweep's cell refuses is tried + // again on that cell scaled to the frame's own spots - the crystal's cell grows with the dose. std::optional DetermineRefineAnalyze(DataMessage &msg, const SpotFindingSettings &spot_finding_settings, - int64_t *spots_on_lattice = nullptr); + int64_t *spots_on_lattice = nullptr, + bool follow_cell_drift = false); void RefineGeometryIfNeeded(DataMessage &msg, IndexingOutcome &outcome); void QuickPredictAndIntegrate(DataMessage &msg, const SpotFindingSettings &spot_finding_settings,