From 49963ab855e1e9bc02247fd6db9ec96716353019 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 14 Jul 2026 15:35:40 +0200 Subject: [PATCH] rugnux: move indexing-ambiguity resolution out of ScaleOnTheFly Indexing-ambiguity detection/resolution is no longer part of scaling: it is a dedicated per-image step for stills (or a post-full-merge step for rotation). Pull the stills per-image resolver out of ScaleOnTheFly into a new ReindexAmbiguityResolver in ReindexAmbiguity.{h,cpp}, alongside the existing rotation free functions (ChooseReindex / ReferenceIntensityCC). Both workflows now share the operator generation and the best-op selection (new file-local PickBestReindex helper); ChooseReindex is rewired onto it. ScaleOnTheFly is now purely a scaling engine (no ambiguity_ops, no resolve_ambiguity flag, no ResolveIndexingAmbiguity). Pure refactor: behaviour is unchanged. IndexAndRefine now holds a ReindexAmbiguityResolver and calls Resolve() at the same call site. Co-Authored-By: Claude Opus 4.8 (1M context) --- image_analysis/IndexAndRefine.cpp | 11 +- image_analysis/IndexAndRefine.h | 3 +- .../scale_merge/ReindexAmbiguity.cpp | 124 ++++++++++++++++-- image_analysis/scale_merge/ReindexAmbiguity.h | 32 +++++ image_analysis/scale_merge/ScaleOnTheFly.cpp | 59 +-------- image_analysis/scale_merge/ScaleOnTheFly.h | 14 +- rugnux/Rugnux.cpp | 2 +- 7 files changed, 153 insertions(+), 92 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index f2178a4b..036e00ae 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -10,6 +10,7 @@ #include "indexing/FFTIndexer.h" #include "indexing/MultiLatticeSearch.h" #include "lattice_search/LatticeSearch.h" +#include "scale_merge/ReindexAmbiguity.h" #include "scale_merge/ScaleOnTheFly.h" IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer, bool retain_outcomes) @@ -485,14 +486,14 @@ std::optional IndexAndRefine::FinalizeRotationIndexing() } IndexAndRefine &IndexAndRefine::ReferenceIntensities(std::vector &reference) { - // An external reference is trusted to be in the correct hand, so let the per-image scaler also break - // the merohedral indexing ambiguity against it (serial stills index each crystal independently). - scaling_engine = std::make_unique(experiment, reference, /*resolve_ambiguity=*/true); + // An external reference is trusted to be in the correct hand, so use it to break the merohedral + // indexing ambiguity per image (serial stills index each crystal independently). + reindex_resolver = std::make_unique(experiment, reference); return *this; } void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) { - if (!scaling_engine) + if (!reindex_resolver) return; // The external reference fixes the cell/space group, breaks the indexing ambiguity and reports CCref, @@ -506,7 +507,7 @@ void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) { return; auto scaling_start_time = std::chrono::steady_clock::now(); - scaling_engine->ResolveIndexingAmbiguity(outcome.reflections); + reindex_resolver->Resolve(outcome.reflections); auto scaling_end_time = std::chrono::steady_clock::now(); msg.image_scale_time_s = std::chrono::duration(scaling_end_time - scaling_start_time).count(); } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 3f912a8b..ed6b85d3 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -17,6 +17,7 @@ #include "lattice_search/LatticeSearch.h" #include "rotation_indexer/RotationIndexer.h" #include "rotation_indexer/RotationIndexerCounter.h" +#include "scale_merge/ReindexAmbiguity.h" #include "scale_merge/ScaleOnTheFly.h" #include "scale_merge/ScalingResult.h" #include "IntegrationOutcome.h" @@ -81,7 +82,7 @@ class IndexAndRefine { const BraggIntegrateFn &integrate, const IndexingOutcome &outcome); - std::unique_ptr scaling_engine; + std::unique_ptr reindex_resolver; void ScaleImage(DataMessage &msg, IntegrationOutcome& outcome); std::optional RotationAngle(int64_t image) const; // mid-exposure angle for the indexer diff --git a/image_analysis/scale_merge/ReindexAmbiguity.cpp b/image_analysis/scale_merge/ReindexAmbiguity.cpp index 929f641e..59e53558 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.cpp +++ b/image_analysis/scale_merge/ReindexAmbiguity.cpp @@ -10,6 +10,42 @@ #include "HKLKey.h" +namespace { + constexpr size_t MIN_REFLECTIONS = 20; + + double SafeInv(double x, double fallback) { + if (!std::isfinite(x) || x == 0.0) + return fallback; + return 1.0 / x; + } + + struct BestReindex { + gemmi::Op op = gemmi::Op::identity(); + bool is_identity = true; + double score = 0.0; + double identity_score = 0.0; + }; + + // Among identity (the baseline, scored `identity_score`) and `ops`, return the operator with the + // highest `score_op(op)`. A non-finite score never wins (too little overlap to decide). Shared by the + // stills per-image resolver and the rotation post-merge ChooseReindex. + BestReindex PickBestReindex(double identity_score, const std::vector &ops, + const std::function &score_op) { + BestReindex best; + best.identity_score = identity_score; + best.score = identity_score; + for (const auto &op : ops) { + const double sc = score_op(op); + if (std::isfinite(sc) && (!std::isfinite(best.score) || sc > best.score)) { + best.score = sc; + best.op = op; + best.is_identity = false; + } + } + return best; + } +} + std::vector ReindexAmbiguityOperators(const UnitCell &cell, int space_group_number, double max_obliquity_deg) { const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(space_group_number); @@ -37,21 +73,18 @@ ReindexChoice ChooseReindex(const std::vector &merged, const UnitCell &cell, int space_group_number, const std::function &)> &score, double max_obliquity_deg) { - ReindexChoice best; - best.identity_score = score(merged); - best.score = best.identity_score; - const auto ops = ReindexAmbiguityOperators(cell, space_group_number, max_obliquity_deg); - best.n_candidates = 1 + static_cast(ops.size()); - for (const auto &op : ops) { - const double s = score(ReindexReflections(merged, op)); - if (s > best.score) { - best.score = s; - best.op = op; - best.is_identity = false; - } - } - return best; + const BestReindex best = PickBestReindex( + score(merged), ops, + [&](const gemmi::Op &op) { return score(ReindexReflections(merged, op)); }); + + ReindexChoice choice; + choice.op = best.op; + choice.is_identity = best.is_identity; + choice.score = best.score; + choice.identity_score = best.identity_score; + choice.n_candidates = 1 + static_cast(ops.size()); + return choice; } double ReferenceIntensityCC(const std::vector &merged, @@ -87,3 +120,66 @@ double ReferenceIntensityCC(const std::vector &merged, const double vy = n * syy - sy * sy; return (vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : 0.0; } + +ReindexAmbiguityResolver::ReindexAmbiguityResolver(const DiffractionExperiment &x, + const std::vector &reference) + : s(x.GetScalingSettings()), + hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)) { + for (const auto &r : reference) + reference_data[hkl_key_generator(r)] = r.I; + if (x.GetUnitCell().has_value() && x.GetSpaceGroupNumber().has_value()) + ops = ReindexAmbiguityOperators(*x.GetUnitCell(), static_cast(*x.GetSpaceGroupNumber())); +} + +bool ReindexAmbiguityResolver::Accept(const Reflection &r) const { + if (r.on_ice_ring) // ice-contaminated intensity would bias the correlation; keep it out + return false; + return AcceptReflection(r, s.GetHighResolutionLimit_A()); +} + +double ReindexAmbiguityResolver::ReferenceCC(const std::vector &reflections, + const gemmi::Op &op) const { + double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; + size_t n = 0; + for (const auto &r : reflections) { + if (!Accept(r) || r.partiality < s.GetMinPartiality()) + continue; + if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || r.partiality <= 0.0f) + continue; + const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); + const auto it = reference_data.find(hkl_key_generator(h[0], h[1], h[2])); + if (it == reference_data.end()) + continue; + const double x = static_cast(r.I) * SafeInv(r.rlp, 1.0) / r.partiality; + const double y = it->second; + if (!std::isfinite(x) || !std::isfinite(y)) + continue; + sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; + ++n; + } + if (n < MIN_REFLECTIONS) + return NAN; + const double nd = static_cast(n); + const double cov = sxy - sx * sy / nd; + const double vx = sxx - sx * sx / nd; + const double vy = syy - sy * sy / nd; + return (vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : NAN; +} + +// Serial stills index each crystal in one of the merohedrally-equivalent hands at random; pick, for this +// image alone, the reindexing whose intensities correlate best with the external reference and apply it. +void ReindexAmbiguityResolver::Resolve(std::vector &reflections) const { + if (ops.empty()) + return; + + const BestReindex best = PickBestReindex( + ReferenceCC(reflections, gemmi::Op::identity()), ops, + [&](const gemmi::Op &op) { return ReferenceCC(reflections, op); }); + if (best.is_identity) + return; + + for (auto &r : reflections) { + const gemmi::Op::Miller h = best.op.apply_to_hkl({{r.h, r.k, r.l}}); + r.h = h[0]; r.k = h[1]; r.l = h[2]; + } +} diff --git a/image_analysis/scale_merge/ReindexAmbiguity.h b/image_analysis/scale_merge/ReindexAmbiguity.h index 300276d3..93f44ce5 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.h +++ b/image_analysis/scale_merge/ReindexAmbiguity.h @@ -4,10 +4,13 @@ #pragma once #include +#include #include #include "gemmi/symmetry.hpp" +#include "HKLKey.h" +#include "../../common/DiffractionExperiment.h" #include "../../common/Reflection.h" #include "../../common/UnitCell.h" @@ -50,3 +53,32 @@ ReindexChoice ChooseReindex(const std::vector &merged, double ReferenceIntensityCC(const std::vector &merged, const std::vector &reference, int space_group_number); + +// Per-image resolver for serial stills. Each crystal is indexed independently in one of the +// merohedrally-equivalent hands at random, so the ambiguity has to be broken per image (unlike rotation +// data, which resolve it once for the whole merged set with ChooseReindex above). Holds the external +// reference intensities and the crystal's twin-law operators, and reindexes one image's reflections into +// the hand best-correlated with the reference. This runs in the per-image pipeline at integration time, +// decoupled from scaling. +class ReindexAmbiguityResolver { + const ScalingSettings s; + const HKLKeyGenerator hkl_key_generator; + std::map reference_data; + // Twin-law reindexing operators (identity excluded). Empty = no ambiguity (holohedral crystal) or the + // cell/space group is unknown, in which case Resolve() is a no-op. + std::vector ops; + + bool Accept(const Reflection &r) const; + // Pearson correlation of this image's partiality/Lorentz-corrected intensities against the reference, + // matched by ASU index under `op`. Scale-invariant (the unknown per-image scale cancels). NaN when + // too few reflections match. + double ReferenceCC(const std::vector &reflections, const gemmi::Op &op) const; + +public: + // The reference must be an external, correctly-handed dataset (not the data's own running merge). + ReindexAmbiguityResolver(const DiffractionExperiment &x, const std::vector &reference); + + // Reindex this image's reflections in place into the hand best-correlated with the reference (done + // once per image, for good). No-op when there is no ambiguity. + void Resolve(std::vector &reflections) const; +}; diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index a2b3c8f0..175b2aac 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -96,67 +96,13 @@ namespace { }; } -ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref, - bool resolve_ambiguity) +ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref) : s(x.GetScalingSettings()), hkl_key_generator(s.GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)) { for (const auto &r: ref) { const auto key = hkl_key_generator(r); reference_data[key] = r.I; } - if (resolve_ambiguity && x.GetUnitCell().has_value() && x.GetSpaceGroupNumber().has_value()) - ambiguity_ops = ReindexAmbiguityOperators(*x.GetUnitCell(), - static_cast(*x.GetSpaceGroupNumber())); -} - -// Serial stills index each crystal in one of the merohedrally-equivalent hands at random; pick, for this -// image alone, the reindexing whose intensities correlate best with the external reference and apply it. -// The correlation uses the partiality/Lorentz-corrected intensity (scale-invariant, so the unknown -// per-image scale G cancels) matched to the reference by ASU index under the candidate hand. -void ScaleOnTheFly::ResolveIndexingAmbiguity(std::vector &reflections) const { - auto reference_cc = [&](const gemmi::Op &op) -> std::pair { - double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; - size_t n = 0; - for (const auto &r: reflections) { - if (!Accept(r) || r.partiality < s.GetMinPartiality()) - continue; - if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || r.partiality <= 0.0f) - continue; - const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); - const auto it = reference_data.find(hkl_key_generator(h[0], h[1], h[2])); - if (it == reference_data.end()) - continue; - const double x = static_cast(r.I) * SafeInv(r.rlp, 1.0) / r.partiality; - const double y = it->second; - if (!std::isfinite(x) || !std::isfinite(y)) - continue; - sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; - ++n; - } - if (n < MIN_REFLECTIONS) - return {NAN, n}; - const double nd = static_cast(n); - const double cov = sxy - sx * sy / nd; - const double vx = sxx - sx * sx / nd; - const double vy = syy - sy * sy / nd; - return {(vx > 0 && vy > 0) ? cov / std::sqrt(vx * vy) : NAN, n}; - }; - - const double identity_cc = reference_cc(gemmi::Op::identity()).first; - const gemmi::Op *best = nullptr; - double best_cc = identity_cc; - for (const auto &op: ambiguity_ops) { - const double cc = reference_cc(op).first; - if (std::isfinite(cc) && (!std::isfinite(best_cc) || cc > best_cc)) { - best_cc = cc; - best = &op; - } - } - if (best != nullptr) - for (auto &r: reflections) { - const gemmi::Op::Miller h = best->apply_to_hkl({{r.h, r.k, r.l}}); - r.h = h[0]; r.k = h[1]; r.l = h[2]; - } } bool ScaleOnTheFly::Accept(const Reflection &r) const { @@ -223,9 +169,6 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { if (integration_outcome.reflections.empty()) return; - if (!ambiguity_ops.empty()) - ResolveIndexingAmbiguity(integration_outcome.reflections); - auto start = std::chrono::steady_clock::now(); ScaleOnTheFlyResult result{ .B = 0.0, .G = 1.0 }; diff --git a/image_analysis/scale_merge/ScaleOnTheFly.h b/image_analysis/scale_merge/ScaleOnTheFly.h index 63d9e587..67bd5623 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.h +++ b/image_analysis/scale_merge/ScaleOnTheFly.h @@ -5,7 +5,6 @@ #include "HKLKey.h" #include "Merge.h" -#include "ReindexAmbiguity.h" #include "../../common/DiffractionExperiment.h" #include "../IntegrationOutcome.h" @@ -33,22 +32,11 @@ class ScaleOnTheFly { const ScalingSettings s; const HKLKeyGenerator hkl_key_generator; std::map reference_data; - // Twin-law reindexing operators for the merohedral indexing ambiguity (empty = no ambiguity, or - // resolution disabled). Populated only when the ref is a trusted external reference (see constructor). - std::vector ambiguity_ops; bool Accept(const Reflection &r) const; [[nodiscard]] std::pair CalculateGlobalCC(const std::vector &reflections) const; public: - // resolve_ambiguity: when the reference is an external, correctly-handed dataset (not the data's own - // running merge), pick per image the reindexing that agrees best with it - serial stills index each - // crystal independently, so the merohedral ambiguity must be broken per image, not globally. - ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref, - bool resolve_ambiguity = false); - - // Reindex this image's reflections into the merohedral-ambiguity hand that best correlates with the - // external reference (stills; done once per image at integration time, decoupled from scaling). - void ResolveIndexingAmbiguity(std::vector &reflections) const; + ScaleOnTheFly(const DiffractionExperiment &x, const std::vector &ref); void Scale(IntegrationOutcome &integration_outcome) const; void Scale(std::vector &integration_outcome, size_t nthreads = 0) const; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 210ff371..72406ddb 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -825,7 +825,7 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // the crystal has an indexing ambiguity (merohedral: lattice symmetry higher than the Laue group), // pick the reindexing that best correlates with the reference intensities and re-merge in it. The // twin-law reindex is metric-preserving, so only the hkl labels change (the cell is unchanged). - // Stills resolve the ambiguity per image in ScaleOnTheFly, not here. + // Stills resolve the ambiguity per image with ReindexAmbiguityResolver, not here. if (rsm && !config_.reference_data.empty() && result.consensus_cell && experiment_.GetSpaceGroupNumber().has_value()) { const int sg_num = static_cast(*experiment_.GetSpaceGroupNumber());