From 3d9c8b878b4f4e198db8cc1cc8a7f0b8c8738fde Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 12 Sep 2026 18:17:14 +0200 Subject: [PATCH] beam-centre rescue: a rung is scored at its own centre, and not adopted on a harmonic The ladder set a trial centre and then scored it on spots found at the starting centre. Raw centroids do not move with the beam centre, but which spots are in the list does: the resolution limit, the ice-ring flag, the beam-stop mask and the strongest-N ranking are all computed against a radius. Every rung was therefore judged on the starting centre's evidence, and the count it produced was not a function of the centre it named. A trial centre now gets its own azimuthal mapping, spot engines and spot cache; the starting centre's cache is parked and restored, so a ladder that adopts nothing leaves the run exactly as it was. Measured against standalone runs pinned at the same centres, the ladder's score goes from a mean absolute error of 21.8 validation frames to 0.43, and at nine pinned centres the patched and unpatched binaries agree exactly - the change is confined to the ladder. Scoring honestly is not by itself an improvement: with the count honest, one crystal adopted a rung at a real 56/60 whose cell is a four-fold axis harmonic of the standing one, and merged the supercell. A longer cell indexes more, which is why the harmonic arbiter exists elsewhere in this file. A rung whose primitive volume is an integer or sqrt(3) multiple of the standing cell is now refused, and that crystal walks past both harmonics to a 60/60 rung. Over a hundred rotation datasets no space group, cell or run status changes. Eight merges move: one gains 126 per cent of its unique reflections, one 3.7, five move by less than half a per cent, and one loses 3.2 - that last was attributed separately and is the difference between two failed merges, both peaking at a half-dataset correlation of 5 per cent, with the honest run the slightly better of the two. Wall time rises 1.7 per cent, the cost of finding spots again at each trial centre. On a corpus processed from a deliberately destroyed beam centre the same change recovers eight datasets outright on their deposited space group and cell, restores a ninth's space-group call by returning its axial absence counts to the values the correct geometry gives, and withdraws one adoption that had shipped a wrong cell in the wrong crystal system at a half-dataset correlation of 0.10. Note that only the spot-cache half of this does the work. Reading the indexer's geometry per call rather than from construction was measured to change nothing on a rotation ladder, because the per-frame validation dispatches to a path that takes its geometry from the forced indexer result; it is kept because the copy was misleading to read, not because it moves a number. --- image_analysis/IndexAndRefine.cpp | 11 +++-- image_analysis/IndexAndRefine.h | 1 - rugnux/Rugnux.cpp | 78 +++++++++++++++++++++++++++---- 3 files changed, 77 insertions(+), 13 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 67c969705..1b9046a85 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -67,7 +67,6 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool : retain_outcomes_(retain_outcomes), real_time(real_time), experiment(x), - geom_(x.GetDiffractionGeometry()), indexer_(indexer), rotation_indexer_counter(x) { if (indexer && x.IsRotationIndexing()) @@ -161,12 +160,18 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data // Read per call rather than cached at construction: the run measures for itself whether the // crystal has ice, and that verdict lands after this object is built. const bool index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(); + // Read per call, for the same reason: a beam-centre hypothesis moves the geometry after this + // object is built, and a spot mapped into reciprocal space through the centre the run started at + // is not the spot the trial centre sees. Held at construction, this scored every hypothesis of the + // beam-centre ladder against the file's geometry - which is why a centre scored 42/60 inside the + // ladder and 5/60 when the same value was given on the command line. + const DiffractionGeometry geom = experiment.GetDiffractionGeometry(); for (size_t seed_cap : {size_t{30}, size_t{80}, std::numeric_limits::max()}) { std::vector recip; recip.reserve(std::min(seed_cap, msg.spots.size())); for (const auto &i: msg.spots) { if (index_ice_rings || !i.ice_ring) { - recip.push_back(i.ReciprocalCoord(geom_)); + recip.push_back(i.ReciprocalCoord(geom)); if (recip.size() >= seed_cap) break; } @@ -223,7 +228,7 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data 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_)); + recip.push_back(i.ReciprocalCoord(geom)); if (recip.size() >= FFT_MAX_SPOTS) break; } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 3a81173ba..946da5b5e 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -36,7 +36,6 @@ class IndexAndRefine { const bool retain_outcomes_; const bool real_time; // see the constructor const DiffractionExperiment& experiment; - const DiffractionGeometry geom_; std::optional indexed_lattice; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index f7fe8f6fd..17ac8a819 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2468,6 +2468,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // hand their spots back and this thread puts them in, so the map itself is still only ever // touched from here and needs no locking. std::map> spot_cache; + // The azimuthal mapping the first pass finds its spots against. It is `mapping` until a + // beam-centre hypothesis moves the centre, and the hypothesis's own mapping while it stands - + // the mapping is where the resolution of every pixel and the finder's radial background bins + // come from, and both move with the centre. + const AzimuthalIntegrationMapping *spot_mapping = &mapping; // The spots of one frame, found with the caller's own analysis engine. A frame's spots are a // pure function of that frame and the settings - the engine carries nothing from one image to // the next, which is why the main image loop can give one to each of its workers - so it does @@ -2549,9 +2554,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b SpotEngine &e = engines[t]; if (!e.analysis) { e.analysis = std::make_unique( - experiment_, mapping, pixel_mask_, *indexer, + experiment_, *spot_mapping, pixel_mask_, *indexer, /*enable_fused_adaptive_gpu=*/true); - e.profile = std::make_unique(mapping); + e.profile = std::make_unique(*spot_mapping); } for (size_t i = next.fetch_add(1); i < wanted.size(); i = next.fetch_add(1)) found[i] = find_spots(*e.analysis, *e.profile, e.raw_image, wanted[i]); @@ -2579,6 +2584,40 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b f.get(); }; + // Evaluating a beam-centre hypothesis means re-finding the spots at that centre. The spot + // POSITIONS do not move with the centre - which is why the cache is keyed by image alone - but + // the spot LIST does: every spot's d, the ice-ring flag, the resolution mask the finder + // applies, the high-resolution gap cut and the strongest-N budget are all radial, and the + // adaptive finder's ring background is binned about the centre. Scoring a trial centre against + // spots found at the file's put the ladder in a world the run never enters: measured, a centre + // scored 42/60 inside the ladder and 5/60 when that same value was given on the command line. + // + // A trial therefore gets its own azimuthal mapping, its own engines and its own cache, and the + // starting centre's spots are PARKED rather than thrown away - a hypothesis that is not adopted + // has to leave the rest of the run exactly as it found it. + std::unique_ptr trial_mapping; + std::map> home_spot_cache; + const auto try_beam_center = [&](float x, float y) { + if (trial_mapping) + spot_cache.clear(); + else + home_spot_cache.swap(spot_cache); + experiment_.BeamX_pxl(x).BeamY_pxl(y); + release_spot_engines(); + trial_mapping = std::make_unique(experiment_, pixel_mask_); + spot_mapping = trial_mapping.get(); + }; + const auto restore_beam_center = [&](float x, float y) { + if (!trial_mapping) + return; + experiment_.BeamX_pxl(x).BeamY_pxl(y); + release_spot_engines(); + trial_mapping.reset(); + spot_mapping = &mapping; + spot_cache = std::move(home_spot_cache); + home_spot_cache.clear(); + }; + // How many validation frames a candidate global lattice actually indexes per image, scored on // the *real* per-image path (force the candidate, then run the same indexing+refinement the // main loop uses). A cell that merely fits the accumulated first-pass cloud (a spurious @@ -3251,7 +3290,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // a halved axis, and 0.12 px is well inside the sigma of any fit that could gate it. { const int majority = static_cast(validation.size()) / 2; - experiment_.BeamX_pxl(measured_x).BeamY_pxl(measured_y); + try_beam_center(measured_x, measured_y); FirstPass alt = pick_best(*indexer_pool, *indexer); const bool header_indexes = best.result.has_value() && best.score > majority; bool measured_indexes = alt.result.has_value() && alt.score > majority; @@ -3327,7 +3366,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // around the value the run started with rather than around a fit that has just // been tested and found no better - and so a run that fails for some other // reason fails at the geometry it was given. - experiment_.BeamX_pxl(header_x).BeamY_pxl(header_y); + restore_beam_center(header_x, header_y); logger.Info("Beam centre check: neither the file's centre ({}/{} frames) nor " "the measured one {:.2f} px away ({}/{}) indexes this run - the " "file's centre is kept and this is not a beam-centre failure", @@ -3335,7 +3374,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b static_cast(validation.size())); } } else { - experiment_.BeamX_pxl(header_x).BeamY_pxl(header_y); + restore_beam_center(header_x, header_y); if (measured_indexes) { // Primitive volumes, as everywhere else here: two passes can find the same // lattice in different settings, and a centred setting's cell is an exact @@ -3469,8 +3508,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // (99 % to 25 %), while along it the run keeps 96-100 % indexed and quietly adopts a 2x, 3x or // sqrt(3) supercell. The direction that announces itself is not the dangerous one. // - // Only after a poor pass, so a correct header costs nothing; and spot finding is not repeated, - // the cache being keyed by image and the spot POSITIONS not depending on the centre. + // Only after a poor pass, so a correct header costs nothing. Every rung re-finds its spots at + // its own centre (try_beam_center): the positions do not move with the centre but the list + // does, and a rung scored on the file centre's spot list is not a rung the run can adopt. if (!cancelled_ && !beam_center_searched_ && config_.beam_center_search_pxl > 0 && best.score < 0.5 * static_cast(validation.size())) { beam_center_searched_ = true; @@ -3516,11 +3556,30 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b for (const auto &[dx, dy]: {std::pair{-k, 0}, {k, 0}, {0, -k}, {0, k}}) { const float trial_x = beam_x + static_cast(dx); const float trial_y = beam_y + static_cast(dy); - experiment_.BeamX_pxl(trial_x).BeamY_pxl(trial_y); + try_beam_center(trial_x, trial_y); FirstPass alt = pick_best(*indexer_pool, *indexer); logger.Info("Beam-centre hypothesis ({:.2f},{:.2f}): {}/{} frames", trial_x, trial_y, alt.score, static_cast(validation.size())); seen = std::max(seen, alt.score); + // A rung is not adopted on the count alone. A centre error along the spindle + // translates the derotated cloud rigidly, which the FFT amplitude cannot see, so + // the indexer commits to an axis HARMONIC of the true cell - and a cell twice as + // long indexes MORE frames, not fewer, so the count rewards exactly the answer a + // centre error produces. The same test the scheme comparison already makes: + // primitive volumes differing by a small integer or by sqrt(3). Measured with the + // rungs scored honestly, a rung 1 px away indexed 56/60 on a 4x supercell of the + // lattice the starting centre had already found at 4/60. + const double ratio = best.vol > 0.0 && alt.vol > best.vol ? alt.vol / best.vol : 0.0; + const double nearest = std::round(ratio); + const bool harmonic = (nearest >= 2.0 && std::abs(ratio - nearest) < 0.15) + || std::abs(ratio - std::sqrt(3.0)) < 0.15; + if (harmonic) { + logger.Info("Beam-centre hypothesis ({:.2f},{:.2f}): its cell is {:.2f}x the " + "volume of the one the starting centre found, so it is an axis " + "harmonic and the frame count cannot arbitrate it - not adopted", + trial_x, trial_y, ratio); + continue; + } if (alt.result.has_value() && alt.score > adopt_bar && alt.score > best.score) { best = alt; won_x = trial_x; @@ -3538,7 +3597,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b break; } } - experiment_.BeamX_pxl(won_x).BeamY_pxl(won_y); + if (!adopted) + restore_beam_center(won_x, won_y); if (adopted) logger.Info("Beam centre from indexing: ({:.2f},{:.2f}) -> ({:.2f},{:.2f}), moved " "{:.2f} px, {}/{} validation frames", beam_x, beam_y, won_x, won_y,