diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index fbde4042b..f56aa0b26 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/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(scaling_end_time - scaling_start_time).count(); } @@ -990,9 +1000,16 @@ std::optional IndexAndRefine::GetConsensusUnitCell() const { std::vector 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); diff --git a/image_analysis/scale_merge/ReindexAmbiguity.cpp b/image_analysis/scale_merge/ReindexAmbiguity.cpp index ef46588e5..61783520c 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.cpp +++ b/image_analysis/scale_merge/ReindexAmbiguity.cpp @@ -50,6 +50,18 @@ std::vector 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(inv.rot[i][0]) / gemmi::Op::DEN) + + old_star[1] * (static_cast(inv.rot[i][1]) / gemmi::Op::DEN) + + old_star[2] * (static_cast(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 ReindexReflections(const std::vector &merged, const gemmi::Op &op) { std::vector out = merged; @@ -254,7 +266,7 @@ double ReindexAmbiguityResolver::ReferenceCC(std::span 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 reflections) const { +void ReindexAmbiguityResolver::Resolve(std::span reflections, CrystalLattice &latt) const { if (ops.empty()) return; @@ -268,4 +280,10 @@ void ReindexAmbiguityResolver::Resolve(std::span 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); } diff --git a/image_analysis/scale_merge/ReindexAmbiguity.h b/image_analysis/scale_merge/ReindexAmbiguity.h index 55282e232..773011342 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.h +++ b/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 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 ReindexReflections(const std::vector &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 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 reflections, CrystalLattice &latt) const; }; diff --git a/image_analysis/scale_merge/SublatticeRelabel.cpp b/image_analysis/scale_merge/SublatticeRelabel.cpp index ed9b3e0e3..299df94e5 100644 --- a/image_analysis/scale_merge/SublatticeRelabel.cpp +++ b/image_analysis/scale_merge/SublatticeRelabel.cpp @@ -68,6 +68,17 @@ std::vector, 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, 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, 3> &M, @@ -155,16 +166,18 @@ SublatticeRelabelResult RelabelOnSublattice(std::vector &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 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(); +} diff --git a/image_analysis/scale_merge/SublatticeRelabel.h b/image_analysis/scale_merge/SublatticeRelabel.h index 72b3c7895..1dbe0c49b 100644 --- a/image_analysis/scale_merge/SublatticeRelabel.h +++ b/image_analysis/scale_merge/SublatticeRelabel.h @@ -25,6 +25,12 @@ std::vector, 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 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 diff --git a/tests/ReindexAmbiguityTest.cpp b/tests/ReindexAmbiguityTest.cpp index e092c464d..1f821295f 100644 --- a/tests/ReindexAmbiguityTest.cpp +++ b/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(h) + latt.Bstar() * static_cast(k) + + latt.Cstar() * static_cast(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()); +}