rugnux stills: a given cell is matched in any setting, and per-image reindexing moves the lattice with the indices

Three defects that stopped known-cell stills runs:

- The consensus cell compared every frame's cell to -C axis by axis. A run whose frames index in another
  setting of the same lattice (axes in another order, another monoclinic choice) had no frame accepted, no
  consensus cell and no output ("No crystal lattice was determined"), although it merged fine. A frame
  whose lattice contains the given cell now counts on the given cell's axes (CellOnReference, the same
  FindSublatticeBases the merge-stage relabel uses).
- The per-image ambiguity resolver relabelled a frame's reflections but left its lattice on the old axes.
  Everything after it predicts from index and lattice: the sublattice relabel put such a frame through a
  second transformation, and on a monoclinic cell with a near-symmetry (a+c nearly as long as a) the
  merge came out at indices finer than anything measured, shells over 100% complete. The lattice now
  follows the indices (ReindexLattice, B' = R^-1 B).
- With -C/-S the stills path replaced the known-cell indexer's lattice, already on the given axes, by
  LatticeSearch's conventional choice, which on that near-symmetric cell is the other description on half
  the frames. The indexer's lattice is now kept when it is close to the given cell.

Measured (cell + space group given, dimple against the deposited model and free set):
- monoclinic MbCO-type still set: was over-complete/broken with --model; now 2.06 A, R/Rfree 0.233/0.314
  at 2.11 A (P1 merge before: 0.275/0.306).
- C-centred orthorhombic set: no MTZ before; now 3.27 A, R/Rfree 0.298/0.403 at 3.68 A (de novo 0.336/0.448).
- monoclinic two-form set, one form: cell P 1 21 1 with beta (was alpha); Rfree 0.310 -> 0.260, 3.05 -> 2.72 A.
- No change on seven other stills sets, nor on two rotation battery sets (5reo, lyso_x10sa_90deg_1).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

image_analysis/IndexAndRefine.cpp

@@ -16,6 +16,7 @@
 #include "lattice_search/LatticeSearch.h"
 #include "scale_merge/ReindexAmbiguity.h"
 #include "scale_merge/ScaleOnTheFly.h"
+#include "scale_merge/SublatticeRelabel.h"
 
 namespace {
     // Iterations the offline geometry refinement is allowed, standing in for the 40 ms the online path
@@ -339,8 +340,17 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data
                 // niggli_class is left unassigned (0): it needs the primitive cell incl.
                 // centering, which LatticeSearch cannot recover from a (possibly centred, e.g.
                 // C2) user cell. A proper primitive-cell indexing path (CrystFEL-style) is deferred.
+                // A lattice FFBIDX returned on the given cell is already on its axes; LatticeSearch's
+                // conventional choice is only one of the descriptions of it, and on a cell with a
+                // near-symmetry (a monoclinic one whose a+c is nearly as long as a) it picks the other
+                // one on half the frames, which the merge then has to tell apart again.
+                const auto settings = experiment.GetIndexingSettings();
                 auto sym_result = LatticeSearch(latt);
-                if (sym_result.system == sg->crystal_system())
+                if (!de_novo && latt.GetUnitCell().is_close(*experiment.GetUnitCell(),
+                                                            settings.GetUnitCellDistTolerance(),
+                                                            settings.GetUnitCellAngleTolerance_deg()))
+                    outcome.lattice_candidate = latt;
+                else if (sym_result.system == sg->crystal_system())
                     outcome.lattice_candidate = sym_result.conventional;
                 else if (!de_novo)
                     outcome.lattice_candidate = latt;
@@ -938,7 +948,7 @@ void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) {
         return;
 
     auto scaling_start_time = std::chrono::steady_clock::now();
-    reindex_resolver->Resolve(outcome.reflections);
+    reindex_resolver->Resolve(outcome.reflections, outcome.latt);
     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();
 }
@@ -990,9 +1000,16 @@ std::optional<UnitCell> IndexAndRefine::GetConsensusUnitCell() const {
         std::vector<UnitCell> accepted;
         accepted.reserve(cells.size());
 
+        // A frame whose lattice is the given cell in another setting (axes in another order, another
+        // monoclinic or centred choice) counts on the given cell's axes - the merge puts it there
+        // (RelabelOnSublattice) - rather than being left out. Compared axis by axis only, a run whose
+        // frames all index in another setting than the one typed after -C has no consensus at all.
         for (const auto &cell: cells) {
             if (cell.is_close(*experiment.GetUnitCell(), dist_tolerance, angle_tolerance))
                 accepted.emplace_back(cell);
+            else if (const auto on_ref = CellOnReference(cell, *experiment.GetUnitCell(),
+                                                         experiment.GetIndexingSettings()))
+                accepted.emplace_back(*on_ref);
         }
 
         return MeanUnitCell(accepted);

image_analysis/scale_merge/ReindexAmbiguity.cpp

@@ -50,6 +50,18 @@ std::vector<gemmi::Op> ReindexAmbiguityOperators(const UnitCell &cell, const gem
                                  /*all_ops=*/false);
 }
 
+CrystalLattice ReindexLattice(const CrystalLattice &latt, const gemmi::Op &op) {
+    const gemmi::Op inv = op.inverse();
+    const Coord old_star[3] = {latt.Astar(), latt.Bstar(), latt.Cstar()};
+    Coord rs[3];
+    for (int i = 0; i < 3; ++i)
+        rs[i] = old_star[0] * (static_cast<float>(inv.rot[i][0]) / gemmi::Op::DEN)
+              + old_star[1] * (static_cast<float>(inv.rot[i][1]) / gemmi::Op::DEN)
+              + old_star[2] * (static_cast<float>(inv.rot[i][2]) / gemmi::Op::DEN);
+    const float v_star = rs[0] * (rs[1] % rs[2]);
+    return {(rs[1] % rs[2]) / v_star, (rs[2] % rs[0]) / v_star, (rs[0] % rs[1]) / v_star};
+}
+
 std::vector<MergedReflection> ReindexReflections(const std::vector<MergedReflection> &merged,
                                                  const gemmi::Op &op) {
     std::vector<MergedReflection> out = merged;
@@ -254,7 +266,7 @@ double ReindexAmbiguityResolver::ReferenceCC(std::span<const Reflection> reflect
 
 // 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::span<Reflection> reflections) const {
+void ReindexAmbiguityResolver::Resolve(std::span<Reflection> reflections, CrystalLattice &latt) const {
     if (ops.empty())
         return;
 
@@ -268,4 +280,10 @@ void ReindexAmbiguityResolver::Resolve(std::span<Reflection> reflections) const
         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];
     }
+    // Everything after this - the sublattice relabel, the partiality post-refinement - predicts a spot
+    // from its index and the lattice, so the lattice has to describe the new indices. Left on the old
+    // axes, a frame relabelled by an operator that is only a near-symmetry of the lattice was put
+    // through the sublattice relabel a second time, and its reflections came out at indices finer
+    // than anything measured.
+    latt = ReindexLattice(latt, best.op);
 }

image_analysis/scale_merge/ReindexAmbiguity.h

@@ -16,6 +16,7 @@
 #include "../../common/DiffractionExperiment.h"
 #include "../../common/Reflection.h"
 #include "../../common/UnitCell.h"
+#include "../../common/CrystalLattice.h"
 
 // Resolving the indexing ("alternative indexing" / reindexing) ambiguity: when the lattice symmetry
 // is higher than the crystal Laue symmetry (merohedral / pseudo-merohedral), a crystal can be validly
@@ -34,6 +35,10 @@ std::vector<gemmi::Op> ReindexAmbiguityOperators(const UnitCell &cell, const gem
 // A copy of `merged` with each Miller index transformed by `op` (the intensities are unchanged).
 // The indices are left where `op` puts them, which is what the scoring callers want - they reduce to
 // the ASU themselves.
+// The lattice an image has once its indices are moved by `op` (h' = h R): the reciprocal axes follow as
+// R^-1 times the old ones, so h' (a*, b*, c*)' is the same reciprocal vector as h (a*, b*, c*).
+CrystalLattice ReindexLattice(const CrystalLattice &latt, const gemmi::Op &op);
+
 std::vector<MergedReflection> ReindexReflections(const std::vector<MergedReflection> &merged,
                                                  const gemmi::Op &op);
 
@@ -130,5 +135,7 @@ public:
 
     // 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::span<Reflection> reflections) const;
+    // Relabels the image's reflections and moves its lattice onto the same axes, so that a spot's new
+    // index still predicts it where it was measured.
+    void Resolve(std::span<Reflection> reflections, CrystalLattice &latt) const;
 };

image_analysis/scale_merge/SublatticeRelabel.cpp

@@ -68,6 +68,17 @@ std::vector<std::array<std::array<int, 3>, 3>> FindSublatticeBases(const Crystal
 }
 
 namespace {
+    // The reference basis in this frame's orientation: the reciprocal axes are rows of M times the
+    // frame's, and the real axes their reciprocal.
+    CrystalLattice OnBasis(const CrystalLattice &frame, const std::array<std::array<int, 3>, 3> &M) {
+        const Coord fs[3] = {frame.Astar(), frame.Bstar(), frame.Cstar()};
+        Coord rs[3];
+        for (int i = 0; i < 3; ++i)
+            rs[i] = fs[0] * M[i][0] + fs[1] * M[i][1] + fs[2] * M[i][2];
+        const float v_star = rs[0] * (rs[1] % rs[2]);
+        return {(rs[1] % rs[2]) / v_star, (rs[2] % rs[0]) / v_star, (rs[0] % rs[1]) / v_star};
+    }
+
     // The frame's reflections at nodes of the lattice M picks out, re-indexed on it: a frame index h'
     // is h = h' M^-1 = h' adj(M) / det(M) there, a node only where all three are integers.
     ReflectionVector NodesOf(const ReflectionVector &reflections, const std::array<std::array<int, 3>, 3> &M,
@@ -155,16 +166,18 @@ SublatticeRelabelResult RelabelOnSublattice(std::vector<IntegrationOutcome> &out
         result.reflections_kept += best_nodes.size();
         o.reflections = std::move(best_nodes);
 
-        // The reference basis in this frame's orientation: the reciprocal axes are rows of M times the
-        // frame's, and the real axes their reciprocal.
-        const auto &M = *best;
-        const Coord fs[3] = {o.latt.Astar(), o.latt.Bstar(), o.latt.Cstar()};
-        Coord rs[3];
-        for (int i = 0; i < 3; ++i)
-            rs[i] = fs[0] * M[i][0] + fs[1] * M[i][1] + fs[2] * M[i][2];
-        const float v_star = rs[0] * (rs[1] % rs[2]);
-        o.latt = CrystalLattice((rs[1] % rs[2]) / v_star, (rs[2] % rs[0]) / v_star, (rs[0] % rs[1]) / v_star);
+        o.latt = OnBasis(o.latt, *best);
         ++result.frames;
     }
     return result;
 }
+
+std::optional<UnitCell> CellOnReference(const UnitCell &cell, const UnitCell &reference,
+                                        const IndexingSettings &settings) {
+    auto bases = FindSublatticeBases(CrystalLattice(cell), reference, settings);
+    if (bases.empty())
+        return {};
+    // The reference in another setting (index 1) before any axis multiple of it.
+    std::sort(bases.begin(), bases.end(), [](const auto &x, const auto &y) { return Det(x) < Det(y); });
+    return OnBasis(CrystalLattice(cell), bases.front()).GetUnitCell();
+}

image_analysis/scale_merge/SublatticeRelabel.h

@@ -25,6 +25,12 @@ std::vector<std::array<std::array<int, 3>, 3>> FindSublatticeBases(const Crystal
                                                                    const UnitCell &reference,
                                                                    const IndexingSettings &settings);
 
+// A frame's cell described on the reference's axes, when its lattice contains the reference (the
+// reference in another setting, or an axis multiple of it): the cell the frame has once it is put on the
+// reference. Empty where the lattice does not contain it.
+std::optional<UnitCell> CellOnReference(const UnitCell &cell, const UnitCell &reference,
+                                        const IndexingSettings &settings);
+
 struct SublatticeRelabelResult {
     size_t frames = 0;            // frames put on the reference cell
     size_t reflections_kept = 0;  // reflections of those frames at nodes of the reference lattice

tests/ReindexAmbiguityTest.cpp

@@ -11,6 +11,7 @@
 #include "../image_analysis/scale_merge/CrystalSetting.h"
 #include "../image_analysis/scale_merge/HKLKey.h"
 #include "../image_analysis/scale_merge/ReindexAmbiguity.h"
+#include "../image_analysis/scale_merge/SublatticeRelabel.h"
 
 namespace {
     const gemmi::SpaceGroup &SG(int number) { return *gemmi::find_spacegroup_by_number(number); }
@@ -239,3 +240,69 @@ TEST_CASE("Reindex: a near-symmetry law moves the cell with the indices", "[rein
         CHECK(moved);
     }
 }
+
+namespace {
+    // A P2(1) cell with a near-symmetry: a+c is within 1% of a, so the lattice is metrically close to
+    // a C-centred orthorhombic one (obliquity ~0.7 deg) without being one.
+    UnitCell PseudoMonoclinic() { return UnitCell{64.0f, 29.0f, 36.0f, 90.0f, 107.0f, 90.0f}; }
+
+    Coord Recip(const CrystalLattice &latt, int h, int k, int l) {
+        return latt.Astar() * static_cast<float>(h) + latt.Bstar() * static_cast<float>(k)
+             + latt.Cstar() * static_cast<float>(l);
+    }
+}
+
+TEST_CASE("Reindex: an ambiguity operator of a pseudo-symmetric monoclinic cell keeps every reflection's resolution",
+          "[reindex]") {
+    const UnitCell cell = PseudoMonoclinic();
+    const auto laws = ReindexAmbiguityOperators(cell, SG(4));
+    REQUIRE_FALSE(laws.empty());
+    const CrystalLattice latt(cell);
+    for (const auto &op : laws)
+        for (int h = -6; h <= 6; ++h)
+            for (int l = -6; l <= 6; ++l) {
+                if (h == 0 && l == 0) continue;
+                // As the resolver applies it: an operator that is a near-symmetry of the metric moves an
+                // index to one of nearly the same length, never to a finer one.
+                const auto m = op.apply_to_hkl({{h, 3, l}});
+                const float d0 = 1.0f / Recip(latt, h, 3, l).Length();
+                const float d1 = 1.0f / Recip(latt, m[0], m[1], m[2]).Length();
+                CHECK(std::fabs(d1 - d0) / d0 < 0.03f);
+            }
+}
+
+TEST_CASE("Reindex: the reindexed lattice predicts every spot where it was", "[reindex]") {
+    const UnitCell cell = PseudoMonoclinic();
+    const auto laws = ReindexAmbiguityOperators(cell, SG(4));
+    REQUIRE_FALSE(laws.empty());
+    const CrystalLattice latt(cell);
+    for (const auto &op : laws) {
+        const CrystalLattice moved = ReindexLattice(latt, op);
+        for (int h = -4; h <= 4; ++h)
+            for (int k = -4; k <= 4; ++k)
+                for (int l = -4; l <= 4; ++l) {
+                    const auto m = op.apply_to_hkl({{h, k, l}});
+                    const Coord before = Recip(latt, h, k, l);
+                    const Coord after = Recip(moved, m[0], m[1], m[2]);
+                    CHECK((after - before).Length() < 1e-5f);
+                }
+    }
+}
+
+TEST_CASE("Sublattice: a frame on the given cell in another setting is described on the given cell's axes",
+          "[reindex]") {
+    // The given C222(1)-style cell, and the same lattice with b and c swapped (and a flipped to keep the
+    // basis right-handed): not close axis by axis, but the same lattice.
+    const UnitCell given{40.0f, 180.0f, 143.0f, 90.0f, 90.0f, 90.0f};
+    const UnitCell frame{40.2f, 143.4f, 180.5f, 90.0f, 90.0f, 90.0f};
+    IndexingSettings settings;
+    REQUIRE_FALSE(frame.is_close(given, settings.GetUnitCellDistTolerance(), settings.GetUnitCellAngleTolerance_deg()));
+    const auto on_given = CellOnReference(frame, given, settings);
+    REQUIRE(on_given.has_value());
+    CHECK(on_given->is_close(given, settings.GetUnitCellDistTolerance(), settings.GetUnitCellAngleTolerance_deg()));
+    CHECK(on_given->b == Catch::Approx(180.5f).epsilon(0.001));
+    CHECK(on_given->c == Catch::Approx(143.4f).epsilon(0.001));
+
+    // An unrelated cell is not put on it.
+    CHECK_FALSE(CellOnReference(UnitCell{55.0f, 66.0f, 77.0f, 90.0f, 90.0f, 90.0f}, given, settings).has_value());
+}