mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity

Merge branch 'regr-5ky6' into rc175

Browse Source
This commit is contained in:
leonarski_f
2026-10-08 22:01:38 +02:00
parent 41cdde68ea a825fb6cba
commit eac0833c88
3 changed files with 66 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/CPU_DATA_ANALYSIS_INDEXING.md
+1 -1
View File
@@ -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.
image_analysis/IndexAndRefine.cpp
+61 -5
View File
@@ -698,7 +698,7 @@ void IndexAndRefine::ProbeSupercellFrame(const DataMessage &msg, BraggPrediction
std::optional<IndexAndRefine::IndexingOutcome>
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<IndexingOutcome> 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<SpotToSave> 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<float>(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<float>(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);
}
image_analysis/IndexAndRefine.h
+4 -2
View File
@@ -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<IndexingOutcome> 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,