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) <noreply@anthropic.com>

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<RotationIndexerResult> IndexAndRefine::FinalizeRotationIndexing()
 }
 
 IndexAndRefine &IndexAndRefine::ReferenceIntensities(std::vector<MergedReflection> &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<ScaleOnTheFly>(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<ReindexAmbiguityResolver>(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<float>(scaling_end_time - scaling_start_time).count();
 }

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<ScaleOnTheFly> scaling_engine;
+    std::unique_ptr<ReindexAmbiguityResolver> reindex_resolver;
     void ScaleImage(DataMessage &msg, IntegrationOutcome& outcome);
 
     std::optional<float> RotationAngle(int64_t image) const;   // mid-exposure angle for the indexer

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<gemmi::Op> &ops,
+                                const std::function<double(const gemmi::Op &)> &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<gemmi::Op> 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<MergedReflection> &merged,
                             const UnitCell &cell, int space_group_number,
                             const std::function<double(const std::vector<MergedReflection> &)> &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<int>(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<int>(ops.size());
+    return choice;
 }
 
 double ReferenceIntensityCC(const std::vector<MergedReflection> &merged,
@@ -87,3 +120,66 @@ double ReferenceIntensityCC(const std::vector<MergedReflection> &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<MergedReflection> &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<int>(*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<Reflection> &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<double>(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<double>(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<Reflection> &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];
+    }
+}

image_analysis/scale_merge/ReindexAmbiguity.h

@@ -4,10 +4,13 @@
 #pragma once
 
 #include <functional>
+#include <map>
 #include <vector>
 
 #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<MergedReflection> &merged,
 double ReferenceIntensityCC(const std::vector<MergedReflection> &merged,
                             const std::vector<MergedReflection> &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<HKLKey, double> 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<gemmi::Op> 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<Reflection> &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<MergedReflection> &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<Reflection> &reflections) const;
+};

image_analysis/scale_merge/ScaleOnTheFly.cpp

@@ -96,67 +96,13 @@ namespace {
     };
 }
 
-ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector<MergedReflection> &ref,
-                             bool resolve_ambiguity)
+ScaleOnTheFly::ScaleOnTheFly(const DiffractionExperiment &x, const std::vector<MergedReflection> &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<int>(*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<Reflection> &reflections) const {
-    auto reference_cc = [&](const gemmi::Op &op) -> std::pair<double, size_t> {
-        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<double>(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<double>(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 };

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<HKLKey, double> 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<gemmi::Op> ambiguity_ops;
 
     bool Accept(const Reflection &r) const;
     [[nodiscard]] std::pair<double, size_t> CalculateGlobalCC(const std::vector<Reflection> &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<MergedReflection> &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<Reflection> &reflections) const;
+    ScaleOnTheFly(const DiffractionExperiment &x, const std::vector<MergedReflection> &ref);
 
     void Scale(IntegrationOutcome &integration_outcome) const;
     void Scale(std::vector<IntegrationOutcome> &integration_outcome, size_t nthreads = 0) const;

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<int>(*experiment_.GetSpaceGroupNumber());