diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 82598783..d436f2a7 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -120,14 +120,16 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data outcome.lattice_candidate = sym_result.conventional; } - // Multi-lattice search for stills: keep extra lattices that share the - // reference unit cell but differ only in orientation. Quick refinement - // (orientation only) is done later in RefineGeometryIfNeeded. + // Multi-lattice search for stills: store rotations that map the reference + // lattice to each accepted extra lattice. Candidates are materialized later + // in RefineGeometryIfNeeded so they're rooted in the refined (and, for + // monoclinic, reordered) main lattice. if (indexer_result.lattice.size() > 1) { auto ml_latt = MultiLatticeSearch(indexer_result.lattice); for (auto &ml : ml_latt) { - if (outcome.extra_lattice_candidates.size() >= experiment.GetIndexingSettings().GetMaxExtraLattices()) + if (outcome.extra_lattice_rotations.size() >= experiment.GetIndexingSettings().GetMaxExtraLattices()) break; + outcome.extra_lattice_rotations.push_back(ml.rotation_vector); RotMatrix rot(ml.rotation_vector.Length(), ml.rotation_vector.Normalize()); outcome.extra_lattice_candidates.push_back(outcome.lattice_candidate->Multiply(rot)); } @@ -186,6 +188,17 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Monoclinic) outcome.lattice_candidate->ReorderMonoclinic(); + // Rebuild extra-lattice candidates from the refined (and possibly reordered) main + // lattice so they share its cell and obtuse-beta convention. + if (!outcome.extra_lattice_rotations.empty()) { + outcome.extra_lattice_candidates.clear(); + outcome.extra_lattice_candidates.reserve(outcome.extra_lattice_rotations.size()); + for (const auto &rv : outcome.extra_lattice_rotations) { + RotMatrix rot(rv.Length(), rv.Normalize()); + outcome.extra_lattice_candidates.push_back(outcome.lattice_candidate->Multiply(rot)); + } + } + // Quick orientation-only refinement of extra lattices (stills path). // Cell, beam center, detector geometry are taken from the first lattice. if (!experiment.IsRotationIndexing() && !outcome.extra_lattice_candidates.empty()) { diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 8df60d21..0f03ee97 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -36,6 +36,7 @@ class IndexAndRefine { struct IndexingOutcome { std::optional lattice_candidate; std::vector extra_lattice_candidates; + std::vector extra_lattice_rotations; DiffractionExperiment experiment; LatticeMessage symmetry{ .centering = 'P', diff --git a/image_analysis/indexing/MultiLatticeSearch.cpp b/image_analysis/indexing/MultiLatticeSearch.cpp index 3080fc31..d40bc3fb 100644 --- a/image_analysis/indexing/MultiLatticeSearch.cpp +++ b/image_analysis/indexing/MultiLatticeSearch.cpp @@ -3,6 +3,9 @@ #include "MultiLatticeSearch.h" +#include +#include + #include namespace { @@ -33,30 +36,42 @@ namespace { } return R; } -} -CrystalLattice Transform(const CrystalLattice &in_latt, int i) { - CrystalLattice latt = in_latt; - switch (i) { - case 0: - break; - case 1: - latt.FlipSign(0); - latt.FlipSign(1); - break; - case 2: - latt.FlipSign(0); - latt.FlipSign(2); - break; - case 3: - latt.FlipSign(1); - latt.FlipSign(2); - break; - default: - break; + struct BasisOp { + std::array perm; + std::array signs; // each +1 or -1 + }; + + // The 24 right-handed (det=+1) sign+permutation reparametrizations of a basis (a,b,c). + // Identity (perm 0,1,2; signs +,+,+) is first so it's preferred when several ops match + // the reference cell. For low-symmetry cells, only a handful pass is_close; for cubic, + // all 24 may pass and the first-match rule keeps behavior deterministic. + const std::vector &RightHandedOps() { + static const std::vector ops = []() { + std::vector v; + const std::array, 6> perms = {{ + {0, 1, 2}, {1, 2, 0}, {2, 0, 1}, // even + {0, 2, 1}, {2, 1, 0}, {1, 0, 2} // odd + }}; + for (size_t p = 0; p < perms.size(); p++) { + const int parity = (p < 3) ? +1 : -1; + for (int s0 : {+1, -1}) + for (int s1 : {+1, -1}) + for (int s2 : {+1, -1}) + if (parity * s0 * s1 * s2 == +1) + v.push_back({perms[p], {s0, s1, s2}}); + } + return v; + }(); + return ops; } - return latt; + CrystalLattice ApplyBasisOp(const CrystalLattice &in, const BasisOp &op) { + const Coord v[3] = {in.Vec0(), in.Vec1(), in.Vec2()}; + return CrystalLattice(v[op.perm[0]] * static_cast(op.signs[0]), + v[op.perm[1]] * static_cast(op.signs[1]), + v[op.perm[2]] * static_cast(op.signs[2])); + } } std::vector MultiLatticeSearch(const std::vector &lattices, @@ -70,9 +85,8 @@ std::vector MultiLatticeSearch(const std::vector MultiLatticeSearch(const std::vector(rod.x()), - static_cast(rod.y()), - static_cast(rod.z())); - - if (found) - continue; - ret.push_back({rotation_vector}); - found = true; + ret.push_back({Coord(static_cast(rod.x()), + static_cast(rod.y()), + static_cast(rod.z()))}); + break; } }